A window to the future: effects of climate change on the distribution patterns of Iranian Zygaenidae and their host plants

Climate change has been suggested as an important human-induced driver for the ongoing sixth mass extinction. As a common response to climate change, and particularly global warming, species move toward higher latitudes or shift uphill. Furthermore, rapid climate change impacts the biotic interactions of species, particularly in the case of Zygaenid moths which exhibit high specialization in both habitat and host plant preferences. Iranian Zygaenidae are relatively well-known and represent a unique fauna with a high endemism rate (46%) in the whole Palearctic; as such they are a good model group to study the impact of climate change on future distributions. In this study, we used species distribution models (SDMs) and ensembles of small models (ESMs) to investigate the impact of climate change on the future distribution of endemic and non-endemic species of zygaenids, as well as their larval host plants. Three different climate scenarios were applied to forecast the probable responses of the species to different climate change intensities. Our results suggest that the central and southern parts of the country will be impacted profoundly by climate change compared to the northern regions. Beyond this, most endemic species will experience an altitudinal shift from their current range, while non-endemic species may move towards higher latitudes. Considering that the regions with higher diversity of zygae-nids are limited to mountainous areas, mainly within the Irano-Anatolian biodiversity hot-spot, the identification of their local high diversity regions for conservation practices has a high priority.


Introduction
Alongside habitat degradation and overexploitation, human-induced climate change has been considered one of the main drivers of the sixth mass extinction, which increased the rate of extinction by 100 times compared to the last five mass extinctions (Maxwell et al. 2016;Pimm et al. 2014;Shivanna 2020).Even optimistic scenarios for future climate change anticipate an increase in temperature, the concentration of CO 2 , and greenhouse gases in the near future (Bellard et al. 2012;IPCC 2022;Ripple et al. 2020).It has been well documented that species commonly respond to climate change through changes in morphology, behavior, phenology, and geographical range shifts (e.g.Bellard et al. 2012;Della Rocca and Milanesi 2022;Howard et al. 2023;Rödder et al. 2021).However, pieces of evidence demonstrate that this response can vary among species (Parmesan and Yohe 2003;Weiskopf et al. 2020).Different studies showed that the response of most taxonomic groups to climate warming is an altitudinal and latitudinal retreat (e.g.Biella et al. 2017;Dieker et al. 2011;Hickling et al. 2006;Parmesan and Yohe 2003;Rödder et al. 2021;Thomas 2010).Undoubtedly, the species with narrower distribution ranges and strong reliance on host plants will be significantly more severe impacted compared to generalists (Bellard et al. 2012;Biella et al. 2017;Hoffmann et al. 2019;Mori et al. 2018;Thomas 2010).
The presence/absence of a species in a certain geographical and ecological space reflects a nexus between abiotic and biotic factors on the one hand and area(s) that species historically have been able to occupy on the other hand (Antonelli 2017;Brown and Carnaval 2019).The potential niche of a species is a geographical intersection between favorable abiotic (A) and biotic (B) factors for the species (Machado-Stredel et al. 2021;Milanesi et al. 2022;Soberon and Peterson 2005).However, due to some factors like dispersal ability, competition, and natural barriers, species would be able to occupy only some parts of the potential niche, named accessible area (Della Rocca and Milanesi 2020;Machado-Stredel et al. 2021).Therefore, any change in the abiotic and biotic factors can impact the habitat suitability of the species and consequently result in geographical shifts in the species range, particularly regarding endemic and rare species with narrow distributions (Bellard et al. 2012;Della Rocca and Milanesi 2020).Furthermore, species that rely heavily on specific interactions with other species, such as zygaenid moths, may be more vulnerable to these changes compared to generalist species.
Species of the family Zygaenidae, particularly endemic ones, spend the majority of their developmental time on a few plant species of the same genus or family (Hofmann and Tremewan 1996;Naumann et al. 1999).Two zygaenid subfamilies, the Zygaeninae and Procridinae, are mainly distributed across the Palearctic and their species can usually be found in open biotopes, areas exposed to the sun and light forests (Naumann et al. 1999).Although some species have a wide distribution range (across the whole Palearctic), most zygaenids have a regional or local distribution, and are extremely restricted to a few adjacent sites or are even unilocal (Hofmann and Tremewan 2017).Their larvae are mono-or oligophagous, their occurrence therefore strongly dependent on their host plant (Naumann et al. 1999).Consequently, future climate change can have a double effect on biodiversity, either directly by changing the habitat suitability of the species, or indirectly by affecting their interactions with their host plants.
With more than 46% of endemism, Iran is a unique hotspot for Zygaenidae diversity across the entire Palearctic (Hofmann and Tremewan 2017;Rajaei et al. 2023a).Iran is characterized by high landscape heterogeneity and a steep climatic gradient and it has been suggested as a transitional zoogeographical region between Palearctic, Oriental, and Saharo-Arabian realms (Holt et al. 2013;Rueda et al. 2013;Yusefi et al. 2021;Fig. 1).Several mountain ranges surround the high-elevation plateau of the country in the southwest of Asia (e.g., Zagros, Alborz, Kopet-Dagh), which provide a wide range of habitat for its biodiversity.Up to now, 73 species of the family Zygaenidae are known from the country, sorted in two subfamilies: Procridinae with 32 species and Zygaeninae with 41 species (Keil 2014;Rajaei et al. 2023a).The Iranian Zygaenidae are one of the best studied groups of non-Papilionoid, particularly over the past 50 years mainly by European zygaenologists (Rajaei et al. 2023b).Most of the species of this family are distributed across the mountainous areas in the north and west of the country, along two main mountain ranges (Mts.),Zagros and Alborz (Keil 2014;Rajaei et al. 2023a).These mountain ranges have historically served as barriers and corridors for gene flow and provided many species with microhabitats and glacial refugia (Ghaedi et al. 2021;Noroozi et al. 2018;Sanmartín 2003).Distribution analyses confirmed that 86% of all the Zygaenidae species are distributed across these two mountain ranges (Hofmann and Tremewan 2017;Keil 2014).Furthermore, more than 60% of the endemic species have been recorded in the Zagros Mts. and especially the central regions of this mountain range (Hofmann and Tremewan 2017;Keil 2014).Besides these mountain ranges, isolated mountains on the sidelines of the central basin provide suitable habitats for several endemic species (e.g.Hofmann and Tremewan 2017;Keil 2014).
Despite the long-term study of the biology and ecology of Zygaenidae species in Iran, there is a big knowledge gap concerning their taxonomy, ecology, distribution patterns, and Fig. 1 Species richness of Iranian Zygaenidae.The map shows distribution of endemic species (red dots), and non-endemic species (gray dots) occurrences across the country at 1-degree geographical unit (≈ 12,400 km 2 ) their conservation status.Therefore, in this study, we aim to (1) shed light on the ecology and the species distribution patterns of endemic and non-endemic species of Zygaenidae and their host plants across Iran.(2) explore potential habitat change in response to climatic changes for endemic and non-endemic species of Zygaenidae using the simulated distribution of species and their host plants under different optimistic and pessimistic climatic scenarios at the end of the current century.To achieve our objectives, we used the most comprehensive existing dataset for the group.We expect that narrowly distributed species suffer more from climate change and particularly global warming than widely distributed, non-endemic species.The outcomes of our study would provide a resource for directing conservation practices toward areas with high priority for conservation under future climate change.

Study area
Iran is the 18th largest country worldwide and is located in southwestern Asia between 25-40° north and 44-64° east (Fig. 1).The country spans across three different macrobioclimatic regions: the Mediterranean (in center and north), the Tropical (in south), and a small part of the Euro-Siberian regions (along the southern coastline of the Caspian Sea; Djamali et al. 2011).In general, the climate of Iran can be considered arid to semi-arid with low annual precipitation (~ 250 mm; Ghasemi and Khalili 2008;Madani 2014).There is a steep gradient of temperature depending on the location between − 3 and 7 °C in the coldest month to 29 to 37 °C for the warmest month (https:// www.irimo.ir/ en/ clima te).

Occurrence data
The checklist of the Zygaenidae of Iran includes the genera Rhagades (2 species), Zygaenoprocris (17 species), Adscita (one), Jordanita (12 species) and Zygaena (41 species; Rajaei et al. 2023a).To generate the dataset of occurrence data for all the Zygaenidae species, all data was gathered in a comprehensive literature review (e.g., Hofmann andTremewan 2017, 2020a, b;Keil 2014;Fig. 1).Additionally, Axel Hofmann provided additional occurrence data for distribution of the endemic and non-endemic species of Zygaenidae worldwide.The locations of all collected occurrences were carefully georeferenced using the software Google Earth Pro (v. 7.3.6.9345;Noori et al. 2023).The final dataset covered all 73 species of the family including more than 2500 occurrences, of which 1710 records remained after removing duplicated and missing values (with minimum 1, mean 23, and maximum 117 occurrences per species).
For this study, we selected 18 species (Table 1), from which 12 are endemic to Iran and the remaining are widely distributed in Central Asia, the Middle East, and the Palearctic.The species were selected if (a) there were more than 10 records for the species in our dataset.(b) The host plant(s) of the species is known, and we were able to gather more than 20 occurrences data for the host plant.(c) The biology of species is well-documented to test the factor of feeding mode (seven species are monophagous, and the remaining 11 species are oligophagous).Furthermore, we used the dataset of the host plants of the examined species, including 10 plant species of the families Apiaceae, Fabaceae, Polygonaceae, and Rosaceae.Two out of these 10 species are endemic to Iran (namely: Prunus eburnea (Spach) Aitch.and Ferulago carduchorum Boiss and Hausskn; see Table 1).The coordinate dataset host plants were gathered from different published datasets, i.e.Flora Iranica (Rechinger 1963(Rechinger -2015)), Revision of genus Eryngium (Wörz 2011), and GBIF database (Global Biodiversity Information Facility; www.gbif.org; the list of host plants and number of occurrences is provided in Table 1; see the reference for occurrence dataset of host plants from GBIF: Table S1).The dataset includes occurrence data for species of host plant across their ranges worldwide.

Environmental variables
We obtained climate variables from the CHELSA dataset (Climatologies at high resolution for the earth's land surface areas, version 2.1) for both presence  and future (2071-2100; https:// chelsa-clima te.org).The dataset includes high-resolution raster files with a spatial resolution of 30 arc sec (WGS84) of downscaled model outputs for different parameters of temperature and precipitation globally (Karger et al. 2017).The future variables of the CHELSA dataset include three out of five simulated socio-economic scenarios of CMIP6 (the international Coupled Model Intercomparison Project 6) at roughly 1 km resolution from the National Oceanic and Atmospheric Administration (from GFDL-ESM4 model; Karger et al. 2017).These scenarios (hereafter: climate scenarios) have been modeled based on increases in temperature and the concentration of CO 2 and greenhouse gases such as methane (Table 2; Karger et al. 2017).Table 2 shows a short description of the climate scenarios we applied in this study (Spp126, Spp370, and Spp585; for more details, see https:// chelsa-clima te.org).
To compare the current and future distribution of the species, we selected 19 bioclimatic variables (including different climatological variables for temperature and precipitation; see https:// chelsa-clima te.org) for the current and three future climate scenarios.We selected a subset of the variables avoiding any high multicollinearity using pairwise Pearson's correlation coefficients (r > 0.75).We also tested multicollinearity using VIF (Variance Inflation Factor) from the usdm package in R (Naimi et al. 2014;R Core Team 2022; see Supplementary Information I, Fig. S1).Furthermore, a Principal Component Analysis (PCA) was performed to explore the contribution of environmental variables to the principal components (ade4 package; Dray and Dufour 2007; Fig. S1).Finally, we used five of the variables to calibrate our models; these variables mainly represent the extreme levels of temperature and precipitation (for more details, see https:// chelsa-clima te.org), as: bio5 (mean daily maximum air temperature of the warmest month), bio6 (mean daily minimum air temperature of the coldest month), bio7 (annual range of air temperature), bio13 (precipitation amount of the wettest month), and bio14 (precipitation amount of the driest month).

Data preparation
Most of the analyses were run in the R environment (version 4.2.1;R Core Team 2022).Using the raster and biomod2 packages, a Presence-Absence Matrix (PAM) for both Zygaenidae and their host plants were generated (Hijmans 2022;Thuiller et al. 2021).Depending on the number of occurrences, we generated 15 times pseudo-absences for each species within a 500 km buffer around the occurrences to consider the enough environmental space and generate informative pseudo-absences (Barbet-Massin et al. 2012).The PAM includes binary data for the presence/absence of the species (1 and 0, respectively) and extracted values from environmental variables for current and future climate scenarios.The PAMs were created for further analysis in the ecospat and biomod2 packages (Broennimann et al. 2022;Thuiller et al. 2021).The ecospat package includes a set of comprehensive functions to study species distribution, and niche qualification and comparison (Di Cola et al. 2017).

Tuning the models
The following setup was applied to tune the computing function in biomod2.Modeling (biomod2) and ecospat.ESM.Modeling (ecospat).The functions were run for 10 replicates each; 80% of the occurrences were put aside to test the model.Table 3  values for evaluation of the model performance for each species based on Kappa (also known as Cohen's kappa coefficient), TSS (True Skills Statistic), and AUC (Area under the ROC Curve) metrics.To ensemble final models, we set a threshold for the models higher than AUC/ROC 0.7.Furthermore, models with higher values for TSS (> 0.7) were selected for projecting the species distribution to future climate scenarios (Table 3).Finally, we report the contribution of environmental variables and host plants for each species of zygaenid moth for both SDM and ESMs approaches (for more details see SI. section III).
To avoid this, we used the Ensemble Small Models (ESMs) approach in the ecospat package (Broennimann et al. 2022).The results of previous works demonstrated that the ESMs could reduce overfitting of the models regarding species with few occurrences (e.g., Breiner et al. 2018;Della Rocca and Milanesi 2022;Herrera et al. 2022).We computed the SDMs for the species with more than 30 unduplicated occurrences using the biomod2 package (Thuiller et al. 2021; Table 1).We executed our SDM and ESMs using maximum entropy (MaxEnt) machine learning algorithms, as they have the best performance, particularly for presence-only data (Elith et al. 2011;Fourcade et al. 2014;Phillips and Miroslav 2008;Phillips et al. 2006).The model performance was evaluated using different metrics such as True Skills Statistic (TSS), Area under the Receiver Operating Characteristic Curve (AUC/ROC), and Boyce Index.The latter has been designed to evaluate the performance of the models regarding presence-only data in the ecospat package (Di Cola et al. 2017; Table 3).

Host plants as abiotic variable
To import the host plant distribution in our model as a predictor, we first modeled the distribution of each host plant as it was described for zygaenids above and then used the result as predictor in combination with other abiotic variables in our models for each climate scenario (Table 2).Depending on the number of occurrences and distribution pattern we used SDM or ESMs to model the species distribution of the host plants and projected it to the future climate scenarios (Table 3).In case of host plants with a high number of occurrences (Securigera varia and Falcaria vulgaris), the coordinate dataset was thinned by two steps using the gridSample function from the dismo package, which thinned the dataset by selecting one record per pixel.Then the result of gridSample was used for further thinning steps using function thin from spThin (Aiello-Lammens et al. 2015).The thin function randomly keeps the occurrences with a user defined distance from each other (we considered a distance of 20 km as suitable).
Finally, to standardize the interpretation of the modeling products, we rescaled the resulting raster files from MaxEnt for habitat suitability.Raster files of habitat suitability were rescaled using omission of 10% of the lowest probabilities at the species records of SDM and ESMs predictions.In the next step, using the mess function in the dismo package, multivariate environmental similarity surfaces (MESS) were computed to assess the occurrence of extrapolation areas when projecting the models outside of their training range (Elith et al. 2010).Then values of MESS were used to evaluate accuracy of resulting habitat suitability in the future scenarios (the full projections of the species distribution under each climate scenario are shown in SI: Figs.S2-S30).In this study, we used a reclassified version of MESS to highlight those regions where at least one of the predictors exceed the training range (Rödder et al. 2013).The positive values of MESS were assigned to 0 and negative values to 1.This modification effectively characterizes the similarity or dissimilarity between the surveyed pixels and all the areas under study.

Area of habitat suitability
The predicted layer of habitat suitability can depict a wide range of areas as suitable habitats for a given species, including areas which are unlikely to be accessible.Therefore, to delimit the real suitable habitat we cropped the predicted layers for the current and future climate scenarios with a buffer of 50 km for a hull polygon of species occurrence using the rangemap package (Cobos et al. 2021;Soberon and Peterson 2005).We considered the area with the higher probability (> 25% and > 50%) within the buffer, as the area where species might be present in each climate-scenario.The area of the raster pixel with presence of the species were calculated in km 2 (Table 4).Furthermore, we only accepted the areas after subtracting the overlap of negative values for MESS which indicate extrapolation (Table 4).Table 4 depicts the area of current species range and the percentage of remaining areas compared with current species range under each future climate scenarios for values of habitat suitability more than 25% and 50%.
Table 4 The areas of current species range of zygaenid species and their host plants and the percentage of remaining areas under future climate scenarios (for higher probability values > 25% and > 50%). (H)

Species range shift
The PAM of species was generated from the raster values with the high probability for each species (> 50%).This PAM was used to look at the species' altitudinal preference for current and future climate scenarios, and to explore the overlap between species and their host plant.The elevation values were extracted for each simulated species' occurrence from the digital model for the earth elevation (Global Digital Elevation Model, ver.3; www.nasa.gov), then density graphs were generated with the ggplot2 package.Finally, we used the raster.overlapfunction from the ENMTools package to measure the overlap between the habitat suitability of moths and host plants (Warren et al. 2008;Warren and Dinnage 2023).The raster.overlap has been developed to measure the niche overlap resulting from species distribution modeling and has several metric values to explore the niche overlap.We used two metrics to explore the overlap between moths and their host plant, Schoener's similarity (D) and similarity statistic (I), of which both will result in a value between 0 (no overlap), and 1 (identical niche prediction).While the former has been used because of its simplicity and long-term use in biological interpretation, the latter, which is a modification of the Hellinger distances, is a measure of the similarity between two probability distributions and was developed to compare the community composition of different sites (Rödder and Engler 2011;Warren et al. 2008)

Results
Overall, considering the most probable habitats of zygaenid species (> 50%), our results show that more than 80% of the studied species will lose around 30%, 70%, and 75% of their habitat suitability under future climate scenarios (Spp126, Spp370, and Spp585) compared to their current distribution, respectively (Table 4; Fig. 2).Although under future climate scenarios some of the non-endemic species of zygaenids and their host plants will experience an expansion in their current species range e.g., Z. loti, F. vulgaris, S. varia, most endemic species and their host plants will lose a dramatic area of their ranges (Table 4, Fig. 2).

Species richness and endemism
Results of the current study suggest that most of the Zygaenidae species in Iran are distributed across mountainous areas in the north and the western half of the country (Fig. 1).However, as evident in Fig. 1 richness of endemic species is more pronounced along the southern parts of Zagros Mts., and mountainous regions of Kerman province in the south.Furthermore, central, and eastern parts of Alborz Mts. and the southern areas of Kopet-Dagh Mts. are other regions with a high number of endemic species.Ghohrud Mts., a chain of segregated high-elevation mountains along the western margin of the Central Basin and parallel to the Zagros Mts., is another hotspot of endemism (Fig. 1).

Habitat suitability under climate change
Results of this study demonstrate that most Zygaenidae species and their host plants will, at least to some extent, experience a shrinking in their current distribution range under both optimistic and pessimistic climate scenarios by the end of this century (Table 4; Fig. 2).The degree of habitat loss does show a significant correlation with endemism.There is a distinct difference between the response of endemic and non-endemic species to climate scenarios.As shown in Table 4 and Fig. 2, all non-endemic species will experience a significantly smaller reduction in their species range under each climate scenario.The only exception is Z. haematina, which is not endemic for Iran but mainly distributed from the southeast Turkey toward the center of Iran (endemic in Zagros Mts.).Except for one of the endemic species (Zygaena pseudorubicundus), other endemic species show a dramatic decline in their ranges.Unlike other endemic species, the rate of habitat loss in Z. pseudorubicundus is not significant, and as shown in Table 4 and Figs. 2 and 3, the area of habitat suitability for the species is shrinking by 40% even under pessimistic (Spp585) climate scenarios.However, the rate of habitat loss for the species in the eastern and southern part of its range (Zagros Mts.) is much higher than the northern distribution across the central Alborz Mts.Our model predicts that the Zagros's population of Z. pseudorubicundus might shrink significantly under pessimistic climate scenarios (Table 4, Fig. 3).
In most areas, habitat loss is much more severe under the pessimistic scenarios (Spp370 and Spp585; Table 4).The models predict the complete vanishing of habitat suitability for species like Zygaena aisha, Z. mirzayansi, and Z. ginnereissi, considering area with higher probability for species distribution (> 50%; Table 4; Figs. 2, 4).
Additionally, the suitable habitat area for species like Z. nocturna, Z. fredi, Z. kermanensis, Z. haematina, and Rhagades brandti, is getting dramatically smaller than the current species range (< 0.1%, < 2%, < 2.9%, and 4%, respectively) under pessimistic scenarios Fig. 2 A comparison between the area of habitat suitability for the Zygaenidae species in current and future climate scenarios (Spp126, Spp370, Spp585).The gradient of color from green to dark red is comparable with the rate of habitat loss, which is significantly higher for endemic and narrow-distributed species compared with non-endemic species.(*) indicates the endemic species (Spp585).However, the rate of habitat loss is lower for non-endemic species, like Z. loti (< 25%; Table 4; Fig. 2).
A bit less pronounced, we detected a similar trend of species-range shifts for the host plants (Table 4).As already discussed for the zygaenid moths, non-endemic species of host plants with wider distribution range will be affected less and even experience an expansion in their ranges e.g., Securigera varia, Falcaria vulgaris.However, the rate of habitat loss even in non-endemic species with narrower species range such as Onobrychis cornuta, Eryngium billardieri is higher (Table 4).The rate of habitat loss for the host plant species is much higher in southern parts of the country compared to its northern areas.
Our results suggest a higher impact of climate change, and particularly global warming, on the southern and central regions of the country.Zygaenoprocris duskei, a monophagous species, has a narrow distribution from the center to the southeast of Iran, and its larvae feed only on the Atraphaxis spinosa (Fig. 5).Under the pessimistic climate scenarios (Spp370, and Spp585), the host plant will lose most of its range in southern Iran and will shift towards higher latitudes, which will result in habitat loss and habitat fragmentation of Z. duskei (Fig. 5).
There are two general trends visible in our results: (a) shifting in the ranges towards higher latitude for mainly non-endemic species, and (b) shifting of species range towards higher elevation in the endemic species (Figs. 6,7,8).Zygaena loti, which lives on the Securigera varia as its host plant (in Iran), is mainly distributed in the Alborz Mts., the Fig. 3 Forecast of endemic species, Zygaena pseudorubicundus and its host plant Falcaria vulgaris for current and under three climate scenarios (Spp126, Spp370, Spp585) by the end of the century in Iran.The yellow points represent occurrences of the species in our dataset; Gradient of red represents habitat suitability for the moth (Z.pseudorubicundus), and gradient of green represents habitat suitability for the host plant (F.vulgaris), the intensity of the color depicts the probability of species presence.Areas of potential extrapolations (MESS) are indicated as grey shading for both moth and host plant Caucasus, the Transcaucasus, Turkey, and southeastern Europe (Fig. 6).In response to the different climate scenarios, the area of habitat suitability of Z. loti will be stable (< 25% reduction under Spp585).Although the species will move towards higher latitudes in general, the population in Alborz Mts. will shift to higher elevations.The same distribution pattern can also be detected in distribution ranges of other non-endemic species like Zygaena tamara, and Zygaena araxis.
On the other hand, the species with narrower distribution range will experience a dramatic decline in their current range under all the future climate scenarios (Figs. 4,5,7,8).The endemic species, Zygaena ecki is one of those species with a small distribution range in the central and eastern part of Alborz Mts.(Figs. 7, 8d).Under pessimistic climate scenarios the species will lose most of its eastern distribution, particularly at the high elevation of Shahkuh Mountain in Semnan Province (> 84% reduction).A similar trend is observed for Onobrychis cornuta, the host plant of Z. ecki.(Fig. 7).
Moving toward higher elevation is more pronounced for endemic species e.g., in Z. kermanensis, and Z. nocturna (Fig. 8e, f), while species with wider distribution across central Asia to western Europe, do not show any significant altitudinal shifts in their species range, e.g., Z. loti (Fig. 8c).The results of the overlapping of habitat suitability show no significant gap in species range between zygaenid species and their host plants.
Results of model evaluation revealed that precipitation amount of the wettest month (bio13) has highest contribution in modeling distribution for both zygaenid moths and their Fig. 4 Forecast of endemic species, Zygaena aisha and its host plant Ferulago carduchorum for current and under three climate scenarios (Spp126, Spp370, Spp585) by the end of the century in Iran.The yellow points represent occurrences of the moth species in our dataset; Gradient of red represents habitat suitability the moth (Z.aisha), and gradient of green represents habitat suitability for the host plant (F.carduchorum), the intensity of the color the probability of species presence.Areas of potential extrapolations (MESS) are indicated as grey shading for both moth and host plant host plants (SI.Section III).On the other hand, while bio7 (annual range of air temperature) and distribution of host plants play an important role on prediction of species range of moths, daily minimum air temperature of the coldest month (bio6), has a greater influence on the species distribution of the host plants (SI.Section III).

Discussion
Shifting towards higher elevation has been documented as a common response of different insect taxa and their host plants to climate change word wide (Biella et al. 2017;Della Rocca and Milanesi 2022;Filazzola et al. 2020;Rödder et al. 2021;Pyke et al. 2016).For instance, Rödder et al. (2021) revealed a constant altitudinal shift of species range for several butterflies in the eastern Alps during the past six decades.According to the results of the present study, zygaenid species of Iran will generally experience altitudinal range shift and a high habitat loss (> 64%) under the most extreme climate scenarios.However, the rate of habitat loss is twice as high for endemic species compared with non-endemics (Table 4; Fig. 2).This might be explained by the fact that non-endemic species with their wider distribution have access to a wider range of habitats and host plants compared to the endemics (Biella et al. 2017;Filazzola et al. 2020;Rödder et al. 2021;Pyke et al. 2016).As example, Zygaena aisha and Z. ginnereissi have an extremely restricted distribution Fig. 5 Forecast of endemic species, Zygaenoprocris duskei and its host plant Atraphaxis spinosa for current and under three climate scenarios (Spp126, Spp370, Spp585) by the end of the century in Iran.The yellow points represent occurrences of the moth species in our dataset; Gradient of red represents habitat suitability for the moth (Z.duskei), and gradient of green represents habitat suitability for the host plant (A.spinosa), the intensity of the color depicts the probability of species presence.Areas of potential extrapolations (MESS) are indicated as grey shading for both moth and host plant range across the higher mountains of Kerman in the southern part of the Iranian central basin (Fig. 4).This region includes some mountains with elevation higher than 4000 m surrounded by the central deserts.Furthermore, zygaenids are not strong and fast fliers and therefore are highly dependent on their habitat (Naumann et al. 1999).Therefore, the response of these species to climate change will be limited to a shift toward higher elevation (Biella et al. 2017;Della Rocca and Milanesi 2022;Filazzola et al. 2020;Rödder et al. 2021;Fig. 4).On the other hand, non-endemic species like Zygaena loti, have more opportunities to move toward the higher latitudes across the Caucasus and Transcaucasia regions under extreme climate scenarios (Filazzola et al. 2020;Rödder et al. 2021;Table 4; Fig. 5).

High species-richness across mountains ranges
Our analyses depict a strong association between species richness and endemism with high-elevation regions across the main mountain ranges in most parts of the country: Zagros Mts. and Ghohrud Mts., Alborz Mts., Kopet-Dagh Mts., Kerman-Yazd Massif and Makran-Taftan Mts., which highlights the effect of the complex topography on distribution pattern of the Zygaenidae species.These results are in line with the previous results of the independent studies, which highlight the important roles of the mountainous areas to shape the biodiversity in Iran (Ghaedi et al. 2021;Noori et al. 2021;Noroozi et al. 2018Noroozi et al. , 2019;;Yusefi et al. 2019).The above mentioned mountainous ranges increase the rate of Fig. 6 Forecast of non-endemic species, Zygaena loti and its host plant Securigera varia for current and under three climate scenarios (Spp126, Spp370, Spp585) by the end of the century in Iran.The yellow points represent occurrences of the moth species in our dataset; Gradient of red represents habitat suitability for the moth (Z.loti), and gradient of green represents habitat suitability for the host plant (S.varia), the intensity of the color depicts the probability of species presence.Areas of potential extrapolations (MESS) are indicated as grey shading for both moth and host plant isolation and at the same time provide a wide variety of microhabitats, which can act as refugia to buffer the effect of climate change for different species (e.g., Albrich et al. 2020;Della Rocca and Milanesi 2022;Djamali et al. 2012;Paknia and Rajaei 2015;Rajaei et al. 2013).Several studies suggested the dual effects of mountain ranges in Iran as corridors and simultaneously as a barrier to the speciation of different taxa (e.g., Ghaedi et al. 2021;Sanmartín 2003).The unique species composition of Zygaenidae in Iran, where many closely related species occur, suggests that Iran has played a significant role in the diversification of this family of moths (Hofmann and Tremewan 2017).Considering that most of above listed mountains fall in the Irano-Anatolian biodiversity hotspot, identifying speciesdiverse regions within this hotspot will help to delineating the areas with higher priority for conservation (Cañadas et al. 2014;Noroozi et al. 2018).

Heterogenous impact of climate change
Although climate change has impact on biodiversity at all levels, the risk of extinction is much higher for the species that occur in smaller and patchier habitats (Della Rocca and Milanesi 2022;Filazzola et al. 2020;Pardini et al. 2017;Rödder et al. 2021).Therefore, access of the species to larger and more diverse habitats may increase the ability of that species to tolerate climate change better (Filazzola et al. 2020;Franzén and Ranius 2004;Rödder et al. 2021).It has been well documented that species with more restricted ranges Fig. 7 Forecast of endemic species, Zygaena ecki and its host plant Onobrychis cornuta for current and under three climate scenarios (Spp126, Spp370, Spp585) by the end of the century in Iran.The yellow points represent occurrences of the moth species in our dataset; Gradient of red represents habitat suitability for the moth (Z.ecki) and gradient of green represents habitat suitability for the host plant (O.cornuta), the intensity of the color depicts the probability of species presence.Areas of potential extrapolations (MESS) are indicated as grey shading for both moth and host plant will suffer much more than species with a broad ecological amplitude from the rapid climate change (Bonelli et al. 2022;Rödder et al. 2021).Furthermore, species which are highly dependent on their host plants might experience a dual impact of climate change.Directly by effecting the species habitat suitability and indirectly by changing the interaction of species and host plants (Bellard et al. 2012;Blois et al. 2013;Filazzola et al. 2020).
In line with the previous studies, our results reveal a higher impact of climate change on the biodiversity in the center and south of the country than northern regions, which are the regions with higher species richness and pronounced endemism across mountainous areas (Figs. 1, 5;Ashrafzadeh et al. 2019;Shamsabad et al. 2018).Different studies suggested a higher impact of climate change in the Middle East and especially in Iran, because of the high level of contribution to the emission of greenhouse gas (Mansouri et al. 2019;Segan et al. 2016;Waha et al. 2017).As a general trend, our model predictions show that the impact of climate change even across a mountain range (e.g., Alborz Mts.) is not homogeneous.The eastern parts of the Alborz Mts. will be affected more than its western parts, which can be seen e.g., for Zygaena ecki (Fig. 7).This might be interpreted by heterogenous topology of the mountainous areas which provide a wide range of habitats with different climatic setups (Albrich et al. 2020;Djamali et al. 2012).Most of the endemic Zygaenidae species are distributed in the small zones above the tree line of mountains (Hofmann and Tremewan 2017;Keil 2014).To this end, different studies confirmed the important role of mountain ranges in the configuration of the country's biodiversity (Noori et al. 2021;Noroozi et al. 2018;Yousefi et al. 2019Yousefi et al. , 2023)).

Threats of Iranian Zygaenidae
Several studies have highlighted a significant gap between the current protected areas of the country and the most diverse regions for different groups of animals and plants, particularly in mountainous areas (Noori et al. 2021;Noroozi et al. 2019Noroozi et al. , 2023;;Yusefi et al. 2019).The network of protected areas of the country is under intense pressure by human activities (Karimi and Jones 2020).For instance, overgrazing has been reported as one of the most important threats to natural habitats and particularly high-elevation biodiversity in Iran (Karimi and Jones 2020;Jowkar et al. 2016).This factor has been suggested as one of the important factors which affect the habitat quality of the burnet moths (Franzén and Ranius 2004;Naumann et al. 1999).Franzén and Ranius (2004) suggested that the high correlation between distribution of burnet moths and butterflies, reflecting the fact they have similar habitat requirements.Therefore, simulation of the future distribution pattern of the species can help to understand the effect of climate change not only on zygaenid moths, but also other groups of Lepidoptera and maybe other insects (Bellard et al. 2012).Consequently, defining higher priority areas for the conservation of vulnerable groups like Zygaenidae under ongoing climate change is inevitable.This will help scientists and decision-makers to estimate the extinction risk of the different species by investing limited resources for highly protecting species-diverse areas efficiently and develop target habitat management plans for habitats that are particularly at risk.

Limitation of the model and potential enhancements
While the insights gained from this study are valuable, it is important to acknowledge the inherent limitations in our approach.As previously discussed, most of the endemic zygaenid moths are in small populations that are highly localized, rendering it impractical to gather a more extensive dataset (Hofmann and Tremewan 2017;Keil 2014;Naumann et al. 1999).This may increase the risk of model overfitting and the impact of spatial autocorrelation (SAC) on our analysis (Dormann et al. 2007).Our examination of SAC within the data reveals the chance for SAC for some of the species, particularly endemic one.Consequently, the results necessitate cautious consideration.To enhance the reliability of our findings, it is imperative to conduct more intensive surveys in the study area, thereby reducing the bias in sampling effort.Furthermore, incorporating other variables (i.e., land cover, topology, etc.) linked to the physiology and phenology of the zygaenid moths and their host plants could bolster the robustness of our results.

Conclusion
Rapid anthropogenic climate change impact on the current biodiversity and accelerates the risk of extinction higher than at any time on the planet earth (Pimm et al. 2014;Settele et al. 2016;Shivanna 2020).Climate change is not only affecting the habitat suitability of the species but also the interactions between species (Bellard et al. 2012;Blois et al. 2013).Estimating the reaction of different species to climate change may help to design more effective conservation strategies.Although the present study was limited to the species of the family Zygaenidae, it provides an example of how climate change will affect biodiversity unevenly at the level of a country.Our results show different responses of the endemic and non-endemic species to future scenarios of climate change.While nonendemic zygaenid species might move poleward, the endemic species may move towards higher elevations, especially due to their high dependence on their host plants/habitats and low flying ability.Our models predicted that higher mountains in the southern and central parts of Iran may be affected more severely than higher latitudes.Considering that the mountainous areas with high biodiversity are under high pressure from human activities due to being close the populated cities, expansion of the current network of protected areas toward regions with higher species diversity is an inevitable solution.However, designing an effective conservation practice depends on improving our understanding regarding the distribution pattern of different species, particularly the mega-diverse group of insects.

Table 1
List of the examined Zygaenidae species and their host plants, with endemism status and number of occurrences for both moth and its host plants after omitting duplicated, missing values, and thinning the coordinate dataset

Table 2
List of the future climate scenarios and short description of their features

Table 3
Results of model evaluation for the studied species and(M)stand for Host plant and Moth, respectively.The table is sorted alphabetically based on species names of moths and host plants