A total of three A. mellifera honeybee hives with mite populations were used in this study. All colonies were bred at the apiary of the Institute of Apicultural Research, Chinese Academy of Agricultural Sciences, Beijing, under the same conditions. No acaricide treatment was applied in this study to avoid pesticide bias. According to the previous sample collection method (Dietemann et al., 2013; McAfee et al., 2017), V. destructor families were collected from a single A. mellifera colony. Eggs, protonymph, deutonymph, and adults were transferred directly to tubes using a soft paintbrush and soft tweezers. Approximately 200 individuals were pooled for each replicate (four developmental stages, n = 3 for each stage). All collected samples were immediately fresh-frozen using liquid nitrogen and stored at −80°C for further processing.

Total RNA of each sample was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) according to the recommendations of the manufacturer. Extracted RNA was treated with Turbo DNA-Free Kit (Ambion Inc., Foster City, CA, United States) to remove contaminating genomic DNA. The quality and quantity of total RNA were evaluated using the A260:A280 ratio. NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States) was used to determine the concentration of RNA, and the integrity was analyzed with the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States) using the Agilent RNA 6000 Nano Kit.

First, before library construction, ribosomal RNA (rRNA) was removed using the Ribo-Zero™ rRNA Removal Kit (Illumina Inc., San Diego, CA, United States). Second, the generation of libraries for sequencing was made using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, United States) by following the instructions of the manufacturer. The products were then purified using a TruSeq RNASample Prep Kit v2 (NEB, Ipswich, MA, United States), and library concentration was measured using the Quant-iT PicoGreen Double-Stranded DNA Assay Kit (Invitrogen, Carlsbad, CA, United States). Fragment size was estimated using the Agilent High Sensitivity DNA Analysis Kit with the Agilent Bioanalyzer 2100 system and quantified by the quantitative real-time PCR (qRT-PCR) using the KAPA Library Quantification Illumina/ABI Prism Kit protocol (KAPA Biosystems, Wilmington, MA, United States). Finally, the libraries were sequenced on an Illumina HiSeq Xten platform. The sequencing data have been submitted to the NCBI Sequence Read Archive (SRA) and are accessible through the BioProject accession number: PRJNA777703.

The raw reads were filtered using the SOAPnuke (v1.5.2, parameters: −l 15 −q 0.2 −n 0.05 −i), and clean reads were obtained by removing contaminating adapter molecules, reads containing poly-N, and low-quality reads. All following analyses were performed based on the high-quality clean data. Sequenced reads in the FASTQ format were aligned to the V. destructor genome (Vdes_3.0¹) by using HISAT2 (v2.0.4, parameters: –phred64 –sensitive –no-discordant –no-mixed −l 1 −X 1000), and then the mapped reads were assembled using StringTie (v1.0.4, parameters: −f 0.3 −j 3 −c 5 −g 100 −s 10000 −p 8).

To identify new transcripts, all the reconstructed transcripts of each sample were aligned to the reference genome by using Cuffcompare; then, the novel transcripts were filtered and assembled to obtain putative transcripts. The putative novel lncRNAs were screened by the following steps: first, transcripts with transcript length >200 bp were selected; second, transcripts with Fragments Per Kilobase of exon per Million fragments mapped (FPKM)≥0.5 expression level were selected; third, transcripts with coverage >1 were selected; fourth, transcripts with one of the five class codes, namely, “i,” “j,” “u,” “o,” and “x,” were included (Sun et al., 2012), which were potentially recognized as novel ones. In addition, to further distinguish lncRNAs and mRNAs, the coding capacity of the novel transcripts was jointly calculated by Coding Potential Calculator (CPC), Coding-Non-Coding Index (CNCI), txCdsPredict, and Pfam-scan. Transcripts revealing a coding potential with CPC score < 0, CNCI score < 0, txCdsPredict score < 500, and Pfam-scan < 0.001 were identified as lncRNAs; others were considered as mRNAs. The transcript was identified as either lncRNA or mRNA only when the results were consistent among at least three judgment methods.

Gene expression was calculated using the normalized expression values as FPKM by RSEM (v1.2.12). Then, the differential expression analysis among the four groups was performed using the DEGseq software. Adjusted P-value ≤ 0.001 and absolute value of fold change ≥2 were set as the threshold for significantly differential expression of the coding genes and lncRNAs.

To predict the cis-target genes of identified lncRNAs, the adjacent protein-coding gene was searched within the range of 10-kb upstream or 20-kb downstream of a lncRNA (Yang et al., 2021). Then, the expressions of the lncRNAs and their candidate mRNAs were calculated by Spearman’s and Pearson’s correlation coefficients, and spearman_cor ≥ 0.6 and pearson_cor ≥ 0.6 were considered as the correlated expressions.

The differentially expressed mRNAs or lncRNA target genes were used for functional enrichment by the GO analysis and the KEGG pathway analysis. For the GO analysis, the differential expression genes (DEGs) were classified according to the principle classification of GO,² and the number of differential transcripts included in each GO category was counted. For the KEGG³ enrichment analysis, the hypergeometric distribution test was used, and pathways with uncorrected P-value ≤ 0.05 were considered as significantly enriched.

Three lncRNAs and their target mRNAs were selected for quantitative real-time PCR (qRT-PCR) to validate the RNA-Seq results. qRT-PCR was performed using SYBR Fast qPCR Kit Master Mix (2×) (KAPA Biosystems, Boston, MA, United States) on a QuantStudio 6&7 Flex Real-Time PCR system (Applied Biosystems, Foster City, CA, United States) according to the instructions of the manufacturer. All primers used in this study are listed in Supplementary Table 1. GAPDH was used as the reference gene. The conditions were set as follows: 95°C for 2 min, followed by 40 cycles (95°C for 3 s and 60°C for 30 s); then, 95°C for 15 s and 60°C for 15s; finally, 95°C for 15 s. The relative expression levels were then calculated based on the 2^–ΔΔCT method.

In this study, a total of 12 complementary DNA (cDNA) libraries were constructed and sequenced on the Illumine HiSeq Xten platform by using total RNA from each sample, and each sample was measured in three biological replicates (Supplementary Figure 1). After quality control, 12.64 Gb of clean reads were obtained, and the Q30 ranged from 83.54 to 97.66% for each sample. In total, 37.78–70.31% of clean reads from each library were mapped onto the V. destructor reference genome, and a total of 42,555 transcriptions were detected (Supplementary Table 2).

After mapping and coding capacity analysis, 9,760 (3,075 known and 6,685 novel) non-coding transcripts and 32,795 (28,091 known and 4,704 novel) protein-coding transcripts were obtained, which were used for downstream analyses. To comprehensively examine the differences between the two transcript species, we compared different characteristics (i.e., transcript number, transcript length, and exon number) of mRNAs and lncRNAs. The average length of mRNA was longer than that of lncRNA. Compared to mRNAs, a smaller number of exons and genes were identified in lncRNAs (Figure 1).

By calculating the gene expression levels from the RNA-Seq data using RSEM, 12,804 mRNAs and 7,954 lncRNAs were quantified in at least one biological replicate (FPKM > 0), and all lncRNAs tended to be expressed at a lower level than that of mRNAs (Figure 2A and Supplementary Tables 3, 4). For each stage, 11,750 mRNAs and 5,880 lncRNAs were identified in eggs, 11,502 mRNAs and 5,605 lncRNAs in protonymphs, 10,612 mRNAs and 3,572 lncRNAs in deutonymphs, and 11,241 mRNAs and 4,409 lncRNAs in adults (Supplementary Figure 2). Based on the expression changes, principal component analysis (PCA) was used to evaluate the distance relationship among the 12 samples. PCA score plot showed a clear separation of the four developmental stages (Figures 2B,C). Then, to understand the expression patterns of mRNA and lncRNA in mite samples, we generated the heatmaps using hierarchical clustering analyses to show the different changes of mRNAs and lncRNAs according to the development of female V. destructor(Figures 2D,E).

To further investigate the expression pattern of mRNAs and lncRNAs during the four developmental stages of female V. destructor, the Short Time-series Expression Miner (STEM) was applied to cluster mRNAs and lncRNAs with similar expression patterns, respectively. A total of 12 models were established for mRNAs and lncRNAs, but only five and three were significantly enriched (Bonferroni-adjusted P-value ≤ 0.05), respectively (Figures 3A,B). Given that the mRNAs with similar expression patterns may also be functionally connected, we subsequently performed functional enrichment for each mRNA model and found intracellular membrane-bounded organelle ATP-binding cassette (ABC) transporters (GO: 0043231), regulation of protein serine/threonine kinase activity (GO: 0071900), the cytosolic ribosome (GO: 0022626), the primary metabolic process (GO: 0044238), and the ribonuclease T2 activity (GO: 0033897) as the leading terms in each model.

To explore the key mRNAs and lncRNAs involved in the developmental regulation of female V. destructor, the data were pairwise compared within three groups (i.e., protonymph vs. egg, deutonymph vs. protonymph, and adult vs. deutonymph) (Supplementary Tables 5, 6). Consequently, in the protonymph vs. egg comparison, 1,085 differentially expressed mRNAs were detected, of which 498 were upregulated and 587 were downregulated (Supplementary Figure 3A). A total of 862 differentially expressed lncRNAs were detected, of which 454 were upregulated and 408 were downregulated (Figure 4A). In the deutonymph vs. protonymph comparison, 2,461 differentially expressed mRNAs were detected, of which 664 were upregulated and 1,797 were downregulated (Supplementary Figure 3B). A total of 1,963 differentially expressed lncRNAs were detected, of which 170 were upregulated and 1,793 were downregulated (Figure 4B). In the adult vs. deutonymph comparison, 5,429 differentially expressed mRNAs were detected, of which 3,992 were upregulated and 1,437 were downregulated (Supplementary Figure 3C). A total of 1,826 differentially expressed lncRNAs were detected, of which 1,112 were upregulated and 714 were downregulated (Figure 4C).

To address how lncRNAs interact with their target genes (mRNAs) to regulate the developmental process of V. destructor, the potential targets of lncRNAs by cis-regulation were predicted. A total of 4,089 target genes within the range of 10-kb upstream or 20-kb downstream were predicted in the cis role (Supplementary Table 7). To better understand the gene regulation involved in the developmental transitions of female V. destructor, especially at the BP, the molecular function (MF), and the cellular component (CC) levels, we clustered the GO annotation of the differentially expressed mRNAs (Supplementary Figures 3D,E) and lncRNAs (Figures 4D,E). Based on the BP in the GO classification, the differentially expressed mRNAs and lncRNAs were classified into different functional categories, and the top five GO terms were enriched in “cellular process,” “metabolic process,” “single-organism process,” “biological regulation,” and “regulation of biological process.” In the CC domain, the top five GO terms were “cell,” “cell part,” “membrane,” “membrane part,” and “organelle.” In the MF, “catalytic activity,” “binding,” “transporter activity,” and “signal transducer activity” were the most highly represented terms (Supplementary Tables 8, 9).

To gain further understanding of the biological functions of the differentially expressed lncRNAs and mRNAs during developmental transitions, the KEGG enrichment analyses of the mRNAs and lncRNA target genes in the three comparison groups were performed (Supplementary Tables 10, 11). In this study, we only selected the significantly enriched terms in the KEGG pathway analysis for each comparison, and the associations with human disease pathways were excluded. For protonymph vs. egg, deutonymph vs. protonymph, and adult vs. deutonymph comparisons, 16, 31, and 25 pathways that were highly enriched with mRNAs were detected, respectively. “ABC transporters” and “Hippo signaling pathway” were significantly enriched in the protonymph vs. egg comparison, “Glycosylphosphatidylinositol (GPI)-anchor biosynthesis” and “Thiamine metabolism” were significantly enriched in the deutonymph vs. protonymph comparison, and “Ribosome biogenesis in eukaryotes” and “RNA polymerase” were significantly enriched in the adult vs. deutonymph comparison (Supplementary Figure 4).

Notably, 5, 37, and 15 pathways were significantly enriched in the three comparison groups (i.e., protonymph vs. egg, deutonymph vs. protonymph, and adult vs. deutonymph), respectively. For example, “Fat digestion and absorption” and “ABC transporters” were significantly enriched in the protonymph vs. egg comparison; “Valine, leucine, and isoleucine degradation” and “D-Glutamine and D-glutamate metabolism” were significantly enriched in the deutonymph vs. protonymph comparison; and “Fat digestion and absorption” and “Riboflavin metabolism” were significantly enriched in the adult vs. deutonymph comparison (Figure 5).

To verify the RNA-Seq results, three predicted lncRNAs (i.e., XR_002672204.1, LTCONS_00014321, and LTCONS_00011419) and their target genes (i.e., cuticle protein 65, ABCG5, and MUC5AC) were selected, and their expression levels were detected by qRT-PCR. As shown in Figure 6, the expression of lncRNA XR_002672204.1-cuticle protein 65, LTCONS_00014321-ABCG5, and LTCONS_00011419-MUC5AC was high in the deutonymph stage compared to the other stages. Moreover, all three lncRNAs and their target genes showed similar expression patterns.