Genetic Profile of the Parasitic Varroan Mite Varroa destructor (Arachnida: Mesostigmata: Varroidae) in Taiwan: a New Taiwanese Haplotype Intermediate Between the Highly Virulent Russian and Less Virulent Japanese Types Identified in the Honey Bee Host Apis cerana

© 2023 Academia Sinica, Taiwan Open Access Genetic Profile of the Parasitic Varroan Mite Varroa destructor (Arachnida: Mesostigmata: Varroidae) in Taiwan: a New Taiwanese Haplotype Intermediate Between the Highly Virulent Russian and Less Virulent Japanese Types Identified in the Honey Bee Host Apis cerana Tsen Hua1,2 , Panuwan Chantawannakul3,4, Cheng-Lung Tsai1,5 , and Wen-Bin Yeh1,* 1Department of Entomology, National Chung Hsing University, Taichung 402, Taiwan. *Correspondence: E-mail: [email protected] (Yeh). E-mail: [email protected] (Hua); [email protected] (Tsai) 2Department of Plant Medicine, National Pingtung University of Science and Technology, Pingtung 912, Taiwan 3Department of Biology, Faculty of Science, Chiang Mai University, Chiang Mai, 50200 Thailand. E-mail: [email protected] (Chantawannakul) 4Research Center of Microbial Diversity and Sustainable Utilization, Faculty of Science, Chiang Mai University, Chiang Mai 50200, Thailand 5Department of Entomology and Nematology, Ft. Lauderdale Research & Education Center, University of Florida, Ft. Lauderdale, FL 33314, USA Received 30 June 2022 / Accepted 3 January 2023 / Published 31 March 2023 Communicated by Benny K.K. Chan The Modern beekeeping industry is being challenged by the varroan mite and its transmitted pathogens. Various types of Varroa destructor exhibit different levels of virulence toward honey bees, but only the Japanese (J) and Russian (R) types were found to infect Apis mellifera. Type R was more highly virulent against A. mellifera in comparison with type J. Examining the genetic profile of Varroa species is therefore of crucial importance in apiary management. In this study, maternally inherited mitochondrial cytochrome oxidase I (COI) and bisexual nuclear internal transcribed spacer (ITS) sequences of V. destructor individuals from Taiwan were determined. All 168 COI sequences observed in populations obtained from A. mellifera were identical and belonged to type J, with one base difference to that of populations collected from A. cerana; the new type is named ‘T type’ (Taiwan type). ITS sequences of V. destructor and its sister species V. jacobsoni were identical. A network analysis based on 611 COI sequences compiled from references indicated the presence of 27 haplotypes in V. destructor. Epidemic history and relationship analyses of V. destructor showed that the basal haplotypes were those from A. cerana and many R-extending haplotypes infesting A. mellifera involving amino acid substitutions. Calibration dating based on COI analysis revealed that V. destructor differentiated from its sibling lineage (occurring in Sri Lanka) prior to 1.3 million years ago (Mya). The ancestral haplotype retention and drift in V. destructor that occurred locally during 0.10–0.64 Mya might be relevant to its host A. cerana, which had been isolated geologically. The highly virulent type R was spreading quickly and could gradually outcompete the common and less virulent type J. Type T, being intermediate between types R and J, ought to be studied to better understand the pathogenic mechanism of V. destructor in A. mellifera. Moreover, for areas where type R does not occur, such as Taiwan, quarantine requirements are crucial for reducing invasion risks. Key words: Mite, Varroa destructor, Bee, COI, ITS. Citation: Hua T, Chantawannakul P, Tsai CL, Yeh WB. 2023. Genetic profile of the parasitic varroan mite Varroa destructor (Arachnida: Mesostigmata: Varroidae) in Taiwan: a new Taiwanese haplotype intermediate between the highly virulent Russian and less virulent Japanese types identified in the honey bee host Apis cerana. Zool Stud 62:11. doi:10.6620/ZS.2023.62-11. Zoological Studies 62:11 (2023) doi:10.6620/ZS.2023.62-11 1 © 2023 Academia Sinica, Taiwan BACKGROUND Colony Collapse Disorder has considerably reduced the honey production and pollination efficiency of honey bees worldwide over the past three decades. Multiple factors can drive honey bee colony losses, the parasitic varroan mites and their transmitted viruses being among the most significant ones, exerting a substantial adverse effect on the modern beekeeping industry (Rosenkranz et al. 2010; Mondet et al. 2014; McMenamin and Genersch 2015). Four Varroa species have been described, namely Varroa destructor Anderson and Trueman, Varroa jacobsoni Oudemans, Varroa rindereri de Guzman and Delfinado-Baker, and Varroa underwoodi Delfinado-Baker and Aggarwal. Each varroan species has its Apis host (de Guzman and Delfinado-Baker 1996; de Guzman and Rinderer 1999; Wang et al. 2019); however, only V. destructor infects Apis mellifera (Anderson and Trueman 2000). For a long time, Varroa destructor was considered a biotype of V. jacobsoni until Anderson and Trueman (2000) described it as a distinct species. The varroan species that infects A. mellifera is larger than the one that infects A. cerana, and it was regarded as a distinct biotype with different levels of virulence toward these two bee species (Delfinado-Baker 1988; DelfinadoBaker and Houck 1989). Currently, we know that biotype A (also known as the Russian, German, or Korean type) is a group of V. destructor that infects A. mellifera and frequently causes its death. Biotype B is a group of V. jacobsoni that infects several Apis species but generally does not cause their death. Biotype C (also known as the Japanese type) is a group of V. destructor that infects A. mellifera but does not cause its death (Delfinado-Baker 1988; Delfinado-Baker and Houck 1989; Anderson and Trueman 2000). The body size and shape of V. jacobsoni and V. destructor are different; their size ratio is ≥ 1.4 (Dietemann et al. 2013). Cytochrome oxidase I (COI) sequences are different between V. jacobsoni and V. destructor (Anderson 2000; Anderson and Trueman 2000). Kraus and Hunt (1995) studied the genetic variability of V. destructor using randomly amplified polymorphic DNA (RAPD) analysis and reported that the Varroa species occurring in the United States has the same RAPD pattern as the one found in Germany. Based on a RAPD analysis, de Guzman et al. (1997) proposed that Eastern Russia and Japan might have been the original territories of the V. destructor infecting A. mellifera. Subsequent genetic studies have revealed that the RAPD pattern is consistent in identification with COI divergences (de Guzman et al. 1998; Anderson and Fuchs 1998; de Guzman and Rinderer 1999; de Guzman et al. 1999; Anderson and Trueman 2000). Anderson and Trueman (2000) recognized V. destructor and V. jacobsoni as two different species based on mitochondrial COI sequences and morphological characters, and revealed that V. destructor but not V. jacobsoni is harmful to A. mellifera; the divergence of their COI sequences is approximately 7%. In subsequent studies, COI sequences were extensively used for identification of these two species of Varroa. It has been established that V. jacobsoni distributed in Malaysia and Indonesia infects several Apis species but does not reproduce in A. mellifera, whereas V. destructor infects both A. cerana and A. mellifera. Before the introduction of the Western honey bee into East Asia, V. destructor was considered to infect only A. cerana, a bee species distributed in South and East Asia, including Pakistan, Nepal, India, Sri Lanka, northern Thailand, Vietnam, China, Korea, the Russian Far East, Japan, and Taiwan (Crane 1978; de Guzman et al. 1999; Anderson and Trueman 2000; Solignac et al. 2005). Based on the presence of sympatric populations of A. mellifera and A. cerana, two independent host switch events, one in Eastern Russia and the other in Japan in the 1950s, were proposed by Crane (1984) and de Guzman et al. (1997), respectively. The Russian type was probably first introduced to eastern Europe in the 1960s via transported brood combs. Subsequently, it spread to Germany and other parts of Europe, later establishing itself worldwide and becoming the main pest that damages A. mellifera (Crane 1978; de Guzman et al. 1997; Rosenkranz et al. 2010). Varroa destructor has been known to infest A. cerana in Japan since 1909, but had not been detected in A. mellifera in Japan until 1957 (Crane 1984). The Japanese type may have been exported to Taiwan and southeastern Asia during 1972 (Lo and Chao 1975). It was transported to South America in the 1970s and then to the United States in the late 1990s (de Jong and Goncalves 1981; de Jong et al. 1982; Delfinado-Baker and Houch 1989; de Guzman et al. 1997 1999). A comprehensive sampling by Anderson and Trueman (2000) enabled haplotype identification of V. destructor based on the mitochondrial COI amplicon. Each haplotype was assigned the name of the area (country or island) in which it was discovered in A. cerana. Of their six defined haplotypes (Nepal, Sri Lanka, Vietnam, China, Japan/Thailand, and Korea), only the Japanese and Korean haplotypes were found to infect A. mellifera. The Korean type was found to be identical with the Russian type (Anderson 2000), but in the present study, the Japanese (J) and Russian (R) types are used to represent haplotypes of the defined types of Japan/Thailand and Korea, respectively, in accordance with the original area of the relevant host-shift event (Crane 1978 1984) and the priorities of synonyms page 2 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan indicated by de Guzman et al. (1997) and Anderson (2000). A new, more complex haplotype naming system, based on newly available sequences of COI and other amplicons, e.g., “pure J”, “J1-2” or “J-like”, has been proposed (Solignac et al. 2005; Navajas et al. 2010; Traynor et al. 2020). Although the reasons why only types J and R are capable of infecting A. mellifera have remained unclear, type R has been found to be highly virulent toward A. mellifera and it can outcompete type J in localities where they are sympatric (de Guzman et al. 1997; Anderson 2000; Garrido et al. 2003; Strapazzon et al. 2009; Dietemann et al. 2019). By contrast, type J is considerably less aggressive and it does not cause mortality (Delfinado-Baker 1988; de Guzman et al. 1998; de Guzman and Rinderer 1999; Carneiro et al. 2007). Colonies of type R were found to be more abundant than type J in sympatric areas (de Guzman et al. 1999; Carneiro et al. 2007). A naturally Varroaresistant population of A. mellifera has also been observed worldwide (Locke 2016; Mondet et al. 2020). Microsatellite analysis indicated that hybridization events were not rare in areas where both types R and J occurred (Solignac et al. 2005; Dietemann et al. 2019) without the presence of any particular asymmetrical introgression, suggesting a postzygotic isolation between the two types. It is important to identify and elucidate haplotypes of Varroa mites since their identification is the basis for species management, which is necessary to protect the bee industry from an invasion of the highly virulent type R. Additionally, the identification of new haplotypes infesting A. cerana is highly significant, because such haplotypes pose a potential threat to A. mellifera, as it might become a new host in the future (Warrit et al. 2006; Navajas et al. 2010; Dietemann et al. 2019; Wang et al. 2019; Traynor et al. 2020). Thus, the haplotype profiles of V. destructor documented in Argentina, Brazil, China, Madagascar, Mexico, Spain, and Thailand in the last two decades are critical for inferring the evolutionary history of the species and developing management strategies against it (Medina and Martin 1999; Warrit et al. 2006; Muñoz et al. 2008; Navajas et al. 2010; Maggi et al. 2012; Rasolofoarivao et al. 2013; Octaviano-Salvadé et al. 2017; Dietemann et al. 2019; Wang et al. 2019; Traynor et al. 2020). Apis cerana started being cultured mainly in the southwestern areas of Taiwan more than 300 years ago, while A. mellifera was imported to the island in 1910 from the United States (Himura 1912; Inoue 1913; Hong 1962). Varroa jacobsoni was first recorded in 1975 in Taiwan, but the population was later identified as V. destructor (Lo and Chao 1975). The present study is based on a survey of V. destructor conducted throughout Taiwan between 2017 and 2020. A total of 14 Varroa colonies from A. mellifera and A. cerana in Taiwan, China, and Thailand were collected and their morphological and molecular characters were evaluated. The 211 sequences acquired in this study were analyzed in combination with 400 additional available COI sequences of V. destructor. Moreover, bisexual nuclear internal transcribed spacer (ITS) sequences of V. destructor were determined. Similar to hybrid detection using microsatellite markers, the identification of ITS sequences can be helpful for evaluating gene flow phenomena or hybridization events between various Varroa lineages (Solignac et al. 2005; Dietemann et al. 2019). Findings regarding haplotype profile and globally differentiated pattern can be informative for the development of V. destructor management strategies as well as quarantine requirement considerations, particularly in areas without severe damage by type R. MATERIALS AND METHODS Sample collection Adult females of V. destructor were collected from the bee brood cells of A. mellifera colonies in 12 apiaries across Taiwan including 2 populations of A. cerana (Fig. S1). During the autumns of 2017 to 2020, 182 Varroa mites (5–25 individuals from each apiary) were analyzed (Table S1). Moreover, 8 and 30 samples from A. mellifera apiaries in China (Jilin Province) and Thailand (Chiang Mai), respectively, were employed to determine the COI and ITS sequence variation among populations. Morphological description of Taiwanese V. destructor Measurements of body widths and lengths of 31, 3, and 23 individuals of Varroa mites collected from A. mellifera in Taiwan, China, and Thailand, respectively, as well as 5 individuals from A. cerana in Taiwan, were taken using a Microsight 4.1.2 (Q-Optics, USA) connected to a Canon EOS 800D camera (Tokyo, Japan) (Fig. 1). Specimen images were captured using a Nikon Coolpix B700 camera (Tokyo, Japan) with a Raynox DCR-250 macrolens (Tokyo, Japan). A statistical analysis of the measurements was carried out using the R software, following statistics proceedings described by R Development Core Team (R version 4.1.2). Analysis of variance (ANOVA) was used to test the measured variables for the aforementioned 4 varroan groups. The nonparametric Kruskal-Wallis test followed by Wilcoxon rank sum test for post hoc comparison page 3 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan (p < 0.05) were applied for the significant test of size measurement (Valdano and Di Rienzo 2007). DNA extraction, amplification, and sequencing Genomic DNA was extracted from adult female V. destructor individuals by using a Quick Extract DNA extraction kit (Epicentre Biotechnologies, Madison, WI). Primer sets were used to amplify the mitochondrial COI gene and nuclear ITS region. The sequences of forward and reverse primers for COI amplification, WX-990422 and ZSF-20060123, respectively, were reported by Wang (2007). Forward and reverse primers for the ITS region, namely ITSvarroa1 (TGARMCTGCGGARGGAWCATTAC) and ITSvarroa2 (CACACTTGATTTCAGATAAACATAAC), respectively, were designed in this study on the basis of the acquired Varroa mite sequences and conserved regions of ribosomal 18S and 28S DNA sequences (Kjer et al. 1994). Polymerase chain reaction (PCR) was performed in a volume of 25 µl containing 0.5 μl of each primer (10 mM), 0.2 μl of each dNTP (25 mM), 2.5 µl 10X Taq buffer, 0.5 µl Super Taq polymerase, and 1 μl of DNA template. The PCR programming conditions were as follows: 94°C for 2 min for denaturation, followed by 35 cycles at 94°C for 50 s, 42°C for 50 s, and 72°C for 50 s. The final extension was 72°C for 10 min. The PCR products were purified from 1% agarose gel by using the QIA Quick Gel Extraction Kit (Qiagen, Hilden, Germany) and then sequenced using an ABI 3730XL DNA Analyzer (Applied Biosystems, Foster City, CA, USA). Phylogenetic and sequence analyses The sequences of forward and reverse strands of each specimen were piled up with the published COI sequences of V. destructor and then aligned using ClustalW multiple alignment in Bioedit 7.0 (Hall 1999). In total, 211 COI sequences were acquired in this study, 168, 25, and 3 of which were collected from A. mellifera from Taiwan, Thailand, and China, respectively; and 15 were from Taiwanese A. cerana. Moreover, the 400 COI sequences of V. destructor obtained from GenBank and data in prominent references were compiled to elucidate the COI divergence (Table S2 and Fig. S2). The putative amino acid sequence was translated from the COI sequences and the pairwise variation of sequences was estimated based on the uncorrected proportional divergence by using MEGAX software (Kumar et al. 2018). The COI sequence divergences were also compared among V. destructor, V. jacobsoni (AF107902-AF107910), V. underwoodi (AF107260), and V. rindereri (AF107261). Phylogenetic analysis based on the COI haplotype sequences of V. destructor and V. jacobsoni was performed using Bayesian inference. The best-fit substitution model for COI was analyzed in jModelTest 0.1 by using the Bayesian information criterion (BIC) and was discovered to be GTR+I+G in MrBayes version 3.2 (Posada 2008; Ronquist et al. 2012). Three hot chains and one cold chain were applied for 1 × 108 generations with sampling every 1000 generations (Ronquist et al. 2012). The analysis was halted when the average standard deviation of the split frequencies was less than 0.002. Fig. 1. Morphological measurements of Taiwanese Varroa destructor. BL: body length, BW: body width. Scale bar in 200 μm was shown. page 4 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan The first 25% of trees were discarded as burn-in, and the remaining 75% of trees were used to construct a consensus tree. Haplotype separation was conducted using DNASP software (Rozas et al. 2003). Haplotype networks were analyzed using TCS version 1.21 with a statistical parsimony criterion (Clement et al. 2000). Based on the phylogenetic inferences, the most divergent type of V. destructor, namely the Sri Lankan type, was used for the outgroup comparisons in the network analysis. Moreover, DELTRAN optimization was applied in the haplotype network to determine the evolved processes of ancestral and derived haplotypes (Agnarsson and Miller 2008). This is because the most recent mutation of V. destructor in the new host, A. mellifera, is likely to occur as far as possible towards the tips in the network process. Calibration dating for COI sequences The divergence time of V. destructor haplotypes was estimated with the help of BEAST 2.6 (Bouchaert et al. 2014) by using the closely related haplotypes of Luzon 1 and Luzon 2, and those of V. jacobsoni as outgroups (Anderson and Trueman 2000). The bestfit model for nucleotide substitution was determined in jModelTest 0.1 by using the BIC (Posada 2008). GTR was the suitable model obtained in BEAST version 2.6 for COI sequences. Two substitution rates of 1.5% and 1.15% per lineage over a million years for COI sequences, commonly used for insects, were applied separately in this study (Brower 1994; Papadopoulou et al. 2010). A Markov chain Monte Carlo analysis was run for 1 × 108 generations with sampling performed every 1 × 103 generation. The effective sample size for the posterior distribution of estimated parameter values was analyzed using Tracer version 1.5 until the suggested value was reached (i.e., 200) (Rambaut and Drummond 2009). Subsequently, the initial 25% of the run was discarded as burn-in (Rambaut and Drummond 2009). RESULTS Morphological measurements of the V. destructor population in Taiwan Body lengths and widths of V. destructor specimens collected from A. mellifera in three different countries and from A. cerana in Taiwan are shown in table 1. Although a similar body length of 1200 m exists among the varroan mites, the body width is found to be significantly different in Taiwanese and Thai populations based on the nonparametric Kurskal-Wallis test. COI and ITS sequence variations of V. destructor populations in Taiwan compared with other varroan species The COI and ITS sequences of 458 bp and 516 bp, respectively, were successfully amplified. COI sequences of 196 and 15 specimens of V. destructor obtained from A. mellifera and A. cerana, respectively, were determined. COI and ITS sequences of V. destructor were deposited into the NCBI GenBank under accession numbers OK335527-OK35737 and OK336567-OK336693, respectively. The other two ITS sequences of Thai V. jacobsoni were deposited under accession number OK335749-OK335750. Individuals of V. destructor from 12 colonies across Taiwan were found to have identical COI and ITS amplicons except for a single base difference in the COI sequence in samples taken from A. cerana. The average genetic distances in the COI sequences between V. destructor samples from Taiwan in relation to V. underwoodi, V. rindereri, V. jacobsoni, the Mindanao varroan type, and the Luzon varroan type were found to be 9.0%, 8.6%, 6.6%, 4.8%, and 4.2%, respectively. In the ITS region, the average proportions in T, C, A, and G in the V. destructor specimens were 34.4%, 13.5%, 29.1%, 23.0%, respectively, and no sequence variation was observed between V. destructor and V. jacobsoni Table 1. Mean and stand deviation (SD) of body size measurement (mm) for Varroa destructor populations from Taiwan, Thailand, and China Measurement Body Length* Body Width* Ind. No1Host2 Population Mean ± SD Mean ± SD Taiwan 1.201 ± 0.051a1.720 ± 0.053b31 A. mel China 1.259 ± 0.126a1.783 ± 0.029b3A. mel Thailand 1.201 ± 0.061a1.695 ± 0.061c23 A. mel Taiwan 1.288 ± 0.071a1.836 ± 0.058a5 A. cer *Significant divergence was determined by Kruskal-Wallis test followed by Wilcoxon rank sum test for post hoc significance (p < 0.05). Different letters mean significantly divergent. 1Number of individuals measured. 2A. mel: Apis mellifera; A. cer: Apis cerana. page 5 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan (EF025475 and two Thai individuals of OK335749OK335750). Relationship of V. destructor haplotypes Twenty-seven haplotypes were identified among the 611 COI sequences obtained from V. destructor, among them six haplotypes previously identified by Anderson and Trueman (2000). All V. destructor COI from A. mellifera in Taiwan belonged to type J; one new haplotype, named here as the Taiwan type (type T), was found in Taiwanese A. cerana. A network analysis revealed a detailed evolved pattern among haplotypes (Fig. 2A). The Sri Lankan haplotype used as outgroup had at least 13 substitution steps to other haplotypes, and substitution steps ranged from 1 to 9 among the remaining 26 haplotypes, except the peculiar type MK482687 from Iran. The network originating from the sibling Sri Lankan haplotype is closely connected to the haplotypes of China3 and J. The newly identified type T from Taiwanese A. cerana is intermediate between types J and R. The network also shows that 13 out of 16 R-derived types found in tropical areas of Asia, i.e., southern India, Iran, and Saudi Arabia, involved amino acid substitutions (Fig. 3). Type R is found to be distributed worldwide except in Australia, Taiwan, and Puerto Rico (Table S2). The Vietnamese haplotype was also observed in China and Thailand. A Bayesian analysis conducted on COI haplotype sequences revealed that the lineages of Luzon 1 and Luzon 2 are distributed between V. jacobsoni and V. destructor (Fig. S3). In V. destructor, the Sri Lanka type is the most basal, followed by the types of J and Nepal. The remaining types either come from A. cerana (China, Vietnam, Taiwan, and Russia) or from A. mellifera, exhibiting a polytomous relationship (Fig. S3). Calibration dating of V. destructor Although a polytomous relationship is found among the shallow divergent haplotypes, calibration dating based on the COI haplotype sequences of V. destructor and V. jacobsoni revealed that the types from Nepal, Japan, Vietnam, and China found in A. cerana are grouped in the same cluster, whereas those found in A. mellifera in India, Iran, and Saudi Arabia are grouped with types R and T (Fig. 4). Calibration dating revealed that the lineage of V. destructor diverged from that Fig. 2. Haplotype network (A) and the possible evolved processes of ancestral and derived haplotypes according to the DELTRAN optimization (B) of Varroa destructor based on the COI haplotype sequences. The haplotype of Sri Lanka in the rectangle was used as the outgroup. Red and black circles represent the sampling mites parasitized in the host of Apis cerana and A. mellifera, respectively. The dotted red circles of Japanese and Russian haplotypes have both hosts of A. cerana and A. mellifera. A hypothetical haplotype is shown in a dot black circle. (A) The small circle corresponds to one individual, and the individual’s number of remaining haplotypes is marked within the corresponding medium and large circles. The nonsynonymous substitution position and its inducing amino acid changes between haplotypes are labeled. (B) The abbreviations of the haplotype are shown based on the haplotype name on panel A and Re means the extending haplotype from the Russian haplotype. C, C2, C2-1, C3, J, N, R, T, and V indicate the haplotypes of China, China 2, China 2-1, China 3, Japan, Nepal, Russia, Taiwan, and Vietnam, respectively. page 6 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan Fig. 3. Variable nucleotide (A) and amino acid sequence positions (B) of the COI of Varroa destructor haplotypes. Sequences identical to the first sequence of the Sri Lankan haplotype are indicated by a dot. Nonsynonymous nucleotide substitutions are shown in oblique and bold. Haplotypes from the parasite host of A. cerana are labeled in italic and those from A. mellifera are labeled in accession number with its acquired country. Fig. 4. Calibration dating of Varroa destructor based on COI haplotype sequences. Outgroups of Varroa jacobsoni and the varroan haplotypes of Luzon were used for comparisons. Divergent time is labeled in million years ago (Mya). Haplotypes from the parasite host of A. cerana are labeled in italic and those from A. mellifera are labeled in accession number with its acquired country. page 7 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan of V. jacobsoni approximately 2.6 million years ago (Mya) and that a subsequent differentiation occurred in the Luzon and Sri Lanka lineages approximately 1.3–1.8 Mya (Fig. 4). The differentiation of the remaining haplotypes found in A. cerana occurred at approximately 0.10–0.64 Mya. DISCUSSION Morphological characteristics of V. destructor populations in Taiwan and other areas Body size measurements and the target Apis host were found to be suitable for identifying four varroan species. Varroa underwoodi is the smallest in length at approximately 716 μm and is easily distinguished from the other three bigger varroan species with body lengths ranging from 1037 to 1180 μm (Delfinado-Baker and Houck 1989; de Guzman and Delfinado-Baker 1996; de Guzman and Rinderer 1999; Dietemann et al. 2013). While V. rindereri only parasites Apis koschevnikovi, its body length, approximately 1180 μm, is similar to both V. destructor and V. jacobsoni, which parasite A. mellifera and A. cerana. Moreover, a remarkably stable bimodality between body length measurements of V. jacobsoni and V. destructor was found, i.e., 1037–1063 and 1106–1167 μm, respectively (Delfinado-Baker and Houck 1989; Anderson and Truemen 2000). Furthermore, the body length of the Argentinan V. destructor is 1178 μm (Maggi et al. 2012). In this study, the body lengths of V. destructor collected from A. mellifera and A. cerana in Taiwan were 1201 and 1288 μm, respectively. These sizes are slightly larger than those reported previously. When varroan specimens had body lengths longer than 1106 mm, they could be recognized as V. destructor rather than V. jacobsoni (which had lengths of 1037–1063 μm), though body length was obviously variable in V. destructor populations. COI sequence differentiation in V. destructor In total, 27 V. destructor haplotypes are found worldwide, while only Jand R-related types could reproduce in A. mellifera and cause serious harm (Trynor et al. 2020). Type R, however, can infect more A. mellifera colonies and it can induce a more severe damage than type J (Delfinado-Baker 1988; de Guzman et al. 1999; Carneiro et al. 2007). Only type J has been identified in A. mellifera in Taiwan, and no colony collapse resulting from a type J infestation has been recorded, although apiary cleaning and pesticide spraying were necessary (Chen 2011). Moreover, a newly discovered haplotype, named type T, an intermediate between types J and R, was found in A. cerana in Taiwan. Currently, the R-related type has replaced type J in certain areas, such as Brazil, Japan, Thailand, and the United States, where mixed infestation has previously been found (Carneiro et al. 2007; de Guzman et al.1999; Dietemann et al. 2019). The genetic drift arising from host-shift events and the population expansion resulting from human activity may be responsible for the consistent genetic structure of types J and R of V. destructor worldwide. Although type R did not exhibit polymorphism in A. cerana, nearly all R tip-extending haplotypes found in A. mellifera in southern India, Iran, and Saudi Arabia contained amino acid substitutions, which may result from either particular environmental adaptations or other host-shift events. Moreover, new host-shift events caused by pre-adaptation or opportunistic events may have occurred because a few V. jacobsoni can survive in A. mellifera (Roberts et al. 2015). The network outgroup comparisons, host-shift events, and artificial transportation of the current 27 haplotypes provide evidence for identifying ancestral haplotypes in V. destructor. The major haplotype typically recognized as an ancestral haplotype in population genetic studies may not be found in V. destructor because the current major J and R types had occurred through recent host-shift events in modern apiculture. The intermediate haplotypes closer to the outgroup, however, could be recognized as ancestral states that evolved prior to their tip-extended haplotypes (Fig. 2B). The local haplotypes derived from the R type found in A. mellifera, however, may have mutated recently. COI sequences have been commonly applied in molecular calibration related to speciation events involving temperate Palaearctic and tropical Amazonian biota (Avise 2000; Knowles 2000; Hoorn et al. 2010; Rull 2011). Figure 4 suggests that V. destructor differentiated from the Sri Lankan type prior to 1.3 Mya. Ancestral haplotype retention and genetic drift subsequently occurred in local populations, as observed in the intermediate haplotypes during 0.10–0.64 Mya. This estimation was consistent with that for the period 0.15 Mya on the basis of microsatellite markers by Solignac et al. (2015). Rapidly arising mutations of V. destructor have occurred in A. mellifera from India, Iran, and Saudi Arabia, which exhibit a likely pattern of varroan mites in Papua New Guinea (Roberts et al. 2015). The differentiation history of V. destructor on its natural parasite host A. cerana was explored. The population structure inferred from morphometric and DNA sequence analyses indicated that A. cerana page 8 of 12Zoological Studies 62:11 (2023) © 2023 Academia Sinica, Taiwan comprises two major clusters and several submajor clusters relevant to refugia formation, with these clusters affected by Pleistocene glacial cycles (Smith and Hagen 1996; Radloff et al. 2010). Ancestral haplotype differentiation and retention in local V. destructor populations might be relevant to the refugium of its host during Pleistocene glacial cycles. Hybridization in Varroa mites Two independent host-shift events occurred in both the Russian Far East and Japan. The possible occurrence of R–J hybrids has been proposed because they have been observed to be sympatric in Brazil, Japan, Thailand, and the United States (de Guzman et al. 1999; Carneiro et al. 2007; Dietemann et al. 2019). A worldwide sampling analysis based on microsatellite markers revealed the occurrence of hybrid events and did not point to the presence of asymmetrical introgression between the J and R types, although many hybrids involving J alleles in the genetic background of the R type have been reported (Solignac et al. 2005; Dietemann et al. 2019). A postzygotic barrier has been proposed to serve as a mechanism for selection against type J, possibly explaining why type R replaced type J and coexisted with it in sympatric areas (Solignac et al. 2005). Examining the genetic content of R–J hybrids in sympatric areas would provide additional information on the selection fitness of the highly virulent type R relative to the less virulent type J. Gene expression levels are key factors that can help elucidate why type R but not type J is highly virulent for A. mellifera. Allozyme patterns used to examine genetic differentiation furnished no recognizable differences for populations of V. destructor in Brazil, China and Europe (Issa 1989; Biasiolo 1992). Zhang et al. (2010) reported that V. destructor exhibited differences in expression of genes related to metabolic functions and nerve transmission signature of A. mellifera versus A. cerana. Oldroyd (1999) indicated that bees have been exposed to type R for the longest period in eastern Russia; thus, tolerance to Varroa is most common in that region, although many naturally resistant populations of A. mellifera have been found (Locke 2016). Moreover, on the basis of genome-wide analyses, Techer et al. (2019) found that divergent selection regimes existed between V. destructor and V. jacobsoni. In this study, network analysis indicated that type T was linked to the most common types R and J. The mechanism underlying the tolerance or virulence difference between the R and J types in A. mellifera could be elucidated using the available genomic sequences of the J, R, and T types. Hybridization events between V. jacobsoni and V. destructor might be also possible. Spillbacks of V. destructor to A. cerana and spillovers of V. jacobsoni to A. mellifera have been found in sympatric populations in Thailand, where reciprocal genetic admixture has occurred between V. destructor and V. jacobsoni (Dietemann et al. 2019). The ITS sequence (EF025470) of V. jacobsoni was obtained from Java, Indonesia, where coinfestation of the two mites has commonly occurred (Anderson and Morgan 2007; Anderson and Fuchs 1998). Therefore, the identical ITS sequence in V. jacobsoni and V. destructor may have resulted from hybridization events. However, the fact that V. destructor and V. jacobsoni have identical ITS patterns may also have resulted from the ancestral retention of an ITS region. Quarantine management for colonies of A. mellifera Identifying the isolates of V. destructor species, especially those of the highly virulent type R, is crucial in the prevention of invasions of Varroa mites. Solignac et al. (2005) and Navajas et al. (2010) examined microsatellite markers and COI sequences in 12 individuals collected from Taiwanese A. mellifera and found pure J and J-like types of the Taiwanese V. destructor. In this study, a sufficient sample obtained across Taiwan revealed the presence of only the less virulent J type in A. mellifera. R-free regions have been recorded in some areas. Honey bee movement in each country should be examined, and quarantine requirements should be implemented for the trade of live honey bees worldwide to prevent the introduction of the R type to R-free countries, islands, or areas (Iwasaki et al. 2015; Traynor et al. 2020; Hall et al. 2021). The haplotype composition and distribution areas shown in the table S2 and figure S2 can support the clarification of existing haplotypes and examination of new haplotypes in V. destructor. For varroan mite management, Oldroyd (1999) also suggested that the virulent R type could be controlled by deliberately releasing the less virulent J type. CONCLUSIONS Two hundred and eleven COI were sequenced and compared with the 400 available COI sequences of V. destructor. A total of 27 haplotypes were discovered worldwide, and only the J and R-derived types were found to be able to reproduce in A. mellifera. In Taiwanese apiaries, the less virulent type J was found, highlighting Taiwan’s success in preventing an invasion of the highly virulent type R. Honey bee page 9 of 12Zoological Studies 62:11 (2023)