- Research article
- Open Access
RNA-sequencing elucidates the regulation of behavioural transitions associated with the mating process in honey bee queens
BMC Genomicsvolume 16, Article number: 563 (2015)
Mating is a complex process, which is frequently associated with behavioural and physiological changes. However, understanding of the genetic underpinnings of these changes is limited. Honey bees are both a model system in behavioural genomics, and the dominant managed pollinator of human crops; consequently understanding the mating process has both pure and applied value. We used next-generation transcriptomics to probe changes in gene expression in the brains of honey bee queens, as they transition from virgin to mated reproductive status. In addition, we used CO2-narcosis, which induces oviposition without mating, to isolate the process of reproductive maturation.
The mating process produced significant changes in the expression of vision, chemo-reception, metabolic, and immune-related genes. Differential expression of these genes maps clearly onto known behavioural and physiological changes that occur during the transition from being a virgin queen to a newly-mated queen. A subset of these changes in gene expression were also detected in CO2-treated queens, as predicted from previous physiological studies. In addition, we compared our results to previous studies that used microarray techniques across a range of experimental time-points. Changes in expression of immune- and vision-related genes were common to all studies, supporting an involvement of these groups of genes in the mating process.
Our study is an important step in understanding the molecular mechanisms regulating post-mating behavioural transitions in a natural system. The weak overlap in patterns of gene expression with previous studies demonstrates the high sensitivity of genome-wide approaches. Thus, while we build on previous microarray studies that explored post-mating changes in honey bees, the broader experimental design, use of RNA-sequencing, and focus on Australian honey bees, which remain free from the devastating parasite Varroa destructor, in the current study, provide unique insights into the biology of the mating process in honey bees.
Mating is a key, and complex social behaviour, which is central to reproductive success across the animal kingdom. While mating behaviour has been the focus of numerous phenomenological studies (reviewed in ), there is growing interest in elucidating its molecular underpinnings in order to characterize its plasticity and evolution across animal taxa (e.g., [2, 3]). Brain transcriptomes can be used to associate complex behaviours like mating with patterns of gene expression at a genomic scale [2, 3], and to identify conserved molecular pathways across taxa [4, 5]. Such studies have revealed some of the specific molecules responsible for courtship or post-mating changes in behaviour in model organisms ranging from Drosophila melanogaster to voles [6–9]. However, the field of behavioural genomics is still in its infancy and further studies are needed to characterize the complexity and plasticity of mating behaviour in non-model organisms (see ) and to investigate the existence of genetic toolkits shared across taxa.
Social insects are likely to play a key role in developing our understanding of the molecular underpinnings of mating behaviour and post-mating changes. In social Hymenoptera (ants, some bees and wasps) queens only mate during a short period early in their life and undergo profound behavioural changes after mating, as they transform into nest-bound egg-laying machines. Mated queens have hugely extended lifespans relative to non-reproductive workers or similar solitary female insects [10–12]. Furthermore, because of the specialisation of reproductive and non-reproductive individuals in insect societies, the mating process may trigger sets of trade-offs that are different to those observed in non-social organisms. For example, the trade-off between reproduction and both immunity and longevity seen in solitary invertebrates [13, 14] may be uncoupled in social insects [10–12]. Consequently, studies of social insects may provide unique insights into the molecular mechanisms at play during and after the mating process. Thanks to modern genomic tools, such insights are emerging for ants [15, 16], wasps [17, 18], termites [19, 20], and honey bees [21, 22].
Honey bees (Apis mellifera) provide an excellent model system for molecular studies, as the behavioural ecology of mating is well understood. Queens mate on the wing when 6–13 days of age . Virgin queens leave their nest during the afternoon , and may fly several kilometres  to mate at leks known as drone congregation areas . At the lek they mate with multiple males in quick succession , before returning to their nest. About 50 % of queens commence oviposition after one mating flight, while the remainder fly on 1–2 subsequent afternoon(s) and again mate with multiple males [23, 27–29]. Physical stimulation of the bursa copulatrix during mating, as opposed to mating flights, contact with drones, or the presence of semen in the spermatheca, is responsible for activating oviposition . Mating is associated with profound behavioural  and physiological [32–34] changes that are unrelated to the age of the queen. In particular, virgin queens are photophilic and eager to fly , and are aggressive towards other virgin queens . In contrast, mated queens are photophobic, seek the protection of clustered workers and are unlikely to engage in fighting with other queens. Microarray studies have shown that these changes in behaviour and physiology are associated with profound differences in gene expression in the brain [22, 31].
Interestingly, similar physiological and behavioural changes to those found in mated queens are observed in queens subjected to double narcosis with carbon dioxide (CO2), providing an elegant tool for dissecting the molecular aspects of post-mating changes [36–40]. Insects possess specialized receptor cells that can detect and measure environmental CO2 [41, 42]. Honey bees use these receptors to help regulate CO2 concentrations in their nest by wing fanning [43–45]. Narcosis, caused by artificial CO2 exposure at much higher concentrations (50 % vs. a maximum of 3 % in the nest), usually triggers oviposition [46–48] by accelerating germ cell differentiation and stimulating the initial differentiation of the vitellarium . It has been hypothesized that CO2 might achieve these results by exploiting the carbonic anhydrase pathways, the hypoxia-induced transcription factors or via activation of transferrin . Narcosis can thus be used experimentally to separate the effects of insemination and mating flights from ovary activation on physiology and gene expression. Such studies have shown that treatment of honey bee queens with carbon dioxide stops mating flights, activates ovaries, alters the chemical profiles of mandibular glands, and affects gene expression both globally and for specific key genes such as dopamine receptors [36, 50–52].
Here we examine the effects of the mating process and CO2-narcosis on gene expression in the brains of eight day old queen honey bees from an Australian population using RNA-sequencing (RNAseq). We expected CO2-treated queens to exhibit gene expression patterns intermediate between virgin and naturally-mated queens. Our study is an important step forward in understanding the relationship between mating and gene expression in honey bees. First, it directly compares the neurogenomics of virgin, CO2-treated and naturally-mated honey bee queens within the same study, enabling us to dissect how the mating process, as opposed to just the stimulation of egg-laying, changes gene expression. Second, it is the first RNAseq study of such gene expression, and thus builds on and complements previous studies [21, 22, 31, 50] that have used the informative but less-powerful microarray approach .
In November 2012 we reared 20 sister queens of standard Australian commercial stock (primarily of Italian A.m. ligustica heritage) using standard beekeeping techniques from a single source colony . The day before the young queens were due to emerge from their pupal cells we transferred them to an incubator at 35 °C and emerged them in individual glass vials . At one day of age we paint-marked the queens and introduced each of them into their own nucleus colony. The entrance of each hive was covered with a queen excluder – a grid that allows workers to exit and enter, but confines the larger queens within the hive.
When the queens were 6 days old we randomly assigned the twenty queens to one of three treatments: 1) Mated queens (n = 7) were queens that successfully participated in a mating flight. Early afternoon of the 6th day, we removed the queen excluders, and then monitored the entrances until mating flights with successful mating (as indicated by a mating plug) were observed. 2) CO 2 -treated queens (n = 7) were subjected to 10 min of CO2-narcosis on day 6 and on day 7. The seven caged queens were placed in a zip-lock plastic bag, which was then flushed with compressed CO2 until all queens were completely immobile, and sealed for 10 min. Queens were released in their individual colonies after treatment on both days. 3) Virgin queens (n = 6) were caged along with the CO2-treated queens, but were returned to their respective colonies without narcosis.
All queens were harvested 2 days after treatment when they were 8 days old, at the same time of day (15:00 h), directly into liquid nitrogen and then stored at −80 °C. We used 4 queens per treatment group for RNAseq (total = 12 samples).
Dissections and RNA isolation
Queens were stored at −80 °C prior to dissection. Abdomens were dissected to examine ovary activation and the presence or absence of semen in the spermatheca. All mated queens had activated ovaries (defined as possessing developing eggs) and semen in their spermathecae, while no eggs or semen were detected in virgin and CO2-treated queens. Head capsules were removed and brains were dissected over dry ice . Brains were placed individually in Trizol and macerated to breakdown tissue. After addition of chloroform, samples were briefly vortexed and then centrifuged at 12,000 g for 15 min at 4 °C to generate an aqueous phase. After addition of isopropyl alcohol, the aqueous phase was again vortexed and centrifuged (10 min) to produce a pellet. The pellet was then washed with 75 % ethanol prior to drying and redissolving in ddH2O. The dissolved RNA was treated with DNAse I buffer (Ambion) to remove gDNA contaminations, prior to centrifugation and collection of the aqueous phase.
The RNA content of the sample and the purity of the extracts were assessed using a NanoDrop™ 1000 spectrophotometer (Thermo Scientific Inc., Bremen, Germany).
We sequenced RNA with an Illumina HiSeq system using 2 lanes of a plate (6 samples per lane) and producing 50 bp single-end reads. Reads were checked with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) for quality control and were subsequently processed with Trimmomatic  to remove adapters and low quality bases: reads less than 36 bases long were filtered out (Additional file 1). We also dropped reads that matched ribosomal RNA sequences (rRNA) by means of SortMeRna . Surviving reads were aligned with TopHat for Illumina  to the latest version of the honey bee genome (Apis mellifera Assembly 4.5) available on BeeBase (http://hymenopteragenome.org/beebase/?q=download_sequences) by using the Galaxy web-based platform (https://usegalaxy.org/). Mapped reads were converted into raw read counts with SAMtools idxstats  and these were used to quantify differential gene expression.
Analysis of gene expression
Raw read counts were analysed with R using the edgeR package from Bioconductor . Read counts were log2 transformed to correct for the skew to zero and large values. Only genes with at least 10 reads per sample were kept in the analysis (12,992 genes). Normalization was performed with Trimmed Mean of M-values (TMM), a method that is implemented in the edgeR Bioconductor package . We detected differential levels of gene expression using a modified Fisher’s exact test that takes into account both dispersion and multiple samples. Finally, raw P-values for each gene were corrected for multiple comparisons setting a false discovery rate (FDR) of 5 %.
For global analyses of gene expression we used hierarchical clustering (Ward method) and principal component analysis in JMP Pro 10.0 (SAS, Cary, NC). To perform Gene Ontology (GO) analyses we obtained Drosophila melanogaster orthologs with BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) for honey bee genes that were significantly differentially expressed between treatments and computed functional annotation clustering in DAVID version 6 [62, 63] with medium stringency and a cutoff of P-value < 0.05. To identify overrepresented biological functions (enrichment analysis) we compared the annotation composition in our list of differentially expressed genes to that of a population background composed of all the honey bee genes with Drosophila orthologs.
We used Venny (http://bioinfogp.cnb.csic.es/tools/venny/index.html) to overlap D. melanogaster ortholog matching lists of significantly differentially expressed genes and a Hypergeometric test (http://nemates.org/MA/progs/overlap_stats.html) to assess whether genes overlapping between studies occurred significantly more often than expected by chance. We compared our study to three microarray studies [22, 31, 50] and one qPCR study  on honey bee gene expression after either mating or CO2-narcosis (for more details see Results and Additional files 2, 3, 4, 5, 6, 7). For comparisons with microarray studies, we first overlapped whole datasets from each study to identify potential candidate genes for mating behaviour and response to carbon dioxide and performed GO analysis on them. We also performed overlap analyses on genes from pairwise comparisons of interest to evaluate the level of agreements across studies.
The sequencing produced 30 million reads per sample on average (min 22,298,993 and max 39,301,839, see Additional file 1). About 7 % of total reads were dropped during the filtering steps and a further 4 % were excluded as ribosomal RNA. On average, 98 % of surviving reads per sample were aligned to single locations in the honey bee genome.
Global analysis of gene expression
A total of 1088 genes were significantly differentially expressed (FDR < 0.05) in at least one of three pairwise comparisons across treatments (see below); this represents 7.10 % of the 15,314 coding sequences present in the honey bee genome. We performed a hierarchical clustering analysis on this set of genes to identify common patterns of gene expression across individual samples. We obtained a clear separation of our samples into two groups: mated queens clustered alone while virgin queens and CO2-treated clustered together forming a separate macro-group (Figure 1). This suggests that the mating treatment was the major driver of gene expression in the honey bee brain.
We also performed a principal component analysis on the same set of 1088 significantly differentially expressed genes (Fig. 2). Again, most of the changes in gene expression (77 %) were associated with the mating process: mated queens clustered on the opposite side of the chart compared to virgin and CO2-treated queens. Only 23 % of the difference in gene expression was associated with CO2-treatment as compared to virgin queens.
Pairwise comparisons between treatments
Mated vs. virgin queens
Pairwise analysis of gene expression between treatments revealed a large set of genes that were significantly differentially expressed between mated and virgin queens: 829 genes at FDR < 0.05 (Fig. 3 and Additional file 8). Of these, 654 genes were up-regulated and 175 were down-regulated in mated queens. Our GO analysis was based on the 475 D. melanogaster orthologs for these genes. Six metabolic pathways and 28 GO terms were significantly overrepresented (Table 1 and Additional file 9). Of particular interest were GO terms associated with the following key biological processes: sensory perception (P-value = 0.011, 15 genes, up-regulated in mated = 5, down-regulated = 10), detection of stimulus (P-value = 0.040, 8 genes, up-regulated in mated = 2, down-regulated = 6), multiple metabolism-related GO terms such as the fatty acid metabolic processes (P-value =0.015, 7 genes, all up-regulated in mated) and the immune-related GO terms defense response (P-value = 0.009, 11 genes, up-regulated in mated = 10, down-regulated = 1), melanisation defense response (P-value = 0.012, 3 genes, all up-regulated in mated) and innate immune response (P-value = 0.015, 7 genes, up-regulated in mated = 6, down-regulated = 1). See Additional file 9 for a list of the genes associated with these GO terms.
Mated vs. CO2-treated queens
There were 629 significantly differentially expressed (FDR < 0.05) genes between mated and CO2-treated queens (Fig. 3 and Additional file 10): 409 genes were up-regulated and 220 were down-regulated in mated queens. GO analysis of the matching 310 D. melanogaster orthologs revealed that 13 GO terms were significantly overrepresented (Table 1 and Additional file 11), including sensory perception (P-value = 0.024, 11 genes, up-regulated in mated = 4, down-regulated = 7) and response to organic substance (P-value = 0.044, 8 genes, up-regulated in mated = 5, down-regulated =3). See Additional file 11 for a list of the genes associated with these GO terms.
CO2-treated vs. virgin queens
A much smaller set of 151 genes were significantly differentially expressed (FDR < 0.05) between CO2-treated and virgin queens (Fig. 3 and Additional file 12): 117 genes were up-regulated and 34 were down-regulated in CO2-treated queens. The 59 D. melanogaster orthologs for these genes produced 6 significantly overrepresented GO terms (Table 1 and Additional file 13), including the biological process cognition (P-value = 0.004, 7 genes, up-regulated in CO2-treated = 1, down-regulated =6). See Additional file 13 for a list of the genes associated with this GO term.
In order to compare our RNAseq analysis of the honey bee mating process and CO2-treatment to previous research addressing similar questions with the microarray technique, we overlapped lists of significantly differentially expressed genes from our study and the three following studies: 1) Kocher et al. , where the authors examined brain gene expression in virgin, mated and egg-laying honey bee queens, providing obvious potential comparisons with our mated queens; 2) Kocher et al. , where they examined the effects of mating and instrumental insemination with saline or semen on gene expression in the brains of honey bee queens, providing useful comparisons with both our mated and CO2-treated queens (as instrumental insemination involves CO2-treatment); 3) Niño et al. , where the brain transcriptomic profile was evaluated in virgin, CO2-treated, and physically manipulated honey bee queens (i.e., exposed to CO2 and sham-inseminated), providing valuable comparisons for our CO2-treated queens.
Very few genes were shared across studies for the focal pairwise comparisons, resulting in a lack of statistically significant overlap (see Figs. 4 and 5 and Additional files 2, 3, 4, 5, 6, 7 for details on these analyses). In contrast, at the whole dataset level, a significant number of differentially expressed genes were shared between our study and Kocher et al. . However, no significant overlap was found in comparisons with data from Kocher et al.  and Niño et al.  (see Additional files 2 and 7). At the coarser GO level of analysis, immune-related GO terms and the GO term ‘response to other organism’ were recurrent across all queen studies (Kocher et al. [22, 31], Niño et al.  and this study).
The mating process and CO2-induced narcosis significantly affect the neurogenomic profile of honey bee queens, by changing the expression of more than a thousand genes within two days of treatment. The majority of these differences were driven by the mating process, maturation than natural mating. Importantly, changes in pathways related to vision, metabolism, chemoreception, and immunity match expectations from known behavioural changes. These next-generation sequencing results partially corroborate previous microarray studies, and provide new insight into the molecular regulation of key behavioural transitions in honey bee queens.
Neurogenomics of mating behaviour
Gene expression in the brains of mated queens differs strongly from that seen in virgins. Strikingly, genes that are associated with vision such as Rhodopsin 2 (Rh2), neither inactivation nor afterpotential A (ninA), Arrestin 2 (Arr2), G protein beta-subunit 76C (Gβ76C), chaoptin (chp) and Calhotin (Cpn) were all down-regulated in mated queens compared to virgins (see Additional file 6). Changes in the expression of these visual perception genes mirrors the transition from photophilic behaviour observed in virgin queens that engage in mating flights, to more photophobic behaviour in mated queens confined within the nest. Queens are required to fly during swarming events , and it would be interesting to see if the vision system is reactivated in queens as they prepare to swarm. It would also be interesting to determine whether queens of open nesting honey bee species like Apis florea, in which the queens are able and ready to fly at all times , show the same decline in vision-related genes after mating.
Another group of genes that are differentially expressed in mated and virgin queens belong to the family of odorant receptors (ORs) and odorant binding proteins (OBPs). Odour and pheromone perception are central to the social life of honey bees [67–69] and play an important role in mate location and efficient mating. Of particular interest is the presence of the Pheromone-binding protein-related protein 2 (Obp19d), which was down-regulated in mated queens. This gene is involved in the detection of pheromones, which are chemical compounds released by one individual to trigger a social response in members of the same species and therefore are widely used among social insects to mediate mating behaviour. This protein may be used by virgin honey bee queens to locate rival virgins. Interestingly, of ten genes encoding predicted pheromone or odorant binding proteins that are affected by mating in Drosophila  only four were also significantly differentially expressed in our study (Obp19d, Obp56g, Or43a and Or49b).
As previously reported [22, 31, 50, 71], genes broadly related to metabolism were also differentially expressed between mated and virgin queens (see Additional file 9). Mating flights are energetically expensive  whereas nest-bound life is not. It is therefore unsurprising that genes related to carbohydrate metabolism, glycolysis and gluconeogenesis were differentially regulated in mated queens relative to virgin queens. In contrast, the requirement for fatty acid metabolism increases after mating, because novel lipids are required for pheromone synthesis and egg production . As expected, all seven genes in this cluster were up-regulated in mated queens.
Finally, the last important group of genes that differ between mated and virgin queens is the immune genes. Three of these genes, Serine Protease Immune Response Integrator (spirit), Peptidoglycan recognition protein SA (PGRP-SA) and Gram-negative bacteria binding protein 1 (GNBP1) are important players in the Toll pathway and their expression is usually triggered by a challenge from Gram-positive bacteria . Three additional genes, Melanization Protease 1 (MP1), Serpin-27A (Spn27A) and Hemolectin (Hml) are involved in the melanisation response [75, 76] and can therefore be triggered by bacterial infection or wound healing reactions. Finally, two other genes in this group, the antimicrobial peptides defensin and serpent, are involved in hematopoiesis . With the exception of defensin, all immune genes were up-regulated in mated queens and this is likely to result in higher immunocompetence as more defense molecules, such as antimicrobial peptides, are produced, and cellular responses or wound healing reactions may be more effective. Increased expression of immune genes post mating has been observed repeatedly in honey bees [22, 31, 50, 71] and other organisms . In Drosophila, for example, 19 immune-related genes respond to mating, including defensin and PGRP-SA . One explanation for up-regulation of immune genes in females after mating might be activation of the immune system by immune elicitors associated with the male reproductive apparatus or the ejaculate, or by a traumatic insemination event that leads to wound healing reactions. However, this does not seem to be the case in the honey bee, where it has been previously observed that non-traumatic instrumental insemination with sterile saline is sufficient to up-regulate the expression of immune genes . Consequently, up-regulation of immune genes may be triggered, not by the physical process of mating, but through molecular cross-talk between reproductive and immune gene-expression pathways. Such cross-talk may have been selected because of adaptive advantages to the queen and to the colony, such as increased protection from horizontally transmitted parasites and pathogens . In contrast, in solitary species immunocompetence is reduced after mating, so that energy can be redirected from the immune system to reproductive activity, i.e., sperm storage or egg-laying [80, 81]. Our results demonstrate that this type of trade-off is not present in honey bees, perhaps because the newly-mated queen can rely on the social immunity conferred by her natal colony, and because queens are fed a near-perfect diet of worker mandibular secretions [82, 83]. It would be interesting to conduct a comparative study across solitary- and swarm-founding social insects, to test this hypothesis. Another explanation for the observed pattern could be that virgin queens reduce their investment in immune defence in preparation for mating flights: this would allow them to allocate a greater portion of their energy resources in two activities (flying and mating) that are energetically expensive. A time-course study on the levels of expression of immune genes in honey bee queens from emergence to full reproduction would address this question.
In addition to the GO analyses, we examined the expression of several genes that might play a role in mating and ovary activation. Genes of the insulin/insulin-like signalling (IIS) and TOR pathways are responsible for regulating growth and nutrition and are fundamental to the process of honey bee caste determination . We found that only one gene in these pathways was significantly differentially expressed between mated and virgin queens, Phosphoinositide-dependent kinase 1 (PDK1). This gene was up-regulated in mated queens. The fact that PDK1 was up-regulated in the only group of queens with activated ovaries in our studies is in line with the hypothesis that PDK1 activity is linked to ovary size in worker bees, where foragers with a bias toward pollen have both larger ovaries and higher levels of expression of PDK1 compared to foragers that are nectar-biased .
We also examined expression of the biogenic amines, as these compounds may be involved in mediating the interactions between brain and ovaries during reproductive activation. Dopamine signalling pathways are positively associated with reproductive status in workers [37, 40, 86]. For example, the gene N-acetyldopamine is positively correlated with ovarian development . In queens, instead, dopamine and reproductive status are negatively correlated and dopamine levels decrease after mating . Our study provides further support for this reversed relationship in queens, as dopamine N-acetyltransferase (Dat), a component of the catabolism of dopamine, was down-regulated in mated queens.
Neurogenomics of CO2-narcosis
The most important differences in gene expression profiles between CO2-treated and virgin queens relate to cognition (see Additional file 13 for other differentially expressed genes). Most of the genes in this group were the same genes found in the GO terms “sensory perception” and “detection of stimulus” for the mated vs. virgin queen contrast; in addition we found neither inactivation nor afterpotential C (ninaC) and no receptor potential A (norpA, see Additional file 6). All cognition genes but one were down-regulated in CO2-treated queens, as they were in mated queens. The similarity in patterns of expression of these genes suggests that CO2-treated queens undergo a process of de-activation of visual perception and eye development genes similar to that seen in mated queens. One potential explanation for this is that CO2-narcosis induces acidosis, as may mating flights , indicating a potential role for body pH as a trigger for these changes in gene expression. However, this must be separate from the processes that link CO2-treatment to oviposition, as Koeniger et al.  have demonstrated that flight is not sufficient to induce oviposition in honey bee queens. This suggests that different aspects of the mating process may trigger distinct sets of changes in gene expression.
Despite such similarities in gene expression in mated and CO2-treated queens when compared to virgins, a direct comparison of the two groups also highlighted interesting differences. Expression of genes related to sensory perception differed strongly across the two groups, as did genes associated with olfactory and gustatory activity. Five ORs (Or13a, Or43a, Or49a, Or85b and Or85c) and Obp19d were up-regulated in CO2-treated queens, while three other OBPs (Obp56d, Obp56g and Obp83b) were up-regulated in mated queens. CO2-treated and mated queens also differ for genes related to the response to organic substances: among these, of particular relevance is the Insulin-like receptor (InR), which was down-regulated in CO2-treated queens. InR is a major player in the insulin signalling pathway and regulates important biological functions such as metabolism, growth, reproduction, and aging . Again, this supports the idea that CO2-narcosis induces only a subset of the changes in gene expression caused by the process of mating and onset of oviposition. In addition, the two-day post-treatment time interval prior to sampling did not allow for complete development of eggs in CO2-treated queens, while it was sufficient for mated queens, as shown by our dissections. This is not surprising, as egg development has been observed after twelve days from narcosis or instrumental insemination [50, 71], which is a much longer time span. Together with our gene expression results, this confirms that mating and CO2-narcosis, despite producing the same final result, follow apparently different pathways of action.
Only two of the three whole dataset comparisons with previous studies produced statistically significant overlaps. However, there is a common pattern in the biological functions that are associated with genes that are shared between studies. GO terms related to immune functions and response to other organism are overrepresented across studies [22, 31, 50]. This makes us confident that these processes are key to the physiological and behavioural changes that take place as a queen transitions from her initial virgin state to that of a mated matriarch.
Interestingly, while the findings of Kocher et al.  are significantly congruent with our study, those from Kocher et al.  and Niño et al.  are not, no doubt as a consequence of the similarities and differences across the studies in experimental treatments. In our study and in Kocher et al. , queens were collected two days after treatment; in contrast, Kocher et al.  and Niño et al.  analysed samples collected 5 and 10 days after treatment, respectively. This suggests both that the neurogenomic state of an individual is repeatable across studies, and that it is temporally dynamic. Future studies are needed to understand how and why gene expression changes over time, both before and after mating.
Finally, a number of other factors may explain the low-level of congruence between our results and earlier studies. An obvious methodological factor is that, while previous studies used microarrays, we used RNAseq to generate transcriptomic data. RNAseq is a more powerful technique than microarrays, and relies on a different experimental/statistical approach (for a comparison between microarray and RNAseq platforms see Guo et al. ), which may hinder comparisons among studies. From a biological perspective, in all studies experimental bees were produced by a single colony. This approach is frequently used in studies of social insects, as it eliminates variance related to inter-colonial variability, resulting in higher statistical power to detect inter-individual differences in gene expression as a result of treatments applied to individual bees. However, such benefits are lost when comparing across studies, as each colony is a unit on its own (the super-organism) characterized by a particular social environment.
This study is an important advance in the molecular characterization of the mating process in the honey bee Apis mellifera, a model system for sociobiology and the most important managed crop pollinator. Our results are partially in line with previous studies, but demonstrate interesting mismatches, denoting the importance of consistent experimental design in genomic studies. By uncovering many important biological functions associated with the mating process, this study also stresses once more the complexity of this behaviour and the need in the future to isolate the single components (e.g., mating flights and copulation) and analyse them separately. In addition, this study focused on Australian honey bees, which comprise one of the last populations free of the invasive parasitic mite Varroa destructor [91–93]. As such, our results provide a biological baseline that can be used as a reference to understand the impact of external challenges on bee decline in the Western world . At the same time, this important difference in parasite load between populations of honey bees could underlie the differences in patterns of gene expression revealed by our comparative analyses. Future studies need to focus on determining the mechanisms behind cross-study variation to isolate key and repeatable differences in gene expression during the mating process in honey bees, as well as the impact of important parasites and pathogens on this process.
Availability of supporting data
RNAseq raw sequence reads and normalized expression values for each gene are available in the NCBI Gene Expression Omnibus repository, Series record GSE65833 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=oludywiunlwblgd&acc=GSE65833. Other data sets supporting the results of this article are included within the Additional files.
For this study honey bees from commercial colonies were used and therefore no ethics statement is required.
- CO2 :
High-throughput RNA sequencing
genomic Deoxyribonucleic acid
False discovery rate
quantitative real-time polymerase chain reaction
Odorant binding protein
Insulin/insulin-like signaling pathway
Target of rapamycin pathway
Kyoto encyclopedia of genes and genomes
Shuster SM, Wade MJ: Mating systems and strategies: Princeton, New Jersey: Princeton University Press; 2003.
Hitzemann R, Bottomly D, Darakjian P, Walter N, Iancu O, Searles R, et al. Genes, behavior and next-generation RNA sequencing. Genes Brain Behav. 2013;12(1):1–12.
Harris RM, Hofmann HA. Neurogenomics of Behavioral Plasticity. In: Ecological Genomics. Springer; 2014. p. 149–68.
O'Connell LA, Hofmann HA. Genes, hormones, and circuits: An integrative approach to study the evolution of social behavior. Front Neuroendocrinol. 2011;32(3):320–35.
Dolezal AG, Toth AL. Honey bee sociogenomics: a genome-scale perspective on bee social behavior and health. Apidologie. 2014;45(3):375–95.
Robinson GE, Fernald RD, Clayton DF. Genes and social behavior. Science. 2008;322(5903):896–900.
Ram KR, Wolfner MF. A network of interactions among seminal proteins underlies the long-term postmating response in Drosophila. Proc Natl Acad Sci. 2009;106(36):15384–9.
Hammock EA, Young LJ. Microsatellite instability generates diversity in brain and sociobehavioral traits. Science. 2005;308(5728):1630–4.
Lim MM, Wang Z, Olazábal DE, Ren X, Terwilliger EF, Young LJ. Enhanced partner preference in a promiscuous species by manipulating the expression of a single gene. Nature. 2004;429(6993):754–7.
Castella G, Christe P, Chapuisat M. Mating triggers dynamic immune regulations in wood ant queens. J Evol Biol. 2009;22(3):564–70.
Page REJ, Peng CY-S. Aging and development in social insects with emphasis on the honey bee, Apis mellifera L. Exp Gerontol. 2001;36(4):695–711.
Heinze J, Schrempf A. Aging and reproduction in social insects–a mini-review. Gerontology. 2008;54(3):160–7.
Lawniczak MK, Barnes AI, Linklater JR, Boone JM, Wigby S, Chapman T. Mating and immunity in invertebrates. Trends Ecol Evol. 2007;22(1):48–55.
Morrow EH, Innocenti P. Female postmating immune responses, immune system evolution and immunogenic males. Biol Rev. 2012;87(3):631–8.
Feldmeyer B, Elsner D, Foitzik S. Gene expression patterns associated with caste and reproductive status in ants: worker‐specific genes are more derived than queen‐specific ones. Mol Ecol. 2014;23(1):151–61.
Wurm Y, Wang J, Keller L. Changes in reproductive roles are associated with changes in gene expression in fire ant queens. Mol Ecol. 2010;19(6):1200–11.
Ferreira PG, Patalano S, Chauhan R, Ffrench-Constant R, Gabaldón T, Guigó R, et al. Transcriptome analyses of primitively eusocial wasps reveal novel insights into the evolution of sociality and the origin of alternative phenotypes. Genome Biol. 2013;14(2):R20.
Toth AL, Varala K, Henshaw MT, Rodriguez-Zas SL, Hudson ME, Robinson GE. Brain transcriptomic analysis in paper wasps identifies genes associated with behaviour across social insect lineages. Proc R Soc B Biol Sci. 2010;277(1691):2139–48.
Korb J, Weil T, Hoffmann K, Foster KR, Rehli M. A gene necessary for reproductive suppression in termites. Science. 2009;324(5928):758–8.
Steller MM, Kambhampati S, Caragea D. Comparative analysis of expressed sequence tags from three castes and two life stages of the termite Reticulitermes flavipes. BMC Genomics. 2010;11(1):463.
Grozinger CM, Fan Y, Hoover SER, Winston ML. Genome-wide analysis reveals differences in brain gene expression patterns associated with caste and reproductive status in honey bees Apis mellifera. Mol Ecol. 2007;16:4837–48.
Kocher SD, Richard F-J, Tarpy DR, Grozinger CM. Genomic analysis of post-mating changes in the honey bee queen (Apis mellifera). BMC Genomics. 2008;9(1):232.
Oertel E. Mating flights of queen bees. Gleanings in Bee Culture. 1940;68:292–3.
Neumann P, van Praagh JP, Moritz RFA, Dustmann JH. Testing reliability of a potential island mating apiary using DNA microsatellites. Apidologie. 1999;30(4):257–76.
Loper GM, Wolf WW, Taylor OR. Honey bee drone flyways and congregation areas - radar observations. J Kansas Entomol Soc. 1992;65(3):223–30.
Koeniger G, Koeniger N, Fabritius M. Some detailed observations of mating in the honeybee. Bee Wld. 1979;60:53–7.
Tarpy DR, Page RE. No behavioral control over mating frequency in queen honey bees (Apis mellifera L.): Implications for the evolution of extreme polyandry. Am Nat. 2000;155(6):820–7.
Roberts WC. Multiple mating of queen bees proved by progeny and flight tests. Gleanings in Bee Culture. 1944;72(6):225–59. 303.
Woyke J. Causes of repeated mating flights by queen honeybees. J Apic Res. 1964;3(1):17–23.
Koeniger G. In welchem Abschnitt des Paarungsverhaltens der Bienenkönigin findet die Induktion der Eiablage statt? Apidologie. 1981;12(4):329–43.
Kocher S, Tarpy D, Grozinger C. The effects of mating and instrumental insemination on queen honey bee flight behaviour and gene expression. Insect Mol Biol. 2010;19(2):153–62.
Fahrbach SE, Giray T, Robinson GE. Volume changes in the mushroom bodies of adult honey bee queens. Neurobiol Learn Mem. 1995;63:181–91.
Kocher SD, Richard F-J, Tarpy DR, Grozinger CM. Queen reproductive state modulates pheromone production and queen-worker interactions in honeybees. Behav Ecol. 2009;20(5):1007–14.
Richard F-J, Tarpy DR, Grozinger CM. Effects of insemination quantity on honey bee queen physiology. PLoS One. 2007;2, e980.
Pflugfelder J, Koeniger N. Fight between virgin queens (Apis mellifera) is initiated by contact with the dorsal abdominal surface. Apidologie. 2003;34(3):249–56.
Vergoz V, Lim J, Duncan M, Cabanes G, Oldroyd BP. Effects of natural mating and CO2 narcosis on biogenic amine receptor gene expression in the ovaries and brain of queen honey bees, Apis mellifera. Insect Mol Biol. 2012;21(6):558–67.
Koywiwattrakul P, Thompson GJ, Sitthipraneed S, Oldroyd BP, Maleszka R. Effects of carbon dioxide narcosis on ovary activation and gene expression in worker honeybees, Apis mellifera. J Insect Sci. 2005;5.
Berger B, Abdalla FC, Cruz-Landim C. Effect of narcosis with CO2 on the ovarian development in queens of Apis mellifera (Hymenoptera, Apini). Socbiol. 2005;261–270.
Ebadi R, Gary NE, Lorenzen K. Effects of carbon dioxide and low temperature narcosis on honey bees, Apis mellifera. Environmental Entomology. 1980;9(1):144–8.
Harris JW, Woodring J, Harbo JR. Effects of carbon dioxide on levels of biogenic amines in the brains of queenless worker and virgin queen honey bees (Apis mellifera). J Apic Res. 1996;35(2):69–78.
Bierbower SM, Cooper RL. The effects of acute carbon dioxide on behavior and physiology in Procambarus clarkii. J Exp Zool A Ecol Genet Physiol. 2010;313A(8):484–97.
Robertson HM, Kent LB. Evolution of the gene lineage encoding the carbon dioxide receptor in insects. J Insect Sci. 2009;9.
Sudarsan R, Thompson C, Kevan PG, Eberl HJ. Flow currents and ventilation in Langstroth beehives due to brood thermoregulation efforts of honeybees. J Theor Biol. 2012;295:168–93.
Southwick EE, Moritz RFA. Social control of air ventilation in colonies of honey bees, Apis mellifera. Journal of Insect Physiology. 1987;33(9):623–6.
Seeley TD. Atmospheric carbon dioxide regulation in honey-bee (Apis mellifera) colonies. J Insect Physiol. 1974;20(11):2301–5.
Kaftanoglu O, Peng YS. Effects of insemination on the initiation of oviposition in the queen honeybee. J Apic Res. 1982;21(1):3–6.
Engels W, Gonçalves L, Engels E. Effects of carbon dioxide on vitellogenin metabolism in unmated queen honeybees. J Apic Res. 1976;15:3–10.
Engels W, Ramamurty P. Initiation of oogenesis in allatectomised virgin honey bee queens by carbon dioxide treatment. J Insect Physiol. 1976;22(10):1427–32.
Berger B, Abdalla F, Cruz-Landim C. Effect of narcosis with CO2 on the ovarian development in queens of Apis mellifera (Hymenoptera, Apini). Socbiol. 2005;45(2):261–70.
Niño E, Tarpy D, Grozinger C. Genome‐wide analysis of brain transcriptional changes in honey bee (Apis mellifera L.) queens exposed to carbon dioxide and physical manipulation. Insect Mol Biol. 2011;20(3):387–98.
Niño EL, Malka O, Hefetz A, Tarpy DR, Grozinger CM. Chemical profiles of two pheromone glands are differentially regulated by distinct mating factors in honey bee queens (Apis mellifera l.). Plos One. 2013;8(11).
Thompson GJ, Yockey H, Lim J, Oldroyd BP. Experimental manipulation of ovary activation and gene expression in honey bee (Apis mellifera) queens and workers: testing hypotheses of reproductive regulation. J Exp Zool A Ecol Genet Physiol. 2007;307(10):600–10.
Guo Y, Sheng Q, Li J, Ye F, Samuels DC, Shyr Y. Large scale comparison of gene expression levels by microarrays and RNAseq using TCGA data. PLoS One. 2013;8(8), e71462.
Harbo JR. Propagation and instrumental insemination. In: Rinderer TE, editor. In: Bee Genetics and Breeding. Orlando: Academic; 1986. p. 361–89.
Wagener-Hulme C, Kuehn J, Schulz D, Robinson G. Biogenic amines and division of labor in honey bee colonies. J Comp Physiol A. 1999;184(5):471–9.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;170.
Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28(24):3211–7.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.
Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Brito R, McHale M, Oldroyd B. Expression of genes related to reproduction and pollen foraging in honey bees (Apis mellifera) narcotized with carbon dioxide. Insect Mol Biol. 2010;19(4):451–61.
Winston ML: The biology of the honey bee: Cambridge, MA, USA: Harvard University Press; 1991.
Oldroyd BP, Wongsiri S: Asian honey bees: biology, conservation, and human interactions: Cambridge, MA, USA: Harvard University Press; 2009.
Forêt S, Maleszka R. Function and evolution of a gene family encoding odorant binding-like proteins in a social insect, the honey bee (Apis mellifera). Genome Res. 2006;16(11):1404–13.
Gary NE. Chemical mating attractants in the queen honey bee. Science. 1962;136(3518):773–4.
Free JB: Pheromones of social bees: London, UK: Chapman and Hall; 1987.
McGraw LA, Gibson G, Clark AG, Wolfner MF. Genes regulated by mating, sperm, or seminal proteins in mated female Drosophila melanogaster. Curr Biol. 2004;14(16):1509–14.
Niño E, Tarpy D, Grozinger C. Differential effects of insemination volume and substance on reproductive changes in honey bee queens (Apis mellifera L.). Insect Mol Biol. 2013;22(3):233–44.
Kunieda T, Fujiyuki T, Kucharski R, Foret S, Ament SA, Toth AL, et al. Carbohydrate metabolism genes and pathways in insects: insights from the honey bee genome. Insect Mol Biol. 2006;15(5):563–76.
Teerawanichpan P, Robertson AJ, Qiu X. A fatty acyl-CoA reductase highly expressed in the head of honey bee (Apis mellifera) involves biosynthesis of a wide range of aliphatic fatty alcohols. Insect Biochem Mol Biol. 2010;40(9):641–9.
Kambris Z, Brun S, Jang I-H, Nam H-J, Romeo Y, Takahashi K, et al. Drosophila Immunity: A Large-Scale In Vivo RNAi Screen Identifies Five Serine Proteases Required for Toll Activation. Curr Biol. 2006;16(8):808–13.
Tang H, Kambris Z, Lemaitre B, Hashimoto C. Two Proteases Defining a Melanization Cascade in the Immune System of Drosophila. J Biol Chem. 2006;281(38):28097–104.
Lesch C, Goto A, Lindgren M, Bidla G, Dushay MS, Theopold U. A role for Hemolectin in coagulation and immunity in Drosophila melanogaster. Dev Comp Immunol. 2007;31(12):1255–63.
Rehorn K-P, Thelen H, Michelson AM, Reuter R. A molecular aspect of hematopoiesis and endoderm development common to vertebrates and Drosophila. Development. 1996;122(12):4023–31.
Rolff J, Siva-Jothy MT. Copulation corrupts immunity: a mechanism for a cost of mating in insects. Proc Natl Acad Sci. 2002;99(15):9916–8.
Schmid-Hempel P: Parasites in social insects: Princeton, New Jersey, USA: Princeton University Press; 1998.
Baer B, Armitage SA, Boomsma JJ. Sperm storage induces an immunity cost in ants. Nature. 2006;441(7095):872–5.
Siva‐Jothy MT, Tsubaki Y, Hooper RE. Decreased immune response as a proximate cost of copulation and oviposition in a damselfly. Physiol Entomol. 1998;23(3):274–7.
Crailsheim K. Interadult feeding of jelly in honeybee (Apis mellifera L.) colonies. J Comp Physiol B. 1991;161(1):55–60.
Crailsheim K. Trophallactic interactions in the adult honeybee (Apis mellifera, L). Apidologie. 1998;29:97–112.
Wheeler DE, Buck NA, Evans JD. Expression of insulin/insulin-like signalling and TOR pathway genes in honey bee caste determination. Insect Mol Biol. 2014;23(1):113–21.
Wang Y, Amdam GV, Rueppell O, Wallrichs MA, Fondrk MK, Kaftanoglu O, et al. PDK1 and HR46 gene homologs tie social behavior to ovary signals. PLoS One. 2009;4(4):e4899.
Vergoz V, Lim J, Oldroyd B. Biogenic amine receptor gene expression in the ovarian tissue of the honey bee Apis mellifera. Insect Mol Biol. 2012;21(1):21–9.
Sasaki K, Nagao T. Distribution and levels of dopamine and its metabolites in brains of reproductive workers in honeybees. J Insect Physiol. 2001;47(10):1205–16.
K-i H, Sasaki K, Nagao T. Depression of brain dopamine and its metabolite after mating in European honeybee (Apis mellifera) queens. Naturwissenschaften. 2005;92(7):310–3.
Harrison JF. Insect acid–base physiology. Annu Rev Entomol. 2001;46(1):221–50.
Wu Q, Brown MR. Signaling and function of insulin-like peptides in insects. Annu Rev Entomol. 2006;51(1):1–24.
Navajas M, Migeon A, Alaux C, Martin-Magniette M-L, Robinson G, Evans J, et al. Differential gene expression of the honey bee Apis mellifera associated with Varroa destructor infection. BMC Genomics. 2008;9(1):301.
Aronstein KA, Saldivar E, Vega R, Westmiller S, Douglas AE. How Varroa parasitism affects the immunological and nutritional status of the honey bee Apis mellifera. Insects. 2012;3(3):601–15.
Alaux C, Dantec C, Parrinello H, Le Conte Y. Nutrigenomics in honey bees: digital gene expression analysis of pollen's nutritive effects on healthy and Varroa-parasitized bees. BMC Genomics. 2011;12(1):496.
Oldroyd BP. What's killing American honey bees? PLoS Biol. 2007;5(6), e168.
This study was funded by a BBSRC ISIS grant BB/J019453/1, a Royal Holloway Research Strategy Fund Grant, and a Leverhulme Grant F/07537/AK to MJFB. BPO was supported by Australian Research Council Discovery grants DP150100151 and DP120101915. FM was supported by a Marie Curie International Incoming Fellowship FP7-PEOPLE-2013-IIF-625487 to MJFB. We would like to thank Dave Galbraight (Penn State) and Alberto Paccanaro (RHUL) for support with analysis of RNAseq data and four anonymous reviewers for providing thoughtful insights that helped to improve the manuscript.
The authors declare that they have no competing interests.
The study was conceived by MJFB, in collaboration with BPO. BPO provided and prepared the honey bee queen samples. MJFB and VV conducted the RNA extraction. FM conducted the analyses of transcriptomic data and prepared the initial draft of the manuscript. All authors contributed to writing the manuscript and approved the final version.
Fabio Manfredini and Mark J F Brown contributed equally to this work.
Summary statistics for Illumina HiSeq, read preparation and alignment to the honey bee genome.
Comparative studies. Overlap analyses with previous microarray studies on the effect of mating and CO2-narcosis on gene expression in honey bees.
Overlap analysis with Kocher et al. [ 22]. Genes that were significantly differentially expressed between mated and virgin queens in both studies.
Overlap analysis with Kocher et al. [ 31]. Genes that were significantly differentially expressed between mated and virgin queens in both studies.
Overlap analysis with Niño et al. [ 50]. Genes that were significantly differentially expressed between CO2-treated and virgin queens in both studies.
Genes involved in eye development and visual perception that were differentially expressed between treatments in this study.
Results of comparative studies.
Pairwise comparisons. List of the genes that were differentially expressed between mated and virgin queens at FDR < 0.05.
Gene ontology analysis for the pairwise comparison mated vs. virgin. Significantly enriched GO terms and KEGG pathways (functional annotation chart, P-value < 0.05) and details on associated genes for terms of interest.
Pairwise comparisons. List of the genes that were differentially expressed between CO2-treated and mated queens at FDR < 0.05.
Gene ontology analysis for the pairwise comparison CO 2 -treated vs. mated. Significantly enriched GO terms (functional annotation chart, P-value < 0.05) and details on associated genes for terms of interest.
Pairwise comparisons. List of the genes that were differentially expressed between CO2-treated and virgin queens at FDR < 0.05.
Gene ontology analysis for the pairwise comparison CO 2 -treated vs. virgin. Significantly enriched GO terms (functional annotation chart, P-value < 0.05) and details on associated genes for term of interests.