**1. Introduction**

The parasite *Varroa destructor* is considered a major pathogenic threat to honey bees [1] and to beekeeping. This mite is an ectoparasite affecting both adult bees and broods. Female mites have two distinct stages: a phoretic stage on adult bees and a reproductive stage, which takes place inside a capped brood during bee metamorphosis. The Varroa threat is not new for the beekeeping community, but with colony importations and the commerce of bees, this threat continues to increase. Indeed, these circumstances favor the Varroa spread throughout territories and the world's apiaries. This threat is all the more important given that the parasite spread is rapid [2]. Thus, bees and beekeepers cannot adapt and respond efficiently; on the contrary, Varroa, with continuous exposure to miticide treatments, responds with mechanisms of resistance [3–5]. Consequently, the current challenge is to develop new methods to limit Varroa numbers inside colonies.

**Citation:** Dechatre, H.; Michel, L.; Soubeyrand, S.; Maisonnasse, A.; Moreau, P.; Poquet, Y.; Pioz, M.; Vidau, C.; Basso, B.; Mondet, F.; et al. To Treat or Not to Treat Bees? Handy VarLoad: A Predictive Model for *Varroa destructor* Load. *Pathogens* **2021**, *10*, 678. https://doi.org/10.3390/ pathogens10060678

Academic Editor: Giovanni Cilia

Received: 21 April 2021 Accepted: 26 May 2021 Published: 30 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Without regular efficient treatment against this mite, honey bee colonies can collapse within a 2–3-year period in temperate climates. Varroa feeding on pupal hemolymph can induce a decrease in adult bee body weight and malformations as well as reducing their life spans, thus weakening their immune systems [6–8]. Thus, it seems logical that infested colonies are less productive and efficient than healthy colonies, which can have significant economic impacts for beekeepers [9,10]. Beyond a threshold of 3 phoretic Varroa mites per 100 bees, the decrease in performance is correlated with the Varroa load [10]. According to this study, a colony with more than 3 phoretic Varroa mites per 100 bees produces, on average, 2.65 kg less honey than a colony below this threshold. Unfortunately, until now, it has not been possible to predict, from the mite population size in the spring, the population load in the summer, despite studies by Arechavaleta-Velasco and Guzman-Novoa (2001), Harris et al. (2003), and Lodesani et al. (2002), confirmed significant correlations between the amount of brood and/or the fertility of the mites [11–13] and population growth [1].

Models of Varroa dynamics have been previously established but mainly carried theoretical descriptions and only allowed for the evaluation of the instantaneous Varroa load. Wilkinson and Smith's model [14] was built from virtual colonies, and DeGrandi-Hoffman and Curry's model [15] was based on the BEEPOP honey bee colony population dynamics model [16]; a BEEHAVE Varroa unit was developed by Becher et al. [17]. These models were based on parametric values available from previous studies [16,18–26]. As these models primarily work with mathematical extrapolation, instead of being dataderived, we assumed that the resulting parametric values could be revised. Additionally, in the twenty years since these models were published, Varroa biology may have coevolved with its host. The coevolution between Varroa and honey bees has been reported by Kurze et al. [27] and includes host resistance behaviors, which involve a decrease in the Varroa reproduction rate as well as perturbations in the biological cycle of the mite. Moreover, previous studies serving as the basis for model construction were based on honey bees with different European origins and on Africanized honey bees [19]. Honey bee origins affect Varroa reproduction [28] and, consequently, Varroa population sizes. To increase its predictability, here, we used a model based on empirical data.

The most important information for a beekeeper is not the Varroa load at the time of honey flow because most treatment compounds, even some labeled "natural" (e.g., formic acid or thymol), are banned or not recommended during honey flow [29]. The aim of this study was therefore to predict the Varroa load one or three months later, from its baseline level in early spring, to anticipate colony performance for honey flow, knowing that the reduced performance threshold is 3 phoretic Varroa mites per 100 bees. Aimed as a useful tool for beekeepers, the model Handy Varload is based on inexpensive, accessible, and quickly measurable data in the field.

#### **2. Results**

#### *2.1. Variable Selection (for Variable Definitions, See Materials and Methods, Statistical Analysis)*

The variable "phoretic Varroa" measured at *t* = 0 was continuous with 25% of zeros, 27% of the data in [0, 1], and 48% of the data in [1,30]. The zero-inflated beta distribution is similar to the beta distribution but allows zeros as response values in which the ν parameter models the probability of obtaining zero. The distribution features of the variable "phoretic Varroa" (*Vpt*) require dividing by 100 in order to fit the data to the interval [0, 1]. We then modeled this new response variable by a zero-inflated beta distribution, with parameter variation depending on covariates. A first model selection was performed to choose the best variables to model μ (see AICc comparisons in Table 1; more details are provided in Supplementary Materials Table S1). At the end of this preliminary selection, two models including the apiary random factor were retained, one for the 1-month adjustment and the other for the 3-month adjustment, noted (\*) and (\*\*) in Table 1, respectively.

**Table 1.** Comparisons of the tested models investigating the influence of phoretic Varroa numbers (per 100 bees) at *t* = 0, capped brood cell numbers, varbrood, and date of predicted phoretic Varroa numbers as a function of the estimation length, using the AICc criterion. N = 867 for data adjustment at one month (x = 1) and N = 93 for data adjustment at three months (x = 3).


The final models (A and B, see below) were obtained after a second variable selection based on AICc comparisons, using the modeled σ and ν added to the (\*) and (\*\*) preliminary models. The number of phoretic Varroa present at *t* was modeled by the following zeroinflated beta models (BEZI in "*gamlss*"):


*For data adjustment at one month*: (A)

<sup>μ</sup> = logit−<sup>1</sup> (α<sup>0</sup> <sup>+</sup> <sup>α</sup>1*Vb*t−<sup>x</sup> <sup>+</sup> <sup>α</sup>2*Cp*t−<sup>x</sup> <sup>+</sup> <sup>α</sup>3*D*<sup>t</sup> <sup>+</sup> *Ap*) σ = exp (β<sup>0</sup> + β1*Vb*t−<sup>x</sup> + β2*D*<sup>t</sup> + *Ap*) <sup>ν</sup> = logit−<sup>1</sup> (γ<sup>0</sup> <sup>+</sup> <sup>γ</sup>1*Vb*t−<sup>x</sup> <sup>+</sup> <sup>γ</sup>2*Cp*t−<sup>x</sup> <sup>+</sup> <sup>γ</sup>3*D*<sup>t</sup> <sup>+</sup> *Ap*)

*For data adjustment at three months*: (B)

<sup>μ</sup> = logit−<sup>1</sup> (α<sup>0</sup> <sup>+</sup> <sup>α</sup>1*Vb*t−<sup>x</sup> <sup>+</sup> *Ap*) σ = exp (β<sup>0</sup> + *Ap*) <sup>ν</sup> = logit−<sup>1</sup> (γ<sup>0</sup> <sup>+</sup> <sup>γ</sup>1*Vb*t−<sup>x</sup> <sup>+</sup> <sup>γ</sup>2*Vp*t−x)

where the α, β, and γ parameters are coefficients used to model μ, σ, and ν, respectively. As a consequence of this second variable selection, the final AICc was −3179.2 (A) and −343.5 (B) (see details in Tables S2 and S3). Varbrood, which was retained by model selection in all cases except for σ of model B, appeared as the most important explanatory variable.

#### *2.2. Goodness of Fit and Prediction Evaluation*

#### 2.2.1. Parameter Uncertainty

For models A and B, the α, β, and γ parameters associated with μ, σ, and ν were estimated and their CI95%s were computed (see Table 2). Based on the intercept, we can note that varbrood had the largest influence on the data adjustment for each parameter of model A. Moreover, the order of influence of model covariates was the same regardless of

the parameter: varbrood > date > capped brood cells. For ν of model B, phoretic Varroa had a larger influence than varbrood. The CI95 range as positively correlated with covariate weights, i.e., the greater the weight, the larger the uncertainty.

**Table 2.** Estimated coefficient and 95% confidence interval (CI95%) of models A and B investigating the influence of varbrood, capped brood cells, phoretic Varroa, and date on the number of phoretic Varroa mites for mu, sigma, and nu parameters.


Moreover, the apiary effect depended on the horizon of prediction. Thus, the mean apiary effect was zero with varying estimated standard deviations depending on the data adjustment; at one month, the estimated standard deviation was 0.285, with a standard deviation of this estimate of 0.798, and at three months, the estimated standard deviation was 0.681, with a standard deviation of this estimate of 0.914 (see Supplementary Materials Tables S2 and S3).
