<figure><img src="image_1.png" /></figure>

HAL
open science

# Population genetics of ectoparasitic mites suggest arms race with honeybee hosts

Alexis Beaurepaire, Arrigo Moro, Fanny Mondet, Yves Le Conte, Peter

Neumann, Barbara Locke

► To cite this version:

Alexis Beaurepaire, Arrigo Moro, Fanny Mondet, Yves Le Conte, Peter Neumann, et al.. Population genetics of ectoparasitic mites suggest arms race with honeybee hosts. Scientific Reports, 2019, 9 (1), 10.1038/s41598-019-47801-5. [https://doi.org/10.1038/s41598-019-47801-5](https://doi.org/10.1038/s41598-019-47801-5)

HAL Id: hal-02626257

https://hal.inrae.fr/hal-02626257

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# Population genetics of ectoparasitic mites suggest arms race with honeybee hosts

Received: 19 February 2019

Accepted: 22 July 2019

Published online: 06 August 2019

Alexis L. Beaurepaire<sup>D 1,2</sup>, Arrigo Moro<sup>D 2</sup>, Fanny Mondet<sup>1</sup>, Yves Le Conte<sup>D 1</sup>, Peter Neumann<sup>2,3</sup> & Barbara Locke<sup>D 4</sup>

The ectoparasitic mite, *Varroa destructor*, is the most severe biotic threat to honeybees (*Apis mellifera*) globally, usually causing colony death within a few years without treatments. While it is known that a few *A. mellifera* populations survive mite infestations by means of natural selection, the possible role of mite adaptations remains unclear. To investigate potential changes in mite populations in response to host adaptations, the genetic structure of *V. destructor* in the mite-resistant *A. mellifera* population on Gotland, Sweden, was studied. Spatio-temporal genetic changes were assessed by comparing mites collected in these colonies, as well as from neighboring mite-susceptible colonies, in historic (2009) and current (2017/2018) samples. The results show significant changes in the genetic structure of the mite populations during the time frame of this study. These changes were more pronounced in the *V. destructor* population infesting the mite-resistant honeybee colonies than in the mite-susceptible colonies. These results suggest that *V. destructor* populations are reciprocating, in a coevolutionary arms race, to the selection pressure induced by their honeybee host. Our data reveal exciting new insights into host-parasite interactions between *A. mellifera* and its major parasite.

Anthropogenic movement of species can threaten biodiversity, agriculture, ecosystem functioning and can facilitate the spread of harmful pathogens<sup>1,2</sup>. The western honeybee, *Apis mellifera*, is a perfect example of a species that has experienced a rapid expansion in geographical distribution due to international trade and globalization and in this process has acquired many novel parasites and pathogens<sup>3–5</sup>.

The most dramatic consequence of the global spread of *A. mellifera* is the propagation of the invasive ectoparasitic mite, *Varroa destructor*. This mite is arguably the most severe threat to *A. mellifera* globally, practically exterminating wild colonies and severely affecting the management and profitability of beekeeping in the wake of its global spread during the 1980's and 1990's<sup>6</sup>. The damage this parasite causes to its new host by feeding on adults and brood is amplified by the multiple viruses it carries and transmits<sup>7–9</sup>.

*V. destructor*, originally restricted to the Asian continent, has a non-lethal relationship with its natural host, the eastern honey bee, *Apis cerana*<sup>10</sup>. However, due to the transportation and introduction of *A. mellifera* in Asia, the mite managed to switch host and has successfully established itself as a harmful parasite in *A. mellifera* honeybee colonies<sup>11</sup>. Without a long-term coadaptive evolution, as is shared between the mite and its natural host, *A. mellifera* colonies are ill-equipped to cope with this new invasive parasite.

Current treatment strategies used in apiculture for the control of *V. destructor* infestations are costly, time-consuming, can harm the host and alter the quality of bee products<sup>6</sup>. Despite these drawbacks, most colonies of *A. mellifera* require treatment to have a chance for surviving the parasite infestation for more than 1–2 years. Yet, a few populations of *A. mellifera* exhibit resistance or tolerance traits that allow them to survive extended periods without treatments<sup>12</sup>.

To understand how these colonies deal with the invasive parasite, surviving honeybee populations have been studied extensively in the past decades. This research has highlighted that a wide range of individual or colony-level mechanisms are involved in their survival and that many of these traits are inheritable<sup>13–17</sup>. One such well-studied isolated honeybee population on the island of Gotland, Sweden, has been living treatment-free for almost two decades<sup>12</sup>. This population was established in the late 1990's as an isolated natural selection

<sup>1</sup>INRA, UR 406 Abeilles et Environnement, Avignon, France. <sup>2</sup>Vetsuisse Faculty/University of Bern, Institute of Bee Health, Bern, Switzerland. <sup>3</sup>Agroscope, Swiss Bee Research Center, Bern, Switzerland. <sup>4</sup>Department of Ecology, Swedish University of Agricultural Sciences, Uppsala, Sweden. Correspondence and requests for materials should be addressed to A.L.B. (email: alexis.beurepaire@vetsuisse.unibe.ch)
<table><caption>**Table 1.** Information on the samples. Table indicating the group names (Groups), period of sampling (Period), GPS coordinates (GPS coord.), number of hives (N<sub>Hives</sub>) and number of individuals for each group (N<sub>Ind</sub>) used in this study.</caption><thead><tr><th>Groups</th><th>Period</th><th>GPS coord.</th><th>N<sub>Hives</sub></th><th>N<sub>Ind</sub></th></tr></thead><tbody><tr><td>Resistant</td><td rowspan="2">Historic 2009</td><td>57°4′7.3″N 18°12′27.0″E</td><td>8</td><td>146</td></tr><tr><td>Susceptible</td><td>57°8′9.4″N 18°18′46.0″E</td><td>11</td><td>108</td></tr><tr><td>Resistant</td><td rowspan="3">Current 2017–18</td><td>57°4′7.3″N 18°12′27.0″E</td><td>3</td><td>54</td></tr><tr><td>Susceptible</td><td>57°22′27.0″N 18°40′24.3″E</td><td>4</td><td>41</td></tr><tr><td>Mainland</td><td>59°49′4.9″N 17°39′22.9″E</td><td>4</td><td>83</td></tr></tbody></table>

experiment, allowing *A. mellifera* bees to naturally adapt to *V. destructor* parasitism¹⁸. This experiment resulted relatively quickly in improved survival rates through naturally adapted traits that limit the reproductive success of *V. destructor* and are genetically inheritable within the population¹⁶,¹⁸,¹⁹.

Host-parasite coevolutionary theory predicts that adapted host resistance traits, in response to parasitic pressure, are expected to induce a reciprocated selection pressure on the parasite, strong enough to drive counter adaptations towards a fitness optima²⁰,²¹. This ultimately results in an arms race with a series of adaptations and counter-adaptations between the host and the parasite²². In most scenarios, the parasite, usually with a shorter generation time than its host, would have an advantage in this arms race²³,²⁴. However, in this particular system, *V. destructor* is at a disadvantage due to an important factor: the invasive mite population in Europe has a low level of genetic diversity²⁵. This is in part due to a founder effect abolishing most genetic diversity during only a few invasion occasions from its original host, but also due to the reproductive biology of the mite with frequent inces- tuous mating²⁵. Within the sealed worker brood cells of developing *A. mellifera* pupae, the mother mite produces a single haploid male offspring followed by 4–5 diploid females²⁶. The adult male copulates multiple times with its adult sisters before the parasitized bee emerges as an adult. At which point, the male dies in the cell and the fertilized daughters along with their mother will enter the honeybee colony’s mite population growth cycle with on average 10–15 generations per year⁶. This incestuous reproductive system generates high inbreeding levels. Nevertheless, occasionally two mites may enter a brood cell together and mating can then occur between lineages introducing genetic admixing²⁷.

To date, it is unclear if or how mites reciprocate with antagonistic adaptations to their adapted host. Antagonistic coevolution (such as that between hosts and parasites) is expected to drive molecular evolution leading to genetic divergence between populations at a much faster rate compared to selection pressures of envi- ronmental change²⁸,²⁹. The aim of this study was therefore to analyse temporal genetic structure changes of *V. destructor* mites infesting the naturally adapted mite-resistant honeybee population on Gotland, Sweden, to investigate whether the parasite population shows signs of reciprocated adaptations in response to the resistant adaptations of its host. The mites from the mite-resistant population on Gotland were compared with a geograph- ically neighboring mite population from managed mite-susceptible honeybees that have not experienced natural selection pressures. The changes in the genetic diversity and genetic structure was compared over time using microsatellites to compare historic mite samples collected in 2009 with current mite samples collected from the same populations in 2017 and 2018.

## Results

A total of 432 *V. destructor* adult females were genotyped at nine polymorphic microsatellite markers to study the temporal genetic structure of the mites in the mite-resistant colonies and local susceptible colonies from 2009 (“historic samples”) and 2017–18 (“current samples”) located on the island of Gotland, Sweden (Table 1). In addi- tion, we used samples from the mainland from the apiary of the University of Uppsala to control the degree of isolation of the Gotland mite populations. Overall, a low number of alleles was found for all markers (N<sub>A</sub> = 2–3, Table S1). An analysis of linkage disequilibrium revealed that none of the marker pair was significantly linked after Bonferroni corrections for multiple testing (all *p*-values > 0.00035), confirming the independence of the markers used.

The comparison of genetic diversity showed that the number of alleles and the level of heterozygosity did not differ significantly among the five *V. destructor* groups (Kruskall Wallis, *p* > 0.05) (Fig. S1). In addition, all five groups significantly deviated from the Hardy-Weinberg equilibrium (*p* < 0.001). Notably, a few private alleles were detected in all groups, but these alleles had a very low frequency (<3%). Finally, a rarefaction analysis to assess the impact of sample size on the estimates of genetic diversity showed that the level of allelic richness reaches a plateau after about 20 mites are included, irrespective of the group considered (Fig. S2).

To assess the level of genetic differentiation between the mite groups, we calculated two distinct pairwise genetic differentiation indexes (F<sub>ST</sub> and D<sub>est</sub>). These tests indicated different patterns between the groups com- pared (Table 2, Fig. 1). Minute and non-significant levels of genetic differentiation were found between mites coming from naturally-surviving and susceptible colonies of the historic collection (F<sub>ST</sub> = 0.002 and D<sub>est</sub> = 0.000, *p* = 0.720). In contrast, the level of differentiation between these two groups was more marked and significant in the current samples (F<sub>ST</sub> = 0.046 and D<sub>est</sub> = 0.015, *p* < 0.01). In parallel, when comparing mites from the same host population between the two sampling dates, more elevated and significant levels of population differ- entiation were obtained. Notably, this difference was about three times higher in the samples from the resistant colonies (F<sub>ST</sub> = 0.178 and D<sub>est</sub> = 0.069, *p* < 0.001) compared to the susceptible ones (F<sub>ST</sub> = 0.067 and D<sub>est</sub> = 0.022,
<table><caption>**Table 2.** Results of the pairwise analysis of genetic differentiation. Table indicating the results of the pairwise population differentiation analysis, showing the groups compared (Groups), the level of genetic differentiation using two distinct estimates (F<sub>ST</sub> and D<sub>est</sub>) and the associated *p*-values based on 9999 permutations.</caption><thead><tr><th>Groups</th><th></th><th>F<sub>ST</sub></th><th>D<sub>est</sub></th><th>p</th></tr></thead><tbody><tr><td>S-historic</td><td>R-historic</td><td>0.002</td><td>0.000</td><td>0.720</td></tr><tr><td>S-current</td><td>R-current</td><td>0.046</td><td>0.015</td><td>&lt;0.01</td></tr><tr><td>S-historic</td><td>S-current</td><td>0.067</td><td>0.022</td><td>&lt;0.01</td></tr><tr><td>R-historic</td><td>R-current</td><td>0.178</td><td>0.069</td><td>&lt;0.001</td></tr><tr><td>S-current</td><td>Mainland</td><td>0.069</td><td>0.031</td><td>&lt;0.001</td></tr><tr><td>R-current</td><td>Mainland</td><td>0.049</td><td>0.020</td><td>&lt;0.001</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Genetic differentiation across populations of *Varroa destructor*. Schematic maps showing the results of the pairwise genetic differentiation among groups. Significant results are represented by full lines, non-significant by dashed lines, the size of the lines is proportional to the values of the estimates, also indicated next to the relevant lines (on top: F<sub>ST</sub> values, and D<sub>est</sub> below and between parentheses).</figcaption></figure>

*p* < 0.01). Finally, the comparison between the samples from the mainland and the two current Gotland populations indicated low but highly significant genetic differentiation (Table 2, Fig. 1).

To investigate how individual mite genotypes were represented in the five groups, we analysed the frequency and the distribution of Multi-Locus Genotypes (MLGs) in our samples. Overall, 64 MLGs were detected (Table 3). The mainland group was the most diverse, with 27 MLGs detected and 19 private MLGs (i.e. MLGs sampled only in one group). In Gotland, the number of MLGs found in the historic samples was higher (N = 33) than in the current samples (N = 22). However, the number of private MLGs was higher in the current samples from the mite-resistant group, with eight unique MLGs accounting for 34.38% of the mites of this group. The comparison of the distribution of the five most prevalent MLGs across populations and periods revealed highly significant differences across groups (Chi² test, *p* < 0.001). Interestingly, the most predominant MLG in the historical samples (>58% in both populations) was not found in the current mite-resistant group but was found in the current susceptible group (18.18%) (Fig. 2).

Finally, to assess the spatio-temporal genetic changes in the four mite groups from the island of Gotland, an analysis of population structure was conducted using the software Instruct³⁰. Here, the samples from the mainland were excluded to focus on the changes due to the host colony phenotype (resistant vs susceptible) and exclude any background noise from a potential island effect. This analysis revealed that the most likely number of genetic clusters in this mite sample was two (ΔK₂ = 269.46). Displaying the individuals according to these two clusters (Fig. 3) revealed that most mites from the two historic groups belonged to the first cluster whereas samples from the two current groups were mostly of the second cluster. Accordingly, the likelihood of the four groups
<table><caption>**Table 3.** Distribution and prevalence of Multi-Locus Genotypes (MLGs). The number (N) and % of individuals of the different categories indicated in the left column.</caption><thead><tr><th>Category</th><th></th><th>N</th><th>% of individuals</th></tr></thead><tbody><tr><td rowspan="6">Number of Samples</td><td>S-historic</td><td>93</td><td>32.07%</td></tr><tr><td>R-historic</td><td>91</td><td>31.38%</td></tr><tr><td>S-current</td><td>22</td><td>7.59%</td></tr><tr><td>R-current</td><td>32</td><td>11.03%</td></tr><tr><td>Mainland</td><td>52</td><td>17.93%</td></tr><tr><td><strong>Total</strong></td><td>290</td><td>100.00%</td></tr><tr><td rowspan="10">Number of MLGs</td><td>S-historic</td><td>19</td><td>32.07%</td></tr><tr><td>R-historic</td><td>23</td><td>31.38%</td></tr><tr><td>S-current</td><td>12</td><td>7.59%</td></tr><tr><td>R-current</td><td>15</td><td>11.03%</td></tr><tr><td>Mainland</td><td>27</td><td>17.93%</td></tr><tr><td><strong>Total</strong></td><td>45</td><td>100.00%</td></tr><tr><td>Historic</td><td>33</td><td>63.45%</td></tr><tr><td>Current</td><td>22</td><td>18.62%</td></tr><tr><td><strong>Total</strong></td><td>55</td><td>82.07%</td></tr><tr><td>Susceptible</td><td>24</td><td>39.66%</td></tr><tr><td>Resistant</td><td>34</td><td>42.41%</td></tr><tr><td><strong>Total</strong></td><td>55</td><td>82.07%</td></tr><tr><td rowspan="6">Private MLGs</td><td>S-historic</td><td>8</td><td>10.75%</td></tr><tr><td>R-historic</td><td>11</td><td>21.98%</td></tr><tr><td>S-current</td><td>3</td><td>18.18%</td></tr><tr><td>R-current</td><td>8</td><td>34.38%</td></tr><tr><td>Island-current</td><td>15</td><td>44.44%</td></tr><tr><td>Mainland</td><td>19</td><td>50.00%</td></tr><tr><td><strong>Total</strong></td><td>32</td><td>20.17%</td></tr><tr><td rowspan="6">Distribution of the most prevalent MLG</td><td>S-historic</td><td>54</td><td>58.06%</td></tr><tr><td>R-historic</td><td>53</td><td>58.24%</td></tr><tr><td>S-current</td><td>4</td><td>18.18%</td></tr><tr><td>R-current</td><td>0</td><td>0.00%</td></tr><tr><td>Mainland</td><td>0</td><td>0.00%</td></tr><tr><td><strong>Total</strong></td><td>111</td><td>38.28%</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distribution of the most prevalent MLGs across groups. Prevalence of the different *V. destructor* Multi-Locus Genotypes (MLGs) across the different groups. The coloured bars represent the five most prevalent MLGs and the white bars represent the rest of the MLGs. The letters represent the populations (S: Susceptible, R: resistant, M: mainland) and dates (H: historic and C: current).</figcaption></figure>

to be affiliated with the two genetic clusters was significantly different (χ² = 142.40, p < 0.001). In addition, when comparing the mites collected in the mite-resistant and the mite-susceptible colonies, significant differences were found in the current samples (χ² = 9.22, p < 0.01) but not in the historic ones (χ² = 0.25, p > 0.05).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Genetic Structure of the Gotland mite populations. The blue and yellow colours indicate the two genetic clusters, the X-axis represents individuals and the Y-axis represents the probability to belong to the two different clusters ($\Delta K_2$ = 269.46).</figcaption></figure>

Discussion

The results of this study show that the genetic structure of *V. destructor* populations in the island of Gotland changed significantly between the time when our historic and current samples were collected, with differing degrees between the resistant and susceptible colonies. In addition, the distribution of the mite genotypes sampled in 2017 in the colonies from the resistant population suggests that genetic divergence in the parasite population of Gotland is ongoing as a result of selection by the host.

Despite marked differences in the number of individuals screened between the historic and the current groups, the markers used were able to accurately grasp the overall diversity in the populations studied here. For instance, when comparing the level of genetic diversity between the five *V. destructor* groups, no differences were found in the number of alleles or in the level of heterozygosity. Moreover, the rarefaction analysis performed on the five groups of mites did not reveal any difference between the groups with more samples (i.e. the historic groups) and the groups with less (i.e. the current groups), with a plateau phase being reached after about 20 individuals are analysed. This result is in line with the findings of another study using the same method³¹. Here, the lowest number of samples analysed in a group was well above this threshold (N = 41 in the S-current group). Additionally, populations of *V. destructor* are homogenous at the apiary level³¹,³². Therefore, clearly neither the difference in the number of colonies nor the difference in sample size used across the five groups had any significant impact on the genetic diversity estimates in our results.

Low levels of genetic diversity have initially been reported from a large number of *V. destructor* populations infesting *A. mellifera* throughout the world, suggesting that strong bottlenecks have taken place after this parasite's host switch and further dispersal of the infested colonies²⁵. Yet, more recent studies have revealed that the diversity of *V. destructor* is not as limited as previously thought in populations of Europe²⁷, the USA³¹,³² and Asia³³,³⁴. The discrepancies between the earlier work on the varroa mite's population genetics and the more recent investigations could be due to (i) the average low number of samples of the populations initially screened by Solignac *et al.* (2005) (avg = 12.44 mites per location), (ii) the degree of polymorphism of the markers screened, and/or (iii) an increased diversification of the parasite populations over the ten years separating the studies. Discriminating between the first two possibilities is not possible given the lack of investigations using the same markers as Solignac *et al.* (2005). Nevertheless, the results of this study suggest that the possibility (iii) is not a valid one, since we did not detect significant changes in the genetic diversity between our historic and current sample groups.

Despite this lack of differences in the level of genetic diversity over time, significant changes in population structure between the historic and current samples were detected. The low number and low frequency of private alleles detected demonstrate that these changes in population structure were mainly driven by variation in the frequency of the alleles initially found in the mite populations and not because of the addition or deletion of alleles. The results of the genetic structure analysis using the software Instruct did not show a pure segregation of the two genetic clusters across locations or dates, which can be explained by the fact that a unique lineage of *V. destructor* invaded *A. mellifera* colonies of Europe a few decades ago only⁷,²⁵,³⁵. Yet, despite this unique origin, significant genetic variation can now be observed, and the number of individuals belonging to the two genetic clusters detected here differed significantly across time and across the resistant and susceptible colonies in the current sampling period but not in the historic one. This temporal population structure was confirmed by the analysis of Multi-Locus Genotypes. The latter analysis permitted to assess the diversity and structure of *V. destructor* populations in a subtler way. Moreover, the differences of magnitude between the two indexes of population differentiation ($F_{ST}$ and $D_{est}$) reported also reflect these changes.

Hence, these independent but complementary analyses clearly indicate that *V. destructor* populations have changed within the eight years separating the collection of the current and historical samples. These changes led to a higher level of differentiation and an increased diversity of mite genotypes in the resistant colonies compared to the susceptible ones. These variations of population structure and diversity may reflect host-parasite interactions between *A. mellifera* and *V. destrutor*, but could also be caused by an influx of foreign mite genotypes in the current samples, the development of resistance to acaricide treatments and/or genetic drift. First, the comparison of the mainland mite population with the two populations located on Gotland suggest that an influx of mites was not responsible for the differences detected. Despite the common origin of the mites infesting both locations, the significant barrier to gene flow observed between the island and the mainland and the high proportion of private MLGs in the mainland and island mites in the current sampling period (50% and 44.44%, respectively) indicate that these populations are isolated. Second, both susceptible populations (in Gotland and Uppsala) are treated with oxalic acid, a compound that does not select specific resistant mites³⁶. Thus, the changes of genetic structure
detected is not likely be caused by the treatment of the susceptible colonies. Lastly, genetic drift, a process causing the allele frequencies of populations to vary from one generation to the next as a result of chance³⁷, is probably not responsible for the integrality of the changes in population structure we observed within the time frame of our study. With a strong genetic drift effect, a reduction of diversity is expected³⁷,³⁸. In this study, there was no significant difference in the number of alleles and heterozysogity levels between the different mite populations from Gotland and from the mainland. In contrast, we found an increased diversity of MLGs in the resistant colonies of Gotland between 2009 and 2018. Moreover, if genetic drift was responsible for the increased differentiation between the two sampling dates, the most prevalent MLGs in the historical sampling period would be more likely to prevail in the current period³⁸. Instead, the patterns of changes of mite MLGs between the two dates show that the most common MLG, representing more than half the mites sampled during the historical collection, was not sampled in the resistant group later on, but was still infesting the susceptible population. This observation suggests that a strong selection against the most common mite genotypes is occurring in the resistant colonies. To conclude, the above-mentioned evidences strongly suggest that the main factor responsible for the temporal changes in the genetic structure of *V. destructor* measured here are caused by the strong selection pressures induced through the co-evolution of the mite with its host, rapidly leading to observable genetic changes in the parasite population, potentially aiding their survival in adapted mite-resistant *A. mellifera* colonies.

The notion of *V. destructor* adaptation in response to its host or environment is not far fetched. Examples such as the ability of the mite to adapt to acaridice treatments²⁷,³⁹ show how quickly this parasite's populations can react to selection pressures, despite it's apparent lack of genetic diversity. In our study, given an average of ten generations of mites per year in the Northern-European climate of Gotland⁴⁰, the time separating the collection of samples from the historic and current group is more than enough to assess the impact of strong selective forces guided by antagonistic coevolution²⁸,²⁹. *V. destructor* genetic variation has been linked to phenotypic variation associated with fertility and virulence differing between the two parasite haplotypes that spread globally from Japan and Korea³⁵,⁴¹. When *V. destructor* was first reported in Brazil, the fertility and the virulence of the parasite were low, causing little damage to *A. mellifera* colonies⁴². However, only a few years later, the virulence of *V. destructor* had increased in the region similar to what was experienced in European populations at the time. Genetic analyses demonstrated that the initial avirulent Japanese *V. destructor* haplotype had been replaced by the more virulent Korean haplotype⁴³. The fact that variability in the consequences of mite infestation so clearly exists within the *V. destructor* haplotypes, having a significant effect on colony survival, suggests that disregarding the possibility of, for example, reduced virulence adaption of the mite influencing host survival, is no longer reasonable. Indeed reduced virulence in mites has been proposed previously to explain the long term survival of another isolated honeybee population in the Arnot Forest in Ithaca, NY⁴⁴.

In 2006, Fries and Bommarco performed a cross-infection experiment with the Gotland resistant honeybee population to test for varying host responses to mites, sourced from either the resistant or susceptible populations⁴⁵. In their study, only 3 years before the historical samples of this study were taken, they found that mite source had no effect on the colonies and that Gotland colonies had significantly reduced mite infestation rates. Since their study was published, there has been a consensus in the scientific community to regard the mite as a "fixed factor" in the studies of underlying mechanisms to explain the long-term survival of naturally adapted honeybee populations, supported by the limited genetic diversity in the mite population²⁵. The historic mite samples of this study confirm the results of Fries and Bommarco's in 2006 that the mite source did not influence the bees, since the mites in 2009, only 3 years later, did not significantly differ genetically between groups. However, by looking over a decade after Fries and Bommarco's study it is evident that a change has occurred in the genetic structure of the mite populations between surviving and treated colonies. Performing a similar cross-infection experiment testing the current mite genotype in the Gotland resistant and susceptible colonies will show whether the genotypic differences observed in the current samples of this investigation are associated with phenotypic differences and will help understanding how the diversity of *V. destructor* populations observed at the genetic level may impact the survivability of their hosts. If the genetic changes on the neutral markers of this study are confirmed by phenotypic changes, the recently updated version of the genome of the mite⁴⁶ could be used to unravel the genetic bases of the parasite adaptations. Inhibiting mite reproductive success is a well defined genetically inheritable trait of the Gotland mite-resistant honeybee population, even if the mechanisms explaining how the bees are capable of this are still not completely clear¹⁶. Whether this host trait could explain the differences in mite genotype distribution we observed here remains to be studied.

In conclusion, the observed changes over time in the genetic structure of *V. destructor* suggest adaptations of the parasite, in a host-parasite coevolutionary arms race, most likely in response to selection pressures applied by the adapted resistant traits of the host. The magnitude of these changes between the historic and current samples of this study demonstrate a relatively fast response. Model systems for experimental host-parasite coevolution research has been largely restricted to microbes and short generation hosts⁴⁷. The natural selection experiment on Gotland has been, and continues to provide, a unique opportunity to study host-parasite coevolution in real time between two larger eukaryote organisms.

## Material and Methods

**Sampling.** In summer 2009, *V. destructor* adult females were collected from the brood of eight honeybee colonies in the naturally adapted mite-resistant honeybee population on the island of Gotland Sweden as well as from eleven mite-susceptible colonies of *A. mellifera* that require regular mite-treatment with oxalic acid to survive, located approximately ten km away from the resistant colonies (Table 1, Fig. 1). In the summer of 2017, mites were collected again from three colonies of the same population of mite-resistant *A. mellifera* colonies. However, control mites were collected in 2018, because parasite loads in treated colonies were very low in 2017 (most likely due to the regular treatments applied by beekeepers). In addition to the samples collected on Gotland, mites were collected on the mainland in Sweden from four mite-susceptible honeybee colonies at the Swedish University
of Agricultural Sciences in Uppsala during the summer of 2017 as a secondary control group to account for a possible island effect. All mite samples were kept at -20 °C until DNA extraction. From hereinafter, the 2017 and 2018 samples will be referred to as the temporally "current" samples while the 2009 samples are referred to as the "historic" samples for streamlining the presentation and discussion of the data (Table 1).

**Genotyping.** Total mite DNA was extracted using a Chelex protocol⁴⁸. Mites were genotyped using a set of microsatellite markers. In all, 20 markers were tested to determine optimal PCR conditions (ie multiplexing) and polymorphism on a subset of 24 individuals from three colonies of the susceptible historic group. From these 20 markers, nine were polymorphic and worked well under our lab conditions (Table S1). PCR products were sent to Genoscreen (Lille, France) to run on a 3730XL sequencer (Applied Biosystems ®). Once received, the genotypes were scored manually using Peak Scanner v 1.0 (Applied Biosystems ®).

**Analyses.** For the following analyses, the mite samples from Gotland were grouped according to their population of origin (resistant or susceptible) and the date of sampling (historic or current). The mite samples from Uppsala were included in the analyses and labeled as "mainland". Only individuals with genotypes available for five or more loci were kept in this dataset. Along these lines, a total of 432 mites were kept in the analyses (Table 1).

First, the linkage disequilibrium between pairs of markers was tested along with the Hardy-Weinberg equilibrium for each group using Fstat v 2.9.3⁴⁹. Then, to compare the genetic diversity across the five groups of mites, the number of alleles ($N_A$) and the level of observed heterozygosity (Ho) were calculated using the Microsatellite Toolkit⁵⁰ and compared using Kruskal-Wallis tests in R v. 3.5.2⁵¹. We used the software ADZE⁵² to conduct a rafefaction analysis on the allelic richness over markers for every mite groups to estimate whether the lowest sample size we used was sufficient to grasp accurately the genetic diversity of the group considered.

The software GenAlex v. 6.5³³ was then used to assess genetic differentiation in the populations of *V. destructor* of the four groups from Gotland and between the current mainland and Gotland samples by calculating two distinct estimates of population differentiation ($F_{ST}$ and $D_{est}$) as recommended by⁵⁴. $F_{ST}$ reflects the probability that two alleles drawn at random from within a group are identical by descent⁵⁵ and is commonly used to compare the genetic differentiation between two or more group of samples⁵⁴,⁵⁶. Generally, $F_{ST}$ values below 0.05 suggest low genetic differentiation across groups, whereas when this estimates ranges between 0.05 to 0.25, moderate genetic differentiation is suggested³⁷. In parallel, $D_{est}$ complements this estimate and better reflects differences in allele frequencies when markers are highly polymorphic⁵⁶,⁵⁷.

In addition, to explore the distribution of *V. destructor* genotypes, the package poppr⁵⁸ for R v. 3.5.2⁵¹ was used to compare the prevalence of Multi-Locus Genotypes (MLGs) in the five mite groups. All individuals with missing data were excluded from this analysis (Table 3).

Finally, an analysis of genetic structure without a priori was performed using the software Instruct³⁰ in order to infer the optimal number of genetic clusters underlying the mite sample (called "K") from Gotland (parameters set as following: admixture model, 20 chains of each K varying from 1 to 8, 50'000 burn-in and 100'000 iterations). The methods from Evanno *et al.* (2005) was then used to determine the most likely K in the overall dataset. The results of the Instruct analysis were pooled using the software CLUMPP⁶⁰ and a figure corresponding to the most likely number of K was prepared with the software Distruct⁶¹. In order to test whether the probabilities of belonging to a cluster differed among the different groups, a Chi-squared test was used to evaluate the relationship between the inferred clusters and the mite groups. To do so, the number of individuals in the different genetic clusters for each group was obtained from the ancestry values of the populations obtained with the CLUMPP software and was compared using a contingency table. In parallel, the differences between naturally-surviving and treated colonies for each sampling date was tested using the same method⁵⁹,⁶²,⁶³.

**Data Availability**

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

**References**

1. Banks, N. C., Paini, D. R., Bayliss, K. L. & Hodda, M. The role of global trade and transport network topology in the human-mediated dispersal of alien species. *Ecol. Lett.* **18**, 188-199 (2015).

2. Tabak, M. A., Piaggio, A. J., Miller, R. S., Sweitzer, R. A. & Ernest, H. B. Anthropogenic factors predict movement of an invasive species. *Ecosphere* **8**, e01844 (2017).

3. Moritz, R. F. A., Härtel, S. & Neumann, P. Global invasions of the western honeybee (*Apis mellifera*) and the consequences for biodiversity. *Ecoscience* **12**, 289-301 (2005).

4. Evans, J. D. & Schwarz, R. S. Bees brought to their knees: Microbes affecting honey bee health. *Trends Microbiol.* **19**, 614-620 (2011).

5. Oldroyd, B. P. Coevolution while you wait: *Varroa jacobsoni*, a new parasite of western honeybees. *Trends Ecol. Evol.* **14**, 312-315 (1999).

6. Rosenkranz, P., Aumeier, P. & Ziegelmann, B. Biology and control of *Varroa destructor*. *J. Invertebr. Pathol.* **103**, S96-S119 (2010).

7. Wilfert, L. *et al*. Honeybee disease: Deformed wing virus is a recent global epidemic in honeybees driven by *Varroa* mites. *Science* **351**, 594-597 (2016).

8. Levin, S., Sela, N. & Chejanovsky, N. Two novel viruses associated with the *Apis mellifera* pathogenic mite *Varroa destructor*. *Sci. Rep.* **6**, 37710 (2016).

9. Tentcheva, D. *et al*. Prevalence and Seasonal Variations of Six Bee Viruses in *Apis mellifera* L. and *Varroa destructor* Mite Populations in France. *Appl. Environ. Microbiol.* **70**, 7185-7191 (2004).

10. de Jong, D., Morse, R. A. & Eickwort, G. C. Mite Pests of Honey Bees. *Annu. Rev. Entomol.* **27**, 229-252 (1982).

11. Anderson, D. L. & Trueman, J. W. H. *Varroa jacobsoni* (Acari: Varroidae) is more than one species. *Exp. Appl. Acarol.* **24**, 165-189 (2000).

12. Locke, B. Natural *Varroa* mite-surviving *Apis mellifera* honeybee populations. *Apidologie* 1-16, https://doi.org/10.1007/s13592-015-0412-8 (2016).
13. Rinderer, T. E., Harris, J. W., Hunt, G. J. & de Guzman, L. I. Breeding for resistance to *Varroa destructor* in North America. *Apidologie* **41**, 409–424 (2010).

14. Büchler, R., Berg, S. & Le Conte, Y. Breeding for resistance to *Varroa destructor* in Europe. *Apidologie* **41**, 393–408 (2010).

15. Locke, B., Le Conte, Y., Crauser, D. & Fries, I. Host adaptations reduce the reproductive success of *Varroa destructor* in two distinct european honey bee populations. *Ecol. Evol.* **2**, 1144–1150 (2012).

16. Locke, B. Inheritance of reduced *Varroa* mite reproductive success in reciprocal crosses of mite-resistant and mite-susceptible honey bees (*Apis mellifera*). *Apidologie* **47**, 583–588 (2016).

17. Kurze, C., Routtu, J. & Moritz, R. F. A. Parasite resistance and tolerance in honeybees at the individual and social level. *Zoology*, https://doi.org/10.1016/j.zool.2016.03.007 (2016).

18. Fries, I., Imdorf, A. & Rosenkranz, P. Survival of mite infested (*Varroa destructor*) honey bee (*Apis mellifera*) colonies in a Nordic climate. *Apidologie* **37**, 564–570 (2006).

19. Locke, B. & Fries, I. Characteristics of honey bee colonies (*Apis mellifera*) in Sweden surviving *Varroa destructor* infestation. *Apidologie*, https://doi.org/10.1007/s13592-011-0029-5 (2011).

20. Schmid-Hempel, P. *Evolutionary Parasitology*. (Oxford University Press, 2011).

21. Van Valen, L. A new evolutionary law. *Evol. theory* **1**, 1–30 (1973).

22. Thompson, J. N. *The coevolutionary process*. (University of Chicago Press, 1994).

23. Hafner, M. S. et al. Disparate Rates of Molecular Evolution in Cospeciating Hosts and Parasites. *Science* (80–). **265**, 1087–1090 (1994).

24. Schmid Hempel, P. Evolutionary parasitology. The integrated study of infections, immunology, ecology, and genetics. *Oxford Biol* (2011).

25. Solignac, M. et al. The invasive Korea and Japan types of *Varroa destructor*, ectoparasitic mites of the Western honeybee (*Apis mellifera*), are two partly isolated clones. *Proc. Biol. Sci.* **272**, 411–419 (2005).

26. Martin, S. J. Ontogenesis of the mite *Varroa jacobsoni* Oud. in worker brood of the honeybee Apis mellifera L. under natural conditions. *Exp. Appl. Acarol.* **18**, 87–100 (1994).

27. Beaurepaire, A. L., Krieger, K. J. & Moritz, R. F. A. Seasonal cycle of inbreeding and recombination of the parasitic mite *Varroa destructor* in honeybee colonies and its implications for the selection of acaricide resistance. *Infect. Genet. Evol.* **50**, 49–54 (2017).

28. Paterson, S. et al. Antagonistic coevolution accelerates molecular evolution. *Nature* **464** (2010).

29. Thrall, P. H. et al. Rapid genetic change underpins antagonistic coevolution in a natural host-pathogen metapopulation. *Ecol. Lett.* **15**, 425–435 (2012).

30. Gao, H., Williamson, S. & Bustamante, C. D. An MCMC Approach for Joint Inference of Population Structure and Inbreeding Rates from Multi-Locus Genotype Data. *Genetics* (2007).

31. Dynes, T. L. et al. Fine scale population genetic structure of *Varroa destructor*, an ectoparasitic mite of the honey bee (*Apis mellifera*). *Apidologie* **48**, 93–101 (2017).

32. Beaurepaire, A. L., Ellis, J. D., Krieger, K. J. & Moritz, R. F. A. Association of *Varroa destructor* females multiply infested cells of the honeybee Apis mellifera. *Insect Science*, https://doi.org/10.1111/1744-7917.12529 (2017).

33. Beaurepaire, A. L. et al. Host specificity in the honeybee parasitic mite, *Varroa* spp. in *Apis mellifera* and *Apis cerana*. *PLoS One* **10** (2015).

34. Navajas, M. et al. New Asian types of *Varroa destructor*: a potential new threat for world apiculture. *Apidologie* **41**, 181–193 (2010).

35. Anderson, D. L. Variation in the parasitic bee mite *Varroa jacobsoni* Oud. *Apidologie* **31**, 281–292 (2000).

36. Maggi, M. D. et al. The susceptibility of *Varroa destructor* against oxalic acid: A study case. *Bull. Insectology* **70** (2017).

37. Freeland, J. R., Kirk, H. & Petersen, S. D. *Molecular Ecology*. (Wiley-Blackwell, 2011).

38. Hamilton, M. *Population Genetics*. (Wiley-Blackwell, 2009).

39. Milani, N. The resistance of *Varroa jacobsoni* Oud. to acaricides. *Apidologie* **30**, 229–234 (1999).

40. Calis, J. N. M., Fries, I. & Ryrie, S. C. Population modelling of *Varroa jacobsoni* Oud. *Apidologie* **30**, 111–124 (1999).

41. Delfinado-Baker, M. Variability and biotypes of *Varroa jacobsoni* Oudemans. *Am. Bee J.* **128**, 567–568 (1988).

42. Rosenkranz, P. Honey bee (*Apis mellifera* L.) tolerance to *Varroa jacobsoni* Oud. In South America. *Apidologie* **30**, 159–172 (1999).

43. Garrido, C., Rosenkranz, P., Paxton, R. & Gonçalves, L. Temporal changes in *Varroa destructor* fertility and haplotype in Brazil. *Apidologie* **10**, 34 (2003).

44. Seeley, T. D. Honey bees of the Arnot Forest: a population of feral colonies persisting with *Varroa destructor* in the northeastern United States. *Apidologie* **38**, 19–29 (2006).

45. Fries, I. & Bommarco, R. Possible host-parasite adaptations in honey bees infested by *Varroa destructor* mites. *Apidologie* **38**, 525–533 (2007).

46. Techer, Maeva, A. et al. Genomic analyses of sibling honey bee ectoparasitic mite species show divergent strategies of adaptation. *bioRxiv*, https://doi.org/10.1101/512988 (2019).

47. Brockhurst, M. A. & Koskella, B. Experimental coevolution of species interactions. *Trends Ecol. Evol.* **28**, 367–375 (2013).

48. Walsh, P. S., Metzger, D. A. & Higuchi, R. Chelex 100 as a medium for simple extraction of DNA for PCR-based typing from forensic material. *Biotechnique* **v. 10**, 506–513 (1991).

49. Goudet, J. FSTAT (Version 1.2): a computer program to calculate F-statistics. *Heredity* (Edinb). **96**, 485–486 (1995).

50. Park, S. D. E. Excel Microsatellite Toolkit. (2008).

51. R Core Team. R: A language and environment for statistical computing. (2018).

52. Szpiech, Z. A., Jakobsson, M. & Rosenberg, N. A. ADZE: A rarefaction approach for counting alleles private to combinations of populations. *Bioinformatics*. https://doi.org/10.1093/bioinformatics/btn478 (2008).

53. Peakall, R. & Smouse, P. E. GenAlEx 6.5: genetic analyses in Excel. Population genetic software for teaching and research - an update. *Bioinformatics* **28**, 2537–2539 (2012).

54. Meirmans, P. G. & Hedrick, P. W. Assessing population structure: FST and related measures. *Mol. Ecol. Resour.* **11**, 5–18 (2011).

55. Wright, S. The interpretation of population structure by F-statistics with special regard to system of mating. *Evolution* (N. Y). **19**, 395–420 (1965).

56. Whitlock, M. C. G'ST and D do not replace FST. *Mol. Ecol.* **20**, 1083–1091 (2011).

57. Jost, L. GST and its relatives do not measure differentiation. *Mol. Ecol.* **17**, 4015–4026 (2008).

58. Kamvar, Z. N., Tabima, J. F. & Grünwald, N. J. *Poppr*: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. *PeerJ* **2**, e281 (2014).

59. Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. *Mol. Ecol.* **14**, 2611–2620 (2005).

60. Jakobsson, M. & Rosenberg, N. A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. *Bioinformatics* **23**, 1801–1806 (2007).

61. Rosenberg, N. A. DISTRUCT: A program for the graphical display of population structure. *Mol. Ecol. Notes* **4**, 137–138 (2004).

62. Cornman, S. R. et al. Genomic survey of the ectoparasitic mite *Varroa destructor*, a major pest of the honey bee Apis mellifera. *BMC Genomics* **11** (2010).

63. Evans, J. D. Microsatellite loci in the honey bee parasitic mite *Varroa jacobsoni*. *Mol. Ecol.* **9**, 1436–1438 (2000).
**Acknowledgements**

The authors wish to thank Åke Lyberg for his long-standing support in maintaining the Gotland mite-resistant honeybee population. Cato Östmoe is also deeply thanked for providing the samples from the mite-susceptible honeybee colonies on Gotland in 2018. The study was financially supported by the Swedish Research Council FORMAS (grant diarienr. 2016-00481), the Ricola Foundation *Nature and Culture*, the Vinetum Foundation and The Persephone Charitable and Environmental Trust. Financial support for maintaining the Gotland mite-resistant honeybee population was provided by the EC funded National Program through the Swedish Board of Agriculture.

**Author Contributions**

A.B., F.M., Y.L.C., P.N. and B.L. conceived and designed the study; B.L. collected samples; A.B. and A.M. extracted and genotyped the samples; A.B. analyzed the data; B.L. contributed reagents/materials/analysis tools; A.B., F.M., A.M., Y.L.C., P.N. and B.L. wrote the manuscript, all authors participated to and approved the submitted manuscript.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-019-47801-5.

**Competing Interests:** The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data, in the writing of the manuscript, and in the decision to publish the results.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019