#### State of Vermont’s Wild Bees - Methods

Modeling species distributions

Bee occurrence records
Presence data were obtained from Global Biodiversity Information Facility (GBIF.org) and consisted of historical museum records, and community science observations submitted to iNaturalist.org. Occurrence data were obtained on 13 April 2022 and available at https://doi.org/10.15468/dl.j53x2g.

Distribution models
We used the presence-only data to generate species distribution models using Maxent (Phillips et al. 2006). We included several biologically relevant bioclimatic variables in our species distribution models. The bioclimatic variables (resolution: 1km x 1km) were obtained from AdaptWest.org and were based on the Coupled Model Intercomparison Project phase 6 (CMIP6) climate models. In addition to the bioclimatic variables, we included a binary land cover surface indicating terrestrial and aquatic environments. We accounted for sampling bias by deriving background locations, or ‘pseudo-absences’ from a probability distribution that accounted for human population density (2020), human modification (2016, road density, powerlines, etc) and all observations in the Class Insecta. Each component of the sampling bias layer was scaled between 0 and 1 so they received equal weighting in the sampling bias layer. We fit several species distribution models using the ENMevaluate package (Muscarella et al. 2014). We allowed the smoothing parameter to vary between 0.5 and 2.5 while fitting linear and quadratic functions within the distribution model. We used ‘block’ partitioning to aggregate data into four bins based on latitude and longitude. We chose ‘block’ partitioning so our species distribution models could be extrapolated to future climate scenarios. ‘Block’ partitioning facilitates model transfer to potentially non-analogous environmental conditions (Wenger and Olden 2012, Muscarella et al. 2014). We fit distribution models for species with 5 or more independent observations from unique locations within the broader New England region. Model parameters for each species are available upon request. The input models for all species were identical.

The model with the lowest AICc value was used to determine a species distribution. We converted the probability surface created from the most parsimonious model into a binary surface that represented the locations where the species is likely to be present across the region. We used a sensitivity threshold to objectively assign whether a species was likely to be present or not at a given location. Locations across the northeast that were greater than the threshold value were deemed locations where a species was likely to be present while locations with predicted values less than the threshold value are locations where the species is unlikely to be present. We used the model with the lowest AICc score to make predictions of future distributions given climate change. We used several climate projection ensemble models to model how species distributions may change given different climate conditions (RCP 2.6, 4.5, 7.0, 8.5). RCP 8.5 is associated with rising carbon emissions resulting in an average global rise in temperature of 4.3ºC by the year 2100, while RCP 4.5 is associated with slowly declining carbon emissions and an average global rise in temperature of 2.4ºC (Zeller et al. 2021).

Locating unique communities

We used the resulting species distribution maps to identify areas within the region that support unique insect communities (Dansereau et al. 2022). To do this, we created a binary species occurrence matrix (1 = present, 0 = absent) for each 1x1km grid cell across the region for each insect order. We transformed species presence data using a Hellinger transformation using the labdsv R package (Roberts 2019) for use in determining beta diversity. We measured total beta diversity as the variance of the community matrix. We determined each cell’s local contribution to beta diversity (LCBD). The LCBD value is a quantitative measure for how much a given cell (1x1 km) contributes to the overall variance in the community (Legendre and De Cáceres 2013).

We identified species that contribute disproportionately to beta diversity by summarizing the LCBD score for each species within an order. However, species with larger ranges would have higher values than species with lower values. Therefore, we modeled the expected relationship between distribution size and species LCBD scores. We assessed the role of each species’ contribution to beta diversity using the residuals from the aforementioned regression. Residual scores that fell more than 2 standard deviations away from zero were considered to be species that play a disproportionate role in beta diversity.

Climate sensitivity

We used the relative variable importance from the most parsimonious species distribution model (lowest dAICc value) as an indicator for species sensitivity to climate change. We categorized the bio climatic and physical explanatory variables into drivers of change, such as temperature, precipitation and physical characteristics. We summarized the relative importance scores that represent drivers or season of change.

Rate and direction of distribution change

We quantified species richness across the region by summarizing the current species distribution models and future predictions for each order. We calculated an analog of climate velocity (km/yr) using the species richness information for each insect order observed in northeastern North America with 3 or more species distributions. We used the dVoCC function in the VoCC R package (García Molinos et al. 2019) to quantify the velocity and direction communities would need to move to maintain current species richness. We summarized the resulting velocity and directional data across the region to make inference about velocity and potential movement directions. Velocity metrics are reported in km/year measured using great circle distance while direction is reported in degrees where 0 indicates northward or pole ward movement while 180 indicates a southward or equatorial movement. Areas with the slowest rate of change were identified as areas of potential refugia. In addition to the metrics described above, we calculated the anticipated rate of change in species richness across the region for the four carbon emissions scenarios for each 1x1km grid cell using linear regression.

State Rankings

S-ranks were initially generated using The NatureServe Conservation Rank Calculator and data available through GBIF. The results were then reviewed by Spencer Hardy and several regional bee experts (John Ascher, Tracy Zarrillo, Michael Veit), leading to refinement of many rankings.

The RSGCN list was developed by the Northeast Fish and Wildlife Diversity Technical Committee in consultation with regional experts, including Spencer Hardy.

The 2015 VT State Wildlife Action Plan (SWAP) designated SGCN species, including several Bumble Bees. These lists also relied on state experts, including Kent Mcfarland who assisted with the ranking of the Bumble Bees.

Watchlist

Using the expert opinion, S-ranks, G-ranks, RSGCN, and SGCN lists, we created a watchlist of species of potential conservation concern. By default, all SGCN and RSGCN species that occur within the state were included, except for Bombus perplexus (SGCN high priority) which we felt should be removed from the SGCN list in 2025. Other criteria for inclusion included; specialized habitat, very few records in the region, and/or some indication of population decline. Not all S1 species were included, since several are at the northern edge of their range in Vermont and are abundant elsewhere.

Threat matrix

We ranked the possible effects of various threats for most genera to create the threat matrix. Rankings were based on various natural history traits and published literature. Nesting location, diet, flight window, habitat specificity, and overwinter stage were considered in most rankings, with parasitic genera usually ranked similarly to their hosts.

For example, bees at the northern edge of their range (i.e. Xylocopa, Melissodes bimaculatus) are likely to benefit from warmer temperatures, while those near their southern distribution limits in Vermont (i.e. Bombus borealis, Melissodes illiatus) will likely be negatively affected. Bees found primarily in semi-open settings (i.e. residential areas) will likely benefit from forest fragmentation, while forest obligate species are likely to be negatively impacted. If multiple species in the same genera had different anticipated responses, we scored the genera ‘mixed impacts’.

Genera represented by a few species or only be introduced species were omitted for clarity unless we suspected strong responses to multiple threats (i.e. Xylocopa).

These rankings are not definitive but represent the possible impact on populations.

Flower Interactions

When known, flower associations for individual bees were included in all datasets published for this project. An additional 7,000+ records were collected through iNaturalist, using the observation field “interaction -> visited flower of:”. This data relies on community scientists for collection and verification, though many records were added and/or checked by project staff. Associated flower data was cleaned prior to use.

Bees in Agriculture

Flower interaction data was filtered to generate a list of bees that visited agricultural crops. Plants were generally assigned as agricultural if they are commercially available as a food product. This definition included wild plants such as Acer saccharum and Allium tricoccum, but did not include plants such as Apios americana and Amelanchier that are occasionally eaten, but rarely sold. Some leniency was given around genera that contained both commercial and wild species, particularly Rubus, Vaccinium, and Fragaria. In many of these cases, the plant was only recorded at the genus level, in which case it was included, as were most wild members of the genera, assuming the fruit and flower were similar to the cultivated species. A few wild Vaccinium (i.e. Vaccinium stamineum) and Rubus (i.e. Rubus odoratus) were excluded given significant differences in flower morphology.

Registered Honey Bee hive locations

We received the location of all registered honey bee hives through a freedom of information request from the Vermont Agency of Agriculture, Food & Markets. Vermont has a mandatory registration program that collects the number and location of all active honey bee hives in the state. Not all hives were reported with coordinates, when they weren’t available we geo-referenced to the extent possible from the address or town provided. Data from 2020 through 2022 were included.

#### Literature cited

Dansereau, G., P. Legendre, and T. Poisot. 2022. Evaluating ecological uniqueness over broad spatial extents using species distribution modelling. Oikos 2022:e09063.
García Molinos, J., D. S. Schoeman, C. J. Brown, and M. T. Burrows. 2019. VoCC: An r package for calculating the velocity of climate change and related climatic metrics. Methods in Ecology and Evolution 10:2195–2202.
Legendre, P., and M. De Cáceres. 2013. Beta diversity as the variance of community data: Dissimilarity coefficients and partitioning. Ecology Letters 16:951–963.
Muscarella, R., P. J. Galante, M. Soley-Guardia, R. A. Boria, J. M. Kass, M. Uriarte, and R. P. Anderson. 2014. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution 5:1198–1205.
Phillips, S. J., R. P. Anderson, and R. E. Schapire. 2006. Maximum entropy modeling of species geographic distributions. Ecological Modelling 190:231–259.
Roberts, D. W. 2019, August. Labdsv: Ordination and Multivariate Analysis for Ecology.
Wenger, S. J., and J. D. Olden. 2012. Assessing transferability of ecological models: An underappreciated aspect of statistical validation. Methods in Ecology and Evolution 3:260–267.
Zeller, K. A., C. A. Schroeder, H. Y. Wan, G. Collins, K. Denryter, A. F. Jakes, and S. A. Cushman. 2021. Forecasting habitat and connectivity for pronghorn across the Great Basin ecoregion. Diversity and Distributions 27:2315–2329.