Climate change may cause the extinction of the butterfly Lasiommata petropolitana in the Apennines

Climate change represents a threat to narrow-ranged mountain species living in low-altitude massifs. We studied the disjunct Apennine population of Lasiommata petropolitana (Lepidoptera, Nymphalidae) in the Gran Sasso and Monti della Laga National Park. We quantified the altitudinal shifts undergone in the last decades (1964–2021) in the Alps and Apennines and estimated the local extinction risk due to climate change. We also sequenced the COI mitochondrial marker of seven Apennine specimens, comparing them with those available across the Palearctic. We projected the probability of presence for the species under a future climatic scenario using an ensemble forecasting approach. We found that, despite geographical isolation, the Apennine population of L. petropolitana displays a single widespread COI haplotype also occurring in most European populations. In the Alps and Apennines, this species has shifted uphill an average of 6.3 m per year since 1964. Accordingly, our model predicted a likely extinction in the Apennines by about 2060, due to a reduction of the climatic suitability in this region. Implications for insect conservation Despite its potential loss in the Apennines would not erode mitochondrial diversity, L. petropolitana characterises the butterfly community of the Gran Sasso massif as an alpine enclave. The loss of the Apennine population, together with those of other orophilous butterflies, could trigger a homogenization of alpha and beta diversity and induce a loss of functional diversity in the impoverished high-altitude biotas. As habitat heterogeneity is a key aspect for populations to endure climate change, the maintenance of varied microhabitats, mainly through grazing management, could address the decline of this population.


Introduction
Changes in climate have produced dramatic alterations in species distributions over geological times (Peterson and Lieberman 2012;Blois et al. 2013). The recent humaninduced climate changes represent a growing threat to biodiversity (Halsch et al. 2021;Habibullah et al. 2022). Their impacts are now widely documented (Parmesan 2006;Renner and Zohner 2018;Shah et al. 2020) and are especially noticeable in European southern mountain ecosystems, where populations that were more widely distributed during the coldest phases of the late Quaternary, are now restricted to the highest altitudes (Grabherr et al. 2010;Ernakovich et al. 2014;Ohlemüller et al. 2008). In these areas, even a moderate increase in temperature could determine significant changes in community composition and species richness (Shah et al. 2020;Viterbi et al. 2020).
Among invertebrates, butterflies are considered good models to study climate change effects because their distribution has been generally well studied across time and because they often have a relatively small climatic niche (Parmesan et al. 1999;van Swaay et al. 2010;Habel et al. 2011;Schweiger et al. 2014). Moreover, butterflies tend to react quickly to environmental changes by virtue of biological features such as their complex demands during life cycle, their mobility and short lifespan (van Swaay et al. 2010;Bonelli et al. 2018). Different consequences of climate change have been documented for mountain butterflies, including distribution shifts towards higher altitudes (Forister et al. 2010;Rödder et al. 2021) and an overall expansion of thermophilous and generalist species, in parallel to a decrease of more specialised ones (Zografou et al. 2014;Cerrato et al. 2019;Macgregor et al. 2019;Bonelli et al. 2022).
Many butterfly species inhabiting relatively cold habitats were presumably more widely distributed during the last glacial maximum and in the current interglacial they shifted towards high latitudes and/or to scattered mountain areas. Examples among Holarctic butterflies are included in the genera Erebia (De Groot et al. 2009;Hinojosa et al. 2018;Minter et al. 2020;Sistri et al. 2022) and Parnassius (Ashton et al. 2009;Todisco et al. 2010).
Upward shifts are associated with the reduction and the fragmentation of suitable habitats for orophilous butterflies, increasing their extinction risk because of fitness lowering and genetic impoverishment (Schmitt and Hewitt 2004;Habel et al. 2011). In the short term, an elevational displacement might have a positive effect on a given butterfly population, as at high altitudes land use intensity is generally lower than in lowland areas, but in the long run climate change is expected to lead to a reduction in habitat area (Hülber et al. 2020).
Indeed, the possibility to shift towards higher elevations is clearly limited for populations inhabiting mountains, since the area of the habitats rapidly decreases with altitude (Habel et al. 2011). In addition to ecological consequences of population extinction, the loss of intra-specificgenetic diversity is an important aspect to consider, since many relict populations are markedly genetically distinct from the other lineages, as a result of their protracted isolation (Schmitt et al. 2005;Konvička et al. 2014;Minter et al. 2020;Sistri et al. 2022). A high extinction risk has been predicted for several species with disjunct populations in the Apennines, such as Erebia gorge (Hübner, [1804]) (Piazzini and Favilli 2020) and Erebia pandrose (Borkhausen, 1788) (Sistri et al. 2022), the latter represented by a genetically differentiated population on these mountains.
In this study, we focused on the Apennine population of the butterfly Lasiommata petropolitana (Fabricius, 1787) (Nymphalidae, Satyrinae). This species occurs from the Pyrenees to Scandinavia and as far as eastern Russia (Bozano 1999). While its distribution is relatively continuous in northern parts of its range, towards the south, particularly in Europe, it becomes fragmented into a series of increasingly isolated populations restricted to mountainous areas (Kudrna et al. 2011). In Italy, the distribution of L. petropolitana is similar to that of other butterfly species such as Polyommatus eros (Ochsenheimer, 1808), Boloria pales ([Denis & Schiffermüller], 1775), Erebia pandrose and E. pluto (Prunner, 1798), which are widely distributed in the Alps and restricted to a few mountain areas in the Apennines. These isolated populations often represent endemic lineages responsible for a strong genetic distinctiveness between the communities of the Apennines and the Alps (Ehl et al. 2020;Menchetti et al. 2021;Schmitt et al. 2021;Sistri et al. 2022). In the Apennines, L. petropolitana only occurs in a few sites in three neighbouring massifs: Monti della Laga, Gran Sasso and Monti Marsicani (Balletto et al. 2007). These areas are included within the National Park of Gran Sasso and Monti della Laga and that of Abruzzo, Lazio and Molise. The presence of the species in the Apennines was confirmed in the 1980s after the first report dating back to the end of the XIX century (Prola et al. 1978). Then, this population was described as a separate subspecies (Pararge petropolitana centralapennina Teobaldelli and Floriani 1983), mainly on the basis of wing pattern. In Italy, L. petropolitana inhabits grassy habitats in woodland (e.g., with Fagus sylvatica L. and Larix decidua Mill.) and shrub areas (e.g., with Pinus mugo Turra), mainly on rocky substrates and in cold and humid sites (Prola et al. 1978;D'Alessandro et al. 2008;Hellmann and Parenzan 2010;Bonato et al. 2014). In the Apennines the species is monovoltine, with a flight-period between late May and early July, occurring above 1200 m of altitude in areas close to the limit of beech forests (Prola et al. 1978;Teobaldelli and Floriani 1983;Balletto et al. 2007).
Our study has two main research goals: (i) to assess the future extinction risk of this population and ii) to evaluate the potential impoverishment of mitochondrial genetic diversity if the Apennine population of L. petropolitana would be lost. To achieve these goals, (a) we applied Species Distribution Modelling (SDM) to high-resolution occurrence data to predict the persistence of L. petropolitana in the Apennines during the next decades; (b) we evaluated the historical elevational shift of the species over the Alps and Apennines; and (c) we obtained cytochrome c oxidase subunit 1 (COI) sequences to evaluate the mitochondrial DNA (mtDNA) genetic identity of the Apennine population with respect to other populations from Europe and Asia.

Data collection
We surveyed the Central Apennines in 2021, recording the geographic coordinates and the elevation (10 m precision) of all the sites where we observed L. petropolitana (Online Resource 1). Occurrence data for the Alps and the Apennines, as defined by Menchetti et al. (2021) (Fig. 1), were retrieved from literature, museum and private collections with reference to the Italian CKmap project (Balletto et al. 2007) and from collection data stored in Roger Vila's laboratory at the Institut de Biologia Evolutiva CSIC-UPF (Barcelona, Spain). Data deposited in the platform iNaturalist have also been included if they were associated with high-resolution coordinates (error lower than 1000 m) and if accompanied by a photograph in order to avoid the risk of confusion with the similar L. maera (Linnaeus, 1758). A subset of occurrences was obtained using only the specimens with year and altitude recorded. We found a total of 147 suitable occurrence data. Due to the different sources of these data, and to reduce clustered occurrences, we then filtered the dataset selecting a minimum distance between the occurrence points of 2 km using the 'thin' function of 'spThin' R package (Aiello-Lammens et al. 2015) in R v. 4.0.3 (R Core Team, 2020. This reduced the final dataset used for the Species Distribution Modelling to 109 records ( Fig. 1).

Species distribution modelling
We modelled the distribution of L. petropolitana in the study areas as a single entity, although Alps and Apennine populations may be differently adapted to local conditions. Recent studies highlighted the importance of including adaptive genetic variation in climate change vulnerability assessment (Razgour et al. 2019). However, given the very limited distribution of L. petropolitana in the Apennines (4 raster cells) this was not possible. Modelling the entire species as a single entity is expected to generate more generous estimates of species tolerance to environmental conditions.
In order to reduce the risk of overfitting through the inclusion of spurious correlates , which can be particularly problematic in climate change projections of species distributions (Santini et al. 2021), we restricted the selection of predictor variables to those for which we could hypothesize a causal relationship with species presence. There are few studies on physiological and life history traits of L. petropolitana and they do not clearly suggest which climatic variables are most influential on the species (Nylin et al. 1996;Gotthard 1998). However, previous research on butterfly physiology and distribution documented some major climatic constraints in Southern Europe. In general, the loss of snow cover strongly affects the survival of butterflies at high altitudes due to exposure to low winter temperatures (Vrba et al. 2012(Vrba et al. , 2017Konvicka et al. 2021) while, in summer, the main influencing factor is water availability (Vodă et al. 2015), which also counteracts the strong seasonal increase in temperature. We obtained 19 climatic variables from the CHELSA database (Karger et al. 2017(Karger et al. , 2018) (Online Resource 3) with a resolution similar to that of the occurrence data (30 s, ~ 1 km 2 ), and clipped them for the study area. Then we carried out a PCA on the variables to obtain their correlation structure (Online Resource 4). Afterwards, based on the available biological knowledge and the correlation structure among available climatic variables, we selected four weakly correlated variables: 1) minimum temperature of the coldest month ('bio6' in CHELSA database); (2) precipitation of the coldest quarter ('bio18'); (3) precipitation of the warmest quarter ('bio19'); (4) temperature annual range ('bio7'). Regarding the climatic variables in the study region, in addition to a generalised increase in temperature, there is evidence of a decrease in precipitation frequency and amounts in the last decades (Valt and Cianfarra 2010;Di Lena et al. 2012;Gobiet et al. 2014;Brugnara and Maugeri 2019).
We evaluated the Variance Inflation Factor (VIF) for each selected variable using the R package 'usdm' (Naimi et al. 2014) and verified that all values were < 3, which indicates very low multicollinearity (Alin 2010; Kock and Lynn 2012) (Online Resource 5). Then, we downloaded and clipped the same variables for future projections for the 2041-2060 period, as well as the Representative Concentration Pathway (RPC) 4.5, described by Moss et al. (2008)  Species distribution modelling was performed using the R package 'biomod2' (Thuiller 2014). In general, overall statistical models tend to be less prone to overfitting than machine learning models, but the latter usually achieve higher accuracy in present predictions (e.g. Merow et al. 2013). As such, it is recommended to use multiple and diverse algorithms (Araújo and New 2007;Norberg et al. 2019) and validate them using a robust blocking approach to estimate the uncertainty around predictions (Bahn and McGill, 2013;Roberts et al. 2017). Accordingly, we used four statistical models to forecast the current and future species distribution: generalised additive model (GAM), generalised linear model (GLM), maximum entropy (MAXENT) and random forests (RF). The first two are statistical linear models, while MAXENT and RF are machine learning algorithms, so the entire set covers a wide range of complexity in distribution models (Merow et al. 2013) allowing appropriate estimates of uncertainty in model projections (Buisson et al. 2010).
Thereafter, we generated 10 datasets of 1000 randomly selected background points, one dataset for each of 10 repetitions used for the model validation. To estimate the relative importance of the predictor variables, we correlated the predicted probabilities of presence, using the full dataset, with the predicted probabilities of presence obtained after permuting the variable of interest. The relative importance of each variable was quantified as one minus the Pearson rank correlation coefficient (Thuiller et al. 2009). We validated the models using spatial block validation with the R package 'blockCV' (Valavi et al. 2019). We split the study area into six spatial blocks and two folds (Online Resource6), iteratively fitted the model on all but one block, and tested against the left-out block. This procedure is known to provide a more objective assessment of model transferability for future projections (Bahn and McGill 2013;Roberts et al. 2017). We measured model performance using the True Skill Statistic (TSS) and the area under the curve (AUC). For both present and future projections, we obtained one occurrence probability raster for each statistical model by calculating the mean of all the projections of models with a TSS > 0.5 and an AUC > 0.7. In order to obtain a single projection for each scenario, we then averaged the four rasters and excluded cells with an occurrence probability < 10%. In the results, the term 'cells' will thus indicate the cells with an occurrence probability > 10%. Then we calculated the difference between the two scenarios by subtracting the current average predictions from the future ones; raster cells with positive delta values indicate a predicted improvement of climatic conditions, whereas raster cells with negative delta values indicate a deterioration of climatic conditions. To estimate the uncertainty in the predictions due to disagreements among the four different algorithms, we assigned − 1 to all cells with negative values, + 1 to all cells with positive values (and 0 otherwise) of the average single-model predictions, and assessed the consensus of model predictions by summing the four binarized maps. This resulted in a raster map with values ranging between − 4 and + 4, with extreme values indicating that all the four statistical models predicted a decrease (− 4) or increase (+ 4) in the probability of occurrence, and intermediate values indicating partial (± 2-3) or high disagreement (− 1 to 1) among the predictions of the algorithms.
To ensure transparency and full replicability of the methods we provide an Overview, Data, Model, Assessment and Prediction (ODMAP) protocol as part of the supplementary material, that details all methodological choices and parameters (Online Resource 7; Zurell et al. 2020). The occurrence data used for species distribution modelling are available in the Online Resource 8.

Changes in elevation over time
Lasiommata petropolitana shifts in elevation were evaluated from 1964 to 2021, assuming 1964 as one of the dates suggested to mark the beginning of the Anthropocene (Lewis and Maslin 2015). When explicitly reported, elevation values were obtained from literature and collections data. For the occurrences without a clear report of elevation, but with precise coordinates (error lower than 1000 m like most citizen science data), the elevations were extracted using a 1 km 2 layer available at https:// www. earth env. org/ topog raphy (Amatulli et al. 2018). The final dataset was composed of 323 records over a period of 57 years. An ordinary least squares (OLS) regression between time and elevation was used to test for the existence of a temporal trend in elevational shifts.

COI analyses
To evaluate the genetic identity of the Apennine population of L. petropolitana we gathered COI sequences for this species from across its distribution range (Eurasia) by downloading data publicly available in GenBank and in BOLD (Barcode of Life Data System), and by adding a series of seven unpublished sequences specifically analysed for this study, representing specimens collected in 2021 in Monti della Laga and Gran Sasso massifs. COI sequences from the Apennine specimens were obtained using LepF1 and LepR1 primers, following the protocol by Platania et al. (2020) at the Institut de Biologia Evolutiva CSIC-UPF (Barcelona, Spain). We aligned all the sequences with MegaX. We removed potential contaminant sequences as well as sequences shorter than 600 bp. This resulted in a final alignment of 55 specimens (Online Resource 2). We constructed maximum parsimony haplotype networks using the COI alignments and TCS 1.21 (Clement et al. 2000) by setting a 95% connection limit (11 steps). The haplotype network was graphically edited with tcsBU (Múrias dos Santos et al. 2016).
Haplotype colours were assigned according to p-distances projected after Principal Coordinates Analysis (PCoA) on the RGB space as done by the 'recluster.col' function of 'recluster' R package https:// cran.r-proje ct. org/ web/ packa ges/ reclu ster/ index. html) (Dapporto et al. 2013). Subsequently, the colours were plotted on a map using the recluster.plot.pie function of the same package. The haplotypes in the network were coloured according to their geographical origin.

Results
We recorded L. petropolitana in four sites in the massifs of Gran Sasso and Monti della Laga within the borders of the homonymous National Park. The species was found in grassy areas on rocky substrate in beech forest clearings, as well as in adjacent pastures with sparse shrubs, in an altitudinal range between 1455 and 1775 m above sea level (Online Resource 1). Among literature data for the Apennines, we excluded the single pre-1980 record for the massif of Monti Marsicani in the Abruzzo, Lazio and Molise National Park (Prola et al. 1978), as it lacks both coordinates and altitude. In the latter area, which constitutes the southern range limit of L. petropolitana in the Apennines, we have not detected the species and therefore its current presence needs confirmation.

Species distribution modelling
After thinning a total of 147 occurrence data, a dataset of 109 records was selected for the species distribution modelling of Alpine and Apennine populations (Fig. 1).The spatial-block validation of the statistical models indicated good predictive performances, with TSS = 0.53-0.82 and AUC = 0.81-0.93 (Online Resource 9). The most important variable in predicting the occurrence of L. petropolitana, according to all four statistical models, was the minimum temperature of the coldest month (values of the mean importance of all variables are given in the Online Resource 10). Accordingly, the smoothed response curves showed a strict dependence of L. petropolitana occurrence probability on the minimum temperature of the coldest month, with a great concordance between different algorithms. Conversely, the three other variables showed weaker relationships and lower concordance between the algorithms (Online Resources 10 and 11).
The projections of each statistical model (Online Resource 12) pictured slightly different results that were averaged in the ensemble model. The ensemble model projections for the present and the future indicated a high probability of occurrence across the Alps, but not over the Apennines (Figs. 2a, b). The predicted change in the probability of occurrence raster highlighted a widespread decline in climatic suitability across the distribution range of the species in the Apennines (Fig. 2c). We also found a high consensus among model predictions, indicating a high probability of decline throughout the current distribution of the species (Fig. 2d). Similarly, many marginal areas of the pre-Alps are predicted to deteriorate their climatic suitability for L. petropolitana (Fig. 2c).
The occurrence suitability for both present and future of all the cells is reported in Fig. 2e. The four cells assessed for the Apennines are indicated by numbers from one to four (Fig. 2e). These cells will shift from a present mean value of 0.47 ± 0.12 mean ± s.d.) to 0.24 ± 0.15 (mean ± s.d.) in the future (Online Resource 13). Only 2,7% of the occupied cells have a probability of 0.24 in the present.This value is much lower than the 10% quantile threshold used in similar studies 1 3 to predict future occurrence at a given site (e.g., Habel et al. 2011).Overall, 75.4% of the cells will decrease their suitability in the future (Online Resource 14). More than half of the cells (54.1% + 19.3%) have a high consensus value between the four algorithms (− 4, − 3; + 3, + 4 respectively; Fig. 2f) and belong to high elevation areas (Fig. 2d).
Regression between elevation and year of records indicated an average increase of 6.3 m (± 1.5 s.e.) in altitude per year since 1964 (t = 4.28, P < 0.001, Fig. 3).

COI structure
We retrieved a total of 48 COI sequences with a length between 600 and 658 bp and obtained seven original sequences (all with 658 bp) for the four Apennine sites. The haplotype network analysis identified few (8) haplotypes, with a star-like distribution (Fig. 4c). One haplotype appeared as highly frequent and widespread, while the other 7 haplotypes are rarer, separated by one or two mutations from the main haplotype.
While the most common haplotype is widely distributed across the range of L. petropolitana, including the isolated Apennine population, rare haplotypes are segregated in different regions. If we exclude the 5 singletons, the 2 remnant haplotypes are endemic to the Alps and the Pyrenees (Fig. 4c).

Discussion
According to species distribution modelling and the analysis of the altitudinal trend, the population of L. petropolitana in Central Italy is strongly threatened by climate change, as there is evidence that the climatic suitability of the areas currently occupied in the Apennines will lower critically in the next decades.
The analysis of the COI mitochondrial marker showed that the Apennine population belongs to a common haplotype occurring all over the distribution range of the species, therefore questioning the validity of its subspecific status attributed by Teobaldelli and Floriani (1983). This pattern is in agreement with recent findings on this genus (Platania et al. 2020), and reflects the tendency of boreoalpine species to show a limited genetic differentiation after postglacial dynamics of populations (Hewitt 2000;Mutanen et al. 2012;Wallis and Arntzen 1989). Among the selected variables the minimum temperature of the coldest month ('bio6' in CHELSA database) explained most of the distribution of the species across the Alps and Apennines (Online Resources 10, 11), with an optimum value around − 10 °C. During the current interglacial these temperatures only occur on high mountain systems in Southern Europe. Presumably, L. petropolitana was more largely and quite homogeneously distributed in non-glaciated areas of Central and Southern Europe before the onset of the current interglacial period that started 11.600 years ago (Otto-Bliesner et al. 2021). During the interglacial period, L. petropolitana should have quite rapidlyexpanded far to the north in areas formerly covered by ice-sheets (Scandinavia) and retreated to scattered mountain areas in Central and Southern Europe (Pyrenees, Alps, Apennines, Balkans, Carpathians) (Mutanen et al. 2012). Likely, this process is currently exacerbated by human-induced climate changes and we found a significant elevational shift of about 6 m per year, in line with results obtained with Lepidoptera in other studies (Forister et al. 2010;Rödder et al. 2021;Sistri et al. 2022). A consistent shift in altitude may force the Apennine population to occupy smaller areas still characterised by locally suitable environmental conditions. Such areas may indeed act as microclimate refugia in the short term, but in the long run their suitability will be reduced due to climate change, probably leading to the local extinction of the species. On the Apennines, the minimum altitude recorded in recent years for the species is 1455 m a.s.l. (site 4, Online Resource 1). The relatively low altitude of the mountains in its range determines that suitable sites are isolated from each other by valleys that do not conform to the species' demands in terms of climate and habitat, thus limiting the possibility of dispersal among the mountain groups.
In these areas L. petropolitana lives near the treeline, a limit around which there is increasing evidence for butterfly decline in the Apennines (Scalercio et al. 2014;Cini et al. 2020). The local morphology of the territory could be also a limiting factor for the dispersal of the species, such as in areas where rocky cliffs enclose the suitable habitats, Fig. 2 Species distribution modelling of Lasiommata petropolitana (cells with an occurrence probability < 10% were excluded). The panels (a-d) are accompanied by insets showing a detail of the Central Apennines area and four yellow arrows indicate the position of the four occurrence records for this area. Panels (a) and (b) show the current  and future (2041-2060) distribution of the species, respectively. They are based on the ensemble models, ensemble present (a) and ensemble future (b), and the colours represent the species occurrence probability in a range between 0 and 1. c Ensemble delta: this panel shows the occurrence probability difference between future projections and current projections (future probability minus present probability for each cell). Negative values (in red) indicate a predicted reduction of suitability in the future, whereas positive values (in blue) indicate a predicted increase of suitability in the future (the values between -0.1 and + 0.1 are excluded). d Consensus change. This panel shows the confidence in directional change of the four statistical models. The colours represent the number of statistical models which agree in the future occurrence predictions: green areas indicate full agreement among all the models (+ 4 + 3 − 3 − 4 values), while red areas indicate an intermediate agreement in which different models predict different results (+ 1 + 2 0 − 1 − 2 values). e Histogram showing the proportion of raster cells depending on the probability of the occurrence of L. petropolitana in present (light blue) and future (red with white stripes). The numbers (1, 2, 3, 4) indicate the occurrence probability values of the four sites in the Apennines both for the present (blue numbers) and for the future (red numbers). These numbers are associated with those in Online Resources 1 and 13. f Frequency of the confidence in directional changes among the four statistical models represented as the percentage of cells for which models agreed. The consensus categories are indicated on the x-axis and their percentages on the y-axis (the colours are associated with (d) ◂ without allowing connections with nearby ones (e.g., site 4, Online Resource 1).
Another important factor responsible for the range shift of many butterfly species is land-use change (van Swaay et al. 2010;Bonelli et al. 2011). We have no information on the historical occurrence of L. petropolitana in the Apennines, as its presence there has been confirmed in the 1980s (Prola et al. 1978). Furthermore, recent records are few and we don't know if anthropic activities, e.g. linked to livestock or forestry activities, have contributed to its rarity. At present, however, the distribution range of the species is included in a National Park, where land use changes are limited thanks to current regulations.
Despite being clearly threatened at a local scale, L. petropolitana is assigned to the Least Concern category in the European Red List of Butterflies (van Swaay et al. 2010). Nevertheless, in some countries it is recognised as a species under extinction risk, being listed as endangered in Germany (Reinhardt andBolz 2011) andPoland (Buszko andNowatzki 2000). In the Czech Republic it is now extinct, following habitat disappearance due to land use changes, including cessation of grazing and afforestation (Spitzer et al. 2017). These examples refer to populations at the margins of their main distribution range, comprising massifs with lower altitudes than the Alps and the Pyrenees and/or areas located at the southern limit of the Northern Europe population. This tendency has been highlighted in other butterflies at their range limit (Wilson et al. 2015;Fourcade and Öckinger 2017). Lasiommata petropolitana is assigned to the Least Concern category also in the Italian IUCN Red List, due to the presence of large and relatively stable populations in the Alps (Bonelli et al. 2018), but the marginal population in the Apennines should be categorised as Endangered according to the B2ab (iii, v) IUCN criteria (http:// www. iucn. it/ scheda. php? id=-16629 73159) and our projections for the next decades confirm this assessment.
Although L. petropolitana does not show a particularly conspicuous variability in the COI marker across its range, some of its biological features, probably having a genetic basis, are known to vary among different geographical areas, such as larval diapause in alpine populations in contrast to pupal diapause reported further north (Nylin et al. 1996;Gotthard 1998). However, mtDNA may show an evolutionary trajectory that may be different than that of nDNA. From this point of view, it would be worth checking in future studies whether the Apennine population represents an Evolutionary Significant Unit (Casacci et al. 2013) based on highthroughput sequencing on the nuclear genome.
Even if the potential loss of the Apennine population of L. petropolitana in the future would not erode the mitochondrial genetic diversity of the species, the prospect of its local extinction is a cause for concern. Indeed, L. petropolitana is a species that characterises the rich butterfly community of the Gran Sasso and Monti della Laga National Park. Together with communities from neighbouring highaltitudes mountains, these massifs are recognised as alpine enclaves in the Apennines, because of both species' richness and composition (Balletto et al. 1982). Among the most  Reissinger, 1990, Lycaena hippothoe (Linnaeus, [1760) and Erebia euryale (Esper, 1805), all fairly widespread in the Alps but more localised in the Apennines (Balletto et al. 2007;Dincă et al. 2011).
The loss of L. petropolitana in the Apennines, together with other orophilous butterflies (Bonelli et al. 2011) could trigger a homogenization of alpha and beta diversity in the region and possibly induce a loss of functional diversity in the impoverished high-altitude biota. In turn, these losses could erode the ecosystem services provided by a restricted number of butterflies species occurring on high altitude grassy habitats (Noriega et al. 2018). Highly charismatic butterfly species that characterise unique communities can also act as attractive components for tourism in protected  Fig. 2a. c Maximum parsimony COI haplotype network. The haplotypes (coloured circles) are separated by segments; each segment represents a single mutation between two different haplotypes. Haplotype colours have been assigned according to the geographical origin of the COI sequences. The circle area is proportional to the number of specimens showing each haplotype, as shown in the legend areas, and their persistence could benefit the touristic offer of the Gran Sasso and Monti della Laga National Park.
It is difficult to hypothesise any action to avoid the climate change impacts on the Apennine population of L. petropolitana. Translocation to other Apennine massifs makes little sense since our predictions indicate a decrease of the climatic suitability for the entire mountain chain in the next decades, and no positive trend is expected anywhere in the region (Fig. 2c). However, the entire range of the species in the Apennines is included within protected areas, where it is possible to implement conservation measures. In general, the maintenance of habitat heterogeneity is a key aspect for the conservation of rich butterfly communities (van Swaay et al. 2012). In the study area, L. petropolitana partly occurs in grassy habitats maintained by cattle and sheep grazing, which prevents the expansion of shrubs and arboreal vegetation near the treeline. The maintenance of a moderate grazing, as long as it is correctly planned and adapted to local conditions (Metera et al. 2010;Jerrentrup et al. 2016;Johansson et al. 2017;Cutter et al. 2022), is thus desirable at the wood margins, considering that the treeline generally tends to advance rapidly both with land abandonment and climate warming (Gerhig-Fasel et al. 2007;Wieser et al. 2014).However, there are also large upland areas in Gran Sasso and Monti della Laga National Park, including some where L. petropolitana occurs (e.g. sites 1 and 2, Online Resource 1), which are subject to intense grazing both by cattle and sheep. A decrease of the treeline elevation is reported on the southern slopes of Gran Sasso massif, probably due to the combination between water shortage in summer and alteration of microclimate conditions when beech cover was reduced for pasturage in the past (Bonanomi et al. 2020). At these sites, some efforts should therefore be dedicated to allow the natural advance of the beech woods, in order to limit the mismatch between the elevational shift of L. petropolitana and the availability of suitable vegetation structure units for the species. Moreover, in such areas it could be possible to implement local strategies to limit the impact of grazing on herbaceous vegetation and on the biodiversity associated with it. Among the insects occurring in the Park, other localised species would also benefit from this measure, such as the endemic bumble bee Bombus konradini Reinig, 1965, for which intensive grazing is recognised as a major threat (Quaranta et al. 2018).
Continuous monitoring will also help in detecting negative population trends and prime for the need to implement changes in management. Italian National Parks, including Gran Sasso and Monti della Laga are called upon to conserve biodiversity, including threatened butterfly populations (Cini et al. 2020;Sistri et al. 2022). The same Parks have been recently committed by the Italian Environment Ministry to monitor pollinators using periodic transects following the EU Pollinators Initiative. The information obtained by transect counts will allow the monitoring of population trends and, in the case of species threatened by climate change such as L. petropolitana, they will enable to evaluate whether a reduction of estimated climatic suitability is associated with a reduction in population size.
Acknowledgements This study has been carried out in collaboration with the Gran Sasso e Monti della Laga National Park within the project "Ricerca e conservazione sui lepidotteri diurni di sei Parchi Nazionali dell'Appennino Centro-Settentrionale". MB and LD thank Giorgio Davini (Gran Sasso and Monti della Laga National Park) for his contribution in the project. VD also acknowledges the Visiting Professor fellowship awarded by the Research Institute of the University of Bucharest. Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement. This study was funded by the Ministero Italiano della Transizione Ecologica within the project "Ricerca e conservazione sui lepidotteri diurni di sei Parchi Nazionali dell'Appennino Centro-Settentrionale". Support was also provided by the Academy of Finland (Academy Research Fellow, decision no. 328895) to VD. RV is supported by Grant PID2019-107078 GB-I00 funded by Ministerio de Ciencia e Innovación and Agencia Estatal de Investigación (MCIN/AEI/https:// doi. org/ 10. 7554/ eLife. 76576).

Conflict of interest
The authors declare no conflict of interest.

Ethical approval N/A
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.