Yaping Ye1†
Xiaoxue Fan1†Qi Long1
Jie Wang1
Wende Zhang1Zongbing Cai1Minghui Sun1Xiaoyu Gu1Peiyuan Zou1Dafu Chen1,2
Rui Guo1,2*- 1College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou, Fujian, China
- 2Apitherapy Research Institute, Fujian Agriculture and Forestry University, Fuzhou, Fujian, China
Ascosphaera apis is a fungal pathogen that exclusively infects bee larvae, causing chalkbrood disease, which results in severe damage for beekeeping industry. Long non-coding RNAs (lncRNAs) are versatile regulators in various biological processes such as immune defense and host-pathogen interaction. However, expression pattern and regulatory role of lncRNAs involved in immune response of bee host to A. apis invasion is still very limited. Here, the gut tissues of Apis mellifera ligustica 4-, 5-, and 6-day-old larvae inoculated by A. apis spores (AmT1, AmT2, and AmT3 groups) and corresponding un-inoculated larval guts (AmCK1, AmCK2, and AmCK3 groups) were prepared and subjected to deep sequencing, followed by identification of lncRNAs, analysis of differentially expressed lncRNAs (DElncRNAs), and investigation of competing endogenous RNA (ceRNA) network. In total, 3,746 A. m. ligustica lncRNAs were identified, including 78 sense lncRNAs, 891 antisense lncRNAs, 1,893 intergenic lncRNAs, 346 bidirectional lncRNAs, and 210 intronic lncRNAs. In the 4-, 5-, and 6- comparison groups, 357, 236, and 505 DElncRNAs were discovered. Additionally, 217, 129, and 272 DElncRNAs were respectively predicted to regulate neighboring genes via cis-acting manner, and these targets were associated with a series of GO terms and KEGG pathways of great importance, such as response to stimulus and Jak-STAT signaling pathway. Moreover, 197, 95, and 356 DElncRNAs were observed to target 10, eight, and 21 DEmiRNAs and further target 147, 79, and 315 DEmRNAs, forming complex regulatory networks. Further investigation suggested that these targets were engaged in several key cellular and humoral immune pathways, such as phagosome and MAPK signaling pathway. Ultimately, the expression trends of nine randomly selected DElncRNAs were verified by RT-qPCR, confirming the authenticity and reliability of our transcriptome data. Findings in this current work not only provide candidate DElncRNAs for functional study, but also lay a foundation for unclosing the mechanism underlying DElncRNA-regulated larval immune responses to A. apis invasion.
Introduction
Western honey bees (Apis mellifera) are widely reared in apiaries worldwide and produce a variety of honey bee products, such as honey, royal jelly, and propolis. A. mellifera pollinate a substantial quantity of crops and wild flowers, playing a critical role in ecological balance and food security (Potts et al., 2010). Ascosphaera apis is a widespread fungal pathogen which can not only infect bee larvae to causes chalkbrood disease, but also lead to a sharp reduction in colony population and productivity (Aronstein and Murray, 2010). In the last 2 decades, researchers have conducted an array of studies on the interaction between honey bee larvae and A. apis at the morphological, ethological, or biochemical level (Garrido-Bailón et al., 2013; Chaimanee et al., 2017; Guo et al., 2019b). Our group previously systematically investigated the response of Apis mellifera ligustica larval and the infection of A. apis during chalkbrood (Chen D. et al., 2017). For example, Chen D. et al. (2017) revealed the cellular and humoral immune responses of western honey bee larvae to A. apis invasion based on deep sequencing combined with bioinformatics; Chen D. F. et al. (2017) deciphered the transcriptomic dynamics of A. apis during the infection process in A. m. ligustica larvae.
Long non-coding RNAs (lncRNAs) are a class of ncRNA with a length of more than 200 nt containing two or more exons (Zhang et al., 2019). Similar to mRNA structures, most lncRNAs are transcribed by RNA polymerase II and hence contain the 5′ cap and the 3′ polyA tail (Li J. H. et al., 2012). According to the location with respect to protein-coding genes, lncRNA can be classified into sense lncRNA, antisense lncRNA, intergenic lncRNA, intron lncRNA, and other lncRNA (Jarroux et al., 2017). LncRNA can be capable of exerting functions in versatile manners, namely, cis-acting regulation, trans-acting modulation, miRNA precursor, competing endogenous RNA (ceRNA) network, and translation into protein (Gil and Ulitsky, 2020). Evidences suggest that as pivotal regulator, lncRNA palys an important role in a large number of biological processes, including dose compensation (Carmona et al., 2018), cell cycle (Ma et al., 2021), and host-pathogen interaction (Carnero et al., 2016).
With the rapid development of high-throughput sequencing, abundant of lncRNAs have been identified in animals, plants, and microorganisms (Guo et al., 2018; Zhou R. et al., 2021; Fan et al., 2021). In insects, lncRNA-associated studies have mainly focused on model insects such as fruit flies (Drosophila) and silkworms (Bombyx mori); however, knowledge of lncRNA-associated in other insects, including honey bees, is still limited at present. Zhou H. et al. (2021) discovered that lncRNA-CR11538 in Drosophila restored Toll immunity homeostasis via interaction with the transcription factor Dif/Dorsal, further participating in host innate immunity. Previous works demonstrated that lncRNAs in honey bees were potentially involved in behavior (Feng R. R. et al., 2022), caste development (Humann et al., 2013), and host-pathogen/parasite interaction (Chen et al., 2019a). Our previous studies deciphered the responses of both A. m. ligustica and Apis cerana larvae to A. apis invasion at the mRNA and miRNA levels (Chen D. et al., 2017; Guo et al., 2019a). Accumulating evidences have demonstrated that molecules containing miRNA response elements (MREs) such as mRNA, lncRNA, circRNA, can competitively bind miRNA and further regulating the expression of neighboring genes and biological processes (Ala, 2020). However, until now, there was no reported study on the expression profile and regulatory role of lncRNAs during the A. m. ligustica larval response to A. apis invasion.
In the current work, 4-, 5-, and 6-day-old A. m. ligustica larval guts inoculated with A. apis spores and un-inoculated larval guts were subjected to strand-specific cDNA library construction and deep sequencing, followed by identification of lncRNAs and investigation of the expression profile. The potential regulatory functions of host DElncRNAs were then analyzed in combination with previously obtained small RNA-seq (sRNA-seq) data. To our knowledge, this is the first report of the lncRNA-regulated response of honey bee larvae to A. apis invasion. Our findings will shed light on the mechanism underlying the lncRNA-mediated larval response to A. apis infection and offer novel insights into the interaction between A. m. ligustica larvae and A. apis.
Materials and method
Bee and fungus
A. m. ligustica colonies were reared in the teaching apiary of the College of Animal Sciences (College of Bee Science) at Fujian Agriculture and Forestry University. A. apis spores were prepared following the method developed by Jensen et al. (2013) and stored in the Honey Bee Protection Lab, College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University.
Experimental inoculation and sample preparation
In our previous study, A. m. ligustica larvae were reared in 48-well culture plates in a constant temperature and humidity incubator (Shanghai Yiheng Scientific Instrument Co., Ltd.) as described by Peng et al. (1992). Briefly, the diet was mixed and frozen in smaller aliquots and preheated to 34°C before feeding; 2-day-old larvae were removed from the combs with a transferring tool to 10 μL diet; 3-day-old larvae (n = 9) in the treatment group were fed 20 μl diet containing A. apis spores with a final concentration of 107 spores/mL and then fed once a day with 30 μl (4-day-old), 40 μl (5-day-old), and 50 μl (6-day-old) diet; 3-day-old larvae (n = 9) in the control group were fed once a day with 20 μl (3-day-old), 30 μl (4-day-old), 40 μl (5-day-old), and 50 μl (6-day-old) diet without A. apis spores; gut tissues of 4-, 5-, and 6-day-old larvae were harvested utilizing our previously developed protocol (Guo et al., 2019b), and then frozen in liquid nitrogen and kept at −80°C until deep sequencing and molecular experiments. Gut samples of 4-, 5-, and 6-day-old larvae in the control groups were named the AmCK1 group, AmCK2 group, and AmCK3 group, while those in the treatment groups were named the AmT1 group, AmT2 group, and AmT3 group, respectively. There was one replica in each group, and each group included three gut tissues.
RNA extraction, cDNA library construction, and deep sequencing
The total RNA of the gut samples in the aforementioned six groups was extracted using the TRIzol method (Promega, United States). Oligo (dTs) was used to isolate poly (A) mRNA, which was subsequently fragmented and reverse transcribed using random primers (QIAGEN, Germany). Next, second-strand cDNA was synthesized with RNase H and DNA polymerase I, and the double-strand cDNA was then purified by the QiaQuick PCR extraction kit (QIAGEN, Germany). After agarose gel electrophoresis, the required fragments were purified using a DNA extraction kit (QIAGEN, Germany) followed by enrichment through PCR amplification (NEB, United States). The reaction conditions were set as follows: 98°C for 30 s, followed by 13 cycles of 98°C for 10 s and 65°C for 75 s, and 65°C for 5 s. Ultimately, the 6 cDNA libraries were subjected to deep sequencing by Guangzhou Gene Denovo Biotechnology Co., Ltd. Using the Illumina HiSeqTM 4000 platform. Raw data generated from strand-specific library-based RNA-seq were deposited in the NCBI Sequence Read Archive (SRA) database and linked to BioProject number PRJNA406998.
Quality control of raw data
Fastp software (version 0.18.0) (Chen et al., 2018) was used to perform quality control on the raw data by removing reads that contained adapters, more than 10% N-base or low-quality reads (Q-value≤20) to obtain high-quality clean reads. The obtained clean reads were mapped to the reference genome of A. mellifera (Amel_HAV3.1) using HISAT2 software (Kim et al., 2015) with default parameters.
Identification of lncRNAs
The transcripts were assembled using a combination of Cufflinks (Trapnell et al., 2012) and TopHat2 software (Kim et al., 2013) and then aligned to the A. mellifera reference genome (Amel_HAV3.1) with Cuffcompare software (Meana et al., 2010) to detect novel transcripts. Transcripts with one of the classcodes “u, i, j, x, c, e, o” were defined as novel transcripts. Next, lncRNAs were filtered out according to the following criteria: (1) length ≥ 200 bp and (2) exon number ≥ 2. Furthermore, both CPC (Kong et al., 2007) and CNCI software (Sun et al., 2013) were used to predict the coding ability of the above lncRNAs, and only those without coding ability were regarded as reliable lncRNAs. The numbers of various types of lncRNAs, including intergenic lncRNAs, sense lncRNAs, and antisense lncRNAs, were calculated.
Analysis of DElncRNAs
Following the method described by Chen D. et al. (2017), the transcript expression level was normalized with the FPKM (fragments per kilobase of transcript per million mapped reads) method, which can eliminate the influence of different transcript lengths and sequencing data amounts on the calculation of transcript expression. EdgeR software (Robinson et al., 2010) was employed to screen differentially expressed lncRNAs (DElncRNAs) in the AmCK1 vs. AmT1, AmCK2 vs. AmT2, and AmCK3 vs. AmT3 comparison groups following the criteria of p ≤ 0.05 (corrected by false discovery rate) and |log2(FC)|≥1. Finally, Venn diagram analysis and expression clustering of DElncRNAs in different comparison groups were conducted using the related tools in the OmicShare platform (www.omicshare.com).
Investigation of the cis-acting role of DElncRNAs
Protein-coding genes located 10 kb upstream and downstream of DElncRNAs were surveyed and then annotated to GO (http://www.geneontology.org/) and KEGG (https://www.kegg.jp/) databases utilizing the Blast tool. Subsequently, gene numbers were calculated for each GO term or KEGG pathway. Furthermore, significantly enriched GO terms by neighboring genes were defined by hypergeometric testing, while KEGG pathway enrichment analysis was performed by KOBAS 2.0 software (Xie et al., 2011), with the A. mellifera reference genome (Amel_HAV3.1) as the background. Only terms or pathways with corrected p values of less than 0.05 were considered enriched.
Source of small RNA-seq datasets
In our previous study, gut tissues of A. apis-inoculated 4-, 5-, and 6-day-old larvae and corresponding uninoculated 4-, 5-, and 6-day-old larval guts were prepared following the aforementioned method, followed by total RNA isolation, cDNA library construction, and sRNA-seq by Guangzhou Gene Denovo Biotechnology Co., Ltd. Using the Illumina MiSeq™ 4000 platform. Raw data derived from sRNA-seq are available in the NCBI SRA database under the BioProject number PRJNA406998. Quality control was previously performed, and the results were suggestive of the high quality of sRNA-seq-derived data (Feng W. et al., 2022).
Investigation of ceRNA regulatory networks
According to the method described by Chen et al. (2019a), target DEmiRNAs of DElncRNAs and target DEmRNAs of DEmiRNAs were predicted by using a combination of three software programs, MiRanda (V3.3a) (Enright et al., 2003), RNAhybrid (V2.1.2) (Krüger and Rehmsmeier, 2006) + SVM_light (V6.01) (Kumar et al., 2008), and TargetFind (Allen et al., 2005), and the intersection was deemed reliable targets. Based on the prediction results, ceRNA regulatory networks were constructed and visualized by Cytoscape software (Smoot et al., 2011) with default parameters.
RT-qPCR validation of DElnRNAs
Nine DElncRNAs were randomly selected for RT-qPCR, including three (MSTRG.1133.1, XR_001704875.2, and MSTRG.11613.1) from the AmCK1 vs. AmT1 comparison group, three (MSTRG.4918.2, MSTRG. 9603.5, and MSTRG.8790.2) from AmCK2 vs. AmT2 comparison group, and three (XR_001705688.2, XR_410074.3, and XR_003304187.1) from AmCK3 vs. AmT3 comparison group. Specific forward and reverse primers for each of the nine selected DElncRNAs were respectively designed using Primer Premier 6 (Singh et al., 1998), and synthesized by Sangon Biotech (Shanghai) Co., Ltd. Gene actin (GeneBank ID: 406122) was used as internal reference. Total RNA from gut samples in each group was respectively isolated using FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, China), followed by reverse transcription with Random primers. The resulting cDNA was used as the template for qPCR reaction, which was performed on QuanStudio 3 Fluorescence quantitative PCR instrument (ABI, United States) The reaction system (20 μl) contained 10 μl of SYBR Green Dye, 1 μl of upstream and downstream primers (10 μmol/L), 1 μl of cDNA template, 7 μl of DEPC water. The reaction conditions were set as follows: 95°C pre-denaturation for 5 min; 95°C denaturation for 30 s, 60°C annealing and extension for 30 s, a total of 40 cycles each group of qPCR reaction and set three times experiment for repeating. The relative expression level of each DElncRNA was calculated using 2−ΔΔCT method (Livak and Schmittgen, 2001). All experiments were run with at least three parallel samples and were repeated three times. Data were shown as mean ± standard deviation (SD) and subjected to Student’s t-test by Graph Prism 8 software (San Diego, United States) (p < 0.05 was considered statistically significant). Detailed information of primers used in this work was presented in the Supplementary Table S1.
Result
Quality control of deep sequencing data
Totally, 85,811,046, 81,962,296, 85,636,572, 79,267,686, 82,889,882, and 100,211,796 raw reads were produced from AmCK1, AmCK2, AmCK3, AmT1, AmT2, and AmT3 groups, respectively (Supplementary Table S2). After strict quality control, 85,739,414, 81,896,402, 85,573,798, 79,202,304, 82,828,926, and 100,128,692 clean reads were obtained, respectively (Supplementary Table S2). In brief, the Q20 and Q30 of clean reads in six groups were above 98.07% and 94.30%, respectively (Supplementary Table S2). In addition, the mapping ratio of clean reads in the reference genome was above 99.92% (Supplementary Table S2). The results showed that the next-generation sequencing data were reliable and could be used for further analysis. Compared with the A. mellifera reference genome (Amel_HAV3.1), the comparison rate ranged from 92.06% to 94.81%, and there were 63.10%–69.80% clean reads were compared to the exon region, 8.56%–9.59% of them compared to the intron region, and 20.98%–27.32% compared to the intergene region (Supplementary Table S3).
Identification and investigation of A. m. ligustica lncRNAs
In the aforementioned six groups, 1,991, 2,031, 1,970, 2,101, 1,955, and 1,944 lncRNAs were identified. After removing redundant lncRNAs, a total of 3,746 A. m. ligustica lncRNAs were discovered, including 3,146 known lncRNAs and 600 novel lncRNAs. Among these, there were 78 sense lncRNAs, 891 antisense lncRNAs, 1,893 intergenic lncRNAs, 346 bidirectional lncRNAs, and 210 intron lncRNAs (Figure 1).
FIGURE 1. Type statistics of A. m. ligustica lncRNAs identified in this work.
Differential expression profile of lncRNAs in larval guts infected by A. apis
Here, 156, 98, and 361 up-regulated and 201, 138, and 144 down-regulated lncRNAs were identified in the 4-, 5-, and 6-day-old comparison groups (Figure 2A). Venn analysis indicated that 32 DElncRNAs were shared by the three comparison groups, and the numbers of specific DElncRNAs were 191, 93, and 322, respectively (Figure 2B). As shown in Figure 2C, the shared DElncRNAs displayed various expression patterns during A. apis infection; such as the expression levels of XR_003305757.1, XR_003304308.1, and XR_001705213.2 were elevated in AmCK1 group but decreased in AmT1 group (Figure 2C).
FIGURE 2. Number and expression of DElncRNAs in A.m. ligustica larval guts responding to A. apis infection. (A) Number of up- and down-regulated lncRNAs in three comparison groups. (B) Venn diagram of DElncRNAs in three comparison groups. (C) Heat map of shared DElncRNAs by three comparison groups.
Cis-acting regulation of host DElncRNAs
In the AmCK1 vs. AmT1 comparison group, 217 DElncRNAs were predicted to regulate 361 neighboring genes, which were enriched in 31 GO terms relative to biological process, cellular component, and molecular function, such as behavior, cell, and binding (Figure 3A); these genes were also involved in 203 KEGG pathways, including ribosome biogenesis in eukaryotes, measles, and viral carcinogenesis (Figure 3D); further investigation suggested that 37 genes were associated with immune pathways such as cell apoptosis, autophagy, and endocytosis (Table 1). Comparatively, 129 DElncRNAs in the AmCK2 vs. AmT2 comparison group were predicted to regulate 217 neighboring genes, which were enriched in 33 GO terms (Figure 3B), such as localization and organelle, as well as 154 KEGG pathways, such as alcoholism and cell adhesion molecules (Figure 3E); additionally, 19 genes relevant to immune pathways, such as phagosome, ubiquitin-mediated proteolysis, and Jak-STAT signaling pathway, were putatively regulated by 129 DElncRNAs (Table 1). In the AmCK3 vs. AmT3 comparison group, 272 DElncRNAs were predicted to regulate 496 neighboring genes, which were enriched in 31 GO terms (Figure 3C), such as immune system process and catalytic activity, as well as 226 KEGG pathways, such as malaria and p53 signaling pathway (Figure 3F); 272 DElncRNAs were found to be engaged in regulating 51 genes related to cellular or humoral immune pathways, such as endocytosis, biosynthesis of insect hormones, and lysosome (Table 1).
FIGURE 3. GO terms and KEGG pathways enrich by neighboring genes of DElncRNAs. (A–C) GO classification of neighboring genes of DElncRNAs in AmCK1 vs. AmT1, AmCK2 vs. AmT2, and AmCK3 vs. AmT3 comparison groups. (D–F) Enriched KEGG pathways by neighboring genes of DElncRNAs in AmCK1 vs. AmT1, AmCK2 vs. AmT2, and AmCK3 vs. AmT3 comprison groups.
TABLE 1. Summary of neighboring genes of DElncRNAs associated with immune pathways.
CeRNA regulatory network of host DElncRNAs
Complex ceRNA networks of DElncRNAs existed in the AmCK1 vs. AmT1 and AmCK2 vs. AmT2 comparison groups, whereas the DElncRNA-DEmiRNA-DEmRNA regulatory network in the AmCK3 vs. AmT3 comparison groups was more complicated (Figure 4). In detail, 197 DElncRNAs in the AmCK1 vs. AmT1 comparison group could target 10 DEmiRNAs and further bind to 147 DEmRNAs (Figure 4A), which were engaged in 24 GO terms, including single-organism process and membrane (Figure 5A), as well as 23 KEGG pathways, such as circadian rhythm-fly and dorsoventral axis formation (Figure 5D). In the AmCK2 vs. AmT2 comparison group, 95 DElncRNAs could target eight DEmiRNAs and further link to 79 DEmRNAs (Figure 4B), which were involved in 22 GO terms, including single-organism process and binding (Figure 5B), as well as 16 KEGG pathways, such as ECM-receptor interaction and beta-alanine metabolism (Figure 5E). In the AmCK3 vs. AmT3 comparison group, 356 DElncRNAs could target 21 DEmiRNAs and 315 DEmRNAs (Figure 4C), which were associated with 28 GO terms, including cellular process and membrane part (Figure 5C), as well as 68 KEGG pathways, such as neuroactive ligand and drug metabolism (Figure 5F).
FIGURE 4. (A) DElncRNA-involved ceRNA network in the AmCK1 vs. AmT1 comparison group. (B) DElncRNA-involved ceRNA network in the AmCK2 vs. AmT2 comparison group. (C) DElncRNA-involved ceRNA network in the AmCK3 vs. AmT3 comparison group.
FIGURE 5. GO terms and KEGG pathways enriched by target DEmRNAs within ceRNA networks in three comparison groups. (A–C) Loop graphs of enriched GO terms by targets in AmCK1 vs. AmT1, AmCK2 vs. AmT2, and AmCK3 vs. AmT3 comparison groups. (D–F) Chord graphs of enriched KEGG pathways by targets in AmCK1 vs. AmT1, AmCK2 vs. AmT2, and AmCK3 vs. AmT3 comparison groups.
Further investigation demonstrated that target DEmRNAs in the AmCK1 vs. AmT1 comparison group were related to two cellular immune pathways (lysosome and endocytosis) and one humoral immune pathway (MAPK signaling pathway); target DEmRNAs in the AmCK2 vs. AmT2 comparison group were associated with two cellular immune pathways (endocytosis and phagosome); target DEmRNAs in the AmCK3 vs. AmT3 comparison group were relevant to six cellular immune pathways (lysosome, endocytosis, phagosome, etc.) and one humoral immune pathway (Toll and Imd signaling pathways). Detailed information about immune-related targets is presented in Table 2.
TABLE 2. Summary of immune-associated target DEmRNAs within ceRNA networks.
Validation of DElncRNAs by RT‒qPCR
The RT‒qPCR results suggested that the expression trends of nine DElncRNAs were consistent with those in the transcriptome data (Figure 6), further confirming the reliability of the sequencing data used in this work.
FIGURE 6. RT-qPCR verification of DElncRNAs. (A–C) DElncRNAs in AmCK1 vs. AmT1 comparison groups. (D–F), DElncRNAs in AmCK2 vs. AmT2 comparison groups. (G–I) DElncRNAs in AmCK3 vs. AmT3 comparison groups. * represents p < 0.05, ** represents p < 0.01, *** represents p < 0.001.
Discussion
Previously, our group identified 6,353 lncRNAs in A. m. ligustica workers’ midguts based on RNA-seq, including 4,749 known and 1,604 novel lncRNAs (Chen et al., 2019a). Here, 3,146 known and 600 novel lncRNAs were identified (Figure 1), further enriching the reservoir of lncRNAs in A. mellifera. In addition, we found that 39 (1.2%) known lncRNAs were shared by the A. m. ligustica larval guts and workers’ midguts, indicative of their key roles in the growth and development of both larval gut and worker’s midgut. It is speculated that those specific lncRNAs exert different functions at different developmental stages of the gut tissue. In animals and plants, the lncRNA expression was suggested to be tissue- and stage-specific (Statello et al., 2021). Hence, the total number of A. m. ligustica lncRNAs should be much greater than the documented ones. It is believed that with the increasing quantity of related studies, more honey bee lncRNAs will be discovered in the future.
LncRNAs were verified to be vital regulators in responses of insects to pathogen/parasite infection (Valadkhan and Valencia-Hipólito, 2016; Meng et al., 2021). For instance, Zhang L. et al. (2020) identified 4,450 DElncRNAs, 66 DEmiRNAs, and 7,448 DEmRNAs in the B. mori BmN cells responding to the B. mori nucleopolyhedrosis virus (BmNPV) infection, and found that DElncRNAs were likely to participate in host response via ceRNA network. Jayakodi et al. (2015) identified 2,470 lincRNAs and 1,514 lincRNAs from Apis cerana and A. mellifera, respectively, in sacbrood virus (SBV) infected A. cerana, 11 lincRNAs are specifically regulated upon viral infection. In this work, 156, 98, and 361 up-regulated lncRNAs were detected in the 4-, 5-, and 6-day-old comparison groups (Figure 2A), which showed that a portion of lncRNAs were induced to activation. Comparatively, 201, 138, and 144 lncRNAs were found to down-regulate in the above-mentioned three comparison groups, indicative of the suppression of these lncRNAs by A. apis. In summary, these DElncRNAs were speculated to be engaged in the host response to A. apis infection and play certain regulatory parts. In addition, three up-regulated lncRNAs (XR_001702513.2 (log2FC = 8.748 2, p = 0.000 2), XR_003304325.1 (log2FC = 1.854 1, p = 0.000 4), and XR_411902.3 (log2FC = 2.667 4, p < 0.000 1)) and two down-regulated lncRNAs (XR_001703713.2 (log2FC = −7.965 8, p < 0.000 1) and XR_003305757.1 (log2FC = −7.409 4, p < 0.000 1)) were shared by the aforementioned three comparison groups (Figure 2C), suggestive of their pivotal roles in host response, therefore deserving additional investigation such as RNAi-based functional dissection.
Accumulating evidences have shown that lncRNAs are capable of regulating the transcription of neighboring genes in a cis-acting manner (Gil and Ulitsky, 2020). Chen et al. (2019a) constructed DElncRNA-miRNA-mRNA networks of the A. m. ligustica response to Nosema ceranae infection and found that a portion of DElncRNAs were likely to participate in regulating host material and energy metabolism as well as cellular and humoral immunity during host responses to N. ceranae invasion. Ropri et al. (2021) reported that SE-lncRNAs (RP11-379F4.4 and RP11-465B22.8) played a potential role in the progression of ductal carcinoma in situ (DCIS) and invasive ductal carcinoma (IDC) by regulating the expression of neighboring genes. Insect only consumes a small amount of energy to maintain basic activities without activation of the immune system; however, once the immune system is activated under pathogen invasion, the energy consumption rate will greatly increase (Kingsolver et al., 2013). Here, we observed that the largest number of neighboring genes in the AmCK3 vs. AmT3 comparison group were involved in energy metabolism-associated pathways such as oxidative phosphorylation and nitrogen metabolism, indicating the continuous proliferation of A. apis within the host larval gut and the A. apis-caused pressure at the later stage of infection, which resulted in the enhancement of host immune defense and thus the elevation of the energy metabolism rate. To cope with the infection of pathogenic microorganisms, insects have evolved an efficient innate immune system including cellular and humoral immune, and the latter is mediated by multiple signaling pathways (Valanne et al., 2011; Hillyer, 2016). In the present study, neighboring genes regulated by DElncRNAs in the 4-, 5- and 6-day-old larval guts infected by A. apis were observed to be involved in several cellular immune pathways namely cell apoptosis, lysosome, phagosome, ubiquitin-mediated proteolysis, and some humoral immune pathways such as Jak-STAT, MAPK as well as Toll and Imd signaling pathways (Figure 3). These results demonstrated that corresponding DElncRNAs were likely to control the transcription of neighboring genes and further participate in regulating host cellular and humoral immune responses to A. apis infection.
Those lncRNAs with MREs could interact with miRNAs and influence downstream gene expression via ceRNA networks (Tang et al., 2017; Wang et al., 2020). An increasing number of lncRNAs have been verified to be pivotal regulators in the occurrence and development of an array of human diseases through ceRNA mechanisms, such as cancer, Alzheimer’s disease, and cardiovascular disease (Wu et al., 2020; Ma et al., 2020). Additionally, lncRNAs were suggested to mediate insect-pathogen interactions via ceRNA regulatory networks (Chen et al., 2019a; Zhang S. et al., 2020; Mao et al., 2022). For example, Mao et al. (2022) reported that aae-lnc-0165 suppressed by Wolbachia induced the expression of REL1 gene in Aedes aegypti through the sequence-specific binding of aae-miR-980-5p, which contributed to the activation of Toll signaling pathway. Previous findings showed that lncRNA-mediated ceRNA regulatory networks were putatively engaged in midgut development and the N. ceranae response of western honey bee workers (Guo et al., 2018; Chen et al., 2019a). In this current work, 197, 95, and 356 DElncRNAs in 4-, 5-, and 6-day-old larval guts could respectively target 10, eight, and 21 DEmiRNAs and further target 147, 79, and 315 DEmRNAs, forming complex ceRNA regulatory networks (Figure 4). This indicated that these DElncRNAs may participated in the regulation larval response to A. apis invasion via ceRNA mechanism. After being ingested by bee larvae, the A. apis spores enter into the midgut and germinate at low level, the diaphragm between midgut and hindgut disappears at prepupal stage (7- and 8-day-old), and the spores rapidly germinate and at the meantime the mycelia grow in abundance in the hindgut when contacting O2, thereafter the mycelia penetrate the gut wall and then the body wall, resulting in chalkbrood mummy (Li T. et al., 2012; Jensen et al., 2013). Intriguingly, it is observed that the number of immune-associated targets within ceRNA networks in the 6-day-old comparison group were more than other two comparison groups. This reflected that with the increased time of A. apis infection the host-pathogen interactions were at a higher level at 6-day-old (3 dpi), a timepoint adjacent to the outbreak of chalkbrood disease.
MiR-1 plays an important role in the pathogenesis of heart disease, Liu et al. (2021) indicated that inhibiting miR-1 may relieve right ventricle hypertrophy and fibrosis in model rats used in their research, which also works significantly in plants. Wang et al. (2018) discovered that Connexin43, an extract of Astragalus root, worked by targeting miR-1 to cure viral myocarditis, and overexpression of miR-1 inhibited endogenous Connexin43 expression significantly. Previous studies have found that the expression of miR-1 in A. m. ligustica workers was significantly downregulated 6 days after N. ceranae inoculation, suggesting its potential involvement in the host immune response (Huang et al., 2015). In the Asian honey bee A. ceranae, Chen et al. (2019b) found that the expression of miR-1-x in the worker’s midgut was significantly downregulated at 7 days post inoculation with N. ceranae spores. Here, 35 DElncRNAs in the 6-day-old larval gut could jointly target miR-1-z (highly homologous to ame-miR-1), which can further target 32 DEmRNAs (Figure 4C). It is speculated that these DElncRNAs potentially regulate the expression of downstream target genes by targeting miR-1-z, further modulating the larval immune response to A. apis invasion. Effective knockdown of lncRNAs in insects such as Drosophila, Helicoverpa armigera, and Plutella xylostella was achieved utilizing the RNAi method (Guan et al., 2020; Zhang L. et al., 2020). Recently, our team conducted dsRNA-based knockdown of lncRNA13164 in the A. ceranae larval guts and found that lncRNA13164 regulated the expression of three immune genes (stk, e3µl and or1) via ace-miR-4968 and further mediated the host immune response to A. apis infection (Fu et al., 2022). In the near future, we will perform a functional study on miR-1-z as well as associated DElncRNAs and explore the mechanism underlying host response mediated by the DElncRNA-miR-1-z-DEmRNA axis.
Conclusion
In a nutshell, 3,146 known lncRNAs and 600 novel lncRNAs were identified in the A. mellifera larval guts; additionally, A. apis infection caused overall change of expression profile of lncRNAs in host guts; DElncRNAs potentially participated in larval immune response to A. apis invasion by regulating the expression of neighboring genes or interacting with DEmiRNAs (Figure 7); corresponding DElncRNAs were potentially engaged in host immune response through ceRNA regulatory networks via absorption of miR-1-z.
FIGURE 7. A working model of DElncRNA-modulated larval immune response of A. m. ligustica bee to A. apis invasion.
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 in the article/Supplementary Material.
Ethics statement
The study was reviewed and approved by Fujian Agriculture and Forestry University.
Author contributions
RG, XF, and JW conceived and planned the experiments. YY, QL, WZ, and ZC carried out the experiments, analyzed and interpreted the data. RG, MS, XG, and PZ designed the figures. DC and RG reviewed and edited the paper. All authors contributed to the review and approval of the manuscript for publication.
Funding
The National Natural Science Foundation of China (31702190), the Earmarked Fund for China Agriculture Research System (CARS-44-KXJ7), the Master Supervisor Team Fund of Fujian Agriculture and Forestry University (RG), the Natural Science Foundation of Fujian Province (2022J01131334), the Scientific Research Project of College of Animal Sciences (College of Bee Science) of Fujian Agriculture and Forestry University (RG), and the Undergraduate Innovation and Entrepreneurship Training Program of Fujian province (202210389114, 202210389131).
Acknowledgments
All authors thanks reviewers and editors for their constructive comments and recommendations. RG appreciates the love from his beloved wife and daughter.
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/fphys.2022.1082522/full#supplementary-material
References
Ala U. (2020). Competing endogenous RNAs, non-coding RNAs and diseases: An intertwined story. Cells 9 (7), 1574. doi:10.3390/cells9071574
Allen E., Xie Z., Gustafson A. M., Carrington J. C. (2005). MicroRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell 121 (2), 207–221. doi:10.1016/j.cell.2005.04.004
Aronstein K. A., Murray K. D. (2010). Chalkbrood disease in honey bees. J. Invertebr. pathology 103 (1), S20–S29. doi:10.1016/j.jip.2009.06.018
Carmona S., Lin B., Chou T., Arroyo K., Sun S. (2018). LncRNA Jpx induces Xist expression in mice using both trans and cis mechanisms. PLoS Genet. 14 (5), e1007378. doi:10.1371/journal.pgen.1007378
Carnero E., Barriocanal M., Prior C., Pablo Unfried J., Segura V., Guruceaga E., et al. (2016). Long noncoding RNA EGOT negatively affects the antiviral response and favors HCV replication. EMBO Rep. 17 (7), 1013–1028. doi:10.15252/embr.201541763
Chaimanee V., Thongtue U., Sornmai N., Songsri S., Pettis J. S. (2017). Antimicrobial activity of plant extracts against the honeybee pathogens, Paenibacillus larvae and Ascosphaera apis and their topical toxicity to Apis mellifera adults. J. Appl. Microbiol. 123 (5), 1160–1167. doi:10.1111/jam.13579
Chen D., Chen H., Du Y., Zhou D., Geng S., Wang H., et al. (2019a). Genome-Wide identification of long non-coding RNAs and their regulatory networks involved in Apis mellifera ligustica response to Nosema ceranae Infection. Insects 10 (8), 245. doi:10.3390/insects10080245
Chen D., Du Y., Chen H., Fan Y., Fan X., Zhu Z., et al. (2019b). Comparative identification of microRNAs in Apis cerana cerana workers’ midguts in responseto Nosema ceranae invasion. Insects 10 (9), 258. doi:10.3390/insects10090258
Chen D. F., Guo R., Xiong C. L., Liang Q., Zheng Y. Z., Jian X. U., et al. (2017). Transcriptomic analysis of Ascosphaera apis stressing larval gut of Apis mellifera ligustica. Acta Entomol. Sin. (in Chinese). doi:10.16380/j.kcxb.2017.04.005
Chen D., Guo R., Xu X., Xiong C., Liang Q., Zheng Y., et al. (2017). Uncovering the immune responses of Apis mellifera ligustica larval gut to Ascosphaera apis infection utilizing transcriptome sequencing. Gene 621, 40–50. doi:10.1016/j.gene.2017.04.022
Chen S., Zhou Y., Chen Y., Gu J. (2018). Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinforma. Oxf. Engl. 34 (17), i884–i890. doi:10.1093/bioinformatics/bty560
Enright A. J., John B., Gaul U., Tuschl T., Sander C., ang Marks D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5 (1), R1. doi:10.1186/gb-2003-5-1-r1
Fan X., Wang D., Shen X., Qiu J., Wu L., Yan J., et al. (2021). Identification of lncRNA expression profiles and analysis of ceRNA in the hippocampus of perinatal glyphosate-exposed mice. Int. J. Dev. Neurosci. 81 (4), 312–323. doi:10.1002/jdn.10102
Feng R. R., Fu Z. M., Du Y., Zhang W. D., Fan X. X., Wang H. P., et al. (2022). Identification and analysis of microRNAs in the larval gut of Apis cerana cerana. Sci. Agric. Sin. 55 (01), 208–218. (in Chinese). doi:10.3864/j.issn.0578-1752.2022.01.017
Feng W., Huang J., Zhang Z., Nie H., Lin Y., Li Z., et al. (2022). Understanding of waggle dance in the honey bee (Apis mellifera) from the perspective of long non-coding RNA. Insects 13 (2), 111. doi:10.3390/insects13020111
Fu Z. M., Gu X. Y., Hu Y., Zhao H. D., Zhu Z. W., Zhang H. Y., et al. (2022). Microbiology. (unpublished).Lnc13164 regulates immune response of Apis cerana cerana larvae to Ascosphaera apis infection via ace-miR-4968-y
Garrido-Bailón E., Higes M., Martínez-Salvador A., Antúnez K., Botías C., Meana A., et al. (2013). The prevalence of the honeybee brood pathogens Ascosphaera apis, Paenibacillus larvae and Melissococcus plutonius in Spanish apiaries determined with a new multiplex PCR assay. Microb. Biotechnol. 6 (6), 731–739. doi:10.1111/1751-7915.12070
Gil N., Ulitsky I. (2020). Regulation of gene expression by cis-acting long non-coding RNAs. Nat. Rev. Genet. 21 (2), 102–117. doi:10.1038/s41576-019-0184-5
Guan R., Li H., Zhang H., An S. (2020). Comparative analysis of dsRNA-induced lncRNAs in three kinds of insect species. Arch. Insect Biochem. Physiol. 103 (1), e21640. doi:10.1002/arch.21640
Guo R., Chen D., Diao Q., Xiong C., Zheng Y., Hou C. (2019b). Transcriptomic investigation of immune responses of the Apis cerana cerana larval gut infected by Ascosphaera apis. J. Invertebr. pathology 166, 107210. doi:10.1016/j.jip.2019.107210
Guo R., Chen D., Xiong C., Hou C., Zheng Y., Fu Z., et al. (2018). Identification of long non-coding RNAs in the chalkbrood disease pathogen Ascospheara apis. J. Invertebr. pathology 156, 1–5. doi:10.1016/j.jip.2018.06.001
Guo R., Yu D., Zhou N. H., Liu S. Y., Xiong C. L., Zheng Y. Z., et al. (2019a). Comprehensive analysis of differentially expressed microRNAs and their target genes in the larval gut of Apis mellifera ligustica during the late stage of Ascosphaera apis stress. Acta Entomol. Sin. (in Chinese). doi:10.16380/j.kcxb.2019.01.006
Hillyer J. F. (2016). Insect immunology and hematopoiesis. Dev. Comp. Immunol. 58, 102–118. doi:10.1016/j.dci.2015.12.006
Huang Q., Chen Y., Wang R. W., Schwarz R. S., Evans J. D. (2015). Honey bee microRNAs respond to infection by the microsporidian parasite Nosema ceranae. Sci. Rep. 5, 17494. doi:10.1038/srep17494
Humann F. C., Tiberio G. J., Hartfelder K. (2013). Sequence and expression characteristics of long noncoding RNAs in honey bee caste development-potential novel regulators for transgressive ovary size. PloS one 8 (10), e78915. doi:10.1371/journal.pone.0078915
Jarroux J., Morillon A., Pinskaya M. (2017). History, discovery, and classification of lncRNAs. Adv. Exp. Med. Biol. 1008, 1–46. doi:10.1007/978-981-10-5203-3_1
Jayakodi M., Jung J. W., Park D., Ahn Y. J., Lee S. C., Shin S. Y., et al. (2015). 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 Genom 16, 680. doi:10.1186/s12864-015-1868-7
Jensen A. B., Aronstein K., Flores J. M., Vojvodic S., Palacio M. A., Spivak M. (2013). Standard methods for fungal brood disease research. J. Apic. Res. 52 (1), 1–20. doi:10.3896/IBRA.1.52.1.1310.3896/IBRA.1.52.1.13
Kim D., Langmead B., Salzberg S. L. (2015). Hisat: A fast spliced aligner with low memory requirements. Nat. methods 12 (4), 357–360. doi:10.1038/nmeth.3317
Kim D., Pertea G., Trapnell C., Pimentel H., Kelley R., Salzberg S. L. (2013). TopHat2: Accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14 (4), R36. doi:10.1186/gb-2013-14-4-r36
Kingsolver M. B., Huang Z., Hardy R. W. (2013). Insect antiviral innate immunity: Pathways, effectors, and connections. J. Mol. Biol. 425 (24), 4921–4936. doi:10.1016/j.jmb.2013.10.006
Kong L., Zhang Y., Ye Z. Q., Liu X. Q., Zhao S. Q., Wei L., et al. (2007). CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic acids Res. 35, W345–W349. Web Server issue). doi:10.1093/nar/gkm391
Krüger J., Rehmsmeier M. (2006). RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic acids Res. 34, W451–W454. Web Server issue): W451–W454. doi:10.1093/nar/gkl243
Kumar M., Gromiha M. M., Raghava G. P. (2008). Prediction of RNA binding sites in a protein using SVM and PSSM profile. Proteins 71 (1), 189–194. doi:10.1002/prot.21677
Li J. H., Zheng Z. Y., Chen D. F., Liang Q. (2012). Factors influencing Ascosphaera apis infection on honeybee larvae and observation on the infection process. Acta Entomol. Sin. 55 (07), 790–797. doi:10.16380/j.kcxb.2012.07.003
Li T., Wang S., Wu R., Zhou X., Zhu D., Zhang Y. (2012). Identification of long non-protein coding RNAs in chicken skeletal muscle using next generation sequencing. Genomics 99 (5), 292–298. doi:10.1016/j.ygeno.2012.02.003
Liu Y., Li Y., Li J., Zuo X., Cao Q., Xie W., et al. (2021). Inhibiting miR-1 attenuates pulmonary arterial hypertension in rats. Mol. Med. Rep. 23 (4), 283. doi:10.3892/mmr.2021.11922
Livak K. J., Schmittgen T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif.) 25 (4), 402–408. doi:10.1006/meth.2001.1262
Ma N., Tie C., Yu B., Zhang W., Wan J. (2020). Identifying lncRNA-miRNA-mRNA networks to investigate Alzheimer’s disease pathogenesis and therapy strategy. Aging 12 (3), 2897–2920. doi:10.18632/aging.102785
Ma Z., Gao X., Shuai Y., Wu X., Yan Y., Xing X., et al. (2021). EGR1-mediated linc01503 promotes cell cycle progression and tumorigenesis in gastric cancer. Cell Prolif. 54 (1), e12922. doi:10.1111/cpr.12922
Mao W., Zeng Q., She L., Yuan H., Luo Y., Wang R., et al. (2022). Wolbachia utilizes lncRNAs to activate the anti-dengue Toll pathway and balance reactive oxygen species stress in Aedes aegypti through a competitive endogenous RNA network. Front. Cell. Infect. Microbiol. 11, 823403. doi:10.3389/fcimb.2021.823403
Meana A., Martín-Hernández R., Higes M. (2010). The reliability of spore counts to diagnose Nosema ceranae infections in honey bees. J. Apic. Res. 49, 212–214. doi:10.3896/ibra.1.49.2.12
Meng X., Li A., Yu B., Li S. (2021). Interplay between miRNAs and lncRNAs: Mode of action and biological roles in plant development and stress adaptation. Comput. Struct. Biotechnol. J. 19, 2567–2574. doi:10.1016/j.csbj.2021.04.062
Peng Y. C., Mussen E., Fong A., Montague M. A., Tyler T. (1992). Effects of chlortetracycline of honey bee worker larvae reared in vitro. J. Invertebr. Pathology 60 (2), 127–133. doi:10.1016/0022-2011(92)90085-I
Potts S. G., Biesmeijer J. C., Kremen C., Neumann P., Schweiger O., Kunin W. E. (2010). Global pollinator declines: Trends, impacts and drivers. Trends Ecol. Evol. 25 (6), 345–353. doi:10.1016/j.tree.2010.01.007
Robinson M. D., McCarthy D. J., Smyth G. K. (2010). EdgeR: A bioconductor package for differential expression analysis of digital gene expression data. Bioinforma. Oxf. Engl. 26 (1), 139–140. doi:10.1093/bioinformatics/btp616
Ropri A. S., DeVaux R. S., Eng J., Chittur S. V., Herschkowitz J. I. (2021). Cis-acting super-enhancer lncRNAs as biomarkers to early-stage breast cancer. Breast cancer Res. BCR 23 (1), 101. doi:10.1186/s13058-021-01479-8
Singh V. K., Mangalam A. K., Dwivedi S., Naik S. (1998). Primer premier: Program for design of degenerate primers from a protein sequence. BioTechniques 24 (2), 318–319. doi:10.2144/98242pf02
Smoot M. E., Ono K., Ruscheinski J., Wang P. L., Ideker T. (2011). Cytoscape 2.8: New features for data integration and network visualization. Bioinforma. Oxf. Engl. 27 (3), 431–432. doi:10.1093/bioinformatics/btq675
Statello L., Guo C. J., Chen L. L., Huarte M. (2021). Gene regulation by long non-coding RNAs and its biological functions. Nat. Rev. Mol. cell Biol. 22 (2), 96–118. doi:10.1038/s41580-020-00315-9
Sun L., Luo H., Bu D., Zhao G., Yu K., Zhang C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic acids Res. 41 (17), e166. doi:10.1093/nar/gkt646
Tang W., Ji M., He G., Yang L., Niu Z., Jian M., et al. (2017). Silencing CDR1as inhibits colorectal cancer progression through regulating microRNA-7. OncoTargets Ther. 10, 2045–2056. doi:10.2147/OTT.S131597
Trapnell C., Roberts A., Goff L., Pertea G., Kim D., Kelley D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7 (3), 562–578. doi:10.1038/nprot.2012.016
Valadkhan S., Valencia-Hipólito A. (2016). “lncRNAs in stress response,” in Long non-coding RNAs in human disease. Editor K. V. Morris (Cham, Switzerland: Springer International Publishing), 203–236.
Valanne S., Wang J. H., Rämet M. (2011). The Drosophila Toll signaling pathway. J. Immunol. 186 (2), 649–656. doi:10.4049/jimmunol.1002302
Wang J., Chen M. Y., Chen J. F., Ren Q. L., Zhang J. Q., Cao H., et al. (2020). LncRNA IMFlnc1 promotes porcine intramuscular adipocyte adipogenesis by sponging miR-199a-5p to up-regulate CAV-1. BMC Mol. cell Biol. 21 (1), 77. doi:10.1186/s12860-020-00324-8
Wang Y., Li J., Xuan L., Liu Y., Shao L., Ge H., et al. (2018). Astragalus Root dry extract restores connexin43 expression by targeting miR-1 in viral myocarditis. Phytomedicine Int. J. phytotherapy Phytopharm. 46, 32–38. doi:10.1016/j.phymed.2018.06.031
Wu X., Sui Z., Zhang H., Wang Y., Yu Z. (2020). Integrated analysis of lncRNA-mediated ceRNA network in lung adenocarcinoma. Front. Oncol. 10, 554759. doi:10.3389/fonc.2020.554759
Xie C., Mao X., Huang J., Ding Y., Wu J., Dong S., et al. (2011). Kobas 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic acids Res. 39, W316–W322. Web Server issue. doi:10.1093/nar/gkr483
Zhang L., Xu W., Gao X., Li W., Qi S., Guo D., et al. (2020). LncRNA sensing of a viral suppressor of RNAi activates non-canonical innate immune signaling in Drosophila. Cell host microbe 27 (1), 115–128.e8. doi:10.1016/j.chom.2019.12.006
Zhang S., Yin H., Shen M., Huang H., Hou Q., Zhang Z., et al. (2020). Analysis of lncRNA-mediated gene regulatory network of Bombyx mori in response to BmNPV infection. J. Invertebr. pathology 170, 107323. doi:10.1016/j.jip.2020.107323
Zhang X., Hong R., Chen W., Xu M., Wang L. (2019). The role of long noncoding RNA in major human disease. Bioorg. Chem. 92, 103214. doi:10.1016/j.bioorg.2019.103214
Zhou H., Li S., Wu S., Jin P., Ma F. (2021). LncRNA-CR11538 decoys Dif/Dorsal to reduce antimicrobial peptide products for restoring Drosophila Toll immunity homeostasis. Int. J. Mol. Sci. 22 (18), 10117. doi:10.3390/ijms221810117
Keywords: honey bee, Apis melliferaligustica, chalkbrood, Ascosphaera apis, ncRNA, lncRNA, ceRNA
Citation: Ye Y, Fan X, Long Q, Wang J, Zhang W, Cai Z, Sun M, Gu X, Zou P, Chen D and Guo R (2022) Comprehensive investigation and regulatory function of lncRNAs engaged in western honey bee larval immune response to Ascosphaera apis invasion. Front. Physiol. 13:1082522. doi: 10.3389/fphys.2022.1082522
Received: 28 October 2022; Accepted: 01 December 2022;
Published: 16 December 2022.
Edited by:
Neal Silverman, Worcester Foundation for Biomedical Research, United StatesReviewed by:
Aili Wang, Shandong Agricultural University, ChinaKui Kang, Zunyi Normal University College, China
Copyright © 2022 Ye, Fan, Long, Wang, Zhang, Cai, Sun, Gu, Zou, Chen and Guo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Rui Guo, ruiguo@fafu.edu.cn
†These authors have contributed equally to this work