Mites (Varroa destructor) were collected in four honey bee colonies (Apis mellifera), from the experimental apiary in Jiangxi Agricultural University. All colonies were maintained within one apiary. The colonies were 1 year old and had been treated with Amitraz during the summer. Twenty mites were collected from each of four honey bee colonies using the powdered sugar method in Spring (13). The body surface of the mites was rinsed with distilled water and twenty mites of each colony were pooled for RNA extraction using TRIzol. One Pacbio Iso-Seq library was prepared for each RNA pool to sequence in a SMRT cell. In total, four Iso-Seq libraries were prepared and sequenced in four SMRT cells.

The long reads were aligned to the V. destructor genome assembly (version V_des 3.0) using Minimap2 with default transcript parameters (14). Single nucleotide variants (SNVs) were identified using Longshot with minimal coverage of 20 (15). The SNVs were further annotated with SNPeff (16). SNVs along the contigs were analyzed with Pearson's Chi-Squared test. To analyze the haplotypes of the transcripts, the reads were aligned to the protein coding sequences of the V. destructor genome using Minimap2 with default parameters. The SNVs shared in all four colonies were concatenated to form haplotypes using SAM4WebLogo and summarized using Whatshap (15, 17, 18).

RNA reads that could not be aligned to the mite genome were retrieved using Samtools and assigned to microbes using Kraken2 (16, 19). The results were viewed with Krona (20). The reads assigned to Deformed wing virus (strain A) were extracted and aligned to a reference DWV genome (NCBI Genome version GCF000852585.1) with Minimap2. SNVs were predicted using Longshot and the nucleotides at SNV positions were clipped with SAM4WebLogo (13, 17, 21), which were concatenated to form a haplotype for each long read. The distribution of haplotypes among the colonies was analyzed with Pearson's Chi-Squared test, using R (22). Additionally, the reads were aligned to DWV-A (GCF000852585.1), DWV-B (AY251269.2) and DWV-C (CEND01000001.1) in parallel to quantify the relative abundance of DWV variants associated with the mites using Minimap2 (13, 23, 24).

To further study the competition of DWV haplotypes in mites and bees, a controlled mite infestation assay was performed. Second-instar larvae were collected from honey bee frames randomly selected from the colonies in the apiary. The larvae were artificially reared to the pupal stage in 48-well microtiter plates maintained in an incubator, with a temperature and humidity setting of 34°C, 70% RH, respectively (25). V. destructor mites were collected from capped cells to minimize any physical damage. During the pupation stage, one mite was transplanted onto each pupa to form a paired infection group (N = 10). Pupae without mites served as a control group to assess mortality due to artificial rearing (N = 10). Pupae, and their associated mites, were collected 4 days after infestation (N_pupa = 5, N_mite = 5), as well as the emerged adults (N_pupa = 5, N_mite = 5). Total RNA was extracted from individual bees and mites with TransZol Up Plus RNA kit (Transgen). cDNA was synthesized using PrimeScript^TM RT reagent Kit (Takara), with 1 μl gDNA Eraser L, 2 μl 5 × gDNA Eraser buffer, 5 μl RNase Free H₂O and 2 μl total RNA incubating at 42°C for 2 min. Additionally, 1 μl PrimeScript RT Enzyme Mix I, 4 μl 5 × PrimeScript Buffer, 1 μl RT Primer Mix and 4 μl RNase Free H₂O were added to the reaction for 37°C, 15 min, 85°C 5 min. The cDNA was stored at −20°C.

Based on the above SNV analysis of the long reads, the PCR fragments containing the highly dense and variable loci at the 218, 325, 555, and 567th nucleotides were amplified in both mites and bees using customized primers (DWV_h_5': ACGCGCGCGATAATGAGT, DWV_h_3′: GATCTCTGGTTTTGCCTGCAC). A 20 μl amplification reaction system was as follows: 1 μl 5′ primer, 1 μl 3′ primer, 10 μl 2x Tag PCR StarMix, 3 μl DEPC water, 5 μl RNA template. The reaction program consisted of 95°C 5 min, 95°C 10 s, 58°C 30 s, for 40 cycles. The 443 bp PCR product was purified and recovered from an agarose gel using a column DNA back kit (TIANDZ). PCR products were sequenced in both directions using the Illumina Nova 6,000 platform. These paired reads were joined and then aligned to the DWV CDS region using Minimap2 (13). The nucleotides at SNV positions (at 218, 325, 555, and 567th nucleotides) was clipped with SAM4WebLogo to form haplotypes (17). The distribution of haplotypes were analyzed with Pearson's Chi-square test, using R (22). The transmission of haplotypes was visualized using Circlize package in R (26).

On average, over 75 million long reads were generated from a SMRT cell, with a mean length of 3.5Kbp (over 150x of the mite transcriptome) (Supplementary Table 1). In total, 4,633, 5,683, 4,214, and 4,803 biallelic SNVs were identified from four honey bee colonies, out of which 892 SNVs were shared (Figure 1). The number of synonymous and non-synonymous SNVs was not significantly different among the four honey bee colonies (Pearson's Chi-square test, P > 0.05). On average, 5.25 SNVs were phased (assigned variants to the maternal or paternal chromosome) in a haplotype block (Supplementary Table 2). The mitochondrial genes cytochrome c oxidase subunit I (NP_758874.1), cytochrome c oxidase subunit III (NP_758878.1), and cytochrome b (NP_758884.1) were highly expressed. However, SNVs were not found in either of the above genes, nor for the mitochondrial gene cytochrome c (XP_022658338.1). The transcript activating signal cointegrator 1 complex subunit 3-like (ASCC3, XP_022662619.1) showed the highest number of SNVs (146 SNVs) in all four colonies. As 20 mites were pooled, a maximum of 40 haplotypes could be identified from this gene in each colony. In theory, 6 biallelic SNVs (2⁶ = 64) were sufficient to distinguish between the 40 haplotypes. Two SNVs at the 3' end of the transcript XP_022662619.1, with the highest coverage, were used as an anchor to construct haplotypes. One SNV was added at a time, until the observed number of haplotypes was less than the expected ones. In total, 30 haplotypes were identified using 5 SNVs and 12 haplotypes were shared among the colonies (Supplementary Figure 1). Surprisingly, a dominant haplotype (VD-1) was found in all four colonies (Figure 2A).

On average, 525,922 reads aligned to microbes associated with mites. Among the DWV variants, 95.41 ± 2.02% (Mean ± SD) reads were assigned to DWV-A, which was orders of magnitude higher than DWV-B (4.58 ± 2.01%) and DWV-C (0.02 ± 0.01%). In total, 381 SNVs were identified from DWV, among which 11 SNVs were commonly shared among the four colonies (at 218, 325, 555, 567, 615, 693, 780, 1809, 1938, 1977, and 8676^th nucleotide, respectively) and 266 SNVs were colony specific. The proportion of shared SNVs was significantly higher in mites compared with DWV along the genome (Pearson's Chi-square test, P < 0.05). Three commonly shared SNVs at the 1809, 1938, and 1977^th nucleotides were used to construct the haplotypes. In total, 11 haplotypes were identified, and 3 haplotypes dominated over 70% of total haplotype counts in all four colonies (Figure 2B; Supplementary Figure 2).

To determine whether the dominant DWV haplotypes in bees originated from high virus abundance in mites, or by viral competition during proliferation, a controlled mite infestation assay was performed. As noted above, for this experiment larvae were reared in a laboratory incubatory until pupation. All larvae survived this incubation. Then, each experimental pupal sample was paired with a Varroa mite. In the control group (no mite), DWV was not detected. In the treatment groups, 2.2 million PCR fragments were sequenced. The 218, 325, 555, and 567^th nucleotides were identified as having the highest diversity, on average, in the DWV genome. Surprisingly, two haplotypes accounted for 95% of DWV in both the mites and the bees. By analyzing the paired pupae and mite, the relative abundance of the DWV haplotypes in mites and bees was consistent, a result that deviated significantly from random (Fisher's exact test, P < 0.001) (Figure 3).