1. Introduction
Long non-coding RNAs (lncRNAs) are a class of linear RNA molecules transcribed by RNA polymerase II with a length of more than 200 nt [1]. Accumulating studies are indicative of the participation of lncRNAs in an array of physiological and pathological processes of great importance, such as cell proliferation and differentiation, material and energy metabolism, and immune defense, as well as the occurrence and progression of disease [2,3,4,5]. LncRNAs were earlier thought to be the byproducts of the transcription process without protein-coding ability; however, with the development of related knowledge and techniques, lncRNAs were suggested to exert pivotal biological functions in diverse manners, including the direct regulation of neighboring genes’ transcription, indirect modulation of downstream gene expression via the competing endogenous RNA (ceRNA) network, and translation into proteins or small peptides [6]. LncRNAs have been studied in depth in a few model organisms, such as humans and Arabidopsis thaliana; however, progress associated with lncRNAs in insects, including honey bees, is limited. Wu et al. [7] conducted a transcriptomic investigation of lncRNAs in Bombyx mori and found that a subseries of lncRNAs were potentially engaged in the regulation of biosynthesis, transport, and secretion of silk proteins by serving as ceRNAs or miRNA precursors. On the basis of a transcriptome-wide analysis of lncRNAs, miRNAs, and mRNAs in Apis mellifera queens during different oviposition periods, Chen et al. [8] revealed that these ncRNAs may be involved in ovary activation and oviposition processes through the modulation of related genes’ expression. Previously, following the next-generation sequencing of the Apis mellifera ligustica workers’ midguts tissues, we identified a total of 6353 lncRNAs and observed that their structural features were analogous to those discovered in other animals and plants [9]. An increasing number of studies have demonstrated that lncRNAs are critical regulators in insect–pathogen interactions; for example Jayakodi et al. [10] who identified 2470 and 1514 long intergenic noncoding RNAs (lincRNAs) in Apis cerana and A. mellifera based on deep sequencing and bioinformatics, detected that lincRNAs exhibited a tissue-specific expression pattern among seven organs or tissues, and they further discovered that 10 host lincRNAs were induced to activation upon SBV and DWV infection. Our group previously analyzed the expression profile of lncRNAs in the A. m. ligustica workers’ midguts challenged by Nosema ceranae, and detected that differentially expressed lncRNAs (DElncRNAs) were likely to participate in modulating the host response via various methods, such as a cis-acting effect, trans-acting effect, miRNA precursors, and ceRNA networks [9]. In addition, we deciphered the expression pattern of lnRNAs in Apis cerana cerana 6-day-old larval gut responding to A. apis invasion, and investigated the putative regulatory role of DElncRNAs in the host response [11].
A. cerana is widely reared in China and many other Asian countries, and plays a crucial role in pollination for considerable crops and wild flowers, as well as production of api-products [12]. A. cerana is the original host for N. ceranae, a widespread fungal parasite causing bee nosemosis, which results in severe losses for the apicultural industry [13]. A. c. cerana, the nominate subspecies of A. cerana, is a major bee species used in beekeeping practices in China. We previously investigated the dynamics of miRNAs in the A. c. cerana workers’ midguts infected by N. ceranae and revealed the potential role of DEmiRNAs in host immune responses [14]. Recently, we analyzed the differential expression pattern of circular RNAs (circRNAs) in the midguts of A. c. cerana workers during N. ceranae infection, and we uncovered that DEcircRNAs were likely to modulate host cellular and humoral immune responses through the regulation of the transcription of source genes or absorption of target miRNAs to affect the downstream gene expression [15]. As a vital layer of post-transcription regulation, lncRNAs are able to not only directly modulate gene expression, but also interact with other ncRNAs such as miRNAs to indirectly regulate the expression of target genes [16]. However, compared with A. mellifera, studies on lncRNAs in A. cerana are scarce at present. Whether and how lncRNAs participate in the interaction between A. cerana and N. ceranae has been unclear until now.
Our team previously conducted deep sequencing of cDNA libraries from the midgut tissues of A. c. cerana workers at 7 days post inoculation (dpi) and 10 dpi (AcT7 and AcT10 groups), as well as corresponding uninoculated workers’ midgut tissues (AcCK7 and AcCK10 groups) [9]. In this current work, on the basis of the obtained high-quality transcriptome datasets, the systematic identification and structural characterization of lncRNAs in the A. c. cerana workers’ midguts were conducted, and the differential expression profile of lncRNAs during the N. ceranae infestation was then analyzed in depth, followed by further investigation of the putative regulatory functions of DElncRNAs in host response, with a focus on immune response. Our findings can offer not only a foundation for identifying the mechanism underlying the DElncRNA-regulated response to N. ceranae infestation, but also a novel insight into the Asian honey bee microsporidian interaction.
2. Results
2.1. Quantity and Structural Properties of A. c. cerana lncRNAs
Here, 2365, 2322, 2487, and 1986 lncRNAs were, respectively, identified in the AcCK7, AcT7, AcCK10, and AcT10 groups; the numbers of known lncRNAs were 1001, 960, 1089, and 799, while those of novel lncRNAs were 1364, 1362, 1398, and 1187, respectively. Among the identified lncRNAs, the three most highly expressed were XR_001765028.1, XR_001765880.1, and XR_001766991.1. Additionally, it was found that 1558 lncRNAs were shared by the abovementioned four groups, whereas the quantities of unique ones were 131, 132, 177, and 70, respectively (Figure 1A). After removing redundant ones, a total of 3496 A. c. cerana lncRNAs (Table S1) were ultimately identified, with a length distribution range of 188 nt~441,800 nt (Figure 1B). As shown in Figure 1C, the length of introns in lncRNAs was smaller than that in mRNAs: 1254 (35.90%) lncRNAs contained only one intron, while the remaining 2242 (64.10%) included two or more introns (Figure 1D). Similarly, the length of exons in lncRNAs was shorter than that in mRNAs (Figure 1E); as many as 2162 (61.84%) lncRNAs contained two exons, and another 1334 (38.16%) lncRNAs included more than two exons (Figure 1F).
2.2. Differential Expression Profile of lncRNAs in the Midguts of A. c. cerana Workers Responding to N. ceranae Infestation
In the AcCK7 vs. AcT7 comparison group, 50 upregulated and 29 downregulated lncRNAs were identified (Figure 2A, see also Table S2); among these, XR_001767038.1 (log2Fold change = 12.47) (log2FC = 12.47), TCONS_00047233 (log2FC = 10.55), and XR_001766127.1 (log2FC = 10.44) were the three most upregulated lncRNAs, while TCONS_00006288 (log2FC = −10.76) was the most downregulated one, followed by XR_001765791.1 (log2FC = −10.59) and XR_001766607.1 (log2FC = −10.37) (Figure 2B). Comparatively, 26 upregulated and 47 downregulated lncRNAs were detected in the AcCK2 vs. AcT2 comparison group (Figure 2C, see also Table S3); among these, three lncRNAs including TCONS_00029732 (log2FC = 13.07), TCONS_00018981 (log2FC = 12.45), and XR_001767090.1 (log2FC = 11.93) were the most upregulated, whereas TCONS_00004222 (log2FC = −15.01) was the most downregulated one, followed by TCONS_00023581 (log2FC = −11.55) and XR_001766877.1 (log2FC = −11.46) (Figure 2D). Further analysis showed that there were six shared DElncRNAs (TCONS_00011721, TCONS_00022015, TCONS_00024562, XR_001765500.1, XR_001765515.1, and XR_001766127.1) by the aforementioned two comparison groups, while the numbers of unique ones were 72 and 67, respectively (Figure 2E). Detailed information regarding DElncRNAs in AcCK7 vs. AcT7 and AcCK7 vs. AcT10 comparison groups are, respectively, presented in Tables S2 and S3.
2.3. Analysis of the Cis-Acting Effect of DElncRNAs on Host Response to N. ceranae Infestation
Accumulating evidence showed that lncRNAs are capable of regulating the expression of one or several linked and continuous genes in close genomic proximity (called the cis-acting effect) [17]. It was suggested that 79 DElncRNAs in the AcCK7 vs. AcT7 in the comparison group potentially regulated 87 up- and downstream genes, which were engaged in 30 Gene Ontology (GO) terms relevant to biological process, cellular component, and molecular function, such as metabolic process, cell part, and transporter activity (Figure 3A, see also Table S4). In contrast, 73 DElncRNAs in the AcCK7 vs. AcT10 comparison group putatively modulated 87 up- and downstream genes, which were enriched in 13 functional terms such as regulation of biological process, muti-organism process, and binding (Figure 3B, see also Table S4).
In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis demonstrated that the up- and downstream genes in the workers’ midguts at 7 dpi were associated with 82 pathways, including apoptosis, adhesion junction, and phagosome (Figure 4A, see also Table S5). Comparatively, the up- and downstream genes in the workers’ midguts at 10 dpi were relative to 29 pathways, such as purine metabolism, RNA degradation, and endocytosis (Figure 4B, see also Table S5).
2.4. Investigation of the Trans-Acting Effect of DElncRNAs on the Host Response to N. ceranae Infestation
In addition, those lncRNAs diffusing from the synthesis site were able to act directly on many genes at a great distance at the same chromosome or other chromosomes (named trans-acting effect) [18]. It was detected that a total of 235 mRNAs were co-expressed with DElncRNAs in the AcCK7 vs. AcT7 comparison group (Table S6); among these, XR_001765926.1 was found to co-express with the largest number of mRNAs (90). For example, the co-expressed zinc transporter ZIP11 protein encoding gene (XM_017059863.1) had a positive correlation with XR_001766833.1, but had a negative correlation with TCONS_00021820; the co-expressed SEC23-like interacting protein isoform X2 encoding gene (XM_017063532.1) was positively correlated with XR_001766919.1, while it was negatively correlated with XR_001765791.1. Meanwhile, 209 mRNAs were observed to be co-expressed with DElncRNAs in the AcCK10 vs. AcT10 comparison group (Table S6); among these, TCONS_00029732, XR_001765691.1, and XR_001766718.1 were observed to co-express with the largest number of mRNAs (88) (Table S6). For instance, the co-expressed DNA mismatch repair protein MSH2 encoding gene (XM_017051143.1) had a positive correlation with TCONS_00007982 and XR_001766995.1; the co-expressed alkaline phosphatase 4 protein-encoding gene (XM_017048382.1) was negatively correlated with TCONS_00047837, TCONS_00004245, and TCONS_00044926 (Table S6); and the co-expressed promoter complex subunit 7 (XM_017063362.1) was negatively correlated with XR_001765313.1 and TCONS_00004222 (Table S6).
2.5. Analysis of DElncRNA-Involved ceRNA Regulatory Networks
Various RNA molecules, including lncRNAs, can competitively bind to miRNAs to modulate the expression of downstream genes and further regulate biological processes such as host response and pathogen infection (called ceRNA mechanism) [19]. Here, nine DElncRNAs in the workers’ midguts at 7 dpi were observed to potentially target eight DEmiRNAs, and further bound to twenty-five DEmRNAs, forming five independent subnetworks (Figure 5A); additionally, TCONS_00022015 targeted both miR-1-x and miR-598-y, TCONS_00022015 targeted both miR-252-y and miR-6717-x, while the other seven DElincRNAs had only one target (Figure 5A). Comparatively, six DElncRNAs in the workers’ midguts at 10 dpi were found to target five DEmiRNAs and further link to eight DEmRNAs, forming five independent subnetworks (Figure 5B); additionally, TCONS_00007982 and XR_001765691.1 both targeted novel_m0003-3p, while the other four DElncRNAs had only one target. (Figure 5B) Detailed information on the targeting relationships among DElncRNAs, DEmiRNAs, and DEmRNAs in the AcCK7 vs. AcT7 and AcCK10 vs. AcT10 comparison groups are shown in Table S7.
2.6. Analysis of DElncRNAs as Putative miRNA Precursors
As a kind of versatile regulator, certain lncRNAs were verified to function as precursor molecules that are processed into miRNAs, greatly enlarging the reservoir of miRNAs within cells [20]. In the workers’ midguts at 7 dpi, TCONS_00024312 and XR_001765805.1 were predicted to be precursors for ame-miR-315 and ame-miR-927, respectively (Figure 6A). As shown in Figure 6B, TCONS_00006120 in the workers’ midguts at 10 dpi was found to be a putative miRNA precursor targeting the mature sequences of both ame-miR-87-1 and ame-miR-87-2.
2.7. Validation of DElncRNAs by RT-qPCR
The primers for RT-qPCR assay were subjected to PCR amplification and agarose gel electrophoresis; the results showed that the expected fragments were amplified using each primer pair (Figure S2), verifying the specificity of these primers. The RT-qPCR results indicated that the expression trends of the randomly selected 12 DElncRNAs were consistent with those in the sequencing data (Figure 7), indicative of the authenticity and reliability of the transcriptome datasets used in this work.
3. Discussion
3.1. Quantity, Structural Property, and Expression Pattern of lncRNAs in the Midguts of A. c. cerana Workers Infected by N. ceranae
Here, a total of 3031 lncRNAs with a length distribution ranging from 188 nt to 441,800 nt were for the first time identified in A. c. cerana midguts. In view of the limited quantity of lncRNAs in A. cerana, the identified lncRNAs could offer valuable resources for the study of lncRNAs in A. c. cerana and other subspecies belonging to A. cerana. Given that the expression of lncRNAs is the tissue, developmental, and stress stages, the identified lncRNAs were speculated to be only a fraction of the total lncRNAs in A. c. cerana. Additionally, the identified A. c. cerana lncRNAs were detected to have shorter lengths of exons and introns, analogous to the characteristics of lncRNAs discovered in other species [21,22].
An array of studies showed that pathogen infections affected the expression pattern of lncRNAs in animals and plants [23,24,25,26]. In this work, 79 and 73 DElncRNAs were screened in the midguts of A. c. cerana workers at 7 dpi and 10 dpi with N. ceranae, respectively, indicative of the overall change of host lncRNAs caused by N. ceranae infestation. Further analysis suggested that six DElncRNAs, such as TCONS_00011721, TCONS_00022015, TCONS_00024562, XR_001765500.1, XR_001765515.1, and XR_001766127.1, are shared by the abovementioned comparison groups, implying that these common DElncRNAs may be vital regulators in host response to N. ceranae infestation, thus deserving further investigation.
3.2. DElncRNAs Were Potentially Involved in Modulating Host Response to N. ceranae Infestation via Both cis- and Trans-Acting Effects
The action modes of lncRNAs are flexible and diverse [3]. One major regulatory manner of lncRNAs is the cis-acting effect, which has been verified in numerous animals, such as Plutella xylostella [27] and A. m. ligustica [9]. LncRNAs located upstream and downstream of the protein-encoding genes may generate partial overlap with promoters, enhancers, or inducible elements, thus exerting regulatory functions in the expression of neighboring protein-coding genes at the transcriptional or post-transcriptional level [2,3]. The Notch signaling pathway is highly conserved in various insects and plays critical roles in controlling organ development and tissue homeostasis [28,29]. The Mitogen-activated protein kinases (MAPK) signaling pathway acts as a signaling hub that can be activated by external signals, such as mitogens, hormones, temperature, and osmolality [30]. Additionally, a series of signaling pathways, including MAPK, JNK, ERK, P53, and PI3K/Akt, jointly regulate diverse insect biological processes, such as eye development, cell shape, cellular apoptosis, neural development, embryonic dorsal plate closure, and immune response [28,31]. The Notch and MAPK signaling pathways could synergistically modulate insect growth, development, and immunity [32,33]. Here, we observed that two neighboring genes of DElncRNAs in the host midguts at 7 dpi, such as the Nicastrin protein isoform X2 encoding gene (NCSTN) (XM_017058633.1) and γ-secretase protein subunit Aph-1 encoding gene (APH-1) (XM_017060355.1)), were upregulated and enriched in the Notch signaling pathway; the vein-like protein encoding gene Vein (XM_017059352.1) in the host midguts at 10 dpi was upregulated and enriched in the MAPK signaling pathway. These results indicate that the corresponding host DElncRNAs (TCONS_00032459, XR_001766455.1, and TCONS_00034390) were putatively engaged in regulating the Notch and MAPK signaling pathways in a cis-acting manner during N. ceranae infestation (Figure 8).
Insects lack adaptive immunity and rely on innate immunity in defense against pathogen invasion [32]. Insect innate immunity includes cellular immune response, such as phagocytosis, colocalization, encapsulation reaction, melanization, ubiquitin-mediated protein hydrolysis, and lysosome, as well as humoral immune responses, represented by the MAPK, Jak-STAT, and NF-κB signaling pathways [28,34,35,36,37,38,39,40]. Here, three, two, one, and one neighboring genes of DElncRNAs relevant to phagosome, endocytosis, autophagy, and lysosome, respectively, were detected in the workers’ midguts at 7 dpi, indicating that the corresponding DElncRNAs were potentially involved in the N. ceranae response of the host by regulating the aforementioned four vital cellular immune pathways with a cis-acting effect.
P53, a tumor suppressor protein that acts as a transcription factor, affected the release of apoptosis-forming factors, such as cytochrome c in mitochondria, by downregulating Bcl-2 expression and, meanwhile, upregulating Bax expression, which, in turn, controlled the cell cycle [41]. The p53 signaling pathway was involved in maintaining genomic stability and regulating DNA damage repair, as well as apoptosis [41]. Here, we observed the upregulation of the cytochrome c protein encoding gene (CytC) (XM_017063771.1), enriched in both apoptosis and the p53 signaling pathway; the universal transcription factor IIH subunit 3 encoding gene (TTDA) (XM_017063798.1), enriched in nucleotide excision repair, was also detected as upregulated. Together, these results demonstrate that corresponding DElncRNAs may be engaged in regulating apoptosis, the p53 signaling pathway, and nucleotide excision repair through activating the expression of CytC and TTDA via a cis-acting manner, further affecting the host response to N. ceranae infestation. However, additional work is required to dissect the underlying mechanisms of the DElncRNA-regulated cis-acting manner. LncRNAs are capable of modulating gene expression via interaction with their co-expressed mRNAs [42,43]. In the present study, co-expressed genes with DElncRNAs in the AcCK7 vs. AcT7 and AcCK10 vs. AcT10 comparison groups were associated with a suite of functional terms, such as stress response, catalytic activity, and biological regulation; six cellular and humoral immune pathways, such as ubiquitin-mediated protein hydrolysis, endocytosis, lysosome, autophagy, and MAPK signaling pathways; 11 signaling pathways associated with growth and development including Wnt, FoxO, Hippo, Hedgehog, mTOR, and p53; and 42 pathways relative to material and energy metabolisms, such as carbon metabolism, sulfur metabolism, and oxidative phosphorylation [9,28,34,35,36,37,38,39,40] (Figure S1). The results showed that DElncRNAs may affect many aspects of host midguts’ growth and development, metabolism, and immune defense in a trans-acting manner, further responding to N. ceranae infestation. More efforts are needed to further ascertain the mechanisms underlying DElncRNA-mediated trans-acting effects in the future.
3.3. DElncRNAs May Play a Role in Host Response to N. ceranae Infestation by Serving as miRNA Precursors
MiRNAs are key regulators in gene expression and considerable biological processes, such as immune response and apoptosis [44,45]. Increasing studies indicate that lncRNAs can act as miRNA precursors to generate abundant mature miRNAs through a series of splicing and catalysis [14]. By using the immunohistochemistry and luciferase reporter system, Yuan et al. [44] reported that miR-315 affected the structure and function of Drosophila neurosynapses by repressing the expression of the target gene dFMR1, followed by the knockdown of miR-315, which led to embryonic death, while the overexpression of miR-315 gave rise to pupation defects. Additionally, it is suggested that miR-315 is mainly expressed in the nervous system of Drosophila and acts as a specific activator of the wingless signaling pathway, which is involved in regulating wing development [45]. MiR-927 was observed to affect the expression of antimicrobial peptides associated with the Toll signaling pathway of Aedes albopictus through the inhibition of serine and to promote DENV infection in humans [46]. He et al. [47] detected that overexpression of miR-927 in Drosophila melanogaster liposomes resulted in the downregulation of the expression of the target gene krl-h1, which resulted in a significant decrease in egg production, increased embryonic lethality, and pupation defects. MiR-87, a key regulator of Drosophila neuronal dendrites C4da, was found to promote neuronal dendrite regeneration by targeting the transcriptional repressor Ttk69 [48]. Here, three DElncRNAs (TCONS_00024312, XR_001765805.1, and TCONS_00006120) were observed to be potential precursors for miR-315, miR-927, and miR-87. It is speculated that these three DElncRNAs may participate in the modulation of gut development and AMP synthesis by generating abundant miRNAs, further regulating host N. ceranae response.
3.4. DElncRNAs May Play a Part in Host Response to N. ceranae Infestation through ceRNA Regulatory Networks
MRE-containing lncRNAs can absorb miRNAs via ceRNA networks, thereby reducing the binding of miRNAs to target mRNAs and attenuating the repressive effect on target genes [49,50]. Liu et al. [51] discovered that ame-miR-279a was highly expressed in Kenyon cells in the mushroom body of the brains of A. mellifera nurse and forager bees, and the overexpression of ame-miR-279a resulted in an enhanced proboscis extension response to sucrose solution, whereas the knockdown of ame-miR-279a led to a diminished proboscis extension response. Here, two DElncRNAs (XR_001766050.1 and TCONS_00024562) in the host midguts at 7 dpi and three DElncRNAs (TCONS_00024562, XR_001765691.1 and TCONS_00042438) in the host midguts at 10 dpi were observed to target miR-279-y, suggesting that these five DElncRNAs, serving as “molecular sponges”, may be involved in the N. ceranae response of A. c. cerana workers.
In Drosophila, miR-315 was a vital regulator in synaptic structure and transmission by targeting dfmr1 [44]. Kitatani et al. documented that miR-87 promoted dendrite regrowth during regeneration, at least in part through suppressing Ttk69 in Drosophila sensory neurons, suggesting that developmental- and injury-induced dendrite regeneration share a common intrinsic mechanism to reactivate dendrite growth [48]. Here, ten and four DElncRNAs in the abovementioned two comparison groups were detected to target miR-315-x, and seven and four DElncRNAs putatively targeted miR-87-y. The results showed that corresponding DElncRNAs may be engaged in modulating signal transmission and midgut growth and regeneration of workers’ midguts during N. ceranae infestation.
The overexpression of miR-927 in the Drosophila larval fat body significantly decreased the expression of Kr-h1 and resulted in reduced oviposition, increased mortality, delayed pupation, and reduced pupal size [47]. In mosquito cell response to viral infection, both the overexpression and knockdown of miR-927 led to the alteration of the expression levels of antimicrobial peptides involved in the Toll signaling pathway [46]. Here, five and three DElncRNAs in the host midguts at 7 dpi potentially targeted miR-927-x and miR-927-y, respectively. This indicates that these DElncRNAs may be involved in the host immune response to N. ceranae infestation.
miR-1, a conserved miRNA family that are abundantly expressed in mammals and insects, are widely involved in growth, development, metabolism, and immunity [52,53]. Tci-miR-1-3p was proved to be closely correlated with the resistance of Tetranychuscinnabarinus to cyflumetofen [52]. In Bactrocera dorsalis, miRNA-1-3p is an important sex determinant in the early male embryo and suppresses the expression of Bdtra [53]. We detected that TCONS_00022015 in the host midguts at 7 dpi putatively targets miR-1-x, which further targets max-binding protein MNT-like transcript variant X3 encoding gene (XM_017051758.1) and uncharacterized protein encoding gene (XM_017062549.1); these targets were enriched in the Hedgehog pathway. We inferred that TCONS_00022015 may be employed by the host to absorb miR-1-x to regulate the Hedgehog signaling pathway and ubiquitin-mediated protein hydrolysis in the N. ceranae response of workers’ midguts.
In Drosophila, miR-252 exerted regulatory functions in metamorphosis, cell cycle, wing formation, and body weight [54]. A. albopictus miR-252 inhibited the replication of DENV by regulating the expression of envelope protein genes [55]. Verma et al. reported that the absence of miR-60 expression in the gut of Caenorhabditis elegans resulted in the upregulation of genes relative to lysosomal protease, exogenous metabolism, and immune defense, as well as the activation of endocytosis, thereby promoting adaptation to chronic oxidative stress [56]. miR-60 in nematode intestines was suggested to participate in controlling the response to heat and cold stress by regulating the expression of zip-10 [57]. Here, we observed that TCONS_00032996 potentially targeted miR-252-y, further targeting the telomerase-binding protein EST1A-like isoform X2 encoding gene (XM_017057307.1), a member closely associated with the mRNA surveillance pathway; TCONS_00032996 putatively targeted miR-252-y, further targeting the enriched telomerase-binding protein EST1A-like isoform X2 encoding gene (XM_017057307.1), enriched as well as the mRNA surveillance pathway. As structure-specific and multifunctional enzymes, flap endonucleases are widely observed in eukaryotes and are involved in DNA replication and repair, thus playing a critical part in maintaining genomic stability. Here, TCONS_00035424 in the host midguts at 10 dpi was found to target miR-60-y, which further targeted the gene encoding the flap endonuclease GEN isoform X1 (XM_017063363.1). In summary, the results demonstrated that TCONS_00032996, TCONS_00032996, and TCONS_00035424 in the workers’ midguts were likely to competitively bind to miR-252 and miR-60 to modulate the mRNA surveillance pathway and DNA replication and repair, further mediating the host response to the infestation by N. ceranae. Recently, we established the technical platform for the functional investigation of honey bee lncRNAs and miRNAs [58,59]. In the near future, we will conduct RNAi of these promising candidate DElnRNAs mentioned above to further explore their ceRNA functions in the response of A. c. cerana workers’ midguts to N. ceranae invasion.
4. Materials and Methods
4.1. Biological Materials
A. c. ceranae colonies were reared in the apiary of the College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou city, China. N. ceranae spores were previously purified [14,60] and conserved at the Honey Bee Protection Laboratory, Fujian Agriculture and Forestry University, Fuzhou city, China.
4.2. Transcriptome Data Source
On the basis of our previously established method, the midgut tissues of A. c. cerana workers at 7 dpi and 10 dpi with N. ceranae spores (AcT7 and AcT10 groups), as well as corresponding uninoculated workers’ midgut tissues (AcCK7 and AcCK10 groups), were, respectively, dissected and then subjected to RNA isolation, cDNA library construction, and Illumina sequencing (Guangzhou Gene Denovo Biotechnology Co., Ltd., Guangzhou, China), followed by strict quality control of the raw data [15], which were deposited in the NCBI SRA database (
4.3. Prediction and Characterization of lncRNAs
Following the method described by Chen et al. [9], a combination of CPC [61] and CNCI [62] software was used to predict lncRNAs in the AcCK7, AcCK10, AcT7, and AcT10 groups, and the intersection of the prediction results was regarded as reliable lncRNAs. The expression level of each lncRNA was normalized to FPKM (Fragments Per Kilobase Million). The OmicShare platform (
4.4. Analysis of DElncRNAs
The DElncRNAs in the AcCK7 vs. AcT7 and AcCK10 vs. AcT10 comparison groups were screened using edgeR software (version 4.2) [63] (
4.5. Investigation of the Cis-Acting and Trans-Acting Effects of DElncRNAs
According to the method described by Ye et al. [64], neighboring genes within 10 kb upstream and downstream of DElncRNAs were surveyed and considered as potential targets. Subsequently, using the Blast tool in the NCBI database (
The co-expressed genes were predicted by correlation analysis or co-expression analysis of lncRNAs and protein-coding genes in the samples. The Pearson correlation coefficient method was used to analyze the correlation between lncRNAs and protein-coding genes between samples, and only the maximum positive correlation result and the maximum negative correlation result were retained. Next, the interactions between DElncRNAs and mRNAs were visualized using Cytoscape (version 3.7.2) (Bethesda, MD, USA) [65]. Further, the co-expressed mRNAs were subjected to GO term and KEGG pathway analyses utilizing the Blast tool.
4.6. Target Prediction and ceRNA Network Analysis
Following the method described by Guo et al. [66], the DElncRNA-targeted DEmiRNAs and DEmiRNA-targeted DEmRNAs were predicted with TargetFinder software (Livermore, CA, USA,
4.7. Analysis of the DElncRNAs Serving as miRNA Precursors
Based on the sequence information and structural features of the miRNA precursors, the sequences of the DElncRNAs were matched to the miRBase database (
4.8. RT-qPCR Validation of DElncRNAs
To verify the reliability of the transcriptome data used in this study, four DElncRNAs in the AcCK7 vs. AcT7 comparison group and eight DElncRNAs in the AcCK10 vs. AcT10 comparison group were randomly selected for RT-qPCR validation. Specific primers (Table S8) for the selected DElncRNAs were designed with DNAMAN software (version 10) (LynnonBiosoft, San Ramon, CA, USA) and then synthesized by Sangon Biotech (Shanghai) Co., Ltd. (Shanghai, China). The total RNA of the midgut samples (n = 3) in the AcCK7, AcT7, AcCK10, and AcT10 groups were isolated using RNA extraction kits (Takara, Dalian, China). Next, reverse transcription was performed using Oligo (dT)23 primers to obtain the corresponding cDNAs (Vazymes, Nanjing, China), which were then used as templates for qPCR detection. The actin gene (GeneBank ID: 107999330) was used as the internal reference. The 20 μL reaction system contained 1 μL of reverse primer (10.0 μmol·L−1), 1 μL of cDNA, 10 μL of SYBR Green Dye (Vazymes, Nanjing, China), and 7 μL of sterile water. The reaction was conducted on a QuanStudio Real-Time PCR System (ThemoFisher, Waltham, MA, USA) following the cycling parameters of 95 °C for 1 min, followed by 40 cycles at 95 °C for 15 s, 60 °C for 30 s, and 72 °C for 45 s. The relative expression level of each DElncRNA was calculated using the 2–∆∆CT method [69]. GraphPad Prism 6.0 software was used for the relevant data analysis and plotting.
4.9. Statistical Analysis
All statistical analyses were performed using SPSS software (IBM, Armonk, NY, USA) and GraphPad Prism 6.0 software. The data are shown as the mean ± standard deviation (SD). The statistical analysis was conducted with the independent samples using t-test and one-way ANOVA.
5. Conclusions
In conclusion, 1795 known lncRNAs and 1701 novel ones were for the first time identified in A. c. cerana workers’ midguts. These A. c. cerana lncRNAs shared similar structural characteristics to those discovered in other animals and plants, such as shorter exons and introns than mRNAs. The overall expression profile of lncRNAs in host midguts was changed due to the fact of N. ceranae infestation. DElncRNAs may play regulatory roles in host N. ceranae response by the regulation of neighboring genes via a cis-acting effect, modulation of co-expressed mRNAs via a trans-acting effect, and control of downstream target gene expression via ceRNA networks.
Conceptualization, R.G. and D.C.; methodology, Z.W. and S.W.; software, Z.W., S.W., X.F. and K.Z.; validation, Z.W., S.W. and J.Z.; formal analysis, Z.W., S.W., H.Z. and X.G.; data curation, Z.W., S.W., Y.Z., S.G. and D.Z.; writing—original draft preparation, Z.W. and S.W.; visualization, Q.L. and Z.N.; supervision, R.G.; project administration, R.G. and D.C.; funding acquisition, R.G. and D.C. All authors have read and agreed to the published version of the manuscript.
Not applicable.
Not applicable.
Raw data generated from RNA-seq are available in the NCBI SRA database under the BioProject number: PRJNA562787 [
All authors thank the reviewers and editors for their constructive comments and recommendations. R.G. genuinely appreciates his beloved father and mother.
The authors declare no conflict of interest. The funders had no role in the design of the study; the collection, analyses, or the interpretation of data; the writing of the manuscript; or the decision to publish the results.
Footnotes
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Enlarge this image.Figure 1. Comparison of the structural characteristics between A. c. cerana lncRNAs and mRNAs. (A): Venn diagram of the lncRNAs identified in four groups; (B): length distribution of lncRNAs and mRNAs; (C): overall lengths of introns in lncRNAs and mRNAs; (D): number statistics of introns in lncRNAs and mRNAs; (E): overall lengths of exons in lncRNAs and mRNAs; (F): number statistics of exons in lncRNAs and mRNAs.
Enlarge this image.Figure 2. Differential expression profile of lncRNAs involved in response of A. c. cerana workers’ midguts to N. ceranae infestation. (A,B): Volcano diagrams of DElncRNAs in the AcCK7 vs. AcT7 and AcCK10 vs. AcT10 comparison groups, DElncRNAs are marked with bule (down), red (up) and grey (not significant), the criteria for a DElncRNAs are |log2FC| > 1 and p < 0.05; (C): Venn diagram of DElncRNAs; (D,E): expression clustering of up- and downregulated lncRNAs in the AcCK7 vs. AcT7 comparison group; (F,G): expression clustering of up- and down-regulated DElncRNAs in the AcCK10 vs. AcT10 comparison group.
Enlarge this image.Figure 3. GO classification of the up- and downstream genes of the top 20 DElncRNAs in the AcCK7 vs. AcT7 (A) and AcCK7 vs. AcT10 (B) comparison groups.
Enlarge this image.Figure 4. KEGG pathways enriched by the up-and downstream genes of the top 20 DElncRNAs in the AcCK7 vs. AcT7 (A) and AcCK10 vs. AcT10 (B) comparison groups.
Enlarge this image.Figure 5. DElncRNA-DEmiRNA-DEmRNA regulatory networks in the AcCK7 vs. AcT7 (A) and AcCK10 vs. AcT10 (B) comparison groups.
Enlarge this image.Figure 6. Secondary structures of four miRNA precursors: (A) two DElncRNAs as miRNA precursors in the AcCK7 vs. AcT7 comparison group; (B) one DElncRNA as miRNA precursor in the AcCK10 vs. AcT10 comparison group.
Enlarge this image.Figure 7. RT-qPCR validation of 12 DElncRNAs. ns, nonsignificant. * p < 0.05, ** p < 0.01, and *** p < 0.001.
Enlarge this image.Figure 8. A hypothetical working model of the DElncRNA-modulated immune response of A. c. cerana workers’ midguts to N. ceranae infestation.
Supplementary Materials
The following supporting information can be downloaded at:
References
1. Kopp, F.; Mendell, J.T. Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell; 2018; 172, pp. 393-407. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29373828]
2. Sinpoo, C.; Paxton, R.J.; Disayathanoowat, T.; Krongdang, S.; Chantawannakul, P. Impact of Nosema ceranae and Nosema apis on individual worker bees of the two host species (Apis cerana and Apis mellifera) and regulation of host immune response. J. Insect Physiol.; 2018; 105, pp. 1-8. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29289505]
3. Wu, Z.; Liu, X.; Liu, L.; Deng, H.; Zhang, J.; Xu, Q.; Cen, B.; Ji, A. Regulation of lncRNA expression. Cell. Mol. Biol. Lett.; 2014; 19, pp. 561-575. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25311814]
4. Murillo-Maldonado, J.M.; Riesgo-Escovar, J.R. The various and shared roles of lncRNAs during development. Dev. Dyn.; 2019; 248, pp. 1059-1069. [DOI: https://dx.doi.org/10.1002/dvdy.108]
5. Ulitsky, I.; Bartel, D.P. LincRNAs: Genomics, evolution, and mechanisms. Cell; 2013; 154, pp. 26-46. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23827673]
6. Choi, S.W.; Kim, H.W.; Nam, J.W. The small peptide world in long noncoding RNAs. Brief. Bioinform.; 2019; 20, pp. 1853-1864.
7. Wu, Y.; Cheng, T.; Liu, C.; Liu, D.; Zhang, Q.; Long, R.; Zhao, P.; Xia, Q. Systematic identification and characterization of long non-coding RNAs in the silkworm, Bombyx mori. PLoS ONE; 2016; 11, e0147147.
8. Chen, X.; Ma, C.; Chen, C.; Lu, Q.; Shi, W.; Liu, Z.; Wang, H.; Guo, H. Integration of lncRNA-miRNA-mRNA reveals novel insights into oviposition regulation in honey bees. PeerJ; 2017; 5, e3881.
9. Chen, D.; Chen, H.; Du, Y.; Zhou, D.; Geng, S.; Wang, H.; Wan, J.; Xiong, C.; Zheng, Y.; Guo, R. Genome-wide identification of long non-coding RNAs and their regulatory networks involved in Apis mellifera ligustica response to Nosema ceranae infection. Insects; 2019; 10, pp. 245-270.
10. Jayakodi, M.; Jung, J.W.; Park, D.; Ahn, Y.J.; Lee, S.C.; Shin, S.Y.; Shin, C.; Yang, T.J.; Kwon, H.W. Genome-wide characterization of long intergenic non-coding RNAs (lincRNAs) provides new insight into viral diseases in honey bees Apis cerana and Apis mellifera. BMC Genomics; 2015; 16, 680.
11. Du, Y.; Feng, R.R.; Wang, J.; Zhu, Z.W.; Zhang, W.D.; Yu, K.J.; Long, Q.; Cai, Z.B.; Xie, Y.L.; Xiong, C.L. et al. Long non-coding RNA response of 6-day-old Apis cerana ccerana larvae to Ascosphaera apis infection. Acta. Microbiol. Sin.; 2021; 61, pp. 1338-1358. (In Chinese)
12. Chen, C.; Wang, H.; Liu, Z.; Chen, X.; Tang, J.; Meng, F.; Shi, W. Population Genomics Provide Insights into the Evolution and Adaptation of the Eastern Honey Bee (Apis cerana). Mol. Biol. Evol.; 2018; 35, pp. 2260-2271. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29931308]
13. Papini, R.; Mancianti, F.; Canovai, R.; Cosci, F.; Rocchigiani, G.; Benelli, G.; Canale, A. Prevalence of the microsporidian Nosema ceranae in honeybee (Apis mellifera) apiaries in Central Italy. Saudi J. Biol. Sci.; 2017; 24, pp. 979-982. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28663691]
14. Chen, D.; Du, Y.; Chen, H.; Fan, Y.; Fan, X.; Zhu, Z.; Wang, J.; Xiong, C.; Zheng, Y.; Hou, C. et al. Comparative Identification of MicroRNAs in Apis cerana cerana Workers’ Midguts in Responseto Nosema ceranae Invasion. Insects; 2019; 10, 258.
15. Zhu, Z.; Wang, J.; Fan, X.; Long, Q.; Chen, H.; Ye, Y.; Zhang, K.; Ren, Z.; Zhang, Y.; Niu, Q. et al. CircRNA-regulated immune responses of asian honey bee workers to microsporidian infection. Front. Genet.; 2022; 13, 1013239.
16. Yoon, J.H.; Abdelmohsen, K.; Gorospe, M. Functional interactions among microRNAs and long noncoding RNAs. Semin. Cell Dev. Biol.; 2014; 34, pp. 9-14. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24965208]
17. Atianand, M.K.; Fitzgerald, K.A. Long non-coding RNAs and control of gene expression in the immune system. Trends Mol. Med.; 2014; 20, pp. 623-631.
18. Lee, J.T. Epigenetic regulation by long noncoding RNAs. Science; 2012; 338, pp. 1435-1439. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23239728]
19. Xie, C.R.; Wang, F.; Zhang, S.; Wang, F.Q.; Zheng, S.; Li, Z.; Lv, J.; Qi, H.Q.; Fang, Q.L.; Wang, X.M. et al. Long Noncoding RNA HCAL Facilitates the Growth and Metastasis of Hepatocellular Carcinoma by Acting as a ceRNA of LAPTM4B. Mol. Ther. Nucleic Acids; 2017; 9, pp. 440-451.
20. Zhang, Y.C.; Liao, J.Y.; Li, Z.Y.; Yu, Y.; Zhang, J.P.; Li, Q.F.; Qu, L.H.; Shu, W.S.; Chen, Y.Q. Genome-wide screening and functional analysis identify a large number of long noncoding RNAs involved in the sexual reproduction of rice. Genome Biol.; 2014; 15, 512. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25517485]
21. Ye, G.; Zhang, L.; Zhou, X. Long noncoding RNAs are potentially involved in the degeneration of virulence in an aphid-obligate pathogen, Conidiobolus obscurus (Entomophthoromycotina). Virulence; 2021; 12, pp. 1705-1716. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34167451]
22. Wang, Y.; Ye, W.; Wang, Y. Genome-wide identification of long non-coding RNAs suggests a potential association with effector gene transcription in Phytophthora sojae. Mol. Plant Pathol.; 2018; 19, pp. 2177-2186.
23. Zhang, S.; Yin, H.; Shen, M.; Huang, H.; Hou, Q.; Zhang, Z.; Zhao, W.; Guo, X.; Wu, P. Analysis of lncRNA-mediated gene regulatory network of Bombyx mori in response to BmNPV infection. J. Invertebr. Pathol.; 2020; 170, 107323. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31926972]
24. Pereiro, P.; Lama, R.; Moreira, R.; Valenzuela-Muñoz, V.; Gallardo-Escárate, C.; Novoa, B.; Figueras, A. Potential Involvement of lncRNAs in the Modulation of the Transcriptome Response to Nodavirus Challenge in European Sea Bass (Dicentrarchus labrax L.). Biology; 2020; 9, 165. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32679770]
25. Cui, J.; Luan, Y.; Jiang, N.; Bao, H.; Meng, J. Comparative transcriptome analysis between resistant and susceptible tomato allows the identification of lncRNA16397 conferring resistance to Phytophthora infestans by co-expressing glutaredoxin. Plant J.; 2017; 89, pp. 577-589.
26. Wang, L.L.; Jin, J.J.; Li, L.H.; Qu, S.H. Long Non-coding RNAs Responsive to Blast Fungus Infection in Rice. Rice; 2020; 13, 77.
27. Zafar, J.; Huang, J.; Xu, X.; Jin, F. Analysis of Long Non-Coding RNA-Mediated Regulatory Networks of Plutella xylostella in Response to Metarhizium anisopliae Infection. Insects; 2022; 13, 916. [DOI: https://dx.doi.org/10.3390/insects13100916]
28. Hoare, M.; Narita, M. Notch and senescence. Adv. Exp. Med. Biol.; 2018; 1066, pp. 299-318.
29. Ikeya, T.; Hayashi, S. Interplay of Notch and FGF signaling restricts cell fate and MAPK activation in the Drosophila trachea. Development; 1999; 126, pp. 4455-4463.
30. Ronkina, N.; Gaestel, M. MAPK-Activated Protein Kinases: Servant or Partner?. Annu. Rev. Biochem.; 2022; 91, pp. 505-540.
31. Shilo, B.Z. The regulation and functions of MAPK pathways in Drosophila. Methods; 2014; 68, pp. 151-159. [DOI: https://dx.doi.org/10.1016/j.ymeth.2014.01.020] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24530508]
32. Adamo, S.A. Why should an immune response activate the stress response? Insights from the insects (the cricket Gryllus texensis). Brain. Behav. Immun.; 2010; 24, pp. 194-200. [DOI: https://dx.doi.org/10.1016/j.bbi.2009.08.003]
33. Erler, S.; Popp, M.; Lattorff, H.M. Dynamics of immune system gene expression upon bacterial challenge and wounding in a social insect (Bombus terrestris). PLoS ONE; 2011; 6, e18126. [DOI: https://dx.doi.org/10.1371/journal.pone.0018126] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21479237]
34. Dostert, C.; Jouanguy, E.; Irving, P.; Troxler, L.; Galiana-Arnoux, D.; Hetru, C.; Hoffmann, J.A.; Imler, J.L. The Jak-STAT signaling pathway is required but not sufficient for the antiviral response of Drosophila. Nat. Immunol.; 2005; 6, pp. 946-953. [DOI: https://dx.doi.org/10.1038/ni1237]
35. Zhang, L.; Xu, W.; Gao, X.; Li, W.; Qi, S.; Guo, D.; Ajayi, O.E.; Ding, S.W.; Wu, Q. LncRNA sensing of a viral suppressor of RNAi activates non-canonical innate immune signaling in Drosophila. Cell Host Microbe; 2020; 27, pp. 115-128. [DOI: https://dx.doi.org/10.1016/j.chom.2019.12.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31917956]
36. Chen, D.; Guo, R.; Xu, X.; Xiong, C.; Liang, Q.; Zheng, Y.; Luo, Q.; Zhang, Z.; Huang, Z.; Kumar, D. et al. Uncovering the immune responses of Apis mellifera ligustica larval gut to Ascosphaera apis infection utilizing transcriptome sequencing. Gene; 2017; 621, pp. 40-50. [DOI: https://dx.doi.org/10.1016/j.gene.2017.04.022]
37. Guo, R.; Chen, D.; Diao, Q.; Xiong, C.; Zheng, Y.; Hou, C. Transcriptomic investigation of immune responses of the Apis cerana cerana larval gut infected by Ascosphaera apis. J. Invertebr. Pathol.; 2019; 166, 107210. [DOI: https://dx.doi.org/10.1016/j.jip.2019.107210]
38. Hillyer, J.F. Insect immunology and hematopoiesis. Dev. Comp. Immunol.; 2016; 58, pp. 102-118.
39. Kingsolver, M.B.; Huang, Z.; Hardy, R.W. Insect antiviral innate immunity: Pathways, effectors, and connections. J. Mol. Biol.; 2013; 425, pp. 4921-4936. [DOI: https://dx.doi.org/10.1016/j.jmb.2013.10.006]
40. Su, G.; Morris, J.H.; Demchak, B.; Bader, G.D. Biological network exploration with Cytoscape 3. Curr. Protoc. Bioinform.; 2014; 47, 8.13.1–8.13.24. [DOI: https://dx.doi.org/10.1002/0471250953.bi0813s47]
41. Robin, M.; Issa, A.R.; Santos, C.C.; Napoletano, F.; Petitgas, C.; Chatelain, G.; Ruby, M.; Walter, L.; Birman, S.; Domingos, P.M. et al. Drosophila p53 integrates the antagonism between autophagy and apoptosis in response to stress. Autophagy; 2019; 15, pp. 771-784. [DOI: https://dx.doi.org/10.1080/15548627.2018.1558001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30563404]
42. Bonasio, R.; Shiekhattar, R. Regulation of transcription by long noncoding RNAs. Annu. Rev. Genet.; 2014; 48, pp. 433-455. [DOI: https://dx.doi.org/10.1146/annurev-genet-120213-092323]
43. Ma, L.; Bajic, V.B.; Zhang, Z. On the classification of long non-coding RNAs. RNA Biol.; 2013; 10, pp. 925-933. [DOI: https://dx.doi.org/10.4161/rna.24604]
44. Yuan, L.; Ren, X.; Zheng, Y.; Qian, J.; Xu, L.; Sun, M. MiR-315 is required for neural development and represses the expression of dFMR1 in Drosophila melanogaster. Biochem. Biophys. Res. Commun.; 2020; 525, pp. 469-476. [DOI: https://dx.doi.org/10.1016/j.bbrc.2020.02.127]
45. Silver, S.J.; Hagen, J.W.; Okamura, K.; Perrimon, N.; Lai, E.C. Functional screening identifies miR-315 as a potent activator of Wingless signaling. Proc. Natl. Acad. Sci. USA; 2007; 104, pp. 18151-18156. [DOI: https://dx.doi.org/10.1073/pnas.0706673104]
46. Avila-Bonilla, R.G.; Yocupicio-Monroy, M.; Marchat, L.A.; Pérez-Ishiwara, D.G.; Cerecedo-Mercado, D.A.; Del Ángel, R.M.; Salas-Benito, J.S. miR-927 has pro-viral effects during acute and persistent infection with dengue virus type 2 in C6/36 mosquito cells. J. Gen. Virol.; 2020; 101, pp. 825-839. [DOI: https://dx.doi.org/10.1099/jgv.0.001441]
47. He, Q.; Zhang, Y.; Dong, W. MicroRNA miR-927 targets the juvenile hormone primary response gene Krüppel homolog1 to control Drosophila developmental growth. Insect. Mol. Biol.; 2020; 29, pp. 545-554. [DOI: https://dx.doi.org/10.1111/imb.12662] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32715555]
48. Kitatani, Y.; Tezuka, A.; Hasegawa, E.; Yanagi, S.; Togashi, K.; Tsuji, M.; Kondo, S.; Parrish, J.Z.; Emoto, K. Drosophila miR-87 promotes dendrite regeneration by targeting the transcriptional repressor Tramtrack69. PLoS Genet.; 2020; 16, e1008942. [DOI: https://dx.doi.org/10.1371/journal.pgen.1008942] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32764744]
49. Sarver, A.L.; Subramanian, S. Competing endogenous RNA database. Bioinformation; 2012; 8, pp. 731-733. [DOI: https://dx.doi.org/10.6026/97320630008731] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23055620]
50. Cesana, M.; Cacchiarelli, D.; Legnini, I.; Santini, T.; Sthandier, O.; Chinappi, M.; Tramontano, A.; Bozzoni, I. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell; 2011; 147, pp. 358-369. [DOI: https://dx.doi.org/10.1016/j.cell.2011.09.028]
51. Liu, F.; Shi, T.; Yin, W.; Su, X.; Qi, L.; Huang, Z.Y.; Zhang, S.; Yu, L. The microRNA ame-miR-279a regulates sucrose responsiveness of forager honey bees (Apis mellifera). Insect. Biochem. Mol. Biol.; 2017; 90, pp. 34-42. [DOI: https://dx.doi.org/10.1016/j.ibmb.2017.09.008] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28941994]
52. Zhang, Y.; Feng, K.; Hu, J.; Shi, L.; Wei, P.; Xu, Z.; Shen, G.; Li, M.; Xu, Q.; He, L. A microRNA-1 gene, tci-miR-1-3p, is involved in cyflumetofen resistance by targeting a glutathione S-transferase gene, TCGSTM4, in Tetranychus cinnabarinus. Insect. Mol. Biol.; 2018; 27, pp. 352-364. [DOI: https://dx.doi.org/10.1111/imb.12375]
53. Peng, W.; Yu, S.; Handler, A.M.; Tu, Z.; Saccone, G.; Xi, Z.; Zhang, H. MiRNA-1-3p is an early embryonic male sex-determining factor in the Oriental fruit fly Bactrocera dorsalis. Nat. Commun.; 2020; 11, 932. [DOI: https://dx.doi.org/10.1038/s41467-020-14622-4]
54. Lim, D.H.; Lee, S.; Han, J.Y.; Choi, M.S.; Hong, J.S.; Lee, Y.S. MicroRNA miR-252 targets mbt to control the developmental growth of Drosophila. Insect. Mol. Biol.; 2019; 28, pp. 444-454. [DOI: https://dx.doi.org/10.1111/imb.12562]
55. Lim, D.H.; Lee, S.; Han, J.Y.; Choi, M.S.; Hong, J.S.; Seong, Y.; Kwon, Y.S.; Lee, Y.S. Ecdysone-responsive microRNA-252-5p controls the cell cycle by targeting Abi in Drosophila. FASEB J.; 2018; 32, pp. 4519-4533. [DOI: https://dx.doi.org/10.1096/fj.201701185RR] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29543534]
56. Verma, P.; Cohen, S.M. MiR-965 controls cell proliferation and migration during tissue morphogenesis in the Drosophila abdomen. Elife; 2015; 4, e07389. [DOI: https://dx.doi.org/10.7554/eLife.07389]
57. Jiang, W.; Wei, Y.; Long, Y.; Owen, A.; Wang, B.; Wu, X.; Luo, S.; Dang, Y.; Ma, D.K. A genetic program mediates cold-warming response and promotes stress-induced phenoptosis in C. elegans. Elife; 2018; 7, e35037. [DOI: https://dx.doi.org/10.7554/eLife.35037]
58. Fu, Z.M.; Gu, X.Y.; Hu, Y.; Zhao, H.D.; Zhu, Z.W.; Zhang, H.Y.; Ji, T.; Niu, Q.S.; Chen, D.F.; Guo, R. Lnc13164 regulates immune response of Apis cerana cerana larvae to Ascosphaera apis infection via ace-miR-4968-y. Acta. Microbiol. Sin.; 2023; 63, pp. 1047-1059.
59. Wu, Y.; Guo, Y.L.; Fan, X.X.; Zhao, H.D.; Zhang, Y.Q.; Guo, S.J.; Jing, X.; Liu, Z.T.; Feng, P.L.; Liu, X.Y. et al. ame-miR-34 Modulates the Larval Body Weight and Immune Response of Apis mellifera Workers to Ascosphara apis Invasion. Int. J. Mol. Sci.; 2023; 24, 1214. [DOI: https://dx.doi.org/10.3390/ijms24021214]
60. Chen, H.; Fan, X.; Zhang, W.; Ye, Y.; Cai, Z.; Zhang, K.; Zhang, K.; Fu, Z.; Chen, D.; Guo, R. Deciphering the CircRNA-Regulated Response of Western Honey Bee (Apis mellifera) Workers to Microsporidian Invasion. Biology; 2022; 11, 1285. [DOI: https://dx.doi.org/10.3390/biology11091285]
61. Kong, L.; Zhang, Y.; Ye, Z.Q.; Liu, X.Q.; Zhao, S.Q.; Wei, L.; Gao, G. CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic. Acids. Res.; 2007; 35, (Suppl. 2), pp. W345-W349. [DOI: https://dx.doi.org/10.1093/nar/gkm391] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17631615]
62. Fan, X.N.; Zhang, S.W. LncRNA-MFDL: Identification of human long non-coding RNAs by fusing multiple features and using deep learning. Mol. Biosyst.; 2015; 11, pp. 892-897. [DOI: https://dx.doi.org/10.1039/C4MB00650J] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25588719]
63. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. EdgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics; 2010; 26, pp. 139-140. [DOI: https://dx.doi.org/10.1093/bioinformatics/btp616]
64. Ye, Y.; Fan, X.; Long, Q.; Wang, J.; Zhang, W.; Cai, Z.; Sun, M.; Gu, X.; Zou, P.; Chen, D. et al. Comprehensive investigation and regulatory function of lncRNAs engaged in western honey bee larval immune response to Ascosphaera apis invasion. Front. Physiol.; 2022; 13, 1082522. [DOI: https://dx.doi.org/10.3389/fphys.2022.1082522]
65. Guo, R.; Chen, D.F.; Huang, Z.J.; Liang, Q.; Xiong, C.L.; Xu, X.J.; Zheng, Y.Z.; Zhang, Z.N.; Xie, Y.L.; Tong, X.Y. et al. Transcriptome analysis of Ascosphaera apis stressing larvel gut of Apis cerana cerana. Acta. Microbiol. Sin.; 2017; 57, pp. 1865-1878. (In Chinese)
66. Guo, R.; Geng, S.H.; Xiong, C.L.; Zheng, Y.Z.; Fu, Z.M.; Wang, H.P.; Du, Y.; Tong, X.Y.; Zhao, H.X.; Chen, D.F. Differential expression analysis of long non-coding RNAs during the developmental process of Apis mellifera ligustica worker’s midgut. Sci. Agric. Sin.; 2018; 51, pp. 3600-3613. (In Chinese)
67. Lavorgna, G.; Guffanti, A.; Borsani, G.; Ballabio, A.; Boncinelli, E. TargetFinder: Searching annotated sequence databases for target genes of transcription factors. Bioinformatics; 1999; 15, pp. 172-173. [DOI: https://dx.doi.org/10.1093/bioinformatics/15.2.172] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10089203]
68. Wu, Y.; Wei, B.; Liu, H.; Li, T.; Rayner, S. MiRPara: A SVM-based software tool for prediction of most probable microRNA coding regions in genome scale sequences. BMC Bioinform.; 2011; 12, 107. [DOI: https://dx.doi.org/10.1186/1471-2105-12-107]
69. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods; 2001; 25, pp. 402-408. [DOI: https://dx.doi.org/10.1006/meth.2001.1262]
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.