Response of Iranian lizards to future climate change by poleward expansion, southern contraction, and elevation shifts

This study explores the relationships between recent Iranian lizard species distributions and the observed climate, as well as potential future distributions of species. For this purpose, an ensemble of seven algorithms was used to forecast the distributions of 30 species for the recent and future (2070) based on the averages of 14 global climate models under optimistic (RCP2.6) and pessimistic (RCP8.5) scenarios. Annual precipitation (n = 16) and annual mean temperature (n = 7) were identified as the most important variables in determining the distribution of 76.66% (23 out of 30) of the species. The consensus model predicts that the ranges of 83.33% of species (n = 25) have the potential to expand poleward at higher latitudes while preserving the majority of their recent distributions (except for four species). Furthermore, the ranges of the remaining species (n = 5) will be preserved at higher latitudes. However, they (n = 22) may contract slightly (n = 13) or excessively (n = 9) in the south of their distribution range at lower latitudes. These results indicate that species (N = 19) situated in mountainous areas such as the Zagros, Alborz, and Kopet Dagh may move or maintain their range at higher elevations as a result of future climate change. Finally, this study suggests that 30% of species (n = 9) may be threatened by future climate change and that they should be prioritized in conservation efforts.


Results
The ensemble models exhibited high quality, with TSS, AUC, and KAPPA values ranging from 0.76 to 1 ( Table 1). The map of suitable habitats based on elevation and bioclimatic variables for recent and future (2070) climate conditions under RCP2.6 and RCP8.5 scenarios for 30 species of lizards that are distributed in Iran is shown in Fig. 2. The mean of variable importance (%) as estimated by the algorithms for the 30 species of lizards are provided in Table 2. For 16 species BIO12, for seven species BIO1, for three species BIO15, for two species BIO5 and one species elevation were found to be the most important variables affecting the distribution of species. Table 3 and Fig. 3 illustrate the range shift of 30 species of lizards in recently suitable habitats (gain/loss) by 2070 under the RCP2.6 and RCP8.5 scenarios.
Future climate change will affect 30 Iranian lizards in different ways, some by expanding their ranges, some by contracting their ranges, and others by remaining relatively unaffected (especially in habitats with 75-100% suitability) (Tables 3 and 4  www.nature.com/scientificreports/ A. pannonicus (AP). BIO12 (49.76%) and BIO1 (12.86%) are the two important variables affecting the distribution of Asian Snake-eyed Skink, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the west and southwest, along with the Zagros, Alborz, and Kopet Dagh ranges, as well as a portion of the northwest and northeast, have 50 to 100% suitability. Habitats with 75 to 100% suitability are concentrated in the middle and southern Zagros, a narrow strip of Alborz, and a small section of northeastern Kopet Dagh, as well as the northwestern and northeastern parts of the country (Fig. 2 AP1). Based on future climate change (especially RCP8.5), suitable habitats (especially 75-100% suitability), in addition to maintaining the recent range, will expand to higher latitudes in the north and northwest. While the eastern and southern margins of the distribution range will contract in the face of future climate change ( Fig. 2 AP2 and AP3 and Fig. 3 AP1 and AP2). Expansion to higher elevations for habitats with 75-100% suitability is also apparent in the Zagros, Alborz, and Kopet Dagh mountains ( Fig. 2 AP2 Table 2). The habitat suitability map for recent climate conditions shows that the west and southwest, along with the Zagros, have 50 to 100% (especially 75-100%) suitability ( Fig. 2 AE1). Based on future climate change (especially RCP8.5), suitable habitats (especially 50-75% suitability), in addition to maintaining the recent range, will expand to higher latitudes in the north and northwest of Iran ( Fig. 2 AE2 and AE3 and Fig. 3 AE1 and AE2). Expansion to higher elevations for habitats with 75-100% suitability is also apparent in the Zagros mountains ( Fig. 2 AE2  www.nature.com/scientificreports/ mountains ( Fig. 2 EA2 and EA3). According to RCP8.5, habitat loss will occur at low elevations in the western and southern margins of the Zagros mountains ( Fig. 3 EA1 and EA2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 0.10% and 4.61%, and new habitat gain at 30.78% and 68.58%, respectively (Table 3).
E. schneideri (EUS). BIO12 (55.41%) and BIO4 (12.61%) are the two important variables affecting the distribution of Schneider's Long-legged Skink, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that most of the northwest, west, and southwest, along with the Zagros range, Alborz and Kopet Dagh mountains, the northeast, and a small portion of southeastern Iran have 50 to 100% suitability ( Fig. 2 EUS1). The deserts of central and northeast Iran are devoid of EUS. Based on future climate change (especially RCP8.5), suitable habitats (especially 75-100% suitability), in addition to maintaining the recent range, will expand to higher latitudes ( Fig. 2 EUS2 and EUS3 and Fig. 3 EUS1 and EUS2). Expansion to higher elevations for habitats with 75-100% suitability is also apparent in the Zagros, Alborz, and Kopet Dagh mountains ( Fig. 2 EUS2 and EUS3). According to both scenarios, especially RCP8.5, habitat loss will occur on the eastern margins of the Zagros mountains, especially toward the south, as well as a small portion of the northeast and southeast (except RCP2.6) of the distribution range in the low latitudes ( Fig. 3 EUS1 Table 2). The habitat suitability map for recent climate conditions shows that the Zagros mountains, south, southeast, northeast and center of Iran have 50 to 100% suitability ( Fig. 2 LN1). Habitats with 75 to 100% suitability are concentrated in the west and southwest, along with the Zagros mountains ( Fig. 2 LN1). Based on future climate change (especially RCP8.5), suitable habitats (especially 50-75% suitability), in addition to maintaining the recent range, will expand to higher latitudes. Expansion to higher elevations for habitats with 50-75% suitability is also apparent in the Zagros mountains ( Fig. 2 LN2 and LN3 and Fig. 3 LN1 and LN2). According to both scenarios, especially RCP8.5, habitat loss will occur on the eastern margins of the Zagros mountains, as well as the south and southeast of the distribution range in the low latitudes ( Fig. 3 LN1 Table 2). The habitat suitability map for recent climate conditions shows that the west and southwest along the Zagros mountains, as well as a small portion of the Mesopotamian plain, have 50 to 100% (especially 75-100%) suitability (Fig. 2 MIH1). Based on future climate change (especially RCP8.5) suitable habitats in addition to maintaining the recent range will expand to higher latitudes ( Fig. 2 MIH2 and MIH3 and Fig. 3 MIH1 and MIH2). Expansion to higher elevations for habitats with 75-100% suitability is also apparent in the Zagros mountains ( Fig. 2 MIH2 and MIH3). According to RCP8.5, habitat loss will occur in the western margins of the Zagros mountains in low elevation, as well as southern margins at low latitudes ( Fig. 3 MIH1 and MIH2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 1.85% and 11.15% and new habitat gain at 47.78% and 95.07%, respectively (Table 3).  www.nature.com/scientificreports/  www.nature.com/scientificreports/ mountains in higher latitudes (Fig. 2 ML2 and ML3 and Fig. 3 ML1 and ML2). According to both scenarios, especially RCP8.5, habitat loss will occur in the eastern and southern parts of the distribution range, especially at low latitudes (Fig. 3 ML1 and ML2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 58.74% and 92.84%, and new habitat gain at 60.15% and 122.23%, respectively (Table 3).

M. persicus (MP)
. BIO15 (50.85%) and BIO1 (29.24%) are the two important variables affecting the distribution of Persian Tiny Gecko, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the southwest, south, and southeast of Iran have 50 to 100% (especially 75-100%) suitability (Fig. 2  HF1). Based on future climate change (especially RCP8.5), suitable habitats, in addition to maintaining the recent range, will expand to higher latitudes in the northern margins of the southern and southwestern ranges. Expansion onto the Mesopotamian plain is also apparent, especially in RCP8.5 (Fig. 2 MP2 and MP3 and Fig. 3  MP1 and MP2). According to RCP8.5, habitat loss will occur in the northern margins of the southeastern range ( Fig. 3 MP1 and MP2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 9.20% and 5.23%, and new habitat gain at 2.47% and 18.31%, respectively (Table 3).
O. elegans (OE). BIO12 (39.50%) and BIO14 (17.69%) are the two important variables affecting the distribution of Elegant Snake-eyed Lizard, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the northwestern, western, and southwestern parts of the country around and along the Zagros range, the Mesopotamian plain, the southern Alborz, a part of the Kopet Dagh mountains, and the southern Iranian plateau have 50 to 100% suitability (Fig. 2 OE1). Based on future climate change in the RCP2.6, suitable habitats (especially 75-100%), in addition to maintaining the recent range, will expand to the northwestern  www.nature.com/scientificreports/ and northern regions of the high latitudes ( Fig. 2 OE2 and Fig. 3 OE1). According to the RCP8.5 scenario, the western margin of the northern distribution, the western and southwestern, especially the Mesopotamian plain, the southern and eastern margins of the distribution area around and along the Zagros Mts., and the southern Iranian plateau at lower latitudes will be lost ( Fig. 2 OE3 and Fig. 3 OE2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 33.20% and 78.89%, and new habitat gain at 66.52% and 51.85%, respectively (Table 3).

P. caucasia (PC)
. BIO1 (32.99%) and BIO5 (25.07%) are the two important variables affecting the distribution of Caucasian Agama, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the Kopet Dagh and Alborz ranges, northwestern of the country towards the central and southwestern Zagros, as well as eastern Iran, have 50 to 100% (especially 75-100%) suitability (Fig. 2 PC1). According to RCP2.6, habitats with 75-100% suitability will be lost at low latitudes, whereas these habitats will remain at higher latitudes, especially in the Kopet Dagh, Alborz, and Zagros highlands and northwest of Iran (Fig. 2 PC2,  Fig. 3 PC1). According to the RCP8.5, the habitats with 50-100% suitability will remain at higher latitudes, especially in Kopet Dagh, the Alborz highlands, northwest, and as well as a small portion of the Zagros highlands (Fig. 2 PC3). The range of species shift in the RCP8.5 scenario reveals a significant loss in the entire range of species distribution (Fig. 3 PC2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 59.32% and 94.67%, and new habitat gain at 9.25% and 4.03%, respectively (Table 3).

P. maculatus (PM)
. BIO12 (36.73%) and BIO5 (17.42%) are the two important variables affecting the distribution of Spotted Toad-headed Agama, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that most of the central Iranian plateau has 50 to 100% suitability (Fig. 2 PM1). Based on future climate change (especially RCP8.5), suitable habitats (especially 75 to 100% suitability), in addition to maintaining the recent range, will expand in higher latitudes (Fig. 2 PM2 and PM3 and Fig. 3 PM1 and PM2). According to both scenarios, especially RCP8.5, habitat loss will occur in the margins of the distribution range, especially in the southeastern parts at low latitudes (Fig. 3 PM1 and PM2) Table 2). The habitat suitability map for recent climate conditions shows that the northwest country towards the central and small portion of the southwestern Zagros mountains, as well as the Kopet Dagh and Alborz ranges, have 50 to 100% suitability (Fig. 2 PP1). Based on future climate change, suitable habitats (especially 75 to 100% suitability), in addition to maintaining the recent range, will expand into higher latitudes (Fig. 2 PP2 and PP3 and Fig. 3 PP1 and PP2). According to both scenarios, especially RCP8.5, habitat loss will occur in the southern margins of the distribution range at low latitudes ( Fig. 3 PP1 and PP2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 10.91% and 37.32%, and new habitat gain at 87.11% and 90.96%, respectively (Table 3).

No. Species Suitability classes (%)
Habitat area suitability (%)/ its coverage (%) within PAs  www.nature.com/scientificreports/ P. scutellatus (PS). BIO12 (44.04%) and BIO1 (22.49%) are the two important variables affecting the distribution of Gray Toad-headed Agama, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that a wide range of the central plateau has 50 to 100% suitability (Fig. 2 PS1). In RCP8.5, habitats with 75-100% will contract as a result of future climate change, whereas habitats with 50-75% suitability will expand, especially at higher latitudes (Fig. 2 PS3). According to both scenarios, especially RCP8.5, new habitats will be gained in the Alborz, Zagros, south, and east highlands, whereas habitat loss will occur at lower elevations in the northern regions and will be restricted to low latitudes in the southern regions (Fig. 3 PS1 and PS2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 59.74% and 88.10%, and new habitat gain at 28.89% and 17.13%, respectively (Table 3).
T. caspius (TC). BIO1 (30.36%) and elevation (28.37%) are the two important variables affecting the distribution of Caspian Thin-toed Gecko, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the west of the Caspian Sea, Hyrcanian forests, the southern range of the Alborz, east of the Alborz and Kopet Dagh mountains, east and also a small part of the central and southern Zagros highlands have 50 to 100% suitability (Fig. 2 TC1). Habitats with 75 to 100% suitability are concentrated in the west of the Caspian Sea, Hyrcanian forests, the southern range of the Alborz, and northeast of Iran (Fig. 2 TC1). Future climate change will significantly decrease the extent of habitats with the suitability of 75-100%, especially in the RCP8.5 scenario, and they will be shifted to the heights of the Alborz and Kopet Dagh ranges (Fig. 2 TC2 and  TC2). According to both scenarios, especially RCP8.5, habitat loss will occur at lower elevations in the northern regions (RCP8.5) and will be restricted to low latitudes (both scenarios) in the southern regions (Fig. 3 TC1  and TC2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 39.90% and 79.30% and new habitat gain at 16.07% and 10.26%, respectively (Table 3).
T. princeps (TP). BIO12 (53.03%) and BIO15 (14.41%) are the two important variables affecting the distribution of Prince Lacerta, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the west and southwest along the Zagros mountains, as well as the Mesopotamian plain, have 50 to 100% (especially 75-100%) suitability (Fig. 2 TP1). According to both scenarios, especially RCP8.5, future climate change will decrease the extent of habitat with 75-100% suitability, especially at lower latitudes ( Fig. 2 TP2 and  TP3). Under both scenarios, especially RCP8.5, habitat loss will occur at the western, southern, and eastern margins of the species distribution range in the Zagros mountains, but new habitats will be gained in the north of the distribution range and at higher latitudes ( Fig. 3 TP1 and TP2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 14.55% and 29.83%, and new habitat gain at 2.63% and 4.15%, respectively (Table 3).
T. septemtaeniata (TS). BIO12 (57.67%) and BIO5 (16.90%) are the two important variables affecting the distribution of Southern Grass Skink, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the northwest, west, and southwest, along with the Zagros mountains, the Mesopotamian plain, north, northeast, along with the Kopet Dagh ranges, and east have 50 to 100% suitability. Habitats with 75 to 100% suitability are concentrated in the middle and southern Zagros mountains (Fig. 2 TS1). Based on future climate change (especially RCP8.5), suitable habitats (especially 75-100% suitability), in addition to maintaining the recent range, will expand to higher latitudes in the northwest, north, northeast, and east of the country (Fig. 2 TS2 and TS3 and Fig. 3 TS1 and TS2). Expansion to higher elevations for habitats with 50-100% suitability is also apparent in the Zagros and Kopet Dagh mountains (Fig. 2 TS2 and TS3). Under both scenarios, especially RCP8.5, habitat loss will occur in the Mesopotamian plains and the southern margins of the Zagros mountains at lower latitudes, whereas new habitats will be gained at higher latitudes ( Fig. 3 TS1 and TS2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 20.19% and 20.77%, and new habitat gain at 58.73% and 107.96%, respectively (Table 3).
T. agilis (TA). BIO1 (24.36%) and BIO4 (21.08%) are the two important variables affecting the distribution of Agile Ground Agama, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that almost all of Iran, except the northwestern part of the Zagros, has 25 to 100% suitability (Fig. 2 TA1). Based on future climate change (especially RCP8.5), suitable habitats, in addition to maintaining the recent range, will expand to higher latitudes (Fig. 2 TA2 and TA3 and Fig. 3 TA1 and TA2). Northwest of the species' distribution range at higher latitudes, especially according to the RCP8.5, may also become habitats with the suitability of 50-75% (Fig. 2 TA2 and TA3). Habitat loss is apparent in the east, south, and southeast of the species distribution range in Iran (Fig. 3 TA1 and TA2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 31.87% and 34.17%, and new habitat gain at 72.15% and 143.79%, respectively (Table 3).
T. ruderatus (TR). BIO12 (58.52%) and BIO5 (12.02%) are the two important variables affecting the distribution of Horny-scaled Ground Agama, respectively ( Table 2). The habitat suitability map for recent climate conditions shows that the northwestern, western, and southwestern parts of the country, along with the Zagros range, Mesopotamian plain, and northeast of Iran, have 50 to 100% suitability (Fig. 2 TR1). Based on future climate change, especially RCP8.5, suitable habitats (especially 75-100%), in addition to maintaining the recent range, will expand to the northwestern regions of the high latitudes. Expansion to higher elevations for habitats with 75-100% suitability is also apparent in the Zagros mountains (Fig. 2 TR2 and TR3, Fig. 3 TR1 and TR2). Additionally, the northeastern part of the country might potentially act as a potential distribution range for the species, especially under the RCP8.5. According to the RCP8.5, the eastern and southern margins of the Zagros www.nature.com/scientificreports/ mountains at low latitudes and elevation, as well as the northwest of the distribution range at low elevation, will be lost (Fig. 2 TR3 and Fig. 3 TR2). Based on the consensus model under RCP2.6 and RCP8.5, habitat loss was estimated at 8.21% and 25.74%, and new habitat gain at 52.36% and 67.24%, respectively (Table 3).
Protected area coverage. The percentage of habitat suitability and the percentage overlap of the PAs network on the habitat suitability maps for each species in recent and 2070 under the RCP8.5 scenario are provided in Table 4. All species have a small area of suitable habitat (especially ≥ 75%) within PAs, both in recent and in 2070 (Table 4). Under RCP8.5, areas of habitat with the suitability of 75-100% would decrease within PAs for 12 species, while they would increase for the remaining 18 species, though the change (decrease or increase) is not substantial (Table 4).

Discussion
The fingerprint of climate change has been reported across a variety of taxonomic groupings that are expected to undergo elevational or poleward shifts in their geographical ranges as a result of global warming in North America, Europe, and Australia 3,5,6,36,[44][45][46][47][48] . However, there is little knowledge about the latitudinal expansion of lizards under future climate change, especially in Asia. According to eSDM results, 83.33% of the 30 lizards studied expanded their range to the north at higher latitudes, while preserving their recent range. On the other hand, the range of 73.33% of species is restricted slightly (43.33%) or excessively (30%) along the southern margins at lower altitudes, but it will also persist or expand to higher latitudes. Following previous research on reptiles [49][50][51] , the findings of this study also indicate that species situated in mountainous areas (N = 19) such as the Zagros, Alborz, and Kopet Dagh may move or maintain their range at higher elevations as a result of future climate change. However, it should be highlighted that although ascending to higher elevations can provide favourable temperatures for threatened species, it can also pose challenges due to factors such as radiation, vegetation cover, and low partial pressure of oxygen (PO2) [52][53][54] , which are characteristics of high elevations and require further investigation for these species. According to the results of this study, annual precipitation for 16 (53.33%); annual mean temperature for 7 (23.33%); precipitation seasonality for 3 (10%); the max temperature of the warmest month for 2 (6.66%) and elevation for 1 (3.333%) of species were the most important variables influencing the Iranian lizards' distribution range ( Table 2). Reptiles have intermediate mobility 36 . Therefore, assumptions of unlimited or null dispersion under climate change are impossible, and future range shifts will probably fall in between 14 . Despite the fact that no research has been conducted on the dispersal ability of Iranian reptiles, particularly lizards, in the face of climate change, few studies in Europe can provide insight on this matter. For example, two southern European squamates, Hierophis viridifl avus and Vipera aspis, have shifted 60 km north in the last 40 years 55 . This is because warming in the colder northern ranges of species may open up new chances for colonization 14,56 . Moreno-Rueda et al. (2011) showed the mean latitude of the Spanish reptiles' ranges as they migrated northward at a rate of 0.5 km/year between 1940-1975 and 1991-2005 37 . They suggest that the rate of species migration to the north is influenced not only by dispersion ability but also by other variables such as geographic barriers and habitat distribution 37 . As a result, for the species under investigation, more research in these areas is required. The present study, which assumed an unlimited dispersion hypothesis, predicted the range loss and gain of 30 Iranian lizards by 2070, as shown in Table 3. In this study, retreat from lowlands or their southern areas was also observed for species (Figs. 2 and 3). Similar results were observed for the Vipera berus that retreated their distribution from the southern range in some regions of France 55 . Another example of range retraction is illustrated by field observations of many populations across the common lizard's distribution region in Europe. According to monitoring in the species' southern range, several lowland populations went extinct in 10 years, or their density was reduced by more than 50% after a warm spell 56 .
Species distribution models based on climatic factors can provide important knowledge on how species will respond to future climate change 14 . Furthermore, the findings of this study may reveal new insights into the fate of mid-latitude lizards as a result of future climate change. On the other hand, elevation can limit species ranges and has been demonstrated to have a role in explaining the distribution of species [57][58][59][60] . In this study, however, climatic factors were shown to be more significant than elevation in the distribution range of the majority of species (N = 29; Table 2), which followed previous research on reptile species richness in Iran 39,61 . Even though the results of this study shed light on how the species may respond to future climate change, it is important to acknowledge that the models in this study do not consider other factors that may contribute to lizard declines, such as anthropogenic pollution, habitat fragmentation, and loss, invasive species predation, disease, and parasitism 62 . For example, several studies have demonstrated that habitat fragmentation negatively impacts the dispersal of lizard species [63][64][65][66][67] . Restricted dispersal can lead to inbreeding, smaller population sizes, and loss of genetic variation [68][69][70][71][72][73][74][75] . However, there are few studies on Iranian lizards in this area, and more research is required, especially in light of climate change. On the other hand, non-climatic factors may have a major role in predicting the ranges of taxa 76,77 , and their inclusion in models, as well as feedback interactions between variables, is expected to improve future estimates of species extinction or decline 14 . Such factors include, for example, habitat management, the spatial distribution of habitats, human disturbance, and nutritional factors. Therefore, to address these complicated relationships, multi-factorial research would be necessary 78 .
The responses of closely related species to environmental conditions are generally similar, but species-specific responses have also been reported [79][80][81][82][83][84][85] . Depending on these two scenarios, the conservation implications may be different 79 . According to this study, climate suitability for some closely related species may be species-specific. For example, P. maculatus and P. persicus, among the three Phrynocephalus species evaluated in this study, will have the potential to expand their distribution range as a result of future climate change. However, the range of P. scutellatus may be significantly reduced (Tables 3 and 4, Fig. 3). As a result, this study suggests further www.nature.com/scientificreports/ investigation into phylogenetic niche conservatism and divergence among Phrynocephalus species, emphasizing the importance of understanding deep-time species history and speciation mechanisms before assuming common responses and conservation strategies delineation. Because it is now known that climate factors play an important role in speciation by promoting range fragmentation that leads to allopatric speciation (through niche conservatism) or promoting parapatric population divergence along climatic gradients (through niche divergence) 83 . Despite this, studies on different species of an Iranian lizard genus are rare, necessitating more research and study in this area. Conservation organizations are being encouraged to adopt proactive efforts to reduce the effects of climate change on biodiversity 26 . Despite the need to conserve Iranian biodiversity, including lizards, from climate change, stakeholders and environmental authorities have issued no specific recommendations for the management of lizards that may be threatened. This study found that a small area of highly suitable habitat exists within the PAs (Table 4). On the other hand, this study suggests that 30% of species (n = 9) may be threatened in the future, particularly along their southern margins (Figs. 2 and 3, Table 3). Additionally, the coverage of suitable habitats (75-100%) within PAs for these species (except OE) would also decrease under future climate change (Table 4). According to the findings of this study, future climate change has resulted in a loss of suitable habitats (e.g. for AC, EP, ERS, IB, OE, PC, PS, TC, and TP) as well as habitat fragmentation (e.g. for EP, ERS, IB, PC, PS, and TC) for these species (Tables 3 and 4; Figs. 2 and 3), which can lead to a reduction in population size 86 . It should also be noted that, although future climate change may result in the expansion of suitable habitats for species (Fig. 3), these new habitats may not be protected or may be less suited than existing habitats 17,87,88 . For example, (i) changing habitat may reduce food intake because new habitats are unfamiliar or of lower quality; (ii) individuals changing social environments may encounter higher aggressiveness from nonfamiliar or nonkin individuals or may prevent the evolution of helping; (iii) individuals may face increased predation risk during the dispersal phase and early in the settling phase in all cases 89 .
Monitoring programs that track lizards' temporal and spatial changes are rare in Iran, and financing such projects should be prioritized as a research priority. Consensus over monitoring schemes and collaboration, as well as monitored species, will be required to achieve these targets. In addition, experiments on the effects of climate change should also be conducted to gain a better understanding of the mechanisms, the causal pathways involved, and nonlinear reactions to future warmer temperatures 56 . This study assessed the effectiveness of the existing PAs network and identified potential conservation areas outside the existing PAs. However, more research into human activities and the presence of natural barriers in the region is required. This new data could support the development of predictive models to define management strategies and prioritize species in Iran. In conclusion, these initial findings can contribute to improving our understanding of the ecology and biology of 30 Iranian lizards, which may be applied to future research and biomonitoring programs, as well as practical conservation actions.

Methods
Study area, species, and occurrence records. This study focuses on Iran, which has a total area of 1.6 × 106 km 2 and is located in southwest Asia between the longitudes of 44° and 63° East and latitudes of 25° and 40° North (Fig. 1). The present study investigated 30 lizard species from 22 genera. These species were chosen for two reasons: (1) they had an adequate number of distribution points, and (2) their distribution range was in the west, east, north, south, center, or the entire country, allowing the response of different species across the country to be investigated under future climate change. Table 1 provides a list of these species, along with their conservation status. There are 13 species with the least concern conservation status, 12 species that are not listed, four species with data deficient, and one species that is vulnerable (Table 1). Of these, four species are endemic to Iran (Table 1). The occurrence points for these species were provided by Global Biodiversity Information Facility (GBIF, http:// www. gbif. org/). To decrease the impact of spatial autocorrelation, duplicate records were removed and occurrence records with a distance of more than 1 km were employed in the analysis 90 . The number of occurrence records used for each species is listed in Table 1. The geographical coordinates of these points are illustrated in Fig. 2. Explanatory variables. Topography and climate are introduced as the most critical factors on reptile richness at the global and regional scales 39,[91][92][93][94] . According to this, lizard's niche models were constructed for recent  and future (2070; the average for 2061-2080) climate change projections. Six bioclimatic variables with 30-s spatial resolution raster grids were downloaded from the WorldClimate (v 1.4) database (https:// www. world clim. org). These bioclimatic variables were annual mean temperature (BIO1 hereafter); temperature seasonality (BIO4 hereafter); the max temperature of the warmest month (BIO5 hereafter), annual precipitation (BIO12 hereafter); precipitation of driest month (BIO14 hereafter); and precipitation seasonality (BIO15 hereafter). BIO1 and BIO12 were chosen because they are the most influential factors for the richness and distribution range of reptiles in Iran 39 . The following four variables were selected because they are likely biologically significant, are weakly associated globally, and might indicate environmental features that limit distributions 95  www.nature.com/scientificreports/ Ensemble species distribution modelling (eSDM). The ensemble of species distribution models (eSDM hereafter) is a suitability-weighted average predicted by multiple algorithms and is one of the best or most powerful techniques for predicting habitat suitability, particularly in the face of future climate change [96][97][98][99][100][101] . Ensemble forecasting helps us to solve the issue of variability in forecasts produced by various modelling approaches or global circulation models 97,102,103 . For this purpose, the "biomod2" package (v 3.4.6) was used to simulate species distribution as an eSDM in the R (v 4.2.0) programming language 104 . The default settings recommended by Guisan et al. (2018) are used in this study 101,105 . The algorithms used in this study for all species were Flexible Discriminant Analysis (FDA); Random Forest (RF, n.trees = 1000), Generalized Boosted Models (GBM, n.trees = 1000, 3 Fold Cross-Validation); Generalized Linear Models (GLM, type = ' quadratic' , interaction. level = 1, the stepwise procedure using Akaike Information Criterion (AIC) criteria); Classification Tree Analysis (CTA, CV.tree = 50, 5 Fold Cross-Validation); Surface Range Envelops (SRE, quant = 0.025); and Maximum Entropy (MaxEnt.Phillips, maximum iterations = 500, https:// biodi versi tyinf ormat ics. amnh. org/ open_ source/ maxent/). These models (except MaxEnt and SRE) require presence and absence data and, therefore, need a set of pseudo-absence background data samples from the landscape of the study area. Since this process involves a random procedure caused by the random selection of the pseudo-absences (possibly stratified), Guisan et al. (2018) suggested establishing several pseudo-absence data sets to avoid sampling bias, especially for a moderate or low number of pseudo-absences. According to the method of Guisan et al. (2018), this study employs random sampling throughout the study area and is repeated three times with an equal number of presence data [106][107][108][109][110][111][112] . For each model, 70% of the data is used to calibrate the model (training set). The Area Under Curve-Receiver Operating Characteristics (AUC hereafter) statistics, Cohen's kappa (KAPPA hereafter), and True Skill Statistics (TSS hereafter) were used to evaluate the remaining 30% predictive capability. However, the final set is constructed with a TSS equal to or greater than 0.70 105 .
To eliminate the splitting of the total record, this process is repeated four times 105 . The TSS value ranges from -1 to + 1, + 1 means perfect agreement, and 0.60 to 0.90 means that the model performance is fair to good 113 . AUC values greater than 0.90 are considered good, those between 0.60 and 0.90 are considered average, and those below 0.60 are considered poor 114 . The importance of the variables is consistent between models that calculate the average importance of the variables used in different sets of pseudo-absences and cross-validation runs 105 .
Species range change (SRC). The species range change (SRC hereafter) was calculated using the "BIO-MOD_RangeSiz function" for each of the 30 species, as the difference between the number of sites lost (that is, the sites where the species may not exist in the future, but currently exists) and the number of sites gained by the species (that is, the number of sites that the species may exist in the future but does not currently exist) compared with the number of sites currently occupied [115][116][117] .