The study domain encompasses the equatorial belt and the Southeast Asian, comprising China, India, Thailand, Laos, Vietnam, Philippines, Cambodia, Bangladesh, Nepal, Singapore, Bhutan, Indonesia, Papua New Guinea, Timor-Leste, Australia, Vanuatu, Fiji, New Caledonia, and Solomon Islands ( Figure 1 ). The region is topographically complex, dissected by numerous rivers, and dominated by tropical rainforest climate, tropical monsoon climate and alpine climates (Wheeler and White, 2002; Vadrevu et al., 2019). The five focal bee species involved are widespread throughout the region, and largely ensure the stability of agro-pastoral production and ecosystems in the region. The occurrence data of the five honey bee species (i.e. Apis dorsata, Apis florea, Apis laborosa, Apis andreniformis and Apis cerana) were obtained from the Global Biodiversity Information Facility (GBIF) database (GBIF, https://www.gbif.org/), China Academic Journal Network Publishing Database (http://www.cnki.net/), El Colegio de la Frontera Sur Database (http://www.ecosur.mx/) and Google Scholar (https://www.google.com/) (Kumar and Stohlgren, 2009; Kumar, 2012; Santana et al., 2019). We checked each occurrence data for the five honeybee species, and eliminated occurrence records with duplicate and unknown coordinates (Tang et al., 2020; Yuan et al., 2020). For species distribution points with only specific place name, Google Earth Pro (https://www.google.com/earth/versions/#earth-pro) was used to query precise coordinates corresponding to the place name, and the outliers were discarded (Li et al., 2020; Tang et al., 2021a; Zhao et al., 2021). In addition, the geographic distribution records with excessively ambiguous labels (e.g., ‘Thailand’) was likewise excluded (Tang et al., 2021b). The actual ranges of some honeybee species may exceed the areas delineated, and the study assumes that species distributed outside the study area (e.g., at higher latitudes) are ecologically and genetically distinct from the populations in the study area (Martínez-López et al., 2021; Wang et al., 2023). After these steps, 20,583 validated occurrences remained for the five bee species.

To mitigate oversampling bias, sampling deviation correction was applied to calibrate the distribution of each species, but different correction methods have certain limitations (Kramer-SChadt et al., 2013). The limited number of distribution points were retained within a certain distance according to the spatial screening method, and systematic sampling was performed again according to the principle of uniformity and randomness to ensure that the geographic autocorrelation of the datapoints was minimized (Tang et al., 2021a, 2021; Wen et al., 2024). Meanwhile, the sampling results were compared with the sampling results of cluster grouping and group modeling (Tang et al., 2021b). While ensuring sufficient modeling data to the greatest extent, the distribution points behave as a true reflection of the bee species distribution during the modeling process (Boria et al., 2014). Using these methods, the 5893 unique occurrences for five bee species were remained to participate in modeling. Finally, the distribution points of five bee species were exported as csv format and mapped with ArcGIS 10.4.1 (Esri, Redlands, California) to visualize spatial coverage ( Figure 1 ). This map identified areas that were better surveyed, guiding possible caveats in our models and highlighting the interpretation of results.

Environmental variables directly affect habitat suitability and are central to forecasting climate responses of species (Hirzel and Lay, 2008; Gogol-Prokurat, 2011). Twenty environmental data were selected to predict habitat suitability of five bee species, including 19 bioclimatic variables and altitude variable (see Table 1 ). The historical bioclimatic variables were downloaded from the WorldClim database (http://www.worldclim.org) with a resolution of 30′′, which were raster data generated by the kriging interpolation method based on the observation data of meteorological stations worldwide from 1970 to 2000 (Fick and Hijmans, 2017; Poggio et al., 2018). The altitude raster with 30 arc seconds was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn/). Future projections were generated for 2050 and 2070 using created nine combinations of three greenhouse gas concentration scenario (RCP 4.5, 6.0 and RCP 8.5) and three Global Climate Models (GCMs: CCSM4, HadGEM2-AO, and MIROC-ESM-CHEM) (Gao et al., 2013; Martínez-López et al., 2021). RCPs is the abbreviation of representative concentration pathways (Drouet et al., 2015; Gebre, 2015). According to the Fifth Assessment Report (AR5) released by Intergovernmental Panel on Climate Change (IPCC), RCP4.5 and RCP 6.0 represent two intermediate greenhouse emission scenarios and RCP8.5 represents the highest emission scenario (Drouet et al., 2015; Martínez-López et al., 2021). The different climate change projections represented by the RCPs take full account of future atmospheric greenhouse gas concentrations, policy factors, and land use change (Meinshausen et al., 2011). For each RCP scenario, different GCMs present different aspects of future climate change because they consider different parameters (Ahn et al., 2021). It is accepted that these models are similar in variables used to develop them, but the models also differ in processes and inflows (Zhang et al., 2016). For example, HADGEM2-AO was established and developed mainly considering troposphere and aerosols, while MIROC-ESM-CHEM was developed based on nutrients from phytoplankton and zooplankton in the ocean (Collins et al., 2011; Gent et al., 2011; Watanabe et al., 2011; Baek et al., 2013). Based on this, the specificity of each model and the differences between models were captured and extended to predict the species distribution under different climatic conditions. We established a uniform standard framework to select the environmental variables used to model the potential distribution of each species. First, twenty environmental variables were selected for each honeybee species to build the model. Secondly, the multicollinearity of variables will lead to over-fitting of the model and reduce the accuracy of prediction results, so Pearson correlation coefficient method was used to calculate the correlation of the paired variables (Padalia et al., 2014; Tang et al., 2019). If the correlation coefficient of two variables is greater than 0.8, it means that the two variables can be substituted for each other (Wang et al., 2010; Khanum et al., 2013). On the basis of fully considering the ecological significance of variables, the less ecologically informative member should be removed (Li et al., 2020). Third, to determine which variables to retain in the final model, we prioritized the variable contributions obtained by the model under default parameters and the studies related to the survival and niche of honey bee species (Tang et al., 2021b).

From the ENM projections under current climate scenario, we quantified the contribution rate of each environmental variable involved in the modeling (Tang et al., 2021a). Based on ENM projections, we further analyzed the prediction results for 2050 and 2070 under different RCP conditions. We first quantified the suitable habitat for the five honey bee species. Following reclassification, only pixels scored as highly or moderately suitable were retained for subsequent area calculations. We quantified the distribution difference between present and projected future distributions to estimate climate-driven habitat losses for each honey bee species in 2050 and 2070 (Martínez-López et al., 2021). In addition, it is assumed that the presence of at least one honeybee species within a given area will basically guarantee the local agricultural and forestry production (Tscharntke et al., 2005; Ssymank et al., 2008). Consequently, we considered the five bee species as a whole in order to determine the distribution of pollinators in the study area. We overlapped the suitable distribution of the five species under same RCP and model at same time periods, and consider high, moderate and low suitable distribution are given decreasing priority. Any predicted site classified as highly suitable for at least one of the five bee species was designated as highly suitable area overall. Any predicted site simultaneously encompassed by highly and moderately suitable area were assigned the highest suitability class. Conversely, concordant unsuitability across all five species at a given site indicates a complete loss of honey bee services under projected climate change.

Climate warming typically drives upslope shifts in species’ altitudinal distribution (Gottfried et al., 1999; Chen et al., 2011). To account for the elevation changes under different RCP conditions, we quantified mean elevation of suitable habitat for the five bee species under each RCP for 2050 and 2070. Under different RCP at current and future, We randomly sampled 10,000 points within the suitable distribution areas of five honey bee species by GCM (Martínez-López et al., 2021). It should be noted again that the randomly sampled suitable distribution here also include only highly and moderately suitable areas. We then extracted the altitude information of the sampling points by ArcGIS 10.4.1 software (Zhou et al., 2010). Based on elevation data obtained of each species, we calculated the mean elevation of the five species in their current and future distribution areas.

We measured the direction and distance of centroid shifts for each bee species under each RCP-time combination by ArcGIS 10.4.1 (Tang et al., 2021b). Meanwhile, we estimated the geographical extent and area of the nature reserves within the study domain. To assess reserve effectiveness, we calculated the proportional overlap between nature reserves and species distribution under each climate scenario (Martínez-López et al., 2021). These results provide us with the minimum and maximum area losses of suitable distribution for bee species. To evaluate the overlap rate between nature reserves and species distribution under climate change, we considered five bee species as a whole, and therefore the suitable areas were overlapped to determine the final species distribution areas (only include highly and moderately suitable areas). The shapefiles of protected natural areas can be obtained from Protected Planet (https://www.protectedplanet.net), which has the most comprehensive information of protected natural areas for each country in the study area.

To assess the relationship between species distribution range and their ability to cope with climate change, we used Pearson correlation coefficients to analyze the relationship between the current distribution and the mean projected change in the future distribution under each climate scenario (Martínez-López et al., 2021). Additionally, we separately evaluated the environmental suitability of each bee species’ distribution records under each RCP condition (Teixeira et al., 2014). Prior to correlation analysis, we normalized current and future suitable distribution of each species using the following formula.

CSD means current suitable distribution and FSD means future suitable distribution.

where x is the current suitable distribution (CSD) or future suitable distribution (FSD) of a species. x ¯ is the mean of all the species’ CSD or FSD. σ is the standard deviation of the CSD or FSD of all species. The standardized procedure effectively avoids the influence of geographical range size on the analysis results when comparing current distribution reductions of each bee species. The resulting dimensionless value, that is the standardized CSD/FSD, was used to test the relationship between mean change of all standardized FSDs and standardized CSD (Martínez-López et al., 2021). For FSDs, we separated the correlation analysis into two time frames (2050 and 2070), and each frame encompassed nine GCM-RCP combinations. To evaluate climate-driven changes in environmental suitability, we extracted the suitability values and FSDs for every occurrence point of each bee species under each GCM-RCP combination. We then calculated the differences between CSD and FSD under climate effects, and determined the suitability changes for each species under different combinations of models and scenarios. Finally, we visualized the suitability changes for each species through box plots.

The current and future suitable distributions of the five bee species are predicted by the optimized ENM based on GCM-RCP combinations. The training and test AUC values exceeded 0.90 for all combinations of climate scenarios, and the highest AUC can reach 0.98. These values significantly exceed the random distribution and have little difference, indicating that the prediction results have high quality performance. Partial AUC analysis also supports this conclusion. With the maximum sensitivity plus specificity, 95.2% of species distribution points can be accurately predicted, and 91.4% of species unsampled or non-distributed areas can be accurately predicted. With the 10 percentile training presence, 93.4% of species distribution points can be accurately predicted, and 89.6% of species unsampled or non-distributed areas can be accurately predicted.

Across all bee species, Bio1, Bio5, Bio7, Bio12, Bio18 and elevation consistently explained the predicted current distribution ( Table 2 ). Bio1, Bio12 and elevation emerged as the primary determinants of suitability across all bee species, with the cumulative contribution of the three variables reaching a maximum to 54.9% for A. andreniformis and a minimum to 43.8% for A. cerana.

Occurrence records for the five honey bee species involved in the modeling exhibit marked spatial heterogeneity across the study region ( Figure 1 ). Occurrence records are markedly clustered on the Korean Peninsula, China’s eastern seaboard, China-Indochina Peninsula and Malay Peninsula, indicating that the above areas were well surveyed. Conversely, the number and density of occurrence records decreased across much of the Malay Archipelago.

The distribution ranges of the five bee species are first considered as a whole. By 2050, the highly suitable distributions of the pollinator are projected to concentrate in the Yangtze River basin in China, all of Japan except Hokkaido, the western side of the Western Ghats, the southern edge of the Qinghai-Tibet Plateau, the Irrawaddy River basin, the northeastern plateau of Thailand, the Truong Son Ra region, all of Cambodia and the Philippines, and the islands of Kalimantan, Sulawesi, and Sumatra ( Figure 2 ). The moderately suitable distributions are spread outward with the highly suitable distributions as the center. By 2070, the highly suitable distribution will expand slightly along the southern edge of the Tibetan Plateau, and the suitable habitats in other places will contract markedly ( Figure 3 ). Notably, almost all highly suitable areas in the China-Indochina Peninsula and the Malay Archipelago will nearly vanish. Climate warming is expected to significantly contract climatically suitable habitat for pollinators across the study region, and the highly suitable habitat will show an obvious fragmented distribution.

ENMs estimate that the current suitable distribution (CSD) of the five studied species range from 0.7342×10⁶ km² (A. laborosa) to 4.3853 ×10⁶ km² (A. cerana) ( Figure 4 ). Based on current distribution area, expert experience and the Labougle (1990), all five species studied are common species (corresponds to rare species). By 2050, all GCM-RCP combinations indicates continued range contraction of the future suitable distribution for pollinators. The loss of honey bee species distribution ranged from 7% (A. floreafor, RCP 6.0, GCM HadGEM2-AO) to 72% (A. laborosa, RCP 8.5, GCM MIROC-ESM-CHEM) ( Figure 4 ). By 2070, the pollinator distribution ranges are projected to contract markedly across all GCM-RCP combinations. The loss of distribution ranged from 14% (A. florea, RCP 4.5, GCM CCSM4) to 87% (A. laborosa, RCP 8.5, GCM MIROC-ESM-CHEM) ( Figure 4 ).

Projections indicate that the median elevation of suitable habitat for all five bee species will continue to shift upslope by 2050 and 2070 ( Figure 5 ). By 2050, mean elevation variation ranged from 35m (A. florea, RCP 4.5) to 360m (A. andreniformis, RCP 8.5) at. By 2070, the effect of climate change on the pollinators elevation change will become more pronounced. By 2070, mean elevation variation ranged from 75m (A. dorsata, RCP 4.5) to 450m (A. laborosa, RCP 8.5). As pollinator habitats shift to higher elevations, most species are projected to experience progressively greater habitat loss under RCP8.5, peaking by 2070. Meanwhile, species whose climatically suitable ranges lie at higher elevations are more significantly affected by climate change.

Centroid analyses reveal that the centroid of all five bee species changed significantly under climate change, primarily in distance and direction ( Figure 6 ). The suitable distribution of A. dorsata, A. laborosa, A. andreniformis and A. cerana will shift southeastward by 2050 and to continue in the same direction by 2070. Conversely, A. florea exhibits a north-eastward displacement. Across all species, the geographic centroids of suitable ranges are projected to shift 65–137 km by 2050 and 2070. Overall, RCP8.5 scenario exerts the strongest influence on geographic-centroid shifts across all species by 2050 and 2070.

We overlapped the suitable distributions of all bee species (only include highly and moderately suitable areas) into one whole. Currently, 17% of the total suitable habitat falls within protected natural areas ( Figure 7 ). By 2050, as habitat contracts, 20% (range 18%-22%) of the remaining suitable area is protected by the nature reserve. By 2070, this proportion rises modestly only 28% (range 24%-32%).

A species’ resilience to climate change scales with its geographic range size. Correlation analyses revealed strong positive relationships between CSD and mean FSD change in 2050 (r = 0.76, df = 16, p < 0.01) and in 2070 (r = 0.86, df = 17, p < 0.01) ( Figure 8 ). The correlation between CSD and mean FSD change will be stronger by 2070, and pollinators will be more strongly affected by global warming. Our models predict that honey bee species with smallest current suitable ranges will be more negatively impacted by projected climate change in the future. This negative impact will be more pronounced by 2070 as warming accelerates. Overall, climate change will diminish pollinators availability across their suitable habitats.

Environmental suitability is projected to decline for all species across every GCM-RCP combination by 2050 and 2070 ( Figure 9 ). Species with smaller ranges face proportionally larger suitability losses over the coming decades. Comparing the median losses of all GCMs, median losses under RCP 8.5 exceed those under RCP 4.5 and 6.0 for both 2050 and 2070.

Ecological niche modelings (ENMs) can accurately predict and analyze the suitable habitat range of species, and reveal the relationship between climate variables and species presence (De Marco et al., 2008; Evans et al., 2015). As the most commonly used species-climate analysis method, ENMs also capture large-scale species distribution information, which provides richer environmental context for establishing a systematic and comprehensive relationship between species and the environment (Elith and Leathwick, 2009; Zhu et al., 2013). Therefore, ENMs are currently the most widely accepted and well-fitted models for species distribution prediction. The International Union for Conservation of Nature (IUCN) has begun to apply the model to species invasion, endangered species conservation and climate-change responses (Cassini, 2011).

To predict the impacts of global climate change on low-latitude pollinator distributions, we constructed ecological niche models under current and future climatic conditions. The modelled ranges and climate envelopes of all five bee species closely match their currently documented distributions within the study region (Hepburn et al., 2005; Radloff et al., 2010; Nagir et al., 2016; Raffiudin et al., 2020). It is worth noting that climate is not the only decisive factor affecting the process of species distribution dynamics (Root et al., 2003; Davis et al., 2010). Intrinsic adaptive capacity, habitat disturbance, reproductive traits, land-use change, pesticide exposure, pathogen pressure and anthropogenic activities will affect the survival and distribution of the species (Otto et al., 2016; Potts et al., 2016b; Heneberg and Bogusch, 2020). These factors, as important external pressures, limit the suitable distribution of pollinators and may potentially aggravate species declines under climate warming. Meanwhile, the construction of the niche models also must be carefully calibrated to approximate each species’ fundamental niche (Owens et al., 2013). Nonetheless, climate is the most important factors that determine the survival and reproduction of species, profoundly shaping the geographic distribution and ecological functions of pollinators (Hegland et al., 2009; Uden et al., 2015). We quantified the climatic niche of five bee species using nineteen environmental variables, and evaluated model accuracy with AUC values from both training and test datasets. We believe that our models provide sufficiently accurate estimates of species distributions. Accuracy analyses also confirmed the robust performance of our SDMs.

Environmental factors are closely related to the geographic distribution of insects (Lobo, 2016). The study area, centered on Southeast Asia at low latitudes, is strongly influenced by maritime climate regimes. The unique temperature-precipitation patterns and geography conditions provide essential environmental information for the survival and distribution of pollinators in the region (Abrahamczyk et al., 2011; Bartomeus et al., 2011). Among the many climatic factors that affect the suitable distribution, only a few exert dominant control over insect ranges (Stiling, 1988). Under current climates, Bio1 (annual mean temperature) and Bio12 (annual precipitation) primarily constrain suitable distributions for all five pollinators, with the lowest cumulative contribution of both exceeding 43%. Reproduction is the basic life activity of insects to maintain population, and temperature is an important influencing factors (Régnière et al., 2012). Annual average temperature can characterize the thermal requirements for each pollinator, and temperature deviation from the optimum impair insect reproduction and physiology (Janowitz and Fischer, 2011). Elevated temperatures impair neural function and affect the short-term memory of bees, which in turn will affect their performance as adults (Jones et al., 2005). In addition, temperature anomalies further destabilize the balance of the population (Savage et al., 2004). Investigation of the correlation between weather conditions and winter mortality of bee species in Austria revealed that warmer and drier weather conditions in the preceding year were associated with elevated winter mortality of bee species (Switanek et al., 2017). To reduce the negative effects of increased temperatures, most species shift upslope to higher elevations where temperatures are more suitable. Altitudinal surveys of the species abundance for bee populations on the Kilimanjaro region also reveal temperature as the primary driver of bee abundance at higher elevations (Classen et al., 2015). Elevation-associated changes in temperature, precipitation, vapor pressure and wind speed reduce both abundance and beta-diversity (Perillo et al., 2017). However, there are interspecific differences in climate change adaptation. Honey bee species inhabiting harsh environments generally exhibit greater tolerance to climate variability and extremes than those in normal climatic conditions (Abou-Shaara et al., 2012). Global warming is an important challenge to pollinators, especially those at low latitudes. The pollination process requires honey bees to fly a certain distance, and the sustained muscle activity raises their body temperature well above ambient (Heinrich and Buchmann, 1986). Bees will stop foraging and pollination in hot environments when their body temperature approaches their critical thermal limit (Willmer and Stone, 2004). Global warming will inevitably reduce the pollination efficiency of most honey bee species, so artificial selection of heat-tolerant bees may be an important way to address the climate challenge (Abou-Shaara, 2016). In conclusion, temperature constrains pollinator diversity by restricting the resource acquisition of ectothermic pollinators under global change.

In addition to temperature, precipitation and ambient humidity also profoundly influence insect physiology and population ecology (Barton and Ives, 2014). Our model indicated that annual precipitation is another key factor limiting suitable distribution of each bee species. Environmental moisture primarily influences the water balance of insects, thereby regulating the growth and development of pollinators (Benoit, 2010; Jactel et al., 2012). The nationwide study of pollinators across Ethiopia found that the precipitation of the driest quarter and the precipitation of the warmest quarter are the best contributors to bee distribution, and fluctuations in precipitation affect the pollination services as well as the ecological and economic roles of honey bee species in the country (Abrha, 2018). As with temperature, altered precipitation patterns in the current and preceding years can alter bee abundance, and plant-resource abundance is likewise affected (Moeller et al., 2012). Both excessively high and extremely low soil-moisture levels greatly reduce the hatching success of honey bee eggs (Dupraw, 1961). Future climatic change will cause irregular shifts of the rainfall belt in equatorial regions, further exacerbating the uneven distribution of water resources (Mamalakis et al., 2021). The survival of honey bee species affected by rainfall will eventually shape the structure of the community. In addition to altering species communities, future changes in rainfall patterns will also influence key physiological and morphological traits of individual insects. Body size increased over time, more so at drier sites (Suni and Dela Cruz, 2021). At the same time, the increase in body weight reduces water loss but, in turn, limits their flight range (Peterson et al., 2006). Pollinators must make developmental decisions without complete information about future conditions, and this uncertainty affects the structure of populations (Forrest et al., 2019). Precipitation may be an important driver of trait shifts in tropical bees, shifts that arise from differential climatic tolerances and thereby motivate investigations of natural selection on color and body size.

Climate change will affect the relationship among human activities, the abiotic environment and organisms, and may further aggravate the loss of pollinator species (Settele et al., 2016). Based on the optimized ENM, we predicted the suitable areas for low-latitude pollinators and analyzed distribution patterns in relation to the key limiting factors. Drawing on multiple databases, we conducted a detailed and comprehensive statistical analysis for each honey bee species in the study region. Therefore, it is reasonable to conclude that the sampling of species distribution points is adequate. Nevertheless, we are still uncertain whether the five bee species studied are transferable across space and time (Martínez-López et al., 2021). This is mainly because the existing distribution point data can provide additional information on the distribution of the species across temporal and spatial scales. Pollinators are widely distributed in the study area, but regions such as Myanmar, Sumatra, Kalimantan and New Guinea islands remain properly and adequately surveyed. We should also be cautious about prediction results in areas characterized by low point density and limited sampling.

Based on an combined analysis of different RCP scenarios and models, we found that the suitable distribution of all honey bee species in the study area will decrease under climate change. The hypothesis that climate change will reduce the suitable range of bee species has been supported by previous studies (Soroye et al., 2020). In a study of regional crops and pollinators, climate change significantly reduced the suitable range of ten Brazilian bee species, particularly Melipona bicolor and Melipona scutellaris (Giannini et al., 2012). It is worth noting that the negative impact of climate change on honey bee populations is inevitable even in climatically suitable sites (Faleiro et al., 2018). Meanwhile, the contraction of suitable pollinator habitat also exhibits clear spatiotemporal patterns. Correlation analysis showed that the loss of pollinator habitat in 2070 was more significantly correlated with the size of the current distribution area than in 2050. Compared to species with large suitable ranges, those with smaller current ranges experience a greater climate-induced reduction in range size. Similar conclusions have also been confirmed in studies of bumble bee species in Mesoamerica (Martínez-López et al., 2021). In addition, environmental suitability for all species is projected to decline, with analyses for 2050 and 2070 indicating that RCP 8.5 could cause the greatest reduction. In the Pacific Northwest, environmental suitability for bumble bees is expected to decline markedly, with that of Bombus. vandykei projected to fall by 63% by 2050 (Koch et al., 2019). The pollinator responses to climate change examined in this study align with projected global range losses of other honey bee species (Le Conte and Navajas, 2008; Willmer, 2014; Pyke et al., 2016).

Our ENMs predicted that the loss of the suitable habitat for low-latitude pollinators will be concentrated in the equatorial region, especially on islands belonging to Indonesia and in the southern regions of Southeast Asia. The spatial mismatch between pollinators and plants caused by habitat loss will disrupt the symbiotic relationship within the ecosystem (Gorostiague et al., 2018). By 2070, only the higher latitudes within the study area are expected to remain highly suitable areas for pollinators. During recent climate warming, many insect species have shifted their ranges to higher latitudes and altitudes (Hill et al., 2011). In addition, Chapman et al. (2012) suggested that poleward migration of insects can promote population increases, highlighting their strong capacity to respond to environmental change (Chapman et al., 2012). Not only are low-latitude insect pollinators affected by climate change, but marine biodiversity near the equator has also declined significantly, with some populations now at risk of extinction (Beaugrand et al., 2015; Davies et al., 2017). Meanwhile, marine organisms tended to migrate poleward, whereas honeybee species showed no comparable shift (Molinos et al., 2016). Centroid migration trends among the five honey bee species could be influenced by variable winds and ocean currents in the study area, but further validation is needed (Jha, 2015). Moreover, the loss of suitable area was mainly concentrated at the edge of suitable range, while the core habitat persisted. Numerous studies have supported this hypothesis, but controversies remain over the underlying causes (Brown, 1984; Martínez-Meyer et al., 2013; Lee-Yaw et al., 2018). Supplementing our models with field investigations at the edge of suitable habitats may clarify this controversial phenomenon of range loss in these species.

Changes in altitude and environmental suitability revealed that future climate change will shift the suitable distribution of insect pollinators to higher elevations. Environmental suitability for all species will decrease in the future. In particular, the most significant decrease in environmental suitability was observed under the RCP8.5 scenario. Many studies have found that most honey bee species shift their altitudinal ranges under climate warming, thereby moving their habitat centroids (Hoiss et al., 2012; Inouye, 2020; Nooten and Rehan, 2020). Insects often migrate to higher elevations to track more suitable habitats (Adedoja et al., 2018). However, this upward shift also results in a loss of suitable habitat for pollinators at higher altitudes. Against the background of marked warming at lower latitudes, pollinators with lower tolerance for warmer or cooler climates are unable to shift toward higher latitudes (Kerr et al., 2015). In addition, the east-west mountains of the continent often block northward movement for pollinators. Kerr et al. (2015) found that the mean elevation of bumble bees’ suitable distribution in Europe rose by 300 m because of the Alps (Kerr et al., 2015). It is worth noting that many studies have shown that vegetation will continue shifting poleward and upslope under climate change (Ashton et al., 2009; Hagedorn et al., 2019). Vegetation shifts at low latitudes may be a key driver of range loss in insects. Although higher altitudes still provide some plant resources for pollinators, plant species richness there is significantly lower than at lower altitudes (Lawton et al., 1987; Becker et al., 2005). Our results suggest that the highly suitable habitat for pollinators is mainly concentrated in the Himalayan region and is projected to shift upslope by up to 500 m in the future. This warrants further investigation, as the upslope shift observed across all bee species may result from their populations being influenced by the east–west-oriented Himalayan mountain range.

Protected terrestrial areas are widespread worldwide and play a key role in protecting endangered species and fragile ecosystems under climate change (Deguise and Kerr, 2006). However, protected areas have often been selected not to protect specific species, but to meet other criteria (Dunn, 2003). For example, Jiuzhaigou Valley in Sichuan and Yushan National Park in Taiwan were both established to protect natural landscapes (Wang et al., 2018). The Wudalianchi Nature Reserve in Heilongjiang Province was established to protect unique geological landforms (Huang et al., 2018). The protected terrestrial areas within our study region are relatively small, spatially fragmented, and poorly connected. Meanwhile, some of these protected areas are themselves at risk from climate change (Hoffmann et al., 2019; Elsen et al., 2020). According to our results, the proportion of suitable pollinator habitat that falls within protected areas remains extremely limited under both current and future scenarios. Therefore, protected-area establishment alone may be insufficient to conserve pollinator diversity and habitat.

Overall, the proportion of suitable pollinator habitat within protected areas is projected to increase, but this does not guarantee effective protection for bee species. If protected areas remain unchanged, the suitable habitat of pollinators will suffer substantial losses under climate warming (Martínez-López et al., 2021). This explains why the proportion of suitable habitat within protected areas is increasing. Although climate change is already so detrimental to pollinators, it is unrealistic to assume that protected areas will remain unchanged. Studies on the negative impacts of climate change on protected areas and forests are frequently reported (Hannah et al., 2007; Seidl et al., 2017). In addition to climatic factors, deforestation, land use change and pesticide abuse all pose significant challenges to vegetation within protected areas (Ashley et al., 2006). Even though certain species can retain ecological functions within protected areas, these areas are often not major agricultural production zones, limiting pollinators’ contribution to agriculture (Mcdonald and Boucher, 2011). Furthermore, our analysis shows that environmental suitability for pollinators is projected to decline in the future. In this context, it remains uncertain whether protected areas will continue to serve as refuges for pollinators.

Predictions indicate that the highly suitable habitat of honey bee species will be significantly fragmented in the future. While the conservation role of protected areas for species needs further study, establishing new nature reserves may enhance pollinator conservation (Martínez-López et al., 2021). In particular, the Himalayan region will become an important habitat for pollinators in the future. Although our projections suggest the Himalayan region is an important area, environmental problems such as habitat fragmentation, air pollution, and deforestation are becoming more pronounced in the region today (Chakraborty et al., 2017; Dhungel et al., 2018). A recent study showed that environmental changes in the Himalayan region have led to significant biomass loss, likely driven by global climate change and intensifying agricultural activities in the region (Chakraborty et al., 2018). However, the significant increase in the average vegetation index during the growing season may indicate improved food availability for pollinators in a warming climate (Shrestha et al., 2012). For the fragmented suitable habitat patches of pollinators in plains and low mountain areas, establishing nature reserves may not be the most effective strategy. For example, the Mekong Plain is the most densely populated area in Southeast Asia, and establishing nature reserves there would exacerbate human–land conflicts, thereby affecting agricultural production (Kontgis et al., 2019). In conclusion, protecting important species under climate change should avoid pure in-situ protection and adopt integrated landscape management strategies, a task that requires the active participation and efforts of land stakeholders. Furthermore, timely environmental education activities may benefit conservation efforts (González-Fernández et al., 2018). The floristic inventories allow us to determine which plants bees pollinate and from which plants they collect pollen and nectar. Conversely, we can grow plant species needed by bees in protected areas to support bee survival. Finally, the floristic inventories in the study area fill an obvious geographical gap in our sampling and lay the foundation for more scientifically and rationally predicting species distributions.