Missing Nurse Bees—Early Transcriptomic Switch From Nurse Bee to Forager Induced by Sublethal Imidacloprid

The environmental residue/sublethal doses of neonicotinoid insecticides are believed to generate a negative impact on pollinators, including honey bees. Here we report our recent investigation on how imidacloprid, one of the major neonicotinoids, affects worker bees by profiling the transcriptomes of various ages of bees exposed to different doses of imidacloprid during the larval stage. The results show that imidacloprid treatments during the larval stage severely altered the gene expression profiles and may induce precocious foraging. Differential expression of foraging regulators was found in 14-day-old treated adults. A high transcriptome similarity between larvae-treated 14-day-old adults and 20-day-old controls was also observed, and the similarity was positively correlated with the dose of imidacloprid. One parts per billion (ppb) of imidacloprid was sufficient to generate a long-term impact on the bee’s gene expression as severe as with 50 ppb imidacloprid. The disappearance of nurse bees may be driven not only by the hive member constitution but also by the neonicotinoid-induced precocious foraging behavior.


INTRODUCTION
With high pollinating efficiency, workers of the honey bee, Apis mellifera, are the most important pollinators for agricultural crops as well as wild plants (Calderone, 2012;van der Zee et al., 2012;Garibaldi et al., 2013Garibaldi et al., , 2016Hung et al., 2018). Although honey bees may be exposed to different types of agricultural chemicals, previous studies evaluating the impacts of pesticides found that a group of insecticides called neonicotinoids may be particularly hazardous to honey bees (Goulson, 2013;Lu et al., 2014;Sandrock et al., 2014;Goñalongs and Farina, 2015). While foragers are more frequently in contact with the pesticides, the whole colony is also under the threat of sublethal neonicotinoids as neonicotinoid-contaminated nectar and pollen are delivered to the colony through foraging activities (Chauzat et al., 2006;Škerl et al., 2009;Codling et al., 2016;Mitchell et al., 2017;Böhme et al., 2018).
A sublethal dose of imidacloprid, one of the major neonicotinoids mostly used around the world, may not cause acute bee death, but it can generate chronic irreversible damage to the bee colony. Chronic exposure to sublethal doses of imidacloprid can generate abnormal behavior or physiological disorders in bees (Goulson, 2013;Lu et al., 2014;Sandrock et al., 2014;Goñalongs and Farina, 2015). A long-term investigation revealed that honey bee colony health and overwintering success are negatively affected by exposure to field doses [20-100 parts per billion (ppb)] of imidacloprid and that colony survival is correlated with the amount of imidacloprid (Dively et al., 2015;Woodcock et al., 2017;Wood et al., 2018).
Larvae and nurse bees are affected by consuming neonicotinoid-contaminated pollen or nectar. Honey bee larvae can tolerate higher doses of imidacloprid than nurse bees [LD 50 = 4.17 µg, LC 50 = 138.84 parts per million (ppm)] (Dai et al., 2017), but larval synapse development is affected by low doses. After exposure to 10 ppb imidacloprid, honey bee larvae pupation and eclosion rates are not affected. However, the developed adults show a reduced synaptic density in the brain region responsible for olfactory and visual function, and the olfaction-related learning behavior is permanently impaired (Yang et al., 2012;Peng and Yang, 2016). With continuous exposure, all these physiological and behavioral disorders induced by the sublethal dose of neonicotinoids can be sustained throughout the lifespan of workers and exhibit in different syndromes with the progress of bee development.
Due to difficulties in observation, only a few studies report the impact of imidacloprid on nurse bee behavior; they report that nurse bees show less activity and social interaction after being exposed to a low dose of imidacloprid (Medrzycki et al., 2003;Forfert and Moritz, 2017). From the molecular perspective, genes with functions related to metabolism, biosynthesis, and royal jelly synthesis are affected after exposure to sublethal doses of imidacloprid or other neonicotinoids (Aufauvre et al., 2014;Shi et al., 2017a,b;Wu et al., 2017;Christen et al., 2018;Wang et al., 2018). Most molecular studies suggest an immediate response after treatment; nevertheless, it is unclear whether the impact is long-lasting. In this report, we comprehensively evaluate the short-and long-term effects of sublethal doses of imidacloprid on honey bees of different ages following exposure during the larval stage. We used next-generation sequencing (NGS) to investigate the transcriptome of 9-day-old larvae and 0-, 7-, 14-, and 20-day-old adults. Different concentrations of imidacloprid solutions, including 1, 10, and 50 ppb, were applied to confirm the dosage effect and determine the lowest effective concentration of imidacloprid that could alter honey bee gene expression.

Honey Bees
Honey bees were purchased from a local bee keeper at Hsinchu area and maintained at National Taiwan University and Miaoli District Agricultural Research and Extension Station (Council of Agriculture, Executive Yuan). Colonies contained at least eight honeycomb frames, a normal egg-laying queen, larvae, and pupae, as well as pollen and honey. Bee stability was continuously monitored. Pollen cake (45% sucrose, 30% pollens, 15% soybean powder, and 10% H 2 O) and 50% sucrose solution were provided during the non-flowering season. The area nearby the colonies was confirmed as not using neonicotinoid pesticides.

Imidacloprid Solution Preparation
Imidacloprid powder (Sigma) was dissolved in 100% DMSO to a final concentration of 100 ppm and then diluted to 1, 10, and 50 ppm with DMSO. All solutions were aliquoted in 1.5-ml tubes and stored at −20 • C until further use. Imidacloprid working solutions were freshly prepared every day before feeding. To prepare working solutions, 1, 10, and 50 ppm imidacloprid stock solutions were diluted with ddH 2 O to final concentrations of 1, 10, and 50 ppb (final concentration of DMSO is 0.1%).

Imidacloprid Feeding Assay
The imidacloprid feeding assay and sample collection were performed during October to December 2016. For sample collection, two colonies were selected. For each colony, at least 300 larvae were used. The queen was restricted in one frame to lay eggs for 16 h to synchronize the age of experimental individuals. After the larvae hatched (defined as 1-day-old larvae), 1 µl of 1, 10, or 50 ppb of imidacloprid solution was added to each cell. This feeding procedure was performed daily for 4 days, from 1-to 4-day-old larvae, to mimic the number of days honey bee larvae consume bee bread. For each experimental group, a total of 4 µl of imidacloprid solution was applied into every cell, resulting in final amounts of imidacloprid consumed per larva of 0.004, 0.04, and 0.2 ng, respectively. Parallel controls were larvae fed with either ddH 2 O or 0.1% DMSO (see Supplementary Figure 1A for the experiment procedure). At 10 days after the cells were capped, pupae were pulled out and kept separately according to the treatments in 96-well plates at 37 • C until eclosion. Bees from two different beehives were pooled as one group for the labeling and sampling procedure. Newly emerged adults were labeled on the thorax with different colors of acrylic paint based on the treatment and then randomly released into two beehives. Bees from different colonies were mixed and randomly assigned for the following experiment. The results represented in this report were the integrated effects from two hives of bee. Sampling was performed on 9-day-old larvae and 0-(newly emerged), 7-, 14-, and 20-day-old adults. The 0-, 7-, 14-, and 20-day-old adults conduct the general tasks of cell cleaning, brood rearing, pollen/nectar processing, and foraging, respectively (Seeley, 1982). For each age of worker bee, three biological replicates were collected, five individuals per replicate to final 15 individuals per age per treatment. Samples were stored in RNAlater (Invitrogen, Carlsbad, CA, United States) at −20 • C for subsequent RNA extraction.

RNA Extraction and RNA-seq Library Preparation
RNA from whole larvae and adults' heads was extracted using an RNeasy mini kit (Qiagen, Hilden, Germany) following the user manual. RNA quality was checked in an agarose gel, with no significant degradation of 28S and 18S bands. RNA quantity was measured using the NanoDrop 2000 UV-Vis spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). For each sample, at least 10 µg of total RNA was used for RNA-seq library construction. The library construction was performed according to Zhong et al. (2011). Briefly, polyA RNA enrichment was performed using Oligo d(T)25 Magnetic Beads (NEB, Ipswich, MA, United States), then eluted and fragmented simultaneously in 2 × RT buffer in the presence of random hexamers (final concentration 6 µM) and d(T) 23 VN (final concentration 5 µM). First-strand cDNA was synthesized using a ProtoScript II First Strand cDNA synthesis kit (NEB). Second-strand cDNA was synthesized using DNA polymerase I (NEB) and RNase H (NEB) with the presence of dUTP mix (dUTP, dATP, dCTP, dGTP; final concentration 1 mM). After end repair and dA tailing, the product was ligated with an Illumina Y-shaped adapter. The final product was treated with UDG (NEB) then PCR amplified with an index primer set. Library size and quality were checked using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, United States). Four lanes of 150-bp pairedend sequencing were performed using a NovaSeq 6000 sequencer (Illumina, San Diego, CA, United States).

Gene Expression, Differentially Expressed Genes, and Gene Ontology Analysis
FastQC was performed to check the quality of the reads, and seqtk 1 was used to remove the reads that failed to achieve the quality score threshold (30). RNA-Seq reads were then aligned to the adaptor, tRNA, and ribosomal RNA reference sequences using Bowtie2 (v.2.3.4.1;doi:10.1038/nmeth.1923) with default parameters. The remaining filtered reads were then mapped to the honey bee genome (GCA_000002195.1, Amel_4.5 assembly 2 ) using HISAT2 (version 2.1.0) (Kim et al., 2015) with parameters [-max-intronlen 50000 -p 6 -dta-cufflinks]. Analysis of the differentially expressed genes (DEGs) was performed based on the number of mapped raw reads and identified using Bioconductor DEseq2 (v.1.12.4; Bioconductor 3.7) with default parameter and Wald test (Love et al., 2014). The DEGs with FDR (false discovery rate) < 0.05 and fold change ≥ 2 (log 2 fold change ≥ 1 or ≤−1) were considered differentially expressed. The functional analysis and GO clustering of remaining DEGs were performed using the online gene functional classification tool, DAVID 3 (Huang et al., 2009a;Huang et al., 2009b). Three tables were generated for each subset of DEGs, including "functional annotation clustering, " "functional annotation chart, " and "functional annotation table."

Gene Expression Estimation and Differentially Expressed Gene Identification
To monitor the effect of sublethal doses of imidacloprid on honey bee gene expression, we collected larvae and adults of worker bee at five different ages. The sampling time was determined based on the relative probability of task performance as described in Seeley (1982). For each age, a total of three biological replicates were collected, with five bees per replicate to a final 15 bees per age per treatment. Honey bees were collected from two beehives and then randomly assigned for each replicate; thus, the results represented in this report may eliminate the colony effects. Four lanes of 150 bp Illumina paired-end sequencing were generated; read yields per sample are shown in Supplementary Table 1A.
Read counts for each gene, as well as FPKM levels at various development stages, are listed in Supplementary Tables 1B-F, 2A-E, respectively. Imidacloprid was supplied in a feeding solution using 0.1% DMSO in ddH 2 O as a solvent; control bees were fed with 0.1% DMSO solution without insecticide.
For each age of worker bee, we compared the gene expression levels between bees fed with imidacloprid and the 0.1% DMSO control to examine the DEGs to understand the impact of imidacloprid on honey bee gene expression at various ages. For each age of bee, three comparisons against the solvent control (bees fed with 1 µl of 0.1% DMSO) were performed: (i) 1 µl of 1 ppb imidacloprid solution; (ii) 1 µl of 10 ppb imidacloprid solution; and (iii) 1 µl of 50 ppb imidacloprid solution (Supplementary Figure 1B). DEGs with FDR values < 5% were considered differentially expressed and selected for further analysis. The DEGs were then filtered based on their log 2 fold change. DEGs with a log 2 fold change < 1 or >−1 were discarded, while genes with a log 2 fold change ≥ 1 or ≤−1 were considered twofold differentially expressed. DEGs with a log 2 fold change ≤−1 were defined as "down-regulated, " and DEGs with a log 2 fold change ≥ 1 were defined as "up-regulated" ( Table 1 and Figure 1A). The list of DEGs is shown in Supplementary Tables 3A,B. The detailed analysis procedure is described in section "Materials and Methods." A cluster analysis was performed using read counts to visualize patterns of DEGs at various developmental stages ( Figure 1C). Complete normalization was applied, and the clustering was performed based on Pearson correlation coefficients. For 9-dayold larvae, honey bees fed with 50 ppb imidacloprid showed distinct expression patterns and were assigned as the outgroup, while patterns of DEGs, between 1 and 10 ppb imidacloprid, and between H 2 O and DMSO, were similar, with only slight Fold change 9-day-old larvae 0-day-old adults heads 7-day-old adults heads 14-day-old adults heads 20-day-old adults heads 1 ppb 10 ppb 50 ppb 1 ppb 10 ppb 50 ppb 1 ppb 10 ppb 50 ppb 1 ppb 10 ppb 50 ppb 1 ppb 10 ppb 50 ppb FIGURE 1 | Numbers of differentially expressed genes at various developmental stages of worker bees. (A) Numbers of differentially expressed genes (DEGs) in 9-day-old larvae, as well as in 0-, 7-, 14-, and 20-day-old adults of honey bees exposed to 1-, 10-, or 50-ppb imidacloprid treatment for 4 days during the larval stage. Original numbers of DEGs are shown on top of each column. X-axis: the developmental stages of honey bees; Y-axis: the number of differentially expressed genes. (B) Venn diagram of shared DEGs of 9-day-old larvae, 0-, 7-, and 14-day-old adults among all concentrations of imidacloprid. The numbers of shared DEGs are shown in the overlapped region. Larvae: 9-day-old larvae; 0d: 0-day-old adults; 7d: 7-day-old adults; 14d: 14-day-old adults. (C) Cluster analysis of DEG expression patterns of 9-day-old honey bee larvae, as well as 0-, 7-, and 14-day-old adults. Read counts were used for cluster analysis. Complete normalization was applied, and the clustering was performed based on Pearson correlation coefficients. The expression trend of each replicate is shown. Controls (H 2 O and DMSO) and different concentrations of imidacloprid treatments (1, 10, and 50 ppb) are labeled in the bottom panel and marked with different colors in the panel of the dendrogram. Red: expression levels higher than average; blue: expression levels lower than average; white: expression levels close to average. differences between one replicate of the DMSO solvent control and one replicate of the 1-ppb treatment compared to the other replicates. For 0-day-old adults, a clearer grouping between controls (DMSO and H 2 O) and imidacloprid treatments (1, 10, and 50 ppb) was observed, although one replicate of the 1-ppb treatment had expression patterns closer to the H 2 O control rather than the other two 1-ppb replicates. For 7-day-old adults, the clustering was chaotic and correlations between controls and treatments could barely be identified. For 14-day-old adults, with vast numbers of DEGs identified at this stage, significantly distinct groups between controls and imidacloprid treatments were observed without ambiguous distributions. Cluster analysis for 20-day-old adults was not performed as less than 30 DEGs were identified. After being exposed to a sublethal dose of imidacloprid, the gene expression of honey bee larvae and adults was significantly affected, and the effects were long-lasting ( Table 1). The number of DEGs increased with the dosage of imidacloprid in 9-dayold larvae (18, 921, and 1,217 DEGs for 1, 10, and 50 ppb, respectively) and 0-day-old adults (381, 903, and 1,141 DEGs for 1, 10, and 50 ppb, respectively). For 7-day-old adults, 600 DEGs were identified in bees fed with 10 ppb imidacloprid, while no DEGs were identified in the other treatments (1 and 50 ppb imidacloprid). Surprisingly, more than 4,800 DEGs were observed in 14-day-old adults (4,871, 5,863, and 5,848 DEGs for 1, 10, and 50 ppb, respectively), and the number of DEGs seems to have increased slightly with the dosage of imidacloprid. For 20-day-old adults, fewer than 100 DEGs were identified for all treatments (Table 1) (20, 29, and 12 DEGs for 1, 10, and 50 ppb, respectively). For 9-day-old larvae and 0-day-old adults, there are fewer upregulated than downregulated genes. In contrast, for 7-, 14-, and 20-dayold adults, a greater number of DEGs were identified in the upregulated groups.
The numbers of unique and shared DEGs among different concentrations of imidacloprid treatment at every age of bee were examined to evaluate the effect of the dosage of imidacloprid (Supplementary Figure 2). DEGs were first divided based on the regulatory trend (i.e., upregulated or downregulated), and then the shared DEGs were examined. The numbers of shared DEGs among the different concentrations of imidacloprid treatments are shown in Supplementary Figures 2A-D. In 9-day-old larvae and 0-and 20-day-old adults, the percentage of shared DEGs was less than 30%; more DEGs were found to be uniquely identified in the different concentrations of imidacloprid treatments. In 14-day-old adults, the percentage of shared DEGs was higher than 50%. The shared DEGs in 9-day-old larvae and 0-to 14day-old adults were then examined among all the imidacloprid treatments, disregarding the regulatory trend. Thirty-five genes were identified as shared among different aged bees and different concentrations of imidacloprid treatments, including gammaaminobutyric acid receptor subunit beta (GB40975), DNA ligase 1-like (GB42264), cuticular protein (GB46297), RNA polymerase II elongation factor EII (GB48521), insulin-like growth factor 2 mRNA-binding protein 1 (GB52056), nicotinamide riboside kinase (Nrk1) (GB53410), trehalose transporter 1 (GB55302), and several other unknown-function-genes ( Figure 1B and Table 2).

Gene Ontology Analysis of DEGs
To further understand the biological meaning of the DEGs, Gene Ontology (GO) analysis was performed to predict their functions. DEGs of every treatment were first divided into two groups, i.e., upregulated and downregulated genes, based on their expression trend. An overview of the results of the "functional annotation chart" among all the treatments is shown in Supplementary Tables 3C,D. A GO term that could be identified in at least three ages of bee was considered as a constantly affected function. In the upregulated genes, GO terms identified in at least three ages of bees were as follows: Coiled coil, ion transport, nucleotide binding, Nucleotide-binding (alpha-beta plait), Ribosome, Ribonucleoprotein, Ribosomal protein, structural constituent of ribosome, Transport, and translation. Among the downregulated genes, most of the GO terms can be only identified in two ages of worker bee (Supplementary Tables 3C,D).
In general, the numbers of functional terms were correlated with the numbers of DEGs. For 9-day-old larvae, the number of downregulated clusters (79 and 35 for 10 and 50 ppb, respectively) was higher than that of upregulated clusters (7 and 14 for 10 and 50 ppb, respectively). The opposite pattern (i.e., a higher number of upregulated clusters) was observed in 0-, 7-, 14-, and 20-day-old adults. The numbers of functional terms of upregulated DEGs for 0-, 7-, 14-, and 20-day-old adults were as follows: 0-day-old adults, 16, 34, and 45 for 1, 10, and 50; 7-day-old adults, 23 for 10 ppb; and 14-day-old adults, 52, 91, and 81 for 1, 10, and 50 (Supplementary Table 3C). The numbers of GO cluster numbers of downregulated DEGs were as follows: 0-day-old adults, 8, 8, and 26 for 1, 10, and 50; 7-day-old adults, 10 for 10 ppb; and 14-day-old adults, 16, 27, and 37 for 1, 10, and 50 (Supplementary Table 3D). For upregulated terms, most of them were identified in 9-day-old larvae or 0-, 7-, and 20-day-old adults can be found in 14-day-old adults, while few functional clusters were only identified in 9-day-old larvae or 0-day-old adults. Terms such as ATPase/dynein-related/AAA domain, Dynein heavy chain, Frag1/DRAM/Sfk1, integral component of membrane, and transferase activity/transferring hexosyl groups were uniquely found in 9-day-old larvae, while 3 ,5 -cyclic-nucleotide phosphodiesterase activity, 3 ,5 -cyclic nucleotide phosphodiesterase/catalytic domain, Biotin/lipoyl attachment, Cyclic nucleotide-binding domain, Developmental protein, DnaJ domain, Groucho/transducin-like enhancer, Neuroactive ligand-receptor interaction, Nuclear hormone receptor/ligand-binding/core, Potassium channel, voltage-dependent, EAG/ELK/ERG, Ras-association, Serinethreonine/tyrosine-protein kinase catalytic domain, Steroid hormone receptor, Synaptotagmin, voltage-gated potassium channel activity, Zinc finger, and nuclear hormone receptor-type were found only in 0-day-old adults (Supplementary Table 3C). Terms of downregulated genes were found to be more diverse among different ages of bee, suggesting that the negatively affected genes may vary or change among different ages or tasks of worker bees. Higher numbers of DEGs were identified in 14-day-old adults than in 9-day-old larvae, yet fewer GO terms were identified in 14-day-old adults than in 9-day-old larvae. Nine-day-old honey bee larvae are transforming from larval to pupal, to finally adult stage, which may be involved in massive transcriptome reprogramming and diverse functional pathways for the differentiation and development of adults' organs. High numbers of GO terms may suggest that various functional pathways could be affect after exposure to sublethal dose of imidacloprid.

High Transcriptome Similarity Between Imidacloprid-Treated 14-Day-Old Bees and 20-Day-Old Control Bees
The upregulation of foraging regulators in 14-day-old imidacloprid-treated adults suggests that the switch from nurse bee to forager was accelerated. To confirm this hypothesis, the transcriptome of 14-day-old bees treated with imidacloprid during larval stage was compared to that of 0-, 7-, and 20-day-old control bees (see Supplementary Figure 1C for a schematic diagram of comparison). A comparison between different ages of control bees was also performed (see Supplementary Figure 1D for a schematic diagram of comparison, Supplementary Figure 3 for results). Compared to 0-and 7-day-old controls, the total numbers of DEGs were 4,073 (14-day-old 1 ppb vs. 0-day-old control), 5,491 (14-day-old 10 ppb vs. 0-day-old control), 5,140 (14-day-old 50 ppb vs. 0-day-old control), 8,298 (14-day-old 1 ppb vs. 7-day-old control), 7,696 (14-day-old 10 ppb vs. 7-day-old control), and 7,436 (14-day-old 50 ppb vs. 7-day-old control) (Figures 2A,B). In the comparison of 14-day-old treatments to the 20-day-old control, the numbers of DEGs were negatively correlated with the concentration of imidacloprid. The numbers of DEGs identified from 14-day-old 1 ppb vs. 20-day-old control, 14-day-old 10 ppb vs. 20-day-old control, and 14-day-old 50 ppb vs. 20-day-old control were 4,743, 1,686, and 447, respectively ( Figure 2C). The decreasing number of DEGs suggests the increase in transcriptome similarity, and the effect was dose-dependent. A cluster analysis of solvent controls and imidacloprid-treated workers in 9-day-old larvae and 0-, 7-, 14-, and 20-day-old adults was also performed to visualize expression patterns of whole transcriptome among different ages of bee. Fourteen-day-old workers with imidacloprid treatment during the larval stage were correlated with 20-day-old solvent control and treatments, suggesting the transcriptome similarity among them (Supplementary Figure 4).

Fluctuation of DEG Numbers Among the Ages of Worker Bees
The number of DEGs decreased with the progress of worker bee development from 9-day-old larva to 7-day-old adult, peaked in 14-day-old adults, and then dramatically decreased to less than 30 in 20-day-old adults ( Table 1 and Figure 1A). The different concentrations of imidacloprid seem to have induced different effects, as the DEGs were diverse ( Supplementary Figures 2A-D). The high similarity of expression patterns between biological replicates in every treatment excludes the possibility of individual variance or artifacts, and distinct DEG patterns were observed between controls and imidacloprid treatments (Figure 1C and Supplementary Figure 5 for sample correlation). The fluctuating DEG numbers suggest that the global gene expression in the honey bee changes with the progress of their development as well as task performance. As nurse bees and foragers contain distinct transcriptomes (Khamis et al., 2015), it is highly possible that the gene expression profile would change with the different tasks performed by the nurse bees. The upregulation of the nursing regulator at 0-day-old adults suggests that those DEGs could be correlated with the development of the nursing behavior (Table 3), which should be performed while the hive bees become 3-day-old or older. A similar trend was also observed for the foraging regulators for 14-day-old adults ( Table 3). The DEGs identified at 9-day-old larvae and 0-and 7-day-old adults may reflect not only the damages but also the early maturation induced by sublethal imidacloprid treatment. The decrease in DEGs from 9-day-old larvae to 7-day-old adults may suggest that the task performance of nurse bees is attributed by a small subset of genes. Further experiments are needed to validate the behavioral change of newly emerged workers and the taskassociated genes of nurse bees.

General Effects of Neonicotinoid Treatment Among Different Species of Insects
We compared our DEG results with previous studies (Forfert and Moritz, 2017;Shi et al., 2017a) and found that 329 DEGs were shared (Supplementary Table 18A). These genes are functionally related to lipid metabolism, pheromone/general odorant-binding protein, cytochrome P450, transmembrane transport, and ATP binding (Supplementary Tables 18B,C,D).
Physiological experiments have confirmed that imidacloprid and clothianidin can impact energetic/nutrient homeostasis of honey bee (Cook, 2019). These genes could be toxification index genes expressed to inspect the existence of neonicotinoids in the honey bee. In bumble bee (Bombus terrestris) workers, exposure to field-realistic imidacloprid concentrations (10 ppb) can induce a total of 405 DEGs with functions related to energy reserve metabolism, synaptic transmission, apoptotic processes, xenobiotic transport, and complement activation, and including P450 genes (Manjon et al., 2018); exposure to a sublethal dose of clothianidin (6.47 ppb) induces the differential expression of 55 genes, including three putative cytochrome P450 genes (Bebane et al., 2019;Colgan et al., 2019). In Bradysia odoriphaga, exposure to the LC 50 dose (4.08 mg/l) of imidacloprid for 24 h induces the differential expression of 674 unigenes with functions related to cytochrome P450, cell adhesion, axon guidance, and carbohydrate metabolic process (Chen et al., 2019). The effects of neonicotinoids on insect gene expression vary among different species of insects, but genes related to metabolism, detoxification (especially cytochrome P450), and neuron development/synaptic transmission seem to be commonly found. It is confirmed from molecular evidence that even at low concentrations, neonicotinoids can have an impact on the neuron development and detoxification of insects.

Forager-Like Transcriptome Induced by Exposure to Sublethal Doses of Imidacloprid During Larval Stage
Honey bee division of labor correlates with the age of the adults, exposure to pheromones, ovary size, biogenic amines, and the proportion of members of the colony currently engaging in each task (Seeley, 1982;Huang and Robinson, 1992;Johnson, 2008Johnson, , 2010Bloch et al., 2009;Sagili et al., 2011;Wang et al., 2012;Reim and Scheiner, 2014). From a physiological perspective, a significant increase in the titers of JH coordinates the transformation from nurse bee to forager (Jaycox, 1967;Jaycox et al., 1974;Robinson, 1987;Fahrbach and Robinson, 1996;Bloch et al., 2009). From a molecular perspective, the switch from nurse bees to foragers is modulated by the differential expression of transcription factors and alternative transcription start sites, resulting in transcriptome changes followed by organized behavioral states (Whitfield et al., 2003;Khamis et al., 2015). In our study, it may not be easy to infer the trend of the JH titer based on the function of the identified DEGs, as the GO term related to JH synthesis/metabolite was not identified. Despite the expression of JH-related genes, the upregulation of the foraging regulators and the low transcriptomic difference between imidaclopridtreated 14-day-old adults and 20-day-old control adults suggest that the precocious worker bee may leave the beehive. Colin et al. (2019) examined the effect of a sublethal dose of imidacloprid by providing the beehive with a sucrose solution containing 5 ppb of the insecticide for 3 weeks and then monitored the foraging behavior of newly engaged foragers; foragers exposed to the insecticide started the orientation flight at an age 1.38 days younger than the control (Colin et al., 2019). Although the treatment methodology was different from the present study, it supports the relationship between the insecticide and precocious behavior. The occurrence of precocious foraging is believed to be triggered by a high mortality of foragers (Beshers and Fewell, 2001;Leoncini et al., 2004;Perry et al., 2015;Ushitani et al., 2016). The young foragers are found to be less efficient and have poor spatial memory and shorter lifespan (Chang et al., 2015;Perry et al., 2015;Ushitani et al., 2016) and cannot sustain the foraging task for long. This accelerated behavioral switch in our study is unlikely to be induced by the low number of foragers, as the experimental individuals were released into a healthy colony with task-performing foragers and food storage. We cannot exclude, however, that other unknown factors, such as pathogens or parasites, could have affected the task transition. Further validation and experiments are required to confirm if this phenomenon is a general effect, as two local colonies were used in this report.
Most of the differentially expressed nAChRs were found in 14day-old treated adults regardless of the dose of imidacloprid. This might be correlated with the development of foraging behavior in adults over 10 days of age (Khamis et al., 2015).
Food reward-induced olfactory stimuli are required to perform food collection and foraging behavior (Menzel, 1999;Dupuis et al., 2011;Giurfa and Sandoz, 2012). It is also possible for this differential expression to be associated with the development of precocious foraging behavior in the 14-day-old treated adults, as previously described due to the transcriptome similarity between these adults and the 20-day-old control. In the 35 shared DEGs among the four ages, we also found nicotinamide riboside kinase (Nrk1) (GB53410). This gene is induced in eukaryotes by the presence of nerve damage to protect damaged neurons from degradation (Belenky et al., 2007). Differential expression of Nrk1 among all ages of treated bees suggests that the neurons were damaged after the sublethal imidacloprid treatment. These results support previous findings that the exposure to sublethal doses of imidacloprid during the larval stage can induce neuron damage and consequently impair honey bee adults' olfaction and related learning behavior (Medrzycki et al., 2003;Forfert and Moritz, 2017).

AUTHOR CONTRIBUTIONS
E-CY and Y-RC designed the experiments. Y-RC, CT, P-SH, and T-HW performed the imidacloprid treatment and honey bee sampling. Y-RC prepared the RNA-Seq library. DT performed the data analysis. Y-RC, DT, E-CY, and SZ wrote the manuscript. All authors contributed to the article and approved the submitted version.