Plant Diversity Patterns and Conservation Implications under Climate-Change Scenarios in the Mediterranean: The Case of Crete (Aegean, Greece)

Climate change poses a great challenge for biodiversity conservation. Several studies exist regarding climate change’s impacts on European plants, yet none has investigated how climate change will affect the extinction risk of the entire endemic flora of an island biodiversity hotspot, with intense human disturbance. Our aim is to assess climate change’s impacts on the biodiversity patterns of the endemic plants of Crete (S Aegean) and provide a case-study upon which a climate-smart conservation planning strategy might be set. We employed a variety of macroecological analyses and estimated the current and future biodiversity, conservation and extinction hotspots in Crete. We evaluated the effectiveness of climatic refugia and the Natura 2000 network of protected areas (PAs) for protecting the most vulnerable species and identified the taxa of conservation priority based on the Evolutionary Distinct and Globally Endangered (EDGE) index. The results revealed that high altitude areas of Cretan mountains constitute biodiversity hotspots and areas of high conservation and evolutionary value. Due to the “escalator to extinction” phenomenon, these areas are projected to become diversity “death-zones” and should thus be prioritised. Conservation efforts should be targeted at areas with overlaps among PAs and climatic refugia, characterised by high diversity and EDGE scores. This conservation-prioritisation planning will allow the preservation of evolutionary heritage, trait diversity and future ecosystem services for human well-being and acts as a pilot for similar regions worldwide.


Introduction
Earth has entered a new geological epoch, the Anthropocene [1], characterised by human-induced temperature increase, thus placing immense impacts upon natural, atmospheric and hydrological processes. Consequently, global ecosystem health is severely challenged, with shifts in biotic composition being among the major threats of human activity worldwide. Species unable to survive to such a novel adaptive matrix are on an inevitable path to extinction [2], often causing a deadlock in nature management and protection practices. Extinction threats might be more prominent in endemic species suffering from habitat loss and fragmented populations [3], since they have narrower geographical ranges and environmental niches. Climate change is projected to rapidly change community structure in the near future, likely surpassing the land-use effects by the 2070s [4] and significantly altering terrestrial biodiversity [3]. In a global analysis of future climate change impacts, ca. 60% of the plant species examined were predicted to lose more than half of their current climatic range in the coming decades, possibly leading to a substantial global biodiversity decrease and ecosystem function degradation [3]. Even though plants might be more resistant to extinction compared to animals (extinction debt- [5,6]), climate change could elevate the estimated extinction rates [7]. Nearly 20-30% of species would face high extinction risks if global temperature rise beyond 2-3 • C above pre-industrial levels [8]. In order to avoid a 1.5 • C rise in global temperatures until 2052, unprecedented changes should be made on a global scale [9]. In this context, it is important to follow a "climate-smart" conservation planning strategy to efficiently conserve as many species/communities as possible [10]. This strategy relies on the identification of climatic refugia, i.e., relatively climatically stable regions for many species of high conservation value that are characterised by high endemic and genetic diversity [11].
The Mediterranean Basin is a biogeographically complex area, with high levels of endemism and diversification, due to its high topographic and climatic heterogeneity [12]. It constitutes the second largest terrestrial biodiversity hotspot in the world [13] with nearly 12,500 plant species endemic to the region, most of which have an extremely narrow geographical range [12]. The Mediterranean islands have high overall and endemic species richness [12]. Plant endemism varies between 9 and 18% on the largest islands and reaches up to 40% in the high-altitude zones of their mountain ranges [12]. Several of these islands constitute biodiversity hotspots and climatic refugia [14]. Crete (S. Greece) is rendered as the hottest endemism hotspot of the Mediterranean Basin, owing to its high species number and endemism levels and long-term geographical isolation and climatic stability, as well as to its high environmental and topographical heterogeneity [12].
The current global mean annual temperature will probably continue to rise [9] and has already forced species to migrate into higher altitudes, changed flowering periods/durations and, in some cases, led species to extinction [15]. These impacts will most likely accelerate [7] and are expected to be more severe for species occurring on mountain ranges, due to the "escalator to extinction" phenomenon [16]. The Mediterranean is a global hotspot of vulnerable species [17] and among the regions expected to experience the largest changes in climate [18], with climate change impacts being more prominent on islands and mountain summits [12]. However, the Mediterranean Basin has experienced few plant species extinctions (22 taxa-0.17% of the total; [19]). Most Mediterranean countries have their own Red List of threatened plants. Nevertheless, regarding Crete, ca. 30% (58 taxa) of the Cretan single island endemics (SIE) have been assessed and only five are considered to be facing imminent extinction [20]. Several SIE occur on the island's mountain ranges and comprise small and isolated populations, thus being prone to climate change impacts.
New light can be shed upon the mechanisms underlying species' possible extinction and their subsequent conservation needs when integrating phylogenetic diversity metrics into conservation-centred analyses [21]. Evolutionary distinct species usually have unusual phenotypes, rare ecological roles and increased functional importance [22]. As such, they have great conservation value, since they cannot be readily replaced. Thus, if evolutionary isolated species become extinct, this would result in great evolutionary loss. Therefore, any conservation plan needs to take into account the advantage of the conservation prioritisation of evolutionary distinct species, as it allows the preservation of evolutionary heritage, trait diversity and potential future services for human well-being [23].
Correlative species distribution models (SDMs) are widely used to identify which species are most vulnerable, as well as to answer how, why, where and when these species are vulnerable [24]. There is ample evidence that the expectations from correlative SDMs have matched recent population trends focusing on birds, mammals and plants, in decreasing order, and they provide appropriate results when the aim is to estimate a species' extinction risk [17]. To our knowledge, even though continent-wide studies exist regarding the expected vegetation changes and the impacts of climate change on European plants [25], none has investigated how climate change will affect the extinction risk of the entire endemic flora of an island biodiversity hotspot, such as Crete. Our aim is to assess, in an integrated manner, the impact of climate change on the biodiversity and conservation patterns of the endemic plants of Crete and to provide a case-study upon which a cost-effective and climate-smart conservation planning strategy might be set.

Study Area
We focus on Crete (Greece- Figure 1), the fifth largest island in the Mediterranean Basin and the richest island hotspot of Europe in terms of endemic plant species [12].
Crete is the largest (8836 km 2 ) and most topographically complex island in the Aegean archipelago, with more than 50 mountain peaks exceeding 2000 m a.s.l. The island is climatically compartmentalised on a W-E axis and is characterised by a sharp altitudinal gradient (0-2456 m a.s.l.). There are three main mountain massifs present on Crete (Lefka Ori, Idi and Dikti), with an entirely different climate to the warm and dry lowland plains.
Geologically, Crete comprises mostly limestones, is part of the Hellenic Arc and was formed in response to the subduction of the African plate beneath the Aegean [26].
Diversity 2020, 12, x FOR PEER REVIEW 3 of 23 change on European plants [25], none has investigated how climate change will affect the extinction risk of the entire endemic flora of an island biodiversity hotspot, such as Crete. Our aim is to assess, in an integrated manner, the impact of climate change on the biodiversity and conservation patterns of the endemic plants of Crete and to provide a case-study upon which a cost-effective and climatesmart conservation planning strategy might be set.

Study Area
We focus on Crete (Greece- Figure 1), the fifth largest island in the Mediterranean Basin and the richest island hotspot of Europe in terms of endemic plant species [12].
Crete is the largest (8836 km 2 ) and most topographically complex island in the Aegean archipelago, with more than 50 mountain peaks exceeding 2000 m a.s.l. The island is climatically compartmentalised on a W-E axis and is characterised by a sharp altitudinal gradient (0-2456 m a.s.l.). There are three main mountain massifs present on Crete (Lefka Ori, Idi and Dikti), with an entirely different climate to the warm and dry lowland plains.
Geologically, Crete comprises mostly limestones, is part of the Hellenic Arc and was formed in response to the subduction of the African plate beneath the Aegean [26].
Its palaeogeographical history is intricate and defined by two main geological events that created important dispersal barriers (for a thorough review of the Aegean's palaeogeographical history, see [27]): (i) the isolation of Crete from Karpathos' and Cyclades' island groups 12 Mya and (ii) the isolation of Crete from Peloponnese after the Messinian salinity crisis [28].  Its palaeogeographical history is intricate and defined by two main geological events that created important dispersal barriers (for a thorough review of the Aegean's palaeogeographical history, Diversity 2020, 12, 270 4 of 22 see [27]): (i) the isolation of Crete from Karpathos' and Cyclades' island groups 12 Mya and (ii) the isolation of Crete from Peloponnese after the Messinian salinity crisis [28].

Environmental Data
Current and future climatic data were obtained from the WorldClim database [29] at a 30 sec resolution. We constructed 16 additional climatic variables at the same resolution via the "envirem" 1.1 [30] R package based on the 19 bioclimatic variables from WorldClim for current and future climate conditions. We selected three Global Circulation Models (GCMs) that are rendered more suitable and realistic for the study area's future climate based on [31] and two different IPCC scenarios from the Representative Concentration Pathways (RCP) family: RCP2.6 (mild scenario) and RCP8.5 (severe scenario). Seven soil variables providing predicted values for the surface soil layer at varying depths were obtained from the SoilGrids database [32]. We obtained elevation data from the CGIAR-CSI data-portal [33] and then aggregated and resampled them using the "raster" 2.6.7 R package [34] to match the resolution of the other environmental variables.
From this initial set of 43 predictors, only seven (isothermality, precipitation of the wettest month, precipitation seasonality, annual potential evapotranspiration, the Thornthwaite aridity index, continentality and soil pH) were not highly correlated (Spearman rank correlation < 0.7 and VIF < 5- [35]) and used in all the subsequent analyses. Multicollinearity assessment was performed with the "usdm" 1.1.18 R package [36].

Climatic Refugia
We located sites in Crete that have been climatically stable (Sensu, [37]) for the past 4 My, using paleoclimatic data obtained from Paleoclim [38] and Oscillayers [39], the framework of [37] and the "climateStability" 0.1.1 R package [40]. Climatic refugia as designated here refer to macro-refugia according to [41] and have a climate stability index equal to or higher than 0.5 (the climatic stability index ranges between 0-1 according to [37] and is derived as described above).

Species Occurrence Data
Crete hosts 2240 native plant species, 395 of which are Greek endemics and 183, single island endemics (SIE- [42][43][44]). Based on the most extensive and detailed database (Flora Hellenica Database, Strid (ongoing)) of plants occurring in Greece (~1.2 M occurrences), we compiled a dataset of all the SIE present in Crete (8773 occurrences). All the SIE were cross-checked for synonyms, following the nomenclature proposed by the Vascular Plants Checklist of Greece and the Atlas of the Aegean flora [42][43][44]. Our final dataset included 172 SIE, since we modelled only those SIE with at least three locations, following [45].
To avoid pseudoreplication and associated spatial sampling biases, we selected occurrences with a minimum distance of 1 km from each other. We removed any records not conforming to this criterion to reduce sampling bias, while keeping as much useful information as possible, according to the procedure described below. By doing so, we spatially disaggregated the occurrences to ensure that the SDMs did not over-represent the environmental conditions associated with over-sampled regions, which hinder the interpretation and application of models [46,47]. This data cleaning and organising procedure followed the protocols as set out in [48], and we used the "biogeo" 1.0 [47] and "spThin" 0.1.0 [46] R packages. We also evaluated whether any geographical sampling bias existed in our species occurrence data by comparing the statistical distance distribution observed in our dataset to a simulated distribution expected under random sampling via the "sampbias" 0.1.1 [49] R package.

Phylogenetic Tree
We used a "supertree" approach to generate our phylogenetic tree, as most SIE do not have available molecular data. A phylogenetic tree was generated for all the SIE included in our analyses, based on the recently published phylogeny of seed plants by [50,51]. We used the GBOTB extended tree (i.e., GenBank taxa with a backbone provided by Open Tree of Life version 9.1), which contains 74,531 taxa and is the largest dated mega-tree for vascular plants [51]. We appended any SIE that were present in Crete but missing from the phylogeny, by adding them next to a randomly selected congener, following [52,53]. We did not add the missing taxa as polytomies to their respective genera, as this approach adds substantial bias to any ensuing analyses [54]. We then pruned the phylogeny to keep only the 172 SIE. We conducted a sensitivity analysis regarding the phylogenetic tree generation and constructed an additional phylogenetic tree following the framework proposed by [51]. Both trees are available as supplementary files. All the subsequent phylogenetic analyses were computed for both trees.
Evolutionary distinctiveness (ED) was then estimated in the phylogenetic tree using the "picante" 1.6.2 [55] R package. Evolutionary Distinct and Globally Endangered (EDGE) scores were calculated according to [56]: where ED is a species' ED value and GE is its weighted International Union for Conservation of Nature (IUCN) threat category (LC = 0; NT = 1; VU = 2; EN = 3; CR = 4) on a log scale. We also estimated for each grid cell the phylogenetic alpha diversity (PD-Sensu, [57]) of the species inhabiting each of the grid cells with the "picante" 1.6-2 R package [54] and the standardised effect size scores (SES) with the "PhyloMeasures" 2.1 R package [58]. We tested for non-random patterns in PD by estimating their SES scores as where X obs is the observed score within each grid cell, and mean (x null ) and standard deviation (X null ) are the mean and standard deviation of a null distribution of scores generated by shuffling the taxa labels of the grid cell-by-species matrix 999 times. We assessed the statistical significance of the SES scores by calculating two-tailed p-values (quantiles) as: p-values = rank obs runs + 1 where rank obs is the rank of the observed scores compared with those of their null distributions, and runs is the number of randomisations [55]. SES scores with p < 0.05 and p > 0.95 indicate a significantly lower and higher than expected for a given PD value, respectively. Positive and negative SES values indicate phylogenetic overdispersion or clustering, respectively. The greater sensitivity of SES PD to a more terminal structure makes it better suited to exploring assembly processes working at finer temporal and spatial scales [59].

Model Parameterisation and Evaluation
We modelled the realised climatic niche of each species by combining the available occurrence data with current environmental predictors with the "biomod" 3.3.7 R package [60]. We used three different modelling algorithms for species with more than ten occurrences-Random Forest (RF), Classification Tree Analysis (CTA) and Multiple Adaptive Regression Splines (MARS)-in an ensemble modelling scheme, as ensemble forecasting integrates the results of multiple SDM algorithms into a single geographical projection for each time period, reducing the uncertainties associated with the use of a single model algorithm [61,62]. Since these algorithms require presence/absence (PA) data, we generated PAs following the recommendations of [63], and pseudo-absences were generated at a minimum distance of 19.8 km from presence locations to reduce the probability of false absences. We chose that minimum distance due to the median autocorrelation of 19.7 km among the non-collinear environmental variables, which we computed with the "blockCV" 1.0.0 [64] R package. PA generation and model calibration were repeated 100 times per species to ensure that the selected pseudo-absences did not bias the final predictions. Regarding species with fewer than ten occurrences, we followed the ensemble of small models (ESM) framework [65], which is suitable for modelling rare species [65][66][67], using the RF algorithm (single-technique ESMs perform equally as good as double ensembles- [66]), which is robust to overfitting [68]. We calibrated ESMs by fitting numerous bivariate models (987 bivariate models), which were then averaged into an ensemble model using weights based on model performances. For all the models, the weighted sum of the presences was equal to that of the PAs. The models' predictive performance was evaluated via the True Skill Statistic (TSS; [69]) based on a repeated (10 times) split-sampling approach in which the models were calibrated with 80% of the data and evaluated over the remaining 20%. We used null model significance testing [70] to evaluate the performance of all the models and estimated the probability that each model performed better than 100 null models. All the models were found to outperform the null expectation at p < 0.001.

Model Projections
Calibrated models were used to project the suitable area for each species in Crete under current and future conditions through an ensemble forecast approach [61]. The contribution of each model to the ensemble forecast was weighted according to its TSS score. Models with a TSS score <0.8 were excluded from building projections, so as to avoid working with poorly calibrated models. The resulting habitat suitability maps were converted into binary maps and then compared to the binary maps obtained for each GCM, RCP and time-period for each SIE.

Area Range Change
To assess whether the 172 SIE would experience range contraction or expansion under future conditions, we used the "biomod" 3.3.7 R package [60]. Species were not assumed to have unlimited dispersal capability, since this assumption could be overoptimistic. Species were assumed to have very limited dispersal ability, based on the results from the "BIOMOD_RangeSize" function of the "biomod" 3.3.7 R package [60].

Hotspots
For each climate scenario, we stacked the final binary maps of all 172 SIE. From this superimposition, we estimated for each 30 s grid-cell the number of taxa that would find suitable environmental conditions there. We defined potential SIE hotspots as the 20% of cells that provided suitable environmental conditions for the highest number of taxa.
To identify sites where the predictions of hotspot locations differed significantly between current and future environmental conditions, we applied the methodology proposed by [71], using the functions from the "SDMTools" 1.1.221 R package [72].

IUCN Measures
We followed the Preliminary Automated Conservation Assessment (PACA) framework [73] and calculated the standard IUCN measures Extent of Occurrence (EOO) and Area of Occupancy (AOO), and we assigned each SIE to a preliminary IUCN threat category according to Criterion A and B under current and future conditions using the "ConR" 1.1.1 R package [74]. Under the same framework, we assessed the SIE under Criterion A, using the R code provided by [73]. To calculate Criterion A with this approach, we used the CORINE land cover (CLC) data v.20 [75]. CLC layers 1-2 [apart from 223 (Olive groves) and 243-244 (land principally occupied by agriculture, with significant areas of natural vegetation and agro-forestry areas, respectively)] are directly linked to the main threats to SIE. After species were assigned to IUCN and PACA categories, we estimated the total number of taxa recorded and the proportion of taxa assessed under each IUCN and PACA category, under Criteria A and B separately and by combining both criteria, i.e., a taxon would, for example, be categorised as Critically Endangered (CR) if it was assessed as CR by at least one of the two criteria [73].

Niche Breadth
Levins' inverse concentration measure of niche breadth [76] was computed for each taxon using the "ENMTools" 0.2 R package [77]. The niche breadth values ranged from 0 (specialists) to 1 (generalists), thus being comparable among taxa [76]. The difference in the niche breadth between the different IUCN threat categories was investigated via a Kruskal-Wallis non-parametric test (KWA), as was the difference in the area range change between SIE with broad or narrow environmental niches (above and below the first quartiles of the niche breadth values of the dataset, respectively).

GDM Analysis
We used Generalised Dissimilarity Modelling (GDM- [78]) to model the pairwise plant community compositional dissimilarity (Sorensen's dissimilarity) within map grid cells across Crete as a response to environmental and spatial variables and aid us in locating the areas with the highest future floristic turnover. We used the same environmental variables as in the SDM analyses, with the significance of all the variables assessed through a Monte Carlo permutation test (1000 repetitions). We quantified the magnitude of turnover for each environmental variable and the relative importance of that variable in driving species turnover, according to [79,80]. To evaluate the unique contributions of environment and space in explaining species turnover, we partitioned the deviance resulting from sets of three GDMs that used environmental variables, geographical distance (a proxy for dispersal limitation and one of the primary drivers of beta diversity in island biotas [81]) or both as predictor variables [82]. We projected models trained on current environmental conditions onto the future climatic models and scenarios described above. The resulting predicted species turnover is a function of the magnitude of climate change at a given location and at a given time, the starting position along climatic gradients and the modelled rate of turnover at that position [83]. All the GDM analyses were performed with the "gdm" 1.3.7 R package [82].

Current and Future Spatial EDGE Patterns
We derived the mean EDGE values for the species present in each grid cell under current and future climate scenarios and calculated their difference by stacking the presences from the individual species models. We used a grid cell resolution of 1 km to match the resolution of the predictor variables. A grid cell was considered occupied if it overlapped any part of a species' projected distribution [84]. Negative and positive values indicate areas that were predicted to become extinction hotspots and coldspots, respectively. The ∆EDGE index can be considered as a proxy of conservation prioritization due to evolutionary distinct and highly threatened species becoming extinct in a given grid cell. Negative values indicate that in a given grid cell several species that are assessed for instance, as Critically Endangered and have high ED value, might be driven to extinction due to climate-change. Thus, the future mean EDGE value for that specific grid cell will be lower than that reported for the current time-period, indicating that probably immediate conservation actions are needed. The more negative the ∆EDGE index, the more urgent the need will be for these actions to take place.

Protected Area Network and Climatic Refugia Overlap
We overlapped current and future hotspot results with the areas recognised as climatic refugia, as well as with the protected areas (PA) network retrieved from the World Database on Protected Areas (WPDA) using functions from the "wpdar" 1.0.0 [85] and the "sf" 0.8.0 [86] R packages. Protected areas exclusively related to marine protection were excluded. Based on each species' current and future EOO, we calculated the irreplaceability of each PA and climatic refugium in Crete for the current and future climate conditions. The irreplaceability index represents biotic uniqueness and quantifies the degree of overlap between each PA/climatic refugium and the range of SIE [87]. We investigated whether the degree of overlap differed between the current and future conditions via a Kruskal-Wallis non-parametric test.

Model Performance
The models for most SIE had sufficient predictive power (TSS ≥ 0.7-mean TSS: 0.94; Table S1). Precipitation-related variables had the highest contribution among the response variables for most SIE (53.5%), while temperature-related variables and soil pH were important for a smaller fraction of SIE (42.4% and 4.1%, respectively). The full details and analyses of the model performance are given in Table S1.

Area Range Change
All SIE species will experience range contraction of varying magnitudes across taxa according to the GCMs and RCPs (11.7-100.0%; Table S1), but the median range contraction was 98.3% (see Table S2 for the median values for each GCM and RCP). The median range contraction differed significantly between species with narrow and broad niches (KWA: H = 20.7, d.f. = 1, p < 0.01; median range contraction was 97.8% and 88.7%, respectively).

Hotspot Analysis
The eastern, central and high-altitude parts of Crete are more likely to lose most of the SIE occurring there (Figure 2 and Figure S1). The SIE hotspots are currently located above 1500 m a.s.l. in Crete (Figure 2), but all of them are expected to shift downwards-even below 1000 m a.s.l. under any GCM/RCP (Figure 2 and Figure S1)-and these shifts are significantly different in absolute or relative terms ( Figures S2 and S3).

IUCN Measures
For the first time, we provide a preliminary assessment of every SIE occurring in Crete according to the IUCN Criteria A and B (Figure 2; Table S1). At present, 58.7% of SIE are facing imminent extinction (Figure 2), and 87.8% are characterised as Likely Threatened (LT) under the PACA categories (Table S1). This phenomenon could greatly deteriorate, since under any GCM and RCP included in our analyses, many of these species are projected to become extinct (median EX%: 48.85%- Figure 2, Table S1). Nineteen taxa are projected to become extinct under any GCM/RCP (Table S1). At the family level, Convolvulaceae and Hyacinthaceae show the lowest and highest median EX% (70.7% and 100%, respectively, Table S3). The distribution of taxa assessed as Threatened under either the IUCN or the PACA categories under both Criteria A and B is not uniform across Crete (Figure 3). These taxa are mostly concentrated at the Cretan mountain massifs, with the highest number of them occurring in Lefka Ori (Figure 3). These high-altitude areas are predicted to become extinction hotspots under any GCM/RCP (Figure 3 and Figure S4).

Hotspot Analysis
The eastern, central and high-altitude parts of Crete are more likely to lose most of the SIE occurring there (Figure 2 and Figure S1). The SIE hotspots are currently located above 1500 m a.s.l. in Crete (Figure 2), but all of them are expected to shift downwards-even below 1000 m a.s.l. under any GCM/RCP (Figure 2 and Figure S1)-and these shifts are significantly different in absolute or relative terms (Figures S2 and S3).

Niche Breadth
The median niche breadth was 0.41 for the SIE, and 25% of the SIE were considered to have a narrow niche breadth (Table S1). The niche breadth differed significantly between all the IUCN categories (KWA: H = 56.6, d.f. = 3, p < 0.05; the median niche breadth for the IUCN categories was 0.602, 0.398, 0.314 and 0.168 for the LC/NT, VU, EN and CR categories, respectively), and species with narrow niches had higher extinction probabilities according to every GCM/RCP (Table S1).

GDM Analysis
GDM analysis helped to disentangle the relative contribution of the geographic and environmental drivers of SIE community composition across time and space. Our model explained 26.6% of the deviance in SIE composition ( Figure S5). The most important gradient for determining SIE turnover was annual potential evapotranspiration (PET), followed by geographical distance (Figure 4). Soil pH was the weakest predictor of SIE turnover. The fitted functions describing the turnover rate and magnitude along each gradient were nonlinear, with the turnover rates varying by position along the gradients (Figure 4). An acceleration zone in SIE turnover was evident for PET and short geographical distances, while a sharp compositional transition was apparent along the precipitation of the wettest month gradient (Figure 4). Under any GCM/RCP, the low-elevation area in W. Crete was predicted to exhibit a greater species turnover than any other region, followed by the high-altitude Cretan mountain ranges ( Figure 5 and Figure S6).

EDGE-Phylogenetic Diversity
The results from both phylogenetic trees were similar, so we report the results based on the phylogenetic tree that was constructed following [52,53].
Based on both Criteria A and B, the SIE have EDGE scores in the range 3.57-8.74, with most species (75%) falling in an EDGE score range of 3.57-6.39. Eleven species (Tables S1 and S4) have an EDGE score exceeding 7.0, belonging to eleven different families, and all are considered as CR taxa (Table S1).
Higher and lower PD than expected was found in 21 and 86 sites, respectively ( Figure S7). These sites did not differ significantly in terms of elevation or pH preferences (both occur above 850 m a.s.l.), but the non-significant sites were located in lower elevation areas with higher pH (KWA: H = 261.93, d.f. = 5, p < 0.05; Table S5).

Current and Future Spatial EDGE Patterns
The patterns of change regarding the EDGE index show that the high-altitude areas of Crete are currently identified as hosting assemblages of great evolutionary distinctiveness facing immediate extinction risks ( Figure S8). These areas are projected to become extinction hotspots under any

EDGE-Phylogenetic Diversity
The results from both phylogenetic trees were similar, so we report the results based on the phylogenetic tree that was constructed following [52,53].
Based on both Criteria A and B, the SIE have EDGE scores in the range 3.57-8.74, with most species (75%) falling in an EDGE score range of 3.57-6.39. Eleven species (Tables S1 and S4) have an EDGE score exceeding 7.0, belonging to eleven different families, and all are considered as CR taxa (Table S1).
Higher and lower PD than expected was found in 21 and 86 sites, respectively ( Figure S7). These sites did not differ significantly in terms of elevation or pH preferences (both occur above 850 m a.s.l.), but the non-significant sites were located in lower elevation areas with higher pH (KWA: H = 261.93, d.f. = 5, p < 0.05; Table S5).

Current and Future Spatial EDGE Patterns
The patterns of change regarding the EDGE index show that the high-altitude areas of Crete are currently identified as hosting assemblages of great evolutionary distinctiveness facing immediate extinction risks ( Figure S8). These areas are projected to become extinction hotspots under any GCM/RCP (Figure 3 and Figure S9). In the future, the mid-altitude areas are predicted to take up their place (Figure 3 and Figure S9).

Climatic Refugia and Protected Area Network Overlap
Six areas were identified as climatic refugia and are located along a W-E axis in Crete ( Figure 6). The largest and most species-rich climatic refugium is located in the wider area of Lefka Ori in W. Crete, while the smallest and most species-poor is near the Asterousia mountain range in south-central Crete ( Figure 6, Table S6). Two climatic refugia are significantly phylogenetically clustered, and a single one is overdispersed (Table S6). In the future, all these areas are predicted to continue to serve as climatic refugia, with different trends: the species-poor and species-rich refugia show trends towards phylogenetic dispersion (PD increases) and clustering (PD decreases), respectively (Table S6).
Diversity 2020, 12, x FOR PEER REVIEW 13 of 23 GCM/RCP (Figure 3 and Figure S9). In the future, the mid-altitude areas are predicted to take up their place (Figure 3 and Figure S9).

Climatic Refugia and Protected Area Network Overlap
Six areas were identified as climatic refugia and are located along a W-E axis in Crete ( Figure 6). The largest and most species-rich climatic refugium is located in the wider area of Lefka Ori in W. Crete, while the smallest and most species-poor is near the Asterousia mountain range in southcentral Crete ( Figure 6, Table S6). Two climatic refugia are significantly phylogenetically clustered, and a single one is overdispersed (Table S6). In the future, all these areas are predicted to continue to serve as climatic refugia, with different trends: the species-poor and species-rich refugia show trends towards phylogenetic dispersion (PD increases) and clustering (PD decreases), respectively (Table  S6). The overlap between the PA network and the climatic refugia ranges between 13.3 and 97.0% (Table S7). The mean irreplaceability index differs significantly between PAs and climatic refugia (KWA: H = 56.6, d.f. = 1, p < 0.01; Figure S10), as well as between PAs and any GCM/RCP (KWA: H = 103.3, d.f. = 6, p < 0.01). Regarding climatic refugia, the mean irreplaceability index is significantly different only between the current and the HadGEM2 RCP 8.5 climate projection (KWA: H = 14.6, d.f. = 6, p < 0.05; Figure S10). When taking out the effect of area, the mean irreplaceability index does not differ significantly between PAs and climatic refugia, except for the current and the HadGEM2 RCP 8.5 climate conditions (KWA: H = 2.79, d.f. = 1, p = 0.09; Figure 7).

Discussion
Climate change is projected to alter biodiversity and biogeographical patterns all over the globe [88]. The Mediterranean Basin is expected to face the largest changes in climate worldwide [18], with The overlap between the PA network and the climatic refugia ranges between 13.3 and 97.0% (Table S7). The mean irreplaceability index differs significantly between PAs and climatic refugia (KWA: H = 56.6, d.f. = 1, p < 0.01; Figure S10), as well as between PAs and any GCM/RCP (KWA: H = 103.3, d.f. = 6, p < 0.01). Regarding climatic refugia, the mean irreplaceability index is significantly different only between the current and the HadGEM2 RCP 8.5 climate projection (KWA: H = 14.6, d.f. = 6, p < 0.05; Figure S10). When taking out the effect of area, the mean irreplaceability index does not differ significantly between PAs and climatic refugia, except for the current and the HadGEM2 RCP 8.5 climate conditions (KWA: H = 2.79, d.f. = 1, p = 0.09; Figure 7).

Discussion
Climate change is projected to alter biodiversity and biogeographical patterns all over the globe [88]. The Mediterranean Basin is expected to face the largest changes in climate worldwide [18], with these impacts being more prominent on islands and mountain summits [12]. Here, we used the hottest endemic Mediterranean hotspot, Crete, as a case-study for conservation-scenario building and decision making, based on an integrated assessment of climate change impacts on biodiversity and conservation patterns. Our results highlight an augmented extinction risk for the majority of the SIE, while pinpointing areas of high conservation and evolutionary value; they should alert conservation practices and management to take measures for biodiversity loss restraint and the halting of further deterioration.

Diversity Hotspots
The high-altitude areas of Crete are identified as species richness hotspots, with Lefka Ori, the western mountain massif, hosting the most SIE ( Figure 2). These areas are predicted to experience a sharp decline in species richness under any GCM/RCP and will no longer constitute biodiversity hotspots, as a result of the "escalator to extinction" phenomenon [16]. The high-altitude Cretan mountain ranges will experience great floristic turnover under any GCM/RCP ( Figure 5 and Figure S6). The most important gradient for determining SIE turnover is annual PET, followed by geographical distance (Figure 4 and Figure S5). Niche-based processes were found to exert some influence on the distribution of species occurring in Crete [89]. Due to upward species' range shifts and the extinction of high-altitude SIE, mid-altitude areas will host an increasing number of species, thus intensifying the observed mid-domain effect of the SIE [90].

Conservation Assessment
By incorporating evolutionary history into biodiversity and conservation analyses, we can decide whether or not PAs are indeed sheltering different aspects of biodiversity [91] and improve conservation management in the unprecedented biodiversity decline setting that is currently underway [92]. By doing so, we were able to identify new areas of high evolutionary and conservation value. In total, 107 sites had PD higher or lower than expected by chance ( Figure S7), scattered across a W-E axis. Phylogenetically overdispersed sites of high conservation importance [93] occur at high altitudes in Crete, and this may be related to a variety of reasons, such as the competitive exclusion of closely related taxa with high niche overlap, the colonisation of phylogenetically distinct lineages (Campanula/Roucela species) and the high topographical heterogeneity of the Cretan mountain massifs that harbour distinctly adapted plant lineages [94].
We evaluated the potential conservation status of SIE using the IUCN Criteria A and B. Most of the SIE were potentially threatened by extinction ( Figure 2; Table S1). These potential extinctions were not randomly distributed among evolutionary groups, when accounting for mean future area loss (C mean = 0.12; p < 0.01), a phenomenon observed in mammals as well [95]. The SIE will be highly vulnerable in the future, since up to 154 taxa are projected to become extinct ( Figure 2; Table S1) and at least 19 taxa are projected to become extinct under any GCM/RCP (Table S1). This is true irrespective of the species' niche breadth (Table S1) and will most probably also lead to severe genetic diversity decline, as has been observed in other taxa occurring in Greece, regardless of their endemicity status (e.g., Cicer graecum and Helleborus odorus subsp. Cyclophyllus- [96,97]). The high-altitude areas that now serve as diversity and conservation hotspots (due to high EDGE scores- Figures S1 and S8) are projected to become extinction hotspots (due to very low ∆EDGE scores- Figure 3 and Figure S9) and should thus be prioritised in terms of conservation efforts. Eleven taxa (Table S4) should also be given high conservation priority, since they have a very high EDGE score. However, none of them is included in the IUCN Top-50 Mediterranean Island Plants initiative [98]. Identifying conservation priorities and implementing effective actions are urgently needed and will become increasingly important, since human activities and land use are exerting unprecedented pressure on natural environments, leading to severe biodiversity changes [9,23]. In this context, locating areas with high EDGE and low ∆EDGE scores could help to prioritise areas of high conservation and evolutionary importance and high extinction risk.

Conservation and Management Implications
Incorporating climate change projections is critically important for conservation strategies. Identifying climatic refugia is an integral part of the "climate-smart" conservation planning framework, since they constitute taxonomic and genetic diversity centres [11], with macro-refugia-areas that sustain climatic suitability along broad spatiotemporal gradients [41]-being generally more easily captured by GCMs. In the Mediterranean basin, where land-use change and the occurrence of fire events are expected to intensify [99], having a resilient PA system should be a priority. However, the current conservation strategy and management agenda for future conservation projects in Crete is based on the obligations of Dir. 92/43/EEC and especially for habitats and species of Annex I and II, respectively. These obligations include reporting for the conservation status and future trends for only a few species identified as under extinction risk by the present study. Recently, a national set of indicators for assessing ecosystem condition and ecosystem services was set and includes biodiversity-related indices [100], including the identification of micro-refugia of floristic and endemic diversity; our work proposes the inclusion of one more indicator related to future hotspot areas of priority importance for protection and/or conservation interest. It is also evident that hesitancy towards anything other than conventional conservation actions persists [101]. Considering the worst-case scenario (i.e., the extinction of the highest number of species by climate change impacts), in situ conservation focused on micro-reserves and ex situ conservation practices should be considered, as an insurance policy against such losses of plant biodiversity [102], which constitute cost-effective conservation measures [103]. Irreplaceability is a measure of the conservation value of an area [104] designated as PA and/or climatic refugium. Both PAs and climatic refugia in Crete have high mean irreplaceability value, with the latter projected to have higher mean irreplaceability value in the future (Figure 7). The post-2020 biodiversity protection agenda [105] will focus on expanding the existing PAs [106,107]. Many PAs are exposed to spillover effects due to land-cover change [108] and are facing increased climate-change-based risks [109]. It would thus be prudent to aim the conservation efforts at areas with overlaps among PAs and climatic refugia, simultaneously characterised by high diversity and EDGE scores (including the taxa with the highest EDGE scores). These areas are qualified as future climatic refugia and may actually constitute Anthropocene refugia [110]. By doing so, this "climate-smart", cost-effective, conservation-prioritisation planning [111,112] will allow the preservation of evolutionary heritage, trait diversity and future ecosystem services for human well-being [23,104]. Thus, in Crete, under this framework, if at least one area should be given immediate conservation priority due to limited funds, this would be the PA in the Lefka Ori mountain massif, since it presents very high plant diversity, EDGE and ∆EDGE scores, while constituting a climatic refugium almost across its entirety (Table S7; Figures 3 and 6; Figures S8 and S9).
Landscape connectivity as a prerequisite for successful migration is an important management issue. The current scale of habitat fragmentation is likely to hamper the potential success of migration, which relies heavily upon the connectedness of populations across suitable environments [113]. Thus, measures to mitigate landscape fragmentation, as an active conservation strategy, should be included in ongoing and future Action Plans for habitats and species based on climate-change projections and management scenarios. To properly assess and predict future projections of landscape fragmentation, changes in demand, use and supply of ecosystem services should be taken into account, using a climate-change, management scenario-based approach [114]. This procedure is important for policy and decision-making related to land and resource use [115].
Future climatic projections can trigger the interest of the scientific community and the conservation society. However, for decision-and policy-making, raw scientific information is not always the appropriate means of communication. All the information produced by the present study should be transformed into tangible material using the ecosystem services approach and, by this, communicate the loss (via plant species extinction) of various ecosystem services and the related impacts on the socio-economic environment. The Mapping and Assessment of the Ecosystems and their Services (MAES) concept [116] encapsulates this idea in the EU and in Greece [117] and acts as a management and dissemination tool among scientists, policy makers and decision makers.

Conclusions
The present study revealed potential future challenges to be faced under scenarios of climate-change impacts on plant diversity in Crete. The high-altitude areas and mountaintops of Crete emerge as the most important and the most threatened regarding their plant diversity composition and the extinction risks of endemic plant species. Present protection and conservation schemes (i.e., the Natura 2000 network) could benefit from the outcomes of this study and integrate its results into protected areas' management and monitoring specifications. The areas identified as climatic refugia that overlap with PAs and that have high plant diversity and EDGE scores constitute Anthropocene refugia. They should thus be prioritised and supported by relevant action plans to maintain and/or ameliorate habitat conditions. By doing so, these areas could effectively support "climate-smart", cost-effective conservation strategies. Future steps may include pilot studies on specific plant taxa for climate change mitigation measures; flagship species (e.g., Horstrissea dolinicola, the sole representative of a Cretan endemic genus) should be selected to easily communicate the results and encapsulate local knowledge, supported by inhabitants, and, by this, trigger public consultation on relevant policy decisions.  Figure S5: The proportion of total deviance explained as attributable purely to climate (grey), purely to geography (black) and purely to soil pH (olive green). Figure S6 Species turnover as calculated through the GDM framework for the SIE assemblages between the current time-period and. Figure S7: Sites with significantly higher or lower than expected phylogenetic diversity (blue and red shades, respectively), Figure S8: Map of Crete showing SIE assemblages with their respective EDGE index. (A) The EDGE scores are based upon the phylogenetic tree that was built under the framework outlined in [52,53]. (B) The EDGE scores are based upon the phylogenetic tree that was built under the framework outlined in [51].  Table S1: The 172 Cretan single island endemics included in the present study, along with information on each taxon's life form, preferred habitat, extinction risk status for every time-period and climate change model/scenario, extent of occurrence, area of occupancy and niche breadth. Information regarding the area loss for each taxon and every climate change model/scenario is also presented, as is the most important variable for each respective taxon's current potential distribution. NB: Niche breadth. ED: Evolutionary Distinctiveness. EDGE: Evolutionary Distinct and Globally Endangered index. HCP: High conservation priority. EOO: Extent of Occurrence. AOO: Area of Occupancy. ER: Extinction risk. PACA: Preliminary Automated Conservation Assessment. PET: Potential evapotranspiration. RCP: Representative Concentration Pathway. TSS: True Skill Statistic. Table S2: Median range contraction values for every Global Circulation Model (GCM) and Representative Concentration Pathway (RCP) included in the analyses. Table S3: Median range contraction values for every plant family included in the analyses. Table S4: The eleven single island endemic plant taxa of Crete that should be prioritised in terms of conservation efforts based on the EDGE index. Table S5: Median altitude and pH for sites having higher or lower than expected phylogenetic alpha diversity (PD) as well as for the non-significant sites. NS: not-significant. Over: sites with higher than expected PD. Under: sites with lower than expected PD. Table S6: Extent (km 2 ), species richness (SR) and the standardised effect scores of phylogenetic alpha diversity (SES PD ) of the areas identified as climatic refugia in Crete for the present and every Global Circulation Model (GCM) and Representative Concentration Pathway (RCP) combination. A: Asterousia (located in south-central Crete). B: East. C: East-Central. D: Eastern. E: West. F: Western. Table S7: Percent overlap (%) between the protected areas (PA) network and the climatic refugia (CR) recognised in Crete. The extent (in km 2 ) of each climatic refugium is also presented. Phylogenetic tree following [51].tre, Phylogenetic tree following [52,53].tre.