**3. Discussion**

#### *3.1. Selected Variables*

The Handy VarLoad (HVL) model allowed for the prediction of the Varroa load at a given moment *t*, as a function of the previously observed Varroa load and of the available area for their reproduction, i.e., the number of honey bee brood cells.

Seasons influence the Varroa load, but only in the short term. This could be explained by the fact that, in one month, a beekeeper management intervention or a particular climatic event can have an effect on one or two Varroa generations, as the generation interval of capped brood is 12 days. The mite population growth rate is exponential during short periods (three months) and when mite populations are low; in contrast, Varroa population growth follows a logistic dynamic over longer periods (covering the entire production period) when density-dependent factors influence population growth [12]. Consequently, an event which increases or decreases Varroa reproduction may change the short-term Varroa load but have an insignificant influence on the long-term Varroa load. For example, disruption of honey bee colony broods could be offset by the Varroa population growth itself. Conversely, if colony brood disruption speeds up Varroa reproduction, the mite population size eventually stabilizes due to density dependence [25,31]. Moreover, during a three-month period, colonies undergo a series of favorable and unfavorable disruptions for Varroa development, particularly climatic, which will balance each other out. Finally, the apiary effect acts regardless of the delay between two phoretic Varroa measurements. Thus, the biological variability between colonies, the differences in management strategy between beekeepers, year, and region (climate) influence the Varroa load of the colony [32–35].

Contrary to previous mathematical models on Varroa load, the HVL model allows one to obtain a prediction with a measure of uncertainty, as well as the associated uncertainty for each parameter. The model uncertainty includes variability at the inter-apiary scale, in beekeeping management strategies, and in year and region effects. The apiary effect included in the model induces a large amount of prediction uncertainty, but, at the same

time, it assimilates the sampling diversity related to apiary characteristics (management, year, and region).

#### *3.2. Beekeepers' Interest*

The model presented here allows one to have a representation of the risk beekeepers take by not treating the apiary, according to the percentage of colonies that exceed the threshold of three Varroa mites. Different quantiles propose different decision-making indicators for beekeepers taking into account trade-offs between cost, time, and environmental effects of treatments, on the one hand, and the risk of losing infested colonies, on the other.

Moreover, beside economic trade-offs, Varroa treatments are not without consequences and, indeed, may induce acaricide resistance in Varroa [3–5], which is why beekeepers should treat only when economic risks are real. It is worth noting that treatments during the beekeeping season are not efficient over the long term [36]. These types of treatment must be used only when the aim is to temporarily decrease the Varroa load to optimize honey flow performance. Thus, this model takes into account integrated pest management.

The model can also be used to determine which apiaries should be given priority on lavender and sunflower fields if the spot number is limited. However, despite the fact that managing colonies at the apiary scale is more efficient, as honey bee colony performances are highly dependent on the characteristics of any apiary (Kretzschmar et al., unpublished data), beekeepers may want to manage Varroa at the colony scale and thus strictly follow the model prediction.

#### *3.3. Limits and Prospects of the Model*

The choice to use only easy-to-measure variables in the field impairs the model's goodness of fit and, consequently, the estimation/prediction accuracy. Taking into account other variables (Varroa foundress density, Varroa infestation rate in the capped brood, natural death of Varroa mites measured on sticky boards, etc.) would have allowed better predicting the Varroa load. Including these additional variables in the present model could have easily improved its prediction power. Nevertheless, it would be far too long and complex to collect that type of data in the actual schedule of a beekeeper. If the sampling plan is unrealistic and impracticable at a large scale, the HVL model will be worthless. However, such improved models could be developed for researchers or technicians who work on a smaller scale and need to have better precision in their experimental frameworks. Another limit of this study is the sampled colony number: the more hives sampled, the better the estimation. In the present study, as the number of repetitions for each factor (management, year, and region) is limited, our sampling variation increased model uncertainty. Nevertheless, the Handy VarLoad model will be improved by the accumulation of data issued from the numerous experiments in which the two handy variables it uses (phoretic Varroa load and capped brood area) are commonly collected. As the database on which the model is based increases, the effect of covariates (apiary, region, season, beekeeping practices, etc.) can be better integrated.

#### **4. Materials and Methods**
