Study area
Galiano Island lies in the rain shadow of the Olympic Mountains and the Vancouver Island Ranges, in southern coastal British Columbia, Canada. This region is defined by its temperate Mediterranean-type climate, with mild, wet winters and warm, dry summers (Klassen et al. 2015). The combined effects of low precipitation, warm temperatures, and high sunshine hours result in an annual moisture deficit during summer months, which varies to some extent depending on precipitation (Moore et al. 2010). These conditions are expected to become more extreme under projected climate change scenarios: increases in winter precipitation and the intensity and duration of seasonal drought are likely in this system (Salathé et al., 2008; Seager et al., 2019; Spies et al., 2010).
Galiano Island is relatively intact ecologically, with about 24% of its land base conserved in protected areas and a high percentage of forest cover, hosting vegetation communities typical of the Coastal Douglas-fir Biogeoclimatic Zone (Nuszdorfer et al., 1991). Although the island has historically been subject to extensive industrial forestry practices, only ~11% of the landscape has been converted for present day human use (Madrone 2008). This largely forested study area thus presents a contrast to the fragmented environments previously studied by pollinator researchers in water-limited ecosystems (e.g., McFederick and LeBuhn 2005; Wray and Elle 2015).
Sampling
This study was based on a 2×2 factorial study design contrasting four conditions of disturbance and soil moisture availability: 1) dry semi-natural environments (woodland and associated rock outcrop communities); 2) wet semi-natural environments (wetlands); 3) dry modified environments (disturbed upland areas such as clear-cuts and hydro-line corridors); and, 4) wet modified environments (gardens, orchards and fields). Sites were stratified based on available terrestrial ecosystem mapping data (Madrone, 2008), with 6 sites selected per condition (for a total of 24 sites). Limitations in existing site conditions and the logistics of site access resulted in an imbalance in the study design, with 4 sites representing the wet seminatural condition versus 8 sites representing the dry seminatural condition (and 6 sites representing the other two conditions).
Terrestrial ecosystem mapping data were used to obtain estimates of land cover falling within a 500 m buffer around each site, orienting around polygon centroids, including natural, modified, and forested matrix habitat (Fig. 1). Natural matrix habitat included coniferous and deciduous forest cover at various structural stages, as well as wetlands, woodlands, meadows, and rock outcrops. Modified matrix included rural, agricultural and developed land cover. Sites ranged in size from 0.21 to 6.3ha (Supplementary Information). While grouped site conditions circumscribed similar habitat types based on common soil moisture regimes, disturbance conditions were diverse, including forestry, fire, landscaping, and other anthropogenic effects. Commonalities between site conditions (henceforth habitat types) are reflected in the similarity of vegetation communities, as shown in ordination plots (Supplementary Information). Additional environmental data were collected along transects at the site level, including estimates of canopy cover, slope, moss cover, coarse woody debris, and bare soil—variables known to correlate with bumble bee nesting habitat (Potts et al. 2005; Wray and Elle 2015).
To sample floral resource availability (FRA), we stratified habitat types and randomly distributed 6–8 (2 × 15 m) belt transects throughout each site, using balanced acceptance sampling methods to ensure a random yet spatially balanced distribution of transects (van Dam-Bates et al., 2018). The number of transects was scaled roughly in proportion to the size of each site, with the aim of capturing variability in floral resource availability across the landscape. Floral resource availability was quantified as counts of flowering shoots, recorded for each plant species at 1-m intervals, with each interval surveyed comprising a 1 × 2 m area spanning both sides of the transect line. Soil moisture was recorded at 5-m intervals as volumetric water content (%VWC) using a Field Scout Time-Domain Reflectometry probe. Sampling was conducted on a monthly basis from April through August 2018, resulting in five samples per site, with each sample period covering 13 days.
To sample bumble bees, three blue vane traps (Stephen and Rao 2005) were placed at roughly equidistant intervals across each site. Traps were positioned to ensure their visibility given surrounding vegetation and landscape features, and to minimize public interference. Bumble bees captured by blue vane traps were pooled into one sample per site, resulting in 24 × 5 (120) samples divided into the four habitat types. Blue vane traps were vandalized in a dry semi-natural site in May, resulting in the loss of one sample (n=120-1=119). Bumble bees were identified following Williams et al. (2014).
Modelling
Statistical analyses were implemented in R Version 3.6.0 (R Core Team 2019). Linear mixed effects models (LMMs) were fitted incorporating log-transformed soil volumetric water content (logVWC) as a response to contrasting habitat types (fixed effects), with transects nested within sites (random effects), to estimate differences in soil moisture between habitat types. To estimate differences in floral resource availability (FRA), we fit negative binomial generalized linear mixed effects models (GLMMs) incorporating counts of flowering shoots as a response to contrasting habitat types (fixed effects) and transects nested within sites (random effects) for discrete sample periods. We also modeled FRA as a seasonal pattern, by fitting a GLMM incorporating time (julian days) as a quadratic term to estimate overall differences in FRA across site conditions. While time emerged as a significant term in these models, the quadratic term did not adequately describe the observed seasonal patterns in FRA, resulting in a poor fit. Hence, our analyses have focused on discrete sample periods instead. Negative binomial generalized linear mixed effects models were also fitted to model mean bumble bee abundance (counts) as a function of habitat type and sample periods (fixed effects), with a random intercept term corresponding to site. Other environmental variables incorporated as predictors included proportions of natural and modified matrix, forest cover, slope and moss cover. Significant terms were evaluated based on model summaries, and top models selected based on AIC scores (see below).
We estimated differences in the total count of bumble bees sampled between habitat types for each sample period separately. Because our focus was on whether certain habitat types supported prolonged foraging under conditions of drought stress, we concentrated primarily on bumble bee worker abundance as a proxy for bumble bee foraging activity, though differences were also considered among drones and queens—especially late-flying queens, whose presence in the environment may be considered a proxy for the reproductive success of colonies. We fitted GLMs to estimate differences in the abundance of bumble bees within isolated sample periods, as temporal autocorrelation was not problematic in these analyses. We also modelled the abundance of each species as a response to habitats and other environmental variables, both across sample periods (GLMMs) and within sample periods (GLMs). Linear mixed effects models were implemented using the R package ‘nlme’ (Pinheiro et al. 2021); GLMs were implemented using the R package ‘MASS’ (Venables and Ripley 2002); and GLMMs implemented using ‘lme4’, and ‘glmmTMB’ in cases where zero-inflation proved problematic (Bates et al. 2015, Brooks et al. 2017). The best models were selected based on AIC test scores, assuming ΔAIC of 2.0 as a threshold for model improvement (Burnham and Anderson 2002). Model effects are reported as Incidence Rate Ratios (IRR), and as marginal R2 values (variance explained by fixed effects) for GLMMs. For GLMs, pseudo (Nagelkerke) R2 values are presented for significant terms. Marginal effects were calculated using ‘sjPlot’ (Lüdecke 2018).