Wild Bee Diversity Increases With Local Fire Severity in a Fire‐Prone Landscape

As wildfire activity increases in many regions of the world, it is imperative that we understand how key components of fire-prone ecosystems respond to spatial variation in fire characteristics. Pollinators provide a foundation for ecological communities by assisting in the reproduction of native plants, yet our understanding of how pollinators such as wild bees respond to variation in fire severity is limited, particularly for forest ecosystems. Here, we took advantage of a natural experiment created by a largescale, mixed-severity wildfire to provide the first assessment of how wild bee communities are shaped by fire severity in mixed-conifer forest. We sampled bees in the Douglas Fire Complex, a 19,000-ha fire in southern Oregon, USA, to evaluate how bee communities responded to local-scale fire severity. We found that fire severity served a strong driver of bee diversity: 20 times more individuals and 11 times more species were captured in areas that experienced high fire severity relative to areas with the lowest fire severity. In addition, we found pronounced seasonality in the local bee community, with more individuals and more species captured during late summer, especially in severely burned regions of the landscape. Two critical habitat components for maintaining bee populations—flowering plants and boring insect exit holes used by cavity-nesting bees—also increased with fire severity. Although we detected shifts in the relative abundance of several bee and plant genera along the fire severity gradient, the two most abundant bee genera (Bombus and Halictus) responded positively to high fire severity despite differences in their typical foraging ranges. Our study demonstrates that within a large wildfire mosaic, severely burned forest contained the most diverse wild bee communities. This finding has particularly important implications for biodiversity in fire-prone areas given the expected expansion of wildfires in the coming decades.


INTRODUCTION
Animal pollinators play a critical role in maintaining biodiversity in natural and managed systems, with most of the world's flowering plants either requiring or benefitting from pollination services for reproduction (Ollerton et al. 2011). Changes in animal pollinator populations can drive changes in wild plant persistence (Biesmeijer et al. 2006, Brosi andBriggs 2013), with implications for flowering plant diversity as well as the broader ecological communities that depend on flowering plants as critical habitat features. In addition to promoting and maintaining biodiversity, an estimated 70% of global crops used for human consumption are enhanced by pollination services (Klein et al. 2007), valued at €153 billion annually (Gallai et al. 2009).
Concerns remain high regarding the economic and ecological consequences of long-term population declines of pollinators (Pollinator Health Task Force 2015, IPBES 2016), yet we lack a foundational understanding of wild pollinator population dynamics in many natural and human-dominated ecosystems (Bartomeus et al. 2018).
Understanding how ecological communities respond to large-scale disturbances is a critical topic for successful pollinator conservation, especially as many disturbance regimes shift worldwide (Turner 2010). Wildfire is a natural disturbance agent that plays a key role in shaping many terrestrial ecosystems (DeBano et al. 1998, Pyne et al. 1996, Bowman et al. 2009), and it can enhance biodiversity via increases in habitat heterogeneity and reductions in species dominance (Bond and van Wilgen 1996, Pyne et al. 1996, DeBano et al. 1998. Historical wildfire regimes have shifted markedly due to extensive suppression efforts of the past (Bowman et al. 2009), and global climate change has led to alterations in the frequency, severity, and extent of wildfire (Westerling et al. 2006, Miller et al. 2009, Flannigan et al. 2013, Jolly et al. 2015. In the western United States, the average wildfire season length increased 64% between the periods of 1970-1986and 1986(Westerling et al. 2006, and total burn area has increased at a rate of 355 km 2 per year (Dennison et al. 2014). Fire activity is projected to continue to change rapidly in the coming decades, potentially faster than many terrestrial species' ability to accommodate (Krawchuk et al. 2009). Given these ongoing changes, a stronger understanding of how wildfire characteristics influence biodiversity is critical, and this is particularly true for pollinators inhabiting fire-prone landscapes.
Despite the important and changing role of wildfire in most biomes, we have a poor understanding of how characteristics of fire influence pollinator diversity within forest ecosystems (Koltz et al. 2018, Rivers et al. 2018a. Fire severity is an important characteristic of wildfire that is measured by the amount of change in organic matter due to burning (Keeley 2009), and wildfires can have markedly different effects on ecological communities depending on the degree of fire severity. For example, high-severity fire typically acts as a stand-replacing event that kills dominant overstory trees, removes the forest canopy, and exposes mineral soil; in turn, this is likely to shift communities back to the early successional pre-forest stage. In contrast, lowseverity fires only burn the most flammable fuels and often have little influence on tree mortality (Perry et al. 2011). In fire-adapted forests, many understory plant species have adaptive strategies that allow them to compete after canopy and/or duff cover is reduced by severe fire (Bond and van Wilgen 1996). Although fire severity influences the composition of understory plant communities by selectively influencing recruitment (Bond and van Wilgen 1996), severe wildfires are followed by a temporary increase in plant species richness and forb cover in fire-adapted regions due to reduced light competition (Keeley 1987, Keeley et al. 2003, Huisinga et al. 2005. The differences in post-wildfire habitats caused by variation in fire severity, in turn, are expected to have important consequences for local pollinator communities. Nevertheless, nearly all studies examining post-fire responses of pollinators in forest ecosystems to date have either used limited-severity prescribed fire (Campbell et al. 2007, Rubene et al. 2015, Rodr ıguez and Kouki 2017 or have focused on small fires (Potts et al. 2003, Bogusch et al. 2014, Lazarina et al. 2017, both of which are unrepresentative of the mosaic of fire severity that is typical of contemporary wildfires in many regions of the world (Strauss et al. 1989, Stocks et al. 2002, Lentile et al. 2005, Halofsky et al. 2011. The lone exception is a recent study that investigated the influence of the diversity of fire histories (pyrodiversity) on pollinators after a mixed-severity fire in forest/scrub habitat (Ponisio et al. 2016), and researchers found positive effects of pyrodiversity on pollinators only within plots that burned at low to moderate fire severity (Ponisio et al. 2016). Given that highseverity wildfires are expected to have stronger effects on pollinator habitat (i.e., increased flowering plant abundance) than low-severity wildfires as outlined above, additional research is needed to understand the response of pollinator communities and their habitat to local fire severity itself. Furthermore, because pollinators play a fundamental role in maintaining plant populations that serve as the basis for terrestrial food webs (Kearns and Inouye 1997), understanding pollinator response to wildfire will help us ❖ www.esajournals.org 2 April 2019 ❖ Volume 10(4) ❖ Article e02668 predict broader ecological shifts as wildfire activity continues to increase in many regions of the world (Jolly et al. 2015) In this study, we took advantage of a natural experiment provided by a large-scale wildfire to evaluate how a key pollinator group-wild, native bees-responded to a mosaic of fire severity within a fire-prone landscape. First, we hypothesized that bee abundance and richness would respond positively to fire severity because more severe wildfires remove canopy to create open areas, promoting wildflowers that provide pollen and nectar needed to feed bee progeny (Bond and van Wilgen 1996, Huisinga et al. 2005, Van Nuland et al. 2013, Bogusch et al. 2014, Burkle et al. 2015, Bassett et al. 2017. Second, because fires can extend the temporal availability of floral resources (Mola and Williams 2018), we hypothesized that the relationship between fire severity and the density of flowering plants would be greatest at the beginning and end of the bee flight season, when bloom is most scarce. Third, fire can promote nesting substrates for ground-nesting species (i.e., exposed soil) while consuming all but the largest available nesting substrates for cavity-nesting bees (dead wood, Potts et al. 2005, Moretti et al. 2009), so we predicted that fire severity would be positively related to the proportion of exposed bare ground and negatively related to the abundance of dead wood (i.e., snags, stumps, and coarse woody debris) in forest patches.

Study area
We focused the study in the Douglas Fire Complex in the Klamath Mountains of the Klamath-Siskiyou ecoregion of southwestern Oregon, USA (Fig. 1), during the spring and summer of 2016 and 2017. The Klamath-Siskiyou ecoregion has a Mediterranean climate with hot, dry summers and wet winters. It is a steeply mountainous area that is characterized by a mixed-severity fire regime (Taylor and Skinner 2003). Prior to European settlement, the fire return interval in the ecoregion was a highly variable function of vegetation type, topography, and elevation (Odion et al. 2004). The median fire return interval in the Douglas-fir-dominated portions of the Klamath mountains was an estimated 14.5 yr pre-settlement, lengthening to 21.8 yr during a century of active fire suppression (1905-1992Taylor and Skinner 1998).
We conducted our study within the Dad's Creek Fire and Rabbit Mountain Fire, which composed nearly the entire Douglas Fire Complex. Lightning ignited these fires in late July of 2013 which burned for >1 month and created a mosaic of fire severity across the mixed-conifer landscape. The Douglas Complex fires burned 9400 ha in private forests and~10,200 ha in public forests managed by the USDI Bureau of Land Management (BLM; Zald and Dunn 2018). BLM holdings within this region, from which we selected our study sites, were typically managed as even-aged Douglas-fir stands before the fire, with some snags and live trees left after harvest for wildlife habitat.
We determined fire severity within the Douglas Complex using the Relative differenced Normalized Burn Ratio (RdNBR; Miller and Thode 2007), a satellite imagery-based metric of pre-to post-fire change that correlates with basal area mortality (Reilly et al. 2017). We selected RdNBR, calculated using data from the 2014 Monitoring Trends in Burn severity database, for ❖ www.esajournals.org 3 April 2019 ❖ Volume 10(4) ❖ Article e02668 several key reasons. First, RdNBR accounts for pre-fire spectral differences associated with variation in vegetation cover (Miller et al. 2009), so the metric allowed us to separate pre-existing site characteristics (i.e., stand age) with the degree of change caused by the wildfire. Second, RdNBR is one of the fire severity metrics commonly used by land managers and researchers to quantify fire severity, as it is based on widely available Landsat data (Miller and Thode 2007 We selected 35 study sites that spanned the fire severity gradient within the Douglas Complex based on published relationships between RdNBR and basal area mortality (Reilly et al. 2017; Appendix S1: Table S1). To select these sites, we randomly generated points on an RdNBR map layer using ArcGIS. We only used sites where >6.5 ha of forest burned within the same severity category, a balance between needed replication and wild bees' typical foraging ranges (Zurbuchen et al. 2010). In addition, we required that sites within the same severity ❖ www.esajournals.org 4 April 2019 ❖ Volume 10(4) ❖ Article e02668 category were >1 km apart to minimize spatial autocorrelation of bee communities, as most wild bee species in conifer forests are small-to moderate-sized bees (Rivers et al. 2018b) whose foraging ranges are typically <1 km (Zurbuchen et al. 2010). We also required that sites were >50 m from the nearest road (mean distance = 158 m) to minimize influences of flowering weeds in ditches. The 35 selected study sites represented nearly the full range of fire severity within the Douglas Fire Complex (RdNBR range = 50-1037). For reference, RdNBR values of <235, 235-649, and >649 correlate with basal area mortality of <25, 25-75, and >75%, respectively (Reilly et al. 2017).

Bee and habitat sampling
We collected bees during four sampling rounds in 2016 and 2017 during the growing season (May-September). In both years, sampling periods were spaced by 3-4 weeks to represent the flight seasons of the regional bee community. Our earliest sampling (mid-May) coincided with bloom of early-season plants attractive to native bees (e.g., Vaccinium ovatum, Gaultheria shallon, Berberis nervosa, Arctostaphylos spp.), and the final sampling coincided with end of bloom for most forbs within the study region. Very little rainfall occurs between mid-May and September in the ecoregion (National Oceanic and Atmospheric Association 2017), so selecting these sampling dates reduced the likelihood that weather would influence our sampling results. Active wildfires prevented us from completing the final sampling round in 2017, so only 22 of the sites could be sampled in that collection period. However, because we sampled sites in a random order, the sites covered in this final sampling period were not biased by fire severity (RdNBR range covered: 60-1037). We refer collection to periods as late spring (late May), early summer (late June/ early July), mid-summer (early August), and late summer (early September).
During each sampling period, we passively sampled bee communities using blue vane traps with yellow collection bottles (Springstar, Woodinville, Washington, USA) whose UV-reflective vanes are known to attract diverse bee taxa in temperate ecosystems (Stephen and Rao 2005) including conifer forest (Rivers et al. 2018b). For each sampling site, we hung a trap each on two separate 1.8 m tall posts so that each trap was at or above the typical height of flowering forbs and shrubs, and thus visible to foraging bees. Traps had no killing agent or preservative. To avoid placement bias, the first trap was placed 10 m from the site center along a randomly selected azimuth, and then, the second trap was placed 10 m from the center in the opposite direction. We avoided placing traps in dense vegetation patches to standardize visibility for bees as much as possible. We left traps for 48 h, after which we sealed them, put them in a cooler with dry ice to kill captured insects, and then froze the contents for later curation.
We assigned each bee to genus using keys from Michener (2007) and Stephen et al. (1969) and to species/morphospecies using both regional synoptic collections and local keys for species such as Agapostemon (Stephen et al. 1969), Anthophora and Ceratina (Discoverlife.org), Bombus (Williams et al. 2014), Halictus (Roberts 1973), and Xylocopa (Hurd and Moure 1963). No species keys were available for several genera in our region (e.g., Lasioglossum (Dialictus) and some Osmia), so we could only group individuals as morphospecies for these genera. Voucher specimens from our study are to be housed in the Oregon State University Arthropod Collection in Corvallis, Oregon, USA (http://osac.oregonsta te.edu/).
To quantify the influence of local habitat variables on bees collected in traps, we established transects that extended 50 m from each trap where we quantified available floral resources, canopy cover, and extent of bare ground in each site within one week of each sampling period. We assessed flowering plant density using ordered distance sampling, which is efficient for estimating flowering plant density in areas with patchy or sparse bloom (Nielson et al. 2004). At every 10-m interval along the two transects per site, we identified and measured the distance to the 4th nearest flowering plant to the nearest cm (Nielson et al. 2004). We then averaged plant densities by collection period and site for the statistical model. In addition, we recorded all flowering plant species seen while collecting the density data as a measure of species richness. Like the BVT samples, we were unable to complete all the floral resource transects during the final round of 2017. As with the bee sampling, ❖ www.esajournals.org 5 April 2019 ❖ Volume 10(4) ❖ Article e02668 active wildfires prevented us from completing the final sampling round in 2017, and we were only able to sample 17 sites floral resources during this final sampling round (RdNBR range covered: 60-1037). We used the same floral resource transects to systematically measure canopy cover and bare ground in each site. We measured canopy cover to the nearest percent at each of the 10 points along the floral resource transect using a spherical densitometer and then averaged the readings by collection period and study site. At each point, we also visually estimated the extent of bare ground at a 1 m radius to estimate the amount of open ground for ground-nesting bee habitat. We considered bare ground to be any area that lacked vegetation and did not include rocks large enough to obscure access to the ground, woody debris, or other materials impenetrable for ground-nesting bees. We recorded bare ground estimates to the nearest 5% and then averaged the data by site.
To estimate potential nesting substrates for cavity-and wood-nesting bees, we quantified the amount of dead wood within two 12.5-m radius plots per site by measuring snags, stumps, and other woody debris. Within each plot, we counted and measured the width (in categories of 5-10 cm, 10-30 cm, and >30 cm diameter) and length (to nearest cm) of all coarse woody debris (CWD) along a line transect through the plot center. We measured the height (m) and diameter (cm) at breast height of all snags and stumps within the plots. From these data, we calculated the average volume of dead wood per plot that was potentially usable by cavity-nesting bees. In addition, we estimated the number of dead wood pieces that had potential nest sites by recording the number of beetle holes in the 1 m section of CWD crossing the transect and in the 1 m section of snags and stumps at breast height. We counted holes that were >1 mm, as cavitynesting bees use a wide diversity of nesting substrates ) and the preferred cavity size is still unstudied for many species.
Finally, we measured abiotic site-level characteristics that could indicate bee or plant habitat conditions. We measured elevation (MASL) at the site center using a handheld Garmin GPS unit. We also measured aspect to the nearest degree at each trap location using a compass and then calculated the average of each per site. Finally, we estimated pre-fire stand age using 2012 Rogue Valley light detection and ranging (LiDAR) data, provided by the Medford District BLM. Though these data are likely to over-estimate the age of older stands with multiple cohorts of trees, they provided us with the best available estimate of pre-fire stand age. Estimates are averaged across stands and binned to 10-yr average age categories.

Analyses
We modeled the effect of fire severity on mean trapped bee abundance, bee species richness, and flowering plant density using generalized linear mixed models (GLIMMIX; SAS 9.4, SAS Institute, Cary, North Carolina, USA). In each model, we included the year (2 levels), collection period (4 levels), and fire severity (site RdNBR, continuous) as fixed effects, as well as fire severity 9 year and fire severity 9 collection period interactions. Including the "year" and "collection" interaction accounts for annual and phenological environmental changes that may influence the magnitude of the effect of RdNBR on bee populations. We also included three additional site-level covariates (i.e., elevation, aspect, and stand age) to control for potential variation between sites in characteristics that may influence bee and/or plant communities. We included two random effects: study site and study site 9 year. Before running our models, we first checked for collinearity of variables (PROC CORR; SAS 9.4) to ensure covariates were independent (R 2 ≤ 0.5). Because traps did not capture bees in some sites during some sampling periods (16% of the 554 traps set), we used a negative binomial distribution for bee abundance and richness models to avoid overdispersion, requiring a log link. Flowering plant density, which we estimated as the log-transformed number of flowering plants per hectare, was modeled with a normal distribution using an identity link.
We present predicted values from the generalized linear mixed models for bee abundance, bee species richness, and flowering plant density as marginal means with elevation, aspect, and stand age held at their means. To explore changes in bee habitat with fire severity, we also plotted non-seasonal habitat variables (canopy cover, bare ground availability, dead wood ❖ www.esajournals.org 6 April 2019 ❖ Volume 10(4) ❖ Article e02668 volume, nesting cavity availability) as scatterplots with RdNBR and then measured the correlation between each of these habitat variables and RdNBR using SAS 9.4 PROC CORR.
In addition to considering the entire bee community, we constructed the same models for the two most abundant bee genera collected to determine whether their numbers responded the same or differently to fire severity. This model included a fixed variable ("group") that represented the two focal genera and a group 9 RdNBR interaction, with the interaction demonstrating the extent of differences in the magnitude of the relationship between abundance and fire severity between the two genera. For this simplified model, we summed all bees from the selected genera by site and did not look at changes over years or collection periods. Due to the reduced degrees of freedom, we only included two additional variables in the model: stand age (fixed effect) and study site (random effect). We illustrate results as a bubble plot showing total observed abundance of species within the two genera along the fire severity gradient.
To illustrate changes in community composition with fire severity, we calculated bee and flowering plant dissimilarity by species and genus across sites using the Bray-Curtis index and visualized clustering of sampling units using nonmetric multidimensional scaling (NMDS; Minchin 1987). NMDS is a method for graphically demonstrating the dissimilarity between variables in a reduced number of dimensions based on the pairwise distances between sites given the metrics of interest. It then positions the sites graphically in an assigned number of dimensions to maximize the rank correlation between the pairwise inter-site distances of the population metric and those of the graphical ordination (Quinn and Keough 2002). To construct our ordination plots, we first calculated the total number of each genus and species observed at each RdNBR value over all collection periods and years, and then implemented NMDS for genus and species abundances separately using the metaMDS function in the R vegan package (version 2.4-5; Oksanen et al. 2013). This function finds a stable solution to the ordination using several random starts. It then adds species scores to the site ordination, using a double square root transformation to down-weight the importance of the most abundant taxa (Legendre and Gallagher 2001). We used stressplots to determine whether the fit was sufficient for plotting (<0.2). We then plotted RdNBR as a continuous environmental vector onto the ordination using function envfit (Oksanen et al. 2013) and selected the 25% of bee taxa and 25% of the floral taxa that were most highly correlated with the environmental vector to include on the plots using species selection for ordination plots (ordiselect) in the goeveg package of R.

RESULTS
We collected a total of 3220 bees representing 26 genera and 105 species/morphospecies (Table 1). More bees were collected in 2017 (n = 1753) compared to 2016 (n = 1467) despite a reduction in trapping effort caused by active wildfires that reduced sampling effort during the last sampling period of 2017. Mid-summer trapping yielded the most bees in each year compared to other collection periods, with >3-fold more bees collected during this period relative to other sampling rounds. Bees collected in the study represented five families, and their occurrence differed by fire severity of the sampling sites (Table 1). The most common genera collected were, in order of abundance, Halictus, Bombus, Lasioglossum, and Xylocopa (Table 1); they collectively accounted for 84% of the total catch. The trapped species of these four genera are all broad floral generalists (polyleges), and except for the one Xylocopa species, all are eusocial or primitively eusocial, producing annual colonies that multiply workers during the growing season (Michener 2007). A single species, H. tripartitus, comprised 1/3 of all bees in the trap catches. The honeybee (Apis mellifera) was the fifth most common genus collected but was excluded from generalized linear mixed model analyses because individuals likely originated from domestic hives.

Bee abundance and richness
Fire severity was a strong predictor of bee abundance in our study sites (RdNBR: F 1, 162 = 33.99, P = <0.001). We observed a positive relationship between fire severity and bee abundance throughout the flight season, which ❖ www.esajournals.org 7 April 2019 ❖ Volume 10(4) ❖ Article e02668 differed in magnitude between collection periods (RdNBR 9 Collection: F 3, 162 = 3.58, P = 0.015), and was greatest in the late summer. More specifically, we trapped 16.9-38.9% more individual bees for every 100-unit increase in site RdNBR, depending on when sampling took place (Appendix S1: Table S2). Based on predicted values from the generalized linear model with all site-level covariates held at their means ( Fig. 2A, B), this translated to an average of 20 times more individual bees collected at the highest severity burn sites compared to the least burned sites across the four collection periods. The relationship between fire severity and bee abundance was consistent across both years of sampling (RdNBR 9 year: F 1, 162 = 0.40, P = 0.528). Bee richness also increased with fire severity (F 1, 162 = 56.67, P = <0.001), though we did not detect an interaction between fire severity and collection period (RdNBR 9 collection: F 3, 162 = 4.62, P = 0.124). We observed a 18.3-26.6% increase in bee richness for every 100-unit increase in RdNBR (Appendix S1: Table S3). Based on predicted values from the model with all site-level covariates held at their means, this translated to an average of 119 more bee species in the highest severity burn relative to the lowest severity burn across collection periods (Fig. 2C,  D). We did not detect an interaction between fire severity and year (RdNBR 9 year: F 1, 162 = 1.18, P = 0.280).
Bee habitat characteristics were variable in the extent to which they were associated with fire severity (Fig. 3). As expected, canopy cover was greatly reduced in stands that burned more severely (R = À0.64, P = <0.001). Regarding nesting substrates, less bare ground remained with higher fire severity (R = À0.18 P = 0.006). The volume of dead wood (i.e., snags, stumps, and coarse woody debris) did not change along the fire severity gradient (R = À0.07, P = 0.701), whereas the number of woody pieces containing potential nesting cavities was positively correlated with fire severity (R = 0.32, P = 0.064).

Bombus vs. Halictus response to RdNBR
The two most common genera in the study, Bombus and Halictus, both increased in abundance along the fire severity gradient with an interaction between fire severity and genus (RdNBR 9 Group: F 1,33 = 3.63, P = 0.037). This interaction stemmed from differences in the rate at which the abundance of the two groups ( Notes: The proportion Relative differenced Normalized Burn Ratio (RdNBR) column shows the proportion of the sampled fire severity gradient where a given species was trapped. In the Family column, values for n are the total number of individuals collected within that family. The minimum and maximum RdNBR represent the most extreme values where each taxon was collected.

Bee and flowering plant community composition
We plotted variation in bee and plant communities along the fire severity gradient using genera only, as the species-level analysis was poorly fit based on high stress scores (>0.2). Fire severity was correlated with the ordination of bee genera (R 2 = 0.5, P < 0.001) and flower genera (R 2 = 0.3, P = 0.003). The bee genera with the best fit were two abundant (Bombus and Xylocopa) and five uncommon genera (Ceratina, Megachile, Melecta, Osmia, and Eucera). Of these genera, the relative abundances of Bombus, Ceratina, and Osmia were negatively correlated with burn severity, and the remaining genera were positively correlated with burn severity (Fig. 5). The floral genera with the best environmental fit were as follows: Berberis  Table S1 for a list of all sampling sites and corresponding RdNBR values.

DISCUSSION
Our study is the first to demonstrate that local wildfire severity serves as a strong predictor of bee diversity, with a 20-fold difference in abundance and an 11-fold difference in species richness across a natural gradient of forest wildfire severity. We detected an unexpectedly large number of wild bee species in recently burned forest patches, with >100 species/morphospecies mostly trapped in regions of the forest altered by severe, stand-replacing wildfire. These bees represented five families and a substantial portion of the estimated 500+ wild bee species in the state of Oregon (A. Moldenke, personal communication; Kincaid 2017), highlighting the role of disturbed patches within mixed-conifer forest in supporting essential habitat for a functionally important group of organisms.
The influence of fire severity we detected on bee populations was especially pronounced relative to fire effects observed in past studies. We detected a six times greater effect on bee species richness than that observed with time since wildfire in Mediterranean forests (Potts et al. 2003), and a 10 times greater effect on species richness than that observed with prescribed fire management in oak forest (Campbell et al. 2007). Based on our findings, the mosaic of habitats left by mixed-severity wildfires could explain why some studies have reported limited changes in bee species richness due to fire (e.g., Lazarina et al. 2017), while others have observed marked increases in bee species richness after fire in similar habitats (e.g., Potts et al. 2003). Without considering the spatial heterogeneity of mixedseverity wildfires, studies could considerably over-or underestimate the overall impact of wildfire on early seral-adapted organisms and their habitat.
In addition to the bee community, fire severity also had a strong effect on the flowering plants on which bees depend. Study sites that experienced the greatest fire severity had more open canopies and a minimum 18 times greater density of flowering plants throughout the bee flight season compared to the least severely burned stands. Spring bloom in our system was dominated by woody perennial shrubs, such as Ericaceous Gaultheria shallon and Vaccinium ovatum and Rosaceous Rubus parviflorus and Rubus ursinus. Late July through early September, bloom shifted to mostly annual or short-lived forbs, such as various Asteraceae species (e.g., Pseudognaphalium thermale, Stephanomeria elata, and Lactuca serriola), and invasive plants like Cirsium vulgare (Asteraceae) and Hypericum perforatum (Hypericaceae). Based on NMDS analyses, many of these plants were relatively more abundant in severely burned patches. Temporal changes in flowering plant availability due to fire disturbance have been previously recorded in grasslands (Mola and Williams 2018), and we demonstrate here that fire severity moderates this relationship within mixed-conifer forest.
The extended temporal availability of flowering plants in severely burned sites led to pulses of bees by the mid-to late summer when the neighboring forest understory had little or no bloom present. Eusocial Bombus and primitively eusocial Halictus produce annual colonies that become populated with workers as the season progresses, so they appear to have benefitted disproportionately from late-season floral resources. Marked increases in population size within these colonies explain why bee abundance-but not species richness-was predicted by the interaction between fire severity and season in our models. The dynamic response of flowering plants and bees to fire severity throughout the season highlights the importance of considering the temporal variability in ecosystem changes driven by mixed-severity fire.
Bee nesting substrates were largely maintained along the fire severity gradient, with only slight changes in bare ground availability and potential cavity-nesting locations. However, two aboveground-nesting genera (Ceratina and Osmia) and one genus composed of species that nest belowand aboveground (Bombus; Williams et al. 2014), comprised a larger proportion of the bee specimens collected at less severely burned stands. Of note, we did not collect any Ceratina at sites that experienced stand-replacing fire severity. The region's Ceratina nest in dead pithy stems of plants (e.g., Rubus ursinus; A. Moldenke, unpublished data, J. Rivers, personal observation), which burn easily and would likely take several years after a fire to regrow. In contrast, ground-nesting species place their nests in areas that are largely protected from heating mortality caused by wildfire Neff 2011, Love andCane 2016) and thus not expected to be impacted by wildfire in the context of our study.
The two most abundant genera we captured (Bombus and Halictus) both increased in abundance along the fire severity gradient, but the response was greater among Halictus species. Both genera primarily nest in the ground, although Bombus often use empty rodent burrows or, less frequently, nest aboveground (Williams et al. 2014) and Halictus burrows into exposed soil surfaces (Cane 1991(Cane , 2015. The greater increase in Halictus abundance in relation to fire severity may be due to the much shorter foraging distance observed in such small bees (Gathmann and Tscharntke 2002), which requires them to nest closer to flower-rich areas. The positive response of Halictus to fire severity thus provides additional evidence that the relationship between fire severity and bee abundance is reflective of enhanced local resources. In contrast, worker Bombus can forage several kilometers from their nest when local bloom is sparse (Redhead et al. 2016, Pope andJha 2018), so it would seem more likely that foragers would occasionally encounter traps placed in foragepoor patches. This helps to explain why Bombus made up a large proportion of bees collected at the less severely burned stands (based on NMDS results), but still increased in abundance with fire severity.
In this study, RdNBR was a useful metric for predicting bee diversity after heterogeneous wildfire. This metric, which is commonly used by researchers and land managers, is a better estimate of fire severity than other Landsat-based measures that do not correct for differences in pre-fire vegetative conditions (Miller and Thode 2007). However, burn severity does correlate with some pre-fire forest characteristics, such as stand age (Zald and Dunn 2018). Because we do not compare our sites to unburned control sites, we are not able to fully address how pre-fire forest characteristics might influence bee habitat. We were able to statistically control for stand age and other site-level characteristics in our model, but more research is needed to explore whether bees and other pollinators respond differently to fire severity within other gradients of stand age and in different forest types. The sites where we sampled in the Douglas Complex covered a broad RdNBR range representative of mixedseverity fires (Miller et al. 2009), so our study holds potential for making comparisons between fires in space and time.
Our study provides additional evidence that many pollinators benefit from disturbances in ❖ www.esajournals.org 15 April 2019 ❖ Volume 10(4) ❖ Article e02668 forest ecosystems (Hanula et al. 2016). Human activities like harvest and prescribed fire can also enhance pollinator diversity by creating patches of open habitat in otherwise closed-canopy forests (Taki et al. 2013, Rubene et al. 2015, Rodr ıguez and Kouki 2017. However, the landscape heterogeneity created by natural disturbances like wildfire is important for supporting pollinator diversity in space and time (Ponisio et al. 2016). The floral resources available to pollinators after wildfires may also be different than those that follow non-fire disturbances, as many plants in fire-adapted regions have adaptations like smoke-or heat-germinated seeds (Keeley 1987). Few studies have examined pollinator dynamics in managed forest systems (Rivers et al. 2018a), and more research is needed to contrast the pollinator response to human disturbances relative to wildfire, especially as disturbance regimes continue to shift worldwide (Turner 2010).
Ours is the first study to demonstrate the positive relationship between local wildfire severity and wild bee diversity in a fire-adapted mixedconifer forest, providing additional evidence that some organisms depend on stand-replacing wildfire to maintain their habitats (e.g., Hutto 2008, Tingley et al. 2016. Large, mixed-severity fires create heterogeneous landscapes that are critical for sustaining biodiversity in fire-adapted ecosystems (Kelly and Brotons 2017), but until this point we had little evidence of the spatial variability of pollinator communities within these disturbed habitats. Wild bees play a critical role in pollinating flowering plants, which in turn support wildlife habitat and biodiversity in many of these same ecosystems. Understanding pollinator response to wildfire characteristics is therefore a necessary step toward predicting changes in ecosystem function as wildfire activity continues to increase in many regions of the world (Westerling et al. 2006, Dennison et al. 2014, Jolly et al. 2015.