Estimating potential range shift of some wild bees in response to climate change scenarios in northwestern regions of Iran

Background: Climate change is occurring rapidly around the world, and is predicted to have a large impact on biodiversity. Various studies have shown that climate change can alter the geographical distribution of wild bees. As climate change affects the species distribution and causes range shift, the degree of range shift and the quality of the habitats are becoming more important for securing the species diversity. In addition, those pollinator insects are contributing not only to shaping the natural ecosystem but also to increased crop production. The distributional and habitat quality changes of wild bees are of utmost importance in the climate change era. This study aims to investigate the impact of climate change on distributional and habitat quality changes of five wild bees in northwestern regions of Iran under two representative concentration pathway scenarios (RCP 4.5 and RCP 8.5). We used species distribution models to predict the potential range shift of these species in the year 2070. Result: The effects of climate change on different species are different, and the increase in temperature mainly expands the distribution ranges of wild bees, except for one species that is estimated to have a reduced potential range. Therefore, the increase in temperature would force wild bees to shift to higher latitudes. There was also significant uncertainty in the use of different models and the number of environmental layers employed in the modeling of habitat suitability. Conclusion: The increase in temperature caused the expansion of species distribution and wider areas would be available to the studied species in the future. However, not all of this possible range may include high-quality habitats, and wild bees may limit their niche to suitable habitats. On the other hand, the movement of species to higher latitudes will cause a mismatch between farms and suitable areas for wild bees, and as a result, farmers will face a shortage of pollination from wild bees. We suggest that farmers in these areas be aware of the effects of climate change on agricultural production and consider the use of managed bees in the future.


Introduction
In recent years, the declining pollinator population has been a major concern globally (Viana et al. 2012), leading to more research on investigating pollinators' threatening factors and analyzing their reduction effects on agricultural and natural systems. Approximately, 88% of angiosperms (Ollerton et al. 2011) and 87 of the 115 most important food products require pollinators (Klein et al. 2007). Global temperatures are predicted to rise by an average of 2 to 4°C by 2050 due to greenhouse gas emissions (Pachauri et al. 2014). These climate changes will shift the geographic distribution of pollinators and jeopardize food security in the future (Imbach et al. 2017). There are three possible scenarios related to species' responses to a regional climate change: (1) adaptation to new conditions, (2) migration to more suitable areas, and (3) extinction (Coope 1995). Although plantspollinators relationships have developed throughout the evolution process, the climate changes are happening at a faster pace, and the plants and pollinators cannot adapt to these changes (Yurk and Powell 2009). Since bees cannot fly during cold, windy, wet, and dark times, wild bees are more adaptable to future climatic conditions (Christmann and Aw-Hassan 2012). Rader et al. (2013) showed that under the worst scenario of climate change, pollination by honeybees would be reduced by 14.5%, while native bee pollination will increase by 4.5% by 2099.
Climate change can affect foraging activity, body size, and longevity of pollinating insects (Scaven and Rafferty 2013). For example, increasing the ambient temperature reduces the body size of bees (Schweiger et al. 2010). It is claimed that the effects of climate change on tropical insects are more significant because these species have narrow niches (i.e., they are specialist species that require very unique resources), and live in an environment that is close to their optimum temperature (Kjøhl et al. 2011). However, species that are at higher latitudes live at temperatures below their optimum temperature, consequently, increasing the temperature enhances the foraging activities of these species (Deutsch et al. 2008). Simulations show that 17 to 50% of pollinating insects are more likely to experience periods that do not access the floral resource for foraging due to the climate change effect (Memmott et al. 2007). The most important effect of climate change on the plants-pollinators relationship is temporal and spatial mismatches (Rafferty 2017;Schweiger et al. 2010). Early flowering plants seem to be less affected by climate change, and even in some cases benefit from climate change (Takkis et al. 2018). However, climate change effects can change the phenology of the early flowering plants (Karlík and Poschlod 2014). Rising temperatures negatively affect agricultural products at low latitudes while having a positive effect on products in high latitudes (Challinor et al. 2009). If bee foraging activity begins earlier than flowering time, survival rates and population size may decrease (Kjøhl et al. 2011). Moreover, the effects of climate change affect the spatial distribution of plants and pollinators, nectar quantity and quality, pollinator energy demand, plant community structure, and pollinator community structure (Schweiger et al. 2010).
To reduce the adverse effects of climate change on pollinators, it is essential to evaluate the effects of climate change on the distribution of pollinators (Polce et al. 2014). In this way, the methods that can correlate the current climate and pollinators' distribution can determine the areas that will be suitable for species in the future (Dew et al. 2019). Species distribution models (SDMs) are numerical tools that combine species presence points with environmental factors. In this type of model, the key assumption is that the species are in equilibrium with their environment (Elith and Leathwick 2009). As input, SDMs require georeferenced individual locations or species' presence as the response or dependent variable, and independent layers of environmental information such as climate, topographic slope, elevation, land cover, and soil attributes. The central idea of SDMs is the niche theory that was introduced by Joseph Grinell and G. Evelyn Hutchinson (Miller 2010). They distinguished two niches, namely, the fundamental niche and the realized niche. The fundamental niche encompasses all abiotic environmental conditions that a species is physiologically able to tolerate, while the realized niche comprises those parts of the fundamental niche where a species can survive despite the presence of competitors (Zurell et al. 2020). Therefore, the realized niche is geographically smaller than the fundamental niche due to negative interspecific interactions. Species niches are not only climate-related niches as many other environmental factors affect the presence of species and their habitat suitability. There are five main steps for SDM studies: conceptualization, data preparation, model fitting, model assessment, and prediction (Zurell et al. 2020).
Species distribution studies have shown that with increasing temperature, species distribution moves toward higher latitudes and polar regions (Hegland et al. 2009). Bumblebees, for instance, are cold-adapted species that have shifted to higher latitudes and mountain ecosystems over the past 100 years (Kerr et al. 2015;Ploquin et al. 2013;Pyke et al. 2016). In North America, the distribution models have shown that the diversity of Bombus species will decrease in southern regions in the future (Sirois-Delisle and Kerr 2018). Several studies have shown that the suitable habitats for wild bees will decrease due to future climate change (Giannini et al. 2012). For example, Suzuki-Ohno et al. (2020) claimed that the distribution range of five species of bumblebees has decreased after 26 years because of the climate change effect. Imbach et al. (2017) indicated that by 2055, suitable areas for coffee will be reduced by 88%. However, other studies have shown that rising temperature increases the suitable areas for bees. For example, in the case of carpenter bee species in Australia, Dew et al. (2019) found that climate change increased the suitability of this species. Besides, climate change reduces the overlap between suitable areas for horticultural and wild bee (Polce et al. 2014). For example, Bezerra et al. (2019) argued that climate change altered the natural range of Xylocopa bees as well as areas suitable for passion fruit production, resulting in a reduction in overlap areas between the two species.
About 14.46% of Iran is covered by agricultural land, of which 54% is irrigated and 46% is rainfed farmlands. Approximately, 30% of Iran's agricultural products need pollinating insects. There are about 800 species of bees in Iran that belong to the families Colletidae, Halictidae, Andrenidae, Melittidae, Megachilidae, Anthophoridae, and Apidae (Mohammadian 2003). A wide range of pollinating bees has been recognized in different parts of Iran. Among them, two species of Apis florea and Apis mellifera meda have a wider distribution (Sanjerehei 2014). A. florea is distributed in southern parts of Iran and Apis mellifera meda has the most widespread distribution among other bees in the country (Rahimi and Mirmoayedi 2013). The economic value of pollination service in agricultural productions was estimated at $6.59 billion in Iran in 2005-2006(Sanjerehei 2014. $5.72 billion was allocated to the honey bees and $0.87 billion to the wild bees. It has been estimated that pollinators affect 25% of total agricultural products in Iran. To the best of our knowledge, no study has been conducted to model the effects of climate change on wild bees in Iran. Therefore, the present study aims to predict the future distributions of five wild bee species under climate change scenarios in northwestern regions of Iran.

Study area
From an ecological point of view, Iran is divided into three phytogeographical regions: the Euxino-Hyrcanian province of the Euro-Siberian region in the north, the Saharo-Sindian region in the south, and the Irano-Turanian region in the western and central sectors of the country (Talebi et al. 2014). The study area in this study is located in two regions of Euxino-Hyrcanian and Irano-Turanian that also include the two mountain ranges of Alborz and Zagros. The area of the study area is 357090 km 2 and mostly has a mountainous climate that is cold in winter and warm in summer ( Fig. 1).

Species occurrence data and environmental predictors
In this study, we used presence-only records of five wild bee species belonging to three families available at the Global Biodiversity Information Facility (GBIF) website (www.gbif.org). These species are Halictus  Fig. 1. Melitta leporina is one of the solitary bees that nest in the ground and the maximum foraging distances of this species in alfalfa fields are estimated to be 25 to 30 cm. Melitta leporina is the most common species among the 17 members of the genus Melitta and is found throughout Europe and around the Caspian Sea, Iran, and Kazakhstan (Celary 2006). Megachile maritima is on the IUCN red list in Germany, the Netherlands, and Ireland. This species is polylectic and collects pollen mostly from the Asteraceae and Rosaceae families (Güler and Bursali 2008). Anthidium florentinum is a solitary bee the size of a honeybee. It is widespread in northern Africa, Europe, and Asia and has been introduced to the USA. Females nest in pre-existing cavities and cover the cavities with leaves (Wirtz et al. 1992). Megachile rubripes Morawitz is a leaf-cutter bee that is found in the northern and southern regions of Iran and is associated with the Asteraceae and Lamiaceae families (Samin et al. 2015). Halictus tetrazonianellus Strand is spread in Iran, Turkey, Azerbaijan, Lebanon, and northern Russia that is associated with the Asteraceae families (Dikmen and Aytekin 2011). We chose these species to compare the effects of climate change on those that nest in the ground (Halictus tetrazonianellus Strand and Melitta leporina) and those that are above-ground nesting and nest in cavities (Megachile maritima, Megachile rubripes Morawitz, and Anthidium florentinum).
The bioclimatic variables as environmental predictors were obtained from the WorldClim database (www. worldclim.org). The WorldClim data layers included 11 temperature and eight precipitation metrics with a spatial resolution of about 1 km 2 . Layers with a correlation greater than 0.8 were excluded from the modeling process, and finally, layers including isothermality (Bio3), temperature seasonality (Bio4), maximum temperature of the warmest month (Bio5), minimum temperature of the coldest month (Bio6), mean temperature of the wettest quarter (Bio8), mean temperature of warmest quarter (Bio10), annual precipitation (Bio12), and precipitation seasonality (Bio15) were selected. We modeled the future possible distribution of all five species under two representative concentration pathway scenarios (RCP 4.5 and RCP 8.5) in the year 2070. Scenarios RCP 4.5 and RCP 8.5 are known as moderate and worst scenarios where the average Earth temperature rises from 1.4 to 1.8°C for the moderate scenario and 2 to 3.7°C for the worst scenario (Stocker et al. 2014).

Species distribution modeling
According to response variables, species distribution models are divided into two categories: discrimination and profiles. Discrimination models require the presence and absence records of the species of interest and perform the modeling process based on the correlation between species records and environmental variables. Discrimination models are divided into two categories: global (parametric) and local (non-parametric) (Tarkesh and Jetschke 2012). Global models include generalized linear models (GLM), multiple logistic regression (MLR), Gaussian logistic regression (GLR), logistic regression (LR), and artificial neural network (ANN). Local models include generalized additive model (GAM), classification and regression tree model (Milano et al., 2019), nonparametric multiplicative regression (NPMR), logistic regression tree (LRT), random forest (RF), and boosted regression tree (BRT). In contrast, profile models use species presence-only data to model the possible distribution of species (Miller 2010). These models including BIOCLIM, DOMAIN, and maximum entropy (MAX-ENT), and ecological niche factor analysis.
In this study, we used the SDM software package (Naimi et al. 2016) in R for modeling potential native bee distribution. To achieve this, we used three classification models including RF, CART, BRT, and two regression models including GLM and GAM to create habitat suitability maps based on climate data. For predicting possible future distributions under RCP 4.5 and RCP 8.5 in the year 2017, we used the RF model because it showed higher performance than other models evaluated in the current study (see the "Results" section).
Random forest (RF) models are among the machine learning methods that provide a single prediction by ensembling several poorly predicted models. This algorithm creates a forest randomly. The constructed forest is a set of "Decision Trees." The construction of the forest using trees is often done by the method of "bagging." The main idea of the bagging method is that a combination of learning models enhances the overall results of the model. Simply, a random forest makes several decision trees and combines them to make more accurate and consistent predictions. This model creates a large number of decision trees and then combines all the trees for the final prediction of the most parsimonious species distribution. The RF model split the predictor variable based on the homogeneity of the responses. Each split point divides the response into two nodes so that there is maximum similarity within the nodes, and minimal similarity between the nodes (Valavi et al. 2020). The number of trees and the number of predictor variables are the key parameters of the random forest model (Jafarian et al. 2019). The high efficiency of the random forest method has been confirmed by many studies (Bradter et al. 2013;Cushman and Wasserman 2018;Rather et al. 2020), especially when the presence points are limited (Strobl et al. 2007).
The CART algorithm creates a decision tree with binary divisions. It tests the input variables to find the best split so that the impurity index (Gini) resulting from the splitting has the lowest value. In the analysis, two subsets are split based on homogeneous values of the dependent variable and each is split into two other subgroups in the next step using the same approach, and the process continues repeatedly (Jafarian et al. 2019). As the tree grows, the nodes become more homogeneous, and more information becomes visible. Eventually, the splitting process stops due to a lack of data.
The boosted regression tree (BRT) is an ensemble of statistical techniques and machine learning. This model is one of the techniques that improve the performance of a single model by using a combination of multiple models. This model combines the strengths of two algorithms: CART and boosting. Boosting is a way to increase the accuracy of a model, and the basis of its work is that building, combining, and averaging a large number of models is better than creating a model (Elith et al. 2008). The basic idea in this model is to combine a set of weak predictive models to achieve a strong prediction. The generalized linear model (GLM) is an example of a linear regression model that estimates the statistical relationship between a response variable and one or more predictor variables (Guisan et al. 2002). This model is one of the parametric models and is used for cases where the observations are not normally distributed and other regression models are not suitable. This model has been developed to reduce limiting assumptions in linear regression and can be fitted for a range of distributions (Binomial, Poisson, Gaussian, Bernoulli, Gamma, etc.).
The generalized additive (GAM) model is a nonparametric model and an extension of generalized linear models. The generalized additive model provides a suitable method for analyzing data and examining the relationship between dependent and independent variables (Kosicki 2020). Generalized additive models (GAM) are superior to generalized linear models (GLM) in several respects and the purpose of using these models is to maximize the predictive quality of the dependent variable, to discover non-linear relationships between the dependent variable, and the set of explanatory variables. Unlike the generalized linear model, in which the relationship between explanatory variables and response is represented by a formula, the data are allowed to determine the shape of the response curve.
Most ecological niche models use presence and absence data for modeling. It is very difficult to obtain species absence data because it requires a lot of sampling effort in the study area. Background points are samples that are randomly extracted from the study area to be used along with presence data for modeling species presence-only data. The central idea of the pseudo-absences points is to allow species distribution models to compare the desired regions for the species to other regions throughout the study area (Valavi et al. 2020). Furthermore, background points are used to quantify environmental conditions in the region of interest, so the number of required points differs depending on the complexity of environments in the region (Renner et al. 2015). Barbet-Massin et al. (2012) recommended the use of a large number of pseudo-absences in using regression models like generalized linear models (GLM) and generalized additive models. For classification techniques such as random forest (RF), boosted regression trees (BRT), and classification trees, they recommend using the same number of available presence data for pseudoabsences. Therefore, we selected 10000 random points from the background as the pseudo-absence points for GLM and GAM models and generated the same number of presence data from the background as pseudoabsences data for classification models like RF, CART, and BRT.

Model assessment
We evaluated the model performance using 10 runs of subsampling replications by randomly splitting our occurrence data into 70% as training data and 30% as test data. Two statistics, true skill statistic and the area under the ROC curve (AUC), were used to verify the model performance. Models with AUC values of 0.7-0.9 are acceptable, whereas the values higher than 0.9 are regarded as models with excellent abilities. TSS ranges from −1 to +1, where +1 indicates perfect classification, and values between 0.4 to 0.6 show a moderate performance of the model. The suitability maps derived from SDMs range from 0 to 1, with the number one indicating the maximum suitability for the species. In this study, suitability maps were divided into three categories: high (0.7-1), moderate (0.3-0.7), and low (0-0.3), so that changes in species distribution range under climate change scenarios would be quantified.

Uncertainty and sensitivity analysis
Uncertainty can be present in different steps of the implementation of species distribution models. Accurate coordinate recording and spatial autocorrelation between absence/presence data in the first step, type and the number of environmental layers, and the correlation between them, spatial extent and grain size of raster layers in the second step, and the selection of different algorithms for modeling habitat suitability maps in the last step are cases that cause uncertainty in SDMs (Fernández et al. 2013). In this study, five models including RF, CART, BRT, GLM, and GAM were used to show the uncertainties associated with species distribution models. Different algorithms produced different results, so we used the mentioned models to provide habitat suitability maps only in the current scenario and compared the results of the models. From these models, we determined the most efficient algorithm according to AUC and TSS statistics and predicted the effects of climate change on the possible future distribution of wild bees under two RCPs in the year 2070.
To perform sensitivity analysis, first, we estimated the relative contribution of the environmental layers for the five mentioned algorithms for species under study. Then, we removed the environmental layer that had the highest contribution in the modeling process and again produced a new habitat suitability map to determine the effect of removing one layer or the most important layer on the results. Sensitivity analysis was performed only for the most efficient algorithm (i.e., random forest). Table 1 shows the results of the model performance using TSS and AUC statistics. In this table, the items with the highest mean of AUC and TSS are bolded. AUC and TSS for RF model are AUC = 0.8 ± 0.04; TSS = 0.53 ± 0.1; AUC = 0.76 ± 0.08; TSS = 0.46 ± 0.1; AUC = 0.85 ± 0.04; TSS = 0.68 ± 0.05; AUC = 0.82 ± 0.05; TSS = 0.61 ± 0.06; AUC = 0.86 ± 0.04; TSS = 0.66 ± 0.07 for Melitta leporina (ML), Megachile maritima (MM), Halictus tetrazonianellus Strand (HT), Anthidium florentinum (AF), and Megachile rubripes Morawitz respectively. AUC statistic shows that the random forest method has the highest performance for all species than other models. However, the TSS statistic shows moderate efficiency for the random forest method in determining the habitat suitability map for Melitta leporina and Megachile maritima species. For this reason, we used only the random forest (RF) model to predict the possible future distribution of bees under climate change scenarios. Table 2 shows estimates of relative contributions of the predictor environmental variables to RF, CART, BRT, GLM, and GAM models. According to this table, RF, CART, BRT, and GAM models show Bio3 as the most effective layer in the habitat suitability modeling process for Melitta Leporina (ML). However, GLM indicated Bio6 as the most effective layer. For Megachile Maritima (MM), Bio6 had the highest contribution in RF, BRT, and GLM models. Bio4 and Bio18 were considered as layers with the highest impact in habitat modeling in GAM and CART models respectively. Halictus Tetrazonianellus Strand (HT), GAM model indicated Bio18 as the most effective layer and in others, Bio3 had the highest importance in the habitat suitability modeling process. For Anthidium Florentinum (AF), RF and BRT models showed Bio12, GLM, GAM, and CART indicated Bio5, Bio15, and Bio18 as the most effective layer respectively. For Megachile Rubripes Morawitz, Bio3 had the highest contribution in RF and CART models. GLM, GAM, and BRT models showed Bio12, Bio5, and Bio8 as layers with the highest importance in the modeling process. Table 3 shows the area and percentage of different classes of habitat suitability maps for different models and species. The results of this table show that different algorithms have different results for each species. For example, for Melitta leporina (ML), BRT has modeled only two classes with low and moderate suitability, but other models have considered an average of 5% of the study area in the high suitability class. RF and GAM models have shown similar results. Megachile maritima (MM), different models have also shown different results, so that the CART model has considered about 10% of the study area as a high-suitability class, while the BRT and GLM models have shown only two classes of low and moderate suitability. For other species, this difference can be seen in the area and proportions of habitat suitability classes obtained from different models. Table 3 shows that there are significant uncertainties in the use of different models that can lead to misleading results.

Uncertainty and sensitivity analysis
Sensitivity analysis of the results of the random forest model was performed based on the removal of the environmental layer, which had the highest relative participation in habitat modeling. For Melitta leporina (ML), Halictus tetrazonianellus Strand (HT), and Megachile rubripes Morawitz, Bio3 was the most effective layer in the habitat suitability modeling. For Megachile Maritima (MM) and Anthidium florentinum (AF), Bo6 and Bio12 were the layers with the highest impact on the modeling process. Table 4 shows that the removal of the most effective environmental layer from the modeling process has a significant impact on the results of the random forest model. For Anthidium florentinum (AF), for example, with the removal of Bio12, about 6% of the study  area is considered in the high-suitability class, while Table 3 showed that in the case where no layer was removed, the RF model only showed two classes with moderate and low suitability. For other species, the removal of the most important layer has less effect on the results of the RF model.

Possible future distributions based on the random forest model
Since the model assessment results showed that the random forest model had the highest performance among all models, we only used this model to predict the possible future distribution of bees under climate change scenarios. Figure 2 shows the habitat suitability maps for the present (row A) and the year 2070 under scenarios RCP 4.5 (row B) and RCP 8.5 (row C). In this figure, the habitat suitability is between 0 and 1, and numbers close to one indicate a high degree of suitability. According to this figure, it can be seen that the species under study can mainly occupy a small area of the study area potentially at present. However, as time goes on and the temperature of the area rises, it can be seen that a large area of the study area will be favorable for wild bees and the distribution of the species will shift to higher latitudes. The highest rate of the shift in distribution can be observed for Halictus tetrazonianellus Strand, which under a scenario of RCP 8.5 covers a large area of the study area. Figure 3 shows the classification maps of habitat suitability in high, moderate, and low classes. This figure shows the current maps (row A), scenarios 4.5 (row B), and 8.5 (row C). According to this figure, it can be seen that with increasing temperature in scenarios RCP 4.5 and RCP 8.5, the range of the suitable areas for the species under study will increase, and their distribution will be shifted to higher latitudes. This increase is more significant for Halictus tetrazonianellus Strand, Anthidium florentinum, and Megachile maritima. For Melitta leporina, the range of suitable areas will decrease with increasing temperature. It is noteworthy that the high suitability class for all species in scenarios RCP 4.5 and RCP 8.5 has disappeared, and only two classes (moderate and low) have remained. Table 5 shows the area and percentage of habitat suitability classes under scenarios 4.5 and 8.5. According to this table, the range of suitable areas for four species will increase with increasing temperature in 2070 under RCP 4.5 and 8.5. For Melitta leporina, with increasing temperature, the range of suitable areas will decrease in 2070. It is noteworthy that for Melitta leporina, under both scenarios, the high-suitability class will disappear and the area of the moderate class will reduce by 20%. For Megachile maritima under both scenarios, the area of parts with moderate suitability will increase by 20% and the class with high suitability will disappear. For Halictus tetrazonianellus Strand, the class with high suitability will disappear in the year 2070 and the area of the class with moderate suitability will increase by 38% under RCP 4.5 and 42% under RCP 8.5. For Anthidium florentinum, no high-suitability class was identified for the current scenario, but in 2070, the area of the moderate class will increase by 12% under RCP 4.5 and 6% under RCP 8.5. For Megachile rubripes Morawitz, the high-suitability class will disappear with increasing temperature in the year 2070, and the area of the moderate class will increase by 8% and 6% under RCP 4.5 and RCP 8.5 respectively. In general, our results showed that with the increase in temperature in 2070 under two scenarios, RCP 4.5 and RCP 8.5, the area of suitable classes for wild bees will increase, but this increase will occur only in the moderate class, and the high-suitability class will disappear. For one species, increasing the temperature will reduce the desired areas for it.

Discussion
Our results showed that with increasing temperature under scenarios RCP 4.5 and RCP 8.5, the possible distribution of wild bees would shifted to higher latitudes. This result has been reported by other studies. For example, bumblebees have been reported to shift to higher latitudes over the past 100 years (Kerr et al. 2015; Pyke  et al. 2016) because these species are more adaptable to cold weather and prefer lower-temperature habitats. In mountainous ecosystems, as the temperature rises, the bumblebees move to higher altitudes (Ploquin et al. 2013). For M. subnitida, Giannini et al. (2017) under scenario RCP 8.5 also found that the position of the best areas for this species will change and the possible distribution of this species will move to higher latitudes. Therefore, it seems that species that live in cold regions move to higher latitudes with a cooler climate as the temperature increases. This shift of wild bees to higher latitudes can jeopardize food security in our study area and deprive agricultural products of pollination by these bees. Decreased overlap between suitable areas of crops and wild bees due to climate change has been reported in various studies (Bezerra et al. 2019;Polce et al. 2014).
Our results also showed that with increasing temperature, the desired areas for the species under study will increase, and these species would have a possible wider range in the future. The result is different for Melitta leporina and the possible distribution ranges of this species will decrease with increasing temperature. Dew et al. (2019) reported that the range of suitable areas for Ceratina australensis under climate change scenarios in 2070 would increase due to climate changes. Rader et al. (2013) also showed that under the worst scenario of climate change, native bee pollination would increase by 4.5% by 2099. Climate change has been reported as the major factor in reducing the population of bumblebee species in Britain (Williams et al. 2007) and declines in bees food in Europe (Dormann et al. 2008). However, Giannini et al. (2012) showed that the optimal habitat of bees decreased due to climate changes. This result has also been reported for several species of bumblebees (Suzuki-Ohno et al. 2020). In addition to species distribution, rising temperatures can also affect species richness. For example, Carrasco et al. (2020) showed that in North America, under scenarios RCP 4.5 and RCP 8.5, the richness of tomato pollinators decreased due to climate changes. Giannini et al. (2020) also showed that out of 216 species under study, 95% of them would experience a decrease in their populations due to climate change. It seems that the effects of climate change on the species distribution ranges are not constant and depending on the type of species and the study area, these effects can be different. Therefore,  (Rader et al. 2013), therefore, if species respond differently to climate change, so do agricultural products.
Although the range of species distribution will increase with increasing temperature, the quality of suitable areas will decrease in terms of climate because for all species under study we found that the high-suitability class disappeared under scenarios RCP 4.5 and RCP 8.5. Therefore, this study emphasizes more on the negative effects of climate change because for one species, these effects would reduce the suitable areas, and for the other four species, the increase in temperature would cause the loss of habitat quality in terms of climate. Therefore, studies that report an increase in wild bee habitat because of rising temperatures should also consider habitat quality changes. It is almost impossible for bees to adapt to climate change because climate change occurs so quickly and bees will have to migrate to higher latitudes or climatically suitable habitats or they will become extinct. In addition, climate change, along with the fragmentation of bee habitats, insecticides, and diseases, will increasingly reduce the population of these species. Different responses of insects and plants to climate change can cause temporal and spatial mismatch for affected species. Mismatches between plants and pollinators reduce pollination on farms by wild bees. If climate change alters the distribution ranges of wild bees in the future, farmers will be forced to use managed bees to pollinate their crops, which could incur additional costs.
The results of the uncertainty analysis showed that different models obtain different habitat suitability maps and as a result, there is significant uncertainty in the use of these models. For different species, the models used in this study obtained different estimates of habitat suitability classes. Therefore, this point should be considered in the species distribution modeling studies. To eliminate uncertainty between these models, the results of all models are usually ensembled to overcome the existing limitations (Jafarian et al. 2019). Using the ensemble method can lead to better predictions when the models have acceptable performance. In this study, because the random forest model showed higher performance than the other models, we preferred to investigate the possible effects of climate change on species distribution based solely on this model. The results of sensitivity analysis also showed that the number of environmental layers used to model the distribution of species has a significant effect on the results. This study showed that the removal of the most effective environmental layer from