Potential climatic and elevational range shifts in the Italian narrow endemic Bellevalia webbiana (Asparagaceae) under climate change scenarios

shifts in the


Introduction
Although soil consumption is a major concern for plant conservation, climate change caused by human activities can further exasperate conservation problems (Corlett 2016(Corlett , 2020. Climate change is already impacting biodiversity and is likely to intensify over the next century (Bellard et al. 2012;Casazza et al. 2014), so that mitigation efforts are urgent to implement. In the case of species endemic to a country, the global chance of survival entirely relies on national policies, so that endemic taxa are key elements for setting national conservation priorities and for assigning conservation tasks (Orsenigo et al. 2018). Italy hosts a high number of endemics, between 19% and 21% of the native flora according to different estimates (Peruzzi et al. 2014(Peruzzi et al. , 2015Bartolucci et al. 2018Bartolucci et al. , 2022, and most of them show a quite restricted distribution range. The Webb's hyacinth (Bellevalia webbiana Parl., Asparagaceae, monocots; Fig. 1) is one of the most relevant Italian narrow endemic species (Chiarugi 1949;Borzatti von Loewenstern et al. 2013;Astuti et al. 2017). The distribution range of this bulbous perennial herb is restricted to two relatively small and disjunct areas of Central Italy, separated by the mountain ridge of Northern Apennine. Despite the small distribution range, this species is linked to relatively common habitat types, since it typically grows in open fields and meadows, wood margins, olive groves, and vineyards (Gestri et al. 2010). Recent studies investigated the population genetics , functional traits (Astuti et al. 2019), and reproductive performance ) of the five largest populations known for this species. These studies highlighted some differentiation in the population located at the upper NE limit of the species range, around Faenza (Ravenna). The International Union for Conservation of Nature (IUCN) listed the Webb's hyacinth in the Red List of Threatened Species (Peruzzi and Carta 2011; see also Orsenigo et al. 2018), as Endangered (EN A2c). During the last century, the main threat to this species was represented by the spread of human settlements, which caused local extinctions from several historical localities (Gestri et al. 2010). Climate change can quickly become a further threat for the persistence of the small and localised populations of this species, magnifying its extinction risk in the near future (Casazza et al. 2021). Therefore, to evaluate the impact of climate change on the conservation of Bellevalia webbiana, we modelled its future potential distribution under different scenarios by using Ecological Niche Modelling (ENM) (Khanum et al. 2013;Assis et al. 2018).

Occurrence data and study area
Distribution data for Bellevalia webbiana were obtained from Gestri et al. (2010), and from numerous recent field observations stored in Wikiplantbase #Italia (Peruzzi et al. 2019 onwards). A total of 122 distribution records retrieved from the abovementioned sources were then filtered to retain only occurrence records in a 1×1 km grid cell. After filtration, 53 1×1 km grid cells remained for modelling the distribution of the target species (Fig. 1). The study area was selected based on the known range of the species, corresponding to two relatively small and disjunct sites of Central Italy. This territory, located inside longitude 10.99-11.94 E and latitude 43.65-44.39 N ( Fig. 1), covers an area of 3,192.37 km 2 and is characterised by a large absence area in correspondence of the Apennine mountain ridge, along the direction NW-SE.

Current and future climatic data
Mean monthly climatic data at the same 1 km 2 grid cell used for the occurrence data were downloaded from CHELSA (Karger et al. 2021) for three temporal ranges (1981-2010 'current times', 2041-2070, and 2071-2100). The climatic data were used to calculate three biologically relevant climatic variables (annual mean temperature, annual precipitation, and annual potential evapotranspiration) with dismo and envirem R packages (Hijmans et al. 2017;Title and Bemmels 2018;Dolci and Peruzzi 2022). Future climatic scenarios were derived from IPSL-CM6A-LR climatic models (Boucher et al. 2020) and related to three Shared Socioeconomic Pathways (SSP): SSP1-2.6 (CO 2 emissions cut to net zero around 2075), SSP3-7.0 (CO 2 emissions around current levels until 2050, then falling but not reaching net zero by 2100), and SSP5-8.5 (CO 2 emissions triple by 2075).

Ecological niche modelling
An ensembled method approach (Araújo and New 2007) was selected to calculate all potential distributions. Maximum entropy algorithm implemented in Maxent soft- ware (version 3.3.4) (Phillips et al. 2022) and run via ENMtools (Warren et al. 2021) was used to calculate 100 models. Each model was calibrated randomly splitting the occurrence data in two subsets: 70% of for training and the remaining 30% for testing the predictive performance. Bootstrap was selected as replicated run type because it provides accurate inferences when sample size is small (Freedman 1981;Chen et al. 2019). Other parameters were kept as the default (regularization multiplier = 1, feature class combinations defined on the number of samples) (Merow et al. 2013). Model quality was assessed with two different evaluation metrics: the Area Under the Receptive Curve (AUROC) and the Continuous Boyce Index (CBI) (Hirzel et al. 2006). The CBI score provides valuable information on the model's quality (robustness, probability of occurrence and deviation from randomness) and is the most appropriate accuracy metric in the case of presence-only data (Hirzel et al. 2006). Albeit AUROC was considered the standard method to assess the accuracy of predictions, several drawbacks have been reported when applied to species distribution modelling (Lobo et al. 2008). Consequently, in this study, AUROC is reported just as complementary information.
The final ensemble maps were constructed by averaging all models in which CBI value calculated with testing subset was equal to or greater than 0.8. In total, predictions generated by 24 models were used to calculate mean habitat suitability maps together to the relative uncertainty maps.
To avoid extrapolation errors (Owens et al. 2013), and to enhance environmental similarity between the calibration and projection regions, the models were calibrated on the area within longitude 6.2-19.1 E and latitude 35.5-47.3 N. Only three non-collinear and biologically meaningful climatic variables were used as predictors (Barbet-Massin and Jetz 2014). All habitat suitability maps were converted to binary presence-absence outputs by applying a threshold value based on equal training sensitivity and specificity (Liu et al. 2005). Variable importance was computed by using the feature importance ranking measure (FIRM) approach (Greenwell et al. 2018;Scholbeck et al. 2019).

Geographic range estimation
The Area of Occupancy (AOO) as defined by IUCN Red List (IUCN Standards and Petitions Committee 2022) was used to assess extinction risk under the impacts of climate change. Accordingly, to estimate the potential AOO for future scenarios, all predicted potential distributions were downscaled to 2×2 km grid cells (IUCN Standards and Petitions Committee 2022) and the number of cells corresponding to predicted areas of potential presence was calculated.

Uncertainty estimation
The model replications allow to estimate the uncertainty associated with predictions. For each future scenario, all habitat suitability maps used to generate the ensemble prediction were also used to calculate variability as standard deviation. Then, the overall quality of predictions was evaluated by reporting values of uncertainty versus habitat suitability in a 2D-scatter plot. Key areas are those showing high or low suitability and low uncertainty (Peterson and Samy 2016), corresponding to areas in which a species can be retrieved or not with a high degree of confidence.

Results
The ensembled Maxent model for the Webb's hyacinth with the selected environmental variables performed particularly well, achieving high evaluation metric scores (AUC-Train = 0.932, AUC Test = 0.928, CBI Test = 0.853). In total, 24 single models were combined into a final ensembled model. The current potential Area Of Occupancy (AOO), made up by 992 2×2 km cells, is forecast to decrease in the range 2041-2100 (Fig. 2).
The feature importance ranking measure (Table 2) points out the great importance of Annual Mean Temperature and Annual Potential Evapotranspiration as predictors, even if Annual Precipitations also give a relevant contribution.
The overall uncertainty of predictions for current climatic condition (Fig. 3), calculated as a range of values, indicated a low level of uncertainty among replicated runs.

Discussion
We modelled the potential changes in distribution range of the Italian narrow endemic Bellevalia webbiana under different climate change scenarios, highlighting a significant reduction in its potential area of occupancy. This confirms the estimations made by Casazza et al. (2021), albeit carried out with a different focus and scope. These authors calculated that the assisted colonization of 13-16 new 2×2 km grid cells would be required to halt the estimated loss of range under a pessimistic scenario.
Several authors already evidenced similar future reductions in other stenochorous endemic plant species (Hong-Wa and Arroyo 2012; Wulff et al. 2013;Dülgeroğlu and Aksoy 2019;Khajoei Nasab et al. 2020;Semu et al. 2021), building up a consolidated body of knowledge on the potential risk of climate changes for the conservation of endemic plant species in different biogeographical contexts. The results of this study also reinforce the findings by Walker (2014) and Darrah et al. (2017), who highlighted the high risk of extinction of bulbous monocots, which often have stable local populations but limited dispersal capacity. In addition, the same authors remarked on the general underrepresentation of geophytes in plant conservation projects. In all future scenarios, the predicted potential distribution tends to shift towards higher elevations, which are located between to the two main areas in which this species currently occurs, but where this species is currently absent for other extra-climatic ecological reasons. We highlight that the Webb's hyacinth typically occurs at low elevations (on average ca. 200 m a.s.l.). Despite this, we predicted a tendency to shift towards higher elevations, as already evidenced for many species (Randin et al. 2009;Feng et al. 2020). Several other studies returned similar results, but concerning montane plants Engler et al. 2011;He et al. 2019). In addition, important resurvey studies of historical data demonstrated the upslope shift of several species, resulting in a complete transformation of the local assemblages in the same mountain range (Morueta-Holme et al. 2015;Steinbauer et al. 2018).
However, the upslope shift here modelled for Bellevalia webbiana and the connected reduction in potential area of occupancy is only predicted on the basis of climatic niche, and does not consider other potential variables affecting the ecological niche of the species. The most relevant factor potentially making the future scenario much worse than those modelled on the basis of predicted climate is represented by the availability of habitat. This species is strongly linked to open habitats, often with traditional management (Gestri et al. 2010), but these habitats are disappearing in most of the hills and submontane ecosystems of peninsular Italy because of socio-economic changes and of the concentration of population in lower elevation areas (Falcucci et al. 2007). Forest communities are occupying previously grazed or cultivated lands in the higher hills and lower mountains (see e.g., Geri et al. 2010;Amici et al. 2013), promoting more natural habitat types but affecting the persistence of species linked to secondary and cultural landscapes (Amici et al. 2015).
In addition to habitat loss, the limited dispersal ability is a further factor which should be considered to assess the actual extinction threat of B. webbiana. Indeed, given its relatively large seeds , the dispersal capacity of this species may be insufficient to cope with the speed needed for the upslope shift requested by the climate change. Then, even admitting that potentially suitable areas at higher elevations could effectively host B. webbiana and that this species could colonise these areas by natural or assisted migration (Brodie et al. 2021), then this would imply an increasing connection between the two currently disjunct population groups. This increasing connection would, in turn, promote a potential loss of the genetic and functional differentiation documented for some populations (Astuti et al. , 2019Peruzzi et al. 2021). Moreover, considering that the actual AOO documented for Bellevalia webbiana amounts to just 47 2×2 km cells (of which only 24 were confirmed in recent times) and that forest expansion in Apennine mountains could further reduce the availability of open habitats (Gargano et al. 2012), this species shows evident conservation concerns. To overcome these problems, in addition to assisted colonization (Casazza et al. 2021), a proper ex situ conservation programme should be planned, as an insurance policy against extinction (Li and Pritchard 2009).