Host Specificity in the Honeybee Parasitic Mite, Varroa spp. in Apis mellifera and Apis cerana

The ectoparasitic mite Varroa destructor is a major global threat to the Western honeybee Apis mellifera. This mite was originally a parasite of A. cerana in Asia but managed to spill over into colonies of A. mellifera which had been introduced to this continent for honey production. To date, only two almost clonal types of V. destructor from Korea and Japan have been detected in A. mellifera colonies. However, since both A. mellifera and A. cerana colonies are kept in close proximity throughout Asia, not only new spill overs but also spill backs of highly virulent types may be possible, with unpredictable consequences for both honeybee species. We studied the dispersal and hybridisation potential of Varroa from sympatric colonies of the two hosts in Northern Vietnam and the Philippines using mitochondrial and microsatellite DNA markers. We found a very distinct mtDNA haplotype equally invading both A. mellifera and A. cerana in the Philippines. In contrast, we observed a complete reproductive isolation of various Vietnamese Varroa populations in A. mellifera and A. cerana colonies even if kept in the same apiaries. In light of this variance in host specificity, the adaptation of the mite to its hosts seems to have generated much more genetic diversity than previously recognised and the Varroa species complex may include substantial cryptic speciation.


Introduction
The Western honeybee Apis mellifera, originally native to Europe, Africa, and the Middle East, has been repeatedly introduced in almost all regions of the world due to its importance for apiculture [1]. Introductions into Eastern Asia have been ongoing for over a century with most drastic negative consequences for global beekeeping [2]. Following its introduction, A. mellifera came into contact with a broad range of parasites and pathogens infecting native Asian we found reproducing mites in both drones and workers cells. The sampling was conducted in the apiary of the University of Los Banos (Philippines) and in Dien Bien and Son La (Vietnam) in 2013. In these locations, both honeybee host species were kept in sympatry on the same or on adjacent apiaries, within honeybee flight distance range (<1000m). Additionally, mites were collected in 2013 on the island of Cat Ba (Vietnam) a natural reserve where only A. cerana occurred. Finally, dead mites were sampled from boards placed at the bottom of three Western honeybee colonies in Lipa city, Philippines in 2015. All mites were directly transferred into 99% ethanol and stored at -20°C shortly after sampling.
The DNA of individual mites was extracted using a standard Phenol-Chloroform protocol [15]. The quality and amount of each DNA extract was determined using a Nanodrop spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, USA). All in all, 372 mites were analysed (263 from 20 A. mellifera colonies and 109 from 14 A. cerana colonies) ( Table 1).

Mitochondrial DNA analysis
A 950bp fragment of the mitochondrial Cytochrome Oxidase I gene (coxI) was amplified and sequenced using the coxI primer set (10KbCOIF1 and 6,5KbCOIR, [14]) for three mites per host species and location. The resulting sequences were trimmed using the software BIOEDIT [16] and subsequently aligned with the software MEGA V. 6.0 [17]. These fragments were compared with the NCBI database using the NCBI-BLAST tool [18] to infer which of the previously described V. destructor haplotypes best matched our samples.
The amount of divergence between all distinct haplotypes generated in this study was calculated using the software MEGA V 6.0 [17]. In parallel, a maximum likelihood tree was built using the same software. Finally, a median-joining network was constructed using the software NETWORK v. 4.6.1.2 [19].

Microsatellite DNA analyses
All Varroa were genotyped at six polymorphic microsatellite DNA loci (VD112, VD125, VD152 from [20], and VJ275, VJ292 and VJ295, from [13]) using the Fragment Profiler software V. 1.2 of the MEGABACE DNA Analysis System (GE Healthcare Life Science, Buckinghamshire, England). The number of alleles (NA), allelic richness (R) and the observed heterozygosity (H O ) were estimated for each sampling location and host species using the software Fstat V. 2.9.3 [21]. Hardy-Weinberg equilibrium tests were performed within samples for each marker using the former software.
A Principal Component Analysis was conducted on the overall microsatellite data using the R package Adegenet [22] to identify the main genetic clusters among the different mite samples based on the individual mites' genotypes. The fixation indexes (F ST ) between and within the main clusters provided by the PCA analysis were estimated using Fstat V. 2.9.3 [21]. In addition, Jost's population differentiation index (D, [23]) was estimated using the software SMOGD [24] between all locations and honeybees host species. Finally, AMOVAs were performed using the microsatellite data to identify the relevant level of Varroa genetic structuring within the PCA clusters (between locations, between colonies within locations and within colonies) using the software Arlequin V. 3.5.1.3 [25].

Mitochondrial DNA analysis
All sequences generated in this study were registered in the NCBI database under accession numbers KR528378 to KR528387 (S1 Table).
Origins of the mites from Vietnam. The mtDNA sequences of the Vietnamese samples clearly segregated according to host species. Both haplotypes were however highly similar ( Fig  1 and S1 Fig).
The sequences of the Vietnamese mites sampled on A. mellifera colonies were all identical to the previously described Korean AmK1-1 haplotype (accession GQ379056; [14]). The mites from A. cerana colonies showed more variability both within and between locations (Fig 1,  Table 2). The previously described haplotype AcV1-1 (accession GQ379061; [14]) matched the haplotype of mites sampled on A. cerana in Dien Bien and Son La. Our samples from Cat Ba were close to the AcC1-1 haplotype from Guangdong province, in Southern China (accession GQ379065; [14]). The haplotypes in Vietnam were also similar to the Japanese haplotype  (accession GQ379074.1; [14]) with only five substitutions between the haplotypes of mites sampled at Son La. This was in the same range as the variance found within the Vietnamese samples where the haplotypes of the mites sampled at Cat Ba differed by five substitutions to the Son La haplotype. The distances within the Vietnamese haplotypes were not significantly larger than those separating the Korean and the Japanese haplotype (t test, p > 0.05) but significantly larger than the haplotypes sampled on the Phillipines (t test, p < 0.001) ( Table 2).
Origins of the mites from the Philippines. The coxI sequences we obtained from the mites sampled in A. cerana and A. mellifera pupae cells in Los Banos were all identical to the Luzon 1 sp. (accession AF106894.1, [7]). In that location, the mites we sampled all shared the identical native mite haplotype irrespective of host species, unlike in Lipa city where all mites sampled in A. mellifera colonies were of the ubiquitous Korean AmK1-1 haplotype (accession number GQ379105.1; [14]).
Differences between the mites from Vietnam and Los Banos. The divergence levels between the mites we sampled in Los Banos and in Vietnam were high, with the sequences differing at the average by 3.50% ± 1.30 SD and 3.60% ± 1.30 SD from the A. mellifera and A. cerana mites from Vietnam, respectively (Table 2).

Microsatellite markers Analysis
The microsatellite markers used in this study were highly polymorphic with an average of 18.83 ± 3.33 alleles ( Table 3). None of the six markers were in Hardy-Weinberg equilibrium, due to a lack of heterozygotes which is expected as a result from obligate brother-sister mating and inbreeding in the Varroa life cycle.
Principal Component analysis. Based on the microsatellite data of all mites, the two first components obtained with the PCA explained together 38.32% of the genetic variation in our samples (first component: 30.67% and second component: 7.45%, Fig 2). When the individual mite genotypes were plotted on these two main axes, three distinct clusters could be observed, matching the three different haplotype described with the mitochondrial DNA analysis. The first one consisted of the mites sampled in A. cerana colonies in the three populations in Vietnam ("Vietnamese cluster"). A second cluster included the Varroa collected in A. mellifera in Vietnam and Lipa city ("Korean cluster"). Finally, the third cluster comprised the parasites from the colonies of A. mellifera and A. cerana from Los Banos ("Philippine cluster"). Comparison of the Genetic Diversity between the different Varroa types. The mites belonging to the Vietnamese type had a significantly higher allelic richness (14.13 ± 2.04) compared to the Varroa mites of the Korean cluster (2.84 ± 1.39) and the Philippine cluster (R = 3.63 ± 0.91, t test: p = 0.001) ( Table 1). The overall heterozygosity in the mites was low.  Genetic structuring of the Varroa haplotypes. The genetic differentiation within the clusters was low and non-significant for the mites sampled in A. cerana and A. mellifera colonies in Los Banos (F ST = 0.052, p > 0.05) ( Table 4). However, medium and significant F ST values were obtained when comparing the different sampling locations where the Korean cluster was found (F ST = 0.172, p < 0.05) and the different sampling locations where the Vietnamese cluster was found (F ST = 0.205, p < 0.05). Much higher and highly significant F ST values were found when comparing the three clusters ( Table 4).
The AMOVA revealed that the geographic location had a highly significant importance for the mites of theVietnamese cluster (15.03%, p < 0.001) and in the Korean cluster (15.12%, p < 0.01) ( Table 5). Notably, the genetic distance in the Korean cluster was lower within the two Vietnamese populations (D = 0.008) than between the mites from Lipa city and the two Vietnamese populations (D = 0.023) (S2 Table). The amount of genetic diversity varied significantly within location for both, the Vietnamese and the Korean clusters (16.72% and 15.85%, respectively, p < 0.001). Finally, the largest source of variation in these two groups resulted from the differences among mites within colonies (68.24% for the Vietnamese and 69.02% for the Korean clusters, p < 0.001). For the mites sampled in Los Banos, only this latter level was significant (80.47%, p < 0.05), but not the differences among hosts or among colonies within host.
Hybrid detection. We found no evidence of direct hybridization between mites of the two host species. Only eight Varroa mites in the whole dataset carried alleles that were also found in mites sampled in the alternative host species. These individuals were exclusively found in A. cerana colonies in Vietnam and carried one or two alleles that were otherwise specific to the Korean haplotype in A. mellifera colonies. However, none of these mites were direct hybrids, as all other microsatellite loci had private alleles of the Vietnamese haplotype. Moreover, these few shared alleles were found in the homozygotic state, suggesting that they were independent homoplasic alleles of the same length as in the A. mellifera sampled mites but not a result of hybridisation.

Discussion
In this study, we found that the Varroa mite shows different patterns of host specificity between A. cerana and A. mellifera in the Philippines and in Vietnam. Whereas the native Philippine mite was found in colonies of both host species in Los Banos, we found strong host specificity and complete reproductive isolation between the Varroa types parasitizing A. cerana and A. mellifera colonies in Vietnam. The Korean haplotype was only found in A. mellifera colonies but never found in any A. cerana colony we sampled in Vietnam and the Philippines.

Origins and diversity of the Varroa mites from Vietnam
Our findings support the suggestion of Fuchs et al. [26], who reported on an almost complete host specificity of the two Varroa lineages in Northern Vietnam, suggesting sympatry of two host specific Varroa types that do not hybridize. Despite only a minute divergence of the Vietnamese coxI haplotype from the Japanese and also the Korean V. destructor haplotype, our nuclear DNA analyses suggest a complete genetic isolation of the mites from the different host species. Not only did we not detect any indication of hybridization, we also failed to sample mites typical for the one host species in colonies of the other host species in Vietnam. Hence the Korean haplotype found in A. mellifera, which may have switched hosts only about 60 years ago [4], is not able to infect different populations of its original host species established in Northern Vietnam.
Our results show that the arm race between Varroa and its hosts has led to the evolution of very specialized mites. Although the underlying mechanisms of this coevolution are not well understood, previous work suggests that the mite is able to mimicry the cuticular hydrocarbons of its host [27][28] to avoid the hygienic behaviour of the honeybees [29]. Even though this trait appears to be plastic [30], the Korean haplotype is apparently not able to overcome the defenses of the populations of A. cerana found in Vietnam.
The level of genetic diversity and genetic structuring among the three sampling locations in Vietnam were higher in the mites sampled in A. cerana colonies than between the two locations where we sampled in A. mellifera colonies. This matches reports of the global spread of very few, genetically almost identical V. destructor lineages in A. mellifera [13].

Origins and diversity of the Varroa mites from the Philippines
The archipelago of the Philippines accommodates distinct and diverse A. cerana host populations that show haplotype variation at the subspecies level compared to mainland Asian populations [31][32]. This clearly sets the stage for independent coevolution between mites and hosts and may explain the large genetic differences between mainland Asia and the Philippine mites previously described by Anderson and Trueman [7]. In that study, Varroa from three provinces of the Philippines were analyzed: A. cerana colonies were sampled in two provinces on the island of Luzon: Batangas (which is the adjacent province located South of our sampling location) and San Fernando (a province located about 100 km to the North of our sampling location), and a third one in Mindanao, a different island. Each of these three sampling locations hosted a distinct mite haplotype in A. cerana colonies, grouping apart from the rest of the V. destructor sequences coming from A. cerana mites sampled in other Asian countries. However, the mites sampled from A. mellifera colonies in the Philippines all shared the Korean haplotype suggesting a separation of native and introduced mite populations as we observed in Vietnam in this study. We also found that the sequences of the mites we sampled in A. cerana in the Philippines differed significantly from the rest of the haplotypes previously described [14]. However, contrary to Anderson and Trueman [7], we also found the Luzon 1 haplotype in the Philippine A. mellifera colonies in Los Banos.
In addition to the lack of mitochondrial DNA variability in Los Banos, we also failed to detect any level of subpopulation structuring in this Varroa population. In fact, the microsatellite markers we analyzed suggested that the mites readily infect both host species. The presence of the native Philippine type in A. mellifera shows that more types than previously thought may be able to infect both Apis species.
A clearer picture on the Varroa genetic and functional diversity By coupling both mitochondrial and nuclear DNA approaches, we were able to infer the origin of the mites, but also to detect more functional mechanism such as the hybridization potential of different Varroa types in their native and non-native range. The substantial differences between the genetic diversity and infestation abilities of the mites sampled in Los Banos and in Vietnam may have far reaching consequences for our understanding of the host parasite biology of honeybees and Varroa.
The mites we sampled in Los Banos show all genetic prerequisites to qualify as a novel Varroa species. Although Anderson and Trueman [7] did not find morphological distinctiveness (based on body size) to completely separate this "Luzon haplotype 1" to the other Varroa species, we found further evidence that the mites from Los Banos differ markedly from the four other previously described Varroa species. Both mitochondrial DNA sequence divergence and the ability to parasitize both A. cerana and A. mellifera render this Varroa type a potential novel species. In addition, however, also the Varroa types we found in Vietnam appear to segregate as if they were distinct species. Although the mitochondrial haplotypes were rather similar to the Japanese haplotype of V. destructor, the microsatellite DNA markers showed a complete separation of these two mite groups even if kept on the same apiary. Thus, despite the fact that the number of substitutions is well below the 2% coxI divergence level typically considered as separating two species [33], there is no indication of any hybridization of nuclear markers. Given the Biological Species Concept of the reproductive isolation between the two sympatric groups [34][35], the mites with the Korean and the Vietnamese haplotypes could also be considered as two distinct species.
Solignac et al. [13] estimated the divergence time between the Korean and the Japanese Varroa type between 5 000 and 15000 years ago. Assuming a constant mutation rate and population size, we can estimate the divergence time among the various haplotypes found in our study. Four substitutions in the cox I sequence separate the Japanese from the Korean haplotype [14] resulting in an estimate 1250 and 3750 years of divergence per substitution. Similarly, the time of divergence between the two types found in Vietnam (which differ in by between five and eight substitutions) from the Japanese type would be between 10 000 and 30 000 years. Considering the short generation time of the Varroa mite [3], this is a rather long speciation period: This result could explain why the ubiquitous Korean haplotype can no longer parasitize A. cerana in Vietnam in spite of ongoing spill overs and spill backs in Korea and in Japan [14]. Following this reasoning, 50000 to 150000 years would separate the mites from the Philippines and the two Vietnamese haplotypes. This falls within the Pleistocene, during which A. cerana may have arrived in archipelago of the Philippines [31][32]. The ancestor of the Varroa mite found in the Philippines nowadays may have evolved in allopatry from the mainland populations since then.

Conclusions
We here provide an example of how host-parasite coevolution can rapidly lead to speciation within a short time span. The initial extreme selection on the Korean mites after the host switch has allowed for only a very limited number of individuals to reproduce in the novel host species less than a century ago. Subsequently, the constant brother-sister mating of Varroa has led to almost clonal population-specific mite types, which have differentiated considerably with time as they were taken away from their native region into allopatry. These extreme characteristics of the mite set the stage for the potential alloxenic speciation [36] of highly specific and most virulent types.
The consequences of keeping the introduced Western honeybee and the native Asian species is a most unfortunate example of transhumance having devastating consequences by promoting the global spread of parasites and associated viruses [8,[37][38][39]. Since apiculture has facilitated the global transmission of Varroa, selection will inevitably favor the most virulent types in A. mellifera as seen for the global spread of the Korean haplotype. Yet, given the tremendous increase of A. mellifera beekeeping in Asia and the wide diversity of Varroa in the native A. cerana populations, it seems possible that more mite types might switch to the western honeybee. At the same time, the spill back of virulent Varroa strains from A. mellifera to A. cerana may also become a risk potential for apiculture for these two economically and ecologically crucial species. Although some Varroa types are apparently strongly specific (Vietnam), others are more generalist (Philippines). If those generalist mites would spread to mainland Asia, it is likely that they would also invade both A. mellifera, but also A. cerana colonies.