Integrative Analysis of lncRNA-mRNA Co-expression Provides Novel Insights Into the Regulation of Developmental Transitions in Female Varroa destructor

Introduction

Varroa destructor has been and still is the greatest biotic threat to the Western honeybee (Apis mellifera) worldwide since its migration from the Eastern honeybee (Apis cerana). As an obligate ectoparasitic mite, V. destructor feeds primarily on the fat body of brood and adult honeybees and thus results in decreased body weight, malformation, as well as a shorter life span of adult bees (Rosenkranz et al., 2010; Ramsey et al., 2019). V. destructor infestation is also intimately coupled to the pathogen transmission within and between honeybee colonies, which impairs immunity defenses at both the individual and colony levels. Great efforts have been made in developing solutions to efficiently control the mite so far; however, the combination of the mite with other biotic and abiotic stressors continues to shake the beekeeping and pollination industries (Noël et al., 2020).

Female V. destructor mites have two distinct stages in their life cycle, the dispersal phase on adult bees, and the reproductive phase, which takes place inside capped brood cells during bee development (Traynor et al., 2020). In the reproductive phase, the mother mite (foundress) completes reproduction by feeding on the fat body of the larva (Li et al., 2019). Depending on the chemical signal stimulus (kairomones), the foundress mites invade appropriately aged larval cells to oviposition, resulting in a single son (about 6.6 days) and up to 7–9 daughters (about 5.8 days) (Martin, 1994; Boecking and Genersch, 2008; Traynor et al., 2020). Once the eggs hatch, the juveniles develop from protonymphs to deutonymphs before molting into adults. Previous studies showed that the development of V. destructor is a complex process, in which many genes/proteins cooperate to regulate developmental transition and sexual differentiation (Cornman et al., 2010; McAfee et al., 2017; Mondet et al., 2018). It has also been shown that the development of V. destructor relies on effective coordination between physiological readiness and responsiveness to host cues (Evans and Cook, 2018). Although these advances have improved our understanding of this complex parasite, our knowledge of V. destructor biology is still insufficient, especially the regulation during the developmental transition from egg to adult remains to be explored.

Long non-coding RNAs (lncRNAs) are an important class of pervasive genes with more than 200 bases that function as regulators with less or no protein-coding capacity (Wang and Chang, 2011; Quinn and Chang, 2016). LncRNAs have been shown to play vital roles in regulating gene expression at transcription, RNA processing, and translation levels through cis (near the site of lncRNA production) or trans (co-expressed with their target gene) ways (Guttman et al., 2011; Batista and Chang, 2013; Li Y.R. et al., 2020). Growing studies suggest that lncRNAs participate in a broad range of biological processes (BPs). For example, in honeybees, lncRNAs have been engaged in the regulation of the oviposition and ovary size of queens (Humann et al., 2013; Chen et al., 2017; Chen and Shi, 2020), neural function (Sawata et al., 2004; Kiya et al., 2012), DNA methylation (Wedd et al., 2016), pathogen infection (Guo et al., 2018a,b), pesticide metabolism (Jayakodi et al., 2015; Chen et al., 2019; Fent et al., 2020), and behavioral transition (Liu et al., 2019). For V. destructor, however, the only report provided the first catalog of lncRNA profile in gravid female mites, suggesting that lncRNAs play a critical role in cellular and biological progresses (Lin et al., 2020). Nevertheless, so far, the potential biological function of lncRNAs during the developmental transitions of female V. destructor and their regulatory roles are still poorly understood.

To address these gaps in knowledge, in this study, we analyzed the expression profiles of lncRNAs and mRNAs across all the developmental stages of female V. destructor (i.e., egg, protonymph, deutonymph, and adult). Using RNA sequencing (RNA-Seq), we identified a large number of new lncRNAs and mRNAs, expanding the genetic database of V. destructor to a new depth. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the differentially expressed target genes of lncRNAs and mRNAs revealed critical pathways that play key roles in the developmental regulation of V. destructor at different stages. Our findings thus indicate significant roles of lncRNAs in orchestrating the development of female V. destructor mites and are potentially useful to develop effective methods for controlling this extremely threatening parasite.

Materials and Methods

Samples Collection and Total RNA Extraction

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.

Complementary DNA Library Construction and Sequencing

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.

Quality Analysis, Mapping, and Transcriptome Assembly

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).

Coding Capacity Analysis

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.

LncRNA and mRNA Expression Analyses

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.

Functional Analyses

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.

Quantitative Real-Time PCR Validation

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.

Results

Overview of RNA-Sequencing

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).

Identification and Characteristics of mRNAs and Long Non-coding RNAs

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).

Figure 1. Identification and characterizationof mRNAs and long non-coding RNAs (lncRNAs). The distribution of transcript number (A), exon number (B), and transcript length (C) of lncRNAs and mRNAs.

Dynamic Expression Patterns of mRNAs and Long Non-coding RNAs

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).

Figure 2. Expression profiles of mRNAs and long non-coding RNAs (lncRNAs). (A) The expression level of mRNAs and lncRNAs in four developmental stages. (B) Principal component analysis (PCA) based on mRNA. (C) PCA based on lncRNA. (D) Heatmap of mRNAs in the expression pattern profiles. (E) Heatmap of lncRNAs in the expression pattern profiles. Blue color indicates low expression, and red color denotes high expression.

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.

Figure 3. Time-series model of long non-coding RNAs (lncRNAs) and mRNAs.Time-series modules of lncRNAs (A) and mRNAs (B). Numbers in the top left corner indicate module number. Numbers in the lower left corner indicate the numbers of mRNAs or lncRNAs in each module. The Gene Ontology (GO) analysis of transcripts within clusters after the Fragments Per Kilobase of exon per Million fragments mapped (FPKM) filtering identified the associated enriched GO terms with corresponding enrichment P-value shown on right.

Differential Expression Analysis of mRNAs and Long Non-coding RNAs

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).

Figure 4. Differentially expressed long non-coding RNAs (lncRNAs) and functional analysis in differently compared libraries. (A) Volcano plot of the lncRNAs between protonymph and egg. (B) Volcano plot of the lncRNAs between deutonymph and protonymph. (C) Volcano plot of the lncRNAs between adult and deutonymph. The left blue points represent significantly decreased lncRNAs; gray points represent lncRNAs without significant changes; the right red points represent significantly increased lncRNAs. (D) The Gene Ontology (GO) analysis of the upregulated differentially expressed lncRNAs. (E) The GO analysis of the downregulated differentially expressed lncRNAs.

Gene Ontology Functional Analysis of mRNAs and Long Non-coding RNAs

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).

Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of mRNAs and Long Non-coding RNAs

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).

Figure 5. The enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the differentially expressed lncRNA. (A) Cis-target genes in the protonymph vs. egg comparison group. (B) Cis-target genes in the deutonymph vs. protonymph comparison group. (C) Cis-target genes in the adult vs. deutonymph comparison group.

Validation of mRNA and lncRNA

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.

Figure 6. Validation of three significantly differentially expressed long non-coding RNAs (lncRNAs) and their target mRNAs by quantitative PCR (qPCR). X-axis represents the four developmental stages. Y-axis indicates the relative expression of each gene; red lines are log₂(FPKM + 1) values of RNA-Seq, and the blue lines show relative expression by qPCR.

Discussion

The parasite V. destructor is an obligate ectoparasite and a highly invasive pest to honeybees (Nazzi and Le Conte, 2016). Due to the limited history of host-parasite coevolution, A. mellifera is suffering from V. destructor infestation. The specific parasitic life and increasing resistance to acaricides make developing efficient treatment strategies for V. destructor even more challenging. To gain a better understanding of the developmental regulation of female V. destructor, in this study, we performed an integrative analysis of lncRNA-mRNA co-expression crossing all the developmental stages (i.e., egg, protonymph, deutonymph, and adult). Consequently, the quantification analysis found 1,085, 2,461, and 5,429 mRNAs, and 862, 1,963, and 1,826 lncRNAs were differentially expressed in pairwise comparison (protonymph vs. egg, deutonymph vs. protonymph, and adult vs. deutonymph). Further functional enrichment revealed that fat digestion and absorption, ABC transporters, mitogen-activated protein kinase (MAPK) signaling pathway, and ubiquitin-proteasome pathway play key roles in the developmental transitions of female V. destructor.

Females of V. destructor show a development-specific pattern of lncRNA and mRNA expression. LncRNAs are an essential transcriptome component involved in the gene expression regulation in nearly all organisms studied (Fatica and Bozzoni, 2014; Gardini and Shiekhattar, 2015; Jarroux et al., 2017). In accordance with the lncRNA screening results of some extensively studied species, such as A. mellifera, Panonychus citri, and Tetranychus cinnabarinus (Zhang et al., 2014; Kim et al., 2018; Chen et al., 2019; Feng et al., 2020; Li Y. et al., 2021), our data showed relatively short length, low exon numbers, low transcript numbers, and generally low expression levels of V. destructor lncRNAs, implying that these features of lncRNAs were conserved in V. destructor. Accumulating evidence suggests that the expression profiles of lncRNAs exhibit tissue-, stage-, oocyte-, and/or development-specific patterns across species (Necsulea et al., 2014; Schor et al., 2018; Wang et al., 2020). In this study, we found that the number of identified mRNAs and lncRNAs gradually decreased from egg stage to deutonymph stage and then increased in the adult stage; however, the identified egg-specific and adult-specific lncRNAs and mRNAs were more than the other two stages, thus suggesting that gene expressions are more actively regulated in egg and adult stages. Furthermore, our quantitative results showed that, in the transition from protonymph to adult, the number of upregulated lncRNAs and mRNAs were gradually increased, indicating that more genes were regulated during developmental maturation.

Specific functional pathways are consolidated by lncRNA regulation at different developmental stages of female V. destructor. LncRNAs play regulatory roles in a wide spectrum of BPs, such as cell proliferation, cell differentiation and development, cell cycle, metabolism, apoptosis, and immune response (Wang et al., 2011; Batista and Chang, 2013; Geisler and Coller, 2013; Delas and Hannon, 2017). Particularly, in insects, lncRNAs have been implicated in mediating developmental processes and behavioral changes (Liu et al., 2019; Li G. et al., 2021; Yang et al., 2021). Our KEGG analysis of upregulated lncRNAs found that the fat digestion and absorption pathway was enriched in all the developmental stages of female V. destructor, which corresponds to the fact that lipid metabolism is fundamentally important to meet the energetic demands in growth and reproduction in insects (Arrese and Soulages, 2010). However, our results showed that energy was generated and utilized diversely at different developmental stages of V. destructor. During embryogenesis, energy was released by lipid mobilization to invest in the formation and development of organs; while, nymph and adult feed directly on honeybee fat bodies to obtain energy for development and reproduction (Zalewski et al., 2016; Ramsey et al., 2019). Noteworthily, in the fat digestion and absorption pathway, more upregulated lncRNAs were identified in deutonymph mites, echoing the notion that the deutonymphs have large energetic requirements to support energy-intensive developmental processes, such as metamorphosis (McAfee et al., 2017). In concert with this, ABC transporters and the Hippo signaling pathway were also found mainly elevated in deutonymphs. ABC transporters have been associated with the mite molting process (Li G. et al., 2021), with ABCGs particularly involved in pesticide transport and resistance (Dermauw and Van Leeuwen, 2014; Gott et al., 2017; Li Z. et al., 2020). We found that lncRNA LTCONS_00014321 and its target genes ABCG5 and ABCG8 were significantly upregulated in deutonymphs, implying that detoxification, metabolism, and especially metamorphosis are prominently reinforced in the deutonymph stage. Similarly, the upregulated Hippo signaling pathway is supposed to boost cell proliferation, tissue growth and differentiation, and organ migration in the deutonymph stage (Meng et al., 2016). In contrast, MAPK signaling was found consolidated in protonymphs. Given the functions of MAPK signaling in governing innate immune responses of mite (Tian et al., 2020), the upregulated lncRNAs and their target genes in this pathway might be correlated with the initial responses of V. destructor to environmental stressors such as pesticides (Lu et al., 2017). V. destructor is a well-known vector for the spreading of different honeybee viruses (Evans and Cook, 2018; Traynor et al., 2020). The ubiquitin-proteasome pathway of the host is usually employed or manipulated by viruses for their own needs (Gao and Luo, 2006); therefore, our results that the upregulated lncRNAs were mainly enriched in the adult stage may reflect the active promotion of virus-vector interactions (Gupta et al., 2019).

More evidence suggests that the deutonymph stage is crucial during the development of female V. destructor. Cuticle development is a fundamental physiological process in insect molting, keeping pace with growth, and metamorphosis (Liu et al., 2020; Li G. et al., 2021). During the maturation of V. destructor, nymph mites transit from possessing a soft cuticle to forming a tough exoskeleton in adults. In our results, both lncRNAs and mRNAs of some cuticle proteins (cuticle protein 65, cuticle protein 6.4, cuticle protein 63, and cuticle protein 38) were upregulated more in deutonymphs, suggesting that the deutonymph stage is critical in cuticle transition, which is in agreement with the findings in spider mites (Li G. et al., 2021). Another important upregulated gene in deutonymphs was mucin, which is primarily expressed in the peritrophic membrane in the midgut, plays important roles in protecting insects from infection by bacteria, viruses, and other pathogens (Hegedus et al., 2009), and has also been connected to keep cuticle integrity during nymph-adult molting (Zhao et al., 2020). We found that lncRNA LTCONS_00011419 and its target gene mucin 5AC (MUC5AC) were significantly highly expressed in the deutonymph stage than that in the adult stage, emphasizing the significant roles of the deutonymph stage in metamorphosis and pathogen resistance during the development of V. destructor.

Altogether, our present study provides the first integrated analysis of lncRNA-mRNA functional synergistic regulation in the developmental transition of female V. destructor. By depicting lncRNA-mRNA expression profiles and identifying differentially expressed lncRNAs and mRNAs crossing developmental stages, we contributed to gaining a foundational understanding of V. destructor biology. Further evidence from both pathway enrichment and analysis of critical genes indicates that the deutonymph is a key stage in the development of female V. destructor. We hope our results will help in developing more effective treatment techniques for this devastating pest by targeting crucial genes, such as cuticle proteins, in the deutonymph stage.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm.nih.gov/ PRJNA777703.

Author Contributions

BH and S-FX conceived and designed the experimental plan. J-LW, R-YH, N-NL, and JT collected the samples. J-LW, C-XZ, and R-YH contributed to reagents, materials, and analysis tools. J-LW, BH, and S-FX wrote the manuscript. All authors approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China (Grant No. 31902221), the Modern Agro-Industry Technology Research System (Grant No. CARS-44-KXJ-6), and the Agricultural Science and Technology Innovation Program (Grant No. CAASASTIP-2021-IAR).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.842704/full#supplementary-material

Supplementary Figure 1 | Workflow of this study.

Supplementary Figure 2 | The distribution of mRNAs (A) and long non-coding RNAs (lncRNAs) (B) in four developmental stages.

Supplementary Figure 3 | Differentially expressed mRNAs and functional analysis in differently compared libraries. (A) Volcano plot of the mRNAs between protonymph and egg. (B) Volcano plot of the mRNAs between deutonymph and protonymph. (C) Volcano plot of the mRNAs between adult and deutonymph. The left blue points represent significantly decreased mRNAs; gray points represent mRNAs without significant changes; the right red points represent significantly increased mRNAs. (D) The Gene Ontology (GO) terms analysis of the upregulated differentially expressed mRNAs. (E) The GO terms analysis of the downregulated differentially expressed mRNAs.

Supplementary Figure 4 | The enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the differentially expressed mRNA. (A) The KEGG analysis of mRNAs in the protonymph vs. egg comparison group. (B) The KEGG analysis of mRNAs in the deutonymph vs. protonymph comparison group. (C) The KEGG analysis of mRNAs in the adult vs. deutonymph comparison group.

Footnotes

1. ^ https://www.ncbi.nlm.nih.gov/assembly/GCF_002443255.1
2. ^ http://www.geneontology.org/
3. ^ http://www.genome.jp/kegg/

References