First Evidence of Presence of Varroa underwoodi Mites on Native Apis cerana Colonies in Primorsky Territory of Russia Based on COX1 Gene

In modern beekeeping, honey bee colonies are moved over long distances to improve crop pollination and increase the honey harvest. Honey bees are introduced to ecosystems beyond their natural distribution range and frequently exposed to new pathogens and parasites that they have never encountered before, including ectoparasitic mites of the genus Varroa (Varroidae) represented by at least four species (Anderson & Trueman, 2000; Rosenkranz et al., 2010). The first species, Varroa jacobsoni Oudemans 1904, originally parasitized on A. cerana on Java island of Indonesia (Oudemans, 1904) for the first time spilled over to Apis nigrocincta in Indonesia (Hadisoesilo & Otis, 1998; Anderson & Trueman, 2000) and to Apis mellifera in Papua New Guinea (Roberts et al., 2015). The second species, Varroa destructor Anderson & Trueman 2000, originally parasitized on Apis cerana in China, Japan, Korea, and Thailand for the first time spilled over to Apis mellifera in Japan (Beaurepaire et al., 2015). The third species, Varroa underwoodi Delfinado-Baker & Aggarwal, 1987, originally parasitized on Apis cerana in Nepal, for the first time spilled over to A. nigrocincta in Indonesia (Delfinado-Baker & Aggarwal, 1987; Anderson et al., 1997; Kuznetsov, 2005). The fourth species, Varroa rindereri de Guzman & Delfinado-Baker, 1996, originally parasitized on Apis koschevnikovi on the Borneo island of Malaysia (de Guzman & Delfinado-Baker, 1996) does not have any spillover events to other host species yet.

The honey bee A. mellifera is predominantly parasitized by V. destructor worldwide (Traynor et al., 2020), and by V. jacobsoni in Papua New Guinea (Roberts et al., 2015). Several mitochondrial V. destructor haplotypes have been described, but only two of them - Korean (K) and Japanese (J) are capable of reproducing in A. mellifera colonies (Anderson, 2000; Anderson & Trueman, 2000; Muñoz et al., 2008). Up to now, two Varroa species - V. rindereri and V. underwoodi have been poorly studied (Oudemans, 1904; Delfinado-Baker & Aggarwal, 1987; de Guzman & Delfinado-Baker, 1996; Anderson et al., 1997; Rath, 1999; Anderson & Trueman, 2000; Wang et al., 2019).

The distribution area of V. underwoodi is constantly expanding and spans populations of A. cerana in Nepal (Delfinado-Baker & Aggarwal, 1987), South Korea (Woo, 1992; Chantawannakul et al., 2016), Indonesia (Anderson et al., 1997; Chantawannakul et al., 2016), Papua New Guinea (Lee, 1995; Anderson et al., 1997; Chantawannakul et al., 2016), Vietnam (de Guzman & Rinderer, 1999; Chantawannakul et al., 2016) and China (de Guzman Rinderer, & 1999; Huang, 2004; Chantawannakul et al., 2016; Wang et al., 2019). Russia was not listed as a proven distribution area for V. underwoodi according to English language publications (Chantawannakul et al., 2016; Wang et al., 2019), even though Russian language publications had been previously described morphometrically from feral A. cerana population in the Primorsky Territory (Kuznetsov, 2005; Kuznetsov & Lelej, 2005).

Along with the expansion of the area of V. underwoodi, the number of its host species has increased: A. cerana in Nepal (Delfinado-Baker & Aggarwal, 1987), Apis nuluensis in Malaysia (Delfinado-Baker & Aggarwal, 1987; de Guzman et al., 1996; Anderson et al., 1997), A. nigrocincta in Indonesia (Anderson et al., 1997; Hadisoesilo, 1997), A. mellifera in Papua New Guinea (Lee, 1995; Anderson et al., 1997; de Guzman & Rinderer, 1999). Although evidence for V. underwoodi reproduction was only found in A. cerana, the frequent reports on this mite species appearing in the colonies of other species suggested abundant opportunities for cross-species transmission. V. underwoodi is especially dangerous to A. mellifera colonies because they are kept near A. cerana colonies in most Asian countries (Zheng et al., 2011, 2018; Chantawannakul et al., 2016; Wang et al., 2019; Roberts et al., 2020).

Varroa species haplotypes have different virulence for host species and only K and J of the V. destructor six haplotypes are capable of reproducing on A. mellifera (Anderson, 2000; Anderson & Trueman, 2000; Muñoz et al., 2008). Similarly, a high level of genetic diversity of V. underwoodi (Navajas et al., 2010; Roberts et al., 2015; Wang et al., 2019) allows several haplotypes to form with various virulences for host species, including honey bee A. mellifera. Besides, the manifestation of similar traits in related species is supported by the law of homologous series (Vavilov, 1920).

Thus, honey bee A. mellifera is a potential host species for V. underwoodi in its further evolution (Anderson, 2000; Anderson & Trueman, 2000; Muñoz et al., 2008). The potential factors for V. underwoodi to become a parasite of honey bee A. mellifera must be investigated to mitigate their further negative effects and prevent future invasions (Thompson, 1994; de Guzman & Rinderer, 1999; Kolar & Lodge, 2001; Woolhouse et al., 2005). In this study, we used morphometry and mitochondrial COX1 gene sequences to prove the presence of V. underwoodi on the northernmost feral A. cerana subspecies Apis cerana ussuriensis inhabiting Russia's Primorsky Territory.

The adult Varroa underwoodi mites were collected in 70% ethanol during summer 2004 from brood cells of two managed Apis cerana ussuriensis colonies in the village Romashka of the Khasansky district, the Primorsky Territory (43.5N, 131.3E). The reproduction ability of V. underwoodi in A. cerana colonies was assumed when 5–6 adult females and 2–3 nymphs were found together in drone brood cells. All collected mites were subsequently stored in 70% ethanol at −20°C until needed for further analyses. V. underwoodi mites were exclusively found in A. cerana colonies, where the colony infestation rate was 50%. The drone brood infestation rate was 2.8% in June, 35% in July, 58% in August. The worker brood infestation rate was 1%.

For a preliminary confirmation of species identity, the morphometrics and size of adult V. underwoodi female mites (N=10) were compared with previous reports (Delfinado-Baker & Aggarwal, 1987; Woo, 1992; Anderson et al., 1997; Huang, 2004; Wang et al., 2019). The sampled V. underwoodi mites were dried with ethanol at room temperature for one minute. The size of the dorsal shield with lateral setae was used to characterize V. underwoodi (Delfinado-Baker & Aggarwal, 1987; Woo, 1992; Anderson et al., 1997; Huang, 2004; Wang et al., 2019). The morphometry of each individual was measured with an EOS Kiss X7 digital microscope (Canon, Japan) with the lens MP-E 65mm f/2.8 1–5× Macro Photo (Canon, Japan) under 150× magnification according to the manufacturer's instruction.

The total DNA of A. c. ussuriensis was extracted from three mites per colony according to Qiagen DNEasy protocol for animal tissue (Qiagen, Valencia, Ca.). A mitochondrial COX1 gene sequence was used to identify mite species as V. underwoodi and to identify its particular mtDNA haplotype. PCR amplified the V. underwoodi COX1 gene according to Wang et al. (2019) using a pair of primers (COX1_821_F: 5′-GGAGTAGGTACAGGTTGAACGG-3′ and COX1_821_R: 5′-ACAACCCCAGCAATAATAG-CAA-3′) with 821 bp product (Wang et al., 2019). The PCR-amplified fragments were sequenced with Sanger's methods for all V. underwoodi samples, with the use of a pair of primers (F-V51: 5′-GTAATTTGTATACAAAGAGGG-3′ and R-V1400: 5′-CAATATCAATAGAAGAATTAGC-3′) (Warrit et al., 2004).

All PCR products were purified with the QIAquick PCR Purification Kit (250) (QIAGEN, Hilden, Germany) according to the instructions of the manufacturer. The nucleotide sequences of the COX1 gene of V. underwoodi samples were determined through the sequencing of the PCR products using the Sanger dideoxy method (Sanger et al., 1977) on the ABI 3730×l (Applied Biosystems, Foster City, CA, USA) with the ABI PRISM BigDye Terminator v3.1 Cycle Sequencing Kit according to the manufacturer's instructions. All PCR products were sequenced from both strands. The 458 bp sequence of the mitochondrial COX1 gene of V. underwoodi was uploaded into the DDBJ/GenBank database with accession number LC532104.

The V. underwoodi COX1 gene sequences - MH205173 (Hangzhou, China), MH205174 (Jinhua, China), MH205175 (Nanchang, China), MH205176 (Jilin, China), MH205177 (Maoming, China)), V. destructor (KJ403739, KJ507740, KJ403742, KJ403744 (Riyadh, Saudi Arabia)), and V. jacobsoni - (MF462134 (Moresby, Papua New Guinea), AF010479 (Canberra, Australia) from GenBank were compared with V. underwoodi from Russia's Primorsky Territory. The samples, which were closely related to species V. destructor and V. jacobsoni were used for outgroup comparison.

The genetic divergences among V. underwoodi, V. destructor, and V. jacobsoni COX1 gene sequences were estimated using p-distance (proportion of nucleotide sites at which two sequences differ) with CLUSTALW alignment as implemented in MEGA 10.0.5 (Kumar et al., 2018). A neighbor-joining (NJ) phylogenetic tree based on the p-distances of the COX1 gene sequences with 2000 bootstrap replications was constructed with CLC Genomics Workbench 20 (Qiagen Inc., Mississauga, ON, Canada). The statistical analysis and the analysis of molecular variance (AMOVA) were performed with the use of ARLEQUIN 3.5.2 (Excoffier & Lischer, 2010), STATISTICA 8.0 (StatSoft, OK, USA), and EXCEL 2010 (Microsoft, CA, USA).

The taxonomic affiliation of Varroa mite samples was provided using morphometry measurements and mitochondrial COX1 gene sequence analysis. The colour of the ellipsoidal body of V. underwoodi females was chestnut brown (Fig. 1). The dorsal shield surface was lightly striated and reticulated with tightly covered setae, which are approximately the same length and a little spiny. The setae on each lateral edge gradually increased in length posteriorly, with the last three pairs decreasing in size again. The body length of adult female V. underwoodi was 767.50±20.5 μm (mean±SD), and the width was 1,300.50±20.5 μm (n=10). For comparison, the body lengths and widths of adult female V. underwoodi were the following: A. cerana - 700–752 μm and 1,089–1,157 μm (n=15); A. mellifera - 700–735 μm × 1,090–1,120 μm (n=6); A. cerana in Irian Jaya - 690–730 μm × 1,050–1,130 μm (n=5); A. cerana in Sulawesi and Java - 720–780 μm × 1,050–1,080 μm (n=2); A. nigrocincta in Sulawesi - 740–760 μm × 1,120–1,220 μm (n=5) (Anderson et al., 1997); A. cerana in Nepal - 741–780 μm × 1,151–1,168 μm (n=2) (Delfinado-Baker & Aggarwal, 1987); A. cerana in South Korea - 703–784 μm × 1,135–1,324 μm (n=2) (Woo, 1992). Due to the congruence of morphology parameters to the previously published morphology of V. underwoodi (Delfinado-Baker & Aggarwal 1987; Anderson et al., 1997; Huang 2004; Wang et al., 2019), the Varroa mite samples were assumed to belong to the species V. underwoodi.

Fig. 1

The COX1 gene sequences of six V. underwoodi samples from both A. c. ussuriensis colonies #2 and #5 were identical to one another and to MH205176 (Jilin, China) (Wang et al., 2019). We called both of them a haplotype China 1 MH205176. The pairwise differences among V. underwoodi, V. destructor, and V. jacobsoni were counted based on the COX1 gene sequences polymorphism. The pairwise number of nucleotide and amino acid differences, p-distances, and percent of genetic divergence based on mitochondrial COX1 gene sequences is presented in Tab. 1.

COX1 gene sequences of mite species differ significantly (p≤0.05). It confirms with a 95% probability that COX1 sequences of V. underwoodi, V. destructor, and V. jacobsoni are quite distinct species. Of the three species, V. destructor, and V. jacobsoni are closest with a 7% value of genetic divergence. The species V. underwoodi differs from both mite species equally with a 10% value of genetic divergence.

The pairwise genetic divergences and sequence differences of nucleotides and amino acids among each sample of V. underwoodi (N=6), V. destructor (N=6), and V. jacobsoni were counted based on a comparison of the mitochondrial COX1 gene sequences. The pairwise number of nucleotide and amino acid differences, p-distances, and percent of genetic divergence are presented in the table (Tab. 2).

In this study, there was a variation in the genetic divergence (0% to 2%), p-distance (0.000 to 0.022), number of nucleotide differences (0 to 10), and number of amino acid differences (0 to 8) in the V. underwoodi samples. Most of them differed from all other V. underwoodi samples MH205177 from China, Maoming. There were no differences between MH205176 from Jilin, China, and LC532104 from Primorsky Territory, Russia. In the V. destructor samples, the genetic divergence varied from 0% to 1%, p-distance varied from 0.000 to 0.022, the number of nucleotide differences from 0 to 4, and the number of amino acid differences from 0 to 1. The samples KJ403744 and KJ403739 from Saudi Arabia, Riyadh are the least different. No differences were found between KJ403742 and KJ403740 from Riyadh, Saudi Arabia. In the V. jacobsoni samples, no differences were found between MF462134 from Moresby, Papua New Guinea, and AF010479 from Canberra, Australia. A NJ phylogenetic tree was constructed based on pairwise p-distances of the mitochondrial COX1 gene sequences among three species of mites V. underwoodi, V. destructor, V. jacobsoni. The COX1 gene sequence LC532104 from Primorsky Territory, Russia (haplotype China 1 MH205176) was clustered with all V. underwoodi sequences MH205173, MH205174, MH205176, MH205177 (China), which was distinct from V. destructor sequences KJ403739, KJ507740, KJ403742, KJ403744 (Riyadh, Saudi Arabia) and V. jacobsoni sequences MF462134 (Moresby, Papua New Guinea), AF010479 (Canberra, Australia) (Fig. 2). On the phylogenetic tree, the sample of V. underwoodi LC532104 (Primorsky Territory, Russia) was located the closest to the northern-most sample MH205176 (Jilin, China) (distance of 450 kilometers) and the farthest from the southernmost sample MH205177 (Maoming, China) (distance of 3000 kilometers). On the tree, V. underwoodi was genetically closer to V. destructor than to V. jacobsoni. Compared to all V. underwoodi samples, only the southern-most sample MH205177 (Maoming, China) was closest to the V. destructor samples.

Fig. 2

The V. underwoodi mites in the Primorsky Territory have probably always parasitized in wild A. cerana colonies living in tree hollows. Conditions are more favourable for the reproduction of V. underwoodi in hives than in tree hollows. Kuznetsov (2005) stated that in 2002 the V. underwoodi mites were found only once in the A. cerana drone brood in the apiary, but in 2004, a massive reproduction of V. underwoodi was observed. In 2004, the infestation rate of A. cerana colonies with the V. underwoodi mites was 2.8% in June, 35% in July and 58% in August in the drone brood cells, and was 1% in worker brood cells in the apiary. Individual A. cerana drone cells contained up to five to six adults and two to three V. underwoodi nymphs. V. underwoodi imagos and nymphs were found on the same pupae of A. cerana drones and were rarely observed in the brood of worker bees. Female V. underwoodi were found mostly on young drones and very rarely on young worker bees of A. cerana. In summer 2004, up to 16% of young drones were infested with V. underwoodi mites, but in autumn, no mites were found on adult bees. In 2004, the high rate of infestation by V. underwoodi mites in the brood contributed to the frequent cleansing swarming of A. cerana colonies (Kuznetsov, 2005).

The Varroa mites sampled in the Primorsky Territory, Russia were assumed as subspecies V. underwoodi based on their morphometrics that was in the range of previously described populations (Delfinado-Baker & Aggarwal 1987; Anderson et al., 1997; Huang 2004; Wang et al., 2019) and mitochondrial gene COX1 matched with the sample available from the province Jilin of China MH205176 (Wang et al., 2019). The body length and width of adult female V. underwoodi from Russia were a little bigger in comparison with other populations, which can be explained by the northernmost distribution. Since the COX1 gene sequence of V. underwoodi from Primorsky Territory, Russia LC532104 was identical to that of Jilin, China MH205176, we defined both of them as a haplotype China 1 MH205176 (Traynor et al., 2020).

The level of genetic divergence based on mitochondrial DNA sequences between insect species within the genera varied from 8% to 17%, and the genetic p-distance based on mitochondrial DNA sequences varied from 0.100 to 0.200 (Tan et al., 2011; Han et al., 2016; Eimanifar et al., 2017; Ilyasov et al., 2018, 2019). In our current study, the genetic divergence of mitochondrial COX1 gene sequence ranges from 7% to 10%, and genetic p-distances ranges from 0.072 to 0.099 between three species of mites V. underwoodi, V. destructor, and V. jacobsoni, which is almost matched with the range of the mtDNA-based intraspecific level of genetic differences in insects (8–17%) (Tan et al., 2011; Han et al., 2016; Eimanifar et al., 2017; Ilyasov et al., 2018, 2019). The levels of divergence described here between V. destructor and V. jacobsoni 7% is a bit higher than 5.9%, which was described, in a recently published paper (Techer et al., 2019). This difference can be explained by averaging a different number of samples of mites used for comparisons; we used only four samples of V. destructor from Saudi Arabia and only two samples of V. jacobsoni from Papua New Guinea and Australia, whereas in the recently published paper we used seven samples of V. destructor from South Korea, France, Vietnam, China, Japan, Nepal, Sri Lanka and eleven samples of V. jacobsoni from Indonesia, Malaysia, Laos, Borneo, Papua New Guinea (Techer et al., 2019). Previous studies based on the morphology and mitochondrial COX1 gene had reported the occurrence of V. underwoodi in A. cerana colonies in Nepal (Delfinado-Baker & Aggarwal 1987), South Korea (Woo, 1992; Kuznetsov, 2005; Chantawannakul et al., 2016), Indonesia (Anderson et al., 1997; Chantawannakul et al., 2016), Papua New Guinea (Lee, 1995; Anderson et al., 1997; Chantawannakul et al., 2016), Vietnam (de Guzman & Rinderer, 1999; Chantawannakul et al., 2016) and China (Huang 2004; Wang et al., 2019). In the previous studies, Russia had not been included in the distribution area of V. underwoodi (Chantawannakul et al., 2016; Wang et al., 2019). To date, this is the first evidence of V. underwoodi distribution in far east Russia in native A. c. ussuriensis colonies using both morphometry and mitochondrial COX1 gene sequencing methods. The distribution of V. underwoodi in Khasansky district, Primorsky Territory, Russia (43.5N, 131.3E) is at a distance of 450 kilometers from the distribution of V. underwoodi in Jilin province, North China (43.1N, 127.1E).

The presence of V. underwoodi in capped worker brood cells in one A. mellifera colony in Papua New Guinea (Roberts et al., 2015) shows that interspecies host switch can occur. To date, V. underwoodi has not been found yet in A. mellifera colonies in far east Russia. This suggests that interspecies transmission to A. mellifera is rare, which considerably limits chances for host switch. The population structure of V. underwoodi mites studied previously suggested that genetic diversity might lead to evolving lineages to reproduce in A. mellifera colonies. As has been the case for V. destructor and V. jacobsoni, rare events can lead to interspecies host switches with devastating effects (Rosenkranz et al., 2010; Roberts et al., 2015). Because the V. underwoodi samples are genetically closest to V. destructor samples that are a common parasite of both A. mellifera and A. cerana, V. underwoodi also has a high probability to host switch and parasitizing on A. mellifera (Roberts et al., 2020).

The infestation rate by V. underwoodi was significantly higher in China's northern provinces than in its southern ones (Wang et al., 2019). The higher infestation rate in the cold region is unsuspected because mass parasite reproduction is not possible when colonies desist brood rearing during winter. The higher infestation rate in northern Asia is affected by drone brood production time that is shorter in a cold climate (Wang et al., 2019). The identity of the COX1 gene sequences of V. underwoodi from Primorsky Territory, Russia LC532104 and from Jilin, China MH205176 located 450 kilometers away can presumably be explained by the fact that both mite populations are predominantly adapted to parasitizing on A. cerana in the cold climate of northern Asia. The identity of the COX1 sequence of the V. underwoodi samples from far east Russia and northern China suggests that this region is inhabited by one single population of V. underwoodi, parasitizing a common host, A. cerana. This indicates that migration between Chinese and Russian A. cerana populations is associated with the exchange of genes and parasites. Thus, haplotype China 1 MH205176 of V. underwoodi can be the first candidate for a cross-species host switch from north Asian A. c. ussuriensis and north European A. m. mellifera, both of which inhabit a cold continental climate. The COX1 gene sequence LC532104 from Primorsky Territory, Russia clustered with all V. underwoodi sequences from China and formed a common subgroup with the northern-most sample MH205176 from Jilin, China, and belonged to one haplotype China 1 MH205176. This suggests that the samples of V. underwoodi from the far east of Russia and North China are representatives of one common big North Asian population containing haplotype China 1 MH205176. This parasite V. underwoodi seems to originate from south Asian areas and later spreads to north Asia together with A. cerana migration. The low level of genetic differences of V. underwoodi samples collected between far east Russia and North China after isolation by distance and geography can be explained by a low rate of molecular evolution due to a low level of natural selection through life in a stable environment inside A. cerana colonies. Similarly, no differences were found between V. destructor samples KJ403742 and KJ403740 from distant regions of Saudi Arabia, between V. jacobsoni samples MF462134 and AF010479 from distant countries Papua New Guinea and Australia. Thus, the rate of molecular evolution in all Varroa mite species is very low, and their genomes are more conservative than the genomes of their host species of genus Apis.

Until recently, the distribution area of V. underwoodi covered almost all countries where cerana is found, with the exception of Russia. On the basis of COX1 gene sequences and morphometry analysis, we expanded the range of V. underwoodi by 450 km to the east to the Primorsky Territory, Russia. The identity of the COX1 gene sequences of the V. underwoodi from northern China and far east Russia indicates the absence of state borders between countries for A. cerana, which freely migrate and spread parasites. The northernmost boundaries of the V. underwoodi range remain unexplored. Last year, we collected northernmost A. cerana samples near the Terney village of the Primorsky Territory (45.06N 136.61E) – the northern border of the species range. We assumed that the northern border of the V. underwoodi range coincides with A. cerana. Furthermore, we are going to characterize V. underwoodi using additional COX3-ATP6 and CYTB markers. It is possible we can find differences between Russian and Chinese V. underwoodi samples and can demonstrate that its biogeography is more complex than previously assumed.