Divergent evolutionary trajectories following speciation in two ectoparasitic honey bee mites

Abstract

Multispecies host-parasite evolution is common, but how parasites evolve after speciating remains poorly understood. Shared evolutionary history and physiology may propel species along similar evolutionary trajectories whereas pursuing different strategies can reduce competition. We test these scenarios in the economically important association between honey bees and ectoparasitic mites by sequencing the genomes of the sister mite species Varroa destructor and Varroa jacobsoni. These genomes were closely related, with 99.7% sequence identity. Among the 9,628 orthologous genes, 4.8% showed signs of positive selection in at least one species. Divergent selective trajectories were discovered in conserved chemosensory gene families (IGR, SNMP), and Halloween genes (CYP) involved in moulting and reproduction. However, there was little overlap in these gene sets and associated GO terms, indicating different selective regimes operating on each of the parasites. Based on our findings, we suggest that species-specific strategies may be needed to combat evolving parasite communities.

Keywords: Coevolution; Comparative genomics.

Fig. 1

V. destructor and V. jacobsoni are morphologically similar sister species originally parasitizing the eastern honey bee (A. cerana). These species were recognized upon the basis of quantitative morphometric and genetic data. Morphological differences are restricted to differences in body size and shape, as shown in a 3D surface models comparison of two fully sclerotized females. b Both mites parasitize A. cerana as their original host and can coexist in sympatry, occasionally even at the same apiary^,, but are mainly parapatric. c 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

Fig. 2

Inter and intraspecies mtDNA variability across the genus Varroa and host-switched lineages. Different mitochondrial markers used in the literature to discriminate the four Varroa species (and unresolved one) based on partial COX1 (a, b) and intraspecies lineages using larger COX1, ATP6, COX3 and CYTB regions^, (a, c). The two cryptic species V. destructor (red) and V. jacobsoni (blue) are genetically divergent as shown by the unrooted phylogenetic tree of Varroa mite partial COX1 sequences (b). Variability of V. jacobsoni is higher than that of any other Varroa species and the reference genome (blue arrow) corresponded to one of the previously detected haplotypes switching on A. mellifera in Papua New Guinea (b). Several haplotypes from the Korean and Japanese lineages successfully jumped on A. mellifera (bold)^, but only K1-1/K1-2 is quasi-cosmopolitan and can even be retrieved in the native range of the Japanese lineage as illustrated by its presence in Okinawa (red arrow) (c). Ac identified from Apis cerana, Am identified from Apis mellifera

Fig. 3

Positive selection and gene duplication of orthologous genes in parasitic and free-living Acari. For the honey bee mite parasites, host ranges are shown as checkboxes and individual records had to meet the following conditions: (1) Mites had to be found reproducing or observed in bee colonies in independent surveys, and (2) the identity of cryptic Tropilaelaps and Varroa species was confirmed by molecular markers (e.g., COX1 barcoding^,). The number of orthologous genes in each lineage is circled at the nodes. The number of positively selected genes for each branch is shown in green, with the number of species-specific genes shown in purple. Gene duplications within honey bee parasitic mites are shown in boxes: white squares = shared genes, black squares = Varroa-specific genes by, red and blue squares = V. destructor and V. jacobsoni specific genes, respectively. The two Varroa mites have similar host ranges, through V. jacobsoni also occurs on A. nigrocincta, a close relative of A. cerana, possibly as a result of a recent secondary host shift^,. Despite sharing the same ancestral host, the Varroa sister species show different patterns of adaptive molecular evolution and exhibit different gene expansion patterns, suggesting different evolutionary trajectories

Fig. 4

Genes and pathways under positive selection in the Varroa sister species. a Red and blue bars (in 5 kb windows) represent locations of genes in V. destructor and V. jacobsoni, respectively. V. jacobsoni data are mapped to V. destructor scaffold positions for comparison. Black asterisks indicated genes under positive selection, shared by both species (n = 12) and detected in the Varroa ancestral lineage prior to the split of the two species (n = 40). b Semantic space analysis of significantly enriched GO terms (BP, CC, and MF) over-represented in genes detected under positive selection for Varroa mites Bubble color indicates the species for which GO terms were enriched (all p-value < 0.05) and size indicates the frequency of the GO term found in the GOA database. There was little overlap in analyses at both gene (2.8%) and functional levels (0.8%), suggesting different selective pressures on the two sister species since they split

Fig. 5

Duplicated genes in V. destructor are found throughout the genome and are involved in different biological pathways. a Chromosomal location and similarity of duplicated genes (arrows) for V. destructor. Different genes and chromosomal regions underwent duplication in the two species. b Cluster analysis of significantly enriched GO terms for biological processes over-represented among duplicating genes in Varroa mites. Not only were most of the GO terms species-specific, but they also comprised non-overlapping categories of biological processes. Neither gene duplication nor selection analysis suggests a substantial degree of parallel evolution in these mites

Fig. 6

The largest IGR repertoire in honey bee parasites show divergent selective trajectories for host and environment chemosensing between V. destructor and V. jacobsoni. The phylogenetic tree was constructed using 153 amino acid sequences from Acari V. destructor (Vdes = genome and Var = protein sequences), V. jacobsoni (Vjac), T. mercedesae (Tmer), M. occidentalis (Mocc), I. scapularis (Isca), and T. urticae (Turt) which were aligned with MAFFT (see list in Supplementary Data 9). Best fit model computed for the tree is VT + R5 using IQ-TREE. Bootstrap values were estimated using an SH-like aLRT with 1000 and bootstraps over 95% are shown by a purple circle

Fig. 7

CYP repertoire in Varroa mites is broader than expected with two functionaly important Halloween genes shade and disembodied under positive selection in V. jacobsoni. The phylogenetic tree was constructed using 174 amino acid sequences from Acari V. destructor (Vdes = genome and Var = protein sequences), V. jacobsoni (Vjac), T. mercedesae (Tmer), M. occidentalis (Mocc), I. scapularis (Isca), and T. urticae (Turt). Additionally, sequences from the fruit fly Drosophila melanogaster (Dmel), mosquito Aedes aegypti (Aaeg) and the honey bee A. mellifera (Amel) and A. cerana (Acer) were downloaded form NCBI and aligned with MAFFT. Best fit model computed for the tree is LG + R5 using IQ-TREE. Bootstrap values were estimated using a SH-like aLRT with 1000 and bootstraps over 95% are shown by a purple circle. V. destructor was previously believed to have only disembodied (CYP302a1), shade (CYP314a1), and spook (CYP307a1) homologs of the seven Halloween genes (background in black)