Introduction 

Understanding the complex interactions between climate change and ecosystem components and their effects on the landscape is the basis for developing effective species conservation strategies. These interactions can affect regional feature and ecosystem services by influencing gradual and abrupt changes within the landscape and vegetation distribution ([23]). In recent years, climate change has become one of the main drivers affecting the distribution range of plant species through its impact on biodiversity conservation ([33], [18]).

Plant species are not equally sensitive to climate change ([35]). The response of plant species to climate change depends on their genetic variability and adaptation to the magnitude of climate change as well as the availability of space to move to local microclimates. Adaptation to new conditions, movement in different directions (higher altitudes or latitudes) to find suitable climate and finally extinction on a local, regional or global scales are among the methods of coping with climate change by plant species ([34]).

In hyper-arid regions the stability of plant populations severely depend on the reproduction of adult individuals ([11]). In fact, no population is permanently viable without effective recruitment and establishment of new individuals into the population ([15]). Therefore, protecting the ability of species to respond to climate change is essential in conservation planning. This requires a comprehensive understanding of the ecological needs of plant species and a functional characteristics of the landscape and predicting the behavior of species in face of different climatic scenarios.

The term environmental heterogeneity is used to describe the number and width of environmental gradients, the diversity of microsites, the complexity of the structure, the availability of resources, and the limiting factors at a given site ([25]). Environmental heterogeneity consists of five classes that include biotic components (vegetation cover) and abiotic components (land, climate, soil and topography - [43]). In general, the effects of environmental heterogeneity components are often interrelated, and all of them affect patterns of biodiversity, population dynamics, and community structure ([3]). Many ecological phenomena are sensitive to spatial heterogeneity and flow within spatial mosaics. There is a positive relationship between landscape heterogeneity and species richness at different spatial scales ([25]). In addition, spatial heterogeneity of climatic factors can significantly affect the pattern of species distribution on a regional and global scales. Spatial heterogeneity at different scales affects ecosystem characteristics such as population structure, community composition, and several ecosystem processes (such as primary and secondary production) in both aquatic and terrestrial ecosystems ([52]). The effect of spatial heterogeneity on biodiversity ([5]), species richness ([51]), phenotypic changes ([26]) has been widely examined. However, scrutinizing published literatures reveals that the heterogeneity effect of environmental factors on the distribution of plant species in climate change scenarios has not been addressed well.

Moringa peregrina (Moringaceae) is a deciduous tree mainly distributed in Baluchistan, Bashagard and Hormozgan in southeastern and southern Iran. Most species of Moringaceae family originate in India and are distributed in tropical and subtropical regions of the world. Globally, the species grows in northeast Africa and southwest Asia. The species is in danger of extinction due to climate change, severe grazing, eradication, over-exploitation and human mismanagement ([12]). Considering the significant importance of M. peregrina in soil and water conservation, fodder production, fuel wood production, edible and industrial oil production and medicinal applications ([39]), and the small area of distribution in southern Iran ([22]), it is necessary to study the extension and fluctuation of heterogeneity of the species in response to climate change.

Species distribution models (SDMs) or ecological niche models (ENMs) can quantify the species’ response to environmental changes across the distribution range and are used to predict the appropriate main occurrence sites for the species in relation to the environmental conditions of its ecological niche. The potential effects of climate change on species distribution have been extensively studied by the ENMs ([27]). Environmental heterogeneity and area are the two important fundamentals of the species diversity. Understanding and predicting species distribution in space and time and consequently community structure and dynamics is an important issue in ecological studies, and particularly in climate change research. If the response of the species to climatic variables in the form of its distribution can be considered as the concept of ecological niche of a species, then the fluctuations that occurs in habitat suitability of the species in a variety of climate change scenarios can be considered heterogeneity. For more than five decades, ecologists have believed that the spatial heterogeneity in habitat can contribute to species richness by creating different ecological niches.

Models that predicts the distribution of species express the probability that it might occupy new ecological niches. If the difference between the habitat suitability in climate change scenarios to be quantified then the attempts that a species makes to occupy new niches or in some way to extend the niche breadth can be clarified. Niche breadth is an important concept for understanding biological adaptation and is relevant to many eco-evolutionary topics including ecological specialization, niche evolution, and speciation.

The approach of this study in evaluation the habitat heterogeneity is different from previous studies. In this study, habitat suitability fluctuations and environmental factors were compared with regard to heterogeneity. Contrary to the studies that have reported habitat heterogeneity as an important factor in maintaining diversity ([8]), in this study the ecological niche heterogeneity and environmental variables were considered as niche fluctuation and species’ efforts to occupy new areas. Therefore, the objectives of this study were to investigate the distribution of Moringa peregrina in two climate change scenarios, to identify the Region of High Heterogenetic (ROHH) of the habitats in those scenarios, and to ascertain the heterogeneity of habitat variables and compare it with current distribution of the species in southern Iran.

  Materials and methods 

Study area and presence points

The study was carried out in an area of 137.244.90 km² in south and southeast of Iran in the Saharo-Sindian bioclimatic region. The study area with an average height of 634 m a.s.l. is located in Sistan and Baluchistan, Kerman and Hormozgan Provinces in southern Iran (Fig. 1). The region is hot deserts with no frosting period and the average annual rainfall is about 150-200 mm. The most important plant species of the region are Acacia ehrenbergiana, Ziziphus spina-christi, Nannorhops ritchiana, Nerium oleander and Dalbergia sisso. Presence data of M. peregrina were recorded by handheld GPS (Global Positioning System) through field surveys and verified with species distribution maps in “Flora Iranica” ([47]) and “Flora of Iran” ([31]). In addition to the distribution areas in these two references, in areas where the presence of the species was high, its presence was also recorded randomly by field survey. The coordinates of the presence points of the species were recorded in a period of five years (2015-2020) in flowering time. To reduce the Spatially Rarefy Occurrence Data (spatial autocorrelation) between presence points, we removed duplicate presence points using only one location record per 1000-m grid cell ([10], [14]). The size of 1000 m is due to the spatial resolution of the habitat variables of the study.

Fig. 1 - Location of the study area and the presence points of Moringa peregrina (in green) in Iran.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Habitat variables

Due to the significant role of topography in creating different microclimates and as a result creating different ecological niches, this variable along with climatic variables were used to determine heterogeneity. Therefore, a set of climate and topography variables called topoclimatic variables were used as inputs for habitat suitability models ([42]). Climate variables of The Community Climate System Model (CCSM4) in the WorldClim (⇒ https:/­/­www.­worldclim.­org/­⇒ https:/­/­www.­worldclim.­org/­) with an accuracy of 30 seconds (almost 1 km) was used for the years of 2050 and 2070, mild to severe scenarios (including RCP2.6, RCP4.5, RCP6.0, and RCP8.5), respectively. To check the current distribution range of the species, the current climate data available in the same database was used. The topographic variables were Altitude, Compound Topographic Index (CTI), Wind Exposure index and Topographical Openness index. Altitude with an accuracy of 30 seconds was also downloaded from the WorldClim. The Wind Exposition index was used to obtain the possible distribution ranges of the species and the role that wind can have on distribution. Since the topographical openness of landscape refers to the relative position of valleys and heights, this variable was used instead of multiple variables which somehow expresses the extent of the landscape ([1]). The output of this analysis is two maps that divide the land surface based on concavities and convexities which have also been used to study plant habitats ([46]). These indexes were prepared in SAGA GIS (⇒ https:/­/­saga-gis.­sourceforge.­io/­en/­index.­html).

The selection of variables has a great impact on the prediction performance of the models ([37]). The SDMs (species distribution models) with many potentially relevant variables may lead to over-fitting and poor prediction performance ([24]). Environmental variables that had a correlation greater than 0.75 ([30]) were excluded from the analysis. Preparation of the variables and checking the correlation between them were done in ArcGIS® v. 10.4.1 (ESRI, Redlands, CA, USA).

Habitat modeling

Since the SDMs have different predictions, different algorithms were used to reduce uncertainties in the predictions. Four modeling algorithms including tree learning (Random Forest and TreeNET Gradient Boosting) and regression (Logistic Regression and Generalized Paths Seeker) models were used ([21], [30]). These models need presence and pseudo-absence points. Therefore, the number of 1000 pseudo-absence points was randomly created with a distance of 50 km from the presence points using the software QGIS® v. 3.16.3 (⇒ https:/­/­www.­qgis.­org/­). Such distance was considered to ensure different conditions from the presence points. 70% of the data was used for model training and 30% for testing.

All the models were run using the default settings and the importance percentage of each variable in the model was also mentioned in the SPM results. Only in the logistic regression, the interpretation of the results is different. In the logistic regression, t-statistics (t-ratio) and p-values were used to measure the effect of influencing variables on distribution. Those variables with high absolute values for the t-statistic, and thus near-zero p-values, are strong predictors. Validation was performed by separating the test and training data using criteria dependent and independent to the threshold, including AUC (area under curve), BER (balanced error rate) and class accuracy. The AUC is mostly used to evaluate habitat models. The value higher than 0.9 is very good, 0.7 to 0.9 is good, and less than 0.7 is suitable ([21]).

Balanced error rate (BER) was chosen for the study due to the imbalance of the input data (the number of presence and pseudo-absence points). BER statistic is the average of the misclassification rates on examples drawn from positive and negative classes (eqn. 1):

where TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives, and FN is the number of false negatives. To calculate the classification accuracy of classifiers, the accuracy metric is primarily used (eqn. 2):

Note that unlike the accuracy metric, for BER the lower its value, the better detection result is ([17]). The models were fitted on the current dataset and then the results were generalized to different scenarios of 2050 and 2070. The models were run with the defaults of each model in the Salford Predictive Modeler v. 8.2. Then, all the models were combined with each other based on the AUC and the Ensemble model was created for each RCP. True Skill Statistic (TSS) was used to identify the suitable habitat ([21]). TSS takes into account the sensitivity and specificity statistics and range from -1 to +1, where +1 indicates a complete match and values of zero and less indicate that the used test did not differ from the random model ([30]). These criteria were calculated using SPSS® (IBM, Armonk, NY, USA). After identifying the threshold using gap analysis, the spatial agreement of habitats in different climate change scenarios was checked with the Kappa index. The Kappa statistic is based on the overall accuracy of the model’s predictions with respect to the expected accuracy in random mode. It ranges from +1 to -1, where +1 indicates a complete match and zero and less values indicate that the model did not perform better than the random mode ([20]).

Calculation of heterogeneity

Two scenarios were considered to calculate the habitat heterogeneity. Based on them the ensemble model of distribution in RCPs 2.6, 4.5, 6.0 and 8.5 of 2050 and 2070 along with the current distribution of the species were entered into the principal component analysis (PCA), i.e., each scenario has five ensemble models. The percentage of eigenvalues of the first 3 components with a 19 × 19 (spatial scale) proximity filter was used to calculate the spatial heterogeneity. In this index, 1.0 represents the least diversity and 100 represents the highest diversity. In addition to habitat heterogeneity, the amount of heterogeneity was also calculated for each bioclimatic variable to determine the fluctuation of each habitat variable in different scenarios. For example, mean diurnal range (Bio2) layer in the current at RCP2.6, RCP4, RCP 6.5 and RCP8.5 scenarios was entered into the PCA analysis with the same proximity filter to determine its fluctuation in different scenarios of 2050 and 2070. The heterogeneity of climatic variables was also calculated with the PCA. The heterogeneity was also calculated for topographic variables. Since these variables are constant during different scenarios, then the topographic heterogeneity tool was used to calculate them. Heterogeneity analyzes were performed using the SDM toolbox ([2]) in ArcGIS^® v. 10.4.1. After creating the heterogeneities of habitat suitability, the bioclimatic and topographical factors of these heterogeneities were compared with the presence points. In this regard, the heterogeneity of habitat suitability was considered more than 50% for scenarios 1 and 2. The heterogeneities greater than 50% of the habitat suitability were done using the Raster Calculator. Then the Fishnet network was designed for this range of these areas and this range was called Region of High Heterogenetic (ROHH) in each scenario. Then, the heterogeneity of bioclimatic and topography variables per ROHH was extracted from each of these maps. The heterogeneity value of bioclimatic and topography variables was extracted for each presence point using the “extract value to point” function. Since different landforms can affect the formation of microclimates and fluctuation of habitat variables, the role of landforms in the ROHH coverage was also investigated in the scope of Fishnet.

Landform can be calculated using the TPI (Topographic Position Index) at two categories, i.e., slope position classification and landform. Since in landform classification, landform category is calculated in scales of smaller neighborhood and larger neighborhood, this method was used to obtain landform categories in this study ([40], [28]).

Landform classification was calculated using the “Topography Tools” toolbox ([50]) in 10 classes: 1 = canyons, 2 = shallow valleys, 3 = upland drainages, 4 = u-shaped valleys, 5 = plains, 6 = open slopes, 7 = upper slopes, 8 = hills in valleys, 9 = small hills in plains, and 10 = mountain tops.

  Results 

Prediction performance of models and participation of effective variables

The evaluation of model performances showed that all models had good predictive power and well separated the presence and pseudo-absence points from each other (Tab. 1). The AUC values of the models used in this study were calculated as 0.99, confirming that the prediction performance of the models was excellent. As sensitivity and specificity increase, value of BER decreases. Low BER values in the test and training data had values close to 0, which indicates that the models have achieved good values based on the average errors of sensitivity and specificity. The amount of this error for the training data in the CART and TreeNET models was lower than the Logistic Regression and GPS models. Class accuracy showed that different models were able to distinguish the presence points from the pseudo-presence points. In the training data of the CART, TreeNET, Logistic Regression and GPS models, at least 90% differentiation was observed. This value was at least 97% for the test data. Both BER and Class Accuracy metrics showed that the Random Forest model performed poorly compared to other models.

The evaluation of the importance of influencing variables showed that climatic variables are more important in species distribution (Tab. 2). Thus, Bio4 and Bio2 variables had the greatest impacts on species distribution. Results of the GPS model revealed that topographic variables did not affect the distribution of M. peregrina in the study area.

Habitat suitability maps

The current and potential distribution maps of M. peregrina in different climate change scenarios are presented in Fig. 2. In these maps, blue color shows the areas with high probability of species presence and yellow and brown colors show the areas with medium and lower suitability. The comparison of current distribution with different scenarios of 2050 shows that by increasing of temperature, the distribution expands towards the north (central parts of Iran) of the study area (yellow and green color). In the scenarios predicted for 2070, the distribution range of the species is extended to the northern parts, and the highest occupation occurs in the RCP8.5 scenario. But the southern parts are the most desirable in all scenarios (blue color). Areas with high habitat suitability in most RCPs are commonly located in the south and southeast of the study area.

Fig. 2 - Potential habitat suitability of Moringa peregrina at current and different climate change scenarios in 2050 and 2070 in southern Iran.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The comparison of the TSS maps of different climate scenarios with the current distribution of the species is presented in Fig. 3. The areas with red color are known as suitable habitats for the species in the current distribution and the compared scenario. The blue and yellow parts are desirable habitats only in the current distribution and desired scenario.

Fig. 3 - Spatial agreement of current distribution of Moringa pergerina with the potential distribution in different climate change scenarios.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Habitats with high suitability in all distribution scenarios are located in the southern parts of Iran (Fig. S1 in Supplementary material). The highest difference in the range of current distribution with the future distribution is related to RCP2.6 in 2050. In other scenarios, suitable habitats that are only in the current distribution do not have occupied large areas. But in RCP8.5 of 2050 and RCP2.6, RCP4.5 and RCP8.5 of 2070, the area of habitats that are only covered by those RCPs is high.

Heterogeneity analysis

The highest value of the heterogeneity index in the distribution range of the species is located in the southwestern parts of the suitable habitat (Fig. S2 in Supplementary material). Areas with more than 50% heterogeneity are located in the north and northwest of Bandar Abbas, south and southeast of Shiraz, and south and southwest of Bushehr, which are marked with red ovals.

In both scenarios 1 and 2, the habitat heterogeneity value is higher than the variables heterogeneity (Fig. S3 in Supplementary material). Investigating the heterogeneity of climatic variables in the ROHH in both scenarios shows that Bio2 variable has the highest and Bio12 variable has the lowest heterogeneity. Among the fixed variables, the Positive Openness variable had the highest diversity in both scenarios. Among the topography variables, the highest value of variation is related to the Wind Index and the lowest value is related to Negative Openness. The correlation between habitat heterogeneity and heterogeneity of other variables in the range of the ROHH showed that the highest correlation was between habitat heterogeneity and height heterogeneity (Fig. S4 in Supplementary material). The correlation value of this variable was 0.57 for scenario 1 and 0.55 for scenario 2. This shows that the fluctuation of height has a significant effect on the fluctuation of habitat suitability.

Comparison of the heterogeneity of habitat topoclimatic variables in the range of the ROHH and current presence points of the species shows that there is a significant difference in terms of heterogeneity of habitat variables in ROHH and presence points (Fig. 4). The diversity of topoclimatic variables in the current habitat of the species is less than the ROHH. In both scenarios, the average heterogeneity value of Bio2 in the current habitat of the species is 17%, while the heterogeneity value of this variable in the range of ROHH is 38%. The changes in the heterogeneity values of Bio3 variable also indicate that the average heterogeneity values in the presence points and the ROHH area are 10 and 20, respectively. In other words, the fluctuation range of this variable in the potential niche of the species under the influence of future climate scenarios is almost twice the current niche of the species. The heterogeneity of the topographical variables also shows that the heterogeneity of Negative Openness in the current presence points is more than that in the ROHH. Therefore, M. peregrina in the ROHH faces a closed range where the distribution range will be more towards the high altitudes. Heterogeneity by landforms indicates that the highest habitat heterogeneity caused by topography in two scenarios occurred in the canyons and mountain tops landforms. Conversely, the lowest amount of habitat heterogeneity belongs to other landforms such as upper slopes, open slopes and plains (Fig. 5).

Fig. 4 - Heterogeneity of climatic and topography variables in the ROHH range compared with the presence points of Moringa peregrina in southern Iran. S1 and S2 indicate scenario 1 and scenario 2, respectively.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Fig. 5 - Comparison of Moringa peregrina habitat heterogeneity by landforms in climate change scenarios.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

This study aimed to investigate the effect of climate change on the distribution of M. peregrina and the heterogeneity of its habitat suitability and habitat topoclimatic variables in southern Iran. Spatial heterogeneity under climate change has the potential to disrupt dispersal patterns, limiting gene flow and constraining the movement of potentially beneficial alleles across a species range ([41]).

The accuracy of the models applied in differentiating the presence and absence of the species was at an acceptable level (Tab. 2). The analysis of Class accuracy and BER criteria shows that CART, TreeNET, Logistic Regression and GPS models are able to predict species distribution with good performance. Only Random Forest did not always provide excellent performances compared to other models, despite the high efficiency reported in some studies ([19]). Our findings show that the non-linear investigation of the relationships between species habitats cannot always replace regression methods (such as logistic regression). The use of BER criterion for the SDMs (based on pseudo-absence points) can be appropriate where the data used is unbalanced.

The results of this study indicate that the potential distribution of M. peregrina covers areas outside its current range in Iran, including southern parts of Kerman province, north of Bandar Abbas, southeast of Bushehr and south of Fars province in southern Iran (Fig. 1). Farahat & Refaat ([9]) investigated the impact of climate change on M. peregrina, finding that isothermality (Bio3), mean temperature of the warmest quarter (Bio10) and the mean temperature of the coldest quarter (Bio11) had the largest impact on the distribution of the species. On the other hand, Moradi et al. ([29]) reported that rainfall had the greatest effect on the distribution of M. peregrina. These contrasting results could be attributed to the difference in the habitat variables considered, the study scale and also the presence points. The abovementioned variables express temperature changes in different ways, which shows the significant effect of this variable on species distribution. We showed that the distribution range of M. peregrina is characterized by low values of Bio4, i.e., the temperature variation over a given year (or averaged years) based on the standard deviation (variation) of monthly temperature mean. Low values of Bio4 indicate that the temperature fluctuation is lower compared to other regions. Therefore, we concluded that M. peregrina tends to establish in areas where the fluctuations in temperature over the year are limited. The same tendency was also revealed by the Bio2 variable. The low value of the monthly temperature ranges means that the distribution range of M. peregrina is relatively stable in terms of temperature, as inferred from its presence near the coastal areas of Oman Sea and Persian Gulf, where the annual temperature fluctuation is low. High fluctuation in these variables can be considered as asynchronous regimes. Recent studies have shown that asynchronous regimes can be effective on the ability of some species to respond to climate change and can lead to the creation of meta-populations ([41]). In the northern parts of suitable habitat patches (Fig. 3) and in Central Iran, the climate is hot and dry, which can limit the distribution of the species northwards. Accordingly, Farahat & Refaat ([9]) also reported that, noticeable changes in the Iranian range of M. peregrina will not occur as a result of climate change, likely due to the low rainfall and low temperature fluctuations.

On the other hand, seed germination studies have shown that flower abortion is high in very dry regions ([16]). This might explain why the species has not been observed in the northern regions of Kerman and Sistan and Baluchistan Provinces. Indeed, the delineation of some parts of the northern strip of Iran (Golestan, Mazandaran and Guilan Provinces) as a suitable habitat for M. peregrina seems mainly due to the climatic stability of the coastal regions of northern Iran (Fig. 2). It has been reported that current habitats of M. peregrina are mainly located on the edge of the sea and ocean, such as the Oman Sea, Gulf of Aden, Red Sea, Gulf of Aqaba, Gulf of Souns ([7], [9], [12]). The comparison of the current distribution with that predicted under different climate change scenarios showed a high spatial overlap (kappa coefficient > 90%). Proximity to coastal areas involves the reduction of extreme temperature fluctuations and the increase of air humidity, which increase the suitability of the habitat for M. peregrina. One of the effects of proximity to oceans in tropical and low latitudes regions is the sea breeze, which can blow up to several kilometers to inlands and transfer humidity and cool air. The penetration range of the sea breeze into inland depends on the temperature difference between inshore and offshore. A greater difference can cause deeper penetration of the sea breeze into inlands and vice versa. It seems that a set of these factors has a special effect in increasing the habitat suitability for M. peregrina in the study area compared to other areas. The consequences of these conditions might be the increase of the suitability of these regions for the species distribution.

Species can face climate changes either by adaptation to changing environmental conditions or migration towards more suitable habitats. M. peregrina has the potential to spread towards suitable habitats at higher altitudes in the current presence areas, and higher latitudes could be favorable for the species under both climate change scenarios. This is in line with the findings of Moradi et al. ([29]) in the study of the distribution range of M. peregrina based on the SDMs.

M. peregrina is known as a species adapted to hot and dry regions with insufficient soil moisture ([45]). One of the important characteristics of arid and semi-arid ecosystems is the discontinuous distribution of vegetation, forming a mosaic of patches separated by areas without vegetation cover. This is the consequence of heterogeneous patterns of soil moisture, which is a driving factor for plant establishment and growth in such regions. Currently, M. peregrina mostly occurs in valleys and waterways, and most occurrences of the species used in this study were from the landform of waterways. Accordingly, our results showed an expected distribution of the species mainly in canyons and mountain tops landforms in the ROHH range under both scenarios (Fig. 5), in particular in areas where suitable environmental conditions for the species occur, i.e., where soil moisture is sufficient ([36]).

In both the climate scenarios considered in this study, the heterogeneity in habitat suitability was far higher than the heterogeneity of environmental conditions. This suggests that the predicted habitat suitability can be the result of a certain degree of fluctuation of climatic variables in the region. Among the analyzed climatic variables, the Bio12 (expressing the annual precipitation) had the lowest fluctuation and its variation is very small in the ROHH under both the climate scenarios considered. Therefore, we concluded that precipitation does not impact on the predicted distribution of the species. On the other hand, the Bio2 variable (the average daily temperature) showed the largest variation.

In the ROHH, the habitat heterogeneity is highly correlated with the elevation heterogeneity (Fig. S4 in Supplementary material). It is well known that the altitude is one of the driving factors in the distribution of the species and in the changes in vegetation along gradients ([48]). The areas characterized by extreme altitudinal heterogeneity can be characterized by a larger differentiation of ecological niches and thus support higher species richness, thereby protecting ecological and evolutionary processes and making population dynamics more resilient to climate change. Indeed, based on the concept of niche differentiation, more heterogeneous environments can support more species through partitioning the niche space ([49]). On the other hand, shifts towards higher elevation due to climate change (as predicted in this study for M. peregrina) could lead to the fragmentation of populations due to the shrinkage of suitable habitats with elevation, making this species more susceptible to the deleterious effects of climate change ([16], [4]).

The maximum altitude in the ROHH area is 2899 m a.s.l., while the optimal altitude range for M. peregrina was reported to be 1000-2000 m by Gebauer et al. ([13]) and 300-500 m by El-Keblawy & Khedr ([7]). It seems that the elevation of the species in southern Iran is higher than in other regions. The distribution range of the species is negatively related to the increase in altitude (Fig. 4). This can seriously limit the distribution of the species by negatively affecting some plant traits such as density, seed production, plant reproduction, and phenology. It has been reported that plants at higher altitudes are more at risk compared to lower altitudes ([7]). The negative effect of altitude and the positive effect of climatic variables on habitat suitability for M. peregrina have been previously reported in several studies ([16], [9], [29]). Moreover, habitat heterogeneity was also correlated with Bio4 and Bio5, which can reflect the role of these variables in plant access to moisture and nutrients, especially in arid and semi-arid regions ([44]).

Comparing current habitat conditions with those predicted from future climate scenarios, it can be inferred that M. peregrina currently grows in a more open landscape than that predicted by our models, which can affect the presence of pollinators and influence the wind speed. On the other hand, the combination of open landscape and wind can facilitate seed dispersal and thus the distribution of plants in suitable habitats. In line with this finding, it has been reported that there is a positive relationship between the openness of the landscape and altitude with the diversity of pollinators ([38]).

Climate change is not a linear process which can solely be described by changes in temperature means ([6]). Our results showed how the physical parameters of the habitat interact with climate and how this interaction could affect the distribution range of the species. A crucial factor in determining the composition and dynamics of populations is habitat heterogeneity, i.e., the number of structural elements in a given locality, which has been proven to support a higher species diversity ([32]). Finally, our findings support the use of the SDMs and topoclimatic variables as viable tools for investigating the fluctuations of the ecological niche of a species.

  Conclusions 

This study aimed to investigate the current distribution of M. peregrina, its response to climate change scenarios, and its niche fluctuation through the comparison of current distribution with the distributions in different climate change scenarios. Our results indicate that the heterogeneity at the habitat level and in the topoclimatic variables differs between the current and the potential range of the species in southern Iran under future climate change scenarios.

Our models did not anticipate drastic changes in the habitat of M. peregrina under different scenarios of climate change, suggesting that the species is well adapted across its current range in southern Iran, with particular regard to temperature fluctations and elevation. However, the latter represents a possible limitation to the shift of the species range to western and northwestern Iran, while the former is a possible limitation to possible shifts to Central Iran.

Identifying the potential habitats of plant species is one of the vital strategies both to cope with current threats due to climate change and to properly manage their populations. Moreover, SDMs are efficient tools to predict possible changes in species distribution range and habitat heterogeneity, in order to identify the limits of its ecological niche and the environmental limitations to its distribution under future climate scenarios.

  Acknowledgements 

Authors have equally contributed in conception and design of the study, data collection, modelling and mapping, drafting manuscript and approval of article. Authors would like to express their gratitude to the Vice-chancellery for Research and Technology, University of Zabol, Iran for funding the project (grant no. IR-UOZ-GR-8694).

  References

  Supplementary Material

Fig. S1 - Overlap between different classes of habitat distribution in climate change scenarios with the current potential distribution of Moringa peregrina.
Fig. S2 - Habitat heterogeneity and distribution range of Moringa peregrina in climate change scenarios in southern Iran.
Fig. S3 - Comparison of habitat heterogeneity values and habitat topoclimatic variables in the ROHH range of Moringa peregrina in climate change scenarios in southern Iran.
Fig. S4 - Correlation between habitat heterogeneity and habitat variables affecting the potential distribution of Moringa peregrina in southern Iran.