Genomic analyses of sibling honey bee ectoparasitic mite species show divergent strategies of adaptation

30 Multispecies host-parasite evolution is common, but how parasites evolve after speciating remains poorly understood. On one hand, their shared evolutionary history and physiology may propel them along similar evolutionary trajectories. Alternatively, they may pursue different strategies to reduce competition with each other. Here, we test these scenarios using the economically important association 35 between honey bees and ectoparasitic mites by sequencing the genomes of the sister species Varroa destructor and Varroa jacobsoni. We also compare them to another honey bee mite (Tropilaelaps mercedesae). We find different sets of genes and gene ontology terms under selection in each of the lineages, indicating distinct selective regimes operating on each of the parasites. Divergent strategies pursued 40 by the parasites may make it harder for the host species to develop tolerance to all of them at the same time. Based on our findings, we suggest that species-specific strategies may be needed to combat evolving parasite communities. Interactions between hosts and parasites shape biodiversity by driving 45 coevolutionary arms races and by regulating populations over ecological timescales. Parasitism is a remarkably successful strategy, occurring in nearly half of recognized animal phyla 1,2. Parasites depend on their hosts, which form their principal ecological niches and evolve rapidly, often with multiple parasite generations for each host generation. This results in convergent selective regimes, 50 leading to reproducible evolution of genomic features such as reduced genomic compaction, functional losses often associated with a reduced metabolism, and adaptive functional gains 3. Most models of parasite coevolution are limited to interactions between single focal parasites and their hosts, although multi-species dynamics are known to be much more complex 4,5. Furthermore, in nature, 55 communities of related species may parasitize individual hosts 6–8. As a result, coevolutionary dynamics involving more than one parasitic species remain poorly understood 9,10. For example, co-infecting parasites occupying the same host (ecological niche) can exclude competitors via direct, or interference competition 11,12. Spatial 60 separation within regions of the host gut is a strategy to reduce interspecies competition adopted by cestodes naturally co-infecting sticklebacks 13. Within the same host, congeneric parasites can co-exist, but spatial segregation or slight differences in niche parameters may explain intrahost speciation, as found in gill parasites (Dactylogyrus) of freshwater fish (Cyprinidae) 14. Niche partitioning across 65 seasons allows two fungal sibling species (Erysiphe quercicola and E. alphitoides) to exploit one host tree (Quercus robur) by reducing direct interferences 15. Alternatively, parasites may shift their strategies to exploit different aspects of the host’s niche. For example, fire ants (Solenopsis spp.) are parasitized by a genus of decapitating phorid flies (Pseudacteon), which are highly host-specific and have the 70 same life cycle, but specialize on attacking ants in different contexts (e.g., while foraging, vs. at the nest.) 16. However, almost nothing is known about the evolutionary trajectories of speciating parasites. On one hand, they share physiology and host-specific adaptations, which may predispose them to evolve along common lines, particularly in allopatry. On the other hand, selective regimes 75 may be less predictable, or even disruptive, in the case of character displacement. Here, we investigated the evolution of speciating parasites by sequencing genomes of two economically important mite species that specialize on honey bees (Apis sp.): Varroa jacobsoni and Varroa destructor. The honey bee colony, which is a densely packed community of genetically similar individuals, can host many 80 diseases and parasites. Brood parasitism by ectoparasitic mites has evolved in all Apis species except for the western honey bee (A. mellifera), the only lineage nonnative to Asia and allopatric with the others 17. The two Varroa sister species are morphologically very similar, and until recently were thought to be the same species 18 (Figure 1A). They originally parasitized the same host, the widespread 85 eastern honey bee (Apis cerana) (Figure 1B). However, as western honey bees were brought into contact with the eastern honey bee in Asia and Oceania, both mites switched repeatedly to this novel host 19,20 (Figure 1C). Since switching hosts to A. mellifera at the turn of the 20th century, V. destructor has spread worldwide, becoming a major driver of global honey bee declines 21,22. V. jacobsoni evolved an 90 ability to parasitize A. mellifera in the past two decades 20,23, but its potential for worldwide spread is not yet known. While both mite species represent potential or actual threats to the survival of A. mellifera, they are tolerated by A. cerana. Both the mites and the original bee host A. cerana have a range of extensively studied adaptations and counter95 adaptations 24–29. In response to Varroa, the eastern honey bee has strategies to tolerate the parasite and to limit its propagation 25. Specifically, the eastern honey bee has impressive host defence traits, such as worker brood incompatibility 26,27, higher toxicity response to Varroa saliva protein than A. mellifera 28,30 and hygienic behavior 31, which stop the mite reproductive cycle and reduce the infection in the 100 brood 25. Western honey been lack these defensive mechanisms, and hives collapse quickly after infestation by Varroa. Figure 1: V. destructor and V. jacobsoni are morphologically similar sister species 105 originally parasitizing the Asian honey bee (A. cerana). These species were recognized by Anderson and Trueman 2000 by quantitative morphometric and genetic data. The morphological differences are slight and rely mostly on body size, as shown on the 3D surface model comparison of two fully sclerotized females (A). Both species parasitize A. cerana as their original host and can coexist in sympatry, including at the same apiary , 110 but are mainly parapatric 47 (B). With the introduction of A. mellifera, V. destructor can be found on this novel host throughout the original V. jacobsoni range. V. jacobsoni has also extended its range into Papua New Guinea on A. cerana followed by a shift to A. mellifera, where V. destructor is currently absent 20 (C).

On the part of V. destructor, which has been most extensively studied as a result of its devastating agricultural effects, behavioral and chemical ecology studies have unraveled highly adaptive strategies to find, select, feed and reproduce on its host 32 . For example, gene expression changes in the chemosensory organs on the female mite legs with the need to find a suitable host within the host colony 120 and detect bee pheromones 33 . Mites also use chemical perception to identify the exact short time window of synchronized reproduction with the host larval development in a protected capped cell 32 . Another crucial part in Varroa parasitic lifestyle is the ability to feed on its host to respond to energy demand linked to egg production and survival 34 . The behavioral and metabolic processes involved have 125 recently been majorly updated as mites assumed hematophagous were shown feeding on fat tissues with the detection of host fat body proteins and transcripts in Varroa 35,36 , microscopic observation 37 , and chemical analysis of parasite waste excretion 38 . While many mite adaptations have already been discovered 32,39 , the recent use of 'omics approaches has quickly broadened our understanding of the 130 mite's hidden specialization to host 29,36,40,41 .
Here, we ask how both mite species have arrived at this equilibrium in A.
cerana -did they follow similar or different evolutionary paths? On one hand, the two species naturally occur in parapatry ( Figure 1B), which is expected to result from divergent local adaptation with the distribution of different A. cerana lineages 135 42,43 . In actuality, what drives the parapatric distribution is unclear. V. destructor, which was historically absent from Southeast Asia, now occurs there on its novel host A. mellifera in sympatry with V. jacobsoni, which does not parasitize A. mellifera in this region ( Figure 1C). Thus, it is possible that the current geographical distribution resulted from secondary contact following allopatric speciation as 140 shown in other closely related species like the alpine rock-jasmine 44 ; and similar evolutionary trajectories and competitive geographic exclusion like in the case of chipmunks (Eutamias dorsalis and E. umbrinus) 45 . This view is further supported by the fact that both species have the genetic capacity to shift to A. mellifera as a novel host, suggesting a level of physiological similarity. We tested these alternative 145 hypotheses (parallel vs. divergent evolution) by generating high-quality genomes of both species and examining them for signatures of adaptive evolution. We also compared them to the recently published genome of a distantly related ectoparasitic mite Tropilaelaps mercedesae, to see whether the patterns appear universal at different timescales.

Materials and Methods
Varroa mite sample collection V. destructor and V. jacobsoni are haplodiploids. Haploid males are only found inside capped brood cells after initiation of reproduction while diploid females can 155 be found in the brood and on adult workers. All V. destructor samples were collected from A. mellifera colonies of OIST experimental apiary in Okinawa, Japan (26.477 N 127.835 E) in August-September 2016. One mature male V. destructor was collected within a worker cell at the red-eye pupa developmental stage using a brush and kept in absolute ethanol at -80०C. Mature adult females were collected 160 from adult honey bee workers. In order to obtain a large number of mites infecting adult workers, we modified the standard powdered sugar method 48 for an entire hive as follows: one hive super was placed on top of a white collecting tray, containing three to four frames with honey bee workers and no queen; approximately 500g of powdered sugar was then applied using sieve filter powder 165 strainer; the sugar and honey bees were shaken in the tray for 2 min and then separated using a sieve filter, placing the sugared-workers directly within hive. This process was repeated on adjacent colonies. Varroa mites were separated from their hosts and trapped in the icing sugar. Following a water rinse on a gauze mesh, a total of 1,207 alive Varroa females were snap frozen at -80०C until laboratory 170 processing (inactive, sluggish or dead mites were discarded). One fully sclerotized mature female V. destructor (VDOKI-01) collected from the same apiary in reproductive phase was processed for X-ray microtomography (micro-CT). Micro-CT outputs were focused on the Varroa idiosoma ventral and dorsal views, and whole body in profile to show the mite anterior and gnathosoma parts. 175 Samples of V. jacobsoni were obtained from its first detection survey where it was found reproducing in the western honey bee 20  Prior to assembly paired-end reads were de-duplicated using ParDRe 49 , and 195 host and microbial contamination was filtered using DeConSeq (parameters: -i 97c 95) 50 , screening against genomes of Apis cerana, Apis mellifera, as well as bacterial genomes and genomes of bee diseases and parasites. The cleaned reads were then merged using PEAR (parameters: --min-overlap 10 --memory 80G -threads 10 -n 200 -m 600 -p 0.0001) 51 and assembled using Newbler (v. 2.6)
Scaffolding was then performed using Phase Genomics' Proximo Hi-C scaffolding pipeline 57 , which builds on the LACHESIS method 58 .
To estimate genome characteristics, K-mers from cleaned reads were counted using jellyfish version 2. Annotation for both genomes was carried out by NCBI using NCBI's automated genome annotation pipeline. This pipeline takes advantage of speciesspecific RNA-seq data, as well as extensive protein homology data stored in 230 GenBank 63 .

Prediction of gene orthology, selection, and duplication in Varroa mites
The Orthonome pipeline (http://www.orthonome.com) 64 was applied to predict orthologues (pairwise and 1:1 orthogroups), inparalogues and lineage-235 specific gene expansions from the two parasitic Varroa mites' annotated genomes and four Acari annotated genomes available from NCBI. The four outgroups were the Asian honey bee ectoparasitic mite Tropilaelaps mercedesae 65  Additionally, tandem duplications in the three parasitic mites were identified from inparalogues in the immediate vicinity of the orthologues.

245
Gene ontology (GO) enrichment analysis of the positively selected genes in honey bee parasitic mites was carried out using the GOstats R package 70 . The same process was applied for gene duplication in Varroa mites. Biological processes associated with those GO terms were summarized and visualized using REVIGO 71 (http://revigo.irb.hr). The semantic similarity threshold used was C = 0.5 (i.e.,

An improved de novo assembly at chromosomal level for V. destructor
We generated 28.6Gb (333 ± S.D. 56 bp) and 56.2 Gb (409 ± S.D. 61 bp) of deduplicated, decontaminated and merged Illumina reads for V. jacobsoni and V.
destructor, corresponding to a coverage of 80× and 152×, respectively.

290
Through the k-mer analysis the genome size estimates ranged from 369.02 to 369.55 Mbp using k = 42 ( Figure S1). The overall final Vdes_3.0 assembly size was close to that estimate reaching 368.94 Mbp. The final size represented a significant improvement of 25.4% in comparison to the first release assembly Vdes_1.0 87 and 11.2% to the later corrected into Vdes_2.0 (Table 1). Not only the overall assembly 295 size, but the contiguity of Vdes_3.0 was highly and significantly improved compared to the previous released assemblies. According to the N50 statistics, more than 50% of the genome was in scaffolds of 58.5 Mbp or more for Vdes_3.0 while the Vdes_2.0 N50 was around 0.1 Mbp. Comparison of genomic contiguity statistics with other invertebrate assemblies revealed that Vdes_3.0 currently has 300 the best reported scaffold N50 in arthropods and best contig N50 for Mesostigmata as reported by 88   The first draft genome for the mite V. jacobsoni, a novel parasite of A. mellifera 315 For V. jacobsoni, the assembly was performed using whole genome shotgun sequence from a Hiseq 2500. Prediction of V. jacobsoni genome size by k-mer frequency was slightly lower than V. destructor with an estimated model of 365.13 Mbp ( Figure S1). The accuracy and heterozygosity were estimated to be 99.6% and 0.06%, respectively. Similarly, for its sister species, the draft genome for V. Given that Cox1 is a standard marker used to identify Varroa mites "haplogroup", the region was extracted and realigned with described sequences from native and

Divergent patterns of selection and gene expansion among three honey bee ectoparasitic mites
We detected a total of 234 and 225 genes under positive selection for V.  Table   435 S7). In V. destructor, GO terms were mainly associated with 1) the regulation of membrane depolarization, 2) retinal cone cell development and myofibroblast differentiation, 3) mRNA cis splicing via spliceosome, 4) cellular response to alcohol (i.e. any process that results in a change or activity of a cell) and more detailed in Figure S3A. In V. jacobsoni, GO terms related to 1) protein processing involved in 440 targeting to mitochondrion, 2) vitamin K metabolic process, 3) response to pH, 4) gonad development l and more processes detailed in Figure S3B.  genes duplicated in V. jacobsoni were involved in biological processes of the skeletal myofibril assembly, Golgi calcium ion transport, striated muscle contraction ( Figure   5B, Supplementary Table S11). mites. Not only were most of the GO terms species-specific, but they were also responsible for non-overlapping categories of biological processes. Neither gene duplication nor selection analysis suggests substantial degrees of parallel evolution in the mites.  Table 2). Interestingly, both species have genes under positive selection that are known to be involved in acaricide resistance. These   (Table 2).

510
Both positive selection and duplication patterns along V. destructor and V. jacobsoni genomes and related GO terms showed that the two sister species underwent dissimilar evolutionary trajectories.

Stress from in-hive temperature
Heat Shock 70 kDa protein-like 93 6

525
In this study, we developed high-quality de novo reference genomes and annotations for economically important honey bee parasites V. destructor and V.
jacobsoni. Because Varroa mites are, to a large extent, responsible for the global honey bee crisis, and research efforts worldwide target these species, these genomic resources respond to a pressing need 29 . Our analysis revealed that the albeit at variable levels of infestation 111 , suggesting that encounters between the two species are indeed likely. If the divergent evolutionary between the two Varroa species resulted from competition, we predict that long-term coexistence within a hive of A. cerana should be unlikely, a hypothesis that we hope will be tested in the future in the zone of sympatry.

575
In addition, or as an alternative to interspecific competition, it is possible that the divergence between V. destructor and V. jacobsoni was driven by local adaptation to different A. cerana populations. For instance, V. destructor was found to infect the "Northern cerana", "Himalayan cerana" and the "Indo-Chinese cerana" morpho-clusters identified by 43 whereas V. jacobsoni was reported on the "Indo-

580
Chinese cerana" and the "Indo-Malayan cerana". These A. cerana populations differ in their nesting behaviour, such as the number of combs or even in their body size 86 . In China alone, A. cerana populations were shown to be as genetically divergent for more than 300,000 to 500,000 years (earlier than the recognized level of subspecies for A. mellifera) 112 . In addition to genetic differences between the hosts, 585 the climactic environments across Asia are highly variable and could result in local adaptation, though, as mentioned earlier, this did not prevent range expansions after switching hosts ( Figure 1) Such knowledge combined with the progress in novel manipulation tools, such as RNAi 94,116 , as well as more experimental tools like gene drive 117 , may provide alternative ways of specifically controlling mite pests without harming honey bees or other insects.

635
(OIST sequencing center) for leading and processing the PacBio RSII sequencing.
We would also like to thank Francisco Hita Garcia (Biodiversity and Biocomplexity OIST) for his support in the preparation and scanning of specimens using the micro-CT. This work was supported by the Okinawa Institute of Science and Technology Graduate University, as well as grants from the Japan Society for the

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable. 670

Competing interests
IL and STS are employees and shareholders of Phase Genomics, a company commercializing proximity-ligation technology.  Table 2: Genes undergoing positive selection or duplication, which were previously implicated in tolerance to external stress sources and stimuli.