An Insight Into The MicroRNA Proles of An Ectoparasite Mite Varroa Destructor (Acari: Varroidae), The Primary Vector of Deformed Wing Virus (DWV) of The Honey Bee

The most important ectoparasite of the honeybee is the remarkably adaptive mite Varroa destructor. The Varroa mite is a competent vector of Deformed Wing Virus (DWV), the Varroa-virus complex is one of the main factors associated with elevated annual honey bee colony mortality. Micro-RNAs (miRNAs) are small, non-coding RNAs of ~22-24 nucleotides, produced by all plants, animals, and viruses that inuence biological processes through post-transcriptional regulation of gene expression. Knowledge of miRNAs and their functional role in mite biology remains limited. This study developed small RNA libraries from male and female V. destructor by utilizing Illumina's small RNA-Seq platform. A total of 101,913,208 and 91,904,732 small RNA reads (>18 nucleotides) from male and female mites were analyzed using the mirDeep2 algorithm. A conservative approach predicted a total of 306 miRNAs, of which 18 were upregulated and 13 were downregulated in female V. destructor compared to males. A qPCR assay validated the expression of selected differentially expressed female Varroa miRNAs. This dataset provides a list of potential miRNA targets involved in regulating vital Varroa biological processes and for deciphering new targets against Varroa and honey bee viruses they carry.


Introduction
Varroa destructor is considered as the most damaging ectoparasite of the European honey bee (Apis mellifera) in every country it has been introduced (Hristove et al., 2020). The signi cance of A. mellifera to world food production (Gallai et  with Varroa will not survive for more than 3 years without treatment or cultural intervention (Martin 1998;von Dooremalen et al., 2012). In the United States, it has been reported that honey bee colonies would have signi cantly reduced chances of survival unless beekeepers use chemical treatments (Webster and Delaplane 2001;Kulhanek et al., 2021). An accurate assessment of loss caused to honey bee colonies by Varroa infestation is a technical challenge. Still, it is safe to assume that Varroa infestation has resulted in the death of thousands of honeybee colonies and losses in billions of dollars (http://entnemdept.u .edu/creatures/misc/bees/varroa_mite.html).
Varroa foundresses vector Deformed Wing Virus (DWV) to immature honey bees when they feed. DWV is the major cause of honeybee colony collapse. Interestingly, DWV infection is maintained naturally in honey bees and rarely pathogenic in the absence of Varroa infestation (Beaurepair et al. 2019). Contrary to low levels of infection maintenance, inoculation of DWV via Varroa feeding signi cantly increases the viral load in the pupal stage (Gusachenko et al. 2020). Bees that are symptomatic for DWV emerge with wing deformities, have reduced weight, and a shortened lifespan or are killed prematurely by other members of the colony. Bees can be asymptomatic but have high viral load that leads to negative effects on lifespan, foraging and ight capability, behavioral maturation, and immunity ( Varroa and infection to honey bees, the underlying genetic mechanisms of vector competence (acquisition, maintenance, and transmission) needs further clari cation. A detailed understanding of vector competence for DWV by Varroa is required in order to develop intervention strategies to prevent Varroa infestation and Varroa-borne viral transmission to honey bees.
The control and prevention of Varroa is primarily based on chemical acaricides. The reliance on synthetic acaricides has resulted in the development of Varroa resistance (Higes et al., 2020; Rinkevich 2020). Varroa is highly e cient in vectoring honey bee viruses and plays a role in driving changes in virus distribution, prevalence, and virulence (Traynor et al. 2020). Clearly, current control strategies are not su cient; more effective and novel approaches need to be adopted to tackle this global problem. A novel microRNA approach has been investigated in this study to nd potential microRNAs (miRNAs) that regulate Varroa vector competence. MicroRNAs are single strands of non-coding RNA about ~ 22-25 nucleotides long derived from larger hairpin RNA precursors. involved in most physiological and pathological processes through post-transcriptional regulation of several target genes. MicroRNAS are being used as therapeutic targets of several diseases (Bartel et al. 2018). In mammals, miRNAs can regulate more than 30% of protein-coding genes and most cellular processes (Filipowicz et al. 2008). Information is limited on arthropod miRNAs and their role in arthropod physiology ( . Therefore, information about Varroa microRNAs could play a signi cant role in understanding how Varroa miRNAs participate in colonizing DWV and its transmission to honeybees.
Our aim was to predict the expression patterns of conserved and novel miRNAs in male and female Varroa. This study provides insight into the physiological roles of miRNAs in gene regulatory networks and their biological function in the ectoparasitic mite. Identi cation of novel miRNAs provide a platform for the fundamental understanding of the role of miRNAs in the transmission of Varroa-borne viruses and other biological pathways.

Materials And Methods
Varroa mite collection Varroa were collected from A. mellifera colonies maintained at the University of Southern Mississippi's Lake Thoreau Environmental Center in Hattiesburg, MS (31° 20' 54.476" N and 89° 25' 9.202" W) during August-October, 2019. Adult female mites were collected from adult honey bee workers using the sugar shake method (Gregorc et al. 2018). Male Varroa were collected from sealed pupal cells by removing the wax cap and positively identifying males in infested cells using microscopy. A total of 10 male and 20 female Varroa were collected separately in tubes containing Trizol for RNA extraction.

RNA extraction
RNA was extracted from male and female mites separately using the TRIZOL extraction method (Chozyzynski et al. 1995) with some modi cations to the original protocol. Brie y, mites were collected individually then pooled male and female mites in separate tubes. Mites were homogenized in 500 µl of Trizol using a plastic pestle. After homogenization, samples were mixed for 10 min at room temperature in a shaker followed by centrifugation at 15000xg for 10 mins at 4ºC. The supernatant was transferred to a new tube and incubated for 5 min at room temperature to permit complete dissociation of the nucleoproteins. One hundred µl chilled chloroform was added to the samples, mixed, and incubated for 10 mins at 4ºC. Samples were centrifuged at 15000xg for 15 mins to obtain an aqueous phase. The aqueous phase was transferred to a new tube and then 600 µl of isopropanol was added and kept overnight at -20ºC. The following day, samples were centrifuged for 15 mins at 4ºC followed by 70% ethanol wash of pellets. RNA pellets were dried appropriately, and RNA samples were resuspended in sterile water and quanti ed using Nanodrop.

Small RNA Sequencing
The RNA concentration of pooled Varroa females was 286.5 ng/ul (260/280 = 1.8) and for pooled Varroa males was 185.6 ng/ul (260/280 = 1.8). Small RNA libraries were prepared using the Illumina Truseq Kit following the manufacturer's instructions. Brie y, short adapter oligonucleotides were ligated to each end of the small RNAs in the samples. Individual cDNA copies were made with reverse transcriptase, and PCR was used to add sample-speci c barcodes and Illumina sequencing adapters. The nal concentration of all NGS libraries was determined using a Qubit uorometric assay. The cDNA fragment size of each library was assessed using a DNA 1000 high-sensitivity chip on an Agilent 2100 Bioanalyzer. Before sequencing, samples passed the essential quality control (QC) test. After puri cation by polyacrylamide gel electrophoresis, the sample libraries were pooled and sequenced on an Illumina Next Seq 500 (single end 36 bases) using TruSeq SBS kit v3 (Illumina) and protocols de ned by the manufacturer. Four small RNA libraries of pooled male and female Varroa mites were sequenced via Illumina small RNA highthroughput sequencing. RNA library preparation and indexing were performed by the University of Mississippi Medical Center's genomics core facility. Small RNA sequencing was performed on Illumina Next Seq 500 high output (300 cycles) as single-end sequences.

Bioinformatics analysis
The miRDeep2 software package (Friedlander et al. 2008), version 2.0.0.8, was used to process the sequencing data. For the novel miRNA prediction step, the reads from all the samples were combined. The mapper function of miRDeep2 trimmed the adapter sequences from the reads and converted the read les from fastq to fasta format. Reads shorter than 18 bases were discarded. Remaining reads were then mapped to the V.destructor reference genome (Techer et al., 2019) using the default miRDeep2 mapper function parameters. Reads that mapped to the genome were used to predict novel miRNAs. The Drosophila melanogaster genome was also provided as a reference genome and mapped reads were aligned to available microRNAs of D. melaongaster in miRbase (Version 22) and quanti ed. The output le included sequence and location of the possible miRNAs, the number of reads that mapped to it, and a score re ecting the likelihood that the predicted miRNA was not due to chance alone. The software aligned the reads to the reference genomes of D. melanogaster and V. destructor and looked for locations where potential miRNA reads accumulated. The regions immediately surrounding the mapped reads were examined for miRNA biogenesis features, including mature miRNAs, Star and precursor reads, and stemloop folding properties. Essentially, the miRDeep2 program models the miRNA biogenesis pathway, using a probabilistic algorithm to score compatibility of the position and frequency of Next-Gen Sequencing (NGS) reads with the secondary structure of the miRNA precursor. BEDtools (Quinlan and Hall, 2010) were used to determine the genomic origin of the precursors with the variant intersect.

Validation of differentially expressed miRNAs by qRT-PCR
Predicted mRNAs that were differentially expressed in small RNA sequencing were validated by qRT-PCR.
RNA samples from male and female Varroa were converted to cDNA, and the miRNA-speci c qRT-PCR reaction was performed for each sample. Mir-X miRNA qRT-PCR TB Green kit from Takara BIO (catalog # 638316) was used for cDNA synthesis and miRNA expression analysis. Conditions used for qRT-PCR were initial denaturation 95º C for 10 mins, then 40 cycles of 95º C for 5 secs, 60º C for 20 secs.
Normalization, Differential expression (DE) and Statistical analysis of miRNA expression between male and female Varroa mites In silico differential expression analysis of predicted miRNAs were performed using the interactive web interface DeApp (Li and Andrade, 2017). DeApp is a web-based, graphical interface developed in R with the shiny package (web application framework for R). Low expression genetic features were removed after alignment if the counts per million (CPM) value was ≤ 1 in less than two samples. Normalization of samples and a multidimensional scaling (MDS) plot have been provided as supplementary data (Fig. S2) and includes details about sample distribution after ltering out the low expression genomic features. Differential expression (DE) analysis was performed by edgeR, the false discovery rate (FDR) adjusted pvalue was set at 0.05, and minimum fold change was set at 1.5. In the end, this interface displays a dispersion plot showing the results of overall DE analysis along with statistical signi cance (p-value, FDR adjusted p-value, Table 1) with a volcano plot (Fig. 3A) corresponding to the speci ed parameters and cutoff values.  . Targets common in all three programs were further considered. In silico target prediction provided a high number of false positives, but cross-species comparisons and combinatorial effects led to a reduction in this number (Min and Yoon, 2010). Lists of target genes were functionally characterized using the STRING webserver (Franceschini et al. 2013). The networks of target genes and the KEGG pathways signi cantly enriched for target genes were extracted using the STRING output (Table 3) The miRNA pro les of Varroa Fig. 2A shows the basic hair-pin-loop structure of a microRNA and other parameters (dicer cut overhangs, total read count, mature read count, loop read count, total read count, randfold score, and total score) used to determine whether a hair-pin-loop structured RNA is a microRNA. In this study, a total of 306 microRNAs were predicted in V. destructor, and several had homologs present in D. melanogaster. The predicted miRNAs were categorized as high (n=50) and low con dence (n=80) (Fig. 2B, Table S1) based on standard criteria (Bartel et al., 2018). Eighteen of the predicted miRNAs were upregulated, whereas 13 were downregulated in female Varroa compared to male Varroa (Fig. 3A-B; Table 1). Four of the predicted mature microRNAs, nDS_019211455.1_10989, nDS_019211455.1_7097, nDS_019211459.1_37116 and, nDS_019211459.1_24860 were detected in female but not male Varroa, and two of the predicted miRNAs, nDS_019211455.1_10989 and nDS_019211455.1_16072 were detected in male but not female Varroa.
Supplementary data (Table S1) includes the consensus mature, star, the total read counts, and the miRDeep2 score and probability for every miRDeep2-predicted miRNA with a miRDeep2 score.
In silico differential expression (DE) analysis of miRNAs predicted in the Varroa mites In silico analysis predicted a total of 306 microRNAs; fty were of high con dence, and eighty were of low con dence (Fig. 2, Table S1). Eighteen of the predicted microRNAs were upregulated, while 13 of them were down-regulated in females compared to male Varroa (Fig. 3-4 and Table 1).

Prediction of target genes and gene ontology (GO) and functional enrichment analyses in the target network
The STRING web analysis (Fig. 6A) showed that target proteins for 31 predicted miRNAs (18 upregulated and 13 downregulated) have more interactions than what would be expected for a random set of proteins of similar size sampled from the V. destructor genome (number of nodes = 42, number of edges = 97, average node degree = 4.62, average local clustering coe cient = 0.476, expected number of edges = 25, PPI enrichment p-value < 1.0e-16). Such enrichment indicates that the proteins are at least partially biologically connected as a group. To minimize the number of false-positive targets, we opted only for those targets which were predicted by all three miRNA target programs (TargetSpy, MIRANDA, and PITA). The predicted miRNAs are of high con dence, and the main KEGG pathways involved were predicted as oxidative phosphorylation (12 out of 80), endocytosis (4 out of 116), protein processing in the endoplasmic reticulum (2 out of 109), and other metabolic pathways (14 out of 810). Many target genes were predicted for the differentially expressed miRNAs in the small RNA-Seq data (Fig. 6B)  Genes controlled by the conserved miRNAs (D. melanogaster, V. destructor) identi ed in this study regulate development, metamorphosis, proliferation, apoptosis, reproduction (Aedes aegypti), insecticide resistance and also antiviral immune responses (Table 2). Interestingly, analysis of GO terms related to genes regulated by the shared miRNAs revealed many speci c processes known to be regulated by the same miRNAs in D. melanogaster, such as developmental processes, metamorphosis, immune response which reinforces its conservative role ( Fig. 7; Table 2).

Validation of predicted microRNAs of female Varroa by qRT-PCR
The expression of 13 predicted differentially expressed miRNAs were validated using qRT-PCR assay (Fig.  8). The qRT-PCR patterns of differentially expressed microRNAs matched the NGS patterns in the majority of evaluated miRNAs. Inconsistencies in the patterns of NGS data and qRT-PCR data were found for nDS_0192211459.1_37129, nDS-miR-6-3p, nDS-miR-4943-3p, nDS_019211456.1_17707, and nDS_019211457.1_26197. These differences could be due to different methodologies of quantifying miRNA expression (Saldana et al., 2017).  (Kloosterman and Plasterk 2006). However, little is known about the expression of miRNAs in parasitic mites such as V. destructor. Therefore, identifying and understanding the expression pro les of miRNAs in male and female Varroa will expand our fundamental knowledge of mite miRNAs and provide insight into the development and vector competence of Varroa. A conservative in silico approach identi ed 306 microRNAs, among which 80 were of low con dence and 50 of high con dence. Among these high con dence microRNAs 18 were up-and 13 were down regulated in female compared to male Varroa. Our qRT-PCR assay con rmed the up-regulation of Varroa miRNAs such as vde-miR-87-3p, vdebantam-3p, vde-miR-375-3p, and vde-miR-34-5p, which have known roles in DENV replication, CHIKV infection, and immune response during DENV-2 infection in mosquitoes (Hussain et  immune response, and reproduction (Fig. 4A, Table 2) in other arthropods. The new and novel microRNAs predicted in this study remain to be characterized using miRNA inhibitory assays.
In another study, it was most abundantly expressed in both pupal and adult male and female mosquitoes, indicating its functional importance (Feng et al., 2018) and its role in Varroa development.
Another conserved miRNA miR-8-3p is signi cantly up-regulated in Aedes aegypti during pupation and has the highest expression at the mid-pupal period. The miR-8-3p is also signi cantly up-regulated in the fat body. Bryant et al., (2010) showed the up-regulation of miR-8-3p in the fat body of a blood-fed female mosquito and suggested a potential regulatory role in Ae. aegypti reproduction. Contrary to Ae. aegypti, miR-8-3p is abundantly expressed in Anopheles stephensi developmental stages (Feng et al., 2018) and found to be equally expressed in uninfected or infected Aedes albopictus saliva upon infection by CHIKV (Maharaj et al., 2015). In Ae. aegypti, miR-8 has been validated to target SWIM (secreted winglessinteracting molecule), thereby regulating reproductive events (Lucas et al.,2015). Additionally, miR-8 shows cell-type speci city and is expressed in S2 cells (a cell line derived from Drosophila melanogaster embryos). Its temporal expression is less restricted in Drosophila, and has been observed across all developmental stages occasionally with signi cant variation in expression (Jin et al., 2012). On the basis of above-mentioned studies, the conserved nature of miR-8-3p also indicates its possible role in Varroa development, reproduction, and virus infection and requires further investigation. Another miRNA predicted in this study that is common to Drosophila is miR-34-5p, which is more pronounced in female midguts in Anopheles gambiae (Winter et al., 2007). In contrast, it is down-regulated in Drosophila during the metamorphosis period (pupa to adult stage transition) (Sempere et al., 2003). Interestingly, in Anopheles gambiae, miR-34-5p expression was found down-regulated in the midgut upon Plasmodium falciparum infection (Dennison et al., 2015). The miR-34-5p is suggested to contribute to anti-pathogen and immune responses during DENV-2 infection in Ae. albopictus (Liu et al., 2015). Conserved nature of this microRNA indicates its possible role towards development regulation in Varroa but requires further studies. Our in silico data showed upregulation, and our qRT-PCR data has shown its upregulation by ~ 5000 fold in DWV-B infected varroa females suggesting its The predicted miRNAs in the Varroa small RNA-seq provide a list of conserved miRNA targets involved in regulating key biological processes. As discussed earlier, conserved Drosophila homolog microRNAs predicted in the Varroa (Fig. 4A, Table 2), including vde-miR-87-3p, vde-bantam-3p, vde-miR-375-3p and, vde-miR-34-5p. These miRNAs have been shown to play a role in inhibiting or replicating dengue (DENV) and chikungunya virus (CHIKV) in mosquito species (Sempere et al., 2003, Maharaj et al., 2015, Hussain et al., 2013. The expressions of these miRNAs were signi cantly up-regulated (100 − 20,000 fold) in DWV-B infected Varroa female mites (Fig. 4B). These results indicate a functional role of these miRNAs in virus load and warrant a detailed study to characterize these miRNAs in Varroa-borne virus transmission studies.
Some of the predicted mature microRNAs such as nDS_019211455.1_10989, nDS_019211455.1_7097, nDS_019211459.1_37116 and, nDS_019211459.1_24860 were female speci c, detected only in the Varroa females but not in males, and 2 of the predicted miRNAs, nDS_019211455.1_10989 and nDS_019211455.1_16072 were male-speci c, detected in males but not in females. These male and female-speci c microRNAs in Varroa need to be dissected on the basis of their functional roles. A previous study in a nematode Ascaris suum has predicted the role of gender speci c miRNAs as a set of elongation factors, heat shock proteins, and growth factors essential for the development of the organism (Xu et al., 2012). Moreover, sperm proteins and sperm cell motility proteins were found as the targets of the male-speci c miRNAs while ovarian message proteins were found as targets of femalespeci c miRNAs (Xu et al., 2012).
In this study, the size distributions of miRNAs in the Varroa male and female data sets established two peaks, with a predominant peak correlated to miRNAs of 24 nucleotides in size, consistent with genomewide identi cation of miRNA study of the Varroa (Fonseca et al., 2021). These miRNAs are major contributors to the total small RNAs from the Varroa males and females. Moreover, it has been suggested that miRNA might be actively involved in the vector competence of the Varroa-borne viruses.
Using publicly available Varroa destructor genome data, we annotated and evaluated possible targets and infer functions of those female Varroa putative miRNAs which showed in silico up-regulation or down-regulation relative to male. The most conserved miRNAs are involved in several biological, cellular, and development processes.
Main KEGG pathways predicted by STRING web analysis were oxidative phosphorylation, endocytosis, protein processing in the endoplasmic reticulum (ER), and other metabolic pathways. In the predicted KEGG pathways, 14 proteins out of 810 have been annotated for metabolic pathways, 12 out of 80 for oxidative phosphorylation, 4 out of 116 for endocytosis, and 2 out of 109 for protein processing in the endoplasmic reticulum (Table 3). Not surprisingly, this data suggests the signi cance of these biological processes. Proteins involved in endocytosis and protein processing in ER could provide signi cant clues This article reports the results of research only. Mention of a proprietary product does not constitute an endorsement or a recommendation by the USDA for its use. The USDA is an equal opportunity provider and employer.

Competing interests
The authors declare that they have no competing interests.  A. Basic stem-loop structures of predicted microRNAs. miRDeep2 software was used to identify potential miRNA precursors based on nucleotide length, star sequence, stem-loop folding, and homology to the Varroa reference genome. It shows the predicted stem-loop structures, star, and mature sequences of predicted miRNA in Varroa mite. B. Annotation of predicted Varroa destructor microRNAs. 306 microRNAs were predicted in varroa samples. 50 of them were categorized as high con dence (n = ~50), low con dence (~80) miRNAs based on the standard criteria A. In silico differential expression of predicted microRNAs in female Varroa destructor relative to males.
EdgeR was used for differential expression analysis. 13 predicted microRNAs were down-regulated, 18 upregulated, while 60 were unaffected. B. Differential expression of individual microRNAs in Varroa females relative to males. miRNAs with a Log2 fold-change expression > |1| and FDR ≤ 0.1 were considered signi cantly differentially expressed with respect to male microRNAs (Refer Table1). Expression of these microRNAs is comparatively higher in female than male Varroa.
Page 27/29  Gene ontology (GO) derived biological processes related to genes targeted by upregulated/downregulated miRNAs in female V. destructor.