Climate Change and Alpine Screes: No Future for Glacial Relict Papaver occidentale (Papaveraceae) in Western Prealps

: Glacial relicts, especially those with very narrow habitat requirements, are particularly a ﬀ ected by global warming. We considered Papaver occidentale , a glacial relict endemic to the Western Prealps, belonging to the alpine poppy complex ( P. alpinum aggr.), as a model taxon to study the actual status and potential future distribution of species restricted to particular microrefugia. For this study, all known localities were visited, each population was georeferenced and the number of individuals was estimated. Species Distribution Modelling (SDM) was used to evaluate the present and future potential distribution range and habitat suitability, taking into account the speciﬁcity of its habitat (calcareous screes). According to our study, there are globally 19 natural populations of P. occidentale , and a total of about 30,000 individuals. The taxon is a highly specialized alpine plant growing in the majority of natural sites between 1900 and 2100 m a.s.l. on north-facing screes. Predictions for the end of the 21st century indicate that a suitable area will signiﬁcantly decrease (0–30% remaining). Under the most severe climatic scenarios (RCP 8.5), the species risks complete extinction. The long-term in situ conservation of P. occidentale , and all other taxa of the P. alpinum complex, is unlikely to be achieved without slowing global climate change. More generally, our ﬁne-scale study shows that local environmental bu ﬀ ering of large-scale climate change in high-mountain ﬂora may be very limited in specialised taxa of patchy environments such as screes.


Introduction
Understanding how species respond to climate change is a major challenge of conservation biology [1].This is especially the case for high mountain plants, and more specifically glacial relicts inhabiting European Alps [2][3][4].Relicts are remnants of past populations that have become fragmented by climate-driven changes and habitat loss [5].These remnants were left behind during past range shifts and can persist today only in enclaves with suitable environmental conditions in areas with inhospitable regional climates [6].The Alps, along with the neighboring mountain ranges, played an important role in forming the biogeographical patterns in Europe and acted as a refugium for many taxa throughout several ice age cycles [7][8][9][10][11][12].During recent decades, the glacial and postglacial histories of arctic-alpine plant species have been extensively investigated [13][14][15][16][17][18].However, many relict taxa of the alpine region still await thorough analysis as only a few punctual studies have explored the conservation status of the glacial relicts in Europe [19,20].Similarly, although glacial relicts are expected to respond dramatically to climate warming, very few studies exist on their reaction to global warming in their current microrefugia [1,21].
In our study, we investigated a typical glacial relict Papaver occidentale (Markgr.)H.E. Hess et al. (Papaveraceae) (Figure 1).The species is a narrow endemic of the Western Prealps, growing in the alpine zone in plant associations typical of calcareous screes (Thlaspion rotundifolii) [22][23][24][25].Papaver occidentale belongs to the monophyletic section Meconella including ca.50 poppy species, that have mainly circumpolar arctic distribution with some geographically isolated species distributed in the alpine zone of temperate regions [26,27].This group hosts some of the most northerly-growing plants in the world, like Papaver radicatum, which grows at latitudes up to 83 • North [28,29].More specifically, P. occidentale belongs to the alpine poppy complex Papaver alpinum s.l., which is distributed in small isolated regions all across the Alps and neighbouring mountain ranges [24,30,31].Because of the highly polymorphic character of P. alpinum complex, its taxonomic division is not fully resolved [31], it illustrates perfectly the complicated evolutionary and biogeographic history of the Alps.Nevertheless, it was demonstrated that P. occidentale forms a natural and probably monophyletic group with three other white-flowered and narrow-lobed alpine poppies, that are confined to calcareous bedrocks: P. alpinum s.str., P. sendtneri (both north-eastern Alps), and P. tatricum (Tatra mountains in the Carpathians) [30][31][32][33][34].The recolonization history of the Prealps by this taxon after the Last Glacial Maximum (about 20 kya BP) was complex and the present genetic structure of P. occidentale is likely a result of a mixture of local survival and immigration [34].Despite its priority status in many regions of Switzerland and France [35] and its symbolic value to some regions [36], the plant was never studied across its entire distribution area.Thus, besides very few local and technical reports [37][38][39], many aspects of its biology and biogeography (e.g., chorology, sensitivity to global change) have not been studied at all.
Papaver occidentale grows exclusively on steep screes and in inaccessible summit parts of the alpine zone [36].Thus, its habitat is usually not disturbed by any direct human activity (Figure 1).However, the plant is perennial and stenothermic, making it vulnerable to climate change [38].Species Distribution Modelling (SDM) combines species' current natural occurrence data and environmental variables to estimate the geographic range of suitable habitats in the past, present, and future [40][41][42][43].It has been widely applied in conservation biology to predict potential shifts in geographic ranges in response to climate changes, e.g., [44].However, the predictions and use of SDM for alpine plants are challenging due to a complex topographic and edaphic pattern of the high-altitude areas [45].This is even more the case for plants inhabiting very specific alpine habitats, such as calcareous screes, which is the case of P. occidentale.
The present study aims at answering the following questions: (1) What is the present global distribution of Papaver occidentale?(2) What are the population sizes in different regions across its distribution area?(3) What are the ecological preferences of the species?Furthermore, we aimed at applying Species Distribution Modelling (SDM) based on the collected occurrence data, to (4) evaluate the present potential distribution range, and (5) predict the location of suitable habitats of P. occidentale in the future, considering not only various climate change scenarios but also its unique habitat specificity.Finally, based on our findings, we discuss the conservation implications of our results for the model taxon as well as for other relict plants, especially those highly specialized for alpine screes.

Study Area and Data Collection
The studied area covered the whole distribution of Papaver occidentale, namely the Western Prealps between Switzerland and France.Information about all known locations of P. occidentale was obtained from Info Flora (www.infoflora.ch), as well as elucidated from the literature [32,36], from the unpublished reports [38,39], and from personal communications of local experts (B.Clément from the Botanical Garden of the University of Fribourg and D. Jordan from France).The populations from the Swiss canton of Lucerne were identified as belonging to P. sendtneri [34] and were not taken into consideration in the present study.
All populations of P. occidentale were visited and investigated between July and September 2018 (Figures 1 and 2, Table S1, Supplementary Materials).The location of each population was georeferenced directly using field GPS receivers.The number of individuals was estimated as follows: for large populations (>100 individuals), several GPS locations were recorded (10 to 50 m distance between points) with each time the estimation of the number of individuals (direct counting or, when the plants are too dense, counting on a small area with an extrapolation).For small populations (<100 individuals) one GPS location was measured and all individuals counted.Because of the field characteristics (moving screes), it was not possible to use a standard procedure (like a grid pattern) to collect the data.The proximity of high cliffs and strong slopes could have decreased the GPS precision, so all GPS locations were checked afterward on a GIS software with orthophotos, and corrected if necessary (final horizontal error < 20 m).For SDM, every single individual is considered as an occurrence (presence points, the methodology described in Table S2, Supplementary Materials).

Study Area and Data Collection
The studied area covered the whole distribution of Papaver occidentale, namely the Western Prealps between Switzerland and France.Information about all known locations of P. occidentale was obtained from Info Flora (www.infoflora.ch), as well as elucidated from the literature [32,36], from the unpublished reports [38,39], and from personal communications of local experts (B.Clément from the Botanical Garden of the University of Fribourg and D. Jordan from France).The populations from the Swiss canton of Lucerne were identified as belonging to P. sendtneri [34] and were not taken into consideration in the present study.
All populations of P. occidentale were visited and investigated between July and September 2018 (Figures 1 and 2, Table S1, Supplementary Materials).The location of each population was georeferenced directly using field GPS receivers.The number of individuals was estimated as follows: for large populations (>100 individuals), several GPS locations were recorded (10 to 50 m distance between points) with each time the estimation of the number of individuals (direct counting or, when the plants are too dense, counting on a small area with an extrapolation).For small populations (<100 individuals) one GPS location was measured and all individuals counted.Because of the field characteristics (moving screes), it was not possible to use a standard procedure (like a grid pattern) to collect the data.The proximity of high cliffs and strong slopes could have decreased the GPS precision, so all GPS locations were checked afterward on a GIS software with orthophotos, and corrected if necessary (final horizontal error < 20 m).For SDM, every single individual is considered as an occurrence (presence points, the methodology described in Table S2, Supplementary Materials).

Environmental Variables
The statistical analyses and species distribution modelling (SDM) was conducted in the Swiss Western Prealps, from the Lake of Geneva to the Lake of Thun (the centre of P. occidentale distribution).Populations in France were excluded because it was not possible to have access to the environmental variables with the same resolution and homogeneity as those of Switzerland.
We selected a set of 14 environmental variables, which are likely to play a role in P. occidentale distribution (Table 1).The most highly reliable data for the studied area, with the finest available resolution, were used.These data were gathered from different sources.The resolution of the environmental variables used for the SDM was 10 m (Table 1).
Climatic and topographic variables are the basis of SDM, but soil or geology plays also an essential role, mainly for plants [46].As P. occidentale is exclusively linked to moving calcareous screes, this azonal factor plays an essential role to model its distribution.As no exhaustive data exist, we developed this layer ourselves, by delimiting area mainly using orthophotos and in-field knowledge.
The variable EASTNESS (west-east slope aspect) was excluded as it was assumed as not relevant and biased as the mountain chains are mainly oriented from SW to NE in the studied area (on clockwise direction).
The variable ELEVATION was not used in the model because it is an indirect variable [42,47].It was substituted by the temperature (TEMP), which is a direct variable and more suitable to make predictions with future climatic scenarios.The original data for temperature (1 km resolution) was highly correlated with ELEVATION (R-squared value = 0.92, p-value < 0.001).TEMP was thus calculated as a function of ELEVATION, to downscale the resolution to 10 m.
We checked for multicollinearity problems between variables, which is assumed as an important source of model uncertainty.It is safer to avoid selecting highly correlated variables for SDM [47][48][49][50].The variance inflation factor (VIF) was calculated and we looked at paired correlation plots and r-squared values (Table S2 [47,[51][52][53]).With the retained variables, all VIF scores were under 5, (recommendation for value below 10, [47]) and paired correlation R-squared values did not exceed 0.6 (Figures S1-S3, Supplementary Materials), which is below an acceptable threshold of 0.8 [54].
Our first tests showed that SCREES (presence/absence of calcareous screes) was the most important variable, which explained above all the distribution of P. occidentale and largely hid other variable effects.As P. occidentale occurs only in calcareous screes, we, therefore, considered only the area with the presence of calcareous screes for the SDM (and SCREES was thus removed from variables in the models).Pseudo-absences (or background data) were generated inside this area.The number of pseudo-absences was selected to equal the number of presences.Pseudo-absences were randomly generated, which is the strategy with the least assumptions that should be used by default if there are no strong arguments in favour of a different approach [47,55,56].
For future predictions under different climate change scenarios, we gathered raster data from Meteoswiss ([57], Table S3, Supplementary Materials).We used the data that includes the evolution of temperature and precipitation (annual, summer, winter) for the end of the century (sometime between the year 2070 and 2099, noted hereafter year 2085).Three different representative concentration pathway (RCP) scenarios were considered: RCP 2.6 (implies a strong reduction of greenhouse gas emissions early in the 21st century), RCP 4.5 (emissions decline after 2050, stabilization of radiative forcing), and RCP 8.5 (continuously increasing radiative forcing) [57,58].
On average, the annual mean temperature is expected to increase in Switzerland [57], from +1.3 • C (RCP 2.6) to +4.4 • C (RCP 8.5).The actual mean precipitation in winter (PREC-Win data) was also replaced by the raster projections for 2085.On average, the mean winter precipitations are expected to increase in Switzerland from +6% (RCP 2.6) to +13% (RCP 8.5).On the contrary, summer precipitations are expected to decrease (but this variable was not retained in the final validated model).Area delimited using GEOCOVER (Geological vector datasets for better subsurface management, [61]) and interpretation of orthophotos (www.swisstopo.ch,Google maps).A 100 m buffer was added to the delimited areas, as the precise limit of the scree is difficult to get.
NA, vector data.Rasterized to 10 m

Model Selection and Calibration
Two logistic regression techniques were used for SDM: the general linear models (GLMs) and the general additive models (GAMs).GLMs are a generalization of ordinary linear regression.They allow the linear model to be related to the response variable via a link function [62].For our purpose, the link function was the logistic function (logit) and we used quadratic functions of all environmental variables to allow for nonlinear relationships between the covariates and response variable.GAMs are GLMs in which the usual linear relationship between the response and covariates are replaced by several nonlinear smooth functions to model and capture the nonlinearities in the data [63][64][65].
All preselected environmental variables were first included in the models.Then, the significance of the variables in the model, the explained deviance, the AIC, Akaike information criterion, and BIC, Bayesian information criterion [66,67] were checked to remove superfluous variables and build the most parsimonious model.

Model Validation
Different cross-validation techniques have been used to validate the model and prevent overfitting, like 10-fold cross-validation [47,68,69] and Monte-Carlo cross-validation [70,71].However, our data were strongly spatially autocorrelated.Consequently, the 10-fold and the Monte-Carlo cross-validations were not suitable to validate our models (for example, the receiver operating characteristic (ROC) curves of validation datasets were always quasi identical to the ROC curve of the full model, see Figure S4, Supplementary Materials).We, therefore, designed a spatial cross-validation [72,73].The studied area was split into two blocks on the W-E axis.The data within one block serves to train the model while the data in the second block serves to validate the model and vice versa.Then the same exercise was done with two blocks on the N-S axis.Presences and pseudo-absences were each time equilibrated.This validation procedure serves us mainly to delete environmental variables from the model that conduct to overfitting.
The results of cross-validation were analysed with a ROC curve [74].The true positive rate was plotted against the false positive rate for each validation dataset, next to the full model ROC curve, which permits to investigate for overfitting problems.With validation datasets, the area under the curve (AUC) and mean squared residuals (MSR) served to evaluate the model performance.
The spatial cross-validation revealed strong overfitting of the models (both GLM and GAM).The models were thus simplified step-by-step (suppression of some environmental variables, decreasing the degree of smoothing in GAM or suppression of some quadratic terms in GLM), to optimize the dataset during the spatial cross-validation.We lastly also tried to add interaction terms, but they never improved the model and were thus not included.All calculations and graphs were done in R ( [75], Table S4, Supplementary Materials).Maps and spatial analyses were created with QGIS [76].

Distribution and Population Size of Papaver occidentale
The species could be confirmed in all sites indicated in the available sources (e.g., old and recent literature, reports, personal communication of local experts), and no new populations were discovered.Papaver occidentale possesses globally 21 populations (Figures 1 and 2, Table S1, Supplementary Materials) growing in six administrative units: 6 populations in the canton of Fribourg, 5 in Bern, 4 in Vaud, and 1 in Valais, as well as 5 populations in Haute-Savoie (France).Two populations of the canton of Vaud from the Jura Mountains are introduced (TEN and CUN, ancient introductions, [38]).Thus, the majority of 19 natural populations (74%) are growing in Switzerland.
The total population size of natural occurrences of P. occidentale was estimated at 30,756 individuals (estimation from 153 GPS locations).The two introduced sites in the Jura Mountains are very poor in individuals (TEN: 6, CUN: 1).The largest population grows in the canton of Bern (Hinderi Fromatt), with 9851 individuals.Generally, the canton of Bern with 22,117 (72%) of all individuals taken together, is the main distribution center of the taxon.Canton of Fribourg has 5167 (17%), Vaud 748 (2.5%) and Valais 141 (0.5%) individuals.Haute-Savoie, and thus France, counts 2547 individuals, and thus merely 8% of all individuals of the species (Figures 1 and 2).
The species grows in the majority of natural sites between 1900 and 2100 m a.s.l. and on the northern slopes (Figure 3, Table S1, Supplementary Materials), usually between 100 and 300 m from the highest summit.The highest mean elevation of P. occidentale occurrence (2400 m a.s.l.) was measured in the population on Pointe Blanche (PTB) in Haute-Savoie as well as in the canton of Fribourg, with 2300 m a.s.l. on Pointe de Paray (PAR) and with 2235 m a.s.l. at Combe Vanil de l'Ecri (COM).The natural populations with the lowest mean altitude occur in the canton of Vaud, with 1760 m a.s.l. in Les Salaires (SAL) and with 1653 m a.s.l. in La Pierreuse (PIE).The two very small and introduced populations of Jura are growing also at very low altitudes (TEN: 1658 m a.s.l and CUN: 1571 m a.s.l.) (Table S1, Supplementary Materials).
Diversity 2020, 12, x FOR PEER REVIEW 9 of 19 individuals taken together, is the main distribution center of the taxon.Canton of Fribourg has 5167 (17%), Vaud 748 (2.5%) and Valais 141 (0.5%) individuals.Haute-Savoie, and thus France, counts 2547 individuals, and thus merely 8% of all individuals of the species (Figures 1 and 2).The species grows in the majority of natural sites between 1900 and 2100 m a.s.l. and on the northern slopes (Figure 3, Table S1, Supplementary Materials), usually between 100 and 300 m from the highest summit.The highest mean elevation of P. occidentale occurrence (2400 m a.s.l.) was measured in the population on Pointe Blanche (PTB) in Haute-Savoie as well as in the canton of Fribourg, with 2300 m a.s.l. on Pointe de Paray (PAR) and with 2235 m a.s.l. at Combe Vanil de l'Ecri (COM).The natural populations with the lowest mean altitude occur in the canton of Vaud, with 1760 m a.s.l. in Les Salaires (SAL) and with 1653 m a.s.l. in La Pierreuse (PIE).The two very small and introduced populations of Jura are growing also at very low altitudes (TEN: 1658 m a.s.l and CUN: 1571 m a.s.l.) (Table S1, Supplementary Materials).
Diversity 2020, 12, x FOR PEER REVIEW 9 of 19 individuals taken together, is the main distribution center of the taxon.Canton of Fribourg has 5167 (17%), Vaud 748 (2.5%) and Valais 141 (0.5%) individuals.Haute-Savoie, and thus France, counts 2547 individuals, and thus merely 8% of all individuals of the species (Figures 1 and 2).The species grows in the majority of natural sites between 1900 and 2100 m a.s.l. and on the northern slopes (Figure 3, Table S1, Supplementary Materials), usually between 100 and 300 m from the highest summit.The highest mean elevation of P. occidentale occurrence (2400 m a.s.l.) was measured in the population on Pointe Blanche (PTB) in Haute-Savoie as well as in the canton of Fribourg, with 2300 m a.s.l. on Pointe de Paray (PAR) and with 2235 m a.s.l. at Combe Vanil de l'Ecri (COM).The natural populations with the lowest mean altitude occur in the canton of Vaud, with 1760 m a.s.l. in Les Salaires (SAL) and with 1653 m a.s.l. in La Pierreuse (PIE).The two very small and introduced populations of Jura are growing also at very low altitudes (TEN: 1658 m a.s.l and CUN: 1571 m a.s.l.) (Table S1, Supplementary Materials).

Species Distribution Modelling
A total of 28,209 individuals (presences) were comprised in the area considered for the SDM (representing 92% of all individuals and 74% of all natural populations).Thus, the number of pseudo-absences randomly generated (in the area with calcareous screes, see methods), was 28,209 as well.
The different steps of model calibration and validation are described in Table S2 (Supplementary Materials).The final validated model chosen to predict P. occidentale habitat suitability was a GAM (which performed slightly better than the GLM during the validation procedure) including the following four environmental variables: SLOPE, TEMP (mean annual temperature), NORTHING (north-south aspect of the slope), PREC_Win (mean winter precipitation).It explained 58.64% of the total variance.The mean AUC of the validation datasets during the spatial cross-validation was 0.9.The four included variables were significant (p-value < 0.001).The best threshold calculated from the ROC curve was 0.55 (Figure S1, Table S5, Supplementary Materials).

Predicted Habitat Suitability and Ecological Preferences
According to our model, in the studied area (2335 km 2 ), the habitat can be considered as suitable (predicted value > 0.55) for P. occidentale on 33 km 2 (1.4%).The species effectively occurs in an area of about 1 km 2 (Figure 4).south to 1: north) as a function of elevation (meters), for all P. occidentale individuals in the studied area (grey dots).The black line shows a smooth approximation of the mean (GAM fit, with a 95% confidence interval).

Species Distribution Modelling
A total of 28,209 individuals (presences) were comprised in the area considered for the SDM (representing 92% of all individuals and 74% of all natural populations).Thus, the number of pseudoabsences randomly generated (in the area with calcareous screes, see methods), was 28,209 as well.
The different steps of model calibration and validation are described in Table S2 (Supplementary Materials).The final validated model chosen to predict P. occidentale habitat suitability was a GAM (which performed slightly better than the GLM during the validation procedure) including the following four environmental variables: SLOPE, TEMP (mean annual temperature), NORTHING (north-south aspect of the slope), PREC_Win (mean winter precipitation).It explained 58.64% of the total variance.The mean AUC of the validation datasets during the spatial cross-validation was 0.9.The four included variables were significant (p-value < 0.001).The best threshold calculated from the ROC curve was 0.55 (Figure S1, Table S5, Supplementary Materials).

Predicted Habitat Suitability and Ecological Preferences
According to our model, in the studied area (2335 km 2 ), the habitat can be considered as suitable (predicted value > 0.55) for P. occidentale on 33 km 2 (1.4%).The species effectively occurs in an area of about 1 km 2 (Figure 4).The most important environmental variables that influence the distribution of P. occidentale are those included in the final model, namely the slope (SLOPE), the mean annual temperature (TEMP) The most important environmental variables that influence the distribution of P. occidentale are those included in the final model, namely the slope (SLOPE), the mean annual temperature (TEMP) and the slope orientation (NORTHNESS).More surprisingly, the mean precipitation during the winter (PREC_Win) also contributes to explaining the species distribution.The response curve of these four variables can be found in Figure S5 (Supplementary Materials).The ecological preferences of P. occidentale concerning some selected environmental variables are shown in Figure 5.
A small percentage of P. occidentale individuals were growing on south-facing slopes, however only above 1900 m a.s.l.Below 1900 m a.s.l.all individuals were only growing on north facing slopes.Above this elevation, the slope aspect becomes less important (Figure 3).
The model was used to predict habitat suitability of P. occidentale under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5, see Methods), for the end of the century (year 2085, Table 2, Figure 6).The suitable area will be drastically reduced, both for all populations of P. occidentale in the study area and for all climatic scenarios.For the entire studied area reduction to 0-30% is predicted, of the current P. occidentale populations.Under the most severe scenario (RCP 8.5), the species risks complete extinction, and this is not only for small isolated populations but also for all five individual-rich and large-surface populations of Spillgerte in the canton of Bern (Table 2, Figure 6).and the slope orientation (NORTHNESS).More surprisingly, the mean precipitation during the winter (PREC_Win) also contributes to explaining the species distribution.The response curve of these four variables can be found in Figure S5 (Supplementary Materials).The ecological preferences of P. occidentale concerning some selected environmental variables are shown in Figure 5. ELEVATION and IRRAD (total irradiance, see Table 1) are shown but were not included in the final model.
A small percentage of P. occidentale individuals were growing on south-facing slopes, however only above 1900 m a.s.l.Below 1900 m a.s.l.all individuals were only growing on north facing slopes.Above this elevation, the slope aspect becomes less important (Figure 3).
The model was used to predict habitat suitability of P. occidentale under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5, see Methods), for the end of the century (year 2085, Table 2, Figure 6).The suitable area will be drastically reduced, both for all populations of P. occidentale in the study area and for all climatic scenarios.For the entire studied area reduction to 0-30% is predicted, of the current P. occidentale populations.Under the most severe scenario (RCP 8.5), the species risks complete extinction, and this is not only for small isolated populations but also for all five individual-rich and large-surface populations of Spillgerte in the canton of Bern (Table 2, Figure 6).ELEVATION and IRRAD (total irradiance, see Table 1) are shown but were not included in the final model.

Discussion
Our study demonstrates that P. occidentale is a highly specialized, cold-adapted alpine plant.It grows in natural sites mainly between 1900 and 2100 m of altitude, and all populations and individuals that grow under 1900 m a.s.l. are strictly confined to north-facing screes.Although some populations (e.g., in the canton of Bern) show locally relatively high number of individuals, the present distribution pattern and the extremely small effectively occupied area, left no doubt: P. occidentale likely reaches today the strongest reduction of its distribution range in the evolutionary history of this taxon.
Papaver occidentale, along with all members of the P. alpinum complex, belongs to the section Meconella, and shows a large number of morphological and ecological similarities with other members of this taxonomic group [30].The Meconella poppies (ca.50 species) are an evolutionarily

Discussion
Our study demonstrates that P. occidentale is a highly specialized, cold-adapted alpine plant.It grows in natural sites mainly between 1900 and 2100 m of altitude, and all populations and individuals that grow under 1900 m a.s.l. are strictly confined to north-facing screes.Although some populations (e.g., in the canton of Bern) show locally relatively high number of individuals, the present distribution pattern and the extremely small effectively occupied area, left no doubt: P. occidentale likely reaches today the strongest reduction of its distribution range in the evolutionary history of this taxon.
Papaver occidentale, along with all members of the P. alpinum complex, belongs to the section Meconella, and shows a large number of morphological and ecological similarities with other members of this taxonomic group [30].The Meconella poppies (ca.50 species) are an evolutionarily young and taxonomically poorly understood branch, mainly due to a Pleistocene or perhaps even later speciation of the majority of its members [77,78].All species of this section are specialized arctic-alpine poppies, adapted to survive in the harshest (edaphically and climatically) ecosystems inhabited by vascular plants [79].However, the P. alpinum group is the most basal in the sect.Meconella [26,27] and some authors proposed that it colonized the Alps rather from the east than from the north [31].
Nevertheless, Meconella poppies were among the very first immigrants after the Last Glacial Maximum in the high Arctic [80].Notably, P. radicatum was probably among the earliest colonizers of deglaciated parts of northern Norway between the 13,000 and 14,000 years BP [28].Papaver radicatum is also today nearly the only vascular plant (along with Saxifraga oppositifolia) inhabiting the northernmost piece of land in the Northern Hemisphere (islet of Kaffeklubben, northern Greenland) [79].Another species of the Meconella section, P. dahlianum, forms floristically poor communities on Svalbard and occupies similar habitats to P. occidentale (scree-covered mountain slopes), but also gravel beds and moraines.The species is also among the very first immigrants growing on moraines left behind the retracting glaciers on Svalbard [81].Furthermore, P. dahlianum makes part of a restricted number of vascular plants capable of growing, flowering and yielding seeds at the average daily temperature slightly over 0 • C [82], and reaching on Svalbard altitude of 940 m a.s.l.[83].
In the northeastern and southeastern Alps, individual plants and even large populations of P. alpinum s.l. were observed in gravel river beds and roadsides at low altitudes [31].Therefore, it is highly probable, that during the last glacial period, ancestors and/or relatives of P. occidentale were frequent in the vast alluvial terraces surrounding the central glaciated part of the Alps.They then followed the retracting glaciers and survived only in very few isolated and extremely small areas [30].However, recent molecular data indicate that both local survival and immigration from edge alpine refugia played an important role during the recolonization of the Western Prealps by P. occidentale, leading to a complex genetic pattern [34].Papaver occidentale is clearly in a refugial situation today.

Current Geographic Range
According to SDM, the habitats that can be considered as suitable for P. occidentale in the studied area reach about 33 km 2 .However, the effective area where the species occurs does not exceed 1 km 2 .There are likely two main reasons why the realized niche is much more restricted than the fundamental niche.Firstly, the history of the recolonization of the Alps by P. occidentale after the last glacial maximum is complex.It probably followed the retracting glaciers in some valleys on the moraines, but not in all of them.With the Holocene climate warming, the species was blocked in these valleys and was unable to reach other parallel valleys because of its limited dispersion potential.It moreover probably became extinct in some areas.The high genetic isolation of some geographically close populations, could support such a scenario or a possible local in situ survival without postglacial expansion [1].
The second reason is the method of the delimitation of calcareous screes, which could lead to a certain bias.Some delimited screes on orthophotos were included in the studied area, but might not be suitable for the species (scree already too stabilized, presence of competitive vegetation, etc.).This led to an over-optimistic suitable area, but the percentage of reduction of those suitable areas under future climate change scenarios is proportional, and thus reliable.

No Future for P. occidentale?
Thuiller et al. [84] estimated that up to 48.5% of mountain plant species in the European Alps are at risk of extinction due to the accelerated climate change.Papaver occidentale belongs without doubts to the group of alpine plants, which are the most exposed to the ongoing global warming.According to our results, the remaining suitable area will be drastically reduced, even by the weakest global warming scenario.By the RCP 8.5 climate change scenario, P. occidentale will face a complete extinction by 2085 across its whole distribution area.The severe future range contraction predicted for P. occidentale reflects similar results for other alpine endemic species, e.g., Berardia subacaulis [85] or Lilium pomponium [86].It was often concluded, that global warming may result in a significant upward shift in species optimum elevation [87], giving the threatened species an additional survival option.However, marginal Western Prealps have relatively low summit heights, making their alpine plant populations particularly vulnerable to extinction [25,88].Additionally, P. occidentale is confined to calcareous scree habitats, which will not enlarge their surface until the end of the 21st century.Although exceptionally at highest elevations, isolated individual plants can be found growing in weathered crevices of calcareous rocks, an altitudinal shift does not appear to be a viable option for this taxon.
Several studies demonstrated that alpine plants' response to climate warming might be much slower than predicted using standard species distribution models [89,90].Additionally, even highly specialized and/or endemic plants can defy unfavorable climatic conditions and persist for millennia in so-called topo-climatic microrefugia, as was demonstrated for Saxifraga florulenta in the Alps [91] or a relict lineage of Ranunculus glacialis in the Carpathians [16].In extreme cases, exiguous spots can ensure long-term survival of isolated relict populations [92].Our study, however, is a rare example of SDM specifically developed for plants inhabiting calcareous screes, enhancing thus its precision and reliability.Plants inhabiting calcareous screes will probably react faster than other plant communities (e.g., subalpine grasslands, [93]).

Conservation Implications
The entire native distribution area of P. occidentale consists of only 19 populations grouped in six isolated regions, with the total number of individuals estimated at slightly more than 30,000.Surprisingly, the region of Spillgerte (comprising all 5 populations of the canton of Bern), harbors more than 72%, and thus a huge portion of the total number of individuals.Such domination in this region is an important discovery of the present study, since the only reports on P. occidentale available to date, attributed regions of Fribourg and Vaud as the main distribution centers of this taxon [38].Furthermore, the five populations of the canton of Bern represent a unique genetic composition, being highly differentiated in comparison with all the other populations [34].These results should be therefore taken into consideration for prioritizing future conservation efforts.The introduced plants in the Jura Mountains grow in a habitat that is very different from the natural ones [38,94,95].The very low number of individuals in the introduced region, detected in our survey, confirms that P. occidentale has very specific habitat requirements and cannot prosper in any other habitat.Thus, no conservation efforts should be undertaken for these non-natural populations.
Despite its quite narrow distribution and its very specialized environmental requirements, the P. occidentale populations are not threatened by any direct human activity in the short term.The steep and mobile calcareous screes are habitats that are not appropriate for any human activity.Merely a few hiking trails get close to populations [38].Wild animals do not put the species under extinction risk either.Ibexes and chamois have been observed in the vicinity but only a few individuals have been reported as grazed [39].Thus, according to our study, global warming will be the main (if not the only) factor threatening the last remaining populations of P. occidentale in the near future.
Our results confront local nature conservation authorities with major challenges.The long-term in situ conservation of P. occidentale, and probably all other taxa of the P. alpinum complex, as well as other glacial relicts in the Alps, is unlikely to be achieved without slowing global climate change.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1424-2818/12/9/346/s1,Table S1.Populations of P. occidentale and estimated population sizes (number of individuals).Table S2.Additional information related to SDM.Table S3.Information on climate change scenarios data.Table S4.R script.Table S5.Cross-validation of the final validated GAM: AUC and MSR of the validated model, with different cross-validation techniques.Table S6.Full-field counting data.Figure S1.Series of maps representing environmental variables on the studied area.The right panels show 3 zoomed areas where P. occidentale occurs (see Figure 4 in the main manuscript).Figure S2.Paired correlation plot with all environmental variables.For the meaning of abbreviations, see Table 1 in the main manuscript.Figure S3.Paired correlation plot with the 8 selected (uncorrelated) environmental variables.For the meaning of abbreviations, see Table 1 in the main manuscript.

Figure 1 .
Figure 1.Papaver occidentale.(a) individual plant forming numerous flowers; (b) close up of the flower; (c) typical population growing in calcareous scree (Spillgerte, canton of Bern); (d) typical habitat, north-facing screes and summits (Spillgerte, canton of Bern); (e) each tuft is formed by one individual, anchored in the mobile scree with one large and long pivotal root.Plants reproduce exclusively sexually (no vegetative reproduction) and reach the maximum age of ca. 8 years.

Figure 1 .
Figure 1.Papaver occidentale.(a) individual plant forming numerous flowers; (b) close up of the flower; (c) typical population growing in calcareous scree (Spillgerte, canton of Bern); (d) typical habitat, north-facing screes and summits (Spillgerte, canton of Bern); (e) each tuft is formed by one individual, anchored in the mobile scree with one large and long pivotal root.Plants reproduce exclusively sexually (no vegetative reproduction) and reach the maximum age of ca. 8 years.

Figure 2 .
Figure 2. Location of the 21 populations of Papaver occidentale, with their corresponding name (acronyms) and number of individuals (size of dots reflect the population size).The shaded zone indicates the area considered for Species Distribution Modelling (SDM).

Figure 2 .
Figure 2. Location of the 21 populations of Papaver occidentale, with their corresponding name (acronyms) and number of individuals (size of dots reflect the population size).The shaded zone indicates the area considered for Species Distribution Modelling (SDM).

Figure 3 .
Figure 3. (a) Radial histogram indicating the number of Papaver occidentale individuals according to the slope aspect.The small surrounding ticks indicate single individuals and the black solid line gives a smooth approximation, calculated with a GAM.(b) Northness (N-S slope aspect, values from −1: south to 1: north) as a function of elevation (meters), for all P. occidentale individuals in the studied area (grey dots).The black line shows a smooth approximation of the mean (GAM fit, with a 95% confidence interval).

Figure 3 .
Figure 3. (a) Radial histogram indicating the number of Papaver occidentale individuals according to the slope aspect.The small surrounding ticks indicate single individuals and the black solid line gives a smooth approximation, calculated with a GAM.(b) Northness (N-S slope aspect, values from −1:south to 1: north) as a function of elevation (meters), for all P. occidentale individuals in the studied area (grey dots).The black line shows a smooth approximation of the mean (GAM fit, with a 95% confidence interval).

Figure 4 .
Figure 4. Modelled current habitat suitability of Papaver occidentale in the Western Swiss Prealps.White areas are less suitable, red and purple are the most suitable.Three areas are zoomed on the right panel (a, Vanil Noir region; b, Spillgerte region; c, Gummfluh region).The black dots indicate the locations of P. occidentale (presences).

Figure 4 .
Figure 4. Modelled current habitat suitability of Papaver occidentale in the Western Swiss Prealps.White areas are less suitable, red and purple are the most suitable.Three areas are zoomed on the right panel (a, Vanil Noir region; b, Spillgerte region; c, Gummfluh region).The black dots indicate the locations of P. occidentale (presences).

Figure 5 .
Figure 5. Ecological preferences of Papaver occidentale concerning some selected environmental variables.P. occidentale individuals are represented by red dots and background (pseudo-absences) are represented by grey dots.The box plots summarize the preference of P. occidentale concerning each variable.ELEVATION and IRRAD (total irradiance, seeTable 1) are shown but were not included in the final model.

Figure 5 .
Figure 5. Ecological preferences of Papaver occidentale concerning some selected environmental variables.P. occidentale individuals are represented by red dots and background (pseudo-absences) are represented by grey dots.The box plots summarize the preference of P. occidentale concerning each variable.ELEVATION and IRRAD (total irradiance, see Table1) are shown but were not included in the final model.

Figure 6 .
Figure 6.Modelled future habitat suitability of Papaver occidentale in the Western Swiss Prealps, under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), at the end of the century (year 2085).White areas are less suitable, red and purple are the most suitable.Two areas are zoomed on the right panel (a, Vanil Noir region; b, Spillgerte region).Black dots indicate the location of present P. occidentale populations.

Figure 6 .
Figure 6.Modelled future habitat suitability of Papaver occidentale in the Western Swiss Prealps, under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), at the end of the century (year 2085).White areas are less suitable, red and purple are the most suitable.Two areas are zoomed on the right panel (a, Vanil Noir region; b, Spillgerte region).Black dots indicate the location of present P. occidentale populations.
Figure S4.Receiver operating characteristic (ROC) curves for different cross-validation techniques (red line, full model ROC curve with all the data, grey lines, ROC curves of validation datasets).
Figure S5.Response curves for the final validated General Additive Model (GAM).Author Contributions: G.K. and Y.F.conceived the idea.L.P., S.B., B.C. and E.G. conducted the fieldwork.Y.F.conducted statistical and modelling analyses.G.K., M.R. and C.P. edited the manuscript.All authors have read and agreed to the published version of the manuscript.Funding: M.R. is financially supported by the statutory funds of the W. Szafer Institute of Botany, Polish Academy of Sciences.

Table 1 .
Environmental variables selected for the Species Distribution Model (SDM) of Papaver occidentale.

Table 2 .
Remaining suitable area (predicted value 0.55) under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), at the end of the century for Papaver occidentale.The remaining future suitable areas are given as a percentage of the actual suitable predicted area.

Change Scenario Remaining Suitable Area (Entire Studied Area)
Diversity 2020, 12, x FOR PEER REVIEW 11 of 19

Table 2 .
Remaining suitable area (predicted value > 0.55) under different climate change scenarios (RCP 2.6, RCP 4.5 and RCP 8.5), at the end of the century for Papaver occidentale.The remaining future suitable areas are given as a percentage of the actual suitable predicted area.