Rear-edge populations are important for understanding climate change risk and adaptation potential of threatened species

Climate change disproportionately threatens alpine species, by reducing available habitat and by isolating their populations. These pressures are particularly relevant for rear-edge populations, which typically occupy more marginal habitat compared to populations at the core of species' ranges. We studied Caucasian grouse Lyrurus mlokosiewiczi in the Caucasus ecoregion, a global biodiversity hotspot where this species is endemic, to understand potential climate change impacts on the species. Specifically, we assessed how climate change impacts rear-edge populations and how important these populations are for understanding range shifts and adaptive capacity under climate change. We used maximum entropy modeling to assess changes in the distribution of climatically suitable habitat under present and 2070 climate conditions for the representative concentration pathways 8.5 (RCP8.5). Our results revealed that ignoring rear-edge populations leads to a significant underestimation of the future range (by about 14,700 km 2 ). Rear-edge populations were better adapted to warmer climates compared to core populations, and ignoring them, therefore, also underestimates adaptive capacity. Preventing the loss of rear-edge populations should, therefore, be a priority for conservation planning in the face of climate


| INTRODUCTION
Climate change is a major driver of biodiversity loss and its impact is expected to rise drastically over the next decades (Leadley et al., 2010).Climate change already affects a wide range of taxa in various ways, including via changes in phenology, physiology and morphology, distributional ranges, and changes in abundances (Bellard, Bertelsmeier, Leadley, Thuiller, & Courchamp, 2012;Scheffers et al., 2016).While climate change can be beneficial for some species, climate change impacts are typically negative and a threat to species' survival (Barrows et al., 2020;Iknayan & Beissinger, 2018;Sweet et al., 2019).Understanding potential climate change impacts for species of conservation concern is therefore critically important.
Climate change is predicted to have particularly strong impacts on mountain biodiversity (Pepin et al., 2015;Urban, 2018).Several factors contribute to this.Alpine species are often adapted to specific climate conditions or vegetation belts, and when these conditions move upslope, available habitat shrinks and population sizes decline.This escalates local extinction risk (Scheffers et al., 2016;Sekercioglu, Schineider, Fay, & Loarie, 2008).Moreover, populations of alpine species often occur naturally isolated (e.g., on mountain tops) and climate change threatens to increase isolation further (Ahmadi, Hemami, Kaboli, Malekian, & Zimmermann, 2019).Finally, evolutionary traits of alpine species, such as narrow thermal tolerance, may further increase their vulnerability to climate change (Foden et al., 2013;Hof, Rodríguez-Castañeda, Allen, Jansson, & Nilsson, 2017).Generally, high vulnerability of mountain specialists to climate change is worrisome, because there is typically a high level of endemism among these species (Rahbek et al., 2019), which increases the risk for global extinctions (Dirnböck, Essl, & Rabitsch, 2011).
Rear-edge populations (i.e., current low-latitude or lowelevation populations within species' ranges; Hampe & Petit, 2005) are often particularly at risk.These populations often have smaller thermal-safety margins and, if they cannot shift their range (Bennett, Wernberg, Arackal Joy, de Bettignies, & Campbell, 2015), are expected to experience higher population loss and habitat shrinkage than core populations (Hampe & Petit, 2005).Across the world, rearedge populations already decline in response to climate change (Marqués, Camarero, Gazol, & Zavala, 2016).This is problematic, because it might be exactly these populations that are better adapted to warmer climates compared to more northern or higher-elevation populations in the core distribution of species.Losing rear-edge populations also means losing future adaptive capacity to cope with warming climates (Razgour et al., 2013).Identifying rear-edge populations at risk is therefore important to manage for their persistence and to maintain species' genetic variability (Ashcroft, Gollan, Warton, & Ramp, 2012).
F I G U R E 1 Study area and Caucasian grouse (Lyrurus mlokosiewiczi) occurrence points used in our distribution models, as well as accessible areas from which background points were sampled.The inset map is the location of the study area.Grouse drawing from Gould (1850) The area hosts an amalgamation of flora and fauna from Europe, Central Asia, and the Middle East, with a very large number of endemic species, many of which are relict populations adapted to high-mountain, alpine conditions (Shatberashvili et al., 2015).While there is limited data on the potential impact of climate change on the unique biodiversity of the Caucasus, alpine species are predicted to undergo serious range contractions due to climate warming (Flousek, Telenský, Hanzelka, & Reif, 2015).In addition to climate change, habitat destruction and overexploitation pose further significant threats to the biodiversity of the Caucasus (Shatberashvili et al., 2015).Despite this, assessments of climate change impacts on species of conservation concern are largely lacking, which is a constraint to conservation planning.
Caucasian grouse is likely sensitive to climate change (Hof & Allen, 2019), given that it occurs at the highest elevations of the Caucasus and demonstrates specific traits such as ground-nesting, a long incubation period (20-25 days), a high-level of habitat specialization, and lekking courtship behavior.Indeed, breeding success of the Caucasian grouse strongly depends on weather conditions during brooding and the first days of life of nestlings (Kotov, 1968) and climate fluctuations have been shown to be a main driver of population fluctuations (Vitovich, 1986).Likewise, a major range contraction has been suggested for the species (Hof & Allen, 2019).However, previous work has focused on the entire distributional range of the species.This is potentially problematic, as ignoring intraspecific variation, such as between core and rear-edge populations, may lead to inaccurate predictions of range shifts (Ikeda et al., 2017;Smith et al., 2019).
The two sub-populations known to exist include one in the northern part of the range (the Greater Caucasus) and another in the southern parts (the Lesser Caucasus and adjacent mountain ranges in Turkey and Iran; Potapov & Pavlova, 2009).Habitat of the rear-edge populations (i.e., populations occurring in the southern part of the range) is characterized by heterogeneous topography, and suitable habitat patches are often isolated from each other (Sultanov, Kerimov, Agaeva, Mammadova, & Talibov, 2004).Rear-edge populations also inhabit warmer areas than core populations (Habibzadeh, Storch, & Ludwig, 2019;Nurtaev & Nurtaev, 2016).Given this, we assumed that the climatic niches of the rear edge and core populations would likely differ.To assess how climate change might impact on rear-edge populations of the species and to assess the importance of these populations for understanding range shift potential, we compared species distribution models (SDMs) based on (a) presence locations from the entire distribution of the Caucasian grouse, including core and rear-edge populations, versus (b) from the core population only.Specifically, we aimed to better understand how the predicted current distribution and potential future range shifts under climate change vary when including or excluding rear-edge populations from species distribution models.

| Occurrence data
We used a total of n = 377 Caucasian grouse's presence locations.These included presence locations from Gavashelishvili and Javakhishvili (2010) (n = 350), which were recorded using radio-telemetry and field surveys in 2004-2005 and were representative for habitat where the bird occurs for most of the year.In addition, we used presence locations compiled by Habibzadeh and Rafieyan (2016) representing Caucasian grouse lekking sites for 2013-2014 (n = 27).Lekking sites are visited during the mating season from mid-April to mid-May (Potapov, 2004).Nesting sites are typically close by lekking sites (Etzold, 2005).
Our predictor variables had a resolution of 1 × 1 km 2 (see below) and where multiple occurrence locations fell in a single grid cell, we kept only one record.To further assess spatial autocorrelation among the remaining occurrence points, we used the Mantel test in the R package "ecodist" (Goslee & Urban, 2007).This suggested 2.5 km as the minimum acceptable distance between occurrence locations (Figure 2).After removing presence points less than 2.5 km apart from each other, we retained 191 points (121 and 70 for the Greater Caucasus and the Lesser Caucasus, respectively; Figure 1).This is considered a reasonably large sample size (Zurell et al., 2020).

| Predictor variables
We used bioclimatic variables from the Climatologies at High resolution for the Earth's Land Surface Areas (CHELSA) data set, version 1.2 at 30 arcsec resolution (Karger et al., 2017).This data set is based on observational data from 1979 to 2013 and has a higher predictive power than other global climatologies (Maria & Udo, 2017).To avoid overly complex models, we selected five environmental predictors, according to the life history and ecology of both black grouse Lyrurus tetrix (Table S1), a closely related species, and Caucasian grouse (Gavashelishvili & Javakhishvili, 2010;Habibzadeh et al., 2019;Hof & Allen, 2019).
Black grouse populations have been shown to be limited by temperature and precipitation, with mild winter temperatures and less snow cover adversely affecting grouse population growth rates (Viterbi, Imperio, Alpe, Bosser-peverelli, & Provenzale, 2015).Heavy rain is detrimental to the breeding success of grouse, a groundnesting bird, and cold, wet weather events shortly after hatching increases chick mortality substantially (Hannon & Martin, 2006;Ludwig et al., 2006).Likewise, the reduction of invertebrate prey densities due to high rainfall forces chicks to forage in locations with higher predation risk (Ludwig, Helle, & Siitari, 2010).To reflect these climate conditions potentially impacting Caucasian grouse, we selected the bioclimatic variables temperature seasonality (bio4), minimum temperature of coldest month (bio6), mean temperature of wettest quarter (bio8), precipitation of wettest month (bio13), and annual mean temperature (bio1; Table S1).
We divided these variables into three distinct sets (Table 1) to avoid overfitting and to avoid collinear variable pairs (Spearman's rank correlation jr j ≥ .70).To represent future climate conditions, we used the same CHELSA variables for the year 2070, using outputs of the Community Climate System Model version 4 (CCSM4) for the representative concentration pathway emission scenario 8.5 (RCP8.5).CCSM4 is developed by the National Center for Atmospheric Research (Gent et al., 2011) and is commonly used to study impacts of climate change on biodiversity (Hof & Allen, 2019;Perktas ¸& Elverici, 2020).

| Species distribution modeling
We delimited the area for calibrating models as the portion of the Caucasus potentially accessible to the species (Figure 1) (Anderson & Raza, 2010).To delimit this area, we used the maximum natal dispersal distance of black grouse, a closely related species (29 km; Caizergues & Ellison, 2002) and considered only elevations between 1,300 and 3,300 m, for which the grouse has been observed (Gavashelishvili & Javakhishvili, 2010).We calibrated two types of species distribution models: one with presence locations from the core population only (hereafter: Core model) and one with the presence locations from all Caucasian grouse populations (hereafter: Full model).To estimate and project the distribution of the Caucasian grouse, we used maximum entropy modeling (Maxent; Phillips, Dudík, & Schapire, 2004).Maxent, the most widely used SDM algorithm (Phillips et al., 2004), is a presence-only modeling technique (Elith et al., 2011), and has been consistently shown to outperform alternative approaches (Elith et al., 2006;Merow, Smith, & Silander, 2013;Srivastava, Griess, & Keena, 2020), especially when the number of presence records is small (Filz & Schmitt, 2015).To identify the best Maxent parameterization, we used the "kuenm" package in R (Cobos, Peterson, Barve, & Osorio-Olvera, 2019).We tested candidate models for all combinations of our three predictor sets, 17 values of the regularization multiplier (from 0.1 to 1.0 at intervals of 0.1, 2-6 at intervals of 1, and 8 and 10), and 29 combinations of variable response types (linear, quadratic, product, threshold, and hinge) (Cobos et al., 2019).We randomly split the occurrence dataset into 70% for model calibration and 30% for evaluation.We first removed non-significant candidate models based on partial receiver operating characteristic (ROC) tests (with 100 iterations and 50% of data withheld for bootstrapping).Next, we dropped all models with an F I G U R E 2 Mantel's correlogram calculated with 100 permutations for evaluating spatial auto-correlation for our occurrence dataset of the Caucasian grouse (Lyrurus mlokosiewiczi).White circles show significant spatial autocorrelation after progressive Bonferroni corrections (α = .05,1, 000 permutations) omission rate > 0.01.Finally, we selected the best models based on the second-order Akaike Information Criterion (AICc, Burnham & Anderson, 2002) with ΔAICc ≤2 as criteria.We generated final models for the best model parameterizations using the complete set of occurrence locations and performing 30 bootstrap replicates.Finally, we projected these models to current and future climate scenarios over a study area including the entire current distributional range of the Caucasian grouse as reported by the IUCN Red List of Threatened Species ((BirdLife International, 2016, Figure 1), preventing extrapolation using the clamping function (Stohlgren, Jarnevich, Esaias, & Morisette, 2011).

| Comparison among models with and without the rear-edge populations
We assessed statistical differences in model predictions of suitable habitat with and without the rear-edge populations using the Welch Two Sample paired t-test.We also used the Shapiro-Wilk test to assess for normality.To convert continuous model predictions to binary suitable/unsuitable habitat maps, we applied the maximum training sensitivity plus specificity threshold (Liu, White, & Newell, 2013).To compare current and future predicted habitat between the Core model and the Full model, we calculated the amount of stable area (defined as the percentage of the current habitat suitability still suitable in future), as well as lost and gained suitable area.We performed these analyses using the R packages "raster" and "stars." To assess whether the potentially gained habitats are accessible to extant populations, we classified these areas into (a) those within the dispersal range of extant populations (<29 km); and (b) those beyond this range, which may need assisted colonization.We also described general habitat patterns in areas predicted to experience range expansion according to land-cover shares (i.e., closed forest, open forest, cropland, shrubland, grassland, sparse vegetation) based on the Copernicus 2019 land-cover dataset (Buchhorn et al., 2020) and mean elevation (Danielson & Gesch, 2011) for the Full and Core models.
To assess differences between core (Greater Caucasus) and rear-edge (Lesser Caucasus), we performed univariate pairwise t-test and derived boxplots on the climate variables.Since pairwise comparisons indicated that substantial climate differentiation existed between the two populations (Figure S1), we used niche similarity tests to more formally analyze niche differentiation between the sub-populations.We used two metrics that measure similarity in environmental space: Warren's I (Warren, Glor, & Turelli, 2008) and Schoener's D (Schoener, 1968).Both metrics range between 0 and 1 (1 = complete niche overlap), are typically used together to analyze niche overlap and are both related to standard measures of distance from probability theory (Warren et al., 2008).We compared two metrics, because niche similarity inference can be impacted by the choice of similarity index (Brown & Carnaval, 2019).
Using these two indices, we performed niche overlap and niche divergence tests.These tests compare the niches of two species or populations (core vs. rear-edge populations in our case) by comparing the actual observed values for Warren's I and Schoener's D with a null distribution created from randomly reshuffling occurrence values (Brown & Carnaval, 2019, Supporting Information Methods).To quantify populations' niches, we first used a principal component analysis to characterize the environmental space spanned by the climate variables in two dimensions characterized by the first two principal components.Then, using the resulting principal components, we created a continuous environmentalspace surface using kernel density functions and estimate the occupied environmental-space for each of the two populations (Brown & Carnaval, 2019).We used 500 replicate runs and the same accessible areas and environmental variables as in our species distribution models.Finally, we compared the realized niches of our two populations against the species' fundamental niche, using the Potential Niche Truncation Index (PNTI) in the R package "humboldt".The higher the truncated proportion, the greater the risk is that the realized niche is a poor representation of the fundamental niche (values from 0.15 to 0.30 indicate moderate risk, values >0.30 high risk; Brown & Carnaval, 2019, Supporting Information Methods).All these analyses were calculated using the R package "humboldt" (Brown & Carnaval, 2019).An advantage of this package is that it allows to calculate the environmental niches of species that are out of equilibrium, as can be expected for species with low dispersal ability and/or where distributions have historically been larger (Brown & Carnaval, 2019).

| RESULTS
We tested a total of 1,479 candidate models (i.e., different parameter combinations) based on current climate conditions for each calibration area relevant to the Core population and Full models.All candidate models were statistically significant compared with a null model.Of these models, 219 (14.8%) and 62 (4.2%) models met our omission rate criteria regarding predictive power (E = 1%) for the Full and the Core population models, respectively.Picking those with the lowest complexity level (AICc), resulted in one and two final models for the Full and the Core models, respectively (Table 2).We chose these final models to predict present and future potential Caucasian grouse distributions.In terms of environmental variables, models performed best with the environmental variables in set 1 (Table 2).Regarding variable contributions in our final models, annual mean temperature with 83.5% (±3.9) and 73.5% (±3.1) contributions to the Caucasian grouse habitat suitability across the best final models was the most important variable, and the habitat suitability was highest at 0.02 C (range: −10 to 13.2) and 0.27 C (range: −11.6 to 10.3) for the Full and the Core models, respectively (Figure 3).The precipitation of wettest month with an average 12.1% (±3.1) and 23.0% (±4.3) was the second most important variable in predicting the climate niche of the species for the Full and the Core models, respectively.Surprisingly, Caucasian grouse responded differently to this climate variable in our two models (Figure 3).The lowest contributing variable was temperature seasonality with 4.3% (±2.2) and 3.6% (±1.4) contributions to the species habitat suitability for the Full and the Core models, respectively (Figure 3).
When using only data from the core population (Greater Caucasus), we found substantially less suitable habitat under current climate conditions (21,501 km 2 ) compared to the Full model (48,891 km 2 ; Table 3 and Figure 4).Likewise, projecting to future conditions yielded a different distribution and area of suitable habitat (t = 11.17,p = <.001;Table 3).For our climate scenario (CCSM4; RCP8.5), we projected 24,576 km 2 (2.89% of the study area) and 10,613 km 2 (1.25% of the study area) for the Full and the Core models, respectively.The Full model predicted suitable habitat to expand on 1.14% (9,650 km 2 ) of the study area, whereas the Core model predicted less expansion of habitat (0.76%; 6,450 km 2 ; Table 3 and Figure 4).The Core model predicted the species will lose a larger portion of  contemporary suitable habitat (81.82% of 21,501 km 2 ; 17,339 km 2 ) than under the Full model (70.60% of 48,891 km 2 ; 33,963 km 2 ).More habitats were predicted to remain stable until 2070 in our Full model (29.40%; 14,927 km 2 ) compared to the Core model (18.20%; 4,162 km 2 ; Table 3 and Figure 4).Our findings also revealed that the most gained habitats (84.0 and 80.0% for the Core and Full models, respectively) are within the dispersal range of the species (i.e., 29 km).Gained areas are largely characterized by grassland and sparse vegetation, with mean elevations >3,150 m (Table S2).
Comparing the niches of the core (Greater Caucasus) and the rear-edge populations (Lesser Caucasus) using our niche similarity tests demonstrated a small but significant difference between the two occupied climatic niches (Niche Divergence Test [NDT]: D=0.381, p=.044, Table 4, and Figures S2 and S3).However, the two populations' niches were not significantly different in relation to the total accessible environmental space occupied by them (Niche Overlap Test [NOT]: D= 0.5, p=.367,Table 4, and Figure S2 and S3).The PNTI, which measures the potential for a species' occupied environmental-space to be truncated by the available environmental-space in its environment, varied between core and rear-edge populations (0.0 and 0.15 respectively), signaling low and moderate values of potential niche truncation; Figure S3) when the shared analogous environmental space integrated in the niche similarity analysis.

| DISCUSSION
Climate change is a rapidly increasing pressure on mountain biodiversity, and rear-edge populations are particularly at risk.Mountain ecosystems in mid-latitudes are strongly affected as they experience faster warming rates near the annual 0 C isotherm due to positive feedback loops where snow melt leads to lower albedo (Pepin & Lundquist, 2008).For species of conservation concern in such mountain systems, it is therefore important to understand how climate change might impact the distribution of rear-edge populations, and how rear-edge populations differ from core populations.Focusing on the Caucasus Mountains, a global biodiversity hotspot, we assessed possible range shifts of Caucasian grouse, a species at risk in from warming.Specifically, we tested range shift estimates based on data from the species' core populations (Greater Caucasus) versus additionally including data from their current rear-edge populations (Lesser Caucasus).Our analyses demonstrate that ignoring rear-edge populations in predictions of range shifts leads to a substantial underestimation of the species' future potential range and overestimates range contraction.Niche differences between core and rear-edge populations further suggest that rear-edge populations harbor important adaptive capacity to cope with warmer climates, which might be lost when these populations vanish.Preventing the loss of rear-edge populations should therefore be a priority for conservation planning.Key conservation measures to reach this goal should include connectivity planning to facilitate range shifts, as well as assisted colonizations.More generally, our study shows how niche modeling can be used to identify populations particularly valuable for maintaining high adaptive capacity in the face of climate change.
Our first main finding was that we detected much less climatically-suitable habitat when rear-edge populations were ignored, both now and in the future.The results of our SDMs together with the niche similarity analysis suggested that rear-edge populations and core population show different responses to climatic conditions.This provides further evidence that the Caucasian grouse's realized climate niche differs markedly between the core and rear-edge populations (Habibzadeh et al., 2019).The core population largely occurs in a single vegetation formation, Caucasus mixed forests (Olson et al., 2001), whereas rear-edge populations utilize a wider range of vegetation formations, including steppes and open woodlands.Our results also indicate that rear-edge populations demonstrate a broader climatic niche with respect to precipitation during the wettest month (Figure 3, Table S1).All this corroborates the assumption that the Caucasian grouse might be able to adapt better to warming climate than would be predicted based on the climate envelop of the core population only.Moreover, it suggests that the impact of climate change on the species as a whole might be overestimated when only considering the core population.
Our second major finding was that climate change is still likely to result in a major habitat suitability contraction of the Caucasian grouse in the 21st century (Fitzpatrick, Gove, Sanders, & Dunn, 2008;Trisos, Merow, & Pigot, 2020).A major shrinkage of the Caucasian grouse's distribution across the entire Caucasus due to climate change has previously been suggested (Hof & Allen, 2019), but we here highlight that without taking into account the rear-edge populations, this shrinkage can be overestimated.Several mechanisms could explain range contraction and collapse under climate change.
The Caucasian grouse is a specialist for the alpine meadow-forest ecotone (Klaus et al., 2003), which shifts upslope when temperatures rise.However, there might be a time lag in meadow establishment, and available area generally shrinks as vegetation belts move upslope, thereby diminishing the suitable habitat of the Caucasian grouse.Such upward shifts in distribution, acting as ecological traps, have been observed for montane species throughout the world (Grytnes et al., 2014;Rumpf, Hülber, Zimmermann, & Dullinger, 2019), including alpine birds (Maggini et al., 2011).A second key mechanisms could be that warmer climates reduce food availability for breeding birds, as demonstrated for related red grouse (Lagopus scoticus) in the Scottish Highlands (Fletcher, Howarth, Kirby, Dunn, & Smith, 2013).
Assessing possible climate-change effects on the timing of food availability for the Caucasian grouse, especially during the critical nesting and brooding phases (Kotov, 1968) is an essential follow-up work of our study.Finally, range shifts toward areas with a large share of grassland and sparse vegetation might increase the exposure of grouse to other anthropogenic threats, such as livestock grazing or disturbances from tourism (BirdLife International, 2016).Unanimously, our study highlights the importance of the rear-edge populations for maintaining the adaptive capacity of Caucasian grouse in the face of climate change.Survival of the rear-edge populations of the Caucasian grouse depends on the species' ability to track their climate niche or to adapt to warmer conditions (Des Roches et al., 2018).The latter, however, is an unlikely scenario for Caucasian grouse considering that similar, but warmer mountain environments exist to the south of the Lesser Caucasus where the species does not occur.The long-term stability of the rear-edge populations of the Caucasian grouse may also indicate a high level of local adaptation, which fits the population-specific climate affinity demonstrated here and elsewhere (Habibzadeh et al., 2019).Together, these considerations highlight that particularly the migration of individuals from the rear-edge populations of Caucasian grouse toward its core distribution might hold important potential for the species to cope with climate change.However, several ecological characteristics of the Caucasian grouse make this questionable, such as naturally low densities, often small and declining population size (BirdLife International, 2016;Gokhelashvili, Kerry, & Gavashelishvili, 2003) and a limited dispersal capability (Hof & Allen, 2019).Together, this suggests that the species is in growing need of conservation attention and active management to foster its survival in the Caucasus.
We used species distribution modeling, using the largest sample of occurrence data collected for the species and testing a wide range of model parameterization.Still, a few limitations should be mentioned.First, we used only one species distribution modeling algorithm (Maxent), one climate model (CCSM4), and one emissions scenario (RCP8.5).Using different algorithms, climate models, and emission scenarios would add further nuance to our analyses and would quantify variability in predictions.However, this would likely not change the main trends we find while increasing the complexity of the modeling exercise considerably.Moreover, a recent study showed no particular benefit of ensemble models based on multiple algorithms compared to predictions from a single algorithm (Hao, Elith, Lahoz-Monfort, & Guillera-Arroita, 2020).Second, more presence locations are always desirable, but are hard to gather for this species, due to the rugged and remote nature of the Caucasus.Yet, we note that Maxent has been shown to be particularly suitable for small sample sizes (Filz & Schmitt, 2015).Third, we used a geographic approach to distinguish between core and rear-edge populations.While our tests for differences in habitat use supported such a differentiation, a more dynamic classification (e.g., based on climate thresholds) could have been used.Finally, although we focused on important climate predictors, it could be important to further explore the effect of snow conditions on the species distribution.Snow burrowing enables grouse species to better cope with cold temperatures (Bocca, Caprio, Chamberlain, & Rolando, 2014;Shipley, Sheriff, Pauli, & Zuckerberg, 2019) and climate change may impact the quantity and quality of available snow cover for the Caucasian grouse (Gottschalk, Ekschmitt, _ Isfendiyaroglu, Gem, & Wolters, 2007).In addition, other ecological factors such as competition, dispersal ability, resource distribution or predation could be important to include in order to fully understand potential range shifts (Engelhardt, Neuschulz, & Hof, 2020).As species distribution models that specifically consider population demographics and dispersal become available, applying such models to predict the persistence of grouse under future conditions would be a useful next step (Benito Garzón, Robson, & Hampe, 2019).This, however, would require substantial biological and ecological data to robustly parameterize such models.

| Conservation planning and management implications
Several clear management implications derive from our work.First, both the uncertain future of the rear-edge populations under climate change, and the uniqueness of these populations, translate into an urgency for lowering prevailing threats to Caucasian grouse.Threats to the rearedge population, such as habitat degradation and poaching, should be lowered to avoid future synergistic impacts with climate change (Brook, Sodhi, & Bradshaw, 2008;Mantyka-Pringle, Martin, & Rhodes, 2012).For instance, communitybased ecotourism and bird watching opportunities can ascribe value to the grouse population and reduce poaching incentives.Moreover, habitat degradation due livestock overgrazing is widespread in the Caucasus, even inside protected areas (Habibzadeh & Rafieyan, 2016).Therefore, sustainable land management schemes, particularly related to livestock grazing, would reduce degradation of grouse habitat.All these measures would reduce overall pressure on grouse populations, thereby helping them to cope with the emerging challenges that climate change brings about.
More directly addressing the threat of the rear-edge populations getting trapped and increasingly running out of suitable habitat requires facilitating the species' migration possibilities.A key measure here can be the protection of potential dispersal corridors, particularly between protected areaswhich should be a priority to reduce the negative impacts of climate change on alpine biodiversity in general.As Caucasian grouse populations are currently highly fragmented and spread across several countries, designing a transboundary protected area network, including corridors, should therefore be a conservation priority.Where corridors are not enough, assisted colonization, duly considering the ethical, policy, and scientific challenges related to the this conservation measure (Schwartz et al., 2012), and reintroductions using captive-bred individuals from the rear-edge populations can be viable options (Heinrichs, McKinnon, Aldridge, & Moehrenschlager, 2019).Importantly, facilitating the range shift of the rear-edge populations would make important contributions to helping the core population to better cope with the warmer climates they will face soon.
T A B L E 1 Three sets of bioclimatic variables used to predict the present and future potential distribution of the Caucasian grouse (Lyrurus mlokosiewiczi)

Note:
Model performance was assessed based on statistical significance (partial ROC), omission rates, and model complexity (corrected Akaike information criterion; AICc).FCs are as follows: l = linear, p = product, q = quadratic, and t = threshold.F I G U R E 3 Response curves for the most influential bioclimatic variables showing the relationships between climate suitability of the Caucasian grouse (Lyrurus mlokosiewiczi) and the predictors for the Full model (top panel) and the Core model (down panel).Red curve shows the average over the 30 replicate runs of the best models and grey shadings represent 95% confidence intervals

F
I G U R E 4 Change in the area of climatically suitable habitat for the Caucasian grouse (Lyrurus mlokosiewiczi) for the (a) Full model and the (b) Core model under Representative Concentration Pathway (RCP) 8.5 and employing outputs of the Community Climate System Model version 4 (CCSM4) by 2070 Comparison of suitable Caucasian grouse (Lyrurus mlokosiewiczi) habitat between the Full and Core species distribution model for 2070, as well as changes compared to the contemporary situation Values shown are averages over 30 replicate runs, with standard deviations in brackets.t signifies the outcome of the Welch Two Sample paired t-tests between the Core and the Full model.
T A B L E 3Note: T A B L E 4 Results of the niche divergence and overlap tests between the rear-edge (Lesser Caucasus) and core populations (Greater Caucasus) of the Caucasian grouse Note: D corresponds to the niche similarity index quantified with Schoener's D metric.E is the observed significance of the equivalence statistic.B 2 ! 1 and B 1 ! 2 correspond to background statistics comparing rear-edge and core populations, respectively.