Bumble bees are important pollinators in temperate forested regions where fire is a driving force for habitat change, and thus understanding how these insects respond to fire is critical. Previous work has shown bees are often positively affected by the post-fire environment, with burned sites supporting greater bee abundance and diversity, and increased floral resources. The extent to which fire impacts variation in bumble bee site occupancy is not well understood, especially in higher latitude regions with dense, primarily coniferous forests. Occupancy models are powerful tools for biodiversity analyses, as they separately estimate occupancy probability (likelihood that a species is present at a particular location) and detection probability (likelihood of observing a species when it is present). Using these models, we tested whether bumble bee site occupancy is higher in burned locations as a result of the increase in canopy openness, floral species richness, and floral abundance. We quantified the impact of fire, and associated habitat changes, on bumble bee species' occupancy in an area with high wildfire frequency in British Columbia, Canada. The burn status of a site was the only significant predictor for determining bumble bee occurrence (with burned sites having higher occupancy); floral resource availability and canopy openness only impacted detection probability (roughly, sample bias). These findings highlight the importance of controlling for the influence of habitat on species detection in pollinator studies and suggest that fire in this system changes the habitat for bumble bees in positive ways that extend beyond our measurements of differences in floral resources and canopy cover.

Our study was conducted on the unceded territories of the Nuxalk and Ulkatcho First Nations, in and around Tweedsmuir Provincial Park in British Columbia, Canada, from June to August 2019 (Fig. 1). We established sites both in and adjacent to four wildfire zones, two of which are recent burns (2017 and 2018) and two of which are older burns (2009 and 2010), though the older burns had not significantly regenerated, as high burn severity and high elevation have limited tree and shrub regrowth. Three of the burns were large in scale (2,800 to over 7,000 hectares), and the most recent burn was smaller (40 hectares). The unburned sites were in forest habitat adjacent to each of these burned areas. 

We sampled a total of 26 circular sites of 100m diameter, 13 in areas impacted by fire (4 in each large fire and 1 in the smaller fire) which we call "burned" sites, and 13 in nearby unaffected forest, which we call "unburned" sites. We selected sites such that edges were a minimum of 1 km from all other site edges to ensure spatial independence, as bumble bee foraging most frequently occurs within 1 kilometre of their nest (Greenleaf 2007, Geib 2015, Kendall 2022). We visited each site twice over the course of the season and, due to logistical constraints, sampled groups of 3–8 sites in spatial and temporal blocks with block composition differing slightly between visits. However, sites were re-sampled in a similar order, such that visits to sites were separated by similar time periods (4–6 weeks between revisits). 

To sample bumble bees, we used blue vane traps (three per site), collecting samples after traps had been out for 2–4 days, in order to ensure minimal negative impacts on bee populations (see Kimoto 2012 and Gibbs 2017 for evidence of negative impacts of long-term trap collecting). We selected blue vane traps as our primary collection method because our site arrangement and sampling structure necessitated the use of passive sampling and because previous work has shown that blue vanes are one of the most effective for per-sample species accumulation (Joshi 2015). Blue vanes are highly attractive to bumble bees (Stephen 2007) and have been shown to collect similar species sets to those obtained by active netting (Rao 2009). In addition, we performed supplementary spot netting surveys, only at burned sites, for 60 person-minutes per visit either during trap setup or trap take down. We did not net bees at unburned sites because, early in the season, many unburned sites had few to no open flowers from which we could collect bees. Netting was conducted as long as temperature was greater than 15 degrees C, wind speed was below 2 m/s, and there was no precipitation. We identified each bumble bee to species using the key in Williams' North American field guide (2014), and follow the bumble bee taxonomy therein, with the updated revision to Alpinobombus for B. kirbiellus (Williams 2019). 

During sampling visits, we also recorded site-level habitat variables, either when blue vane traps were set up or when they were collected. To quantify canopy openness and floral resource availability, we established two 100m transects in N-S and E-W directions at each site. For canopy openness, we took six evenly-spaced upward-facing photographs per transect (using a Canon 5D MK I with Sigma 8mm f/3.5 EX DG Circular Fisheye Lens), for a total of 12 photos per site. We counted all open flowers from a total of 76 different species (we include a full list of floral species and information about attractiveness to bumble bees as a supplementary file, Online Resource 1) along each of the 3mx100m transects, identified to species or genus using local field guides and online resources (Pojar 1994, Parish 1999, e-flora-BC). We calculated floral abundance and species richness by pooling open flower counts (later, floral abundance was log-transformed) and number of flowering plant species across transects. We measured canopy openness at a site level (once per site) and floral resource information at a visit level (twice per site).

To determine canopy openness at each site, we analysed the upward fisheye photos in Gap Light Analyser (GLA), a program designed for analysis of hemispherical canopy cover photos (Frazer 1999). We used default settings (Registration: Geographic North, Location: none added, Orientation: horizontal,  Topographic shading: Use topographic mask data, Solar time step: 2 minutes, Azimuth regions: 36, Zenith regions: 9, Data source: modelled, Solar constant: 1367 Wm^-2, Cloudiness index: 0.5kt, Spectral fraction: 0.5, Units: Mols m^-2 d^-1, Beam fraction: 0.5, Sky-region brightness: UOC Model) along with a custom projection distortion specific to our lens. GLA relies on contrast between sky and foliage to determine percent canopy cover. This required that we sometimes draw boundaries manually and then set local thresholds accordingly in order to ensure correct classification. We used a blue colour plane, as recommended, to enhance contrast between canopy cover and sky. In some cases (e.g., when the sun reflected off trees) it yielded a "canopy" section that was brighter than sky. To ensure correct classification, we manually traced the canopy cover and applied the colour fill tool. To calculate canopy cover at the site level, we calculated the mean cover across the 12 photos for each site. 

All data are in comma-separated variable (.csv) format which can be read by any spreadsheet program (Microsoft Excel, Google Sheets, Apple Numbers) or any text editor. All code is in .R format, which can be read by any text editor and interpreted by the R program.