Influence of a local bioclimatic gradient on wild bees in boreal forest clearcuts

Abstract

Variation in wild bee assemblages along climatic gradients is still poorly known, particularly in eastern Canada. In this study, we assessed wild bee communities in forest stands that were recently clearcut along a latitudinal bioclimatic gradient. The analyses were conducted at two taxonomic levels: first at the genus level for bees as a whole (Anthophila) and then at the species level for bumble bees (Bombus spp. Latreille) (Hymenoptera: Apidae). Coloured pan traps were used to catch bees and a point-intercept method was used to survey floristic composition at each site. Abundance and richness of Anthophila genera and abundance of bumble bee species were inversely related to latitude, but richness of bumble bee species was not, being instead associated with sections of the gradient. These results indicate that the relationship between wild bee diversity and latitude varies among taxa and that climate and floristic resources probably each play a role. Bombus terricola, a species of concern in eastern Canada, was found in relatively high numbers throughout the study area. To our knowledge, this study provides the first assessment of wild bees in clearcuts in the province of Quebec, showing a clear differentiation of the communities along a spatial bioclimatic gradient.

Résumé

Les variations dans les assemblages d’abeilles sauvages le long de gradients climatiques sont encore mal connues, en particulier dans l’est du Canada. Dans cette étude, nous avons évalué les communautés d’abeilles sauvages dans des peuplements forestiers ayant été récoltés récemment le long d’un gradient bioclimatique latitudinal. Les analyses ont été menées à deux niveaux taxonomiques : d’abord au genre pour l’ensemble des abeilles (Anthophila), puis au niveau de l’espèce pour les bourdons (Bombus spp. Latreille) (Hymenoptera : Apidae). Des pièges bols colorés ont été utilisés pour capturer les abeilles et une méthode par points d’interception a été utilisée pour étudier la composition floristique de chaque site. L’abondance et la richesse spécifique des genres d’Anthophila et l’abondance des bourdons étaient inversement liées à la latitude, mais pas la richesse spécifique des espèces de bourdons, ceux-ci étant plutôt associées à des sections du gradient. Ces résultats suggèrent que la relation entre la diversité des abeilles sauvages et la latitude varie selon les taxons, et que le climat et les ressources floristiques jouent probablement tous les deux un rôle. Bombus terricola, une espèce dont le statut est préoccupant dans l’est du Canada, a été trouvé en nombre relativement élevé dans toute la zone d’étude, ce qui pourrait s’expliquer en partie par le fait que les sites étudiés sont situés loin des activités humaines autres que l’exploitation forestière. À notre connaissance, cette étude fournit une première évaluation des abeilles sauvages dans les coupes à blanc de la province de Québec, démontrant une nette différenciation des communautés selon un gradient bioclimatique spatial.

Introduction

Insect pollinators are crucial for the successful sexual reproduction of more than 80% of the world’s flowering plants (Ollerton et al. 2011; Kremen 2018). The pollination services provided are vital to maintaining biodiversity in natural ecosystems and also for human food security (Ollerton 2017). Insect pollinators are known to increase the production of a considerable number of crop species (Klein et al. 2006; Vanbergen and the Insect Pollinators Initiative 2013; Kremen 2018). Bees (Apoidea: Anthophila) are considered the most important pollinators, with 20 000 species (Michener 2007) that vary in their use of floral resources and are distributed across many ecosystems worldwide.

Bees are facing major threats, such as climate change and habitat loss caused by disturbances (Sánchez-Bayo and Wyckhuys 2019; Wagner 2020). Evidence of bee decline has been recorded in recent decades across Europe and North America (Williams and Osborne 2009; Cameron et al. 2011, 2016; Bartomeus et al. 2013; Mathiasson and Rehan 2019; Soroye et al. 2020). Multiple climatic variables impact bee physiology (Wagner 2020) and could be responsible in part for these declines. Temperature is one determining factor in bees’ development rate, survival, reproductive success, flight activity, and foraging patterns (Rader et al. 2013; Whitfield et al. 2013). Climatic variables also influence the productivity and dispersion of flowering plants (Pyke et al. 2016). Variations in temperature and precipitation can reduce phenological synchrony between bees and plants, resulting in declines in the abundance of some species and in shifts in species’ ranges (Williams et al. 2014; Pyke et al. 2016).

Human activities such as the exploitation of natural resources can drastically reshape natural environments in ways that affect bees (Memmott et al. 2004). Landscape disturbances leading to changes in the availability of flower resources, nesting sites, and overwintering habitats particularly can impact bees, eventually resulting in species loss (Zurbuchen et al. 2010a; Williams et al. 2014). Because bee species vary in their nesting requirements, these impacts are expected to be species specific.

The boreal forest is the northernmost forested ecosystem on the planet. It is one of the widest spanning biomes on earth, forming a mostly uninterrupted belt across Canada, Alaska, Russia, Finland, and Scandinavia (Henry 2002). It is characterised by a cold climate that allows the growth of mostly coniferous trees. Disturbances such as fires and defoliating insect outbreaks are common in the boreal forests and contribute to their natural cycle (Shugart et al. 2005). Logging is also a major anthropogenic disturbance type affecting boreal ecosystems. Even though bees are often studied in agricultural or urban environments, open or recently disturbed forests that are rich in understorey flowering plants constitute an important habitat for many bee species. Anthropogenic disturbances such as logging create early succession habitats that are potentially suitable for many species (Romey et al. 2007; Pengelly and Cartar 2010; Taki et al. 2013; Zitomer et al. 2023). Several studies have shown that vegetation development following logging activities can positively impact wild bee communities in temperate forests (Romey et al. 2007; Taki et al. 2013; Rubene et al. 2015; Rodríguez and Kouki 2017; Chase et al. 2023a, 2023b; Davies et al. 2023). Specifically, the creation of canopy openings increases light penetration, which stimulates understorey herbaceous flowering plants (Hanula et al. 2015). Because they feed exclusively on pollen and nectar at every developmental stage (Whitfield et al. 2013), bees are highly dependent on the availability of floral resources. The increased amounts of woody debris, dead plant materials, and exposed soil after forest harvesting also create specific nesting habitats that can be used by ground-, cavity-, and wood-nesting bees (Hanula et al. 2016; Chase et al. 2023a). The presence in clearcuts of a wide variety of nesting habitats and floral resources benefits bees because they must return to their nest several times a day while foraging (Zurbuchen et al. 2010a).

Studies conducted on wild bees in managed boreal forests have mainly focused on the effect of elevation, floral and nesting resources, or management practices (Rubene et al. 2015; Rodríguez and Kouki 2017; Westerfelt et al. 2018; Luna Santa-María 2021). How climate drives the assemblage of bee communities along latitudinal gradients in forest environments remains unknown, particularly at the regional scale. Specifically, bee communities present in recent clearcuts may vary depending on climatic and environmental conditions present in different regions. A better understanding of how bee communities vary along spatial bioclimatic gradients could help forecast potential impacts of climate change. Indeed, current climate change models predict that an increase in both average temperatures and their variability, as well as decreased seasonality, will be more pronounced at higher latitudes in the Northern Hemisphere (Intergovernmental Panel on Climate Change 2014). Wild bees, which are nonmigratory organisms, may be strongly affected by these changes (Pyke et al. 2016), especially given their limited heat tolerance (Gonzalez et al. 2024). Negative consequences of climate change have been documented at the community and population levels (Rafferty 2017; Kammerer et al. 2021), especially for bumble bees (Jackson et al. 2022; Liczner et al. 2023). The consequences include reduction of species-distribution ranges (Sirois-Delisle and Kerr 2018), reductions in floral-resource availability (Phillips et al. 2018), and phenological mismatches between bees and their feeding resources, thereby reducing the breadth of their diet (Memmott et al. 2007; Forrest 2015). Assessing wild bee communities’ potential response to spatial climate variation would facilitate evaluation of their potential sensitivity to climate change (Lebuhn et al. 2013).

The overall aim of the present study was to, respectively, assess the effects of climate and floristic composition on wild bee diversity in clearcuts in the province of Quebec along a single local latitudinal bioclimatic gradient. We used recent clearcuts as a study ecosystem to uniformise the disturbance history of the different sites along this gradient and to facilitate the analysis of bee community responses to the variables of interest, climate and flowering plants.

We first assessed Anthophila communities at the genus level and then pursued our evaluation to the species level specifically for bumble bee communities (Bombus spp. Latreille) (Hymenoptera: Apidae). We conducted species-level analyses for this genus because it contains species that are potential indicators of climate change (Williams et al. 2014; Kerr et al. 2015). We hypothesised that wild bee diversity would be inversely related to latitude; according to McCabe et al. (2019), cold environments can act as environmental filters on bee communities. Moreover, as the seasonal activity period is variable among Anthophila taxa (Packer et al. 2007), a higher diversity of genera is expected in regions where the growing season is longer (Ogilvie and Forrest 2017). To test our hypothesis, we first compared bee abundance, richness, and diversity for three forest areas located in different bioclimatic domains in Quebec, Canada. Second, we highlighted indicator, dominant, and exclusive taxa using bee distribution patterns along the same gradient. Third, we explored the influence of local floristic composition on bee community assemblages. These steps also allowed us to identify which Anthophila genera and which bumble bee species were the most abundant in boreal and mixed forests in the province of Quebec, a region where post-disturbance bee communities have been little studied to date.

Material and methods

Study sites

We conducted our study during the summer of 2020 in three study areas positioned in three adjacent bioclimatic domains located at the interface between the boreal and temperate forest zones in the province of Quebec, Canada (Fig. 1): the balsam fir–yellow birch domain (southern area), the balsam fir–white birch domain (central area), and the black spruce–moss domain (northern area; Berteaux et al. 2014). Besides the dominant tree species representative of each bioclimatic domain – Abies balsamea (Pinaceae), Betula alleghaniensis (Betulaceae), B. papyrifera (Betulaceae), Picea mariana (Pinaceae) – boreal species aspen, Populus tremuloides (Salicaceae), and jack pine, Pinus banksiana (Pinaceae), tend to be sparsely present across the three study areas, whereas relatively thermophilic species such as red maple, Acer rubrum (Sapindaceae), yellow birch, B. alleghaniensis, and sugar maple, A. saccharum (Sapindaceae), are found only in the southern area. Forests in these regions are managed, generally using clearcutting, followed by natural regeneration of the harvested stands, with tree planting being used only when natural regeneration is insufficient. Partial cuts are also used, primarily in tolerant hardwood stands of the southern area.

Figure 1. Distribution of study sites (black triangles; n = 33) across three bioclimatic domains in the province of Quebec, Canada.

Thirty-three sampling sites were established within 50 km of three towns, each located in a different bioclimatic domain: 12 sites around La Tuque (47.43° W, 72.78° W), 10 sites around Dolbeau-Mistassini (48.88° N, 72.23° W), and 10 sites around Chibougamau (49.91° N, 74.36° W). These towns, located at least 150 km apart, were chosen because of their proximity to recent logging activities and their geographical location along the gradient. Each sampling site was positioned inside a forest stand that had been clearcut 3–10 years before. The sites were positioned at least 5 km from power lines and agricultural lands to avoid any bias related to their floral resources or pesticides. The sites were also positioned at least 0.5 km apart (generally more than 1 km) to ensure a relative spatial independence in wild bee communities (Zurbuchen et al. 2010b). In general, all sites were located on glacial tills with loamy textures and that were more than 50 cm thick over bedrock. On average, for the three areas, from south to north, mean growing season lengths were 215.86 days, 203.43 days, and 179.67 days; mean annual temperatures were 2.93 °C, 1.78 °C, and –4.2 °C; and mean annual precipitations were 958.9 mm, 925.79, and 921 mm, respectively, for the 2014–2020 period (Environment and Climate Change Canada; https://climat.meteo.gc.ca). Mean elevation for the sites of the southern, central, and northern areas was 354 m, 257 m, and 403 m, respectively. Further details on site characteristics, including clearcut year, clearcut size, and ecological type, are provided in Supplementary material, Table S1.

On each site, we established a 2000-m² plot in which two parallel 40-m-long transects were positioned 20 m apart. To minimise a potential edge effect, a buffer of approximately 50 m was kept around plots from the nearest active secondary road and from the forest margin. Following plot installation, areas were surveyed every two weeks from mid-June to the end of August. During each survey, all the sites of an area were visited within 2–3 consecutive days. At each visit to a plot, two sampling protocols were performed, starting with the gathering of captured insects and followed by a survey of floral resources. In the southern and central study areas, five sampling surveys were performed over the course of the whole summer. In the northern area, only four surveys were carried out because of the shorter growing season.

Wild bee sampling

Bees were collected using coloured pan traps (Fig. 2). We used blue, white, and yellow 12-oz plastic PartyCity party bowls (7.5 inches diameter × 1.5 inches deep; ∼19 cm × ∼4 cm). These colours are known to attract bees and are the most used in North American pan trap sampling (Droege et al. 2010; Whitfield et al. 2013). Three pan traps (one of each colour) were set 10 m apart along the two transects, for a total of six coloured pan traps per plot. For each trap, binder clips were used to secure a bowl approximately 20 cm above the ground in a tomato cage. We prevented rain from filling the bowl and spilling its contents by fastening a transparent plastic saucer on top of the cage, about 20 cm above the bowl. Approximately 120 mL of propylene glycol was used to catch and preserve specimens in the bowls because it does not evaporate like soapy water and slows the decomposition of dead insects (Whitfield et al. 2013). This also allowed us to leave the traps operational in the field between visits, thereby accounting for the fact that active periods are variable between bee species and can be affected by daily weather (Packer et al. 2007). At each visit, the traps were reinstalled and refilled with propylene glycol after the caught specimens had been collected from the bowls. This way, traps were continually deployed for the entire growing season, from mid-June to the end of August.

Figure 2. Pictures of the pan-trapping setup shown side-by-side, showing blue and yellow 12-oz plastic bowls secured in a tomato cage with binder clips. Transparent plastic saucers fastened on top of the cages prevented rain from filling the bowls and spilling their contents.

Before identification, captured specimens were stored in 40% ethanol. They were then washed and pinned. All specimens were identified at least to the genus level. Identification was pursued to the species level for the Bombus spp. specimens. Packer et al. (2007), Williams et al. (2014), Wilson and Carril (2015), Ascher and Pickering (2020), and Milam et al. (2020) were used as keys. Difficult bee genera and Bombus species were confirmed with the assistance of expert taxonomists (see Acknowledgements). Voucher specimens were deposited at Université Laval (laboratory of Valérie Fournier, Québec, Quebec, Canada) and at the Collection d’insectes du Québec (Ministère de l’Agriculture, des Pêcheries et de l’Alimentation du Québec, Québec, Quebec, Canada).

Floral-resource sampling

The floral resources were surveyed following the FIREMON point-intercept sampling method (Caratti 2006). Two transects of 18 sampling rods (3 mm diameter) were established along those of the pan traps, with 2 m between each, for a total of 36 rods per site. Along the transects, we considered each flowering plant species touching a pole as a “hit.” The point-intercept method allowed us to statistically consider the seasonal species turnover in floral composition.

Relative species abundance $Re{l_{abond}}$ of floristic species $i$ at a site $j$ was obtained using the rule of proportion:

$$Re{l_{abon{d_{ij}}}} = {\mu _{Ob{s_{ij}}}}\sum\nolimits_i {{\mu _{Ob{s_{ij}}}}} \times 100$$

where the mean number of observations ${\mu _{Obs}}$ of each floristic species $i$ encountered within the site $j$ was calculated as:

$${\mu _{Ob{s_{ij}}}} = \sum O b{s_{ij}}\left( {36\;rods\; \times \;{k_j}} \right)$$

where $Ob{s_i}$ corresponds to the number of hits for species $i$ within site $j$ and ${k_j}$ corresponds to the number of surveys that could be done on that same site.

Flowering plants were identified in the field to the lowest possible taxonomic level using the keys of Frère Marie-Victorin et al. (2002), Larivière and Villeneuve (2014), Leboeuf (2016), and Brouillet (2024). The recurring visits allowed us to confirm the identification of some annual species that were insufficiently developed in early spring.

Data analyses

All analyses were first performed at the genus level for the whole Anthophila group, including bumble bees, and then performed a second time for bumble bees only. Data from the whole sampling season were pooled to best assess the factors influencing spatial variation in wild bee populations. Sites are considered as replicates for the three bioclimatic areas. The effectiveness of the sampling effort was assessed by computing species accumulation curves in the Anthophila data and in the bumble bee data using the function specaccum from the vegan package, version 2.5-7, in R (Oksanen et al. 2020; R Core Team 2021). Accumulation curves assess the increase in the number of species observed according to the sampling effort. Mean species accumulation curves and their standard deviations were obtained from repeated richness accumulations in random order from 999 permutations or subsampling without replacement (Gotelli and Colwell 2001).

Assessing community diversity

Abundance (number of individuals) and richness (number of genera for Anthophila; number of species for bumble bees) were calculated for each site and were used to compare community diversity between the three areas. Shannon’s H’ diversity index (Hill 1973) and Pielou’s J’ evenness index (Pielou 1975) were similarly calculated for each site and were used as diversity indexes to quantify how diversified and equally distributed the communities of wild bees were across sites and along the bioclimatic gradient. The two indexes are complementary, which ensures thorough interpretation of the results. Indexes were computed using the function diversity from the vegan package (Oksanen et al. 2020). Shannon’s diversity index considers the total number of species (S) and the individual abundance within each species. Compared to a community dominated by a sole species (minimal value of 0), a community with perfectly codominant species will have a maximum index value (lnS, LogS, or Log₂S, depending on the choice of the logarithm base) and will be considered the most diversified (Hill 1973). Pielou’s evenness index indicates the extent to which species are proportional in numbers on each site along the gradient; a value of 1 indicates a perfectly balanced community (same number of individuals for each species), whereas a value of 0 indicates an unbalanced community where only one species surpasses all the others in numbers (Pielou 1975).

Abundance and richness were compared between areas using generalised linear models with the glm function from the stats package in R (R Core Team 2021), using a log-link function in both models. Overall significance (P < 0.05) of the categorical factor “area” was assessed by performing a likelihood ratio test by comparing the model containing that variable with the null model. For cases where the difference was significant, pairwise comparisons were computed by changing the base level. Shannon’s diversity and Pielou’s evenness indexes were compared between areas using one-way analyses of variance with the lm function from the stats package, version 4.1.1, in R. When results were significant, multiple pairwise comparisons were computed using the joint_tests and emmeans functions from the emmeans package and the cld function from the multcomp package in R (Hothorn et al. 2008; Lenth 2021). Shannon’s diversity data were log-transformed to be normally distributed.

Indicator, dominant, or exclusive taxa

The Dufrêne–Legendre indicator species analysis (IndVal index) was used to determine if any and, if so, which Anthophila genera and bumble bee species were significantly associated with any of the studied bioclimatic domains along the latitudinal gradient. The IndVal index calculates the indicator value of species within clusters of sites by combining their specificity (highest mean abundance within a group compared to the other groups) and fidelity (presence in most sites of that group; Dufrêne and Legendre 1997). The index values were computed using the indval function from the labdsv package, version 2.0-1, in R (Roberts 2019). The IndVal index was also applied to the floristic data to corroborate the observed floristic assemblage within each area with published findings.

To highlight the differences in Anthophila abundance among sites and along the gradient, we also plotted heatmaps with the pheatmap package in R (Kolde 2019). The heatmaps allowed us to visualise the distribution of highly dominant taxa and to identify Anthophila genera and bumble bee species that were exclusively captured in specific parts of the gradient.

Influence of floristic composition

We used a redundancy analysis, a constrained ordination (Legendre et al. 1998), to determine the effect of floristic composition on Anthophila and bumble bee communities along our gradient. The Hellinger transformation was applied to our bee matrix using the decostand function from the vegan package in R (Oksanen et al. 2020). The Hellinger transformation is particularly well suited for species abundance data because it gives low weights to variables with low counts and many zeros (Legendre and Gallagher 2001). The tested explanatory variables were floristic species, based on their relative abundance on each site. Colinearity of our redundancy analysis model was assessed using variance inflation factor scores. The step function (from the vegan package) was used with forward selection, based on both the Akaike information criterion and the Bayesian information criterion, to determine which floristic species best explain bee community composition. We performed an analysis of variance–like permutation test on our redundancy analysis to assess the significance of constraints and axis using the anova.cca function from the vegan package. Anthophila genera and bumble bee species with fewer than 10 individuals were not used in the analysis, given their scarcity. Floristic species with a relative abundance less than 1% were also omitted.

Results

General results

Overall, we retrieved a total of 328 samples from the pan traps in the southern area, 263 from the central area, and 264 from the northern area. The number of consecutive trap days for each area was 76, 71, and 55, respectively. We collected a total of 4255 specimens belonging to the Anthophila group, corresponding to 21 genera in five families (Table 1). The most abundant genus was Lasioglossum spp. Curtis (Hymenoptera: Halictidae), with a total of 2760 individuals, representing 65% of all collected bees. Bumble bees (Bombus spp.) represented the second most abundant genus (20%), with 863 individuals. All other Anthophila genera combined represented only 15% of all collected individuals. Within bumble bees, 10 species were found (Table 2). Bombus ternarius Say was the most abundant, with 339 individuals, representing 39% of all caught bumble bees. Bombus terricola Kirby was the second most abundant, with 173 individuals (20%), followed by B. sandersoni Franklin, with 170 individuals (19.7%), and B. vagans Say, with 93 individuals (10.78%).

Table 1. List of the 21 bee genera collected along a bioclimatic gradient in the province of Quebec, Canada, from early June to late August 2020, with their abundance in each studied area. The genera are sorted according to their corresponding family.

Table 2. List of the 10 bumble bee species collected along a bioclimatic gradient in the province of Quebec, Canada, from early June to late August 2020, with their abundance for each studied area. Species are presented in alphabetical order.

Sampling effort

The accumulation curves for both taxonomic levels, Anthophila genera (Supplementary material, Fig. S1A) and bumble bee species (Supplementary material, Fig. S1C), reached saturation, which indicates a sufficient sampling effort. Accumulation curves by area reached a plateau only in the southern area for Anthophila genera (Supplementary material, Fig. S1B) and in the northern area for the bumble bee species (Supplementary material, Fig. S1D), which suggests that richness in the other areas might have been slightly higher if more sampling sites were established.

Community diversity

Anthophila communities

For the Anthophila communities, diversity measures decreased from the south to the north of the gradient. According to the results of our generalised linear model analyses, the effect of area was significant in both the models for abundance and richness. Abundance (P-value < 0.0001; Fig. 3A) and genera richness (P-value < 0.0001; Fig. 3B) both differed significantly between all three areas. Total abundance was higher in the southern area, with 2616 individuals, followed by the central area, with 1027 individuals, and the northern area, with 612 individuals. A mean abundance of 218 ± 79.09 (standard deviation) was obtained for the southern area, 102.7 ± 44.16 (standard deviation) for the central area, and 55.64 ± 23.59 (standard deviation) for the northern area. Total genera richness was higher in the southern area (20 genera), followed by the central area (18 genera), and the northern area (13 genera). A mean richness of 11.58 ± 2.84 (standard deviation) was obtained for the southern area, 8.2 ± 2.57 (standard deviation) for the central area, and 5.9 ± 1.3 (standard deviation) for the northern area.

Figure 3. A, Mean abundance; B, mean genera richness; C, mean Shannon’s diversity index; and D, mean Pielou’s evenness index of Anthophila from clearcut sites in three studied areas along a bioclimatic gradient in the province of Quebec, Canada. Letters represent the post hoc result (P-value < 0.0001), showing significant differences between the three areas for diversity measures (abundance and richness) and for Pielou’s evenness index. Shannon’s diversity index data were log-transformed to be normally distributed and to respect homogeneity of variance.

Regarding diversity indexes, Shannon’s diversity for Anthophila communities did not differ significantly between the areas (F _2,30 = 1.354, P-value = 0.2736; Fig. 3C). Mean Shannon’s diversity was 1.07 ± 0.31 (standard deviation) in the southern area, 0.95 ± 0.4 (standard deviation) in the central area, and 1.18 ± 0.22 (standard deviation) in the northern area. However, Pielou’s evenness index differed significantly between the northern area and the two other areas (F _2,30 = 24.029, P-value < 0.0001; Fig. 3D), which were themselves not significantly different from one another. Anthophila community evenness was highly unbalanced along the whole gradient (J’ < 0.5). Mean Pielou’s evenness was 0.19 ± 0.05 (standard deviation) for the southern area, 0.2 ± 0.08 (standard deviation for the central area, and 0.36 ± 0.06 (standard deviation) for the northern area.

Bumble bee communities

For the bumble bee communities, the generalised linear model analyses indicate that the effect of area was significant only for the abundance model (P-value = 0.0247; Fig. 4A). Bumble bee species richness showed no relationship to area (P-value = 0.6756; Fig. 4B). Abundance was higher in the southern area than in the northern area (P-value = 0.0070). A total abundance of 345 individuals was obtained in the southern area, 265 individuals in the central area, and 253 individuals in the northern area. Mean abundance of bumble bees was 28.75 ± 17.14 (standard deviation) for the southern area, 26.5 ± 34.97 (standard deviation) for the central area, and 23 ± 9.24 (standard deviation) for the northern area. Total species richness was similar between areas, with nine species in the southern area, eight species in the central area, and seven species in the northern area. A mean richness of 4.75 ± 1.42 (standard deviation) was obtained for the southern area, 4 ± 1.89 (standard deviation) for the central area, and 4.64 ± 0.81 (standard deviation) for the northern area.

Figure 4. A, Mean abundance; B, mean species richness; C, mean Shannon’s diversity index; and D, mean Pielou’s evenness index of bumble bees (Bombus spp.) from clearcut sites in three studied areas along a bioclimatic gradient in the province of Quebec, Canada. Letters represent the post hoc result (P-value < 0.0001), showing significant differences between the southern and northern areas for abundance.

Diversity indexes showed no relation between bumble bee communities and area. Neither Shannon’s diversity index (F _2,30 = 1.598, P-value = 0.219; Fig. 4C) nor Pielou’s evenness index (F _2,30 = 2.864, P-value = 0.0727; Fig. 4D) differed significantly between areas. Mean Shannon’s diversity was 1.27 ± 0.21 (standard deviation) for the southern area, 1.06 ± 0.38 (standard deviation) for the central area, and 1.13 ± 0.24 (standard deviation) for the northern area. Mean Pielou’s evenness was 0.46 ± 0.77 (standard deviation) for the southern area, 0.49 ± 0.12 (standard deviation) for the central area, and 0.39 ± 0.09 (standard deviation) for the northern area.

Indicator, dominant, and exclusive taxa

Anthophila communities

Results from the IndVal analysis of Anthophila identified 10 genera as indicators in the southern area (Ceratina Latreille, Augochlorella Sandhouse, Andrena Fabricius, Lasioglossum, Megachile Latreille, Halictus Latreille, Melissodes Latreille, Osmia Panzer, Nomada Scopoli, and Augochlora spp. Smith), with Ceratina spp. having the highest value (Table 3). The central and northern areas had a single indicator genus each, Colletes spp. Latreille and Sphecodes spp. Latreille, respectively.

Table 3. List of indicators obtained from the Dufrêne–Legendre indicator species analysis for Anthophila genera and bumble bee species from clearcuts along a bioclimatic gradient in the province of Quebec, Canada. Studied areas were clustered. Taxonomic groups are sorted according to their indicator value and area.

According to the abundance heatmap (Supplementary material, Fig. S2), Lasioglossum and Bombus spp. were predominant genera along the whole gradient. Lasioglossum had the highest abundance at every site except one, where it was surpassed by the genus Bombus. Lasioglossum spp. was particularly abundant in the southern area and the first half of the central area (Supplementary material, Fig. S2A). Distribution patterns of the 19 other genera (Supplementary material, Fig. S2B) showed that, overall, the three study areas shared 11 genera. Six genera were found in the southern and central areas but not in the northern area; these were Anthophora Latreille, Ceratina, Dufourea Lepeletier, Epeolus Latreille, Melissodes, and Stelis spp. Panzer). Two genera – Augochlora and Augochlorella spp. – were collected exclusively in the southern area, and the only specimen of Pseudopanurgus spp. Cockerell was caught in the central area. Nomada spp. was collected in limited numbers in both in the southern and northern areas but was absent from the central area.

Bumble bee communities

For the bumble bee communities, results from the IndVal index analysis identified B. terricola as the only indicator species in the southern area (Table 3). Bombus frigidus Smith and B. flavidus Eversmann were indicator species in the northern area. No indicator species was detected in the central area.

According to the abundance heatmap of bumble bee species (Supplementary material, Fig. S2C), B. ternarius was a dominant species along the latitudinal gradient, as were B. sandersoni and B. terricola, which were also well distributed along the whole gradient. Bombus terricola was generally dominant in the southern area, whereas B. sandersoni tended to dominate in the central and northern areas. Overall, the three areas shared six species. Two species were collected only in a specific area – B. frigidus in the northern area and B. rufocinctus Cresson in the southern area. Bombus perplexus Cresson and B. vagans were found in both the southern and central areas but not in the northern area.

Influence of floristic composition

Floristic composition along the gradient

A total of 51 floristic species were recorded during our sampling season (Supplementary material, Table S2). The most abundant species along the whole gradient were Kalmia angustifolia Linnaeus (Ericaceae) (10.35%), Vaccinium angustifolium Aiton (Ericaceae) (10.12%), Prunus pensylvanica Linnaeus fils (Rosaceae) (9.76%), Diervilla lonicera Miller (Caprifoliaceae) (9.44%), and Rubus idaeus Linnaeus (Rosaceae) (8.61%).

In the southern area, 36 floristic plant species were recorded. Floristic composition of sites within the area was dominated by Prunus pensylvanica (16%), Rubus idaeus (15.56%), and Diervilla lonicera (11.41%). Five floristic species were significant indicators (IndVal index, P-value < 0.05) of the southern area: Prunus pensylvanica, Rubus idaeus, Eurybia macrophylla (Linnaeus) Cassini (Asteraceae), Maianthemum canadense Desfontaines (Asparagaceae), and Fallopia cilinodis (Michaux) Holub (Polygonaceae). See Supplementary material, Table S3 for the complete indicator species list with values.

In the central area, 36 floristic plant species were recorded. Floristic composition of sites within the area was dominated by Kalmia angustifolia (13.55%), Diervilla lonicera (11.39%), Vaccinium angustifolium (11.11%), and Cornus canadensis Linnaeus (Cornaceae) (10.84%). Nine floristic species were significant indicators (IndVal index, P-value < 0.05) of the central area: Viburnum cassinoides Linnaeus (Viburnaceae), Cornus canadensis, Linnaea borealis Linnaeus (Caprifoliaceae), Diervilla lonicera, Clintonia borealis (Aiton) Rafinesque (Liliaceae), Ilex mucronata (Linnaeus) M. Powell, V. Savolainen and S. Andrews (Aquifoliaceae), Aralia hispida Ventenat (Araliaceae), Epigaea repens Linnaeus (Ericaceae), and Viola spp. Linnaeus (Violaceae).

In the northern area, 20 floristic plant species were recorded. Floristic composition of sites within the area was dominated by Rhododendron groenlandicum (Oeder) Kron and Judd (Ericaceae) (23.83%), Kalmia angustifolia (23.74%), and Vaccinium angustifolium (23.74%). Three species were significant indicators (IndVal index, P-value < 0.05) of the northern area: Rhododendron groenlandicum, Gaultheria hispidula (Linnaeus), Muhlenberg ex Bigelow (Ericaceae), and Chamaedaphne calyculata (Linnaeus) Moench (Ericaceae).

Anthophila communities

Using the redundancy analysis, we found that only three of the 51 floristic species sampled along our gradient significantly described the relative abundance of Anthophila genera: Diervilla lonicera (DIE, F _1,29 = 13.649, P-value < 0.001), Rhododendron groenlandicum (LEG, F _1,29 = 5.722, P-value = 0.001), and Vaccinium angustifolium (VAA, F _1,29 = 2.178, P-value = 0.076). Only the first axis of the redundancy analysis was significant (Axis 1: F _1,29 = 19.202, P-value < 0.001). Our model explained 42.6% of Anthophila genera relative abundance, with an adjusted R ²-value of 0.36, leaving approximately 57.4% unexplained. The genera omitted from the analysis due to their scarcity in the dataset (n < 10) were Augochlora, Coelioxys Latreille, Dufourea, Epeolus, Hoplitis (Klug), Nomada, Pseudopanurgus, and Stelis spp.

As presented in Figure 5A, the genera Lasioglossum and Bombus spp. tend to be associated with different floral communities. The relative abundance of Lasioglossum spp. seems related to the presence of Diervilla lonicera. The genera Bombus and Sphecodes spp. are positively correlated with sites from the northern area, and their relative abundance appears, respectively, correlated with Vaccinium angustifolium and Rhododendron groenlandicum.

Figure 5. Redundancy analysis (RDA) triplots of Hellinger-transformed bee data explained by A, three significant plant species for Anthophila, and B, eight significant plant species for bumble bees. Anthophila genera code (panel A) appears in black capital letters: A, Andrena spp.; B, Anthophora spp.; C, Augochlorella spp.; D, Bombus spp.; E, Ceratina spp.; F, Colletes spp.; G, Halictus spp.; H, Hylaeus spp.; I, Lasioglossum spp.; J, Megachile spp.; K, Melissodes spp.; L, Osmia spp.; and M, Sphecodes spp. The floristic species code appears in red capital letters: ARN, Aralia nudicaulis; ASA, Oclemena acuminate; ASM, Eurybia macrophylla; ERR, Acer rubrum; DIE, Diervilla lonicera; KAA, Kalmia angustifolia; LEG, Rhododendron groenlandicum; LIB, Linnaea borealis; MAC, Maianthemum canadense; PRP, Prunus pensylvanica; and VAA, Vaccinium angustifolium. Site icon colour and shape reflect the three areas.

Bumble bee communities

Using redundancy analysis, we found that eight of the 51 floristic species sampled along the gradient significantly described the relative abundance of bumble bee species: Maianthemum canadense (F _1,22 = 9.216, P-value < 0.001), Eurybia macrophylla (F _1,22 = 6.450, P-value < 0.001), Prunus pensylvanica (F _1,22 = 5.761, P-value < 0.001), Oclemena acuminata (Michaux) Greene (Asteraceae) (F _1,22 = 4.902, P-value < 0.001), Kalmia angustifolia (F _1,22 = 3.481, P-value = 0.009), Aralia nudicaulis Linnaeus (Araliaceae) (F _1,22 = 3.289, P-value = 0.013), Linnaea borealis (F _1,22 = 2.475, P-value = 0.042), and Acer rubrum Linnaeus (Sapindaceae) (F _1,22 = 2.382, P-value = 0.053). The first three axes of the redundancy analysis were significant (Axis 1: F _1,24 = 17.839, P-value < 0.001; Axis 2: F _1,24 = 11.485, P-value < 0.001; Axis 3: F _1,24 = 6.062, P-value = 0.001). Our model explained 61.2% of bumble bee species relative abundance, with an adjusted R ²-value of 0.48, leaving approximately 38.7% unexplained. The bumble bee species omitted from the analysis due to their scarcity in the dataset (n < 10) were B. rufocinctus and B. insularis Smith.

As presented in Fig. 5B, a strong negative correlation is suggested between B. terricola and B. ternarius. They are, respectively, linked with Prunus pensylvanicus and Kalmia angustifolia. A strong negative correlation is also suggested between B. sandersoni and B. vagans. Bombus sandersoni, B. frigidus, and B. flavidus are positively correlated with sites from the northern area. Also, because it is the bumble bee species correlated with the most floral species, including two significant indicators of the southern area (i.e., M. canadense, E. macrophylla, and A. nudicaulis), B. vagans seems positively correlated with floristic species from temperate forests.

Discussion

To our knowledge, our study is the first to assess the influence of a regional bioclimatic gradient on wild bee communities in managed boreal forests. We found that climate and floristic resources are key factors driving community assemblages for both Anthophila and bumble bees along this spatial gradient. Our hypothesis that wild bee diversity would be inversely related to latitude was partially supported at the genus level for Anthophila communities and at the species level for bumble bee communities. Communities at both taxonomic levels were unbalanced along the gradient, which could be explained by the dominance of a small number of competitive taxa. We were able to identify taxa that were exclusive to, or representative of, specific areas. Finally, several floristic species were found to play a determinant role in the relative abundance of some taxa within wild bee communities.

Anthophila communities

Our hypothesis was partially supported in the assessment of Anthophila genera. Abundance and richness were both inversely related to latitude, but Shannon’s diversity index and Pielou’s evenness index were not. The north–south latitudinal gradient in forest characteristics in the province of Quebec is driven by climatic conditions that become colder towards the north, with a corresponding decrease in forest productivity (Berteaux et al. 2014). Indeed, in the study areas, post-clearcut vegetation development tended to be more advanced in the southern part of the gradient (L.C.G., unpublished data). In most study sites oin the southern area, a young tree layer had already established, which often overshadowed the understorey herbaceous and shrub layers. In contrast, in the northern area, the vegetation was low and consisted mainly of shrubs and herbaceous plants. The arborescent layer was not yet well developed because of the slow growth of the coniferous species typical of this climate. These observations in clearcuts indicate that the growth speed of regeneration and time-lapse between ecological successions depend on the latitudinal situation along the spatial gradient (De Frenne et al. 2012). It is possible that wild bee communities observed in clearcuts in the south of our gradient are typical of more advanced forest succession than those observed in clearcuts to the north.

The floristic composition of the sites of the northern area had up to 1.8 times fewer floristic species than those in the other two, more southerly areas. Moreover, it was largely dominated by three Ericaceous shrubs: Rhododendron groenlandicum, Kalmia angustifolia, and Vaccinium angustifolium. The more diverse floral communities in the southern and central areas of our gradient may have played a determinant role in Anthophila genera richness and abundance. As shown previously, a sufficient variety of floral species must be present to support diverse and abundant wild bee communities (Venjakob et al. 2016). Floristic diversity is known to be critical to ensuring overlapping blooms throughout the entire season (Fortuin and Gandhi 2021). Anthophila taxa have variable seasonal activity periods, from early season to full season, and sometimes emerge and forage within only a short 4- to 6-week window (Packer et al. 2007). Bloom patterns are therefore one of the environmental factors limiting Anthophila taxa within open forest habitats (Fortuin and Gandhi 2021).

The Pielou’s evenness index results indicate that, in the southern and central areas, Anthophila communities are composed of a high number of rare genera but are dominated by one genus – Lasioglossum spp. In the northern area, Anthophila communities are also composed of a great number of rare genera, but the genus Lasioglossum spp. is codominant with the genus Bombus. Lasioglossum spp. was by far the most abundant genus captured in clearcuts in our study. This genus is very common throughout temperate and boreal forests (Packer et al. 2007) and was also the most abundant genus caught in other wild bee studies in clearcuts (Romey et al. 2007; Gibbs 2010; Proctor et al. 2012; Fortuin and Gandhi 2021). The genus Lasioglossum has been described as indicative of recently disturbed landscapes and open forested areas (Proctor et al. 2012), mainly because most Lasioglossum species are ground nesters (Wilson and Carril 2015). Proctor et al. (2012), Rubene et al. (2015), and Fortuin and Gandhi (2021) observed that regenerating clearcuts tend to favour groundnesting bee groups because the percent cover of warm bare soil is consistently higher there than what is found in undisturbed forests.

It is noteworthy that results of the Anthophila analyses would have been different if performed at the species level. For example, the lower evenness levels observed in the southern and central areas are probably explained in part by the disproportionate abundance of Lasioglossum spp. in those areas. However, because this is known to be a speciose genus (Gibbs 2010), the evenness index results could have changed dramatically and even flipped (i.e., they could have been higher in the southern area) if the analyses were made at the species level.

Bumble bee communities

Regarding the bumble bee communities, our hypothesis was only partially supported at the species level for the genus Bombus. Neither richness nor Shannon’s diversity index decreased from south to north, revealing no significant relationship with latitude. The results of Pielou’s evenness index, which also did not differ significantly between areas, showed that bumble bee communities along the entire gradient are generally composed of a few competitive species that dominate a greater number of rare species. Only abundance decreased from the southern to the northern area.

Many factors could explain the absence of the expected decreasing relationship of richness to latitude. For example, bumble bees are known to include species that are particularly well adapted to colder climates (Williams et al. 2014), and some have even been seen foraging at temperatures close to the freezing point (Goulson 2010). Their large bodies are covered by thick and dense coats of hairs. They can shiver their flight muscles, which allows them to produce their own body heat (Heinrich and Esch 1994). The flight muscles can disconnect from the wings, thereby saving energy during shivering (Whitfield et al. 2013; Wilson and Carril 2015).

Despite these genus-level adaptations, individual bumble bee species tend to be limited to relatively narrow climatic niches (Kerr et al. 2015). This may explain how, instead of finding decreasing bumble bee diversity as latitude increased, we observed relatively distinct species assemblages in each bioclimatic domain. Different bumble bee species are known to be differentially distributed in eastern Canada (Laverty and Harder 1988).

All bumble bee species that we observed were already recognised as part of the boreal fauna, but we lacked a clear understanding of how these boreal species sort themselves within the boreal habitat. Our study reveals specific niche requirements for bumble bees within this important zone, which is one of the most vulnerable to climate change. The distribution of bumble bee abundances in our captures suggests that B. sandersoni, B. terricola, and B. ternarius may possess larger climatic niches than the other bumble bee species within our study areas: those three species were well distributed along the entire gradient. Only a few observations of B. sandersoni have been documented to date in the province of Quebec (Bumble Bee Watch 2022; Global Biodiversity Information Facility 2022); our results provide substantial new information on the distribution of this species and its potential response to climate, B. sandersoni being the third most abundant species among our bumble bee captures.

Bombus terricola was the second most abundant species among our bumble bee specimens. This was unexpected because the species has been reported to be declining in parts of its range in North America (Grixti et al. 2009; Cameron et al. 2011; Bartomeus et al. 2013). NatureServe (2022) and the IUCN Red List (Hatfield et al. 2015) register B. terricola as “vulnerable.” The status of B. terricola in the province of Quebec is undetermined because a large part of the species’ distribution range remains sparsely surveyed. Pathogen spillover from domestic bumble bees used in greenhouses and pesticides such as neonicotinoids are correlated with the species’ decline throughout its North American range and appear to be some of the most significant threats to this species (Hopwood et al. 2012; Szabo et al. 2012). However, the remote environment examined in the present study may still be exempt from such influences. The important presence of B. terricola in clearcuts in our study also implies that this kind of altered environment could support the species, at least in the short term, by providing an adequate habitat that is free of those threats. Competition with other bumble bees has also been shown to impact this vulnerable species: Butler et al. (2021) observed that the array of floral species visited by B. terricola in Michigan, United States of America, consistently narrowed when B. ternarius was also present. As a well-known generalist, B. ternarius has a wide dietary breadth, which includes most of the floral species favoured by B. terricola (Butler et al. 2021). In our study, the abundance of B. terricola was greater than that of B. ternarius on most sites in the southern area, indicating that dense regeneration in clearcuts might attenuate the competitive impact of B. ternarius. Furthermore, B. ternarius is a common bumble bee species in the northeastern boreal forest and is known to inhabit many open woodland and wetland habitats within natural and disturbed landscapes (Williams et al. 2014). Our results suggest that although B. ternarius has great versatility in terms of choice of habitat and floral resources, B. terricola can outcompete it in clearcuts with a dense arborescent floral composition.

Effect of sampling method

Concerns about potential biases from using only the pan-trapping method for bee collection are raised throughout the literature (e.g., Popic et al. 2013; Fortuin and Gandhi 2021). Captures of bumble bees and other large-bodied bee groups are limited when using pan traps (Romey et al. 2007), and thus, the dominance of small bees within our samples (∼75% of captures) may be related to this bias. A combination of pan trapping and netting is suggested to obtain samples with the greatest diversity of bees in flower-rich areas (Tucker and Rehan 2017). However, hand netting was impractical in most of our settings, considering the density and height of the vegetation on multiple sites. Also, unpredictable weather during our tight-scheduled surveys shortened the time frame during which hand netting was feasible because this method cannot be performed when the vegetation is wet (Popic et al. 2013). Another inadvertent sampling bias is mentioned in the literature. Because bees use visual colour cues to navigate to flowers, some studies warn that traps visible from a greater distance are more likely to catch bees than traps that have visual obstructions (Song et al. 2015). Bees might also be underrepresented in traps that are situated in understorey vegetation when foraging resources are available in the canopy (Urban-Mead et al. 2021). Even though visual obstruction and canopy resources did not seem to affect our catch rates, caution suggests accounting for possible bias, and therefore, the wild bee diversity metrics documented here should not be considered absolute.

Conclusion

For the two taxonomic levels investigated, our study demonstrates that latitudinal dissimilarities in climate and floristic composition are determinant factors of wild bee community assemblages within boreal forest clearcuts. In general, wild bee communities were unbalanced in that they tended to be dominated by a few taxa. The results of the present study also show that disturbed forests are a critical habitat for bees and that this information could be used in forest planning, in particular by using interventions that might affect floral diversity.

As revealed by the results, abundance and richness of Anthophila genera and abundance of bumble bee species were inversely related to latitude in clearcuts along a bioclimatic gradient in the province of Quebec. However, this relationship was not observed for the richness of bumble bee communities. Divergent results between the two taxonomic levels suggest an underlying complexity that may be masked by the predominance of a sole taxon. Bees of each genus will need to be analysed at the species level because wild bee taxa are highly diverse in terms of their life histories, social structures, and nesting biology (Michener 2007; Wilson and Carril 2015). Also, to maintain pollination services and advise policymakers regarding the migration capabilities of these taxa in the context of climate change, further research is needed to deepen basic knowledge on wild bee dispersion ranges within the boreal forest. To our knowledge, the present study provides the first assessment of wild bees in clearcuts in the province of Quebec and could be used as a baseline for the study of environmental changes on these organisms.