An Evaluation of Habitat Uses and Their Implications for the Conservation of the Chinese Bumblebee Bombus pyrosoma (Hymenoptera: Apidae)

Bumblebees are important pollinators for many wild plants and crops. However, the bumblebee populations are seriously declining in many parts of the world. Hence, the bumblebee conservation strategy should be urgently addressed, and the species distribution modeling approach can effectively evaluate the potentially suitable areas for their conservation. Here, one of the most abundant and endemic species of bumblebee in China, Bombus pyrosoma, was selected to assess current and future climates’ influence on its distribution with MaxEnt. Nine high-resolution bioclimatic/environmental variables with high contribution rates and low correlations were used. Four of the nine bioclimatic/environmental variables, min temperature of the coldest month (bio_06), annual mean temperature (bio_01), precipitation of wettest month (bio_13) and radiation of warmest quarter (bio_26), were found to be the most critical factors influencing the distribution of B. pyrosoma. The modeling results showed that the areas with high and moderate suitability for B. pyrosoma covered 141,858 and 186,198 km2 under the current climate conditions. More than 85% of the sampling sites in 2019 were found to be suitable under the current scenario. Under the future A1B and A2 scenarios in 2050 and 2100, the areas with low and moderate suitability for B. pyrosoma increased. However, alarmingly, the high suitability areas decreased under the future A1B and A2 scenarios in 2050 and 2100. Furthermore, regions covering seven provinces of northern China were the most crucial for developing nature reserves for B. pyrosoma, with the following order of suitable areas: Gansu, Shanxi, Ningxia, Qinghai, Shaanxi, Hebei and Beijing. Our study highlights the impact of future climate changes on the distribution of B. pyrosoma, and conservation strategies should mitigate the threats posed by environmental changes, particularly in the current high suitability areas.


INTRODUCTION
As important pollinators of many wild plants and crops, bumblebees are essential for maintaining the ecological balance and agricultural food production worldwide (Frankie et al., 2009;Arbetman et al., 2017;Bartomeus and Dicks, 2019). However, the distribution and abundance of bumblebees have decreased significantly due to habitat loss, environmental pollution, pathogens and pesticides, and non-native species invasion (Potts et al., 2010;Goulson et al., 2015;Cameron and Sadd, 2020). It has been reported that some bumblebee species have already exhibited apparent declines in Europe and North America, while the related information is limited in Asia (Cameron et al., 2011;Colla et al., 2012;Kerr et al., 2015). China has approximately 50% of the world's bumblebee species, representing the highest number of bumblebee species of all countries in the world Naeem et al., 2018). As China's terrain and climate are diverse, it is difficult to obtain detailed distribution information for bumblebee species and monitor their population dynamics (Peng et al., 2009;An et al., 2011). Therefore, evaluating the distribution of bumblebees in China is vital to their conservation, which will contribute to the sustainable development of ecosystem management and agricultural food production (Williams et al., 2012;Naeem et al., 2020).
The bumblebee Bombus pyrosoma Morawitz (Hymenoptera, Apidae) is an endemic and abundant species in China . This species is distributed from western Liaoning Province to eastern Qinghai Province and has been abundantly reported in Shanxi, Hebei and Gansu Provinces (An et al., 2008;Wu et al., 2009;Ma et al., 2011). Varieties of natural and agricultural plants have been visited by this species, such as Helianthus annuus L., Leonurus sibiricus L., and Medicago sativa L. (Ma et al., 2011). Due to the advantages of large colonies that are easy to rear, this species has become a potential candidate species for artificial rearing (Peng et al., 2009). Although the geographic distribution and differentiation of B. pyrosoma have been explored, detailed distribution models and information on future climate influences are still lacking (Ma et al., 2011;Naeem et al., 2018;Liu et al., 2020).
Due to the vast area and diverse geographical environment, it is difficult to conduct a large-scale evaluation of bumblebee distribution by traditional methods . Species distribution models (SDMs) have become a prominent tool to predict and evaluate species distributions (Pearson et al., 2007;Elith and Leathwick, 2009;Xu et al., 2012). Different types of SDMs are developed according to the ecological niche, among which maximum entropy (MaxEnt) is widely used (Elith et al., 2006;Wang et al., 2017). MaxEnt acquires the maximum entropy of species distribution and constructs an SDM to predict species' potential distribution based on the species distribution coordinates and environmental variables (Araujo et al., 2005;Phillips et al., 2006). MaxEnt is an effective tool for bumblebee prediction and evaluation in the EU, America, and Asia (Penado et al., 2016;Naeem et al., 2018;Naeem et al., 2019;Krechemer et al., 2020). In particular, the potential distribution of the invasive bumblebee species B. terrestris was researched in South America, Japan and China (Kadoya et al., 2009;Marshall et al., 2018;Naeem et al., 2018;Sirois-Delisle and Kerr, 2018;Aizen et al., 2019). Hence, MaxEnt will be used here to evaluate the distribution and conservation assessment of an endemic and abundant bumblebee, B. pyrosoma.
Species distribution modeling will help to improve the precision of conservation (Ma and Sun, 2018). Here, we aim to predict the distribution of B. pyrosoma and evaluate the effects of environmental variables on the habitat changes of this species. The selected current and future climatic variables will be used to assess the influence on the distribution of this species. The prediction and evaluation of the B. pyrosoma distribution will help to improve bumblebee conservation strategies.

Study Area and Species Data
The bumblebee species B. pyrosoma is an endemic species in China. A total of 6317 specimens of B. pyrosoma from 1905 to 2016 were collected from the museum of the Institute of Apicultural Research, Chinese Academy of Agriculture Science and the Institute of Zoology, Chinese Academy of Sciences. The specimens were recorded from seventeen provinces (Beijing, Chongqing, Gansu, Guizhou, Hebei, Henan, Hubei, Inner Mongolia, Liaoning, Ningxia, Qinghai, Shaanxi, Shanxi, Sichuan, Tianjin, Tibet, and Yunnan). Furthermore, in 2019, we collected 2119 B. pyrosoma individuals across its historical distribution in China. Specimens were sampled, and their location information was acquired with a handheld GPS (Garmin cs60).

Bioclimatic/Environmental Variables
Twenty-three bioclimatic/environmental variables, including nineteen bioclimatic variables at a 2.5 min resolution (∼5 km) and four environmental variables for the current period , were downloaded from the WorldClim-Global climate dataset 1 (Santana et al., 2019;Supplementary Table 1). The four environmental variables, related to bumblebee distribution, were growing degree days (30 arc min (∼50 km), Atlas of the Biosphere, The Nelson Institute Center for Sustainability and the Global Environment, University of Wisconsin-Madison) 2 , rainy_days (5 arc min (∼10 km), Climate Change Agriculture And Food Security, CCAFS) 3 , global reference evapotranspiration [Global-ET0, 30 arc sec (∼1 km), Consultative Group for International Agriculture Research (CGIAR)] 4 , and radiation of warmest quarter (bio_26) (10 arc min (∼20 km), Global climatologies for bioclimatic modeling, CliMond) 5 (Williams et al., 2014). The conversion tool in ArcGIS software was used to adjust the resolution of these four environmental variables similar to those of the 19 bioclimatic layers.
Future climatic data for 2050 and 2100 were downloaded from CliMond (Kriticos et al., 2012) based on the global climate model (GCM), CSIRO-Mk3.0 (Australia) with the Intergovernmental Panel on Climate Change (IPCC) IV SRES scenarios (A1B and A2) (Kriticos et al., 2014). The A1B scenario describes the balance between fossil and non-fossil resource use (Paterson et al., 2015). However, the A2 scenario depicts a changing world with high population growth but slow economic development and technological change (Paterson et al., 2015).
In the case of multicollinearity among those 23 bioclimatic/environmental variables, the Pearson correlation coefficient was used to select only those with low correlations (Wang et al., 2019). The Pearson correlation coefficient r of pairwise bioclimatic/environmental variables was calculated with the R project (Version 3.6.2). According to the contribution percentage, only those variables with a pairwise correlation efficiency less than 0.8 were retained. However, in the case of a correlation coefficient greater than 0.8, the variable with the highest contribution was selected (Casey et al., 2015;Naeem et al., 2018;Thapa et al., 2018;Mwakapeje et al., 2019).

Modeling Procedure and GIS Analyses
The maximum entropy model was used to evaluate the potential distribution of the species. First, we conducted an initial model with the thinned sample records and twentythree bioclimatic/environmental variables with MaxEnt (Version 3.4.1) 6 . In the initial model, we set "Random Test Percentage" to 25 and chose "Make pictures of predictions" and "Do Jackknife to measure variable importance." The rest of the model settings were the defaults. Second, the relevant variables were assessed using the abovementioned method in section Bioclimatic/Environmental Variables. Third, the distribution of B. pyrosoma was simulated with the thinned sample records and filtered environmental variables. Response curves were chosen, and the rest of the model settings were consistent with the initial model.
All the bioclimatic/environmental layers were converted into ASCII files from raster grid format using ArcToolbox 2.0 in ArcGIS (Brown et al., 2017). The species distribution data were also converted into comma-separated value (CSV) format in Excel. The predicted distributions of B. pyrosoma were assessed and reclassified into different categories. Furthermore, the distribution area was assessed in km 2 with zonal statistical analysis in ArcGIS.

Model Validation and Potential Habitat Prediction
The MaxEnt software is used to perform the Jackknife test to determine the relative importance and contribution of different variables in the distribution of species. Three approaches are used for this analysis: (1) Models created using all variables except excluding one, (2) models created using one variable, and all the other variables were excluded, (3) models created using all variables (Phillips et al., 2006;Santana et al., 2019). The response curves were calculated to determine the influence of climate change on the distribution of B. pyrosoma. These curves show how each environmental variable have an impact on the prediction of the MaxEnt model. Similarly, those response curves also show how the logistic prediction is changed with the environmental variable's change (Phillips et al., 2006). The receiver operating characteristic (ROC) curve and the area under the curve (AUC) values were used to evaluate the accuracy of the species distribution model (Wang et al., 2019). The theoretical value of AUC set between 0.5 and 1. When the value of AUC is close to 1, the prediction accuracy of the model is high. The model classification index was as follows: "failure" is 0.50 < AUC < 0.60, "poor" is 0.60 < AUC < 0.70, "fair" is 0.70 < AUC < 0.80, "good" is 0.80 < AUC < 0.90, and "excellent" is 0.90 < AUC < 1.00 (Araujo et al., 2005;Wang et al., 2017). The distribution classification was divided into four categories: (I) unsuitable (0-0.50), (II) low suitability (0.50-0.75), (III) moderate suitability (0.75-0.90), and (IV) high suitability (0.90-1.00).

Spatial Conservation Assessment
We compared the historical, current, and future distribution models with the sampling sites of B. pyrosoma in 2019. The distribution modeling accuracy and fitting degree between the models of different periods and all sampling sites were also assessed. We combined the current distribution of B. pyrosoma with the four future potential distributions in 2050 and 2100 (A1B and A2 scenarios). The raster calculator was used to calculate the suitable overlapping areas (Naeem et al., 2018;Zhang et al., 2018). The distribution area and proportion of different categories of habitat suitability were determined and assigned to China's corresponding administrative division. The suitable area of each province was sorted to indicate the priority of conservation.

Environmental Variable Selection and Model Performance
Pearson correlation of pairwise variables showed that the coefficients between 28 pairs of the 23 bioclimatic/environmental variables were greater than 0.8. Twenty-three bioclimatic/environmental variables were sorted according to the initial model's percentage contribution (Table 1). Nine variables were selected based on the Pearson correlation coefficient and contribution rate. The selected variables were annual mean temperature (bio_01), mean diurnal range {bio_02 [mean of monthly (max temp-min temp)]}, isothermality [bio_03 (bio_02/bio_07) (× 100)], maximum temperature of warmest month (bio_05), minimum temperature of coldest month (bio_06), annual precipitation (bio_12) precipitation of wettest month (bio_13), precipitation seasonality [bio_15 (coefficient of variation)], and radiation of warmest quarter (bio_26). The value of the Pearson correlation coefficient between every two variables was < 0.8, and the cumulative contribution rate of the nine remaining variables was 100% ( Figure 1A). The curve of omission on the training samples exhibited a similar trend as that on the test samples. The AUC values of the training data and test data were 0.952 and 0.941, respectively, in the current distribution model (Figures 1B,C). Therefore, the precision of the current distribution model was excellent according to the model evaluation index.
Relationship Between the Occurrence of B. pyrosoma and Bioclimatic/Environmental Variables Jackknife analysis showed that the minimum temperature of the coldest month (bio_06), precipitation of the wettest month (bio_13), annual mean temperature (bio_01) and annual precipitation (bio_12) were the most important environmental variables that affected the distribution of B. pyrosoma, with a training gain of more than 0.6 (Figure 2A). The radiation of the warmest quarter (bio_26) and the maximum temperature of the warmest month (bio_05) were also important bioclimatic/environmental variables, with training gains of more than 0.45. The contribution of the last three variables (bio_02, bio_03, bio_15) was more than 0.14. The response curve thresholds of the nine bioclimatic/environmental variables were obtained ( Figure 2B). Where the presence probability was more than 0.5, the bioclimatic/environmental variable threshold values were as follows: min temperature of the coldest month (bio_06) ranged from −16.78 to −4.10 • C, precipitation of wettest month (bio_13) ranged from 81.33 to 151.65 mm, annual mean temperature (bio_01) ranged from 3.16 to 11.57 • C, annual precipitation (bio_12) ranged from 353.75 to 764.47 mm, radiation of warmest quarter (bio_26) ranged from 181.00 to 208.00 W/m 2 , the maximum temperature of the warmest month (bio_05) ranged from 17.91 to 29.36 • C, the mean diurnal range (bio_02) ranged from 8.95 to 13.07 • C, isothermality (bio_03) ranged from 26.79 to 34.25, and precipitation seasonality (bio_15) ranged from 70.10 to 123.44.

Modeled Distribution of B. pyrosoma
The current distribution of B. pyrosoma shows that the suitable habitat is concentrated from central to northern China (Figure 3). The high suitability areas were located in eastern Qinghai, southern Gansu and southern Ningxia, parts of Shaanxi and Shanxi provinces. Moderate suitability was located in eastern Qinghai, southern Gansu, southern Ningxia, northern and southern Shaanxi, and parts of Shanxi, Hebei, Beijing, Inner Mongolia, Hubei, Hunan, and Chongqing provinces. The low suitability areas were located in eastern Inner Mongolia, western Liaoning, southern Hebei, parts of Shanxi, north and south of Shaanxi, southern Ningxia and southern Gansu, eastern Qinghai, parts of Sichuan, Chongqing, Hunan and Hubei provinces. According to the percentage of the distinct categories of suitability, the high, moderate, and low suitability areas represented 1.48% (141,858 km 2 ), 1.95% (186,198 km 2 ), and 2.72% (259,878 km 2 ) of the total area of China, respectively.

Conservation Assessment
The dynamics of the suitable habitat for B. pyrosoma were obtained (Figure 4). The total suitability area trend decreased, then increased, and then slightly decreased from the 1900s to 2100s. The suitability area of the current distribution decreased by 189,375 km 2 and increased by 150,365 km 2 compared to the historical distribution. Under the A1B scenario, the low suitability area increased by 95,400 km 2 in 2050 and 80,400 km 2 in 2100, and the moderate suitability area increased by 12,135 km 2 in 2050 5,191 km 2 in 2100 compared with the current suitability area. However, the high suitability area decreased by 44,636 km 2 in 2050 and 36,580 km 2 in 2100. Under the A2 scenario, the low suitable area increased by 93,178 km 2 in 2050 and 83,178 km 2 in 2100, and the moderate suitable area increased by 13,524 km 2 in 2050 and 12,413 km 2 in 2100 compared with the current suitability area. However, the high suitability area decreased by 42,414 km 2 in 2050 and 34,636 km 2 in 2100.
The overlapping suitability areas were also calculated by overlapping the current and future (2050 and 2100 A1B and A2 scenarios) potential distributions of B. pyrosoma (Figure 5). The total suitable area was 200,278 km 2 , which accounted for 2.10% of China's total area. The areas of high, moderate, and low suitability were 42,778, 53,889, and 103,611 km 2 , respectively. Seven provinces were identified as important for the conservation of B. pyrosoma ( Table 2). Gansu Province inhabited the most extensive distribution of high suitability areas, with 26,389 km 2 , which is more than six percent of the total area of the province. The second-largest high suitability area was 8,056 km 2 (5.05%) in Shanxi. The Ningxia Hui Autonomous Region was the third most important area, with 4,722 km 2 (8.97%). The highly suitable areas in Qinghai covered 2,222 km 2 (0.31%), while Shaanxi covered 1,389 km 2 (0.68%). Although there were no high suitability areas in Hebei and Beijing, the percentages of . The symbols of "***", "**," and "*" represent significance level of p-value "0-0.001," "0.001-0.01," and "0.01-0.05". moderate and low suitability were large, covering 39,722 and 3,611 km 2 , respectively.

Validation of the Distribution Model
Comparing sample records in 2019 and the modeled habitat showed a better fitness of the current and future habitat models than the historical habitat model. There were 37 sampling sites in 2019 (37.76% of all sampling sites) located in the unsuitable areas, 21 sampling sites (21.43%) in the low suitability areas, 20 sampling sites (20.41%) in the moderate suitability areas and 20 sampling sites (20.41%) in the high suitability areas under the historical habitat model (Figure 6). Under the current habitat model, there were 14 sampling sites (14.29%) located in the unsuitable areas, 19 sampling sites (19.39%) in the low suitability areas, 31 sampling sites (31.63%) in the moderate suitability areas and 34 sampling sites (34.69%) in the high suitability areas. Under the future habitat model, the fitness of the models was similar in all four scenarios (2050 and 2100 A1B and A2). Specifically, 5-8 sampling sites (5.10-8.16%) were located in the unsuitable areas, 18-24 sampling sites (18.37-24.49%) were located in the low suitability areas, 25-26 sampling sites (25.51-26.02%) were located in the moderate suitability areas, and 25-28 sampling sites (25.51-28.57%) were located in the high suitability areas.

DISCUSSION
Bumblebees contribute significantly to maintaining the natural ecological balance and crop production (Brown, 2011;Cameron and Sadd, 2020). To conserve these important pollinators, their potential distribution has often been evaluated by modeling tools (Elith et al., 2006;Naeem et al., 2018). In this study, the current distribution of the bumblebee B. pyrosoma was modeled based on historical sample records and validated with samples collected in 2019. Furthermore, the potential distribution of B. pyrosoma was simulated under future climate scenarios. The potential distribution and future suitability dynamics of B. pyrosoma will help to develop a conservation strategy.
Nine bioclimatic/environmental variables were selected in this study, and they are slightly different from the variables selected in previous studies (Casey et al., 2015;Naeem et al., 2019;Krechemer et al., 2020). The difference in the variable selection may be due to the different bumblebee species and their sampling sites. Certain variables play an important role  in the distribution of certain bumblebees (Casey et al., 2015;Naeem et al., 2018). Currently, the minimum temperature of the coldest month (bio_06), precipitation of the wettest month (bio_13) and annual mean temperature (bio_01) made significant contributions to the model. Such conditions are present in the meadow of Wutai Mountain and Genting Mountain in North China, which also have grasslands with abundant flowers (An et al., 2008;Richardson et al., 2018;Naeem et al., 2019). Other than bioclimatic/environmental conditions, the abundance of foraging plants together with the cold environment and humidity will benefit bumblebees ( Table 1). Most bumblebees are concentrated in temperate and subfrigid zones in the Northern Hemisphere, and these species are suitable for cold and humid climates (Penado et al., 2016;Streinzer et al., 2019).
Furthermore, the contribution of radiation in the warmest quarter (bio_26) cannot be ignored. When it was omitted, the radiation of the warmest quarter was the environmental variable that decreased the gain the most; therefore, this variable had the most information that was not available in the other variables (Figure 2A). The results are consistent with those of Williams's study since extreme values reduce food-plant pollen and nectar production (Williams et al., 2014). The modeled distribution of B. pyrosoma was mainly distributed in the north of the Qinling Mountains and south of the Inner Mongolia Plateau in northern China Ma et al., 2011;Naeem et al., 2018). The suitable areas (high and moderate suitability collectively) covered 328,605 km 2 (3.43%) under the current climate scenario, 297,361 km 2 (3.09%) in 2050 on average, and 301,250 km 2 (3.13%) in 2100 on average. The suitable area found in our current study was smaller than that found in a previous study (844,057 km 2 ) (Naeem et al., 2018). The difference may be due to the use of thinned species records with different distances in our study. We thinned the records by distances of at least 30 km, whereas previously, the samples were thinned by distances of less than 10 km. Additionally, different environmental variables may lead to different model results. We used nine bioclimatic and environmental variables, while Naeem et al. (2019) used eight variables, five of which were the same in our study. The remaining three, temperature seasonality, precipitation of driest month and elevation, were different variables.
A1B and A2 are two different scenarios for the future climate. The future projection years of 2050 and 2100 were selected to provide a reasonable snapshot of two time periods to represent the near future and distant future. Under these two scenarios, the suitable habitat of B. pyrosoma will increase in the future. However, the high suitability area will decrease under both future scenarios (Figure 4). These trends were consistent with those found in the previous study. However, the fluctuations in the previous study were stronger than those found in ours. This difference may be due to the climatic scenarios and the sampling size. Naeem et al. (2019) used different climatic scenarios with two levels of representative concentration pathways (R) (2050s-RCP2.6, 2050s-RCP8.5, 2070s-RCP2.6, and 2070s-RCP8.5) from the Coupled Model Intercomparison Project Phase 5 (CMIP5) in the 5th assessment report of the IPCC. Besides, more than 700 sampling sites in our study were one hundred more sampling sites than in the previous study (Naeem et al., 2019). Furthermore, the B. pyrosoma population dynamics are similar to those in other reports of bumblebees, who also suffer from an existential crisis in the future (Kerr et al., 2015;Sirois-Delisle and Kerr, 2018).
Finally, the suitable areas for B. pyrosoma will extend in the future compared with the current distribution. However, high suitability areas will be lost, which cannot be neglected. The shrinkage in the high suitability area may be caused by future climate changes, agrochemicals, parasites and pathogenic factors Bartomeus and Dicks, 2019). In addition, it is expected that the interaction between climate change and land use, including agricultural lands, developed land, cropland, wetlands, or in the form of crop and grazing lands will exacerbate species range losses (Newbold, 2018;Sirois-Delisle and Kerr, 2018). We should also pay attention to the influence of those factors on the distribution of the species. In this study, the overlaid potential distribution under the current and four future climate scenarios improves the accuracy of conservation area identification. According to the overlaying habitat, conservation of B. pyrosoma should be established in the following administrative divisions: Gansu, Shanxi, Ningxia Qinghai, Shaanxi, Hebei, and Beijing provinces ( Table 2). Gansu, Shanxi, and Ningxia provinces have the most significant proportions of high suitability areas, accounting for more than 5% of the province's total area. These areas are recommended to focus on for primary conservation.
Our study evaluated the consistency between the distribution models of B. pysoroma and the resampling in 2019. More than 85% of the resampling sites were located in potentially suitable areas, and more than 40% of the resampling sites were highly suitable. More than 70% of the resampling sites were situated in future predicted distribution areas. Furthermore, we optimized the bioclimatic/environmental variables and thinned the sample sites to improve modeling accuracy. Through the data screening, modeling verification and analyzing previous studies, the species distribution model established by MaxEnt software was reliable (Elith et al., 2006;Naeem et al., 2018;Mwakapeje et al., 2019;Naeem et al., 2019).

CONCLUSION
In conclusion, the potential distribution of an endemic and abundant bumblebee species in China, B. pyrosoma, was modeled under the past, current, 2050 and 2100 climates (A1B and A2 scenarios). The temperature, precipitation and radiation of the warmest quarter were recognized as the most significant bioclimatic/environmental variables that affect bumblebees' distribution. With climate change, the main potential distribution areas may increase, while the high suitability areas will decrease. Therefore, it is essential to develop a conservation strategy according to the suitability dynamics for bumblebees. Bioclimatic/environmental factors, potential distribution and conservation strategies should be combined to protect bumblebees in Asia.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.