RNA-seq reveals disruption of gene regulation when honey bees are caged and deprived of hive conditions

The honey bee Apis mellifera is an important eusocial insect pollinator (Giannini et al., 2015; Morandin et al., 2001; Sampson and Cane, 2000). Although the debate on the overall impact of pesticides on honey bee health is ongoing (Stokstad, 2017), honey bee exposure to agricultural pesticides, particularly the neonicotinoids, is one of many factors that contribute to bee population decline (Sanchez-Bayo and Goka, 2014; Tosi et al., 2017; Williamson et al., 2013; Woodcock et al., 2017). Measuring pesticide toxicity for bees under field conditions is difficult because of the complexity of honey bee biology and foraging behaviors (Alburaki et al., 2015, 2017b; Cutler and Scott-Dupree, 2007; Stewart et al., 2014). Therefore, cage experiments are commonly used to test bees (e.g. toxicological and behavioral assessments) under more controlled conditions (Alburaki et al., 2017a; Gregorc et al., 2017). Despite providing more controlled conditions, little is known about the physiological and molecular responses of bees to caging stress.

Honey bees live in highly organized colonies in which each worker performs a very specific function and changes tasks throughout her physiological development (Winston, 1987). Absence of the queen or any disruption in the pheromone communication inside the bee colony is considered fatal for the whole population. When bees are caged for experimental purposes, significant pheromone disruption occurs, altering bee social behavior (Grozinger et al., 2003). For example, absence of the queen mandibular pheromone delays honey bee behavioral maturation (Robinson et al., 1998), and its presence prevents the rearing of new queens (Winston et al., 1991) and inhibits worker ovary development (Hoover et al., 2003). The ability of worker bees to execute their appropriate tasks and communicate with their hive mates is crucial for colony wellbeing.

Gene regulation in honey bees is often studied in relation to pathogen infection (Gregorc et al., 2012), exposure to abiotic stressors (Alburaki et al., 2017b), and other physiological and genetic factors among honey bee castes and developmental phases (Evans and Wheeler, 1999). These efforts have led to a better understanding of honey bee gene regulation and the identification and annotation of an important number of antioxidant and immune genes evolved in honey bee response to various stressors (Corona and Robinson, 2006; Evans et al., 2006).

Imidacloprid, a broadly used neonicotinoid, is acutely toxic to bees and can impair honey bee performance at sublethal doses (Chakrabarti et al., 2015; Williamson et al., 2014). Imidacloprid has a high agonistic affinity with nicotinic acetylcholine receptors (nAChR), particularly in brain tissue (Decourtye et al., 2004). Bees are equipped with complex detoxification mechanisms, including networks of enzymatic antioxidants, to help reduce potential oxidative damage provoked by biotic and abiotic stressors, including pesticides (Corona and Robinson, 2006). Interestingly, caged bees administrated ad libitum sugar syrup containing 5–10 times their median lethal concentration (LD₅₀) of imidacloprid have shown the ability to survive for a relatively long period of time (4–7 weeks) (Alburaki et al., 2017a; Meikle et al., 2016). It is yet unclear whether those bees survive by minimizing their contaminated-syrup intake or have a metabolic capacity allowing them to quickly process the lethal molecules, or by the two scenarios together. It is also conceivable that the honey bee has the ability to sense the presence of those chemicals within the proposed syrup and refrains from feeding on such contaminated nutrients. It has been demonstrated that honey bees fed on neonicotinoids (thiamethoxam), 1 ng per bee for 12 days, significantly decreased their affinity for sucrose stimulation (Aliouane et al., 2009), which can partially explain the rejection of syrup with optimal sugar concentrations for bees. However, acetamiprid which is another neonicotinoid, orally administrated to bees (0.1 μg per bee) increased their responsiveness to water (Aliouane et al., 2009).

This study was conducted on sister honey bees exposed to two different conditions: (1) honey bees were removed from their natal hive and placed in cages in the laboratory and (2) honey bees were labelled and kept in their natal hive. We addressed two major and interrelated hypotheses: (1) compared with an untreated control, sister honey bees caged under similar conditions and administrated ad libitum sugar syrup tainted with lethal and sublethal concentrations of imidacloprid would exhibit higher mortality and upregulation in their detoxification genes; and (2) depriving honey bees of hive conditions by caging them would trigger significant disruption of gene regulation compared with their sister-mates of the natal hive. This study provides new insights into how imidacloprid affects honey bee physiological responses and diet behavior as well as its implications for the regulation of the major bee antioxidant and immune genes. Furthermore, we characterized potential genes involved in the response to cage stress by studying the differentially expressed genes of caged bees and their sister-mates operating under natural hive conditions.

One well-established honey bee colony, the ‘mother hive’, headed by a Carniolan queen (Apis mellifera carnica Pollman 1879) was used as the source of all worker samples of this study. A single hive was used as a source of worker bees in this study to minimize variability in hive conditions and genetic make-up of the target organisms. A total of eight capped worker brood frames ready to hatch were removed from this hive and placed in an incubator at 35°C with 70–80% relative humidity. The following day, several thousand 1 day old sister bees were collected into a sterile plastic box for this study.

Newly hatched worker bees collected in the plastic box were gently mixed and 500 of them were randomly picked and marked with a white dot on their dorsal tergite using a collective marking box (Alburaki et al., 2017a). These workers were placed in a new plastic box to allow the paint to fully dry, then these marked 1 day old bees were introduced back into their natal hive. Another set of ∼150, 1 day old bees were stored at −80°C and were designated as time 0 reference bees (Fig. 1). A third set of 1800, 1 day old bees were randomly divided into 12 groups (150 bees per group); each group was weighed separately and introduced into 12 separate cages. The cages used in this study (dimensions of 11.4×6.3×15.2 cm width:depth:height; Fig. S4) were specifically designed for feeding experiments and are fully described in Gregorc et al. (2018). The 12 bee cages were randomly assigned to four treatment groups (three treatments and one control) in triplicate (Fig. 1). Bees in each cage were provided with distilled water and 1:1 sugar syrup using 30 ml syringes, and 10 g of protein patty (Mann Lake Bee Pro Patties). Imidacloprid was administrated to bees through the sugar syrup at four different dosages (0, 5, 20 and 100 ppb; Fig. 1). Patty, water and syrup consumption and bee mortality were documented daily. Dead bees were collected daily from the cages and stored at −80°C. The protein patties, which were placed into rubber plugs, were weighed using a ±0.01 g sensitive scale and gently placed back in the cages. Both syrup and water consumption were recorded visually from syringes at 0.5 ml sensitivity.

Worker bees were sampled at four time points: 25 workers were sampled weekly from each cage as well as 75 marked workers from the mother hive. Each cohort of workers was weighed and stored at −80°C for subsequent molecular analyses. In addition, dead bees collected daily from each cage during the experiment were counted and weighed at the end of the experiment. A representative set of dead bees from each treatment (0, 5, 20, 100 ppb imidacloprid) was sent for chemical pesticide residue analysis. Four neonicotinoid molecules were screened (imidacloprid, imidacloprid olfen, thiamethoxam and clothianidin) by liquid chromatography-mass spectrometry (LC-MS) (Barnett et al., 2007; Walorczyk and Gnusowski, 2009). Chemical analyses for pesticide residue detection were processed at the USDA National Scientific Laboratories (Gastonia, NC, USA).

Hydrogen peroxide and protein carbonyl content assays were conducted to assess potential honey bee physiological stress and protein damage resulting from exposure to imidacloprid and caging stress. The level of H₂O₂ was quantified in the hydrogen peroxide assay using the bee hemolymph of samples collected in week 1. Bees were individually crushed in 1.5 ml tubes with 300 μl ultra-sterilized water and centrifuged at 11,000 g for 3 min. In order to eliminate proteins, the supernatant containing the bee hemolymph was filtered through a 10 kDa filter and the assay was conducted using a BioVision kit (Milpitas, CA, USA) as per the manufacturer's instructions. The protein carbonyl content assay was carried out on samples from weeks 1 and 4. Proteins were solubilized from honey bee thorax in a protein extraction buffer consisting of 20 mmol l⁻¹ Tris-HCl pH 8.0, 30 mmol l⁻¹ NaCl and 10% glycerol. The tissues were crushed using a pestle and sonicated using a Bioruptor Pico (Diagenode) sonication device for 10 cycles of 30 s pulse and 30 s rest at 4°C. Homogenates were centrifuged at 5000 g for 10 min at 4°C and the supernatants were collected. Quantification was conducted using a kit from Sigma-Aldrich (St Louis, MO, USA) as per the manufacturer's protocol.

Bee brains were dissected under a binocular microscope in an RNase free, sterile work environment. Total RNA was extracted from a pool of 15 bee brains per sample following the TRIzol^® Reagent protocol (Invitrogen) (Chomczynski, 1993) with some modifications (see Alburaki et al., 2017a). Briefly, dissected brains were added to 1 ml TRIzol with 5 mg of acid-washed glass beads and gently mixed for 2 min. Then, 200 µl of phenol–chloroform was added, and the total mixture was incubated at room temperature for 15 min followed by a centrifugation at 10,000 g for 15 min at 4°C. The integrity of the RNA was determined by using a Nanodrop (260 nm/280 nm absorbance) and RNA concentration was brought to ∼200 ng µl⁻¹ and stored at −80°C.

In total, 45 honey bee genes involved in various functions (detoxification, immune defense and chemosensory roles, physiological development and nervous system regulation) were studied (Table S1). Approximately 200 bp of each exon was processed using the online primer picking program Batch Primer 3 (https://probes.pw.usda.gov/batchprimer3/) under the default setting for general primers to select forward and reverse primers that produce amplicons of 60–70 bp in length (You et al., 2008). The resulting primers were then tested in a single multiplex mixture using bee genomic DNA as the template. The genomic DNA was extracted using a MagJET system according to the manufacturer's directions (ThermoFisher). All amplifications from genomic DNA or cDNA were accomplished by Floodlight Genomics (Floodlight Genomics, Knoxville, TN, USA) through their no-cost Educational Outreach Program using an optimized Hi-plex approach (Nguyen-Dumont et al., 2013). Briefly, primers were mixed into a single multiplex mixture and amplified according to the parameters previously outlined for the Hi-plex approach. The resulting amplicons have a sample-specific barcode sequence incorporated during the PCR and all amplicons were pooled and a dual-index Illumina library constructed and quantified using a KAPA PCR-free kit according to the manufacturer's directions (Roche Sequencing and Life Sciences). Sequencing was accomplished on an Illumina HiSeq device running a 2×150 configuration (Novogene). The resulting sequences were mapped to the original target sequences using CLC Genomics Workbench version 9.5.2 (Qiagen) at default settings and the sequence coverage assessed manually to determine targets with no amplification.

Total RNA was converted to cDNA using the Promega GoScript reverse transcription system (Madison, WI, USA) according to the manufacturer's directions. The resulting cDNA was used as a template for multiplex amplifications as described above. Each amplification was accomplished twice and an RNA template (no cDNA) reaction was included as a third reaction for each sample. The resulting amplicons were sequenced and mapped as described above. To normalize the data and determine the relative gene expression, the total number of reads mapping to each target was divided by the average sequence coverage of four housekeeping genes (Actin, CaMKII, GAPDH and E2F) (Scharlaken et al., 2008). The sequences of the gene-specific primers have been published in Scharlaken et al. (2008) and Alburaki et al. (2017a).

Statistical analyses and figure generation were conducted within the R environment (http://www.R-project.org/). Analysis of variance (ANOVA) was performed to study the difference between variables regarding the treatment at a 95% confidence level. Sequences were checked for their quality using FastQC software and normalized using four housekeeping genes as mentioned above. A few genes that were very weakly amplified were discarded from the dataset. Bioconductor package EdgeR (Robinson et al., 2010) version 3.6 was used in R environment version 3.4.3 for differential expression analyses of read counts arising from our RNA-seq data (Robinson and Oshlack, 2010). Simple list-based data objects (DGEList) were created using the function readDGE, and genes with very low counts were filtered not by direct count but with counts per million (cpm). RNA composition was normalized using the function CalcNormFactors by finding a set of scaling factors that minimize the log-fold change (logFC) between samples for most genes, and these scale factors were computed using a trimmed mean of M-values between each pair of samples (Robinson and Oshlack, 2010). Heatmaps were carried out using the library ComplexHeatmap of the package bioclite Limma. Other plots such as multi-dimensional scaling (MDS), hierarchical clustering and differential gene expression (DEG) were generated using their appropriate and respected functions in EdgeR. DEG data were calculated using EdgeR's Fisher's exact test after sample normalization. False discovery rates (FDRs) in DGEExact object were calculated at three different P-cutoffs (0.05, 0.01, 0.001) and their values were retained for each gene found to be upregulated or downregulated at 5%.

Average syrup intake per week showed no differences among the groups until week 4, when bees in the 100 ppb imidacloprid treatment significantly reduced their syrup consumption (P<0.05, Fig. 2). Similar results were obtained when accounting for overall syrup consumption (P<0.01, Fig. 2). Interestingly, all groups treated with imidacloprid consumed a significantly (P<0.01) higher quantity of water than the control treatment group (0 ppb, Fig. 2), while protein patty consumption did not differ among treatments (data not shown). The correlation coefficients of syrup consumption and imidacloprid concentration indicated significant negative correlations in all treatment categories, particularly for 100 ppb (R=−0.15, P<0.001), with the exception of the 5 ppb treatment, in which no correlation was found (Fig. 2).

The 1 day old sister bees used in cage and hive experiments had similar body mass at time 0, at week 3 and for overall mass (Fig. 3). However, significant differences between the mass of caged and hive bees were recorded at weeks 1, 2 and 4, with no differences among caged bees (Fig. 3). In the cage experiment, there were no significant differences in the total number of dead bees collected throughout the experiment (P=0.7) amongst the treatments. However, dead bees of the control treatment were significantly (P<0.05) greater in mass than those in all other treatment groups (Fig. 3). LC-MS pesticide residue analysis of the dead bees showed trace levels of imidacloprid and imidacloprid olfen only in bees of the 100 ppb treatment (Fig. 3).

At week 1, caged bees of the control treatment (0 ppb) showed a significantly higher level of hydrogen peroxide (P<0.001) than all other groups including the hive bees (Fig. 4). Protein carbonyl content was significantly higher in caged bees given the highest concentration of imidacloprid (100 ppb; P<0.05) at week 1. However, the same assay applied on samples at week 4 showed higher carbonyl content in bees of the 5 ppb group (Fig. 4). Similar results were obtained for whole-bee body protein conducted on samples at week 4 (Fig. 4).

In order to test the effect of both factors (imidacloprid and caging stress) on honey bee gene regulation, DEG analyses were carried out among caged bees only (imidacloprid versus control) and between caged and hive bees (caged versus hive). To gain an initial overview of the data, heatmaps for the whole dataset, including all samples, were generated. Major variation in color pattern was visually identifiable between the hive bee group and the rest of the samples (Fig. S1). We refined our analysis by running it on a weekly basis (Fig. 5). In the caged bees, imidacloprid showed no DEG at any time in the four studied weeks (Fig. S2). The only variation in gene regulation was found between caged and hive bees for both per-week (Fig. 5) and overall week data (Fig. S3). The biological coefficient of variance constantly distinguished two major groups of variables across the weeks with DEGs (in-hive bees versus caged bees, Fig. 5A–C). Heatmaps carried out weekly produced similar findings; genes with potentially different regulation are marked with asterisks in Fig. 5A–C. Based on their function (Table S1), antioxidant genes were the major DEGs in caged bees compared to hive bees in all weeks. At week 1, caged bees upregulated seven antioxidant genes and downregulated nine with no upregulation in the developmental and hormonal genes (Fig. 5A). Relatively similar regulation was observed in week 2 (Fig. 5B), with an increasing downregulation in all DEG categories up to week 4 (Fig. 5D). At week 4, upregulation of both immune defense and hormone genes depleted completely and a significant downregulation in all DEGs took place in caged bees. Gathering together the weekly DEG data in the caged bees, we identified constant upregulation of two antioxidant genes (Rsod and Trx-1) and downregulation of six others (Trx-2, Nmdar1, Vg, CSP3, AChE-2 and MsrA; Table S2; Fig. 6) Rsod is a poorly known gene that encodes an uncharacterized protein of 1123 amino acids, while Trx-1 is a gene encoding a major honey bee antioxidant of 136 amino acids and belongs to the thioredoxin family.

Honey bees are social insects, living together in highly organized colonies (Winston, 1987). Worker bees carry out different tasks throughout their developmental cycle from nurses to foragers, a behavioral maturation that is mainly governed by gene regulation (Grozinger et al., 2003; Robinson, 2002). Cage experiments offer a more controllable environment and are widely used in honey bee behavioral and toxicological studies. Nevertheless, little is known about the physiological and gene regulation changes in bees when caged and deprived of hive conditions. In this study, we found that imidacloprid induced a slight effect on honey bee diet behavior; however, an overwhelming difference in gene regulation was observed in bees deprived of hive conditions. Cohorts of caged sister bees consumed the same amount of syrup during the first 3 weeks of the study, regardless of the imidacloprid dose (0, 5, 20, 100 ppb) in the syrup. There were negative correlations between syrup ingestion and time for all concentrations except 5 ppb. This correlation was more pronounced at the highest dose (100 ppb, r=−0.15, P<0.001; Fig. 2), and during week 4, bees consumed significantly less syrup in the 100 ppb imidacloprid group. Consequently, overall consumption was less for syrup at the 100 ppb dose (P<0.05). It is not clear whether bees sensed the presence of imidacloprid in the tainted syrup and avoided it (Meikle et al., 2016) or whether a post-ingestive aversion response previously described in other invertebrates occurred (Behmer et al., 2005). Bees exposed to imidacloprid showed similar survival to the control group (P=0.7; Fig. 4). Despite surviving the ingestion of imidacloprid, symptoms of insecticide toxicity such as increased grooming and ventilation, disorientation, slower activity and difficulty in flying (Williamson et al., 2014) were clearly observed in the treated bees. Although only traces of imidacloprid and its metabolite molecule (imidacloprid olfen) were detected in bees fed 100 ppb, imidacloprid did induce significant protein damage in honey bees at weeks 1 and 4 (Figs 3 and 4). Surprisingly, and presumably as a mechanism to minimize the effects of exposure to pesticides, treated bees consumed more water than the control group (Fig. 2). Dead bees in the control group weighed more than those exposed to imidacloprid (P<0.05; Fig. 4), but there were no differences in the mass of living bees (Fig. 3). Pairing these results together, it appears that the evaporation of water content from the bodies of dead bees explains this contradictory finding.

While caged bees showed no weekly differences in mass, hive bees exhibited completely different patterns (Fig. 3). The mass patterns of hive bees clearly reflect their physiological development from nurse to forager bees. Hive bees were significantly heavier than caged bees at week 1 as well as week 2 (P<0.05), until becoming equal in mass at week 3. When they switched to foraging behavior at week 4, hive bees lost significant mass (P<0.001) compared with the relatively inactive caged bees (Fig. 3). Overall, our data indicate that caging bees induces a profound alteration in their physiological development compared with their sister-mates operating under hive conditions.

The DEG analysis showed a clear divergence in gene regulation of both bee sets (caged versus hive). Based on the genes with the five highest FDR values found in our DEG study, gene upregulation is the predominant response during the first 2 weeks in caged bees followed by overwhelming downregulation at weeks 3 and 4 (Fig. 5C,D), which points to potential depletion in caged bees compared with that of bees in the hive. Upregulation of the gene group involved in the antioxidative process was recorded in caged bees as compared with hive bees at all times and particularly during the first 3 weeks (Fig. 5A–D). Some of those genes belong to the thioredoxin family, such as Trx-1, Trx1-like2 and Trxr-1, as well as the glutaredoxin and glutathione families, all known to be involved with Apis cerana’s response to oxidative stress (Meng et al., 2014). Fewer genes coding for other functions were found to be differentially expressed in caged bees, such as Mrjp1, Fed1 and Vg (Table S2). Vitellogenin, a glycolipoprotein, which is often more abundant in nurse bees (Amdam et al., 2005; Ihle et al., 2009; Nelson et al., 2007), was among the six genes that were constantly downregulated (Vg, CSP3, Trx2, Nmdar1, AChE-2, MsrA) in caged bees throughout the 4 week experiment (Fig. 6). The honey bee genome contains 21 genes encoding odorant-binding proteins (Obps; The Honeybee Genome Sequencing Consortium, 2006) which, beside chemosensory proteins (CSPs), mediate both bee perception and release of chemical stimuli and probably have other functions (Iovinella et al., 2011). Obp21 was only upregulated in caged bees at week 3, while CSP3 was constantly downregulated, which suggest different functions of these genes. The constant downregulation of AChE suggests an absence of toxicity stress in caged bees as this gene was described in previous studies to be well linked to neonicotinoid toxicity (Alburaki et al., 2015; Badiou-Bénéteau et al., 2012). In this study, we identified two genes that were constantly upregulated in caged bees; Rsod, an uncharacterized gene, and thioredoxin (Trx1) (Fig. 6; Table S2). Rsod was described in Drosophila as an atypical member of the Cu/ZNSOD family with a highly unusual number of introns (18) and a duplicated SOD domain. Homologous genes were described in Apis, protozoa and fish but not found in mammals (Corona and Robinson, 2006). Nonetheless, the function of this gene (Rsod) remains unknown in insects. The thioredoxin genes (Trx) encode small and highly conserved oxidoreducatse proteins and are required to maintain redox homeostasis of the cell (Holmgren et al., 2005). Thus, we assume that the caging stress challenge may have triggered the bee biological system to mitigate this stress and prevent further cellular damage. This assumption is supported by the significant upregulation of the antioxidant genes and higher protein carbonyl content in caged bees compared with hive bees.

In conclusion, our study provides new insights into honey bee gene regulation when bees are exposed to cage stress compared with a typical hive environment. We showed that major honey bee antioxidants were constantly upregulated when bees are caged, including genes with uncharacterized function such as Rsod and Trx1. Furthermore, to the best of our knowledge, this is the first time that eight genes have been reported to potentially characterize the honey bee caging stress while conducting cage experiments. Our results add a significant contribution to the body of knowledge related to the effect of stressors on honey bee gene regulation.

We are grateful to Dr William Meikle and Milagra Weiss (USDA-ARS Carl Hayden Bee Research Center, Tucson, AZ, USA) for providing material support for this study.