Influence of the Neonicotinoid Insecticide Thiamethoxam on miRNA Expression in the Honey Bee (Hymenoptera: Apidae).

MicroRNAs (miRNAs) are small endogenous noncoding single-stranded RNAs regulating gene expression in eukaryotes. They play important roles in regulating caste differentiation, behavior development, and immune defences in the honey bee, Apis mellifera (Linnaeus) (Hymenoptera: Apidae). In this study, we explored the effect of the neonicotinoid insecticide, thiamethoxam, on miRNA expression in this species using deep small RNA sequencing. The results showed that seven miRNAs were significantly differentially expressed (q-value <0.01 and |log2(fold-change)| >1) upon exposure to 10 ppb thiamethoxam over 10 d. Some candidate target genes were related to behavior, immunity, and neural function. Several miRNAs, including ame-miR-124, ame-miR-981, ame-miR-3791, and ame-miR-6038, were selected and further validated using real-time quantitative PCR analysis. The findings expand our understanding of the effects of neonicotinoid insecticides on honey bees at the molecular level.

Neonicotinoid insecticides are widely used in agriculture for controlling insect pests (Simon-Delso et al. 2015). However, these insecticides also simultaneously affect nontarget insect species including important pollinators such as the honey bee, Apis mellifera. Neonicotinoids are systemic insecticides that mainly act on the insect nervous system through an agonistic action on nicotinic acetylcholine receptors (nAChRs) (Brown et al. 2006). Sublethal doses of neonicotinoids can impair the brain and mushroom body of honey bees (Oliveira et al. 2014, Peng and Yang 2016) and disrupt learning and memory capacity (Aliouane et al. 2009, Mengoni and Farina 2015, Alkassab and Kirchner 2016). In addition, these insecticides also negatively affect secondary targets such as the midgut and Malpighian tubule (Catae et al. 2014, Oliveira et al. 2014) and decrease immunocompetence (Brandt et al. 2016). Recently, the sublethal effects of neonicotinoids on mRNA transcripts have been investigated in honey bees (Christen et al. 2016, 2017), but the impact on miRNAs is poorly understood.

MiRNAs are a class of 18−24 nucleotide (nt), endogenously initiated small noncoding single-stranded RNAs that regulate target gene expression via either translational repression or mRNA degradation by combining with target mRNAs in a complementary manner (Bartel 2004, Pillai et al. 2005). To date, more than 250 miRNAs in the honey bee have been identified and registered (http://www.mirbase.org/index.shtml). A series of studies on honey bee miRNAs showed that they play important roles in regulating caste differentiation (Guo et al. 2013, 2016), behavior development (Behura and Whitfield 2010, Greenberg et al. 2012, Liu et al. 2012, Nunes et al. 2013), and immune defences (Lourenço et al. 2013). Although the effects of neonicotinoids on miRNAs have been largely overlooked, Lourençoet et al. (2013) found that bacterial infection can dysregulate miRNA expression in honey bees.

In this study, we investigated the effects of the neonicotinoid insecticide thiamethoxam on expression of miRNAs in A. mellifera. In total, seven known miRNAs were significantly differentially expressed upon exposure to thiamethoxam. We also predicted the mRNA targets to investigate the functions of the differentially expressed miRNAs.

Materials and Methods

Rearing of Honey Bees

Two frames with sealed broods (near of adult emergence) were taken from a healthy colony located in the Institute of Apiculture Research of Anhui Agriculture University (Hefei, China) with no pesticides use to a conditioned incubator (35 ± 1°C, relative humidity [RH] of 50 ± 10% and in darkness) from July to September 2016. Then, we obtained newly emerged honey bees in 12 h and put them into cages (11 × 10 × 8 cm3). They were fed with enough bee bread (stored at −20°C) collected from the colonies of the same apiary and sucrose–water solution (1:1 w:v), and maintained at controlled incubator (30 ± 1°C, RH of 70 ± 10% and in darkness) for 3 d.

Thiamethoxam Prepared and Exposure

The residues of thiamethoxam are range from 0.6 to 53.3 ppb and 2.5 to 17.2 ppb in trapped pollen (Krupke et al. 2012, Mullin et al. 2010, Stoner and Eitzer 2012, Codling et al. 2016, Sánchez-Hernández et al. 2016) and honey (Codling et al. 2016, Sánchez-Hernández et al. 2016), respectively. On the basis of this, field-realistic level of thiamethoxam, 10 ppb, was selected as the sublethal concentration in this study. A thiamethoxam (>99% purity, J&K, Shanghai, China) stock solution (1,000 ppm) was prepared using acetone as a solvent. About 10 ppb thiamethoxam was prepared in the 50% sucrose–water solution (1:1 w:v) with the final concentration of acetone was 0.03%, 50% sucrose–water solution added the same acetone as control check. Four-d-old bees were used for subchronic exposure bioassays. The experiment was carried out with 60 bees per replicate and three replicates per treatment. Honey bees in each treatment received enough contaminated 50% sucrose–water solution and bee bread every day, and maintained at controlled incubator (30 ± 1°C, RH of 70 ± 10% and in darkness). After 10 d, we collected all remaining bees and stored them at −80°C. Fifteen bees from each replicate were pooled to obtain one RNA sample.

RNA Extraction

Pools of 15 bees from each sample were prepared for total RNA extraction using the Trizol protocol (Invitrogen, Carlsbad, CA). Qubit 2.0 Fluorimeter (Life Technologies, Carlsbad, CA) was used to measure RNA concentration with Qubit RNA Assay Kit (Invitrogen). The RNA Nano 6000 Assay Kit (Invitrogen) was used to assess RNA integrity in Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA).

Small RNA Library Construction and Sequencing

Three micrograms of total RNA from one RNA sample of each treatment were used as input material for the small RNA library. Sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Set for Illumina (Illumina, NEB, San Diego, CA) following manufacturer’s recommendations. An Illumina Hiseq 2500/2000 platform (Illumina) was used to sequence sequencing libraries and 50 bp two single-end reads were generated.

Data Analysis and sRNA Annotation

First, clean reads were obtained by removing reads containing insert tags, ploy A reads, low-quality reads, and oversized insertion from raw reads. Then, all the clean reads were classified into different categories according to the sequences with the noncoding RNAs (rRNA, tRNA, scRNA, snRNA, and snoRNA) in NCBI (http://www.ncbi.nlm.nih.gov/). MiRBase 20.0 was used to confirm the known miRNAs and novel miRNAs were predicted using the miREvo and miRdeep2 software.

Differential Expression of Known miRNA

First, we used transcript per million to calculate the miRNA expression levels. Then, the DEGseqR package was used to perform differential expression analysis of two treatments (thiamethoxam-treated bees and control bees). And, significantly differential expressed miRNAs between two treatments were defined as those with an adjusted P-value <0.01 (q <0.01) and |log2(fold-change)| >1 by default. We predicted the target genes of the significantly differential expressed miRNAs by miRanda (Enright et al. 2003).

Real-Time Quantitative PCR Analysis

To validate the data of deep sequencing, we randomly selected four significantly differential expressed miRNAs (ame-miR-124, ame-miR-981, ame-miR-6038, and ame-miR-3791) for qPCR assays. We started with 1 μg total RNA for each sample, and U6 was used as a reference gene. All the primers used were given in Table 1. SYBR PrimeScript miRNA RT-PCR Kit (Takara, Dalian, China) was used to obtain the cDNA and perform RT-qPCR according to the manufacturer’s protocol. Reverse transcription and RT-qPCR were performed using PCR system (Bioer Technology, Honshu, Japan) and StepOnePlus Real-Time PCR System (Life Technologies), respectively. The qPCR data were analyzed by 2−△△CT method, and the differences in the four miRNAs between thiamethoxam-treated and control bees were significant by independent-samples t-test.

Table 1.

Primer sequences

Genes	Forward primer (5′–3′)
ame-miR-124	TAAGGCACGCGGTGAATG
ame-miR-981	CGTTGTCAACGAAACCTG
ame-miR-3791	CACCGGGTAGGATTCATC
ame-miR-6038	GTATGTTTCTGTCTTATTTCATT
U6	TGCGGGTGCTCGCTTCGGCAGC

Results

Sequencing and Analysis of Small RNAs

Using the Illumina Hiseq 2500/2000 platform, we sequenced two small RNA libraries from thiamethoxam-treated and control honey bees. After removing low-quality reads, adaptor reads, and reads <18 nt, 9,960,504 and 9,149,499 small RNA reads of 18–35 nt remained for thiamethoxam-treated bees and control bees, respectively, which were used for subsequent analysis (Table 2). The length distribution of the small RNAs for both thiamethoxam-treated and control bees peaked at 21 nt, with 18.46% and 20.54% of total clean reads, respectively. The number of small RNAs of 18–20 nt and 23–28 nt was greater in thiamethoxam-treated bees than in control bees, whereas the number of small RNAs of other lengths was greater in control bees than in thiamethoxam-treated bees (Fig. 1).

Table 2.

Summary of ncRNA annotation for small RNAs

Type	Thiamethoxam-treated bees		Control bees
	Total sRNAs (%)	Unique sRNAs (%)	Total sRNAs (%)	Unique sRNAs (%)
All reads	9,960,504	759,970	9,149,499	509,663
Mapped reads	1,877,880 (18.85)	54,276 (7.14)	2,471,155 (27.01)	58,967 (11.57)
Known miRNA	1,667,138 (88.78)	1713 (3.16)	2,204,070 (89.19)	1721 (2.92)
rRNA	3777 (0.20)	1014 (1.87)	4596 (0.19)	963 (1.63)
tRNA	19,970 (1.06)	1157 (2.13)	16,411 (0.66)	798 (1.35)
snRNA	1279 (0.07)	520 (0.96)	1357 (0.05)	507 (0.86)
snoRNA	1231 (0.07)	88 (0.16)	2421 (0.10)	112 (0.19)
Unannotated reads	184,485 (9.82)	49,784 (91.72)	242,300 (9.81)	54,866 (93.05)

Fig. 1.

Distribution of the length of small RNA reads in thiamethoxam-treated bees and control bees.

All clean reads were mapped to the honey bee genome (NCBI: Amel_4.5), resulting in 1,877,880 (18.85%) and 2,471,155 (27.01%) mapped reads in thiamethoxam-treated and control bees, respectively (Table 2). We also divided the mapped reads into different categories. In each library, 88.78–89.19% of small RNA reads belonged to annotated miRNAs, 1.06–1.35% were tRNAs, 0.2–1.63% were rRNAs, 0.07–0.86% were snRNAs, and 0.07–0.19% were snoRNAs. The remaining 9.81–9.82% did not belong to any annotated transcribed reads (Table 2).

Expression of Known and Novel miRNAs

Our experiments detected 164 known mature miRNAs, of which 155 were present in thiamethoxam-treated bees and 148 were present in control bees (Table 3). In all, the expression levels of these miRNAs were to some extent different in thiamethoxam-treated bees compared with control bees (Supp Table 1 [online only]), of which seven were significantly differentially expressed based on a q-value <0.01 and |log2(fold-change)| >1 (Table 4). Among the significantly differentially expressed miRNAs, ame-miR-124, ame-miR-971, ame-miR-981, and ame-miR-6038 were 1.3445-, 2.6075-, 1.245-, and 1.0726-fold more highly expressed in thiamethoxam-treated bees than in control bees. However, ame-miR-279b, ame-miR-3791, and ame-miR-6051 were downregulated in thiamethoxam-treated bees by 1.5921-, 1.38-, and 1.9082-fold compared with control bees.

Table 3.

Mapping statistics of known miRNAs

Types	Total	Thiamethoxam- treated bees	Control bees
Known miRNAs	164	155	148
Mapped hairpin	176	166	160
Mapped unique sRNA	3434	1713	1721
Mapped total sRNA	3,871,208	1,667,138	2,204,070

Table 4.

MiRNAs significantly differentially expressed in thiamethoxam-treated bees compared with control bees

miR name	Fold-change (log2 thiamethoxam- treated bees/control bees)	q-value
ame-miR-124	1.3445	0.0028049
ame-miR-971	2.6075	0.0092951
ame-miR-981	1.245	6.30E−07
ame-miR-6038	1.0726	0.000243
ame-miR-279b	−1.5921	9.54E−05
ame-miR-3791	−1.38	2.32E−05
ame-miR-6051	−1.9082	0.000474

In total, 14 novel miRNA candidates were predicted using miREvo and miRdeep2 software, of which two only belonged to thiamethoxam-treated bees and two only belonged to control bees while the other 10 novel miRNA were expressed in both groups (Suppl Table 2 and 3 [online only]).

Target Prediction

Using the miRanda software (Enright et al. 2003), we predicted the target genes of the seven significantly differentially expressed miRNAs based on the 3′-untranslated region sequences of honey bee genes. The results showed that each miRNA can target several or even hundreds of genes (Supp Table 4 [online only]). For example, ame-miR-6038 targets 435 genes, whereas ame-miR-3791 targets just 55 genes.

As shown in Supp Table 5 (online only), all seven miRNAs except ame-miR-124 were predicted to target the transcription factor mblk-1-like (Mblk-1), whereas 18-wheeler (18w) was only regulated by ame-miR-124 and ame-miR-3791. The ecdysone receptor (Ecr) appeared to be modulated by ame-miR-3791, ame-miR-124, and ame-miR-981, and the nAChRα2 subunit, dopamine receptor 1 (Dop1), and cyclic guanosine monophosphate (cGMP)-dependent protein kinase foraging (For) were all targeted by ame-miR-6038 (Supp Table 5). Finally, the 5-HT1 serotonin receptor, worker-enriched antennal transcript (Wat), and the nAChRα3 subunit were predicted targets of ame-miR-981, ame-miR-971, and ame-miR-6051, respectively (Supp Table 5).

To confirm the results of deep sequencing, four significantly differentially expressed miRNAs were selected for qPCR validation. As shown in Fig. 2, all four selected genes were significantly differentially expressed in thiamethoxam-treated compared with control bees (ame-miR-124: t = 5.108, df = 4, P = 0.007; ame-miR-981: t = 3.875, df = 4, P = 0.018; ame-miR-6038: t = 2.320, df = 4, P = 0.029; and ame-miR-3791: t = −3.909, df = 4, P = 0.017).

Fig. 2.

Quantitative real-time PCR analysis of ame-miR-124, ame-miR-981, ame-miR-6038 and ame-miR-3791 in thiamethoxam-treated bees and control bees. Data are means ± SEM. * indicates that miRNAs were significantly differentially expressed in thiamethoxam-treated bees compared with control bees (Independent samples t-test, p<0.05). The expression levels of all miRNAs were scaled relative to that of ame-miR-124 in thiamethoxam-treated bees, which was set as 1.

Discussion

In this study, deep sequencing of small RNAs and bioinformatic analysis resulted in the identification of seven miRNAs that are significantly differentially expressed in honey bees upon exposure to 10 ppb thiamethoxam. Four selected miRNAs (ame-miR-124, ame-miR-981, ame-miR-3791, and ame-miR-6038) were further investigated by qPCR, and the results were consistent with those of deep sequencing. The functions of some miRNAs have been studied, and miRNA-124 is required for neural development (Wang et al. 2014) and control the phase of circadian locomotor behavior in Drosophila (Zhang et al. 2016).

We also predicted the target genes to further investigate the functions of the seven significantly differentially expressed miRNAs. nAChRs are targets of neonicotinoid insecticides and two nAChRs subunits, nAChRα2 and nAChRα3, were identified as candidate target genes of ame-miR-6038 and ame-miR-6051, respectively. Induction of nAChRs can result in obvious effects; for example, exposure of honey bees to sublethal thiamethoxam doses can affect memory and locomotor capacity (Aliouane et al. 2009, Charreton et al. 2015). The Toll-related receptor gene 18w is potentially involved in antimicrobial immune defences in honey bees (Aronstein and Saldivar 2005), and this gene was a predicted target of ame-miR-124 and ame-miR-3791.

The ame-miR-6038 miRNA had the largest number of candidate target genes (435), including nAChRα2, For, and Dop1. The For gene, which encodes a cGMP-dependent protein kinase (PKG), has been associated with the development of foraging behavior in honey bees (Ben-Shahar 2005). Activation of PKG was shown to induce nurse bees to forage precociously (Ben-Shahar et al. 2003), and increasing PKG activity pharmacologically increased sucrose responsiveness in nurse bees to the level of foragers (Thamm and Scheiner 2014). Biogenic amines, including dopamine, serotonin, octopamine, and tyramine, are candidate regulators of foraging behavior (reviewed in Page et al. 2006). Some appear to modulate sensory sensitivity; for example, dopamine reduces gustatory responsiveness (Scheiner et al. 2002) and serotonin decreases the direction-specific visual antennal reflex (Erber et al. 1993a, 1993b). Interestingly, the 5-HT1 serotonin receptor was a predicted target gene of ame-miR-981. In addition, the wat gene, which was predicted to be modulated by ame-miR-971, is involved in preferential chemosensory signaling in the honey bee (Kamikouchi et al. 2004). Expression of ame-miR-6038, ame-miR-981, and ame-miR-971 was upregulated in thiamethoxam-challenged bees (Table 3), suggesting their target genes might be downregulated. This indicates that the neonicotinoid insecticide thiamethoxam might have negative effects on behavior development and sensory sensitivity in honey bees.

Conclusions

This study investigated miRNA expression patterns in honey bees upon exposure to 10 ppb thiamethoxam over 10 d. Seven miRNAs (ame-miR-124, ame-miR-981, ame-miR-3791, ame-miR-6051, ame-mir-6038, ame-miR-279b, ame-miR-3791, and ame-miR-6051) were significantly differentially expressed following thiamethoxam challenge. The target genes of these miRNAs were predicted to probe their functions. The results provide valuable information for exploring the molecular mechanisms of neonicotinoid insecticides on honey bees.

Supplementary Data

Supplementary data is available at Journal of Insect Science online.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.