Antimicrobial resistance (AMR) is recognized by the World Health Organization (WHO) as one of the top ten public health threats facing humanity1. Globally, AMR contributes to approximately 4.95 million deaths each year2, and without intervention, this number is projected to exceed 10 million by 20501. The overuse of antibiotics, especially within agricultural sectors (for example, veterinary use of medically important antimicrobials accounts for ~80% of total antibiotics sold in the USA3), is considered a major driver of AMR and is directly associated with increased prevalence of drug-resistant infections in both human and livestock populations4,5,6,7. In 2015, the WHO in collaboration with the Food and Agriculture Organization and World Organization for Animal Health developed a global action plan for addressing AMR8. One major advancement was the introduction of the WHO’s AWaRe (Access, Watch, Reserve) classification system, which categorizes antibiotics into three tiered groups to facilitate their prudent use and minimize future emergence of AMR. This framework has since been adopted by many federal agencies, including the US Food and Drug Administration in 20179 and Public Health Agency of Canada in 201810, eliminating over-the-counter access to livestock antibiotics and mandating veterinary supervision via prescription or Veterinary Feed Directive issued within the context of a veterinarian–client–patient relationship11. The United Nations further released a political declaration12, signed by 193 countries, committing the world to take collective action through a ‘One Health’ approach by reducing medically unnecessary use of antimicrobials in humans and in animals13. Such stringent regulatory measures have greatly restricted the use of antibiotics for growth promotion and feed efficiency in livestock animals.

Honeybees epitomize One Health initiatives through their inseparable symbiosis with environmental health determinants14 while also presenting unique challenges in the context of AMR15. In the beekeeping industry, oxytetracycline is the most common antibiotic used for the prevention of infectious diseases such as American foulbrood (AFB) and European foulbrood (EFB) caused by Paenibacillus larvae and Melissococcus plutonius, respectively16. There have been widespread reports of tetracycline-based resistance in both pathogens since the early 2000s17, and in the case of oxytetracycline-resistant AFB, lincomycin and tylosin represent the main alternatives. Concerningly, recent reports suggest that there is also emerging resistance to both of these antibiotics in North America18,19. Fumagillin is another common antimicrobial used in beekeeping for the prevention of a serious diarrhoeal disease known as nosemosis, caused by two microsporidia, Vairimorpha (formerly Nosema) apis and V. ceranae—the latter being more virulent and prevalent worldwide20,21. A recent systematic review of 50 controlled field trials found no evidence to indicate that resistance has developed in either species thus far22. However, there may be inherent differences between the species in terms of baseline susceptibility to fumagillin23, suggesting that labelled dosages should be revised to appropriately reflect treatment thresholds and to achieve proper usage, maximum efficacy and safe administration in alignment with antibiotic stewardship goals24.

Beyond the issue of AMR, the impact of antibiotics on the microbiota of honeybees is another area of concern. Broad-spectrum antibiotics indiscriminately affect both harmful as well as beneficial microorganisms, which can disrupt gut microbiota homeostasis and lead to a scenario of ‘missing microbes’25. This disruption compromises immune function (for example, ~30% reduction in haemolymph killing capacity after 2-week oxytetracycline exposure26) and can paradoxically increase honeybee mortality by causing increased susceptibility towards other types of opportunistic infection by bacterial27, fungal28 and viral29 pathogens. Another overlooked area with regard to AMR is exposure to agrochemicals as well as other hive medications (for example, miticides for the treatment of ectoparasitic Varroa destructor), which, despite being registered as ‘non-antimicrobial’ products, can nonetheless have potent off-target deleterious effects on both honey-bee-associated symbionts and pathogens30. Importantly, while the short-term impacts of antibiotics on individual worker bee health have been well documented, the long-term impacts on overall colony-level health remain vastly unknown31.

Here, we analyse national survey data from the Canadian Association of Professional Apiculturists (CAPA) with the aim of gaining novel insight on longitudinal antibiotic usage trends and characterizing the associations between antibiotic exposures and overwintering mortality rates. This study is timely given recent antibiotic regulatory changes and, in particular, captures the transition from antibiotics used in the beekeeping industry being accessible over the counter to requiring a prescription or Veterinary Feed Directive. Moreover, the CAPA survey is unique among large-scale surveys (for example, the USDA–APHIS Bee Informed Partnership32 in the USA and COLOSS in Europe33) in that it not only collects annual colony loss data but also tracks annual usage of all relevant hive medications at the province scale, including miticides as well as antimicrobials such as oxytetracycline, tylosin, lincomycin and fumagillin. In addition, this study investigates the potential influence of air pollution and climate factors on honeybee mortality, recognizing that environmental stressors may interact with disease dynamics and antimicrobial use to impact colony health. By examining these data, we can better understand the real-world impacts of regulatory changes on antibiotic use in beekeeping and their potential effects on colony survival. Ultimately, this study offers an evaluation of how policy shifts are influencing antibiotic practices and health outcomes in apiculture at a national scale; this insight therefore offers much needed guidance on future policy development in the agricultural and environmental health sectors.

Overview of surveyed beekeeper population from 2015 to 2023

The number of managed honeybee colonies across Canada ranged from 714,144 to 817,926 during 2015–2023. A total of 443–626 beekeepers responded to the survey each year (Supplementary Data 1). At the province scale, there were significant differences in the proportion of respondents (one-way analysis of variance (ANOVA), F9,78 = 9.19, P = 2.34 × 10−9), with Prince Edward Island (90.7 ± 3.0), New Brunswick (75.1 ± 4.5) and Newfoundland and Labrador (73.8 ± 7.5) showing the highest participation (mean ± standard error of the mean) and Saskatchewan (38.6 ± 3.0), Manitoba (46.4 ± 2.6) and Ontario (49.8 ± 4.5) showing the lowest participation (Fig. 1a). Collectively, the respondents managed between 362,949 and 502,764 colonies annually, accounting for 48.0–63.9% of Canada’s total number of registered bee colonies covered by the survey each year (Fig. 1b).

a, A visual representation of the proportion of survey respondents by province from 2015 to 2023, highlighting regions with the highest and lowest participation. The map was generated via the sf (v1.0.15) package in R, with spatial data derived from GADM (v3.6) and the coordinate reference system transformed to ‘EPSG: 3348 – NAD83 (CSRS) / Statistics Canada Lambert’. b, The annual number of managed honeybee colonies operated by survey respondents compared with the total number of colonies across Canada. c, Top factors contributing to overwintering colony mortality as reported by beekeepers. The circle plot generated via the packcircles (v0.3.6) package in R. The radius of circles was scaled to the RCF index, with larger circles representing factors more commonly reported than smaller circles. d, A ridgeline plot highlighting the distribution of overwintering colony losses for each province. The plot was generated via the ggridges (v0.5.6) package in R. e, A line plot highlighting overwintering colony losses for each province over time. The total line (red) represents the weighted estimate for all of Canada. AB, Alberta; BC, British Columbia; MB, Manitoba; NB, New Brunswick; NL, Newfoundland and Labrador; NS, Nova Scotia; ON, Ontario; PE, Prince Edward Island; QC, Quebec; SK, Saskatchewan.

Full size image

In terms of reported mortality rates, nationwide overwintering colony losses (mean ± standard error of the mean) were determined to be 27.6 ± 3.0% with a significant upward trend apparent during the 2015–2023 period (linear regression, F9,78 = 7.86, P = 2.6 × 10−2). There were significant differences in overwintering losses between provinces (one-way ANOVA, F9,78 = 2.75, P = 8.0 × 10−3), with Prince Edward Island (34.5 ± 5.2%), Alberta (31.4 ± 4.2%) and Ontario (30.3 ± 4.0%) showing the greatest overall losses. Uniquely, Nova Scotia demonstrated the lowest overall loss (15.9 ± 0.9%), which was significantly lower than the national mean (Fig. 1d,e). The top three contributing factors to overwintering mortality as cited by beekeepers were determined to be ‘Poor queens’ (ranked causal factor (RCF) 0.69), ‘Weather’ (RCF 0.46) and ‘Varroa’ (RCF 0.41; Fig. 1c).

Trends in usage of antibiotics and other hive medications

To assess the impact of recent regulatory changes on antibiotic usage in the beekeeping industry, we compared the cumulative usage of oxytetracycline during years before (2015–2018) and after (2019–2023) the new policies took effect (Fig. 2a–c). There was a significant decline in the median (± median absolute deviation) usage of oxytetracycline at the national level (46.9 ± 2.2% before versus 24.1 ± 5.0% after, Wilcoxon rank-sum test, W = 20.0, P = 1.6 × 10−2; Fig. 2d). Similar trends were also observed at the province level with few exceptions (Fig. 2e). Newfoundland and Labrador and Quebec reported low oxytetracycline use (<25%) both before and after regulatory changes, which may be explained by differences in recommended best management practices for AFB and EFB in these provinces34. Among provinces in which the majority of beekeepers (>50%) reported using antibiotics during 2015–2018, Ontario was the only province that did not show a significant decline (W = 17, P = 1.1 × 10−1) in antibiotic usage following regulatory changes.

a–c, Line plots depict the relative proportion of beekeepers who reported using oxytetracycline (a), fumagillin (b) or miticides (c) in each province over time. The total line (bolded) for each represents the weighted estimate for all of Canada. Vertical dashed red lines highlight the time at which Public Health Agency of Canada (PHAC) antibiotic regulatory changes took effect. d–i, Comparisons for oxytetracycline (d and e), fumagillin (f and g) and miticide use (h and i) before (2015–2018) and after (2019–2023) PHAC regulatory changes at the national (d, f and h) and province (e, g and i) levels. Data points represent the annual usage reported for each province (e, g and i) or the annual weighted estimate for all of Canada (d, f and h), while box plots depict the respective median (line in box), interquartile range (box) and minimum and maximum (whiskers) values, with comparisons shown for two-sided Wilcoxon rank-sum tests. All plots were generated in R using the ggplot2 (v3.5.0) package. Statistics determined via the ‘wilcox_test’ function of the rstatix (v0.7.2) package and annotated with the ‘stat_pvalue_manual’ function of the ggpubr (v0.6.0) package. AB, Alberta; BC, British Columbia; MB, Manitoba; NB, New Brunswick; NL, Newfoundland and Labrador; NS, Nova Scotia; ON, Ontario; PE, Prince Edward Island; QC, Quebec; SK, Saskatchewan; NA, not available.

Full size image

Fumagillin (bicyclohexylammonium fumagillin) is an antifungal agent and the only registered form of treatment for nosemosis disease in honeybees. Although this antimicrobial falls within a class of drugs not requiring a veterinarian prescription in most areas of Canada, closure of the sole supplier company (Medivet Pharmaceuticals) that coincided with the 2018 regulatory changes resulted in widespread shortages of fumagillin in Canada35. Reported median usage of fumagillin during 2015–2018 (47.5 ± 3.8%) was found to be significantly higher compared with the 2019–2023 period (19.1 ± 1.3%) across Canada (Wilcoxon rank-sum test, W = 20, P = 1.59 × 10−2; Fig. 2f). Similar trends were also observed at the province level albeit with substantial variation in overall usage rates between provinces (Fig. 2g). Alberta and Saskatchewan demonstrated the greatest decline in usage from ~80% and ~50%, respectively, during the 2015–2018 period to less than 30% by 2020 (Fig. 2f). Corresponding with the reintroduction of fumagillin to Canadian markets in 2021 (with approval of Drug Identification Number 02231180, allowing production by Vita Bee Health36), there was a slight rebound (~15% increase), especially in the prairie provinces (Alberta, Saskatchewan and Manitoba), whereas many eastern provinces exhibited persistent declines in usage of fumagillin—with Ontario, Quebec and Prince Edward Island reaching 5% or less reported usage by 2023 (Fig. 2f).

As a relative comparison with non-antimicrobial hive medications used in beekeeping, we assess usage trends of products registered for the management of ectoparasitic Varroa destructor mites, including amitraz and fluvalinate miticide products as well as formic acid and oxalic acid fumigation products (Fig. 2c). We found no observable differences in the reported usage of these products over time with consideration to the pre- and postregulatory changes relevant to antibiotic accessibility during periods of 2015–2018 compared with 2019–2023 (68.8 ± 6.6% versus 74.2 ± 6.4%, respectively, Wilcoxon rank-sum test, W = 4, P = 0.190). Similar trends were seen at the province level, with little variation and >75% of beekeepers in most provinces reporting the usage of at least one product. One exception was Newfound and Labrador, in which none of the products was reported to be used as V. destructor mites are not present in this area of Canada.

Linking antibiotic use and environment with mortality

While diseases such as AFB and nosemosis are known risk factors that negatively impact colony survival37,38, recent evidence also suggests that off-target side effects of antibiotic exposure can have deleterious long-term health effects in honeybees26,27. Here, we examine the relationship between antibiotic usage trends and overwintering mortality based on CAPA survey data (Fig. 3). Oxytetracycline use at the national level showed a significant downward trend over time during the 2015–2023 period (linear regression, F1,7 = 47.82, P = 2.3 × 10−4) and a strong negative trending association with overwintering mortality (Pearson correlation coefficient, r = −0.62, P = 0.0730; Fig. 3a,d). Similarly, fumagillin usage also showed a significant downward trend over the 2015–2023 period (linear regression, F1,7 = 31.21, P = 8.3 × 10−4) and a strong negative association with overwintering mortality (Pearson correlation coefficient, r = −0.67, P = 0.049; Fig. 3b,e). Finally, as a point of comparison with non-antibiotic hive medications, trends in mite management were evaluated for the cumulative usage amitraz-, fluvalinate-, formic acid- and oxalic acid-based products. There was an increase in the usage of these miticides during the 2015–2023 period (linear regression, F1,7 = 6.40, P = 3.9 × 10−2), although the positive association with overwintering mortality did not reach significance (Pearson correlation coefficient, r = 0.58, P = 9.8 × 10−2; Fig. 3e,f).

a, Oxytetracycline usage over time and its negative association with overwintering mortality. b, Fumagillin usage over time and its negative association with overwintering mortality. c, Trends in mite management product usage showing no significant changes over time. Temporal overlay plots were generated in R using the ggplot2 (v3.5.0) package, with mean trendlines calculated using local polynomial regression fitting (geom_smooth function, method = ‘loess’) and shaded areas indicating 95% CIs. d–f, Scatter plots show the relationship between overwinter mortality and the usage of oxytetracycline (d), fumagillin (e) or miticides (f). Two-tailed Pearson correlation statistics shown for each comparison were determined with the ‘sm_statCorr’ function of the smplot2 (v0.2.1) package. g, Mixed linear model results highlighting significant predictors of overwintering mortality for model 6, including antibiotic use and air pollution factors. Data points depict the structure coefficients, and lines depict the 95% CIs. The plot was generated via the ‘forestplot’ function of the partR2 (v0.9.2) package. h, A bar plot showing the relative variance explained by predictors in model 6. Variance statistics determined from partial R2 values using the ‘partR2’ function of the partR2 package (R2_type = ‘conditional’ nboot = 100) and plotted with the ggplot2 package. TmWinter, winter temperature; TmSpring, spring temperature; TmSummer, summer temperature; TmFall, fall temperature; P, precipitation; S, snowfall.

Full size image

In addition to simple regression analysis, we conducted multivariate analysis to account for the potential effects of environmental confounders, such as climate variability39 and air pollution40, on observed overwintering mortality. A series of mixed linear models were tested incorporating multilevel fixed effects (antibiotic use, miticide use, temperature, precipitation, snowfall, nitric oxide (NO), nitrogen dioxide (NO2), sulfur dioxide (SO2), fine particulate matter (<2.5 μm) (PM2.5), carbon monoxide (CO) and ozone (O3)) and random effects (province) to assess the impact of environmental covariates and control for baseline geographical differences, respectively (Table 1). Antibiotic use (combined estimate of oxytetracycline, tylosin, lincomycin and fumagillin) consistently emerged as a significant predictor of overwintering survival across all models evaluated. The inclusion of random effects in model 2 improved the fit of coefficients compared with model 1 (−270.123 versus −272.451 log-likelihood, respectively), whereas in model 5, including all fixed- and random-effect terms demonstrated overall better performance (log-likelihood −206.945, Akaike information criterion corrected (AICc) 419.434) with a low random effect variance (σ2 = 51.79) indicating optimal resolution of between-province heterogeneity and effective capture of the underlying data structure. Inclusion of temperature, precipitation and snowfall as predictors in model 3 resulted in all three being identified as significant. However, subsequent inclusion of air pollution predictors in model 4 revealed that NO, NO2, SO2 and CO had much greater effects on predicted mortality, which remained significant after accounting for co-linearity of all predictors together in model 5 (Fig. 3g). Given the known importance of seasonality and temperature shifts in determining successful active foraging periods for honeybees41, we further assessed the influence of meteorological predictors specifically based on their estimated values during winter, spring, summer and fall seasons. Inclusion of seasonal temperature predictors in model 6 increased the explainable variance and overall fit (conditional R2 = 0.719, AICc 388.759), whereas negligible changes were seen for inclusion of seasonal precipitation (conditional R2 = 0.457, AICc 418.725) and snowfall (conditional R2 = 0.442, AICc 421.225).

Among the 71.9% of total variance explained by model 6, antibiotic use (partial R2 = 0.126, P ≤ 0.0001), NO2 air pollution (partial R2 = 0.175, P = 0.0031), winter temperature (partial R2 = 0.154, P ≤ 0.0001) and snowfall (partial R2 = 0.124, P = 0.0021) emerged as the top contributors and together explained the majority (57.9%) of total variance (Fig. 3h). Based on the fixed part of model 6, structure coefficients (the contribution of a predictor to model prediction irrespective of other predictors) indicated antibiotic use (95% confidence interval (CI) −0.436 to −0.118) as having the strongest negative correlation with mortality, whereas NO2 (95% CI 0.242 to 0.483) and CO (95% CI 0.194 to 0.444) air pollutants showed the strongest positive correlations (Fig. 3g). After adjusting for the variance uniquely explained by each fixed effect, standardized beta weights (standardized estimates of effect size42) indicated that antibiotic use and winter temperature had the greatest magnitude of negative effects on mortality (that is, these factors increased survivorship), while air pollutants generally showed the greatest magnitude of positive effects. Given the vast number of possible model combinations, we also used MuMIn43 to systematically evaluate all 1,048,576 models, which supported the best-fitting model based on AICc weights while considering biological relevance (Supplementary Data 2). Overall, antibiotic use, NO2 and winter temperature stood out as having notably better predictive performance—regardless of the direction of effect—compared with the remainder of variables assessed (Table 1).

This study demonstrates that antibiotic use in the Canadian beekeeping sector has drastically decreased in recent years and that this change in usage tightly coincides with the rising rates of overwintering mortality seen nationwide over the same period. Our analysis also reveals that several air pollutants are significant predictors of overwintering mortality in honeybee colonies—notably nitrogen dioxide. Altogether, the findings provide new insights into the complicated relationships between antibiotic regulation, best management practices and honeybee health outcomes, emphasizing the importance of balanced approaches in mitigating AMR while safeguarding pollinator health. The multivariate and sometimes non-additive nature of these relationships suggest that no single-factor change to regulations or local best-management practice will provide a panacea to beekeeper challenges. While reduced antibiotic usage was associated with increased colony losses, this relationship may not reflect a simple cause–effect dynamic. Years of frequent antibiotic use may have altered honeybee microbiota and disease susceptibility, potentially creating a dependency on antimicrobials for maintaining colony health in the absence of ecological resilience. We therefore advocate for systemic approaches to bee health that are in line with the broader One Health philosophy.

Beyond direct antibiotic use in beekeeping, honeybees are also exposed to antimicrobial residues and AMR genes from environmental sources. Their extensive foraging behaviour makes them effective bioindicators for tracking environmental AMR within the One Health framework44,45. Residues of various antibiotic classes—including tetracyclines, macrolides, aminoglycosides and sulfonamides—have been detected in honeybee products, often originating from environmental contamination46. Furthermore, multiple studies have shown that AMR genes are widespread in honeybee colonies, with some associated with human pathogens such as Salmonella enterica, Escherichia coli and Staphylococcus aureus47. This highlights the potential for honeybees to act as vectors of AMR between agricultural landscapes and urban environments. Given the growing evidence that climate change exacerbates pathogen burden48 and AMR gene prevalence49, future efforts should take a holistic approach to honeybee health50, considering both direct antibiotic use and broader environmental influences on resistance dynamics.

Our analysis showed that national usage of antibiotics such as oxytetracycline, tylosin, lincomycin and fumagillin was effectively halved, decreasing from approximately 50% to 25% following the regulatory changes implemented in 2018 (Figs. 2 and 3). This indicates a successful outcome with regard to the goal of minimizing indiscriminate antibiotic use in agriculture and also mirrors the findings from a recent report by the World Organisation for Animal Health (based on data collected from 72 countries) demonstrating that overall antibiotic use in animals decreased by 27% in response to global regulatory changes between 2016 and 201851. Notably, British Columbia and Quebec consistently reported low usage of bacterial antibiotics (<25%) both before and after regulatory changes (Fig. 2d,e), which is probably explained by the provincial recommendations of metaphylactic as opposed to prophylactic use of antibiotics to control AFB and EFB in these regions34. By contrast, Ontario demonstrated the highest overall usage of bacterial antibiotics and reported only a marginal decrease (~10%) in their usage, which may be related to increased disease pressure in this region (for example, advisory issued for EFB outbreaks in 202352) as well as the enhanced antibiotic access resources made available through collaboration of the provincial government (Ontario Ministry of Agriculture, Food and Rural Affairs) and Ontario Beekeepers’ Association53,54.

The upward trend in overwintering mortality rates alongside reduced antibiotic usage presents an intriguing dilemma (Fig. 3). While the regulatory intent is to mitigate the risk of AMR, the data suggest that other factors may be influencing colony health and survival rates. For instance, anthropogenic air pollutants such as NO2 (commonly emitted in diesel exhaust) can react with floral odours, which are commonly volatile organic compounds, altering their detection by honeybees55, which ultimately reduces foraging efficiency by 83–90% (ref. 56) and, in the case of repeated short-term exposures, can negatively impact overall colony fitness57. Corroborating this, our multivariate analysis demonstrated that NO2 air pollution was the strongest predictor of overwintering mortality among all other environmental factors assessed (Table 1). Timing and fluctuations in atmospheric NO2 levels may also impact overwintering mortality, given that NO2 increases over the summer before peaking in the fall, while dry deposition of NO2 peaks in late summer58,59. These fluctuations are partially driven by the ability of NO2 to act as an oxidizing agent of volatile organic compounds as ultraviolet radiation catalyses these reactions60 and ultraviolet levels are highest during July and August in North America. The timing of high atmospheric NO2 and its dry deposition coincides with periods when honeybees decrease brood rearing and rely more on stored food in the hive61; therefore, any interference with the honeybees’ ability to find flowers and, conversely, nectar and pollen could negatively impact their ability to build up enough food stores to survive overwintering62. In addition, NO2 has potential to disrupt long-term memory through interfering with glycerophospholipid metabolism63,64, although mechanisms of the involved pathways require further characterization in insects. Nonetheless, in bumblebees (Bombus terrestris), glycerophospholipid metabolism is promoted by symbiotic Lactobacillus apis and Gilliamella spp., which is known to influence memory retention65,66. One study also demonstrated that honeybees treated with tetracycline had poorer olfactory memory and significantly decreased glycerophospholipid pathway activity compared with their conventional counterparts65. Together, these findings suggest that antibiotic usage and exposure to NO2 may be synergistically reducing foraging fitness in bees.

Carbon monoxide was also consistently associated with overwintering mortality across all models tested (Table 1 and Supplementary Data 2). Interestingly, carbon monoxide is known to alter detoxification mechanisms such as cytochrome P450 enzyme-mediated metabolism of insecticides67, with a recent study from the Middle East showing that carbon monoxide is negatively associated with the abundance and diversity of a broad range of insect species68. This suggests that co-exposure to certain air pollutants and pesticides could have synergistic negative effects on honeybees and other insects, warranting further investigation in future studies. Higher snowfall also showed a positive trend with overwintering mortality, altogether depicting a scenario whereby harsher winter conditions combined with air-pollution-mediated impairments to foraging capacity during summer months cumulatively impact the likelihood of honeybee survival. Moreover, recent multicontinental data link poor air quality to increased honeybee stress, immune suppression and pathogen loads, with higher temperatures further exacerbating these effects40,69. Our analysis identified air pollutants and temperature as key predictors of colony mortality, suggesting that the combined impact of air pollution and climate change may progressively weaken honeybee disease resistance. As these stressors intensify in the Anthropocene, managing apicultural diseases will become increasingly complex, particularly given the imperative to limit AMR8. The rising mortality rates observed across Canada probably result from a combination of antibiotic use and broader environmental challenges, many of which remain unmonitored. Notably, 67 pollutants have been modelled to disrupt honeybee olfactory perception of floral cues, potentially impairing foraging efficiency70. These findings underscore the urgent need for comprehensive management strategies that address the interplay of multiple environmental stressors affecting honeybee health.

Nonetheless, antibiotic use remained as one of the top predictors and showed a negative relationship with overwintering colony loss after controlling for confounding environmental factors and other covariates. This potentially suggests that antibiotics are still playing a critical role in managing bacterial and fungal diseases effectively, where alternative strategies are not yet fully compensating for their protective effects. Considering fumagillin as an example, representatives of the Canadian Honey Council estimated in 2019 that beekeepers would expect colony losses to double from around 25% to 50% without a stable supply or an alternative solution to control nosemosis infections71. Here, we show that, while losses remained below 50% in most areas of Canada, nearly all provinces that experienced a major drop in fumagillin use over the past several years also showed a corresponding increase in overwintering mortality (Figs. 2f,g and 3b,e). These trends align with other recent Canadian studies demonstrating that fumagillin can enhance colony survival during the winter months72 and further highlight the context-dependent effects of fumagillin across different geographic regions73.

In terms of other hive medications, a noteworthy observation was the positive trend between miticide use and overwintering mortality (Fig. 3f and Table 1). One conceivable possibility is that the observed miticide use and overwintering mortality were simply responses driven by fluctuating V. destructor infestation levels. A previous study identified V. destructor as a primary driver of colony losses in some Canadian regions74. However, Canadian beekeepers attribute only 41% of losses to V. destructor (Fig. 1c), underscoring the importance of considering environmental covariates when assessing colony health and mortality. The relationship between miticides and mortality is complex, given that miticides can bioaccumulate in wax75 and impair immune signalling during early developmental stages76. It is also plausible that various exposures coupled with declining antibiotic use could cumulatively exacerbate honeybee susceptibility to bacterial, fungal and viral infections. Furthermore, the insidious effects of sublethal pesticide exposure77 as well as unintended side effects of miticides and other insecticides on microbial life30 may play important roles. However, these interactions are vastly underexplored and urgently require further investigation, ideally through leveraging the untapped potential of cell culture models to unravel the complexities underlying honeybee–microorganism–xenobiotic interactions78.

The challenges associated with antibiotic restrictions and disease management in the beekeeping industry clearly emphasize the need for more sustainable therapeutics options that align with One Health initiatives. A probiotic approach, such as the tailored application of beneficial lactic acid bacteria with pathogen-excluding properties, represents a promising solution79,80,81,82, as does selective breeding programmes for disease-resistant bee stocks that aim to bolster intrinsic tolerance to V. destructor mites83. Transgenerational immune priming approaches such as AFB-‘vaccinated’ queen bees (exposed to a bacterin of inactivated P. larvae cells) have been introduced to US and Canadian markets recently, although long-term efficacy remains unclear as a honeybees lack adaptive immunity84,85,86. A sustained commitment to funding development of these and other innovative therapeutic strategies in honeybees is warranted.

There are several limitations to note. Our interpretation is constrained by the structure of the national survey data, which do not distinguish between preventive (prophylactic) and therapeutic (curative) antibiotic use—an important consideration for understanding disease management practices. In addition, factors such as queen quality, mite infestation and starvation—often cited by beekeepers as contributing to mortality—were not included in the multivariate models, as these variables were not quantitatively monitored across years. While reduced antibiotic usage was associated with increased colony losses, this relationship may not reflect a direct causal link. Years of persistent antibiotic use may have deteriorated honeybee microbiota and immune function25, potentially creating a dependency on antimicrobials for maintaining colony health. This counterintuitive trend may suggest that antibiotics, despite their drawbacks, continue to offer short-term protection against bacterial and fungal pathogens in some contexts. However, long-term reliance may have undermined natural defences, making sudden withdrawal more harmful in the absence of effective alternatives. Furthermore, the possibility that antibiotic use and NO2 exposure may interact to reduce foraging success warrants further investigation, as such synergistic effects remain to be fully elucidated.

In conclusion, our study highlights the complex interplay between antibiotic regulation, environmental factors and honeybee health. The significant reduction in antibiotic use following the 2018 regulatory changes, while successful in curbing indiscriminate use and potentially mitigating AMR, coincided with a concerning increase in overwintering mortality rates. This paradox underscores the multifaceted challenges in beekeeping, where environmental pollutants such as NO2 and climatic conditions such as snowfall also play critical roles in colony survival. Our findings emphasize the need for a balanced approach that considers both the prudent use of antibiotics and the broader environmental context to ensure the sustainability and health of honeybee populations. This comprehensive understanding is essential for guiding future policy developments in agricultural and environmental health sectors, ultimately aimed at protecting both pollinators and public health.

Description of survey data

CAPA began surveying beekeepers in 2007 in an effort to monitor overwintering colony loss rates across Canada. Notable exceptions are regions above the 60° north latitude, such as Yukon, Nunavut and the Northwest Territories, which are not included in surveys as environmental conditions are largely unsuitable for honeybees. There is inherent scarcity of data for Newfoundland and Labrador in some years due to the small number of registered colonies present. Statistics for all other provinces including British Columbia, Alberta, Saskatchewan, Manitoba, Ontario, Quebec, New Brunswick, Nova Scotia and Prince Edward Island are compiled at the end of each year and made publicly available in the form of Annual Colony Loss Report summary documents. In 2015, additional criteria were added to the CAPA surveys to track the usage of antibiotics and other hive medicants such as miticides. Consequently, statistics are available for overwintering colony losses spanning from 2007 to 2023 and for trends in antibiotic usage from 2015 to 2023.

Survey data acquisition and processing

Annual survey results (‘FINAL – CAPA Statement on Colony Losses in Canada (YYYY)’) were downloaded in PDF format from the CAPA portal87, and then relevant data (colonies per province, respondent demographics, winter loss, RCFs as cited by beekeepers, and usage of oxytetracycline, tylosin, lincomycin, fumagillin and miticide) from each year were manually extracted and combined into a single wide-formatted TSV file. Notably, the various products used for management of V. destructor mites were not differentiated in the CAPA survey, and thus, miticide usage in this study refers to the percentage of beekeepers using one or more of the following: Apiguard (thymol), Apivar (amitraz), Bayvarol (flumethrin), Apistan (tau-fluvalinate), formic acid products or oxalic acid products. Several calculations were performed to derive additional composite statistics. For variables that included multiple summer and fall data points (for example antibiotic and miticide usage), annual estimates were determined by calculating the mean values within each respective year at the province level. Total antibiotic usage statistics were determined by calculating the sum of oxytetracycline, tylosin and lincomycin usage within each respective year at the province level. A RCF index was developed to systematically quantify the top four possible causes of colony loss as determined by beekeepers; the first-ranked cause was assigned a value of 1, the second a value of 0.75, the third a value of 0.50 and the fourth a value of 0.25, ultimately resulting in four normalized rankings per year for each province with a total of eight ranked categories including ‘Weather’, ‘Starvation’, ‘Weak colonies in fall’, ‘Poor queens’, ‘Nosema’, ‘Varroa’, ‘Don’t know’ and ‘Other’. The final dataset containing the combined survey statistics used in this study is publicly available via GitHub at https://github.com/bdaisley/CAPABEESURV.

Climate data acquisition and processing

The Canadian Centre for Climate Services provides authoritative meteorological records from 8,254 stations nationwide with annual release of hourly, daily and monthly datasets88. Monthly summary datasets (CSV format) for temperature, precipitation and snowfall during the years 2014–2023 were accessed via the Government of Canada Historical Climate Data portal (https://climate.weather.gc.ca). The datasets were subsequently processed and merged with CAPA survey data in R (comprehensive details provided in Supplementary Data 1). In brief, 120 monthly summary CSV files were batch imported into a new R session with the read_csv function of the ‘readr’ package (v2.1.5) and the filter function of the ‘dplyr’ package (v1.1.4) to remove data points from northern provinces and territories not relevant to this study. After filtering, each of the individual datasets consisted of 1,172.6 ± 3.8 observations with 140,714 observations in total, for each of the variables of interest. Calculation of seasonal means for temperature, precipitation and snowfall was achieved through using the ‘group_by’ and ‘mutate’ functions of the dplyr package to summarize winter (December, January and February), spring (March, April and May), summer (June, July and August) and fall (September, October and November) data points within each province and for each year assessed. The resulting data frame, which consisted of 90 rows (10 provinces × 9 years) with four columns (winter, spring, summer and fall) for each of the relevant variables, was merged with the CAPA survey data via the ‘merge’ function in R.

Air pollution data acquisition and processing

The National Air Pollution Surveillance programme tracks ambient air quality measurements at 260 stations across Canada and publishes results annually through the Canada-Wide Air Quality Database89. Continuous monitoring datasets (XLSX format) for nitric oxides (NOX, including nitric oxide and nitrogen dioxide), O3, SO2 and PM2.5 were accessed for the years 2015–2022 via the Environment and Climate Change Canada Data Catalogue. Larger-sized particulate matter (<10 μm; PM10) was not evaluated as the majority of provinces lacked monitoring data for this metric. The datasets were subsequently processed and filtered before merging with CAPA survey data in R (comprehensive details provided in Supplementary Data 1). In brief, the datasets were downloaded using the ‘download.file’ function of base R in binary mode (mode = ‘wb’) and then imported into a new R session using the ‘readxl’ package. The province details were then appended to the dataset based on index matching to the ‘City’ columns via the ‘recode’ function of the dplyr (v1.1.4) package, followed by the calculation of per-province values of NOX, O3, SO2 and PM2.5 for each year (‘na.rm = TRUE’ was set to address missing values). The resulting data frame was then merged with the CAPA survey and climate data via the ‘merge’ function in base R and then exported in TSV format before downstream analysis.

Statistical analysis and data visualization

All statistical analyses were performed using R (v4.3.2) software. General descriptive statistics were generated with plotrix (v3.8.4) and dplyr (v1.1.4) package functions. Datasets with unique values were tested for normality using the Shapiro–Wilk test via the ‘shapiro_test’ function of the rstatix (v0.7.2) package. Normally distributed data were statistically compared with two-tailed t-tests, Pearson correlations, one-way ANOVAs or two-way ANOVAs as indicated. Non-parametric datasets were statistically compared using Wilcoxon or Kruskal–Wallis tests. Multiple comparisons were assessed using the Benjamini–Hochberg procedure to control the expected proportion of false discoveries.

Bar plots, scatter plots and line plots were generated using the ggplot2 (v3.5.0) package, with statistics annotated via the ggpubr (v0.6.0) package where applicable. Ridge plots visualizing the overwintering colony loss and boundary limit of sustainability were generated using the ggridges (v0.5.6) package. Geographic information system mapping of CAPA survey statistics was based on spatial data derived from GADM (v3.6; https://gadm.org) and was accessed in simple feature format via the ‘ne_states’ function of the rnaturalearth (v1.0.1) package. Maps were plotted on a Cartesian two-dimensional coordinate system, transformed to the geodetic coordinate reference system of ‘EPSG: 3348 – NAD83 (CSRS)/Statistics Canada Lambert’ according to Canadian census boundary files90 using the ‘st_transform’ function of the sf (v1.0.15) package, followed by boundary line simplification using the ‘st_simplify’ function (‘dTolerance = 5000’).

Mixed linear regression model analysis was performed to assess the impact of antibiotics as a predictor of overwintering mortality while controlling for environmental covariates. Our analysis involved a multilevel approach, where fixed effects included a range of covariates such as antibiotic use (quantified as combined estimates of oxytetracycline, tylosin, lincomycin and fumagillin), pesticide use, climate variables (temperature, precipitation and snowfall) and air pollution parameters (NO, NO2, SO2, PM2.5, CO and O3). Province was incorporated as a random effect to control for geographical differences, with Quebec and Newfoundland and Labrador excluded from analysis owing to deviation in antibiotic regulations compared with other provinces before 201834 and lack of antibiotic use due to absence of relevant diseases, respectively. Extreme outliers were detected and omitted from the analysis using the ‘is_extreme’ function of the rstatix (v0.7.2) package. Simple linear regression was performed with the ‘lm’ function of the stats (v4.3.2) package, while mixed linear regression was performed with the ‘lmer’ function of the lme4 (v1.1.35.1) package. A series of models were tested with varying complexity. The AIC and log-likelihood metrics were used to determine the best-fitting regression model. To further evaluate best-fitting models, we used an automated model selection approach using the MuMIn package (version 1.48.4) in R. Specifically, we used the dredge() function to systematically evaluate all possible subsets of predictor variables from the global model, ranking models based on AIC corrected for small sample sizes (AICc). Final summary statistics for regression models were generated via the ‘tab_model’ function of the sjPlot (v2.8.16) package.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.