Effects of Dinotefuran on Brain miRNA Expression Profiles in Young Adult Honey Bees (Hymenopptera: Apidae).

Honey bees are important pollinators of wild plants and crops. MicroRNAs (miRNAs) are endogenous regulators of gene expression. In this study, we initially determined that the lethal concentration 50 (LC50) of dinotefuran was 0.773 mg/l. Then, the expression profiles and differentially expressed miRNAs (DE miRNAs) in honey bee brains after 1, 5, and 10 d of treatment with the lethal concentration 10 (LC10) of dinotefuran were explored via deep small-RNA sequencing and bioinformatics. In total, 2, 23, and 27 DE miRNAs were identified after persistent exposure to the LC10 of dinotefuran for 1, 5, and 10 d, respectively. Some abundant miRNAs, such as ame-miR-375-3p, ame-miR-281-5p, ame-miR-3786-3p, ame-miR-10-5p, and ame-miR-6037-3p, were extremely significantly differentially expressed. Enrichment analysis suggested that the candidate target genes of the DE miRNAs are involved in the regulation of biological processes, cellular processes, and behaviors. These results expand our understanding of the regulatory roles of miRNAs in honey bee Apis mellifera (Hymenopptera: Apidae) responses to neonicotinoid insecticides and facilitate further studies on the functions of miRNAs in honey bees.

Honey bees are economically important insects that provide humans with many products and are important pollinators of most common crops and wild plants worldwide (Ollerton et al. 2011). However, neonicotinoid insecticides have been clearly indicated as potential causes of declines in managed and wild bee populations (Lundin et al. 2015). Previous studies have shown that physiological and behavioral processes are significantly impaired after chronic exposure to neonicotinoid insecticides (Peng and Yang 2016, Li et al. 2017). In addition, these insecticides affect immunocompetence (Brandt et al. 2016) and have been found to promote the replication of a viral pathogen in honey bees (Di Prisco et al. 2013). Although three neonicotinoids, clothianidin, thiamethoxam, and imidacloprid, have been banned by the European Union, other neonicotinoid insecticides are still widely used for the protection of crops because of their low mammalian toxicity and broad-spectrum activity (Simondelso et al. 2015). With time and changes of the environment, insecticides are gradually decomposed or diluted in the field, resulting in sublethal effects on the behavior, physiology, biochemistry, and tissue of insects (Haynes 1988). Therefore, treatment with the LC1–LC50 of insecticides may be a better way to identify the physiological and immune responses involved (Desneux et al. 2007, Kim et al. 2018, Gao et al. 2020).

MicroRNAs (miRNAs) are a class of noncoding RNAs that are approximately 22 nt in length and exist in almost all multicellular organisms. Previous studies have shown that miRNAs regulate gene expression through translational inhibition or target gene degradation by binding to the 3′-untranslated regions (UTRs; Bartel 2009), 5′-UTRs (Zhou and Rigoutsos 2014), or Coding Sequences (CDSs) (Brümmer and Hausser 2014) of mRNA sequences. Many other studies have shown that miRNAs participate in developmental pathways (Zondag et al. 2012), caste determination (Ashby et al. 2016), memory regulation (Cristino et al. 2014), segmentation (Freitas et al. 2017), and immune defenses (Lourenço et al. 2013) in honey bees. Although research on neonicotinoid insecticides and honey bees has been conducted, there have been few studies on the effects of common neonicotinoid insecticides on honey bee miRNA expression (Shi et al. 2017).

In the present study, we determined the lethal concentration 50 (LC50) of dinotefuran, which is a third-generation neonicotinoid (subclass: furanicotinyl compounds), and investigated the effects of dinotefuran (LC10) on miRNA expression profiles in the brains of honey bees using RNA-seq. The objectives of this study were to provide novel insights into the effects of neonicotinoid insecticides on pollinators and to facilitate further studies on the functions of miRNAs in the honey bee brain.

Materials and Methods

Acute Oral Toxicity Assay

A dinotefuran (Sigma–Aldrich, St Louis, MO) stock solution (1,000 ng a.i./l) was prepared using 50% sucrose:water solution as the solvent. Six concentrations ranging from 1,000 to 0.01 ng a.i./l were prepared directly in a sucrose:water (1:1, w:v) diet after the pretest. In total, 30 newly emerged honey bees from each treatment group were divided into three cages of 10 honey bees, and a control group was set. The bees were fasted for 2 h, and then a plastic container filled with a diet (average of 10 µl diet/bee) and a cotton swab soaked in distilled water were placed at the cage bottom (Oliveira et al. 2014). After 48 h, the number of dead bees in each treatment group was recorded, and the LC10 and LC50 were determined by Data Processing System (DPS) software 18.01 (Tang and Zhang 2013).

Honey Bee Rearing and Dinotefuran Exposure

On the basis of the acute toxicity assay results, the lethal concentration 10 (LC10), 0.465 mg/l, was selected as the experimental dose of dinotefuran in this study. Honey bees with sealed brood (near-adult emergence) were collected from a healthy colony located at the Zhejiang Academy of Agricultural Sciences (Hangzhou, China) and kept in darkness in a climate-controlled incubator (34 ± 1°C, relative humidity [RH] 60 ± 10%). The newly emerged honey bees were put into six cages (11 × 11 × 7 cm3) with 60 bees in each cage, and maintained at 32 ± 1°C with 60 ± 10% RH in darkness. In the dinotefuran-treated group, the honey bees were treated with a 2-ml LC10 dose of dinotefuran sucrose solution in three replicates (three cages). In the control group, the honey bees were fed 2-ml 50% sucrose:water solution every day. The dead bees were removed, and the solution was replaced daily throughout the experiment. Fifteen honey bees from each cage were collected at 1, 5, and 10 d after oral exposure to dinotefuran, respectively. The bee samples were immediately frozen in liquid nitrogen and stored at −80°C until brain dissection.

RNA Extraction, Library Preparation, and Sequencing

The total RNA was isolated from fifteen pooled brains per sample using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s protocol. The RNA purity and quantity were analyzed with a NanoPhotometer spectrophotometer (Implen, Carlsbad, CA) and a Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA). The criteria used to select the RNA for the following analysis were an A260/A280 ≥ 1.8, an A260/A230 ≥ 2.0, and an RNA integrity number > 7.0. Small-RNA (sRNA) libraries were generated according to the instructions of an NEBNext Multiplex Small-RNA Library Prep Set for Illumina (Illumina, San Diego, CA). The purified library products were evaluated using an Agilent Bioanalyzer 2100 system and sequenced on the HiSeq 2500 platform (1 × 50 bp) after cluster generation.

Bioinformatic Analysis

The clean reads were obtained by removing the 5′ and 3′ adapter sequences, reads containing poly-N sequences, and low-quality reads (length < 18 nt) from the raw data. All clean reads were mapped to miRbase 20.0 (www.mirbase.org) to identify the known miRNAs. Novel miRNAs were predicted using the miREvo (Wen et al. 2012) and miRdeep2 (Friedländer et al. 2012) software programs through exploring the secondary structure, the Dicer cleavage site, and the minimum free energy of the sRNA tags unannotated in the former steps. The miRNA expression levels were calculated as transcripts per million reads (TPM) (Zhou et al. 2011) values, as follows: . Differential expression analysis was performed with the DEGseq R package (1.18.0). A |log2 (fold change) | value > 1 and a q value < 0.05 were set as the thresholds for significance. The target genes of the significantly differentially expressed miRNAs (DE miRNAs) were predicted by miRanda (Enright et al. 2003) and RNAhybrid (Krüger and Rehmsmeier 2006). Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the candidates target genes of the DE miRNAs were performed using the GOseq R package (Young et al. 2010) and the KEGG Orthology-Based Annotation System (KOBAS; Mao et al. 2005), respectively. Terms and pathways with a corrected P < 0.05 were considered significantly enriched.

Quantitative Real-Time Polymerase Chain Reaction Analysis

To validate the miRNA expression levels of sequencing, six DE miRNAs (ame-miR-6001-3p, ame-miR-279c-3p, ame-miR-306-5p, ame-miR-3049-3p, ame-miR-3049-5p, and ame-miR-3720-5p) were selected for quantitative real-time polymerase chain reaction (qRT-PCR) validation. RT-PCR was performed on an ABI 7500 Real-Time PCR System (Applied Biosystems, Carlsbad, CA) using a miRcute Plus miRNA qPCR Kit (SYBR Green; Tiangen, Beijing, China) with specific primers (Supp Table 1 [online only]). U6 was chosen as an endogenous control for miRNA expression (Mount et al. 2007). Three technical replicates and three biological replicates were used to analyze the effects of dinotefuran treatment on miRNAs. The 2−ΔΔCt method was used to calculate fold changes in miRNA expression (Livak and Schmittgen 2001). The differences in the relative expression levels of the miRNAs obtained between dinotefuran-treated and control bees were analyzed using an independent-sample t-test with SPSS 22.0 software (IBM, Armonk, NY).

Results

Acute Toxicity of Dinotefuran

Treatment of honey bees with serial dilutions of dinotefuran significantly affected motor coordination and survival. The bees showed irritability, initial overexcitation, and subsequent body curling, wing spreading, and mouthpart extension. Higher oral doses of dinotefuran affected motor coordination more quickly and more clearly compared with lower doses. The LC50 of dinotefuran for newly emerged honey bees was 0.773 mg/l (Fig. 1). The LC10 (0.465 mg/l diet) was calculated for use in subsequent experiments.

Fig. 1.

Toxicity regression curve of dinotefuran on Apis mellifera.

Overview of the sRNA Libraries

To identify the miRNA profiles of dinotefuran-treated and control bees of different times, 18 sRNA libraries were constructed. After low-quality reads and adapter sequences were filtered out, more than 92.49% of clean reads were mapped to the Apis mellifera L. reference genome. Comparisons with the NCBI and Rfam databases revealed that 58.28–71% of the mapped sRNA reads in each library were annotated to known miRNAs, whereas 7.3–14.59% were rRNAs, 0.24–0.85% were tRNAs, 0.03–0.10% were snRNAs, and 0.27–0.62% were snoRNAs (Supp Table 2 [online only]). The lengths of the sRNA reads were between 18 and 35 nt. A peak in frequency was observed at 22 nt for both dinotefuran-treated and control bees (Fig. 2), which accounted for 45.14% of the clean reads in both groups; the lengths with the next highest frequencies were 21 and 23 nt.

Fig. 2.

Distribution of the lengths of small RNA reads in dinotefuran-treated bees and control bees.

Expression of Known and Novel miRNAs

In total, we detected 179 known mature miRNAs (Fig. 3A), 42 novel mature miRNAs, 190 known hairpin miRNAs, and 44 novel hairpin miRNAs in the dinotefuran-treated and control bees (Supp Table 3 [online only]). In total, 172 miRNAs were common to all three sampled age groups of bees during normal growth (Fig. 3B). Among them, ame-miR-276-3p, ame-miR-7-5p, ame-miR-184-3p, ame-miR-1-3p, ame-miR-9a-5p, ame-bantam-3p, ame-miR-8-3p, ame-miR-2796-3p, ame-miR-3477-5p, and ame-miR-87-3p were the most abundant known miRNAs in both groups at three test times. (Supp Table 4 [online only]). The predicted novel miRNAs were all expressed at low levels (Supp Table 5 [online only]). Calculation of the log2 expression ratios with thresholds of a |log2 (fold change)| > 1 and a P < 0.05 revealed that there were two downregulated miRNAs at 1 d of age (Fig. 4A), 13 upregulated and 10 downregulated miRNAs at 5 d of age (Fig. 4B), and 13 upregulated and 14 downregulated miRNAs at 10 d of age (Fig. 4C) in dinotefuran-treated bees compared with control bees. Some of the abundant miRNAs, such as ame-miR-375-3p, ame-miR-281-5p, ame-miR-3786-3p, ame-miR-10-5p, ame-miR-6037-3p, and ame-miR-306-5p, were extremely significant DE miRNAs (P < 0.01). The significantly DE miRNAs also changed during continuous treatment with dinotefuran (Fig. 4D).

Fig. 3.

Distribution of known mature miRNAs. (A) Dinotefuran-treated group versus control group; (B) distribution of miRNA expression species of honey bees at 1, 5, and 10 d old under the normal growth period.

Fig. 4.

Volcano and Venn diagram of DE miRNAs between dinotefuran-treated bees and control bees. Red dots indicate upregulated miRNAs. Green dots indicate downregulated miRNAs. (A) DE miRNAs between dinotefuran-treated bees and control bees at 1 d; (B) DE miRNAs between dinotefuran-treated bees and control bees at 5 d; (C) DE miRNAs between dinotefuran-treated bees and control bees at 10 d; (D) DE miRNAs among dinotefuran-treated bees at three test times.

Target Prediction and Functional Enrichment Analysis

Based on the intersections of the results from miRanda and RNAhybrid, 2,699 target genes were predicted for 45 DE miRNAs among the 1-d dinotefuran-treated group (1dDT) versus 1-d control group (1dC), 5-d dinotefuran-treated group (5dDT) versus 5-d control group (5dC), and 10-d dinotefuran-treated group (10dDT) versus 10-d control group (10dC) comparisons. As shown in Supp Fig. 1 (online only), the putative target genes of the DE miRNAs between the 1-d dinotefuran-treated and control groups were associated mainly with biological regulation processes. The integral component of membrane and intrinsic component of membrane terms was the most enriched terms in the cellular component category. In addition, the biological regulation, regulation of biological process, and regulation of cellular process terms were the three most significantly enriched biological process terms in both the 5dDT versus 5dC and 10dDT versus 10dC comparisons (Supp Figs. 2 and 3 [online only]).

The putative target genes of the DE miRNAs in 1dDT versus 1dC were enriched in 16 pathways, and inositol phosphate metabolism, the Hippo signaling pathway, and ribosome biogenesis in eukaryotes were the most enriched pathways (Supp Fig. 4 [online only]). The target genes of the DE miRNAs in 5dDT versus 5dC were involved in 88 pathways, the most enriched of which were the mucin-type O-glycan biosynthesis, extracellular matrix–receptor interaction, and phototransduction-fly pathways (Supp Fig. 5 [online only]). In the 10dDT versus 10dC comparison, the Hippo signaling pathway-fly was the most enriched pathway, followed by the mucin-type O-glycan biosynthesis pathway, the Wnt signaling pathway, the endocytosis pathway, and the FoxO signaling pathway (Supp Fig. 6 [online only]).

qRT-PCR Validation of DE miRNAs

To validate the reliability of the RNA-seq results, six DE miRNAs were selected from the 5dDT versus 5dC and 10dDT versus 10dC comparisons and subjected to qRT-PCR with three biological replicates (Fig. 5A) All six selected miRNAs were significantly differentially expressed in dinotefuran-treated bees compared with control bees (P < 0.05), and the relative fold change in expression were consistent with the RNA sequenceing data (Fig. 5B). These results indicated that the RNA-seq data for the DE miRNAs were reliable.

Fig. 5.

Validation of the RNA-seq data by qRT–PCR. (A) The relative expression levels of six randomly selected miRNAs between dinotefuran-treated bees and control bees. (B) Fold change of the dinotefuran-treated bees versus the control bees was verified by RT–PCR. The expression data are presented as the mean ± SEM. *P < 0.05, **P < 0.01.

Discussion

Dinotefuran is a third-generation neonicotinoid insecticide with broad-spectrum and systemic insecticidal activity that causes neurotoxicity by strongly binding to nicotinic acetylcholine receptors (nAChRs) in the insect brain (Wakita et al. 2003, Goulson 2013). In this study, we determined that the LC50 of dinotefuran was 0.773 mg/l. This result is similar to that of a previous study reporting that dinotefuran is extremely toxic to adult honey bees (LC50 = 1.29 mg/l; Badawy et al. 2015). Therefore, dinotefuran should be used with caution to avoid very negative impacts on nontarget insects such as honey bees. The incomplete lethal effect of insecticides on living individuals may affect the biological characteristics of insects. It is meaningful to study the effects of incomplete lethal doses of insecticides on insect physiology and immune responses. Studies have shown that treatment with the LC25 of dinotefuran caused effects on survival, and developmental duration of Bemisia tabaci (Qu et al. 2017) and the differential gene expression profiles in Frankliniella occidentalis were affected when treated with LC10 of dinotefuran (Gao et al. 2020).

miRNAs play important roles in the genetic regulation of insect growth, development, and caste differentiation (Ashby et al. 2016, Ninova et al. 2016, Epstein et al. 2017). In this study, the miRNA expression profiles varied during growth and dinotefuran treatment (Fig. 3). This finding indicates that different types of miRNAs expressed at different levels might exert specific functional effects during growth and dinotefuran exposure. The 10 most abundant miRNAs (ame-miR-276-3p, ame-miR-7-5p, ame-miR-184-3p, ame-miR-1-3p, ame-miR-9a-5p, ame-bantam-3p, ame-miR-8-3p, ame-miR-2796-3p, ame-miR-3477-5p, and ame-miR-87-3p) were completely consistent between groups, although the rankings based on expression were different between the dinotefuran-treated and control groups. Among these miRNAs, ame-miR-9a was found involved in male development in bumblebees (Liu et al. 2019). In addition, a previous study has shown that ame-miR-1 affects host and fungal pathogen interactions by regulating the host proliferation, apoptosis, and immune response (Chen et al. 2019). Guo et al. (2016) has found that ame-miR-7 was higher expression in the abdomen of developing (pupal) workers than queens and might be related to the differential development of worker and queen bee larvae. Ame-miR-184 has been found to be associated with caste-independent ovarian activity in queen and worker honey bees (Macedo et al. 2016).

In this study, 2 DE miRNAs in 1dDT versus 1dC, 23 DE miRNAs in 5dDT versus C_5d, and 27 DE miRNAs in 10dDT versus 10dC were investigated. Among these DE miRNAs, 7 miRNAs (ame-miR-6001-3p, ame-miR-6041-3p, ame-miR-279c-3p, ame-miR-306-5p, ame-miR-3049-3p, ame-miR-3049-5p, and ame-miR-3720-5p) that were identified in both 5dDT versus 5dC and 10dDT versus 10dC were investigated. The expressions of ame-miR-276 and ame-miR-3477 in the dinotefuran-treated group were higher than those in the control group. Moreover, these miRNAs were highly expressed. It has been reported that ame-miR-276 is highly expressed in the optic lobes and in the small-type Kenyon cells of the mushroom bodies in honey bee brains, suggesting its involvement in neural function (Hori et al. 2011). These results indicate that these miRNAs might play important roles in mediating the effects of chronic dinotefuran treatment in honey bee brains and demonstrate that the expression levels of miRNAs are likely affected by dinotefuran treatment.

Pipsqueak, which was predicted to be regulated by ame-miR-375-3p, is an important developmental regulator in Drosophila that exhibits pleiotropic functions during oogenesis, embryonic pattern formation, and adult development (Lehmann et al. 1998). Another predicted target gene of ame-miR-375-3p, Ultraspiracle, influences honey bee social behavior and behavior-related gene expression by mediating responses to an endocrine regulator that is juvenile hormone (Ament et al. 2012). The CAT gene, which encodes a binding enzyme with iron porphyrin as a cofactor, is one of the key enzymes involved in the protection of organisms from oxidative stress and is considered an enzymatic biomarker of exposure to xenobiotics such as thiamethoxam in the honey bee, A. mellifera (Badiou-Bénéteau et al. 2012). Hbg3, which is the potential target gene of ame-miR-3477-5p, is an isoform of α-glucosidase (HBGase) that catalyzes the cleavage of α-glycosidic linkages of polysaccharides from the nonreducing ends and produces α-glucose in honey (Na Ayutthaya et al. 2018). In addition, broad-complex (Br-c), which was predicted to be a target of ame-miR-3049-5p, is a kind of ecdysteroid-regulated gene that is selectively expressed in the large-type Kenyon cells of the mushroom bodies in the worker bee brain and is involved in the regulation of brain function (Paul et al. 2006).

The potential target genes of the DE miRNAs were significant enriched in 19 pathways (P < 0.05), including Hippo signaling pathway, Wnt signaling pathway, mucin-type O-glycan biosynthesis, endocytosis, neuroactive ligand–receptor interaction, and so on. The Hippo signaling pathway plays a critical role in organ size control and tissue regeneration through cell proliferation and apoptosis (Halder and Johnson 2011, Zhao et al. 2011). Wnt signaling has an important role in the development and function of the nervous system involving the protection and subsequent nuclear import of β-catenin in insect (Twyman 2009). Endocytosis is responsible for the cellular immunity of honey bees (Strand 2008). These findings suggest that DE miRNAs might play a regulatory role in the development, nerve conduction, and immune defense of honey bees under the treatment of dinotefuran.

Conclusions

This study presented the first description of miRNA profiles in the honey bee brain affected by the LC10 of dinotefuran. In total, 179 known mature miRNAs were identified and 42 novel mature miRNAs were detected. Among them, 2, 23, and 27 miRNAs were significantly differentially expressed on the first day, the fifth day, and the tenth day of dinotefuran treatment, respectively. DE miRNAs may play important roles in development, immune response, and behaviors through the potential target genes. Our findings contribute to a better understanding of miRNA-mediated regulation and provide a theoretical basis for the detection of the molecular mechanisms of dinotefuran in honey bees.

Supplementary Data

Supplementary data are available at Journal of Insect Science online.

Fig. S1. GO enrichment analysis of potential target genes of DE miRNAs in 1dDT versus 1dC.

Fig. S2. GO enrichment analysis of potential target genes of DE miRNAs in 5dDT versus 5dC.

Fig. S3. GO enrichment analysis of potential target genes of DE miRNAs in 10dDT vs 10dC.

Fig. S4. KEGG pathway enrichment analysis of potential target genes of DE miRNAs in 1dDT vs 1dC.

Fig. S5. KEGG pathway enrichment analysis of potential target genes of DE miRNAs in 5dDT vs 5dC.

Fig. S6. KEGG pathway enrichment analysis of potential target genes of DE miRNAs in 10dDT vs 10dC.

Table S1. Primers used for the reverse transcription quantitative PCR (RT-qPCR) analysis.

Table S2. Summary of ncRNA annotation for small RNAs.

Table S3. Identification of known miRNAs and prediction of novel miRNAs.

Table S4. The most abundant known miRNAs in the dinotefuran-treated and control honey bees at three test times.

Table S5. The most abundant predicted miRNAs in the dinotefuran-treated and control honey bees at three test times.

Click here for additional data file.

Click here for additional data file.