Significant transcriptional changes in mature daughter Varroa 1 destructor mites during infestation of different developmental stages 2 of honeybees 3

22 BACKGROUND: Varroa destructor is considered a major cause of honeybee ( Apis 23 mellifera ) colony losses worldwide. Although V. destructor mites exhibit preference 24 behavior for certain honeybee lifecycle stages, the mechanism underlying host finding 25 and preference remains largely unknown. 26 RESULTS: By using a de novo transcriptome assembly strategy, we sequenced the 27 mature daughter V. destructor mite transcriptome during infestation of different stages 28 of honeybees (brood cells, newly emerged bees and adult bees). A total of 132,779 29 unigenes were obtained with an average length of 2,745 bp and N50 of 5,706 bp. About 30 63.1% of the transcriptome could be annotated based on sequence homology to the 31 predatory mite Metaseiulus occidentalis proteins. Expression analysis revealed that 32 mature daughter mites had distinct transcriptome profiles after infestation of different 33 honeybee stages, and that the majority of the differentially expressed genes (DEGs) of 34 mite infesting adult honeybees were down-regulated compared to that infesting the 35 sealed brood cells. Gene Ontology and KEGG pathway enrichment analyses showed 36 that a large number of DEGs were involved in cellular process and metabolic process, 37 suggesting that Varroa mites undergo metabolic adjustment to accommodate the 38 cellular, molecular and/or immune response of the honeybees. Interestingly, in adult 39 honeybees, some mite DEGs involved in neurotransmitter biosynthesis and transport 40 were identified and their levels of expression were validated by qPCR. 41 CONCLUSION: These results provide evidence for transcriptional reprogramming in 42 mature daughter Varroa mites during infestation of honeybees, which may be relevant 43 to understanding the mechanism underpinning adaptation and preference behavior of 44 these mites for honeybees. 45

its host. 14 However, there is a lack of data on the global dynamic transcriptome of 80 mature daughter Varroa mite during its transition from young bees to newly emerged 81 bees, to adult bees. 82 In this study, we performed a de novo transcriptome assembly and annotation of  Adult female V. destructor mites with the same age were equally obtained from 6 102 honeybee colonies with a high level of parasitism. These mites were delivered into 103 brood cells in non-infested colonies shortly before capping. Mature daughter mites were 104 collected from sealed brood cells in the comb (S1), newly emerged bees from brood 105 cells within one day (S2) and adult bees (emerged from the cell after 7 days, S3). For 106 S1 group, mature daughter Varroa mites were harvested from soon-to-emerge bees. For 107 S2 and S3, mother mites were removed from brood cells at pupal stage and the cell 108 sealed with melted beeswax. Each group contained 3 replicates in 3 non-109 infested colonies. Each replicate consists of 100 adult mites. All the mite samples were 110 frozen in liquid nitrogen and stored at -80℃ until RNA extraction. Using in-house Perl scripts, raw reads were filtered by removing the low-quality reads 137 and reads that contains adapter or poly-N. Also, Q30, GC-content and sequence 138 duplication level were calculated based on the clean reads in order to evaluate the 139 sequencing quality. Sequence contaminants from other sources, such as bacteria and 140 fungi, present in the samples were removed prior to functional annotation analysis. The   were identified using the DESeq R package based on the negative binomial 161 distribution. 19 The P values were adjusted using the Benjamini and Hochberg's 162 approach for controlling the false discovery rate (FDR). An adjusted P-value ＜0.05 163 along with at least two-fold change was used to identify significantly differential  (Table S1), and the qPCR was performed on LightCycler 480 (Roche 184 Diagnostics, Tokyo, Japan). The amplification reactions were performed with the 185 following conditions: 2 min at 95°C, 40 cycles of 95 °C for 5s, 60°C for 30s. The 186 experiment was repeated three times using three independently isolated RNA samples.

187
The glyceraldehyde-3-phosphatedehydrogenase (GAPDH) gene was used as a 188 reference, and relative expression levels were calculated using the 2 −ΔΔCT method.  Table 1, we acquired > 46 million 150bp paired-197 end-seq raw reads from each cDNA library. After eliminating adapters, ambiguous 198 nucleotides and low-quality sequences, > 44 million clean reads (Q20>95%) were 199 retained, which accumulated to > 6.75 Giga bases (Gb) read length with a GC

213
In order to identify the putative functions of the unigenes, BLAST programs (e value< 214 1.0E-5) was employed to search against public databases (Nr, Nt, KO, Swiss-Prot, 215 PFAM, GO and KOG), which were used for gene annotations. As shown in Table 3, 216 the results showed that 82,068 (61.8% of 132,779) unigenes were matched to one or 217 more databases. A total of 68,940 unigenes were found to have homologs in the NR 218 database with an e-value < the cutoff (e-value=1E-5). The e-value distribution analysis 219 of the hit unigenes showed that 52.7% of V. destructor unigenes have the highest 220 homology with an e-value cut-off < 1E-100 ( Fig. 2A). Likewise, the similarity 221 distribution showed that 56.9% of all the unigenes had a similarity > 80%, whereas 222 42.9% of unigenes had similarity that ranged from 40% to 80% and only 0.1% had 223 similarity below 40% (Fig. 2B). As anticipated, the top unigene hit was found in the while only 0.8% of contigs were fragmented (9 BUSCOs) and 0.4% were missing (4 232 BUSCOs) (Fig. 3). 233 All unigenes were aligned to the Cluster of Orthologus Groups (COG) database for 234 functional prediction and classification. A total of 44,672 unigenes were assigned to 235 appropriate COG clusters, which could be classified into 25 functional categories. As 236 shown in Fig. 4, 'General function prediction only' was the largest category (7,698 237 unigenes); followed by 'signal transduction mechanisms' (7,604 unigenes), and 238 'posttranslational modification, protein turnover, chaperones' (4,545 unigenes).

280
To better understand the biological regulatory mechanisms underlying mite infestation, 281 we performed GO annotation analysis of the DEGs identified in this study. Significantly 282 enriched GO terms were identified using an adjusted P-value based on hypergeometric 283 distribution. There were 60 significantly enriched GO terms between S2 and S1, and 284 49 between S3 and S1. However, there were no significantly enriched GO terms 285 between S3 and S2. In terms of S2 group vs. S1 group, the GO terms with the maximum 286 number of DEGs in biological process, cellular component and molecular function were 287 'cellular process', 'cell', and 'binding', respectively (Fig. 8A). When comparing S3 288 group and S1 group, the maximum number categories of DEGs enrichment in the GO 289 three categories also were 'cellular process', 'cell', and 'binding', respectively (Fig.   290   S3). 291 Next, we performed the KEGG pathway enrichment analysis based on these DEGs.

346
In our study, 155 DEGs were found between S2 and S3. Also,1,492 DEGs were 347 previously detected in mites collected from capped brood cells containing developing 348 bees ready to emerge compared to mites infesting adult bees. 16 The majority of the 349 DEGs in V. destructor mites infesting newly emerged bees and adult bees were down-350 regulated compared to mites infesting the sealed brood cells, suggesting that infestation 351 of the more mature honeybee stages seems to be associated with suppression rather than 352 activation of V. destructor genes. These results agree with a previous study, which 353 showed that mites have a preference to adult honeybees, but seems inconsistent with 354 the fact that more mites are found in the sealed brood cells than on adult bees. 25

355
There were a small number of genes whose expression was discordant between 356 mature daughter mites infesting the newly emerged bees and the adult bees. These 357 DEGs were involved in cellular and metabolic processes, and were largely up-regulated 358 between S2 vs. S3, suggesting increased dysregulation of biological functions and 359 processes that were already dysregulated in mites infesting newly emerged bees.

360
Preference for specific honeybee development stage can influence the mite's food 361 intake and metabolism. A previous study showed that V. destructor mites consume fat 362 body tissue rather than honeybee's hemolymph. 2 In our study, active fatty acid 363 metabolism was detected, and numerous lipid metabolic pathways were affected during 364 V. destructor infestation of different stages of honeybees.

365
The GO enrichment analysis showed that DEGs identified in the comparison 366 between S2 vs. S1 and S3 vs. S1 were strongly associated with metabolic process and 367 binding activities, suggesting that transcriptional changes observed in mites infesting 368 different development stages of honeybees are related to metabolic processes.

369
Oxidative phosphorylation is another pathway that was significantly affected, 370 indicating that energy metabolism is involved in mite-honeybee interaction. Compared 371 to S1, DEGs in both S2 and S3 were significantly involved in the proteasome pathway.         Table S4. KEGG pathway analysis of the DEGs identified between S2 and S1.  Table S5. KEGG pathway analysis of the DEGs identified between S3 and S1.