https://orcid.org/0000-0002-6023-0719
Abstract
Managed honey bee (Apis mellifera L.) colonies in North America and Europe have experienced high losses in recent years, which have been linked to weather conditions, lack of quality forage, and high parasite loads, particularly the obligate brood parasite, Varroa destructor. These factors may interact at various scales to have compounding effects on honey bee health, but few studies have been able to simultaneously investigate the effects of weather conditions, landscape factors, and management of parasites. We analyzed a dataset of 3,210 survey responses from beekeepers in Pennsylvania from 2017 to 2022 and combined these with remotely sensed weather variables and novel datasets about seasonal forage availability into a Random Forest model to investigate drivers of winter loss. We found that beekeepers who used treatment against Varroa had higher colony survival than those who did not treat. Moreover, beekeepers who used multiple types of Varroa treatment had higher colony survival rates than those who used 1 type of treatment. Our models found weather conditions are strongly associated with survival, but multiple-treatment type colonies had higher survival across a broader range of climate conditions. These findings suggest that the integrated pest management approach of combining treatment types can potentially buffer managed honey bee colonies from adverse weather conditions.
Introduction
Honey bees (Apis mellifera Linnaeus, Hymenoptera: Apidae) are critical pollinators for agricultural systems (Calderone & Smagghe, 2012, Hung et al., 2018), but beekeepers have reported high losses among their colonies in recent years in North America and Europe (Bruckner et al. 2022, Gray et al. 2023). Beekeepers in the United States typically report losing 30% of their colonies on average each winter despite intensive management inputs (Kulhanek et al. 2017, Bruckner et al. 2022). A better understanding of the drivers of colony loss and their interactions could benefit beekeepers by informing in-hive management and apiary placement decisions. Additionally, wild bees and other pollinators are often exposed to similar environmental stressors as honey bees (Woodcock et al. 2017, Evans et al. 2018, Kammerer et al. 2021), though the extent to which honey bees can serve as an indicator of wild bee declines is an ongoing area of debate (Wood et al. 2020, Kammerer et al. 2021, Jackson et al. 2022).
The seasonal growth and associated behaviors of honey bee colonies in temperate regions are well understood (Winston 1987, Grozinger et al. 2014, Döke et al. 2015). Colony population size increases in spring, and colonies reproduce by colony fission (i.e., swarming) once the colony reaches a density threshold. Colony size continues to grow over the summer as foragers collect and store pollen and nectar. In the fall, colonies taper brood production and shift to producing long-lived winter bees, whose extended lifespan carries the colony through the winter until early spring when the queens start laying eggs again (Mattila et al. 2001).
In studies modeling colony loss at large geographic scales, weather variables consistently rank among top predictors (Switanek et al. 2017, Van Esch et al. 2020, Calovi et al. 2021, Insolia et al. 2022). Relationships between weather and honey bee health are often nonlinear or seasonally specific and likely reflect both direct physiological and behavioral effects on honey bees, as well as indirect effects through floral resource availability (Corbet 1990) and host–parasite dynamics with Varroa destructor, the preeminent pest of honey bees in temperate regions. For example, increased rates of winter colony loss have been associated with high fall temperatures (DeGrandi-Hoffman and Curry 2004, Smoliński et al. 2021, Overturf et al. 2022), fewer summer foraging days (Becsi et al. 2021), and extremes of summer temperature and precipitation (Calovi et al. 2021). Colony performance has also been negatively associated with above-normal spring precipitation (Quinlan et al. 2022, 2023). Uncoupling the direct and indirect effects of weather on colony performance and survival can be challenging in these large-scale studies.
Several other factors also catalyze winter loss, including large populations of V. destructor and nuanced interactions between Varroa, weather, and colony nutritional status. Winter is costly and high risk for honey bees as the colonies must cluster to maintain a steady temperature (Gates 1914) and are reliant on their stored resources and the longevity of winter bees. Varroa populations grow with colonies over the spring and summer months and can overwhelm an untreated colony by reducing the lifespan of winter bees (Dainat et al. 2012), and Varroa has been shown to contribute significantly to winter colony loss in various geographic regions (Dahle 2010, Genersch et al. 2010, Guzmán-Novoa et al. 2010, Molineri et al. 2018, Hernandez et al. 2022). Furthermore, weather and landscape conditions can exacerbate Varroa infestations. Longer growing seasons can lead to higher Varroa loads within honey bee colonies due to the extension of the brood-rearing period (DeGrandi-Hoffman and Curry 2004, Nürnberger et al. 2019). Both Varroa infestation and lack of forage in the spring and summer reduce population size and small colonies struggle to thermoregulate and have lower overwintering success (Owens 1971). Furthermore, poor nutrition disturbs honey bee colonies’ ability to withstand Varroa infestations; several studies describe a connection between protein supplemented diets, Varroa levels, and deformed wing virus (a virus vectored by Varroa mites) (Janmaat and Winston 2000, DeGrandi-Hoffman et al. 2010) as well as between landscape quality and Varroa parasitization (Dolezal et al. 2019, Hernandez et al. 2023).
Mite infestations may also alter the way bees are experiencing their environment. Varroa parasitization impacts foraging ability in honey bees by impeding learning (Kralj et al. 2007), flight duration and homing ability (Kralj and Fuchs 2006), landing ability (Muijres et al. 2020), and waggle dance communication (Pusceddu et al. 2021), thus reducing honey bees’ ability to take advantage of the available forage in their landscape. Varroa also increases bees’ vulnerability to temperature (Aldea-Sánchez et al. 2021) and water stress (Annoscia et al. 2012). Thus, there are several mechanisms whereby Varroa infestation and consequently Varroa management strategies may interact with landscape and weather conditions to influence colony survival, but few studies have demonstrated this at a landscape scale (but see Dolezal et al. 2019).
Examining how weather and landscape impact colony survival under different Varroa management regimes can provide insights into whether and how these different factors interact, potentially informing treatment decisions. Surveys of beekeepers provide a valuable tool to assess environmental drivers of colony loss alongside management (Bruckner et al. 2022, Gray et al. 2023). We leveraged a unique dataset, an annual survey of Pennsylvania beekeepers from 2017 to 2022, and combined the survey results with satellite-derived remote sensing data and seasonal forage availability datasets using a Random Forest model. We built upon previous analysis of this survey data (Calovi et al. 2021), which highlighted the importance of summer weather conditions to winter colony survival. Our analysis uses 3 additional years of data and investigates the effect of management by separating colonies into 2 models based on their level of Varroa treatment. This approach allowed us to holistically predict loss in Pennsylvania, first asking how management decisions impact winter colony loss, and then investigating whether these management practices interact with landscape and weather factors in multiplicative or additive ways.
Materials and Methods
Survey Data
An annual survey on winter honey bee colony loss and beekeeping management has been sent to Pennsylvania beekeepers each spring for over 15 years; we used 6 years of data from the 2017 to 2022 surveys as these were conducted online and included apiary locations for some sites. The survey consists of 36 questions about beekeeping management and colony loss including pre-winter and post-winter colony numbers (November and April, respectively), beekeeper years of experience, Varroa treatment (binary), type of Varroa treatment, supplemental feed (binary), type of supplemental feed, and an option to submit apiary coordinates. Anonymized survey data, after preprocessing (see below), is provided in Supplementary Table S1 in Supplemental Materials.
Preprocessing of Survey Data
From the survey responses, we calculated apiary survival rate (number of colonies post-winter/number of colonies pre-winter), Varroa treatment category, Varroa treatment types, and number of treatment types used. Mite treatment type was based on beekeeper responses to the question “what mite treatment(s) did you use?” for which participants could select multiple choices from a list. We categorized treatment into 3 categories of increasing intervention and toxicity: mechanical, soft chemical, and hard chemical (Underwood and López-Uribe 2019). Mechanical treatments described nonchemical interventions (i.e., drone brood removal, breaking of the brood cycle, screened bottom boards, powdered sugar). Soft chemical treatments were designated as naturally derived, organic compounds (i.e., thymol, formic acid, oxalic acid, or hop beta acids), which have been shown to have fewer lasting effects on bee health and less accumulation in honey and wax (Lindberg et al. 2000, Imdorf and Charrière 2003, Rosenkranz et al. 2010). Hard chemicals were synthetic chemicals (i.e., tau-fluvalinate, amitraz), which have historically been widely used and effective but are becoming less utilized due to the increasing prevalence of resistant Varroa populations (Hillesheim et al. 1996, Milani 1999) and concerns about negative impacts on bees due to persistence in wax (Bogdanov et al. 1998, Martel et al. 2007). Number of Varroa treatment types captures how many distinct treatment types a beekeeper used; this does not describe how many times a beekeeper treated for Varroa, which was not a question on the survey. For example, if a beekeeper responded that they treated with oxalic acid and formic acid, the number of treatment types is 2, even though a beekeeper may have used each treatment multiple times. Multiple-treatment type types within the same category, such as multiple soft chemicals, were considered as separate treatments, but multiple types of application for the same chemical, such as Oxalic (drip) and Oxalic (vapor), were not.
There were initially 3,230 beekeeper surveys with 23,922 total colonies recorded over the survey years of 2017–2022. Responses from beekeepers who reported using their colonies for migratory pollination services and beekeepers who reported having multiple apiaries were excluded. In the survey that was used, each beekeeper could only report 1 GPS location, and thus georeferenced locations of additional apiaries or migratory locations were not available. The single reported GPS location would not reflect the landscape and weather conditions, as they were exposed to throughout the year. Thus, the exclusion of this cohort ensures that the geospatial data corresponding to the points is accurate and that the statistical analysis of treatment and survival is reflective of apiaries within Pennsylvania, which was the focus of the study. Additionally, any records where the number of colonies in the spring exceeded the number of colonies in the preceding fall were excluded, as this could represent either incorrect record keeping or the addition of colonies from other apiaries. After these quality control measures, 3,047 responses with 11,799 colonies remained for analysis. The remaining apiaries ranged in size from 1 colony to 102 colonies, with a median of 3 colonies and a mean of 4.53. The survey data are anonymized, so while the same beekeepers are likely to repeat the survey in multiple years, we had to treat each entry as a unique apiary within each year.
We divided survey responses into subsets appropriate for 2 analyses: a nonspatial model examining the effects of beekeeper management and a spatial model including management and spatial variables (landscape and weather). Our nonspatial model included all 3,047 response records remaining after the filtering described above, while the spatial model included the 708 records for which the beekeeper treated for Varroa and provided apiary location information. For these georeferenced points, we standardized the coordinates and eliminated any that were outside Pennsylvania. To increase the strength of the model, each colony was used as an individual record in this analysis (N = 3,072), and survival was binary, where 0 was mortality and 1 was survival, as in Calovi et al. (2021).
Climate Variables
For each apiary location in the georeferenced dataset, we derived climate variables from the PRISM climate dataset, which provides gridded daily coverage of temperature and precipitation for the conterminous United States at approximately 4-km resolution (PRISM Climate Group 2014). To understand the specific seasonal effects of weather on colony survival, we calculated total precipitation and mean temperature for the summer (June, July, August) and fall (September, October, November) preceding the survey year, and the winter (December, January, February) and spring (March, April, May) of the survey year. We additionally included the BIOCLIM weather variables (Fick and Hijmans 2017) that were both relevant for plant and bee phenology and that had sufficient variability across Pennsylvania. These were growing degree days, annual precipitation, and temperature seasonality (SD × 100).
Landscape Variables
We used a 5 km buffer around each apiary to generate landscape metrics, since this contains 95% of a colony’s typical foraging (Visscher and Seeley 1982, Couvillon et al. 2015). We combined data on land cover with estimates of forage quality per land cover class to represent landscape-total forage available to honey bees (Lonsdorf et al. 2009, Koh et al. 2016). We used relative values of forage quality per land cover class per season from Koh et al. (2016) to convert the Cropland Data Layer (USDA NASS 2021) land cover to available forage. Prior to summing, we weighted values with an exponential decline function (Lonsdorf et al. 2009). We used CDL from the preceding year (since this was the time period when the bees foraged) and applied distance-weighting calculation before summing forage quality within 5 km of apiary locations. We calculated forage quality in spring, summer, and fall, with seasonal date ranges defined by Koh et al. (2016).
To capture developed land area, which presumably has low forage availability, we calculated the mean impervious surface value within 2.5 km of each apiary from the National Land Cover Dataset (NLCD) impervious surface map for 2019 (Dewitz and US Geological Survey 2021). The NLCD is only available every 5 years with the most recent data being from 2019, but we would not expect considerable changes over the study period. All model variables are presented in Table 1
Weather, landscape, and management variables included in RF models of winter colony survival. Seasonal weather variables as well as several annual BIOCLIM variables are included. Summer (June, July, August) and fall (September, October, November) variables are for the year preceding the survey; Winter (December, January, February) and Spring (March, April, May) are for the survey year. All forage variables are for the preceding year
| Category | Variable | Source |
|---|---|---|
| Annual weather | Annual precipitation | PRISM |
| Growing degree days | PRISM | |
| Seasonality | PRISM | |
| Seasonal weather | Spring total precipitation | PRISM |
| Summer total precipitation | PRISM | |
| Fall total precipitation | PRISM | |
| Winter total precipitation | PRISM | |
| Spring mean temperature | PRISM | |
| Summer mean temperature | PRISM | |
| Fall mean temperature | PRISM | |
| Winter mean temperature | PRISM | |
| Landscape | Spring forage score | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) |
| Summer forage score | NLCD & and estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Fall forage | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Impervious surface | NLCD | |
| Beekeeping management | Treatment type | Beekeeper survey |
| Supplemental feed type | Beekeeper survey |
| Category | Variable | Source |
|---|---|---|
| Annual weather | Annual precipitation | PRISM |
| Growing degree days | PRISM | |
| Seasonality | PRISM | |
| Seasonal weather | Spring total precipitation | PRISM |
| Summer total precipitation | PRISM | |
| Fall total precipitation | PRISM | |
| Winter total precipitation | PRISM | |
| Spring mean temperature | PRISM | |
| Summer mean temperature | PRISM | |
| Fall mean temperature | PRISM | |
| Winter mean temperature | PRISM | |
| Landscape | Spring forage score | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) |
| Summer forage score | NLCD & and estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Fall forage | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Impervious surface | NLCD | |
| Beekeeping management | Treatment type | Beekeeper survey |
| Supplemental feed type | Beekeeper survey |
Weather, landscape, and management variables included in RF models of winter colony survival. Seasonal weather variables as well as several annual BIOCLIM variables are included. Summer (June, July, August) and fall (September, October, November) variables are for the year preceding the survey; Winter (December, January, February) and Spring (March, April, May) are for the survey year. All forage variables are for the preceding year
| Category | Variable | Source |
|---|---|---|
| Annual weather | Annual precipitation | PRISM |
| Growing degree days | PRISM | |
| Seasonality | PRISM | |
| Seasonal weather | Spring total precipitation | PRISM |
| Summer total precipitation | PRISM | |
| Fall total precipitation | PRISM | |
| Winter total precipitation | PRISM | |
| Spring mean temperature | PRISM | |
| Summer mean temperature | PRISM | |
| Fall mean temperature | PRISM | |
| Winter mean temperature | PRISM | |
| Landscape | Spring forage score | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) |
| Summer forage score | NLCD & and estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Fall forage | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Impervious surface | NLCD | |
| Beekeeping management | Treatment type | Beekeeper survey |
| Supplemental feed type | Beekeeper survey |
| Category | Variable | Source |
|---|---|---|
| Annual weather | Annual precipitation | PRISM |
| Growing degree days | PRISM | |
| Seasonality | PRISM | |
| Seasonal weather | Spring total precipitation | PRISM |
| Summer total precipitation | PRISM | |
| Fall total precipitation | PRISM | |
| Winter total precipitation | PRISM | |
| Spring mean temperature | PRISM | |
| Summer mean temperature | PRISM | |
| Fall mean temperature | PRISM | |
| Winter mean temperature | PRISM | |
| Landscape | Spring forage score | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) |
| Summer forage score | NLCD & and estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Fall forage | NLCD & estimates of per-land cover class flower cover (Lonsdorf et al. 2009, Koh et al. 2016) | |
| Impervious surface | NLCD | |
| Beekeeping management | Treatment type | Beekeeper survey |
| Supplemental feed type | Beekeeper survey |
Statistical Analysis
To test the effect of Varroa treatment on survival, we used Mann–Whitney tests and Kruskal–Wallis 1-way analyses of variance in R Statistical Software v4.1.2 (R Core Team 2022). We estimated mean percent survival based on Varroa treatment (binary), Varroa treatment type, Varroa treatment category, number of Varroa treatment types, feed type, beekeeper experience (binned in 4 categories), supplemental feed (binary), supplemental feed type, and apiary size. We compared treatment type and category only among beekeepers who reported using only 1 type of treatment.
Random Forest Model
We adapted the Random Forest (RF) classification model from Calovi et al. (2021) using the ranger package in R (Wright and Ziegler 2017). Random Forest models are strong in characterizing nonlinear relationships with many correlated predictor variables (Berk 2008, Shoemaker et al. 2018). Due to the high degree of variability between years, a single year does not act as a good testing dataset for the model, so instead we used 10-fold cross-validation to estimate error, ensuring approximately equal stratification of training and testing datasets across years. We first used 10-fold validation to tune our model to the best number of trees and variables to sample in each tree (mtry) and then used 10-fold cross-validation to generate model accuracy estimates (Wong and Yeh 2020). We stratified the groups by year to ensure a similar proportion of each year represented in each group because the weather was substantially different between years.
To determine variable importance for each variable in the model, the RF model reports premutation variable importance, which is the mean decrease in prediction accuracy across all trees when a predictor variable is randomly ordered. In our full model of all treated colonies, growing degree days and diversity of treatment types had the highest importance. We then separated our data into single-treatment type and multiple-treatment type colonies. Again, number of treatment types refers to the number of different types of Varroa treatment reportedly used, not to the number of applications of treatment. For example, colonies within the single-treatment type model may have been treated once or multiple times with the same type of treatment, as data on treatment quantity and frequency was not included in the survey. Thus, in comparing variable importance between these 2 models, we did not assess the effect of treatment quantity but rather treatment diversity. We included 16 variables in the single-treatment type model and 17 in the multiple-treatment type model as number of treatment types was only relevant for the latter. We repeated calculations of variable importance within a 50-fold cross-validation on each of 10 different random seeds to ensure the stability of the results (Strobl et al. 2007) and took the mean permutation importance value to rank the variables. The variance between the importance values was very low among model runs (see Supplementary Fig. S1). We generated the plots using the ggplot2, cowplot, and pdp packages (Wickham 2016, Greenwell 2017), with additional refinements in Adobe Illustrator 26.2.1.
Partial Dependence
To examine the predicted relationship between individual variables and colony survival, we generated partial dependence plots (Greenwell 2017) for top variables for both the single-treatment type and multiple-treatment type models. We generated both single variable and bivariate partial dependence plots, to visualize how the predicted colony survival rate changes when either 1 or 2 variables, respectively, are varied while holding all other predictors constant.
Results
Beekeeping Management Variables
The average loss among nonmigratory, small-scale Pennsylvania beekeepers included in the survey was 50.3% with a standard deviation of 40.7; mean annual losses ranged from 43.4% in the 2017 survey to 55.3% in the 2019 survey. Apiaries that were ten colonies or fewer accounted for 91.9% of survey responses (2,801) and 79% (2,406) had 5 or fewer. Among the nonmigratory and nonmoved apiaries, the mean number of respondents to the survey was 508 per year. Most beekeepers in the survey (64.9%) had 1–5 years of experience once those with multiple apiaries and those with migratory operations were excluded.
Apiaries treated to control Varroa mites had higher survival (mean = 56.0%, SD = 39.2, n = 2,342) than untreated apiaries (31.1%, SD = 39.8, n = 705, P < 2.2e-16, Fig. 1A). We next separated apiaries into those that received a single type of Varroa treatment and those that received multiple types of treatment. The survey did not request information on how frequently apiaries were treated, and thus “single-type” treatment apiaries may have received that treatment multiple times in the year. The survival rate for apiaries treated with multiple types of Varroa treatment (mean = 67.1%, SD = 34.7) was significantly higher than apiaries treated with 1 type of Varroa treatment (mean = 51.5%, SD = 40.0, P = 0.0008, Fig. 1B). There was no difference between colonies treated with 2, 3, and 4 types of treatment (n = 581, 77, 5, χ2 = 0.93, P = 0.62, data not shown).
Box and whisker plots showing the statistical summaries of the winter colony survival rates relating to Varroa management practices. Boxes represent the interquartile range (IQR); the horizontal line represents the median and the whiskers represent minimum and maximum values 1.5 times the IQR. There was a significant increase in colony survival among beekeepers who treated for Varroa A) and beekeepers who used multiple types of treatment had higher colony survival than those who used no treatment or one type of treatment B). There was no significant difference between beekeepers who used hard and soft chemicals but both had significantly higher survival than untreated C). All treatment types had higher survival than none except hop beta acids D). *Degree of significance between variables connected with lines. (**** indicates p<0.0001, *** indicates p<0.001, ** indicates p <0.01 , * indicates p<0.05).
Among those apiaries where only 1 type of chemical treatment was used, there was a significant difference in survival rate between categories of treatment from a Kruskal–Wallis test (χ2 = 10.28, P = 0.01); however, the pairwise comparisons between categories showed no significant differences in survival between groups. This discrepancy is likely due to the small sample size of mechanically treated colonies (n = 21). There was no significant difference in survival between apiaries treated with soft chemicals and those treated with hard chemicals (P = 0.37; Fig. 1C). Apiaries treated with chemical treatments (hard and soft) had a mean survival of 52% and those treated with mechanical treatments had a mean survival of 36.6% though this difference was marginally nonsignificant (P = 0.117). There was no difference (P = 0.41) between colonies treated with mechanical treatment and untreated apiaries, though few apiaries were mechanically treated (n = 21).
There was a significant difference in survival rate among types (χ2 = 19.882, P = 0.0005) of treatment for singly treated colonies (Fig. 1D). Apiaries treated only with hop beta acids (n = 62) did not have significantly different survival than untreated apiaries (P = 0.26) and had significantly lower survival than those treated with formic acid (n = 649), thymol (n = 68), amitraz (n = 260), or oxalic acid (n = 596, P < 0.05). Apiaries treated with only formic acid had lower survival than those treated with only oxalic acid (P = 0.02), and there was no significant difference between the remaining treatment types. Significance values for all pairwise comparisons of both Varroa treatment and supplemental feed are presented in the Supplementary Materials.
There was a significant positive relationship between apiary size and survival (P = 0.04, data not shown) and no relationship between beekeeper years of experience and survival (P = 0.95, data not shown). There was a significant difference in survival among type of supplemental feed (χ2 = 69.053, P < 0.001, Supplementary Fig. S2), and in particular, survival of apiaries fed Fondant or Sugar Candy (mean = 55.8%, SD = 40.3) was higher than unfed apiaries (mean = 47.6%, SD = 41.4, Supplementary Fig. S2). There was no difference between unfed colonies and all fed colonies however (mean = 50.5%, sd = 40.7, P = 0.35). The analysis of supplemental feed (binary) was highly unbalanced, as only 185 records were associated with apiaries that were not fed supplemental feed, while 2,862 records were. There was no significant relationship between survival and number of supplemental feed types (P = 0.65). Significance values for all analyses are presented in Supplementary Material.
Effects of Landscape Conditions, Weather, and Management on Predicted Survival
Using a RF model, we examined the effect of 17 landscape, weather, and management variables on winter survival of apiaries managed using a single Varroa treatment type and apiaries managed using multiple types of Varroa treatment. For both models, we used the apiaries for which beekeepers reported that it was their only apiary, they were not part of migratory practice, and they volunteered their apiary location. The final single-treatment type model trained on 1,856 colonies used 2,500 trees with an mtry of 3 and the multiple-treatment type model trained on 1,054 colonies used 4,000 trees with an mtry of 3. The model for the single-treatment type apiaries predicted colony survival with a cross-validated out-of-bag (OOB) error of 19.5% and a prediction accuracy of 70.5%. The model for the multiple-treatment type apiaries predicted colony survival with a cross-validated OOB accuracy of 17% and a 76.3% prediction accuracy.
The most important variables for the single-treatment type model were growing degree days, seasonality, fall precipitation, spring precipitation, summer precipitation, and winter temperature (Fig. 2A). The most important variables for the multiple-treatment type model were summer forage, summer precipitation, growing degree days, annual precipitation, spring forage, and fall forage (Fig. 2B). The landscape variables had higher ranking importance in the multiple-treatment type model (3 of the top 6 variables) than in the single-treatment type model, where the top 7 variables were all weather-related.
Variable importance of input variables using the permutation importance metric for A) single-treatment type and B) multiple-treatment type RF models of winter colony survival.
The single-treatment type and multiple-treatment type models show similar trends for many of the weather variables including a positive relationship between spring, fall, and winter precipitation and survival and a negative relationship between summer precipitation and survival (Fig. 3). Growing degree days and annual precipitation have nonlinear trends with survival, showing positive associations at intermediate values and then negative relationship at the low or high extremes. Predicted survival also increases with seasonality, suggesting that areas with larger differences in mean temperatures between each season have higher survival. All forage variables were positively associated with winter survival.
Partial dependence of model prediction on seasonal weather variables, select annual variables, and summer forage. RF predicted survival for single-treatment type colonies (smoothed in orange) and multiple-treatment type colonies (smoothed in blue) have similar trends, but survival for colonies treated with multiple types of treatment is higher. Rug diagrams on the x axis show the distribution of data for each variable for the single-treatment type (orange) and multiple-treatment type (blue) training datasets.
Bivariate partial dependence plots show predicted survival based on both growing degree days and summer precipitation. The plots demonstrated that multiple-treatment type apiaries have higher survival rates across a greater range of both growing degree days and precipitation (Fig. 4). Survival is lowest for the single-treatment type model when precipitation is over 450 mm and growing degree days are outside of a middle range of approximately 2,700–3,200 degree days. Survival is lowest for the multiple-treatment type model when summer precipitation and growing degree days are both low. Multiple-treatment type apiaries had a higher survival across a geographic range in Pennsylvania than single-treatment type apiaries (Fig. 5).
Bivariate partial dependence plots showing the combined influence of summer precipitation and growing degree days on A) single-treatment type (A) and B) multiple-treatment type predicted survival. High survival is in blue, and low survival is in yellow, and data outside the range of the training data are cropped from the predicted results. Middle range precipitation and growing degree-day values predict higher survival for both models, and the multiple-treatment model has higher predicted survival rates across the range.
Predicted survival for 2022 across Pennsylvania for A) single-treatment type and B) multiple-treatment type treated apiaries from RF models based on exclusively climate variables. High predicted survival is in dark green, and low predicted survival is in white. The multiple-treatment type model has higher predicted survival across the state.
Discussion
In an analysis of 3,047 survey results from beekeepers in Pennsylvania, we found that beekeepers who used Varroa treatment in their apiaries had significantly higher survival than those that did not treat. This underscores the importance of this parasite and is consistent with previous studies (Imdorf and Charrière 2003, Rosenkranz et al. 2010, Jacques et al. 2017, Beyer et al. 2018). We also found no significant difference between “hard” and “soft” chemical treatments, which again is consistent with previous studies supporting the efficacy of both treatment types (Imdorf and Charrière 2003, Vandervalk et al. 2014, Sabahi et al. 2020, Qadir et al. 2021, Underwood et al. 2023). Soft chemicals were as effective at increasing winter survival as hard chemicals, which have been shown to leave long-term residue in the colony (Martel et al. 2007) and continuous exposure to both pyrethroids and organophosphates has generated resistance of Varroa to these chemicals (Rinkevich 2020). We also found that beekeepers who used multiple Varroa treatment types had significantly higher colony survival than those who used a single Varroa treatment type. This finding underscores the potential power of using treatments in combination, which has been shown to effectively diminish Varroa populations and is often recommended to beekeepers as an Integrated Pest Management strategy (Delaplane et al. 2005, Roth et al. 2020, Jack and Ellis 2021).
In order to understand how treatment type diversity may influence survival, we divided our data into 2 models, based on whether the beekeeper used 1 or multiple types of treatment. We found no evidence of interactive or multiplicative effects of treatment diversity and weather variables on survival. The partial dependence trends were similar across both models, demonstrating that all treated colonies have similar relationships between these weather variables and survival. However, the multiple-treatment type colonies had higher survival on average and higher predicted survival across a greater range of temperature and precipitation conditions. These results suggest that the relationship between Varroa treatment diversity and weather is additive rather than synergistic, as the variation in survival across conditions is similar among partial dependence plots for both models. Colonies treated for Varroa have been shown to have higher survival than untreated ones across geographies with a range of climatic conditions (Genersch et al. 2010, Guzmán-Novoa et al. 2010, Dainat et al. 2012, Nazzi et al. 2012, van Dooremalen et al. 2012, Molineri et al. 2018, Hernandez et al. 2022). Varroa infestations suppress immunity (Gregory et al. 2005, Zaobidna et al. 2017, Annoscia et al. 2019), reduce nutrition and foraging abilities (Kralj and Fuchs 2006, Dolezal and Toth 2018, Yang et al. 2021), and diminish the lifespan of honey bees (Kovac and Crailsheim 1988, Amdam et al. 2004, van Dooremalen et al. 2012), all of which can also be exacerbated by adverse weather conditions (Le Conte and Navajas 2008, Clarke and Robert 2018, Szentgyörgyi et al. 2018). With this in mind, our results suggest that colonies adequately treated to prevent Varroa may be buffered against some of these deleterious effects of weather. While we did not find a significant interactive relationship between Varroa treatment and seasonal weather, previous work has demonstrated that timing of Varroa treatments influences their efficacy (Strange and Sheppard 2001, Currie and Gatien 2006). This is likely largely due to seasonal variations in the Varroa reproductive cycle and development of honey bee brood (Rosenkranz et al. 2010), but Beyer et al. (2018) also found an interactive relationship wherein weather conditions influenced treatment efficacy and Desai and Currie (2016) found an interaction between season and wintering method on concentration of deformed wing virus.
Forage availability had higher variable importance in the multiple-treatment type model than in the single-treatment type model. This suggests that in colonies where Varroa was managed in a more integrated manner, using multiple-treatment types, the relative importance of weather versus landscape variables may have shifted, allowing the relationships between landscape and colony survival to be illuminated, although it should be noted that the treatment quantity and frequency is not known, and that treatment diversity may not be related to treatment intensity. Predicted survival increased with forage score for spring, summer, and fall in our study, which is consistent with previous work establishing links between forage availability in the landscape improving honey bee colony health (Donkersley et al., 2014, St Clair et al. 2020, Ochungo et al. 2022). Increased nutrition can mitigate the impacts of Varroa infestation as well (Aronstein et al. 2012, DeGrandi-Hoffman and Chen 2015, Dolezal and Toth 2018). Indeed, DeGrandi-Hoffman et al. (2020) found that survival of fed colonies was higher than unfed colonies, despite fed colonies having higher Varroa populations and deformed wing virus levels. Furthermore, Hernandez et al. (2023) demonstrated that increased forage quality led to higher populations of honey bees and consequently higher resistance in the population against Varroa. Future work should investigate these links, particularly the influence of landscape quality and supplemental feed on the ability of a colony to withstand Varroa. However, in our study, forage availability was less important for the survival of single-treated colonies than the weather was, suggesting that management strategies that keep Varroa below a certain threshold are necessary to obtain the benefits of landscape-scale floral resources.
Our results also support evidence found in other studies that the effect of weather on colony survival varies seasonally (Switanek et al. 2017, Beyer et al. 2018, Becsi et al. 2021) and highlights the value of including seasonal variables in models of colony survival. While spring, fall, and winter precipitation increased survival, summer precipitation decreased colony survival, consistent with a previous analysis of the first several years of this data (Calovi et al. 2021). Our finding of a detrimental effect of summer precipitation is consistent with some previous studies (Beyer et al. 2018, Becsi et al. 2021), but not others (Quinlan et al. 2023 found no effect of summer weather on winter cluster size). Many consecutive days of summer rain may decrease survival by reducing bee foraging, leading to reduced brood production and less food stored for winter. Our finding that spring, autumn, and winter precipitation supports higher survival is consistent with Switanek et al. (2017), which found that across all seasons, precipitation was linked to higher winter survival, but contrasts with other recent work showing that above-normal spring precipitation is linked to higher mortality (Quinlan et al. 2023, Quinlan and Grozinger 2023) possibly due to spring rain encouraging grasses, but not flowering plants. Furthermore, we found that winter temperature and survival had a nonlinear relationship, with high and low winter temperatures increasing colony loss. This may reflect that honey bees can succumb to overly cold winters (Southwick and Heldmaier 1987), but that there is a necessity for the temperature to stay consistently low enough for the bees to remain in their winter cluster to minimize energy consumption (Currie et al. 1992). Differences in findings about the effects of weather across these studies may be partially due to geographic differences between studies as different ecoregions can vary in overall baseline climatic conditions as well as plant communities and land use types, which in turn mediate the impacts of weather (Quinlan et al. 2022).
In diverse agricultural systems, deliberate and proactive management can reduce the impacts of climate change (Lipper et al. 2014, Malhi et al. 2021). In many examples of climate-smart agricultural adaptations, the extent of the management intervention contributes to the resilience of the system (Zilberman et al. 2018). For managed honey bee systems, controlling the introduced pest V. destructor augments overwintering survival (Rosenkranz et al. 2010, Jacques et al. 2017, Sabahi et al. 2020) and most beekeepers in the United States use some kind of treatment regime for the pest (Underwood et al. 2019). To our knowledge, our study is the first to show that not only treatment, but using multiple-treatment type types, bolsters winter survival of honey bees. The adaptive management capacity of a beekeeper has been demonstrated as integral to their success in the face of uncertain or variable conditions (Kouchner et al. 2019), and while weather conditions (which are outside of a beekeepers’ control) impacted colony survival considerably in our models, Varroa treatment is an effective way to increase survival across a range of otherwise-deleterious weather conditions.
Our study is limited in that we lacked data on some key aspects of beekeeping management and colony health, including mite counts, queen genotypes, honey harvest quantities, and number and timing of Varroa treatments applied. While the weather variables had higher predictive power on survival in the single-treatment type model than in the multiple-treatment type model, comparison of variable importance between models should be cautiously interpreted as the 2 models are trained on separate sets of data. Without data on treatment quantity, we cannot determine whether the results we report are effects of combining treatments or are caused by increased treatment quantity or time spent treated in colonies with multiple types of treatment. Beekeepers in the single-treatment type group may have treated repeatedly or heavily, while those with multiple-treatment types may have applied treatment less often, or vice versa. Furthermore, variable importance scores can be inconsistent in RF models, but we addressed this by calculating variable importance using 50-fold cross-validation on 10 random seeds and by using the unbiased permutation importance metric (Strobl et al. 2007). Our data are self-reported by beekeepers across Pennsylvania, which creates room for error and some responses, such as number of treatment types, may be confounded by general beekeeping effort and investment. Additionally, due to regulatory changes in the United States which have limited the amount of data that is publicly available on pesticide use in recent years (Hitaj et al. 2020), we were not able to estimate insecticide exposure. Finally, seasonal floral resources values are derived using expert opinion and based on the Cropland Data Layer and thus may not always capture habitat at a suburban and urban scale; future studies should incorporate finer resolution landscape data.
Overall, our study demonstrates the role of beekeeping management in buffering honey bee colonies against adverse weather conditions in temperate regions. In the context of the growing evidence that increasing climate variation may impact bee health (reviewed in Le Conte and Navajas 2008, Cunningham et al. 2022, de Jongh et al. 2022), these results underscore the ability of Integrated Pest Management strategies to increase the survival of honey bees even in adverse weather conditions.
Acknowledgments
We greatly appreciate the continuous efforts made by the beekeepers of Pennsylvania to document trends in winter survival through the annual survey, and we would like to especially thank Ken Hoover for his longstanding dedication to collecting these data and generously sharing it with us. We would like to thank Dr. Gabriela Quinlan for valuable discussions and critical reading of the manuscript.
Funding
This research is supported by the Food and Agriculture Cyberinformatics grant no. 2021-67021-34146 project from the USDA National Institute of Food and Agriculture to C.M.G. and S.G. This research was supported in part by the US Department of Agriculture, Agricultural Research Service. The findings and conclusions in this manuscript are those of the author(s) and should not be construed to represent any official USDA or US Government determination or policy. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. USDA is an equal opportunity provider and employer.
Author Contributions
Darcy Gray (Data curation [Equal], Formal analysis [Lead], Investigation [Equal], Methodology [Equal], Software [Equal], Validation [Equal], Visualization [Lead], Writing—original draft [Lead]), Sarah Goslee (Data curation [Equal], Formal analysis [Supporting], Investigation [Equal], Methodology [Equal], Software [Equal], Validation [Equal], Visualization [Supporting], Writing—review & editing [Equal]), Melanie Kammerer (Data curation [Equal], Methodology [Equal], Writing—review & editing [Equal]), and Christina Grozinger (Conceptualization [Equal], Data curation [Supporting], Funding acquisition [Equal], Investigation [Equal], Methodology [Equal], Project administration [Equal], Resources [Equal], Supervision [Lead], Writing—review & editing [Lead])
