Identifying conservation priority areas for North American bumble bee species in Canada under current and future climate scenarios

Many bumble bee species are declining globally from multiple threats including climate change. Identifying conservation priority areas with a changing climate will be important for conserving bumble bee species. Using systematic conservation planning, we identified priority areas for 44 bumble bee species in Canada under current and projected climates (year 2050). Conservation priority areas were identified as those that contained targeted amounts of each species predicted occurrence through climate envelope models, while minimizing the area cost of conserving the identified conservation priority areas. Conservation priority areas in the two periods were compared to established protected areas and land cover types to determine the area of current and future priority sites that are protected and the types of landscapes within priority areas. Notably, conservation priority areas were rarely within established protected areas. Priority areas were most often in croplands and grasslands, mainly within the mountain west, central and Southern Ontario, Northern Quebec, and Atlantic Canada under all climate scenarios. Conservation priority areas are predicted to increase in elevation and latitude with climate change. Our findings identify the most important regions in Canada for conserving bumble bee species under current and future climates including consistently selected future sites.


| INTRODUCTION
Reversing biodiversity loss requires recovering species at risk of extinction by protecting species from on-going threats. Efforts to protect species at risk in Canada have generally been unsuccessful with only a minority of species' statuses improving following listing, while most remain unchanged (Favaro et al., 2014). This lack of effective conservation could be because current protected areas (e.g., parks and conservation reserves) were not designed with the intention of conserving biodiversity, but often were created in regions unsuitable for or far away from development (Bolliger et al., 2020;Kerr & Deguise, 2004).
Systematic conservation planning can be used to solve the issue of ineffective protected area design by identifying conservation priority areas (Margules & Pressey, 2000;Schwartz et al., 2018). Conservation priority areas are locations that would effectively and efficiently conserve species through the establishment of new protected areas or other targeted conservation actions (Margules & Pressey, 2000). The conservation objectives (i.e., species presence, genetic diversity, etc.) and targets (i.e., how much area or how many species, etc.) is determined by the use of the prioritization models. The models efficiently identify conservation solutions by considering the cost of conserving the priority areas. Examples of cost can be the price of purchasing land and the size of the priority area, since larger areas would have higher costs to purchase, maintain and be more difficult to acquire, etc., the cost of restoration and/or loss of income from the land (Ball et al., 2009;Naidoo et al., 2006). This allows managers to adapt models to the objectives, targets, and costs relevant to their targeted species or group of species in need of conservation.
Bumble bees (Apidae: Bombus spp.) are currently experiencing multiple threats in North America including climate change, habitat loss, pesticides, pathogen spillover, and competition from non-native species Szabo et al., 2012;Williams et al., 2009Williams et al., , 2014. These threats have led to the declines of 26% of the approximately 50 North American bumble bee species that were assessed by the IUCN Red List (Cameron & Sadd, 2020) and are need of conservation action. Despite being a relatively well studied insects, there are still some knowledge gaps including the range and occurrences of bumble bee species. Climate envelope models (i.e., environmental niche models) can be used to predict the suitability of a landscape by comparing climatic conditions of known occurrences to areas where occurrence is uncertain. Previous research has shown that climate change has already contributed to observed range declines and is predicted to continue to have a negative impact on bumble bees (Lee et al., 2019;Martínez-L opez et al., 2021;Martins et al., 2015;Pyke et al., 2016;Sirois-Delisle & Kerr, 2018). Despite southern range contractions noticed for many North American bumble bee species, there has not been an equivalent northern expansion (Kerr et al., 2015). Accordingly, conservation actions need to consider how to mitigate the negative impacts of climate change. Systematic conservation planning can be used to identify conservation priority areas for bumble bees under multiple climate scenarios to help conserve declining species now and into the future by combining current and future climate envelope models. Bumble bees are habitat generalists and can be found in many different land cover types. Once conservation priority areas are identified the types of land covers within the conservation priority areas can be identified to help guide management actions. For example, if agricultural areas are common within priority areas, then actions to mitigate threats associated with agriculture (e.g., pesticides use, habitat loss) should be a focus.
The objective of our study was to identify conservation priority areas for bumble bees under current and future predicted climate scenarios. Using climate envelope models, we estimated the predicted change in distribution for bumble bees in Canada for the year 2050 under the best (Representative Concentration Pathway-RCP 2.6)and worst-case (RCP 8.5) greenhouse gas emissions scenarios according to the IPCC (Intergovernmental Panel on Climate Change) and using multiple range targets (i.e., the amount of each species' range that should be within conservation priority areas). We compared the geographic area of established protected areas within bumble bee identified priority areas to determine if already established protected areas are adequately protecting bumble bee conservation priority areas. We also measured the land cover classes within bumble bee conservation priority areas to identify land cover targets for conservation. The results of this study are an important step for the effective conservation management of critical pollinator communities under future predicted climate change.

| Bumble bee occurrences
We obtained bumble bee occurrences for North America for a recent decade of observations (2008-2018; 227,735 records, 47 species) from the Bumble Bees of North America Database (Richardson, 2020, Appendix S1). Much of the records in this database come from citizen science databases (Richardson, 2020 for a list of contributors to the database). These records are included after screening for identification and location information accuracy. Naming conventions followed the Natural History Museum of London's project "Bombus" bumble bees of the world checklist (Williams, 2022), which contains updated species lists from Williams (1998). We removed observations at the generic level, without geospatial information and duplicate records (i.e., the same latitude and longitude for the same species). Any species with fewer than 10 occurrences (i.e., Bombus cockerelli, B. distinguendus, and B. variabilis with 5, 6, and 9 occurrences, respectively) were removed because these species could not be reliably modeled in subsequent analyses. After refinements there were a total of 70,164 occurrences for 44 bumble bee species ( Figure 1a, Appendix S2).

| Climate envelope models
All analyses were performed in R version 4.0.4 (R Core Team, 2021). The species distribution models used 19 BIOCLIM  climate variables (Fick & Hijmans, 2017) at a 2.5 arcmin ($5 km) resolution for North America (https://www.worldclim.org). Bumble bee species distributions were initially modeled across their full North American range only identifying conservation priority areas in Canada (see next section) to allow species to respond to changes in climate without limiting by geopolitical borders.
Climate envelope models were created using MaxEnt (using the ENMevaluate function in ENMeval package F I G U R E 1 Methods flow chart for the climate envelope models and conservation prioritization models. (a) Bumble bee occurrence data for 44 species and climate data for current and two future climate scenarios were acquired. Future predicted climate was limited to a $ 300 km buffer around each species' distribution. (b) Ensemble models were created using MaxEnt across 7 Global Circulation Models and the three climate scenarios. (c) Ensemble models were cropped to Canada, cost data (equivalent pixel area across Canada), and three conservation targets for each species' distribution were parametrized in the conservation priority models. (d) Conservation priority areas were compared to MODIS land cover and current Canadian protected area to compare coverage across land cover classes and current protected areas. [Kass et al., 2021]) across North America for each of the 44 modeled bumble bee species under current and future RCP 2.6 and RCP 8.5 climate scenarios. MaxEnt creates climate envelop models by taking a set of climate variables to predict the suitability for species presence in other areas by comparing conditions where the species occurs to conditions in other areas using pseudoabsences or background points (Phillips et al., 2006).
We tested different combinations of feature classes and regularization parameters for each species in a fullfactorial design. Feature classes are different types of functions used to represent the relationship between the environmental variables and probability of occurrence in MaxEnt (Phillips & Dudík, 2008). Regularization parameters are used to control overfitting, with low values allowing more parameters while high values allow fewer parameters in the model (Phillips et al., 2006;Phillips & Dudík, 2008). The feature classes we tested were L, LQ, H, LQH, LQHP, LQHPT (L = linear, Q = quadratic, T = threshold, P = product, H = hinge) and the regularization parameters tested were 0.5,1,2,3,4,5.
Pseudoabsences (the same used to assess collinearity, also referred to as background points) were generated by adding 10Â the number of observations of each species and randomly sampling across North America. Pseudoabsences were not spatially restricted to allow broad sampling of suitable and non-suitable climatic conditions to increase model predictability. Future predicted distributions were restricted to 2.7 (approximately 300 km) of the historical distribution of each species. This distance was selected as it is the product of years predicted into the future (30 years) and the estimated maximum dispersal rate for bumble bees (10 km/year). Dispersal estimates for bumble bees are likely species specific however, current data on bumble bee dispersal is extremely limited. The 10 km/year maximum dispersal rate is based on the rate of expansion of an introduced bumble bee and this rate has been used in previous studies (Hingston, 2006;Sirois-Delisle & Kerr, 2018;Stout & Goulson, 2000). Thus, predicted distributions are restricted within a plausible range expansion for each species.
Bumble bee occurrences are highly concentrated near urban areas, likely due to biased sampling near dense human settlements. To adjust for this, we generated a bias file for all bumble bee occurrences from the full Bumble Bees of North America Dataset (years 18,052,019, 47 species, 585,425 observations) using a 2D-kernel density function (kde2D function in KernSmooth package [Wand, 2020]). This bias file downweighs heavily sampled areas (Fourcade et al., 2014) in the MaxEnt model.
Model evaluation was performed by withholding 20% of the occurrences and pseudoabsences to test the model (k-fold cross-validation). Model performance was compared using Area under the Curve (AUC), and the Continuous Bayes Index (CBI) (Fawcett, 2006;Hirzel et al., 2006). B. natvigi was not included in further analyses as the climate envelope models had predicted occurrences outside of the known range, low predicted values at occurrences and poor model performance (Williams et al., 2014, Appendix S3). The median AUC value for the remaining 43 bumble bee species was 0.90 (0.64 0.90), and the median CBI was 0.94 (0.120.94) Appendix S3. To estimate the percent range change associated with climate change, the current range of each species was compared to the predicted future ranges under both RCP 2.6 and RCP 8.5. This was done by comparing differences in pixel suitability (>0 was considered suitable) between current and RCP 2.6 and RCP 8.5 climate scenarios.

| Conservation prioritization
The conservation objective for the bumble bee priority areas was to conserve a targeted amount of each species' distribution while minimizing the cost of conserving the priority areas under current and two future climate scenarios. Conservation priority areas were only identified for Canada, so each climate envelop model was cropped to Canada (Figure 1c). We focused on a single country because conservation legislation is often limited by political boundaries and Canada has been identified as a country with high conservation responsibility particularly in terms of range shifts with climate change (Coristine et al., 2019). Species models were stacked by climate scenario into three separate raster stacks (one each for current, RCP 2.6 and RCP 8.5, Figure 1c).
The conservation priority areas were identified using prioritizr (Hanson et al., 2022). Prioritizr is an R package that can identify conservation priority areas using mixed integer linear programming. We selected three targets (17%, 30%, and 50%, Figure 1c) of each species' distribution determined from the climate envelope models to be in the conservation priority areas. These targets were selected as they match other international conservation goals: 17% (Aichi target, Convention on Biological Diversity [CBD], 2010), 30% (UN convention on biodiversity 2030 target, CBD, 2020a), and 50% (nature needs half, Locke, 2015). The proportion of distributions within conservation priority areas for at-risk and least-concern bumble bee species was calculated to determine if at-risk species and least concern species are adequately represented in conservation solutions. Species were classified as at-risk based on IUCN Red List classifications Critically Endangered, Endangered, Vulnerable, Near Threatened, and Data Deficient (Appendix S2).
Within priotitizr conservation priority areas were identified using the add_min_set_objective function. This function works by minimize the cost of the conservation priority areas while meeting set conservation targets. The "cost" that was minimized in our solution was area, as we assumed that a larger area would be more difficult and costly to acquire as well as manage.
The area of conservation priority for bumble bees within current protected areas and MODIS land cover class (Figure 1d) was calculated using the zonal function in the raster package (Hijmans 2020). The protected areas shapefile was obtained from Environment and Climate Change Canada (ECCC 2019). The MODIS land cover data used was MODIS Land Cover Type (MCD12Q1) version 6 (2018) with a 500 m resolution. The land cover data was re-projected to match the resolution of the species distribution models (2.5 arcmin) by the nearest neighbor method using projectRaster in the raster package.

| RESULTS
Bumble bee conservation priority areas in Canada under all climate scenarios were identified primarily in the mountainous regions of Western Canada (Pacific Maritime, Montane Cordillera, Boreal Cordillera, Boreal Plain), regions of the tundra and taiga (particularly Baffin Island and northern Quebec), northern softwood forests (Boreal Shield) of Ontario, Southern Ontario and Atlantic Canada, Mixed Wood Plain of Southern Ontario, and the Atlantic Maritime (Figure 2, Appendix S4). As the targets for conserving each species' distribution were increased from 17% to 50%, more area within the already identified priority regions was prioritized (Appendix S5) but some new areas were identified including the Prairies, as well as more area of taiga, western forests and western coastal regions (Figure 2, Appendix S4). With climate change, conservation priority areas shifted further north, into coastal regions (especially in the far north) and into areas with higher elevation (Rocky Mountains in the West and Laurentian Mountains in Quebec). Within the conservation priority areas, both at-risk and stable species were well represented (Appendix S6).
Predicted changes to species' range with climate change were species-specific. Our models of bumble bee species climatic suitability under predicted future climate change suggest about half are expected to decrease (22/43), one is expected to remain unchanged, five had mixed responses (i.e., increase in one RCP and decreased in the other) and 17 species were expected to increase in distribution (Figure 3). Most at-risk species (11) are predicted to experience a decrease in range with climate change under both RCP scenarios, with only five expecting an increase, and one each remaining stable or mixed. Generally, many of the species predicted to decline in range occur in western North America or are far northern species and those that are predicted to have an increased distribution are mainly eastern North American species or those with a wide distribution, although there are exceptions. No species were predicted to lose all climatically suitable habitat by 2050.
Identified bumble bee conservation priority areas were not well represented within already established protected areas under current or future climate scenarios (Table 1). On average, 13.7% of bumble bee conservation priority areas are within already established protected areas with a minimum of 12.3% (current climate, 30% target) and a maximum of 15.5% (future RCP 2.6, 30% target).
The most common land cover classes within bumble bee conservation priority areas are shown in Figure 4. In order of most common to least common the land cover classes are cropland, grassland, open shrubland, mixed forest, permanent snow, and ice and woody savannah. Under future climate change scenarios, grasslands, mixed  snow, and ice will occur less often in bumble bee conservation priority areas while grasslands, mixed forest, and woody savannah will occur more often which may affect where conservation efforts should focus. These results highlight regions across Canada where bumble bee conservation efforts should be focused to ensure effective and efficient conservation of North American bumble bee species. A small proportion of conservation priority areas for bumble bees were within already established protected areas. This result is consistent with previous studies examining the extent of priority areas and species at-risk within established protected areas (e.g., Bolliger et al., 2020;Deguise & Kerr, 2006). Bumble bee conservation cannot rely on current protected areas to prevent continued declines. Instead, either new protected areas need to be established within conservation priority areas or private lands within the conservation priority areas may need to be managed to support bumble bees. Although there is international support for increasing protected areas, actual establishment of new protected areas in recent years has failed to meet international targets (CBD, 2020b;Perino et al., 2022;Xu et al., 2021). Establishing protected areas is costly (Wise et al., 2012), and for these protected areas to be effective there are management costs (Le Bouille et al., 2022;Wauchope et al., 2022) and there can be costs for local communities (Green et al., 2018) highlighting the challenges of relying on protected areas for conservation. It may be prudent then to better manage private lands for bumble bee conservation. Actions on private lands that could benefit bumble bees could include improving and expanding forage, nesting and overwintering habitat including forests and early flowering plants that may be particularly beneficial for reproductive success, establishing monitoring programs, limiting exposure to managed bees which could spread disease, and eliminating the use of harmful insecticides (Gill et al., 2016;Goulson et al., 2011;Rundlöf et al., 2022;S anchez-Bayo & Wyckhuys, 2019). To encourage beneficial actions on private lands incentives could be provided to help alleviate the financial burden of creating pollinator habitat or enacting pollinator-friendly management actions (Bommarco et al., 2021;Breeze et al., 2014;Ouvrard et al., 2018).
With climate change, bumble bee priority areas were identified in higher latitudes and in more mountainous regions, a pattern which has been observed previously (Biella et al., 2017;Kerr et al., 2015;Martins et al., 2015). Regions with overlap among conservation priority areas for each examined climate scenario (current, RCP 2.6 and RCP 8.5) may be especially important for conservation because they will provide suitable climate areas now and under multiple possible emission scenarios (Rose & Burton, 2009;Vos et al., 2008). These overlapping conservation priority areas may be particularly important to concentrate conservation efforts for bumble bees as current research suggests they are not tracking shifting climates (Kerr et al., 2015;Sirois-Delisle & Kerr, 2018). If bumble bees continue failing to track climate change, these consistent conservation priority areas are predicted to be suitable throughout a range of climate scenarios and may be important future refuge areas.
Climate change will impact bumble bees in a speciesspecific manner, consistent with previous findings that bumblebee conservation management needs to be tailored instead of a "one size fits all" approach (Liczner & Colla, 2020). Most North American bumble bee species, and even more concerningly, most at-risk species, are predicted to be negatively impacted by climate change. Many eastern species or widely distributed species were predicted to have an increase in their area of habitat suitability across Canada while species that are predicted to have lower areas of climate suitability generally are western species. This may be due to predicted patterns of climate change which are expected to cause the West Coast to become warmer and drier, potentially exceeding the thermal tolerances of bumble bees (Martinet et al., 2021) and limit food availability (Aldridge et al., 2011;Phillips et al., 2018). Additionally, montane bumble bee species would have restricted regions of suitability once they can no longer move to higher elevations or if alternate suitable mountainous regions are not within range (Biella et al., 2017;Inouye, 2020;Lee et al., 2019). However, eastern species or more distributed species may simply have more area to expand as the climate warms allowing for greater predicted gains without as many geographical barriers. Bumble bee conservation and monitoring plans should consider the species-specific impacts of climate change particularly for at-risk species by incorporating conservation actions that can mitigate the impacts of climate change into species recovery documents. Actions that may mitigate the impact of climate change include identifying and protecting climate refugia, improving habitat while reducing ongoing threats, enhancing connectivity to allow dispersal to new regions (Carrasco et al., 2021;Choe et al., 2021;Heller & Zavaleta, 2009;McGuire et al., 2016).
It should be noted that in our models we assumed current climatic conditions (specifically the recent past  to be optimal for species now (nonoverlapping more recent past 2008-2018) and into the future, yet this may not be true in the future. As our climate, models are predictions as to what could occur in the future, continued monitoring for distributional changes would be recommended to validate predicted models. Similarly, our method does not include a species' ability to adapt to changing conditions, whether current conditions will continue to be favorable in new habitats, the thermal tolerances of species, or the role of biotic interactions in promoting or inhibiting species range expansion (i.e., distribution of host species for Psithyrus bumble bee social parasites [Fourcade et al., 2021;Suhonen et al., 2016]), although this is a criticism of many climate change predictions and developing methods of improving predictions is an active area of research (Littlefield et al., 2019). Also unknown are species-specific estimates of dispersal capacity and barriers to dispersal and the synergistic effects of other threats with climate change. This information is generally lacking for bumble bees (but see, Martinet et al., 2021) and future research to address these knowledge gaps would greatly improve the accuracy of bumble bee climate models and conservation effectiveness. Similar to findings in previous work with other species (e.g., Coristine et al., 2018), most bumble bee occurrences are biased to southern Canada and there are a fewer surveys the Far North and in Western Canada, and an overall lack of absence data. This limits our ability to accurately assess suitable habitat for these species as habitat suitability models often perform poorly with low sample sizes and if large portions of a species' range are not sampled (Phillips et al., 2009).
Croplands were the most common land cover identified within conservation priority areas for bumble bees. While this is good news for safeguarding pollination services provided by bumble bees, it poses considerable conservation challenge. Agricultural areas, particularly intensive agricultural areas, may lack high quality habitat, often use pesticides which can be harmful, may be a source of pathogen exposure and sites of high competition for floral resources from managed bees (Cameron & Sadd, 2020;Colla et al., 2006;de Pereira et al., 2021;Grixti et al., 2009;Siviter et al., 2023;Szabo et al., 2012). Agricultural lands within conservation priority areas should be managed for more bee-friendly practices including incorporating semi-natural habitats, reducing or eliminating pesticide use, restricting the spread of disease from managed bees and providing enough natural floral resources to mitigate competition (Crone & Williams, 2016;Hemberger et al., 2021;Klatt et al., 2020;Miller et al., 2023;Purvis, Meehan, & Lindo, 2020). Grasslands are also commonly represented within conservation priority areas although they are some of the most at risk land cover types in the world (Carbutt et al., 2017;Hoekstra et al., 2005). Current efforts to restore and conserve grasslands should also benefit bumble bee conservation by providing high-quality habitat (Purvis, Vickruck, et al., 2020;Redpath-Downing et al., 2013;Rosenberger & Conforti, 2020). Mixed forests are expected to be more represented in conservation priority areas in the future with climate change. Forests are important sources of early-forage plant species which may be critical for the reproductive success of colonies, and forests may be important nesting and overwintering habitat (Malfi et al., 2019Mola, Hemberger, et al., 2021;Pugesek & Crone, 2021). Although often overlooked in bumble bee conservation, the role of forests in providing high quality habitat particularly for queens should be considered in bumble bee conservation plans. We did not model land cover change with future climate change, so we cannot assess whether landscape conditions will be suitable in future modeled priority areas. However, we do not expect drastic land cover change across Canada by 2050, althought this may be optimistic given increasing wildfire risk and pressure for development. Indeed, at least three land cover types (forest, cropland, and rangelands) have remained stable since 1960 (Winkler et al., 2021). In addition, land cover change can be difficult to predict due to sometimes rapidly changing local conditions such as development projects and urbanization. Future research studies would benefit from modeling predicted land cover change and how that would impact conservation priority areas (Kremen et al., 2007) as interactions between climate and land cover change can negatively impact biodiversity (Betts et al., 2019;Mantyka-Pringle et al., 2012).
Here, we identify priority areas across Canada needed for the preservation of bumble bees under climate change. These areas are not consistent under current and future climate scenarios, suggesting protected areas need to be expanded to capture the wide range of conditions species can occupy. Supporting bumble bee species under climate change will become increasingly important in the next few decades, especially since these pollinators are already at-risk of extinction due to other anthropogenic threats such as pathogen spillover, pesticide use, and habitat loss. Identification of priority areas, as we have done here, is a crucial step toward the conservation and management of at-risk species, native biodiversity and ecosystem services for long-term sustainability.

ACKNOWLEDGMENTS
We would like to thank Dr. Jeffrey Hanson for assisting with the prioritizr function. We would also like to thank Genevieve Rowe, Sarah MacKell, Hadil Elsayed and Dr. Victoria MacPhail for providing input on methodology. We would like to thank NSERC for providing funding for ARL's salary.