Statistical Analysis
Analyses were performed using the statistical program R (version 4.2.1; R Core Team 2020). First, we evaluated how pollinator visitation behavior varied among honeybees, bumblebees, and other pollinators with a separate model for the following response variables: 1) number of visits to flowers per 30 min video observation period, 2) total duration per visit of pollinator visits to flowers (seconds/visit), 3) duration pollinators spent on petals only per visit (seconds/visit), 4) duration pollinators spent on pollen only per visit (seconds/visit), and 5) duration pollinators spent simultaneously on pollen + nectar per visit (seconds/visit; model output in Appendix S1: Table S5). Each model was a zero-inflated generalized linear mixed effects model (GLMM) using a negative binomial distribution with a log link function and pollinator group (honeybees, bumblebees, and other pollinators) as the main predictor for both the zero-inflated and GLMM portions of the model (glmmTMB package; Brooks et al. 2017). The pollinator visitation data was aggregated by each flower observed for each pollinator group; therefore, a nested random effect of Flower ID (8 flowers/visit/site) within visit to a site (2 visits/site) within site (6 sites) was used in each model. To model the duration per visit, we used the duration of behaviors in seconds as the response variable with an offset of the log of the number of pollinator visits + 1 to correct for flowers with zero visits. We removed one outlier point from the number of visits per 30 min data, where other pollinators visited a single flower over twice as many times as the next most visited flower in our study. Each model was checked for overdispersion, zero-inflation, and spatial autocorrelation; none of these tests were significant (DHARMa package; Hartig 2020). Then we followed up each model with a post-hoc test to evaluate significant differences among honeybees, bumblebees, and other pollinators’ visitation behaviors (Appendix S1: Table S6; emmeanspackage; Lenth et al. 2020).
To evaluate V. ceranae prevalence in honeybees and bumblebees, we initially calculated the total apparent V. ceranae prevalence in each host species (epiR package; Stevenson et al. 2021) and used a Chi-squared test of two proportions to determine if there was a significant difference among the two host species. Then we ran two sets of models: the first to evaluate how the number of visits to flowers (per 30 min) at each site influenced V. ceranaeprevalence in honeybees and bumblebees, and the second to test how the duration per visit of specific behaviors on the flowers was correlated with V. ceranae prevalence in each host species. All models were generalized linear mixed effects models with V. ceranae prevalence in either honeybees or bumblebees as the response variable, a binomial distribution, and a logit link function (package lme4 ; Bates et al. 2015). For all models, we used a random effect of each visit to a site nested within site to account for both the variation between the two different dates on which each site was visited and the variation among different sites, though the random effects were singular for many models (model outputs in Appendix S1: Table S7 and S8). In the first set, models included the average number of honeybee visits, bumblebee visits, and all other pollinator species visits to flowers during the 30-min observation period as main effects. In the second set, we ran a series of models to evaluate the average total duration per visit (seconds/visit) and the durations per visit of petal-only, pollen-only, and pollen+nectar interactions of honeybees, bumblebees, and all other pollinators that visited the flowers. To deal with zeros in the data, all main effects had ‘1’ added to the value before log transforming the variable, and then they were scaled and centered to generate standardized estimates from the models. The variance inflation factor (VIF) for all main effects in all models was < 2.2, indicating that there was no multicollinearity in our models. Additionally, none of the models were over-dispersed. There was no evidence of spatial autocorrelation in the model residuals, indicating that V. ceranae prevalence was not correlated among sites based on their spatial proximity (DHARMa package; Hartig 2020). Finally, we used a Bonferroni Correction of four comparisons to adjust our alpha significance threshold from 0.05 to 0.0125 to account for four separate analyses, one for each of four different pollinator duration behavioral parameters (total duration per visit, duration on petals per visit, duration on pollen per visit, duration on pollen + nectar per visit). The number of visits model was not included in the Bonferroni Correction and was evaluated with the usual 0.05 alpha threshold.