Adaptive population structure shifts in invasive parasitic mites, Varroa destructor

Abstract Comparative studies of genetic diversity and population structure can shed light on the ecological and evolutionary factors governing host–parasite interactions. Even though invasive parasites are considered of major biological importance, little is known about their adaptative potential when infesting the new hosts. Here, the genetic diversification of Varroa destructor, a novel parasite of Apis mellifera originating from Asia, was investigated using population genetics to determine how the genetic structure of the parasite changed in distinct European populations of its new host. To do so, mites infesting two categories of hosts in four European regions were compared: (a) adapted hosts surviving through means of natural selection, thereby expected to impose strong selective pressure on the mites, and (b) treated host populations, surviving mite infestations because acaricides are applied, therefore characterized by a relaxed selection imposed by the host on the mites. Significant genetic divergence was found across regions, partially reflecting the invasion pattern of V. destructor throughout Europe and indicating local adaptation of the mite to the host populations. Additionally, varying degrees of genotypic changes were found between mites from adapted and treated colonies. Altogether, these results indicate that V. destructor managed to overcome the genetic bottlenecks following its introduction in Europe and that host‐mediated selection fostered changes in the genetic structure of this mite at diverse geographic scales. These findings highlight the potential of parasites to adapt to their local host populations and confirm that adaptations developed within coevolutionary dynamics are a major determinant of population genetic changes.


| INTRODUC TI ON
Coevolution is a process of reciprocal evolutionary changes between interacting species (Ehrlich & Raven, 1964;Thompson, 2005). Adaptations developed within host-parasite interactions strongly shape the genetics of the coevolving species (Anderson & May, 1982;Thompson, 2005;Woolhouse et al., 2002). Notably, the rate of genetic changes is expected to be particularly swift in parasite populations compared to their hosts because of their shorter generation time (Paterson et al., 2010;Schmid-Hempel, 2011). Yet, the selection of specific parasite adaptations is not uniform for parasites infesting distinct populations of hosts, because different selective forces may act on these populations (Gandon & Zandt, 1998). These forces, composed of diverse elements of environmental factors or unique host population adaptations, generally result in the genetic diversification of parasite populations at varying geographic scales, thereby generating geographic mosaics of coevolution (Thompson, 2005).
Although adaptations of parasites to their host traits have been studied extensively in silico (Gandon & Michalakis, 2002) and in vitro (Brockhurst & Koskella, 2013), evidences of mosaics of coevolution in invasive parasites infesting novel hosts remain scarce. The relatively recent coevolving system between the Western honey bee, Apis mellifera, and the obligate ectoparasitic mite, Varroa destructor, represents an ideal opportunity to investigate the impact of host adaptations on parasite evolution in real time (Dietemann et al., 2012;Oldroyd, 1999). Originally, V.
destructor infested colonies of Eastern honey bees, Apis cerana, in Asia (Oldroyd, 1999;Rath, 1999). However, the introduction of A. mellifera colonies in the native range of the parasite resulted in the mite shifting host in the middle of the 20th century, and in its spread to almost all locations where beekeeping with the Western honey bee is practiced (Traynor et al., 2020). The spread of these parasites was particularly swift in Europe. Having been first detected in eastern regions of the continent in the beginning of the 1970s, it had dispersed throughout most of the continent just two decades later (Rosenkranz et al., 2010).
Out of the many V. destructor haplotypes found in A. cerana colonies (Navajas, 2010), only two are known to have switched to A. mellifera and to have emigrated from Asia (Anderson & Trueman, 2000). In addition to this original bottleneck, and coupled with the incestuous mating behavior of the mite (Rosenkranz et al. 2010), the host shift and subsequent spread of the pest has resulted in a highly homogenous genetic population structure within the invasive populations of V. destructor (Solignac et al., 2005).
Despite this low genetic diversity, V. destructor has flourished as an invasive parasite, with devastating consequences for its new host (Neumann & Carreck, 2010;Potts, Biesmeijer, et al., 2010) and the quasi-eradication of wild honey bee populations in Europe (Jaffé et al., 2010;Moritz et al., 2007). In A. mellifera colonies, mite population growth is exponential (Martin, 1998) and, due to its ability to vector lethal honey bee viruses Traynor et al., 2020), survival of susceptible colonies is dependent on annual mite control treatments implemented by beekeepers (Boecking & Genersch, 2008;Le Conte et al., 2010).
These regular mite treatments may remove the selective pressure imposed by the parasite (May & Anderson, 1990;Schmid-Hempel, 2011), thereby limiting the full potential of host adaptations (Fries & Camazine, 2001;Neumann & Blacquière, 2017) as well as the selection pressure of the host on V. destructor. Only occasionally were populations of European Western honey bees in France, Sweden, The Netherlands, and Norway left untreated and had the opportunity to adapt by means of natural selection (Fries et al., 2006;Kruitwagen et al., 2017;Le Conte et al., 2007;Locke, 2016;Oddie et al., 2017;Panziera et al., 2017). These adapted small honey bee populations, in France, The Netherlands, Norway, and Sweden, are considered resistant as they are surviving without the need for acaricide treatments by expressing a wide variety of traits that enable them to interfere with V. destructor population growth (Locke, 2016;Mondet et al. 2020).
Given that populations of V. destructor were genetically homogenous upon their introduction in Europe (Solignac et al., 2005), the adaptation potential of this invasive parasite has received little attention (Eliash & Mikheyev, 2020). Yet, the reduced diversity and common origin of the recently introduced mites represents a perfect starting point to study how adaptation to novel host populations and different environments can induce genetic diversification in the parasite. More specifically, when coevolving with adapted A. mellifera colonies, mites would be expected to swiftly develop counteradaptations against the resistance traits of the hosts in order to survive. Such adaptations should result in changes in the levels of genetic diversity and population structure of the mites (Schulte et al., 2010(Schulte et al., , 2013. For instance, the length of brood production during the season, which is a trait known to enable honey bee survival to V. destructor (Locke, 2016), may influence the yearly number of mite reproductive cycles and the level of recombination of inbred mite lineages (Beaurepaire et al., 2017). In contrast, mites infesting susceptible colonies that require regular treatments may face different selective forces. In this case, treatments may relax the selective pressure of mite infestation on the hosts but are instead a strong selective pressure on the parasite population, which could also lead to the development of resistance toward treatments (González-Cabrera et al., 2018;Martin, 2004;Milani, 1999). Consequently, the genetic structure of V. destructor populations infesting honey bee colonies that are adapted to the mite versus those that are susceptible and require mite control for survival may follow different evolutionary paths and result in a geographic mosaic of coevolution (Thompson, 2005). To investigate this hypothesis, we performed an analysis of the genetic diversity and population structure of mites infesting five adapted and five treated honey bee populations located in four European regions. Our results show significant genetic diversification across V. destructor populations in the studied regions and host populations, thereby suggesting that mites have overcome the initial bottleneck of their introduction and are adapting to their local host populations as predicted by the geographic mosaic of coevolution theory (Thompson, 2005).

| Sampling
In this study, adapted honey bee populations are defined as groups of colonies that have survived V. destructor infestations without the need for treatments for more than ten years, and do not require treatments to survive for extended periods. In contrast, the susceptible host populations are defined as treated since they require frequent treatments and/or management practice in order to survive mite infestations. Such treated populations generally show significantly lower expression levels of mite resistance traits compared to adapted ones . In 2017, adult female V. destructor mites were sampled from 32 treated and 28 sympatric adapted A. mellifera colonies, in various locations across four different regions: France, Netherlands, Sweden, and Norway (Tables 1 and 2). Mite samples were initially divided into 11 groups (Table 1), reflecting the apiary location and type of colony they were infesting (i.e., adapted or treated).
The treatments of the treated colonies differed across the studied regions ( Table 1). The adapted and treated colonies were located in the same apiary at two locations (Avignon, France; Lelystad, Netherlands). At all other locations, the distance between treated and surviving colonies ranged from 2.5 to 325 km ( Table 2). The mites (N = 1,310) were collected on adult workers using standard methods (i.e., powdered sugar, Dietemann et al., 2013) during the summer of 2016 and 2017 and were immediately transferred into 95% EtOH and stored at −20°C until DNA extraction.

| DNA extraction and genotyping
Established protocols were followed to isolate total mite DNA (Beaurepaire et al., 2017). In brief, mites were washed twice in ddH 2 0 to remove the EtOH and then individually distributed in 96-well plates filled with 100 µl of Chelex™ solution. Individual mites were then crushed with sterile pipette tips, 5 µl of 10 mg/ml Proteinase K was added, and their DNA was extracted following Walsh et al. (2013).
Initially, 20 microsatellites (Beaurepaire et al., 2017;Cornman et al., 2010;Evans, 2000) were tested on 12 individual mites from each location (N = 132) to assess the genetic diversity and population structure of the V. destructor samples. PCRs were conducted as detailed in Beaurepaire et al. (2019). Seven of the tested markers were polymorphic over all regions and were thus chosen for genotyping all samples (N = 1,310, Table 3). Twenty-four mites were genotyped for each colony whenever the infestation level allowed to.
PCR products were sent to Genoscreen (Lille, France) to run on a 3730XL sequencer (Applied Biosystems®, Carlsbad, CA). All mites were genotyped using the Peak Scanner TM software v 1.0 (Applied Biosystems®, Carlsbad, CA).

| Microsatellite DNA analyses
In total, 1,310 mites were individually genotyped at seven polymorphic loci (Table 1). Notably, the data from Sweden (165 mites) have already been published in Beaurepaire et al. (2019) and are used here to compare the amplitude of genetic changes across distinct mite groups. To verify the independence of the markers used, all locus pairs were tested for linkage disequilibrium using the software Fstat V 2.9.3 (Goudet, 1995).
To assess genetic differences between V. destructor infesting adapted versus susceptible A. mellifera populations in the different locations, the mites were initially grouped depending on apiary location and the type of host colonies they infested (i.e., adapted or treated; Table 1). To confirm this a priori sample clustering, the levels  (Goudet, 1995). These results were statistically compared using Kruskal-Wallis tests using R v. 3.5.1 (R Core Team, 2018). In order to understand whether the lowest sample size used in this study (i.e., Gotland treated mites, N = 38) was suitable to accurately represent the genetic diversity of the mite groups, rarefaction analyses were conducted using the ADZE software (Szpiech et al., 2008), using the allelic richness averages obtained across all markers in each mite group.

TA B L E 1 Sampling overview
To test the allelic divergence of V. destructor across regions and between honey bee groups within locations, estimates of D est (Jost, 2008) and pairwise tests of population differentiation were obtained for each pair of location (i.e., grouping individuals from the different host populations in each location) and within location for each possible adapted versus treated comparison using the software GenAlex v. 6.5 (Peakall & Smouse, 2012). The statistical significance of the pairwise population divergence indexes was obtained using Fstat v. 2.9.3 (Goudet, 1995) after 55,000 permutations, as this software allows to correct for multiple comparisons.
Additionally, the level of genetic differentiation (D est ) between mite groups treated with different acaricides was compared using a Student t test. Finally, to investigate the relationships between the genetic distance obtained with D est and the spatial distance separating the mite groups (see Table 2), a Mantel test of correlation (Mantel, 1967) was performed with GenAlex v. 6.5. For mite groups that were present on the same apiary (i.e., Avignon, France, and Lelystad, The Netherlands), a distance of 10 meters was considered.
Finally, the diversity and the prevalence of mite genotypes infesting the different honey bee populations were compared. To do so, the distribution of multilocus genotypes (MLGs) in each mite group was computed using the R package poppr (Kamvar et al., 2014). Comparison of the diversity of MLGs across adapted and treated populations was done by estimating the 95% confidence intervals of Shannon (Shannon & Weaver, 1949) and Simpson (Simpson, 1949) diversity indexes using 1,000 bootstraps. In parallel, the differences in the distribution of the most common MLGs between the adapted and treated host populations within each region were tested. To do so, only the mite genotypes with a frequency of at least 5% in the considered locations were compared using Fisher exact test conducted using contingency tables in R. As tests were conducted pairwise between the populations of any given region, the p-values obtained were corrected for multiple comparisons using Holm's method (Holm, 1979). In this analysis, all individuals with missing data were excluded from the dataset, resulting in a total of 863 individuals.

| RE SULTS
No significant linkage disequilibrium between pairs of markers was detected after correction for multiple comparisons (all pvalues > .05). The AMOVA confirmed our a priori grouping, indicating that sampling regions and host groups (i.e., treated versus adapted) were significantly (p <.01) structuring V. destructor populations (i.e., responsible for 41.9% and 7.6% of the total genetic variance observed, respectively; Table 4).
Variable, but low levels of observed heterozygosity (H o = 0.002-0.065) and number of alleles (NA = 2-4) were found across the seven markers over all samples (Table 3). For all mite groups, the rarefaction analyses showed that allelic richness only increased marginally after ~ 40 mites were analyzed ( Figure 1).
The comparison of allelic divergence of V. destructor populations in Europe revealed diverging patterns across the sampled regions. This analysis revealed very low allelic divergence between France and the Netherlands (D est = 0.01, p <.001), while mites from Sweden were more markedly differentiated from mites of these two regions (D est = 0.12, p <.001). Surprisingly, the genetic divergence levels between samples from these three locations and mites from Norway were thrice as high (D est = 0.32-0.38, p <.001; Figure 2).
The comparison of genetic diversity across mite groups infesting the different host populations revealed that the mean number of alleles and the observed heterozygosity levels did not differ significantly between these groups (Kruskal-Wallis tests, p >.05; Figure 3). However, the levels of genetic differentiation between mites sampled from adapted and treated populations within each region resulted all significant, with diverse D est levels ( Figure 2).
Notably, the degree of genetic differentiation weakly but significantly correlated with the distance separating the groups (Mantel test, R 2 = 0.14, p =.001; Figure 4). However, looking more closely at the pairwise comparisons between mite infesting different host groups revealed that mite genetic differentiation is not clearly determined by geographic distance. For instance, the level of divergence between mites in the two treated French populations was very low (D est = 0.01) and nonsignificant (Table 5) despite the distance separating these two groups (45 km). In contrast, higher and significant differences (D est = 0.029-0.048, p <.001) were found when these two groups were compared with the mites infesting the French adapted population, despite the fact that one of the treated host group was in the same apiary as the adapted honey bees. Also, the lowest level of divergence in the Netherlands occurred between the two adapted groups (D est = 0.013, p <.001), in spite of the ~100 km separating them. In parallel, the pairwise allelic comparisons revealed that mites from the adapted host population in France and from Tiengemeten in the Netherlands were not significantly genetically different (D est = 0.002, p =.056, Table 5), while mites infesting susceptible populations in these two regions were (D est = 0.07-0.08; p <.001, Table 5). Moreover, the average level of population divergence found between mites treated with Amitraz and mites groups treated with oxalic acid (average D est = 0.21, ±0.21 SD) did not differ significantly from that obtained when mite groups treated with oxalic acid were compared with each other (average D est = 0.27, ±0.14 SD, t test, p =.63). varied between groups, but only represented a minor proportion of the sampled mites (0.5%-4.1%, Table 6).

TA B L E 4 AMOVA results
Overall, contrasting patterns of MLG diversity could be observed when comparing mites from adapted and treated populations across the different locations ( Figure 5). Notably, in France, significantly higher levels of the Shannon index (H) were found in mites infesting the adapted host population (95% CI = 2.9-3.3) compared to the two mite groups sampled from the treated populations of the same region (95% CI = 2.0-2.5 and 2.4-2.8). The same pattern was found in one adapted population of the Netherlands (95% CI = 3.1-3.5) when compared to the adapted (95% CI = 2.4 -2. 8) and treated populations (95% CI = 2.4-2.9) from the same region. In addition, the distribution of the dominant MLGs varied significantly between adapted F I G U R E 2 Varroa destructor population structure across locations and groups of colonies. Results of the tests of allelic divergence (D est ) between mites infesting colonies of A. mellifera at the four locations and between treated and adapted colonies within each location. The thickness of the arrows linking populations shows the level of allelic divergence between them, from low (D est < 0.05) to high (D est > 0.25), while the dashed and solid lines represent statistical nonsignificance and significance, respectively (*** indicates p <.001). For every region, significant population structure divergences were found when adapted and treated mites were compared (p <.001). Notably, Norwegian mites were genetically highly isolated from all the other mite groups (D est = 0.32-0.38)

F I G U R E 3
Estimates of genetic diversity. Mean (± standard error) number of alleles and observed heterozygosity for each group of mites. Mites were grouped by location (countries are distinguished by color and locations within countries by color shades) and by the type of colony they infested (designated by symbols). No significant difference was found when mite groups were compared across regions and between populations of the same region (Kruskal-Wallis test, p > .05) and treated colonies in each location (Fisher exact tests corrected for multiple testing, p <.05, Figure 6).

| D ISCUSS I ON
Our results suggest that the genetic structure of European V. destructor is shaped by their interactions with their local host populations. These findings reveal that sympatric and allopatric mite populations are experiencing significant levels of genetic divergence that are probably caused by both unique host population traits and environmental differences across locations.
A weak and significant interaction (R 2 = 0.14, p =.001) was found between the genetic and geographic distance separating the mite groups, indicating that geographic isolation can explain a small part of the genetic differences found between distant mite populations.
Additionally, the pattern of genetic differentiation documented here appears to be further explained by the invasion history of the parasite in Europe. Historical reports of this invasion showed that V. destructor arrived in Europe through two routes. First, the mite is believed to have been introduced in 1971-72 in Eastern European countries and to slowly spread to the Northern regions of the continent, including Sweden (Griffiths & Bowman, 1981;Rosenkranz et al., 2010). Secondly, the parasite was introduced in Germany in 1977, a source from which mites from France and Netherlands derived (Nixon, 1983;Potts, Roberts, et al., 2010;Ruttner & Ritter, 1980). After its introduction in Europe, the trading of honey bee colonies between regions became restricted and subjected to strict regulations (European Commission, 1992), so to limit the spread of the parasite across countries. Our results match these previous reports of introduction as the low genetic differentiation between mites infesting French and Dutch honey bees suggests that these populations derived from a common origin, while more elevated levels of differentiation between these two locations and Sweden suggest a distinct introduction. However, the level of differentiation between mites from these three regions compared to Norwegian mites was much higher, pointing to a third introduction event. This hypothesis matches reports on the arrival of V. destructor in Norway, which was first found in the spring of 1993 in Oslo, while the natural spread of the mite was still confined to the southernmost part of Sweden. V.
destructor was in fact found for the first time in Sweden on the island of Gotland in 1987 (Fries, 1987) and later in Malmö, on the southern Swedish mainland, in the spring of 1991 (Fries, 1991  Notably, the distribution of MLGs significantly differed across the populations studied, and strong and significant differences were also found when performing pairwise allelic differentiation analyses between regions. Altogether, these results suggest that the mites are adapting to their local host populations, despite the relatively recent genetic bottlenecks caused by both the host shift and introduction of the mite in Europe (Solignac et al. 2005).
In addition to the diversification of V. destructor across Europe, the comparison between mites infesting sympatric adapted and treated host colonies revealed intriguing patterns of genetic structure. Although the specific mechanisms of horizontal large-scale transmission of the parasite currently remain unknown, the mite can easily spread within and between apiaries of a given region (Frey et al., 2011;Fries & Camazine, 2001). In molecular terms, this high transmission may prevent genetic isolation of mites infesting honey bee colonies within and across apiaries (Beaurepaire et al., 2015, F I G U R E 4 Correlation between geographic distance (GGD) and genetic distance (D est ). Results of the Mantel analysis of the spatial and genetic distance separating the mite groups. A weak but significant correlation (p =.001) was found between the distances 2017; Dynes et al., 2017). Despite this, we here observed significantly different allelic patterns and MLG distribution between V. destructor samples infesting adapted and treated honey bee colonies in all the four regions studied. These differences may be explained by numerous factors. First, neutral processes such as genetic drift may cause isolated mite populations to diverge (Freeland et al., 2011).
However, the differences reported here do not seem to result entirely from this factor as the patterns of genetic differentiation have been found to be consistently dependent on mite-host association in every location (Figures 2 and 5). Although a weak isolation by distance at the continental scale was found, the patterns of genotypic differences between V. destructor infesting adapted and treated populations did not vary according to the distance separating the groups in a given region. For example, differences remained significant even in locations where the two host populations were located at the very same apiary (i.e., Lelystad, the Netherlands, and Avignon, France).
Second, acaricide treatments of the susceptible colonies could have affected the population structure of the mites. The application of pesticides is known to lead to the development of resistance in pest populations (Georghiou, 1972), V. destructor being no exception (Martin, 2004;Milani, 1999;Spreafico et al., 2001). Although we here did not test directly for the presence of acaricide resistance, the variability of microsatellites markers in response to pesticidedriven population genetic changes can provide indirect evidence for changes in genotypic diversity and structure caused by pesticides (de Meeûs et al., 2007;Osakabe et al., 2009;Pascual-Ruiz et al., 2014).
Notably, the use of acaricides should reduce population sizes and diversity levels. As a consequence, an increase of the level of genetic divergence between treated and nontreated populations may also TA B L E 5 Pairwise population divergence indexes (D est )

Norway Tr Hilton
Ad Sørumtangen 0.009*** occur Uesugi et al., 2009). Here, the number of alleles and the levels of heterozygosity did not differ significantly between mites from adapted and treated colonies, but the diversity of mite MLGs and/or the amplitude of genetic divergence varied between these groups across populations. Notably, three of the four treated populations examined in this study (Netherlands, Sweden, and Norway) have been regularly treated with oxalic acid to control V. destructor infestations. This organic acaricide has high and consistent efficacy (Gregorc & Planinc, 2001) and works through contact by killing mites by means of high acidity (Nanetti, 1999). Given this very general mode of action, oxalic acid is not expected to select for any particular lineage of mites, as it has been demonstrated in a previous study, where even after a repeated and prolonged exposure to this compound, mites remained susceptible (Maggi et al., 2017). In contrast, the synthetic acaricide (i.e., Amitraz) used to treat the susceptible colonies in Avignon and Solérieux can foster the development of resistance in V. destructor (Kamler et al., 2016) Conte et al., 2007, 2. Panziera et al., 2017, 3. Locke, 2016, 4. Beaurepaire et al., 2019, and 5. Oddie et al., 2017, while the percentages represent the proportion of individuals out of the 863 mites included in this analysis.

F I G U R E 5
Comparison of MLG diversity (Shannon Index). Graph representing 95% confidence interval of the Shannon index, illustrating the diversity of multilocus genotypes (MLGs) across mite groups. The diversity analysis was conducted after excluding all mites with missing data. The sample size after this exclusion is reported for every group. Overall, contrasting patterns of MLG diversity could be observed between adapted and treated populations across the different locations. Significantly higher levels of MLG diversity were found in mites infesting adapted host populations in Avignon (France) and Tiengemeten (the Netherlands) compared to the mite infesting treated hosts from their respective region. Notably, for Norwegian mites, a significantly lower level of MLG diversity was found in the adapted group compared to treated one, whereas for Swedish mites, no significant difference was found Another factor that may explain the results documented here lies in the natural adaptations of honey bee colonies to V. destructor.
Most strikingly, the pairwise comparisons revealed that mites from adapted colonies from France and one location of the Netherlands (Tiengemeten) were not significantly genetically different, while the mite infesting susceptible colonies in these two regions were.
This may suggest parallel evolution similar to their honey bee hosts (Locke et al., 2012;Oddie et al., 2018). On the other hand, the patterns of change in the diversity of MLGs between mites infesting the adapted and treated host populations were not consistent across regions. For instance, in the populations located in France and the Netherlands, a higher diversity of MLGs was observed in the adapted colonies compared to local treated ones. In contrast, the level of MLG diversity was higher in the treated colonies from Norway and was not significantly different between the three groups located in Sweden. Possibly, in some surviving populations, host-mediated selection may promote mite genotypes expressing a decreased level of reproduction (i.e., selection for lower parasite virulence; Seeley, 2007), while in others, the hosts may select for mite genotypes having specific chemical mimicry abilities (Kather et al., 2015;Le Conte et al., 2015). These results suggest that different selective forces may be acting on the various mite groups, thereby representing hot spots and cold spots of evolution as postulated by the geographic mosaic of coevolution theory (Thompson, 2005). Although the particular host traits shaping the population structure of V. destructor remain to be discovered, these results confirm previous findings documenting significant temporal changes of population structure between parasites sampled in adapted and susceptible honey bee colonies over nine years (Beaurepaire et al., 2019). Altogether, these results provide empirical evidence that honey bee selective pressure influences V. destructor population structure, as previously observed in other systems tem (Eliash & Mikheyev, 2020). Half a century after its introduction into Europe, the ongoing genetic diversification of V. destructor in Europe illustrates well the evolutionary potential of parasites and represents a prominent example of ongoing coevolution between hosts and parasites.

ACK N OWLED G M ENTS
The authors of this work wish to express their sincere gratitude to COLOSS for providing the platform for the conception and realization of this project and to the beekeepers maintaining the experimental colonies. This study was financially supported by the Ricola

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The dataset used in this study can be found here: https://datad ryad.