Sex-specific molecular specialization and activity rhythm-dependent gene expression in honey bee antennae

Honey bees (Apis mellifera) and other eusocial insects, which show an elaborate division of labour involving morphologically and physiologically specialized phenotypes, provide a unique system to study the molecular underpinnings of behaviour (Brockmann et al., 2007; Liang et al., 2012; Whitfield et al., 2006). Honey bee drones (males) and workers (females), for example, strongly differ in their behaviour, morphology and physiology. Drones have to find and mate with a virgin queen at the same time as outcompeting other drones (Brockmann et al., 2006; Gary, 1962; Koeniger et al., 2005; Ruttner, 1985), whereas workers do all the tasks needed to maintain the colony and organize the underlying division of labour; and as foragers, they have to learn and memorize odour mixtures that indicate different rewarding flowers (Frisch, 1967). The different behavioural functions have led to different evolutionary adaptations in the visual and olfactory systems (Esslen and Kaissling, 1976; Menzel et al., 1991).

Honey bee antennae are multimodal sense organs housing olfactory, gustatory and mechanosensory sensilla, as well as CO₂-, humidity- and temperature-sensitive sensilla. However, Esslen and Kaissling (1976) estimated that in workers, 96.9% and in drones, 99.6% of the antennal sensory neurons are olfactory. Regarding olfaction, drones are considered olfactory specialists and workers generalists (Masson and Mustaparta, 1990). Drone antennae have about 7 times more olfactory poreplate sensilla, each containing about 15–30 olfactory sensory neurons (Esslen and Kaissling, 1976). Furthermore, in drones, most of these sensory neurons are specialized to detect minute amounts of the queen's sex pheromone (Brockmann et al., 1998; Vareschi, 1971; Wanner et al., 2007). Electrophysiological recordings indicated that the sex pheromone-responsive sensory neurons have lower action potential thresholds and higher firing frequencies compared with sensory neurons responsive to floral odours (Vareschi, 1971). Compared with drones, workers have a more generalist olfactory system, probably with a broader odorant response profile. For example, worker antennae house hair-like olfactory sensilla (s. basiconica) that are almost absent on drone antennae, and worker antennal lobes are composed of about 170 isomorphic olfactory glomeruli, whereas drones have only 100 normal-sized and 4 macro glomeruli (Arnold et al., 1985; Brockmann and Brückner, 2001; Esslen and Kaissling, 1976; Flanagan and Mercer, 1989; Galizia et al., 1999; Kamikouchi et al., 2004; Kropf et al., 2014; Sandoz, 2006; Wanner et al., 2007). Given the different behavioural functions and the strong sexual dimorphism, drone and worker antennal sensory systems may exhibit different sensory reception strategies and molecular adaptations (Masson and Mustaparta, 1990; Sandoz et al., 2007; Zayed et al., 2012). Previously, Wanner et al. (2007) reported differences in the expression of odorant receptor (Or) genes between drones and workers and showed that one of the drone-expressed Or genes (Or11) binds 9-ODA, the major sex pheromone compound. However, the differences in the overall molecular equipment for odorant reception might be larger and more complex between drone and worker antennae, particularly if one considers the possibility of temporal and behavioural state-dependent changes in gene expression.

In this RNA-seq study, we compared the antennal transcriptomes of sexually mature drones and time-trained foragers collected at different times of the day and under different activity states (Jain and Brockmann, 2018; Naeger and Robinson, 2016). The goals of our project were 3-fold. First, we wanted to provide a more comprehensive description of gene expression differences between drone and worker antennae. Most studies on the insect antennal transcriptome still focus on identifying and reporting genes involved in odorant binding and detection (Liu et al., 2016; McKenzie et al., 2016; Missbach et al., 2020). In contrast, we also aimed at identifying so far unrecognized molecules, not directly involved in odorant detection, but likely to play an important role in peripheral olfactory processing. Second, global gene expression differences between drone and worker antennae might provide information about differences in peripheral sensory processing strategies between a specialist and a generalist peripheral olfactory system (Amin and Lin, 2019; De Bruyne and Baker, 2008; Renou, 2014). Finally, we wanted to explore whether daily changes in gene expression and correlations between gene expression levels and behavioural activity might be a fruitful approach to identify additional genes involved in odorant reception.

In all experiments we used Apis mellifera Linnaeus 1758 colonies of naturally mated queens, which consisted of about 8000 workers (i.e. 8 frames with approximately 1000 workers) and hundreds of drones. Colonies were acquired from a local beekeeper and maintained on the campus of the National Centre for Biological Sciences – Tata Institute of Fundamental Research (NCBS-TIFR), Bangalore, India. All experiments including sample collections were performed under natural light–dark conditions during the months of October and November in Bangalore. During this time of the year, sunrise was at around 06:10 h and sunset was at around 18:00 h, which is close to 12 h of light and 12 h of darkness (LD 12:12) cycle (timeanddate.com).

Daily drone flight activity was determined for three colonies on three different days, i.e. three measurements for each colony, during the months of October and November. On the experimental days, the number of drones leaving the hive entrance was counted every half an hour for 10 min from morning (07:00 h) to evening (19:00 h) (LD 12:12; sunrise time 06:10 h, sunset time 18:00 h). It should be noted that drones and foragers of A. mellifera do not fly during the dark. On these days, we also recorded temperature and humidity changes every minute using a data logger (EQ-172, Equinox, Valli Aqua And Process Instruments, Chennai, India).

During daily mating flight activity under natural conditions, drones were caught at the hive entrance and colour marked on the thorax. The next day, colour-marked drones were collected at two different time points – 09:00 h (inactive) and 14:00 h (active/mating flight time) (see Naeger and Robinson, 2016) – from three different colonies (five drones per time point per colony). At 09:00 h, drones were collected from inside the colonies and at 14:00 h they were collected from the entrance before they started their mating flights. In addition, we collected colour-marked drones from one of the three colonies at six different time points: 06:00 h, 10:00 h, 14:00 h, 18:00 h, 22:00 h and 02:00 h (10 bees per time point) to determine daily expression changes of four major clock genes: period (per), cryptochrome2 (cry2), cycle (cyc) and clock (clk). Night collections were done under dim red light. All collected drones were immediately flash frozen in liquid nitrogen.

An A. mellifera colony was transferred in an enclosed outdoor flight cage to entrain the foraging activity of the workers to a distinct time of the day. First, the colony was allowed to adjust to the new environment for 10 days. During this period, the sugar and pollen feeders were presented for the whole day. The sugar feeder was a yellow plastic plate surrounded by four filter papers containing a 5 µl drop of 100 times diluted linalool racemic mixture (Sigma-Aldrich, St Louis, MO, USA). Then, for time training, the sucrose reward (1 mol l⁻¹ sucrose solution) was presented either from 08:00 h to 10:00 h (morning training) or from 13:00 h to 15:00 h (afternoon training) for 10 consecutive days (LD 12:12). Time for the afternoon training was chosen according to the drone flight time based on our behavioural data. Two different colonies were used for morning and afternoon training. Every day after the training time the feeder was cleaned with ethanol and linalool-scented filter papers were replaced with fresh unscented filter papers. This cleaned empty feeder was available for the rest of the day. On the 8th, 9th and 10th day of training, foragers visiting the feeder were marked on their thorax with different colours, one type of colour each day, to identify the frequently visiting foragers. On the 11th day, the feeder was not presented and the foragers that had all three colour marks were collected at 09:00 h and 14:00 h from both the colonies (10 bees per time point per colony). Collected foragers were immediately flash frozen in liquid nitrogen.

Collected honey bees were transferred from liquid nitrogen onto dry ice and the entire antennae (i.e. scape, pedicel and flagellum) were cut off. We pooled 10 antennae from five bees for RNA-seq samples and four antennae from two bees for qPCR samples. Total RNA was isolated using the Trizol^® method (Invitrogen, Carlsbad, CA, USA). The tissue (i.e. antennae) was homogenized in Trizol, a solution of phenol and guanidinium isothiocyanate, at room temperature and then chloroform was added to separate the aqueous and organic phase. The aqueous phase, which contains RNA, was put into a new tube and total RNA was precipitated using isopropanol and sodium acetate. Samples were treated with DNaseI (Invitrogen) for 10 min to remove any possible DNA contamination. Final RNA concentration was measured using a Nanodrop and RNA quality was assessed by running samples on an agarose gel.

To quantify the expression levels of clock genes using qPCR, around 800 ng of total RNA from each sample was converted into cDNA using Reverse Transcriptase Superscript III enzyme and oligo d(T)16 primers (Invitrogen). Then KAPA SYBR FAST qPCR Master Mix (Kapa Biosystems, Wilmington, MA, USA) in 7900HT Fast Real-Time PCR system (Applied Biosystems, Carlsbad, CA, USA) was used to perform qPCR. An internal control gene for normalization was selected after comparing the expression profiles of six different housekeeping genes: tubulin (XM_396338), RpS18 (XM_625101), actin (NM_001185145), ef-1a-f1 (NM_001011628), gapdh (XM_393605) and ribosomal protein49 (rp49) (NM_001011587). rp49 showed the highest stability (NormFinder; Andersen et al., 2004) in our experimental conditions (Fig. S3) and was used for normalization of clock gene expression data. All the primer sequences used in this study are available from figshare (https://doi.org/10.6084/m9.figshare.12089370.v1, Table S1). Triplicate reactions (10 μl reaction mix) for all the biological replicates of all six time point samples (n=5 per time point) were run in parallel on the same 384-well plate. This restricted us to analyse just one of the clock genes and rp49 (internal control gene) per plate. We also ran a standard curve for both primers on the same plate using a separate stock cDNA. Final gene expression calculation was based on the linear values interpolated from the standard curves. Efficiency of all the primers was between 95% and 100%. Reactions with a qPCR dissociation curve that did not support specificity of the PCR reaction were discarded from the analysis.

Antennal transcriptomes of drones (n=3 per time point), morning-trained foragers (n=2 per time point) and afternoon-trained foragers (n=2 per time point) were sequenced for two different time points (09:00 h and 14:00 h). Total RNA from all the samples was shipped on dry ice to AgriGenome Labs (Kochi, India). RNA quality was further checked on Agilent Tapestation and Qubit. Libraries were prepared using TruSeq Ribo-Zero Gold+TruSeq Stranded mRNA Library Prep kit (Illumina, San Diego, CA, USA). Sequencing was performed on an Illumina HiSeq2500 platform and around 120 million (∼7–10 million per sample) of 75 bp paired-end reads were generated.

We used the cosinor package (https://github.com/sachsmc/cosinor) in R (http://www.R-project.org/) to fit a 24 h cosine model {y=intercept+amplitude×cos[2×π(x−acrophase)/24]} (Nagari et al., 2017) in the circadian gene expression data. We also performed a non-parametric Jonckheere–Terpstra–Kendall (JTK) cycle analysis (Hughes et al., 2010; Patton et al., 2014) to detect the phase and amplitude of clock gene expression profiles.

Approximately 7–10 million pairs of 75 bp reads per sample were mapped to the A. mellifera genome (NCBI Assembly Amel_HAv3.1; Wallberg et al., 2019) using STAR RNA-seq aligner (Dobin et al., 2013). The alignment rate was more than 75% (75.15–86.82%) for all the samples. The number of reads aligning to each gene was counted using featureCounts, a program specifically designed for read summarization on genomic features (Liao et al., 2014). The read counts for each gene were normalized with the respected library size and differences in gene expression between the samples were calculated using DESeq2 (Love et al., 2014). Genes showing expression differences with adjusted P-values (P_adj) less than 0.05 (Wald test) were called differentially expressed genes (DEGs). The Pathview package (Luo and Brouwer, 2013) in R was used to map the above differential gene expression data to relevant pathway graphs from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and to visualize. In addition, the GAGE package (Luo et al., 2009) in R was used for gene-set analysis (GSA) using normalized count data from featureCounts and the KEGG pathway database (Kanehisa and Goto, 2000). Gene ontology (GO) enrichment analysis was done using g:Profiler (Raudvere et al., 2019) with an alpha of 0.05 as the cut off for significance.

Drone and worker antennae showed characteristic differences in their transcriptomic profiles (Fig. 1), and the 14 RNA-seq samples (6 drones and 8 workers) separated into two distinct sex-specific clusters (Fig. 1B). We could detect and compare the expression levels of 10,635 genes between drone and forager antennae using our RNA-seq data. Overall, we identified 4004 DEGs (P_adj<0.05, Wald test), of which 1814 were more highly expressed in drones and 2190 were more highly expressed in foragers (Fig. 1A; see also figshare: https://doi.org/10.6084/m9.figshare.12089370.v1, Table S2).

We were able to detect expression of most of the annotated olfaction-related genes, i.e. we detected 146 out of 167 Or genes, 15 out of 21 odorant-binding protein (Obp) genes, all 6 chemosensory protein (CSP) genes, 13 out of 17 carboxylesterase (CEst) genes, 8 out of 10 gustatory receptor (Gr) genes, 34 out of 45 cytochrome P450 (CYP) genes, 7 out of 8 glutathione S-transferase (Gst) genes and 42 out of 51 other cytochrome-related genes. Of the 146 detected Or genes, 91 showed significant expression differences (P_adj<0.05, Wald test) between drone and forager antenna, of which 24 (13 with log₂ fold-change>1) were more highly expressed in drones and 67 (54 with log₂ fold-change>1) were more highly expressed in foragers (Fig. 1C; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S3). In addition, we found significant expression differences (P_adj<0.05, Wald test) for 13 Obp genes (4 drone biased and 9 forager biased), 4 CSP genes (all higher in foragers), 12 CEst genes (5 drone biased and 7 forager biased), 5 Gr genes (2 drone biased and 3 forager biased), 22 CYP genes (8 drone biased and 14 forager biased), 15 Gst genes (all higher in drones) and 26 other cytochrome-related genes (22 drone biased and 4 forager biased) between drone and forager antennae (Fig. 1C,D; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S3).

Regarding the differentially expressed Or genes, 21 out of the 24 drone biased Or genes belong to the prominent clade X (hymenopteran subfamily L) which comprises the 9-ODA-sensitive Or11 (https://doi.org/10.6084/m9.figshare.12089370.v1, Table S4; see Karpe et al., 2016; Robertson and Wanner, 2006). The 67 forager-biased Or genes belong to 13 different Or clades, but most of them are members of 4 clades: clade XXI (n=20, hymenopteran subfamily J, a bee expanded clade), clade X (n=16, the ‘9-ODA clade’), clade XI (n=11, 9 exon – putative CHC receptors) and clade XVIII (n=5, hymenopteran subfamily H, putative floral scent receptors) (see https://doi.org/10.6084/m9.figshare.12089370.v1, Table S4).

Finally, we found almost all olfaction-related genes that were previously reported to be differentially expressed between drones and workers were among the DEGs with the highest expression differences: Or10, Or11 (the 9-ODA odorant receptor), Or18, Or170 and CEst1 (higher in drones; Wanner et al., 2007); Or63, Or81, Or109, Or150, Or151, Or152, Obp2, Obp4, Obp11, Obp16, Obp19, Obp21, CSP6 and Cyp6BE1 (higher in foragers; Forêt and Maleszka, 2006; Wanner et al., 2007) (see Fig. 1D; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S3).

Besides the genes that are directly involved in odorant detection, we found a number of differentially expressed genes (log₂ fold-change>1, P_adj<0.05, Wald test) that are probably involved in olfactory sensory neuron processing or other sensory functions of the antennae (Fig. 1A; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S2). Drone antennae showed higher expression of the two major sex-determining genes, complementary sex determiner (csd) and feminizer (fem), as well as nitric oxide synthase (NOS), and the potassium channel Shaw (Shaker cognate w; see Hodge, 2009).

Forager antennae showed higher expression of several genes involved in biogenic amine signalling (Dop1, Dop3, DopR2, 5-HT2alpha, 5-HT2beta, TyrR), glutamate signalling (ionotropic glutamate receptor, glutamate receptor 1, vesicular glutamate transporter 1, metabotropic glutamate receptor 1 and glutamate decarboxylase 1 genes), enzymes of the tyrosine/dopamine metabolic pathway (tyrosine hydroxylase and tyrosine aminotransferase genes), as well as several genes involved in neuropeptide signalling [adipokinetic hormone receptor (Akhr), tachykinin (Tk), prohormone-2 and neprilysin-4 genes].

Further, the TRPV channel nanchung showed a higher gene expression in forager antennae. nanchung was reported to be expressed in the Johnston's organ and involved in hearing and gravity perception (Ai et al., 2007; Sun et al., 2009). Mondet et al. (2015) previously showed a higher expression of nanchung in forager antennae compared with nurse antennae.

Finally, vitellogenin (Vg), several genes of the major royal jelly and yellow proteins (Mrjp1, Mrjp3, Mrjp4, Mrjp6, Mrjp8, Mrjp9, Y-h, Y-y, Y-e3, Y-f), all likely to be involved in sex- and caste-specific behaviours (Buttstedt et al., 2013; Drapeau et al., 2006), as well as the immune response genes Def1, Def2, abaecin, Apid1 and Tsf1 were more highly expressed in worker antennae.

GSA using the KEGG pathway database revealed significant (q<0.05, GAGE) upregulation of 67 biological pathways in drone antennae and 3 biological pathways in forager antennae (see https://doi.org/10.6084/m9.figshare.12089370.v1, Table S5). Two of the most significant pathways (lowest q-values) in drone antennae were oxidative phosphorylation (ame00190) and protein processing in the endoplasmic reticulum (ame04141) (Fig. S1). In contrast, in worker antennae, the two most significant pathways were neuroactive ligand-receptor interaction (ame04080) and tyrosine metabolism (ame00350) (Fig. S2).

GO enrichment analysis using DEGs (P_adj<0.05, g:Profiler) with more than 2-fold expression differences (539 DEGs in drones and 973 in foragers) showed significant enrichment of 57 and 95 GO terms in drones and foragers, respectively (P<0.05; see https://doi.org/10.6084/m9.figshare.12089370.v1, Table S6). Some of the top GO terms in drones were catalytic activity (GO:0003824), oxidoreductase activity (GO:0016491), metabolic process (GO:0008152), protein folding (GO:0006457), mitochondria (GO:0005739) and mitochondrial protein complex (GO:0098798). In foragers, some of the significantly enriched GO categories were signalling receptor activity (GO:0038023), molecular transducer activity (GO:0060089), regulation of cellular process (GO:0050794), cell communication (GO:0007154), integral component of membrane (GO:0016021) and extracellular region (GO:0005576).

Drones of all three A. mellifera colonies performed mating flights between 13:00 h and 15:00 h on all three observation days (Fig. 2A). The drone flight activity did not differ among colonies and experimental days. During these 2 h, the temperature was about 25°C and the relative humidity was around 60–70% (Fig. 2B).

The antennae of drones performing mating flights showed significant 24 h daily rhythms (JTK cycle analysis, P_adj<0.05) in the mRNA levels of major clock genes (n=5 per time point; Fig. 3 and Table 1). cry2 and per mRNA levels peaked during the early morning, while the cyc mRNA level was highest during the afternoon. clk mRNA levels did not change. cry2 levels oscillated with higher amplitude than those of per and cyc. Interestingly, this expression pattern of clock genes is similar to that of afternoon-trained foragers (Jain and Brockmann, 2018; Spangler, 1972).

Here, we analysed the changes in antennal transcriptome with activity state and time of day. In drone antennae, we found 78 genes that were differentially expressed (P_adj<0.05, Wald test) (Fig. 4A; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S7) between morning and afternoon. As drones fly out only in the afternoon, these 78 DEGs represent a combined list of genes affected by time of day and activity state. In foragers, morning and afternoon feeder time training allowed separation of time of day and activity state effects on gene expression differences. First, by comparing gene expression differences between the active forager group (i.e. morning-trained foragers collected in the morning and afternoon-trained foragers collected in the afternoon) and in-active forager group (i.e. morning-trained foragers collected in the afternoon and afternoon-trained foragers collected in the morning), we identified 17 DEGs (P_adj<0.05, Wald test) that are likely to be regulated by activity state (Fig. 4A; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S7). Second, by comparing gene expression differences between morning and afternoon irrespective of activity state or training time, we identified 50 DEGs (P_adj<0.05, Wald test) that were associated with the time of day (Fig. 4A; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S7). Finally, we also compared the combined transcriptomes of the morning-trained bees (collected in the morning and afternoon) against the combined transcriptomes of the afternoon-trained bees (collected in the morning and afternoon), but did not find any DEGs, suggesting there were no colony differences or differences due to morning and afternoon training itself.

Comparing the identified DEGs of these workers with those of the drones, we found that 2 out of the 17 activity state-regulated worker DEGs (per and LOC107966102) also showed a difference in expression between the drone samples. Likewise, 10 out of the 50 time of day-regulated worker DEGs (including cry2) showed a difference in expression between the drone samples (Fig. 4A,B). Furthermore, 8 out of these 12 common DEGs showed expression changes in the same direction in drones and foragers (Fig. 4B). per was more highly expressed during the inactive state and cry2 was more highly expressed during the morning in both drones and foragers. Similarly, LOC107966102, CUBN, DNAH3, LOC100576934, low-density lipoprotein receptor-related protein 2 and RNA-binding protein 1 showed changes in the same direction in both drones and foragers (Fig. 4B). In contrast, the expression changes in the remaining 4 common DEGs (hsp60A, hsp90, ebony, AAEL008004) were in the opposite direction in drones and foragers (Fig. 4B).

Despite these similarities, most of the DEGs showing differences with respect to activity state or time of day were either drone or worker specific. Sixty-six DEGs showed changes in expression only in the drone antennae (see https://doi.org/10.6084/m9.figshare.12089370.v1, Table S7). Among these genes were jun-related antigen (Jra), Hr38, endoplasmin, SIFamide receptor (SIFR), foraging (for), dopa decarboxylase (Ddc), kruppel homolog-1, semaphorin-2A, prohormone-2 and few heat shock protein (hsp) genes. In foragers, we found 15 and 40 unique DEGs associated with time of day and activity state, respectively (DEGs associated with activity state include: bruchpilot, semaphorin-1A, translation initiation factor 2 (IF-2) and histone H1; DEGs associated with time of day include: takeout, ataxin-2 homolog, Cdk4, netrin receptor gene UNC5B and non-coding nuclear RNA gene Ks-1; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S7).

GSA revealed significant enrichment (q<0.05, GAGE) of 3 pathways (Fig. 4C; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S8): oxidative phosphorylation (ame00190), protein processing in endoplasmic reticulum (ame04141) and ribosome (ame03010). The oxidative phosphorylation pathway (ame00190) showed significant enrichment only in foragers and was strongly associated with their activity state (Fig. 4C). ame04141 and ame03010 were found to be significantly enriched in both foragers and drones (Fig. 4C). In foragers, both gene-sets were up-regulated during active states (during their feeder training time) while in drones, ame04141 was up-regulated during the active state (in the afternoon during their mating flight time) and ame03010 was up-regulated during the inactive state (in the morning).

So far, most studies on sexual dimorphism in antennal gene expression have focused solely on olfaction-related genes (Liu et al., 2016; McKenzie et al., 2016; Missbach et al., 2020) and have only used antennae collected at a single time of the day, neglecting that antennal gene expression might be temporally dynamic. Here, we performed an extensive comparison of gene expression levels between the antennae of honey bee drones and foragers collected at different activity states, i.e. resting versus mating flight activity and resting versus foraging flight activity. Altogether, we detected 4004 DEGs (P_adj<0.05), with 1814 more highly expressed in drones and 2190 more highly expressed in foragers (Fig. 1A; see also https://doi.org/10.6084/m9.figshare.12089370.v1, Table S2).

Interestingly, none of the canonical olfaction-related genes showed significant expression changes with activity state or time of day. Thus, drone and worker antennae show stable sex-specific molecular specialization which nicely correlates with the already described morphological and physiological differences (Brockmann and Brückner, 2001; Esslen and Kaissling, 1976; Vareschi, 1971). We identified 20 new drone-biased Or genes, all but one of which belong to the clade X that comprises the 9-ODA-sensitive Or11 (Karpe et al., 2016; Robertson and Wanner, 2006; Wanner et al., 2007). This finding suggests that more than just 4 Or genes (as previously suggested) might be involved in honey bee sex-pheromone communication. For workers, we found 67 more highly expressed Or genes. Intriguingly, most of these genes belong to the bee expanded clade (clade XXI) as well as Or clades that are associated with the perception of mandibular gland pheromone (clade X) cuticular hydrocarbons (clade XI), and floral scents (clade XVIII; see Karpe et al., 2016; Robertson and Wanner, 2006). Regarding the other olfaction-related genes (Obp, CSP, CEst, Gr, CYP, Gst and cytochrome-related genes), most were more highly expressed in foragers with the exception of the Gst and cytochrome-related genes, which were generally more highly expressed in drones. Gst and cytochrome-related genes are assumed to be involved in signal (pheromone) inactivation, which is very important for fast pheromone-guided flight navigation (Leal, 2013; Rogers et al., 1999).

Among the non-olfaction-related genes showing a higher expression in drone antennae were csd, fem, NOS and, most interestingly, the potassium channel gene Shaw (Shaker cognate w). It is not entirely clear why the two major sex-determining genes csd and fem were more highly expressed in the drone antennae, but Naeger and Robinson (2016) also found a higher expression of these two genes in the mushroom bodies of drones compared with those of workers. There is some evidence that NOS modulates responses of olfactory sensory neurons in insects (Wasserman and Itagaki, 2003). In contrast, as far as we know, there is no report regarding a possible function of Shaw in insect olfaction. However, there is evidence that the sensillum lymph contains high K⁺ concentrations and the membrane potential of olfactory sensory neurons is regulated by the movement of K⁺ ions (Küppers and Thurm, 1979; Mohapatra and Menuz, 2019; Stengl et al., 1999). Also, Shaw was demonstrated to play a role in modulating the activity of clock neurons (Hodge and Stanewsky, 2008).

Although worker antennae have only one-seventh the number of sensory neurons of drones, their antennal tissue shows a higher expression of many genes involved in neural modulation (e.g. biogenic amine signalling, glutamate and neuropeptide signalling). In addition, our GO enrichment analysis suggested higher expression of genes involved in secondary messenger cascades, cell signalling and extracellular matrix in forager antennae. Both findings suggest a higher degree of signalling plasticity in the worker antennae. Plasticity in antennal processing might play a role in division of labour and in learning and memory. Previous studies in honey bees (McQuillan et al., 2012; Vergoz et al., 2009) demonstrated that the expression of dopamine and tyramine receptors is age and task dependent and can be modulated by social pheromones.

Further, we found a higher expression of the TRPV channel gene nanchung in worker antennae. nanchung was reported to be expressed in the Johnston's organ and involved in hearing and gravity perception (Ai et al., 2007; Sun et al., 2009). Mondet et al. (2015) previously showed a higher expression of nanchung in forager antennae compared with nurse antennae. The higher expression of nanchung might indicate a more sensitive or even more elaborate functional organization of the worker Johnston's organ, which is involved in several worker-specific tasks and communication (Ai et al., 2007, 2009; Frisch, 1967; Sun et al., 2009; Towne and Kirchner, 1989).

Finally, our GSA showed that genes involved in oxidative phosphorylation are predominantly more highly expressed in drone antennae. Correspondingly, we found that in workers, metabolic genes are regulated with respect to activity state. Drone antennae have a much higher number of sensory neurons than worker antennae and coding in sensory neurons is based on generating action potentials, which are energy expensive (Attwell and Laughlin, 2001; Niven and Laughlin, 2008). Furthermore, a higher expression of genes involved in protein folding and protein processing might be a consequence of a higher protein turnover rate associated with higher neural activity in drone antennae (Dörrbaum et al., 2018).

Similar to worker antennae, drone antennae exhibit robust daily oscillation of the two major clock genes per and cry2. More interestingly, the daily oscillations of per and cry2 expression in the drone antenna were similar to those in the antennae of afternoon-trained foragers (Jain and Brockmann, 2018; Sasaki, 1990; Spangler, 1972). Given that drone mating flight times are species specific and innately activated but foraging time of workers can be shifted, it is highly likely that these daily oscillations are differently regulated. Furthermore, the comparison of expression levels between time of day and activity state suggest that the expression of the two clock genes is sensitive to different environmental factors. Expression of per seems to be sensitive to activity state (rest–activity cycle), whereas cry2 is more stably associated with time of day (light–dark cycle). Unfortunately, not much is known about the plasticity of the circadian clock with respect to individual activity and changing light–dark cycles (Saunders, 2002).

In addition to the clock genes, we identified a total of 133 genes that show changes with time and activity state in the antennae of drones and foragers. Twelve of these genes (Fig. 4B) were common and 121 genes were different between drones and foragers, suggesting a strong sexual dimorphism in the set of genes that are regulated by the circadian clock and behavioural activity. Certainly, a more detailed (e.g. more collection time points) study that identifies which genes are innately regulated by the circadian clock and which genes can come under the control of the circadian clock by time training will be very instructive.

Finally, the findings of our RNA-seq study corroborate that the functional organization of insect antennae might be more complex than just detecting odorants and transmitting sensory signals to the brain (Andersson et al., 2010; Couto et al., 2005; Getz and Akers, 1994, 1995). Drone antennae are optimized to detect small amounts of the queen's pheromone and quickly respond to changes in pheromone concentration (Brockmann et al., 1998; De Bruyne and Baker, 2008), whereas forager antennae are predominantly involved in context-dependent detection and discrimination of complex odour mixtures, which might require flexible filtering of sensory signals sent to the brain (Conchou et al., 2019; Gadenne et al., 2016; Getahun et al., 2013). Previous extracellular recordings from the olfactory poreplate sensilla indicated that there might be physiological interactions between the olfactory sensory neurons within one sensillum (Getz and Akers, 1994, 1995). Furthermore, inhibitory interactions between olfactory sensory neurons had been suggested to sharpen and filter the sensory signal sent to the brain (Andersson et al., 2010; Couto et al., 2005). Our finding that worker antennae show higher expression of genes involved in neural modulation nicely corresponds with the many reports of context-dependent plasticity in peripheral sensory processing in worker honey bees (Bigot et al., 2012; Gadenne et al., 2016; Grosmaitre et al., 2001; McQuillan et al., 2012; Vergoz et al., 2009; Watanabe et al., 2014).

We are grateful to W. Roessler and J. Spaethe for constructive discussion regarding the project and experimental procedure during a research stay of R.J. at the University of Wuerzburg. We thank S. Unnikrishnan, W. Roessler, J. Spaethe, T. Tanimura and F. Marion-Poll for valuable and instructive comments on an earlier version of the manuscript.