Predicting negative Effects of Climate Change on Taiwan’s endemic Bumblebee Bombus formosellus

Introduction: Climate change is threatening biodiversity hotspots. Notably, alpine bumblebees, which are mostly associated with a cold ecological niche, face a higher risk of extinction. Bombus formosellus is one of the alpine bumblebees endemic to Taiwan. Aims/Methods: In this study, we use ensemble ecological niche modeling for the first time to predict spatial and temporal dynamics for habitats suitable to B. formosellus under current and future climate scenarios (representative concentration pathway, RCP 2.6, 4.5, and 8.5 in the 2070s). Results: This model identified that the cool temperature with low variation, a specific range of precipitation and presence of coniferous forest and grasslands were the key factors affecting the distribution of B. formosellus. Using modeling to predict suitable habitats under various scenarios, we discovered that, compared with the current climatic conditions, the predicted suitable habitat area in the future decreased regardless of which climate change scenario was applied. In particular, RCP 8.5 appeared to be the most significant, with an area loss of nearly 87%, and fragmentation of the landscape with poor connection. Discussion: In summary, our analyses indicate that cool environments are suitable for B. formosellus. However, Taiwan’s warming is more significant in the high mountains than in the plains. The climate change trajectory may become increasingly unfavorable to B. formosellus. Consequently, this species may face the risk of extinction in the future. Implications for insect conservation: We predict that many suitable habitats of B. formosellus will disappear or become fragmented in the future. Therefore, the remaining patches have become important refuges, and protection measures in these areas should be strengthened.


Introduction
The worldwide distribution and populations of bumblebees reduced considerably in the past several decades, and many species in Europe, North America, and East Asia face extinction . The main threats include insecticide contact, habitat loss, pathogen infection, invasive Climate change is a severe global problem. Due to overreliance on fossil fuels, the concentration of greenhouse gases in the atmosphere and the global average temperature have increased significantly. The Intergovernmental Panel on Climate Change (IPCC) predicts that by the end of the 21st century (2081-2100), global average temperatures will have increased from between 0.3 and 4.8 °C (IPCC 2013(IPCC ), compared to 1986(IPCC -2005. Taiwan is an East Asian Pacific island heavily impacted by climate change. The temperature of Taiwan rose by approximately 1.4 °C between 1911 and 2005, which was twice as the average of the Northern Hemisphere overall (~ 0.7 °C) (IPCC 2007;Shiu et al. 2009). It is estimated that Taiwan's temperature will continue to rise, perhaps by more than 3 °C by the end of the 21st century (Lin et al. 2015).
The study of species distribution is a prerequisite for drafting conservation policies and may help in alleviating the impact of global climate change on species distribution and range boundaries (Thomas 2010). Ecological niche modeling (ENM) integrates species occurrence records and environmental variables to construct mathematic models for predicting habitat suitability, and these models can be predict spatial and temporal dynamics of species for future climate change scenarios (Elith and Leathwick 2009). To improve the accuracy and reliability of ENM, "ensemble" ecological niche modeling (EENM) is widely used, using different algorithms to construct multiple models and combining the predictions of each model to achieve a collective consensus (Araújo and New 2007). Recently, it has been used to study the impact of climate change on bumblebees and other bee species (Pradervand et al. 2014;Martins et al. 2015;Elias et al. 2017;Dew et al. 2019;Giannini et al. 2020).
Although Taiwan is located between the tropics and subtropics, and therefore representing a connectivity zone for the regional biodiversity, its terrain is steep and mountainous, and the temperate climate of high altitude areas provides suitable habitats for bumblebees. Nine species of bumblebees have been identified on the island to date, including Bombus formosellus (Frison 1934), which is found on both sides of the Central Mountain Range but is unknown in almost all the southern half. Though other details are lacking, it is an endemic and relatively rare species (Starr 1992;Williams et al. 2020). Taiwan's island location limits species migration, and with the continued global warming, the habitat suitable for alpine bumblebees may be limited. This is particularly true for endemic species already been confined to a specific geographic areas. Therefore, it is critical to study the impact of climate change on the distribution of B. formosellus.
This study assesses future threats to B. formosellus based on the species survey data and EENM. The investigation seeks to determine: (1) what are major climate and landscape factors affecting species distribution; (2) how does the suitable habitat change under different climate change scenarios; and (3) if the landscape structure of habitat changes, how will the species be affected?

Study area
Taiwan lies on the boundary between Eurasian Continent and the Pacific Ocean, separated from mainland China by the Taiwan Strait to the west (Fig. 1a). The land area is approximately 36,000 km 2 , and the human population on the island is nearly 23 million. The collision of the Eurasian Plate and the Filipino Plate makes Taiwan a mountainous island, with mountains (elevation ≥ 1000 m asl) covering approximately 32% of the island. More than 200 mountains have an elevation greater than 3,000 m, and the highest peak has an elevation of 3952 m (Guan et al. 2009) (Fig. 1b). The annual average temperature is about 21 °C, and annual precipitation is around 2,500 mm. The Tropic of Cancer divides the island into two climate zones: tropical monsoon to the south and subtropical monsoon to the north (Central Weather Bureau 2009). At higher elevations, the central inner alpine zone has a temperate and cool environment similar to temperate zones. The land cover is predominantly forest, accounting for roughly 61% of the island (Forestry Bureau 2014), with urban and agricultural land areas mainly located on the coastal plains.

Species occurrence record and environment data
Occurrence records of bumblebees (genus Bombus) were collected from databases of the Forest Arthropod Collection of Taiwan (https://fact.tfri.gov.tw/), the Global Biodiversity Information Facility (https://www.gbif.org/) and previously published literature (Starr 1992). To maintain precise geographic coordination, we data with less than three decimal digits were removed. Further analysis was conducted on the R 3.6.3 platform, using package "spThin" for spatial thinning (Aiello-Lammens et al. 2015). To ensure consistency with resolution of the predictors, the neighboring distance between occurrence records had at least 1 km radius, reducing the impact of spatial autocorrelation and geographic sampling bias. After these treatments, field surveys of species occurrence were carried out by our research team from March 2020 to February 2021 to confirm these records. Finally, we retained a total of 203 unique occurrence records for 6 species, including 33 occurrence records of B. formosellus. (Fig. 1b).
For the current climate scenario, we extracted climatic data for 19 bioclimatic variables with 30 s spatial resolution from the WorldClim database (Hijmans et al. 2005), and all layers were re-sampled at a 1 × 1 km grid. We also collected land cover vector data from the National Land Surveying and Mapping Center (https://www.nlsc.gov.tw/), including six land cover categories: coniferous forest, broadleaf forest, mixed forest, grasslands, built-up areas and agricultural lands. Data were transformed to continuous variables by calculating the relative proportion of each land cover category within 1 × 1 km. The R package "usdm" (Naimi et al. 2014) was used to calculate the variance inflation factor (VIF) to overcome multicollinearity, retaining only the predictive variables with VIF < 3.

Climate change scenario data
For future scenarios, land cover categories were not included because the rate of land use transitions in Taiwan has slowed down over the past three decades (Chen et al. 2019). Hence, we assumed that climate factors would follow the climate change scenario "representative concentration pathway (RCP)," published in the 5th IPCC report, with other variables remaining constant. We chose three carbon-emission scenarios: RCP 2.6, 4.5, and 8.5, with each representing the radiative forcing to increase by 2.6, 4.5, and 8.5 Wm − 2 between 1750 and 2100, respectively. The perspective ranges of change in average temperature were 0.3-1.7 °C, 1.1-2.6 °C, and 2.6-4.8 °C, corresponding to the mild, moderate, and severe warming scenarios (IPCC 2013). The WorldClim database also provided future climate data, which were first interpolated from the current actual observations to obtain the baseline climate data. Climate data were downscaled and corrected with future climate data estimated by the Coupled Model Intercomparison Project Phase 5 (Taylor et al. 2012) to generate a climate map with 1 × 1 km resolution. We downloaded the bioclimatic data of 2070s (mean value of 2061-2080) under the three scenarios, and chose only the global climate models (GCM) suitable for predicting Taiwan's climate, including CESM1-CAM5, CSIRO-Mk3-6-0, GISS-E2-R, and MIROC5 (Lin and Tung 2017). The future bioclimatic variables generated by each of these four GCMs were used to construct models for B. formosellus, and the average was calculated as the final  (0) habitats. Such maps can be used to observe changes in the range of suitable habitats at current and in the future. The importance of each predictive variable for model fitting was evaluated by the "getVarImp" function of the R package "sdm". In addition, the response curves revealed the relationship between the species occurrence probability and the variables.

Spatial pattern analysis
We used the Fragstats software (version 4.2) to analyze the spatial pattern of suitable habitats for B. formosellus. The splitting index (SI) was used to understand the status of fragmentation, with higher SI values indicating more severely fragmented habitat patches (Jaeger 2000). The proximity index (PI) evaluated connectivity of habitat patches, where a high PI value indicated the presence of large and adjacent patches surrounding the target patch to maintain a reasonably good connection between populations (Gustafson and Parker 1994). The index required a searching radius, which was at 2 km based on the forage distance of bumblebees, currently determined as 1-3 km (Kreyer et al. 2004;Osborne et al. 2008;Hagen et al. 2011).

Results
Under the criterion of VIF < 3, the 25 predictive variables we obtained earlier were screened down to nine variables for constructing the prediction models of B. formosellus. Table 1 shows the average AUC and TSS calculated from the test dataset by eight algorithms, with 50 replicates each. Overall, BRT, Maxent, MLP, RF and SVM had excellent accuracy; of these, RF showed the best performance, indicating that our models had high prediction power. Only 248 high-performance models (TSS > 0.8) were retained for the following weighted averaging to ensure highly accurate predictions.
Considering the importance of each variable to the model fitting, climate variables were clearly more important than land cover variables (Table 2). Mean temperature of the driest quarter was the most important variable result to evaluate habitat suitability. This method eliminated uncertainty among the GCMs (Araújo and New 2007).

Ensemble ecological niche modeling
The R package "sdm" (Naimi and Araújo 2016) was used to execute EENM to predict a suitable habitat for B. formosellus. Eight algorithms were applied, including generalized additive model (GLM), generalized additive model (GAM), multivariate adaptive regression splines (MARS), boosted regression tree (BRT), random forest (RF), maximum entropy (Maxent), multi-layer perceptron (MLP), and support vector machine (SVM). Since all algorithms required background data (pseudo-absence points), pseudo-absences could only be selected from areas where other bumblebee species had been recorded, a more objective approach to avoid considering under-sampled areas as unsuitable for B. formosellus (Marshall et al. 2018;Ghisbain et al. 2020). From the occurrence records of B. formosellus 70% were randomly chosen as the training dataset, and the remaining 30% were used as the test dataset. The model performance of each algorithm was evaluated by bootstrap sampling. The 50 replicates constructed with 8 algorithms yielded 400 models.
The area under the curve (AUC) generated by the receiver operating characteristic was used to evaluate accuracy, where the range of AUC was 0.5-1 and higher values indicated better prediction accuracy of the model. The true skill statistics (TSS), were calculated for the range of TSS between + 1 and − 1, where values approaching 1 indicated a near-perfect model (Allouche et al. 2006). The standard for excellent models was set as values of AUC and TSS higher than 0.9 and 0.8, respectively (González-Ferreras et al. 2016). The ensemble method was performed by retaining models with TSS > 0.8, and then calculating the weighted average based on their TSS values to generate an ensembled occurrence probability map of species. On the map, greater the probabilities, indicate more likelihood for a suitable habitat for the species. Therefore, the map can indicate habitat suitability. The probability corresponding to the maximum TSS was considered the threshold (Liu et al. 2013) to convert the map into a binary map of suitable (1) and unsuitable  Table 1 The area under the curve for the receiver operating characteristic (AUC) and true skill statistic (TSS) generated by the generalized additive model (GLM), generalized additive model (GAM), multivariate adaptive regression splines (MARS), boosted regression tree (BRT), random forest (RF), maximum entropy (Maxent), multi-layer perceptron (MLP) and support vector machine (SVM) This difference was especially clear with RCP 8.5 (Fig. 4a), while the PI showed an inverse trend with more warming associated with less habitat connectivity (Fig. 4b). The overall spatial pattern suggested that climate change may make B. formosellus habitat more fragmentary and widely separated.

Discussion
EENM is an effective tool for predicting the impact of climate change on the suitable habitats of species, and an ensemble of model types improves prediction accuracy (Araújo and New 2007;Marmion et al. 2009). The present study shows that such modeling yielded good results in predicting the distribution of B. formosellus. In terms of single algorithms, this study found that BRT, Maxent, MLP, RF and SVM showed excellent performance, with RF being the best (Table 1). This result was similar to previous studies investigating other bee species (Biella et al. 2017;Elias et al. 2017).
These models indicated that climate change variables strongly impacted the prediction of suitable B. formosellus habitat. Among the temperature-related variables, mean temperature of driest quarter was most important ( Table 2). The main reason for this is that the relatively sparse rainfall season in Taiwan usually occurs in winter, consistent with the general perception that alpine bumblebees are suitable for relatively cooler habitats. We also noticed that mean temperature of driest quarter higher than 10 °C was unsuitable (Fig. 2a), consistent with alpine bumblebees living in other regions having a narrow ecological amplitude (Penado et al. 2016;Suzuki-Ohno et al. 2017). In addition, temperature seasonality is related to the presence of species, thus B. formosellus may be more associated with cooler habitats with lower temperature variation. Other critical variables of importance would be precipitation of the driest wettest months, which should be associated with the interaction with plants. Plant rehydration increases nectar production, related to the species' occurrence in all algorithms and ensemble prediction (EP). Based on EP, the variables with an importance > 10% included mean temperature of driest quarter (39.85%), precipitation of driest month (20.29%), temperature seasonality (16.59%), precipitation of wettest month (14.11%), coniferous forest (10.17%) and grassland (10.03%).
The primary contribution variables from EP (Fig. 2) indicate that the probability of occurrence was highest at lowest mean temperature of driest quarter, temperature seasonality and precipitation of wettest month, and there was a strong decrease onwards. In contrast, precipitation of driest month was positively correlated with habitat suitability. Analysis of the above response curves indicated that B. formosellus not only prefers a cool environment, it but also requires a niche with a specific range of precipitation. For land cover variables, both coniferous forest and grassland show a trend of higher probability with higher cover.
The constructed ensemble models were used to simulate suitable habitats of B. formosellus, showing that such habitats are generally at higher elevations. The current suitable habitat area was estimated to be 4954 km 2 , which decreased substantially in the climate change scenario for the 2070s, indicating that habitat suitability decreases with climate warming. The estimated suitable habitat areas for RCP2.6, RCP4.5, and RCP8.5 were 2566, 2418, and 668 km 2 , respectively, decreases of 41%, 51%, and 87%, respectively, over the current scenario. The spatial range showed that, although suitable habitats may be reduced, large patches could remain in the northern half under the RCP2.6 and RCP4.5 scenarios. However, for the RCP8.5 scenario, the overall suitable habitat area was reduced dramatically, and almost disappeared in the southern half (Fig. 3).
These results suggested that the spatial pattern for suitable habitats of B. formosellus will change under various future climate change scenarios. We calculated the SI and noted that future climate change would lead to increased habitat fragmentation, with the SI being higher in the 2070s than in the current condition, regardless of the scenario. of Taiwan, the results of all GCM concurred that Taiwan is heavily affected by global warming since average temperature in all predictions showed an upward trajectory (Hsu and Chen 2002;Lin et al. 2015. The spatial pattern suggested that warming is more significant in high mountains than in plains, profoundly impacting ecosystems and species distributions (Lin et al. 2015). In other areas, climatic warming has negatively affected alpine bumblebees with restricted geographic ranges (Rasmont et al. 2015). Our prediction models show that suitable habitats of B. formosellus would reduce significantly in climate warming scenarios. The situation is particularly dire under the severe warming scenario presented by RCP8.5 (Fig. 3), with the few remaining habitat patches becoming the only refuges. Aside from the risk of vast area loss of suitable B. formosellus habitats, predictions for the RCP8.5 scenario also revealed that the remaining suitable habitat patches would be fragmentary with poor connectivity (Fig. 4). Habitat fragmentation is one of the major causes that species becoming threatened. If the habitat breaks into small which in turn increases the rate of bumblebee visits (Burkle and Runyon 2016;Descamps et al. 2018;Höfer et al. 2021). This suggests that although dry and cool environments are suitable for B. formosellus, appropriate rainfall is still required. The results suggest that optimal precipitation is 180-250 mm ( Fig. 2b and d). The Bombus alpinus of the European Alps also falls into a similar niche (Biella et al. 2017). Furthermore, coniferous forest and grasslands are the main land cover categories at high altitudes in Taiwan, which beneficial foraging and nesting habitats for bumblebees (Liczner and Colla 2019;Clake et al. 2022). Thus, the coverage of these two land cover categories was also influential ( Fig. 2e and f).
Bumblebees are sensitive to climate change, making them an effective environmental indicator (Herrera et al. 2018). Taiwan has been experiencing climate warming in the past century (Hsu and Chen 2002), with days of drought increasing in many regions ). Accordingly, the climate of Taiwan will become increasingly less favorable for B. formosellus. In predicting climate change Fig. 2 Response curves of the primary contribution variables according to the ensemble prediction due to their stronger ability to fly and spread (Greenleaf et al. 2007;Wright et al. 2015). However, B. formosellus is a small-body species among the local bumblebees (Starr 1992), so the predicted habitat loss and fragmentation will be severe challenges.
Pollination services provided by bumblebees are an essential component for the maintenance of biodiversity in alpine environments (Bingham and Orthner 1998). Taiwan may be and scattered patches, population size and genetic diversity would reduce, along with the species' ability to adapt to environmental changes (Haddad et al. 2015). Mounting evidence also indicates that habitat fragmentation may directly affect bumblebee's gene flow, reproduction, and interaction with plants (Xiao et al. 2016;Jackson et al. 2018;Maurer et al. 2020). It is also noteworthy that large-body bee species have a better chance of survival in a fragmentary landscape in mountain regions covered by many nature reserves (Tang et al. 2006). However, relying on the in situ conservation one of the potential biodiversity hotspots in the world due to its species-rich ecosystems (Myers et al. 2000), particularly alone is not sufficient to protect B. formosellus. Our projections found a critical impact on B. formosellus under the severe warming scenario, but effects on this species were less severe under the mild or moderate warming scenarios. Consequently, the key to genuinely solving the core problem lies in whether humanity can cooperate to adopt a low-carbon economy and lifestyle. Many anthropogenic activities with high carbon footprints remain around the mountain areas of Taiwan, such as intensive orchard management and flower-viewing tourism, which may accelerate warming of the suitable habitat for B. formosellus. Therefore, we suggest that the possible impacts should be assessed to further mitigate these human disturbances.

Data Availability Available from authors upon request.
Code Availability Program R code used for this study is available upon request.

Conflicts of interest/Competing interests
The authors have no relevant financial or non-financial interests to disclose.

Ethics approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.