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.