Identification of long noncoding RNAs reveals the effects of dinotefuran on the brain in Apis mellifera (Hymenopptera: Apidae)

Dinotefuran (CAS No. 165252–70-0), a neonicotinoid insecticide, has been used to protect various crops against invertebrate pests and has been associated with numerous negative sublethal effects on honey bees. Long noncoding RNAs (lncRNAs) play important roles in mediating various biological and pathological processes, involving transcriptional and gene regulation. The effects of dinotefuran on lncRNA expression and lncRNA function in the honey bee brain are still obscure. Through RNA sequencing, a comprehensive analysis of lncRNAs and mRNAs was performed following exposure to 0.01 mg/L dinotefuran for 1, 5, and 10 d. In total, 312 lncRNAs and 1341 mRNAs, 347 lncRNAs and 1458 mRNAs, and 345 lncRNAs and 1155 mRNAs were found to be differentially expressed (DE) on days 1, 5 and 10, respectively. Gene set enrichment analysis (GSEA) indicated that the dinotefuran-treated group showed enrichment in carbohydrate and protein metabolism and immune-inflammatory responses such as glycine, serine and threonine metabolism, pentose and glucuronate interconversion, and Hippo and transforming growth factor-β (TGF-β) signaling pathways. Moreover, the DE lncRNA TCONS_00086519 was shown by fluorescence in situ hybridization (FISH) to be distributed mainly in the cytoplasm, suggesting that it may serve as a competing endogenous RNA and a regulatory factor in the immune response to dinotefuran. This study characterized the expression profile of lncRNAs upon exposure to neonicotinoid insecticides in young adult honey bees and provided a framework for further study of the role of lncRNAs in honey bee growth and the immune response.


Background
Neonicotinoid pesticides are the most widely used and effective insecticides against multiple herbivorous insects worldwide [1]. Neonicotinoid pesticides are systemic insecticides that act as a nicotinic acetylcholine receptor (nAChR) agonists, which make insects extremely prone to death, mainly by binding to nAChRs [2]. However, neonicotinoid insecticides can have negative effects on the survival and health of nontarget beneficial insects, such as honey bees, which are important pollinators. Forager honey bees are exposed to neonicotinoid pesticides by collecting contaminated pollen and nectar, which are then stored in beehives [3]. Therefore, nurses, young in-hive honey bees and larvae are also exposed to neonicotinoid pesticides through the acquisition of nectar, pollen and bee bread. Sublethal doses of neonicotinoid pesticides have been shown to impair the development of the brain and mushroom bodies [4], affect olfactory learning and memory abilities [5,6], and disrupt navigation by forager honey bees [7]. Neonicotinoid accumulation in the bee brain disrupts circadian rhythmicity and impairs sleep in many bees [8]. Dinotefuran belongs to the third -generation of neonicotinoid insecticides and shows broad-spectrum, systemic insecticidal activity [9], which is highly toxic to A. mellifera [10]. Dinotefuran residues have been detected in nectar, pollen and bee bread, with measured residues in pollen of approximately 147 ng/g [11]. Studies have shown that treatment with a sublethal dose of dinotefuran affects olfaction, octopamine concentrations, learning and homing ability in honey bees [12,13].
Long noncoding RNAs (lncRNAs) are a class of nonprotein-coding RNAs of more than 200 nucleotides in length that have been found in numerous species [14]. It is increasingly clear that lncRNAs act in association with other molecules as regulators in several different physiological processes, immune responses, the pathogenesis of various human diseases and other biological processes [15][16][17]. A series of studies have shown that lncRNAs play important roles in growth, development, caste differentiation, and innate immunity in honey bees. The noncoding RNA Nb-1 in the worker brain is involved in the synthesis and secretion of juvenile hormone, indicating a potential role in social behaviors [18]. Fernanda et al. first observed that lncov1 was overexpressed in worker ovaries and demonstrated an expression peak coinciding with the onset of autophagic cell death [19]. LncRNAs in honey bees show tissuespecific expression and are preferentially expressed in ovarian and fat body tissues, suggesting that they are associated with biological and hormone signaling pathways and various diseases of honey bees [20]. Previous studies have shown the involvement of lncRNAs in the regulation of host-pathogen interactions, such as the responses to Nosema ceranae infection [21], or Ascospheara apis infection [22] and virus infection in honey bees [23]. However, the potential role of lncRNAs in the response of A. mellifera to dinotefuran has not been reported.
In the present study, we performed RNA-seq to detect the profile of lncRNAs in the brains of young adult bees exposed to sublethal doses of dinotefuran and identified differentially expressed (DE) transcripts. To further investigate the function of DE lncRNAs in the honey bee brain, the distribution of the highly expressed and DE lncRNA TCONS_00086519 was determined through subcellular localization. The results not only lay a foundation for elucidating the molecular mechanisms underlying the effects of exposure to neonicotinoid insecticides, but also offer a beneficial resource for the functional study of key dinotefuran-responsive lncRNAs in further studies.

Characterization of the brain tissue transcriptome
Approximately 100 million raw reads per sample were obtained from the 18 cDNA libraries (Additional file 1). After removing adaptor sequences and low-quality reads, more than 82.35% of clean reads were mapped to the A. mellifera reference genome using Hisat2 (Table S1), and approximately 74% of clean reads were aligned with unique loci. Most reads were aligned to exon regions while nearly 2.37-6.54% and 12.52-23.73% originated from introns and intergenic regions, respectively (Additional file 1). However, approximately 10% of the clean reads were mapped outside of annotated regions.

Known mRNA profiling in honey bee brains
In total, 19,837 known mRNAs were identified according to the reference honey bee genome. Most of these mRNAs (19,336) were expressed in both the dinotefuran-treated and untreated control groups, while 275 were specifically expressed in the dinotefurantreated group, and 226 were specifically expressed in the control group (Fig. 1a). However, most of the treatment specific mRNAs exhibited a low expression level. Transcripts from the major royal jelly protein family (Mrjp1, Mrjp2, Mrjp3), chemosensory proteins (CSP1, CSP3), apidermin (APD-2, APD-3) and NADH dehydrogenase subunit 1(ND1) were the most abundant at the three examined time points (Additional file 2). Among the identified transcripts, immune-related transcripts such as defensin 1 (Def1), cyclooxygenase 1 (COX1), abaecin (LOC406144) and hymenoptaecin (LOC406142) were highly expressed in bee brains. Gene ontology (GO) analysis showed that the highly expressed transcripts (FPKM > 1000) were annotated with cellular component-related terms such as extracellular region and ribosome, biological process-related terms such as defense response to bacterium and innate immune response, and molecular function-related terms such as structural constituent of ribosome, suggesting that dinotefuran treatment may cause high energy consumption and an immune response.

Identification of lncRNA
After transcript assembly, lncRNAs were defined according to a series of filtering criteria (Additional file 3) using Cuffmerge and Cuffcompare software. Through coding potential analysis with Coding Potential Calculator 2 (CPC2), PfamScan and Coding-Non-Coding Index (CNCI), 6824 lncRNA transcripts were identified in the honey bee brain transcriptome. Among these transcripts, 186 were specifically expressed in the dinotefurantreated group, and 177 were specifically expressed in the control group (Fig. 1b). Although approximately 90% of lncRNA transcripts exhibited a low expression level (PFKM < 1), transcripts such as TCONS_00024915 and TCONS_00086477 were highly expressed in the honey bee brain (Additional file 4), suggesting that these lncRNAs may play a specific role in in the effects of neonicotinoid insecticides on honey bees.

Comparative features of mRNAs and lncRNAs
Most lncRNAs identified in this study contained fewer exons (Fig. 2a) and were shorter in length than mRNAs (Fig. 2b), which is consistent with the results of previous studies [24,25]. More than 63.34% of the lncRNAs contained two to three exons, while only 11.32% of the mRNAs contained two to three exons. Approximately 40% of the lncRNAs ranged from 300 to 1200 bp in length, which was much higher than the percentage of mRNAs (14%) identified in this size range. Based on reference gene structures and predictions from EMBOSS explorer (http://emboss.bioinformatics.nl/), we found that the open reading frames (ORFs) of the mRNAs were much longer than those of the lncRNAs in both annotated and novel transcripts (Fig. 2c). In addition, the expression levels of the lncRNAs were found to be much lower than those of the mRNAs in each library (Fig. 2d).

GO and KEGG analysis of target genes
In the category of cis regulation, a total of 5175 target genes of DE lncRNAs were predicted among the results of the DT_1d vs. C_1d, DT_5d vs. C_5d, and DT_10d vs. C_10d comparisons. GO analysis demonstrated that the target genes of the DE lncRNAs identified in DT_1d vs. C_1d were enriched in 42 GO terms, including nucleic acid metabolic process, RNA metabolic process, hydrolase activity and RNA polymerase II transcription factor binding transcription factor activity (Fig. 4). The results of KEGG pathway enrichment analysis showed that the target genes of the DE lncRNAs identified in DT_1d vs. C_1d were enriched in terms associated with metabolism (e.g., purine metabolism, nicotinate and nicotinamide metabolism, and fructose and mannose metabolism), RNA information processing (e.g., RNA polymerase, RNA degradation, RNA transport and mRNA surveillance pathway) and signaling pathways (e.g., mitogen-activated protein kinase signaling pathway, TGF-β signaling pathway and FoxO signaling pathway) (Fig. 5). The putative target genes of the DE lncRNAs were associated mainly with neurological system processes and immune effector processes in the comparison  of the 5-d and 10-d dinotefuran-treated groups (Additional file 5). In addition, the target genes were associated with RNA polymerase, purine metabolism, the MAPK signaling pathway, neuroactive ligand-receptor interaction, ECM-receptor interaction, and DNA replication in both the DT_5d vs. C_5d and DT_10d vs. C_10d comparisons (Additional file 6). Based on the trans regulation of lncRNAs [28], 245 target genes were predicted. We found that histone acetylation, protein acetylation and amino sugar metabolic processes were the most common significantly enriched terms in the three comparison groups (P < 0.05) (Additional file 7). Additionally, ribosome biogenesis in eukaryotes and DNA replication were the most enriched pathways in each of the comparisons (Additional file 8).

Subcellular localization of lncRNA TCONS_00086519
To explore the potential function of lncRNAs in the honey bee brain, the differentially and highly expressed lncRNA TCONS_00086519, which showed a full-length sequence of 471 bp (Fig. 7), was selected for the next study. Four protein-coding genes were found close to TCONS_0008615 through bioinformatic analysis (Additional file 10). The subcellular localization of lncRNAs is closely associated with their biological mechanism [30]. To further investigate the precise mechanism of the lncRNAs, an RNA-FISH assay was performed to detect the subcellular localization of TCONS_00086519. As shown in Fig. 8, TCONS_00086519 was mainly distributed in the cytoplasm in honey bee brain cells. This finding provided evidence that TCONS_00086519 may act as a competing endogenous RNA (ceRNA) or a molecular sponge to modulate the expression of its target mRNA or miRNA.

Quantitative real-PCR validation
To validate the RNA-seq results, four DE mRNAs (Mrjp1, LOC408790, LOC113219351 and Hbg3) and three DE lncRNAs (TCONS_00086519, TCONS_ 00024915 and TCONS_00023619) were subjected to qPCR. As shown in Fig. 9, Hbg3, TCONS_00086519, LOC408790 and LOC113219351 exhibited significantly higher expression in dinotefuran-treated bees than in control bees (P < 0.05), which was consistent with the RNA-seq data. The expression levels of TCONS_ 00024915 and TCONS_00023619 decreased significantly after persistent exposure to dinotefuran (P < 0.05). Moreover, Mrjp1 expression showed a 14.9-fold reduction relative to controls after dinotefuran treatment for 5 d.

Discussion
Studies on lncRNA expression in honey bees have been increasing in recent years, focused mainly on caste determination [19,31], the behavioral transition from nurse to forager [25], oviposition [32,33], and the response to pathogen infestation [20,21,34]. In the current study, RNA-seq was performed to identify the lncRNA and mRNA profiles of honey bees exposed to 0.01 mg/L dinotefuran for 1, 5, and 10 d to assess the potential regulators of honey bee growth and immune responses related to neonicotinoid insecticides. Finally, 6824 putative lncRNA transcripts and 19,837 known mRNA transcripts were obtained in the honey bee brain transcriptome. A total of 312 lncRNAs and 1341 mRNAs, 347 lncRNAs and 1458 mRNAs, 345 lncRNAs and 1155 mRNAs were identified as DE in the DT_1d vs. C_1d, DT_5d vs. C_5d, and DT_10d vs. C_10d comparisons, respectively. Among these transcripts, we found that several DE lncRNAs were actively associated with honey bee growth and immune responses after dinotefuran treatment. The lncRNA TCONS_00064387 exhibited lower expression in the dinotefuran-treated group than in the control group and targeted the Mrjp family genes in cis. Mrjp family genes have been validated and determined in honey bees, and encode major royal jelly proteins which are the major organic components of royal jelly [35]. The royal jelly master proteins Mrjp1-Mrjp8 have been repeatedly reported in the brain of honey bees [36][37][38] and play crucial neurobiological roles in the nervous system [39].The expression of Mrjp1 was detected in the mushroom body [36], optic lobe [38] and antennal lobe [40], and reduced expression of the Mrjp1 gene in the mushroom body decreased learning ability in worker honey bees [41]. The effect of sublethal concentrations of dinotefuran on reducing the expression of Mrjps suggested that the learning ability and nervous system of honey bees may be affected. The nAChRs are the targets of neonicotinoid insecticides [2], and five nAChR subunits (nAChRα4, the putative target of TCONS_00040953; nAChRα7, the target of XR_ 001705519.2; nAChRα8, the target of TCONS_ 00031755, TCONS_00031757 and TCONS_00031758; nAChRα9, the target of XR_003306308.1 and XR_ 003304779.1; and nAChRβ4, the target of XR_ 001705519.2) were separately investigated for their effects in trans. The induction of nAChRs can result in distinct effects; for instance, exposure to sublethal dinotefuran doses can affect olfaction, basic motor function and postural control in honey bees [12,42,43]. Defensin genes, which play an important role in the innate humoral immune system of honey bees and are associated with responses against G-and, G+ bacteria and fungi, showed significantly lower expression in the dinotefuran-treated group [44,45]. Defensin1 and Defensin 2 were identified as candidate targets of the annotated DE lncRNAs XR_003304788.1 and XR_ 003306226.1, respectively. This result is consistent with a previous study that showed that thiamethoxam, another neonicotinoid insecticide, was able to suppress defensin gene expression [46]. In general, our results may provide insights into the effect of dinotefuran on the immune system of A. mellifera from the perspective of lncRNAs. These findings suggest that DE lncRNAs may participate in the response to dinotefuran through cis and trans effects. Through GSEA with mRNA expression data sets and the target genes identified from DE lncRNA enrichment analysis, we found that the TGF-β signaling pathway was Fig. 8 The location of RNA TCONS_00086519 in the honey bee brain. (a) RNA fluorescence in situ hybridization (FISH) signal location using probes against TCONS_00086519; (b) the nucleus was stained with DAPI; (c) lncRNA TCONS_00086519 is distributed predominantly in the cytoplasmic region in the honey bee brain the most enriched pathway in both mRNAs and the target genes of lncRNAs. This result is in agreement with another neonicotinoid insecticide on adult honey bees that were exposed to realistic field doses of clothianidin [47]. A previous study showed that Varroa destructor activated the TGF-β signaling pathway to suppress wound healing and the immune response in honey bees [48]. Similarly, TGF-β signaling interferes with innate immunity in Drosophila larvae infected with symbiotic nematodes, and phenoloxidase activity in the hemolymph of infected larvae is regulated by activin signaling via TGFβ pathways [49]. This suggests that transforming growth factor-β (TGF-β) may play a critical role in the honey bee immune response to dinotefuran. In this study, mitogen-activated protein kinase 1 (LOC409523), which was enriched in the TGF-β signaling pathway, was indicated to be a potential target of DE lncRNA TCONS_ 00018497 and showed significantly lower expression in dinotefuran-treated bees. MAPK was associated with decreased developmental time and increased titer of juvenile hormone (JH) [50], a hormone that is essential for development and affects the structure and function of the adult honey bee nervous system [51,52]. The MAPK pathway has been reported to be involved in the antiviral immune responses of honey bees [53]. Moreover, pathway analyses demonstrated that a variety of mRNA expression data sets and target genes of DE lncRNAs in both comparison groups were enriched in material and energy metabolism-associated pathways, such as glycine, serine and threonine metabolism, and pentose and glucuronate interconversion. Metabolism plays an important role in the regulation of sustenance to innate and adaptive immune responses [54]. These results indicated that DE lncRNAs may act as regulators in the development, nerve conduction, and immune response of honey bees treated with dinotefuran.
Accumulating evidence indicates that lncRNAs participate in gene expression regulation by modulating chromosomal architecture [55,56], binding to transcription factors for direct recruitment [57], and RNA polymerase II [58] in the nucleus. However, several lncRNAs were found to control mRNA stability and modulate translation and posttranslational modifications in the cytoplasm [14,59,60]. Consequently, investigation into the subcellular localization of lncRNAs will be helpful to further illustrate the mechanisms and functions of lncRNAs [14,61]. Through RNA sequencing, we obtained the highly and significantly DE lncRNA TCONS_ 00086519, with a length of 471 bp, which potentially targets four genes in trans, including general odorantbinding protein 71 and smoothelin-like protein 1. The results indicated that TCONS_00086519 may be associated with odorant binding, development and immunity [62]. RNA FISH analysis revealed that TCONS_ 00086519 was mainly located in the cytoplasm, suggesting that it is capable of forming complexes with diverse structural and regulatory functions related to mRNA stability, mRNA translation, and signaling pathway modulation, serving as a competing endogenous RNA, or functioning as a precursor of microRNAs [60,63]. The results of this work provide a foundation for further study of the function and mechanism of lncRNA TCONS_00086519 in the honey bee brain.

Conclusion
This study revealed the expression patterns of lncRNAs in honey bee brains following exposure to a 0.01 mg/L dose of dinotefuran. Bioinformatic analysis showed that several lncRNAs and mRNAs were likely to participate in important biological processes associated with honey bee material and energy metabolism, neuroactive ligandreceptor interaction and odorant binding. Additionally, we found that the DE lncRNA TCONS_00086519 was distributed mainly in the cytoplasm and may act as a regulator by serving as a miRNA precursor or ceRNA involved in development and the immune response by regulating gene expression in trans in honey bees. Our results provide a foundation for understanding the lncRNA-mediated regulation of growth and immunity and contribute to a better understanding of honey beeinsecticide interactions and honey bees themselves.

Methods and materials
Honey bee rearing Frames containing sealed broods (near adult emergence) were collected from three healthy colonies located at the Zhejiang Academy of Agricultural Sciences (Hangzhou, China) and maintained in darkness in a climatecontrolled incubator (34°C ± 1°C, relative humidity [RH] 60 ± 10%). Then, we obtained newly emerged honey bees and placed them in cages (11 × 11 × 7 cm 3 , N = 60 in each cage).

Dinotefuran preparation and exposure
The residues of dinotefuran recorded in pollen, honey or syrup range from 0.03 ng/g to 147 ng/g [11,[64][65][66]. On this basis, a field-realistic level of dinotefuran of 0.01 mg/L was selected as the sublethal concentration in this study. A dinotefuran (> 99% purity, Sigma-Aldrich, St Louis, MO, USA) stock solution (1000 ng a.i./L) was prepared using a 50% sucrose:water solution as the solvent. The bees were treated with 2 mL of 0.01 mg/L dinotefuran solution in the experimental group (DT). Fifteen honey bees were collected on the 1st day (DT_1d), 5th day (DT_5d) and 10th day (DT_10d). The bees in the control group were fed 2 mL of a 50% sucrose: water solution, and samples were collected at the same time on the 1st day (C_1d), 5th day (C_5d) and 10th day (C_ 10d). The honey bees were fed ad libitum. Dead bees were removed and both treatment solution and untreated sucrose solution were replaced daily throughout the experiment. The bee samples were preserved in liquid nitrogen and stored at − 80°C until brain dissection. The in vivo portion of the study was carried out at 34°C ± 1°C under 60 ± 10% RH in darkness.

Library preparation and sequencing
Total RNA was isolated from 15 pooled brains per sample using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA was checked for purity and integrity by using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and a Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total of 3 μg RNA per sample was used to generate a complementary library with the NEB-Next® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations. After cluster generation using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, USA) according to the manufacturer's instructions, 18 libraries were sequenced on the Illumina HiSeq 2500 platform, and 150 bp paired-end reads were generated.

Identification of lncRNA and novel mRNA
Clean data were obtained after discarding the reads containing adapters, reads containing poly-Ns sequences (over 0.2%) and reads of low quality (over 50%) from the raw data. Q 20 and Q 30 values and GC contents were calculated to evaluate sequencing quality. The remaining clean data were aligned to the reference A. mellifera genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_ 003254395.2) with Hisat2 (https://ccb.jhu.edu/software/ hisat2/index.shtml) [67]. All mapped reads were assembled by StringTie [68] using the default parameters. The transcripts obtained from splicing were merged by using Cuffmerge, and transcripts of less than 200 nt in length with an uncertain strand orientation were removed. The obtained sequences were subjected to BlAST searches against the reference database, and known transcripts were filtered out by Cuffcompare [69]. The remaining transcripts were considered candidate lncRNAs, and those with coding potential were considered novel mRNAs.

Differential expression analysis
The expression levels of mRNAs and lncRNAs were calculated as fragments per kilobase of exon per million mapped fragments (FPKM) values using the Cuffdiff tool in Cufflinks v2.1.1 [26]. A P-value < 0.05 (Benjamini and Hochberg's false discovery rate) and a |log 2 (fold change) | > 1 were set as the thresholds for significantly differential expression.
ACGTCGAAGTTGGAGTACCCATGATCGACCGT TAG-3′ was designed and synthesized. RNA FISH was performed using an RNA FISH Kit (Gefanbio, Shanghai, China) according to the manufacturer's instructions. Briefly, honey bee brains were dissected and fixed with RNA-free paraformaldehyde. The paraffin sections were pretreated for hybridization with a 30% H 2 O 2 -methanol mixture (1:9) for 10 min at room temperature (RT). Then, the sections were treated via proteinase K digestion (20 min, 37°C) in 25% HCl, followed by a 1-min wash in 0.2% glycine irrigation solution. After fixation in 4% formaldehyde for 10 min, the sections were washed twice in acetic anhydride (pH = 8.0) at RT for 5 min each time. The sections were incubated in a humidified chamber (1 h, 65°C) covered with the prehybridization solution. The tissue sections were hybridized with the probe at 0.1 μM for 48 h at 65°C. The tissues were washed twice on coverslips with 2 × SSC (pH = 7.5). Finally, the sections were stained with DAPI, embedded in mounting agent, and then observed under a fluorescence microscope (IX-71, Olympus, Japan).

Target gene prediction
Differentially expressed lncRNAs were selected to predict potential cis and trans effects. The cis actions of lncRNAs affect neighboring target genes. Based on the positional relationship between lncRNAs and mRNAs, target genes were predicted in the regions 100 kb upand downstream of the DE lncRNAs. Based on the expression correlation between lncRNAs and mRNAs, the target genes were evaluated for a trans effect of lncRNAs with an absolute Pearson correlation ≥0.95.

Enrichment analysis
Gene Ontology (GO) enrichment analysis of the target genes of identified DE lncRNAs was performed using the GOseq R package [70], and Kyoto Encyclopedia of Genes and Genomes (KEGG, www.kegg.jp/kegg/kegg1. html) pathway analysis of the target genes was implemented by using KOBased Annotation System (KOBAS) v2.0 [71], considering a corrected p-value < 0.05 to indicate significant enrichment. GSEA software (http:// software.broadinstitute.org/gsea/index.jsp) was applied to interpret gene expression data [29,72]. Briefly, the genes were metrically ranked according to the signal-tonoise ratio ( μa−μb σaþσb ). Then, it was determined whether the genes in the gene set database (KEGG and GO databases) were ranked at the top or bottom of the list in terms of enrichment. The significance of the enrichment score (ES) was calculated with an empirical phenotypic-based permutation test. The normalized enrichment scores (NES) were used to compare the analysis results across gene sets, accounting for differences in gene set size and in the correlations between the gene sets and the expression dataset. The false discovery rate was calculated to control the false positive rate. IN GSEA, NES are determined as follows: NES ¼ actual ES meanðESs against all permutations of the dataset

Quantitative PCR analysis
First strand cDNA was obtained with the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). Quantitative real-time PCR was performed on an ABI 7500 Real-Time PCR system (Applied Biosystem, Carlsbad, CA, USA) using TB Green® Premix Ex Taq™ II (TaKaRa, Dalian, China) with specific primers ( Table 2). Relative gene expression levels were quantified based on β-actin expression by using the 2 -ΔΔCt method with three independent biological replicates. The differences in the relative expression levels of the mRNAs and lncRNAs identified between the dinotefuran-treated and control bees were analyzed using one-way ANOVA and an independent-samples t-test with SPSS 22.0 software (IBM, Armonk, NY, USA). p < 0.05 was considered statistically significant.