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.
Home | Current Knowledge | Threats | Conservation | Next Steps