The study area included Khuzestan Province (KP), Bushehr Province (BP), Fars Province (FP), and Kohgiluyeh and Boyer-Ahmad Province (KBP) with a total area of about 223000 km ² in the southwest of Iran (Fig. 1). In this area, there were two relatively distinct topographic situations: (i) vast arid plains in the south of KP and west of BP with high temperatures in summer averaging ∼37°C and mild winters averaging ∼17°C ( IRIMO, 2010 ), and (ii) mountainous areas in the north of KP, east of BP as well as the other two provinces (FP and KBP) with moderate temperatures in summer averaging about 25°C and cold winters averaging about 6°C ( IRIMO, 2010 ).

Field survey was done in the study area during March-May of 2015-2017 and presence points of C. pabotii were obtained. This species was classified by the Flora Iranica ( Wagenitz, 1980 ) as well as examination of relevant taxonomic literature ( Negaresh & Rahiminejad, 2015b , 2016 ). Presence points were recorded by a GPS with an error of less than 10 m. The spring season was chosen to conduct field surveys because growing and flowering of C. pabotii in this season facilitated the identification of the species. Altogether, we obtained 25 presence points (Fig. 1).

Habitat quality map was obtained to show areas with high quality for C. pabotii using MaxEnt v3.4.1 ( Phillips et al. , 2017 ). MaxEnt only needs a few presence points to create pseudo-absence points (here we selected 10,000 points), and compares environmental variables in presence cells with values of the same variables in pseudo-absence points to find the discrimination between them ( Phillips et al. , 2006 ; Almasieh et al. , 2016 ). Finally, the software creates a continuous quality map. The area under the receiver operating characteristic (ROC) curve (AUC) in MaxEnt (0.5-1, in which '0.5' indicates random and '1' indicates complete discrimination of presence points from pseudo-absence points) was used to evaluate model performance for both habitat quality modelling and climate change predictions. Jackknife test within MaxEnt was used to evaluate the relative importance of each variable and response curves within MaxEnt were used to show the probability of species presence for each variable considering the correlation between each variable and other variables. We randomly selected 75% of the presence points as the training data set and the other 25% as test data ( Pearson et al. , 2007 ).

Environmental variables including topographic, climatic and land-cover were used for habitat modelling (Table 1). Digital Elevation Model (DEM) at about 1-km spatial resolution was used to create the slope variable using Spatial Analyst tool. Solar radiation variable (kJ m ^-2 day ^-1 ), which was classified for each month of the year (averaged period of 1970-2000) ( Fick & Hijmans, 2017 ), was used for modelling. We averaged solar radiations of April-May using Raster Calculator tool, because this quarter of the year was necessary for growth of C. pabotii . Regarding C. pabotii 's dependence on wheat fields, we extracted agricultural lands, the majority of which were wheat fields, from the land-cover variable ( FRWMO, 2010 ) and distance from agricultural lands was created using Euclidean Distance tool. All the tools mentioned are available in ArcGIS 10.2.

Among 19 climate variables ( Hijmans et al. , 2005 ), we chose annual mean temperature (BIO1), tempe­rature annual range (BIO7), annual precipitation (BIO12), precipitation of the wettest month (BIO13) and precipitation of the wettest quarter (BIO16). We included those climate variables which we assumed had the greatest physiological impact on C. pabotii ( Bradley, 2014 ); Annual temperature and precipitation (BIO1 and BIO12), maximum temperature of the warmest month minus minimum temperature of the coldest month (BIO7), as well as precipitation of the wettest periods (BIO13 and BIO16) which were important for growth of C. pabotii were chosen. We removed BIO1 due to its high correlation with DEM. All variables had a correlation of less than 0.85 ( Yang et al. , 2013 ).

We categorized continuous habitat quality map based on maximum training sensitivity and specificity logistic threshold in MaxEnt ( Jimenez-Valverde & Lobo, 2007 ). Therefore, 0.248 was considered as a threshold of suitable and non-suitable categories. We also categorized the suitable category into three sub-categories: 0.246-0.4 as low quality, 0.4-0.6 as moderate quality, and 0.6-1 as high quality ( Hu et al. , 2015 ).

We used MaxEnt to predict the effect of climate change on the distribution of C. pabotii in the future. MaxEnt can predict future distribution of species more reliable than other prediction models ( Hijmans & Graham, 2006 ). We performed modelling projection in MaxEnt using current bioclimatic variables ( i.e. , BIO1, BIO7, BIO12, BIO13 and BIO16) as well as future bioclimatic variables. We downloaded the projection model of CCSM4 for the year 2070 (average 2060-2080), which was from the fifth phase of coupled model inter-comparison project (CMIP5). Four representative concentration pathways (RCPs) were available based on prediction of greenhouse gases concentration ( i.e. , RCPs 2.6, 4.5, 6.0 and 8.5). Each pathway explained the radiative forcing value as the amount of energy entering the atmosphere and the quantity of reflex. By the end of the year 2100, the radiative forcing levels are estimated to reach 2.6, 4.5, 6.0 and 8.5 W/m ² of surface ( Van Vuuren et al. , 2011 ; Yousefi et al. , 2015 ). We chose RCPs 2.6 and 8.5 as the lowest and highest radiative forcing levels for this study.

We assessed changes in C. pabotii distribution due to climate change for the future under both scenarios of RCPs 2.5 and 8.5. For that, we categorized habitat quality map of the current and future distribution into two categories ( i.e. , non-suitable and suitable) based on maximum training sensitivity and specificity logistic threshold in MaxEnt. Then, we calculated the suitable area for the present and the future, as well as the changes in the suitable area in the future under both scenarios. Also, to assess the trend of shifts in the suitable area, we calculated the centroids of merged suitable polygons for present and future maps using Hawths tools ( Beyer, 2004 ) and finally, we compared the spatial change of centroids of suitable polygons ( Hu et al. , 2015 ).

The AUC of ROC curves were 0.984 and 0.978 for training and test data respectively, which indicated the high performance of the model. The high quality category was mostly located in KP and also several small ranges of the high quality category were located in KBP and BP. FP was the province with the smallest area of high quality category (Fig. 2).

BIO13 variable was the most important variable that solely provided the most useful information (Fig. 3A). C. pabotii preferred precipitation of more than 100 mm in the wettest month (Fig. 3B). Also, DEM, BIO16 and BIO12 variables were the second, third and fourth most important variables, respectively (Fig. 3A). C. pabotii inhabited areas with an elevation of less than 1000 m and preferred precipitation of more than 200 mm in the wettest quarter and 300 mm as annual precipitation (Fig. 3B). Distance from agricultural lands was the variable that, if omitted, the highest reduction in gain of the model would appear. Therefore, distance from agricultural lands provided the greatest volumes of information, which were not present in other variables (Fig. 3A). Due to the dependence of C. pabotii to agricultural lands, the probability of the presence of this species decreased as distance from agricultural lands increased (Fig. 3B).

The AUC of ROC curves were 0.963 and 0.884 for training and test data of scenario of RCP 2.6, and 0.934 and 0.928 for training and test data of scenario of RCP 8.5 respectively, which showed the high performance of model. Our results revealed that the suitable area of the C. pabotii distribution will increase in the future in all provinces. Nevertheless, KP showed the highest increase in suitable area in the future under both scenarios (Table 2 & Fig. 4). We observed a shift for C. pabotii toward lower latitudes in the future based on centroids of polygons, as well as a shift toward higher elevations (averaging 770 m to 1406 m in both scenarios). The suitable area in the future will approximately double under the two scenarios and will expand to the surrounding areas, especially toward the center of KP and north of BP. Also, a new relatively large and isolated suitable polygon was observed in FP under both scenarios (Fig. 4). We observed more expansion of the suitable area in RCP 8.5 (Fig. 4B) with respect to RCP 2.6 (Fig. 4A), particularly more expansion toward the center of KP.