BACKGROUND
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 , there is growing interest in elucidating its molecular underpinnings in order to characterize its plasticity and evolution across animal taxa .
brain transcriptomes can be used to associate complex behaviours like mating with patterns of gene expression at a genomic scale  <cit> , and to identify conserved molecular pathways across taxa  <cit> .
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 .
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  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  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 .
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  <cit>  may be uncoupled in social insects .
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  <cit> , wasps  <cit> , termites  <cit> , and honey bees  <cit> .
honey bees  provide an excellent model system for molecular studies, as the behavioural ecology of mating is well understood.
queens mate on the wing when 6– <dig> days of age  <cit> .
virgin queens leave their nest during the afternoon  <cit> , and may fly several kilometres  <cit>  to mate at leks known as drone congregation areas  <cit> .
at the lek they mate with multiple males in quick succession  <cit> , before returning to their nest.
about 50 % of queens commence oviposition after one mating flight, while the remainder fly on 1– <dig> subsequent afternoon and again mate with multiple males .
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  <cit> .
mating is associated with profound behavioural  <cit>  and physiological  changes that are unrelated to the age of the queen.
in particular, virgin queens are photophilic and eager to fly  <cit> , and are aggressive towards other virgin queens  <cit> .
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  <cit> .
interestingly, similar physiological and behavioural changes to those found in mated queens are observed in queens subjected to double narcosis with carbon dioxide , providing an elegant tool for dissecting the molecular aspects of post-mating changes .
insects possess specialized receptor cells that can detect and measure environmental co <dig>  <cit> .
honey bees use these receptors to help regulate co <dig> concentrations in their nest by wing fanning .
narcosis, caused by artificial co <dig> exposure at much higher concentrations , usually triggers oviposition  by accelerating germ cell differentiation and stimulating the initial differentiation of the vitellarium  <cit> .
it has been hypothesized that co <dig> might achieve these results by exploiting the carbonic anhydrase pathways, the hypoxia-induced transcription factors or via activation of transferrin  <cit> .
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 dopaminereceptors .
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 .
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  <cit>  that have used the informative but less-powerful microarray approach  <cit> .
methods
biological material
in november  <dig> we reared  <dig> sister queens of standard australian commercial stock  using standard beekeeping techniques from a single source colony  <cit> .
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  <cit> .
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  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  were observed.
2) co2-treated queens  were subjected to 10 min of co2-narcosis on day  <dig> and on day  <dig>  the seven caged queens were placed in a zip-lock plastic bag, which was then flushed with compressed co <dig> 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  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 , directly into liquid nitrogen and then stored at −80 °c.
we used  <dig> queens per treatment group for rnaseq .
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  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  <cit> .
brains were placed individually in trizol and macerated to breakdown tissue.
after addition of chloroform, samples were briefly vortexed and then centrifuged at  <dig>  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  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  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™  <dig> spectrophotometer .
rna-sequencing
we sequenced rna with an illumina hiseq system using  <dig> lanes of a plate  and producing 50 bp single-end reads.
reads were checked with fastqc  for quality control and were subsequently processed with trimmomatic  <cit>  to remove adapters and low quality bases: reads less than  <dig> bases long were filtered out .
we also dropped reads that matched ribosomal rna sequences  by means of sortmerna  <cit> .
surviving reads were aligned with tophat for illumina  <cit>  to the latest version of the honey bee genome  available on beebase  by using the galaxy web-based platform .
mapped reads were converted into raw read counts with samtools idxstats  <cit>  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  <cit> .
read counts were log <dig> transformed to correct for the skew to zero and large values.
only genes with at least  <dig> reads per sample were kept in the analysis .
normalization was performed with trimmed mean of m-values , a method that is implemented in the edger bioconductor package  <cit> .
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  of 5 %.
for global analyses of gene expression we used hierarchical clustering  and principal component analysis in jmp pro  <dig>  .
to perform gene ontology  analyses we obtained drosophila melanogaster orthologs with blast  for honey bee genes that were significantly differentially expressed between treatments and computed functional annotation clustering in david version  <dig>  <cit>  with medium stringency and a cutoff of p-value <  <dig> .
to identify overrepresented biological functions  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.
comparative studies
we used venny  to overlap d. melanogaster ortholog matching lists of significantly differentially expressed genes and a hypergeometric test  to assess whether genes overlapping between studies occurred significantly more often than expected by chance.
we compared our study to three microarray studies  <cit>  and one qpcr study  <cit>  on honey bee gene expression after either mating or co2-narcosis .
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.
RESULTS
rna-sequencing
the sequencing produced  <dig> million reads per sample on average .
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  <dig> genes were significantly differentially expressed  in at least one of three pairwise comparisons across treatments ; this represents  <dig>  % of the  <dig>  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 .
this suggests that the mating treatment was the major driver of gene expression in the honey bee brain.fig.
<dig> patterns of gene expression in individual queens.
the heatmap was obtained after hierarchical clustering analysis of log2-transformed and normalized read counts that corresponded to  <dig> genes significantly differentially expressed across treatments.
the analysis shows that mated queens cluster separately as compared to virgin and co2-treated queens.
genes in red are significantly up-regulated while genes in blue are significantly down-regulated.
v = virgin, m = mated, c = co2-treated
we also performed a principal component analysis on the same set of  <dig> significantly differentially expressed genes .
again, most of the changes in gene expression  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.fig.
<dig> multivariate analysis of gene expression after mating and double narcosis with co <dig>  the principal component analysis of  <dig> significantly differentially expressed genes shows that the difference was primarily due to the mating process  and to a minor extent the co <dig> treatment
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:  <dig> genes at fdr <  <dig>  .
of these,  <dig> genes were up-regulated and  <dig> were down-regulated in mated queens.
our go analysis was based on the  <dig> d. melanogaster orthologs for these genes.
six metabolic pathways and  <dig> go terms were significantly overrepresented .
of particular interest were go terms associated with the following key biological processes: sensory perception , detection of stimulus , multiple metabolism-related go terms such as the fatty acid metabolic processes  and the immune-related go terms defense response , melanisation defense response  and innate immune response .
see additional file  <dig> for a list of the genes associated with these go terms.fig.
<dig> pairwise comparisons of gene expression between treatments.
all = genes that were significantly differentially expressed between treatment groups.
up = total numbers of unique genes that were up-regulated across all comparisons for each treatment group.
down = total numbers of unique genes that were down-regulated across all comparisons for each treatment group
<dig> e-05
<dig> e-04
<dig> e-03
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig> e-05
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
mated vs. co2
<dig>
<dig>
<dig>
<dig>
<dig>
<dig> e-03
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig> e-03
<dig>
<dig>
<dig>
<dig> e-04
mated vs. co2-treated queens
there were  <dig> significantly differentially expressed  genes between mated and co2-treated queens :  <dig> genes were up-regulated and  <dig> were down-regulated in mated queens.
go analysis of the matching  <dig> d. melanogaster orthologs revealed that  <dig> go terms were significantly overrepresented , including sensory perception  and response to organic substance .
see additional file  <dig> for a list of the genes associated with these go terms.
co2-treated vs. virgin queens
a much smaller set of  <dig> genes were significantly differentially expressed  between co2-treated and virgin queens :  <dig> genes were up-regulated and  <dig> were down-regulated in co2-treated queens.
the  <dig> d. melanogaster orthologs for these genes produced  <dig> significantly overrepresented go terms , including the biological process cognition .
see additional file  <dig> for a list of the genes associated with this go term.
comparative studies
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.
<cit> , 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.
<cit> , 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 ; 3) niño et al.
<cit> , where the brain transcriptomic profile was evaluated in virgin, co2-treated, and physically manipulated honey bee queens , 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 .
in contrast, at the whole dataset level, a significant number of differentially expressed genes were shared between our study and kocher et al.
<cit> .
however, no significant overlap was found in comparisons with data from kocher et al.
<cit>  and niño et al.
<cit>  .
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 .fig.
<dig> best candidate genes for mating process and reproductive activation.
patterns of expression for  <dig> genes that were significantly differentially expressed in our study and in kocher et al.
<cit>  and kocher et al.
<cit>  .
the table shows the total read counts for each gene in virgin, mated and co2-treated queens from our study plus the direction of expression in mated vs. virgin queens in the other two studies: up = up-regulated in mated queens; down = down-regulated; bold red = same pattern of expression between studies; n/a = difference in gene expression occurs in a comparison which is not relevant for our study
fig.
<dig> response to carbon dioxide.
patterns of expression of  <dig> genes that were included in brito et al.
<cit>  and were also significantly differentially expressed in our study.
lines above bars indicate the pairwise comparisons where the genes were significantly differentially expressed .
arrows show direction of expression for brito et al.
<cit>  : green = same trend as in our study; red = opposite trend.
see additional file 7 for more information on this comparative study
discussion
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  <dig> , neither inactivation nor afterpotential a , arrestin  <dig> , g protein beta-subunit 76c , chaoptin  and calhotin  were all down-regulated in mated queens compared to virgins .
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  <cit> , 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  <cit> , 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  and odorant binding proteins .
odour and pheromone perception are central to the social life of honey bees  and play an important role in mate location and efficient mating.
of particular interest is the presence of the pheromone-binding protein-related protein  <dig> , 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  <cit>  only four were also significantly differentially expressed in our study .
as previously reported  <cit> , genes broadly related to metabolism were also differentially expressed between mated and virgin queens .
mating flights are energetically expensive  <cit>  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  <cit> .
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 , peptidoglycan recognition protein sa  and gram-negative bacteria binding protein  <dig>  are important players in the toll pathway and their expression is usually triggered by a challenge from gram-positive bacteria  <cit> .
three additional genes, melanization protease  <dig> , serpin-27a  and hemolectin  are involved in the melanisation response  <cit>  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  <cit> .
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  <cit>  and other organisms  <cit> .
in drosophila, for example,  <dig> immune-related genes respond to mating, including defensin and pgrp-sa  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
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  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  and tor pathways are responsible for regulating growth and nutrition and are fundamental to the process of honey bee caste determination  <cit> .
we found that only one gene in these pathways was significantly differentially expressed between mated and virgin queens, phosphoinositide-dependent kinase  <dig> .
this gene was up-regulated in mated queens.
the fact that pdk <dig> was up-regulated in the only group of queens with activated ovaries in our studies is in line with the hypothesis that pdk <dig> 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 pdk <dig> compared to foragers that are nectar-biased  <cit> .
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  <cit> .
for example, the gene n-acetyldopamine is positively correlated with ovarian development  <cit> .
in queens, instead, dopamine and reproductive status are negatively correlated and dopamine levels decrease after mating  <cit> .
our study provides further support for this reversed relationship in queens, as dopamine n-acetyltransferase , 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 .
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  and no receptor potential a .
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  <cit> , 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.
<cit>  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  and obp19d were up-regulated in co2-treated queens, while three other obps  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 , 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  <cit> .
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  <cit> , 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.
comparative studies
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  <cit> .
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.
<cit>  are significantly congruent with our study, those from kocher et al.
<cit>  and niño et al.
<cit>  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.
<cit> , queens were collected two days after treatment; in contrast, kocher et al.
<cit>  and niño et al.
<cit>  analysed samples collected  <dig> 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 , 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  characterized by a particular social environment.
CONCLUSIONS
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  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 .
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  <cit> .
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 gse <dig> http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=oludywiunlwblgd&acc=gse <dig>  other data sets supporting the results of this article are included within the additional files.
ethics
for this study honey bees from commercial colonies were used and therefore no ethics statement is required.
additional files
additional file 1:
summary statistics for illumina hiseq, read preparation and alignment to the honey bee genome.
additional file 2:
comparative studies.
overlap analyses with previous microarray studies on the effect of mating and co2-narcosis on gene expression in honey bees.
additional file 3:
overlap analysis with kocher et al.
[
22].
genes that were significantly differentially expressed between mated and virgin queens in both studies.
additional file 4:
overlap analysis with kocher et al.
[
31].
genes that were significantly differentially expressed between mated and virgin queens in both studies.
additional file 5:
overlap analysis with niño et al.
[
50].
genes that were significantly differentially expressed between co2-treated and virgin queens in both studies.
additional file 6:
genes involved in eye development and visual perception that were differentially expressed between treatments in this study.
additional file 7:
results of comparative studies.
additional file 8:
pairwise comparisons.
list of the genes that were differentially expressed between mated and virgin queens at fdr <  <dig> .
additional file 9:
gene ontology analysis for the pairwise comparison mated vs. virgin.
significantly enriched go terms and kegg pathways  and details on associated genes for terms of interest.
additional file 10:
pairwise comparisons.
list of the genes that were differentially expressed between co2-treated and mated queens at fdr <  <dig> .
additional file 11:
gene ontology analysis for the pairwise comparison co
2
-treated vs. mated.
significantly enriched go terms  and details on associated genes for terms of interest.
additional file 12:
pairwise comparisons.
list of the genes that were differentially expressed between co2-treated and virgin queens at fdr <  <dig> .
additional file 13:
gene ontology analysis for the pairwise comparison co
2
-treated vs. virgin.
significantly enriched go terms  and details on associated genes for term of interests.
abbreviations
rnaribonucleic acid
co2carbon dioxide
rnaseqhigh-throughput rna sequencing
gdnagenomic deoxyribonucleic acid
ddh2odouble-distilled water
rrnaribosomal rna
fdrfalse discovery rate
qpcrquantitative real-time polymerase chain reaction
gogene ontology
orodorant receptor
obpodorant binding protein
iisinsulin/insulin-like signaling pathway
tortarget of rapamycin pathway
keggkyoto encyclopedia of genes and genomes
fabio manfredini and mark j f brown contributed equally to this work.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
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.
this study was funded by a bbsrc isis grant bb/j019453/ <dig>  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 dp <dig> and dp <dig>  fm was supported by a marie curie international incoming fellowship fp7-people-2013-iif- <dig> to mjfb.
we would like to thank dave galbraight  and alberto paccanaro  for support with analysis of rnaseq data and four anonymous reviewers for providing thoughtful insights that helped to improve the manuscript.

 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, there is growing interest in elucidating its molecular underpinnings in order to characterize its plasticity and evolution across animal taxa. brain transcriptomes can be used to associate complex behaviours like mating with patterns of gene expression at a genomic scale, and to identify conserved molecular pathways across taxa. 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. 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 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 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. 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 may be uncoupled in social insects. 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, wasps, termites, and honey bees. honey bees provide an excellent model system for molecular studies, as the behavioural ecology of mating is well understood. queens mate on the wing when 6– 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– subsequent afternoon and again mate with multiple males. 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 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. interestingly, similar physiological and behavioural changes to those found in mated queens are observed in queens subjected to double narcosis with carbon dioxide, providing an elegant tool for dissecting the molecular aspects of post-mating changes. insects possess specialized receptor cells that can detect and measure environmental co. honey bees use these receptors to help regulate co concentrations in their nest by wing fanning. narcosis, caused by artificial co exposure at much higher concentrations, usually triggers oviposition by accelerating germ cell differentiation and stimulating the initial differentiation of the vitellarium. it has been hypothesized that co 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 dopaminereceptors. 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. 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 that have used the informative but less-powerful microarray approach. methods biological material in november we reared sister queens of standard australian commercial stock 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 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 were observed. 2) co2-treated queens were subjected to 10 min of co2-narcosis on day and on day the seven caged queens were placed in a zip-lock plastic bag, which was then flushed with compressed co 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 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, directly into liquid nitrogen and then stored at −80 °c. we used queens per treatment group for rnaseq. 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 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  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 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 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™ spectrophotometer. rna-sequencing we sequenced rna with an illumina hiseq system using lanes of a plate and producing 50 bp single-end reads. reads were checked with fastqc for quality control and were subsequently processed with trimmomatic to remove adapters and low quality bases: reads less than bases long were filtered out. we also dropped reads that matched ribosomal rna sequences by means of sortmerna. surviving reads were aligned with tophat for illumina to the latest version of the honey bee genome available on beebase by using the galaxy web-based platform. 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 log transformed to correct for the skew to zero and large values. only genes with at least reads per sample were kept in the analysis. normalization was performed with trimmed mean of m-values, 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 of 5 %. for global analyses of gene expression we used hierarchical clustering and principal component analysis in jmp pro. to perform gene ontology analyses we obtained drosophila melanogaster orthologs with blast for honey bee genes that were significantly differentially expressed between treatments and computed functional annotation clustering in david version with medium stringency and a cutoff of p-value < . to identify overrepresented biological functions 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. comparative studies we used venny to overlap d. melanogaster ortholog matching lists of significantly differentially expressed genes and a hypergeometric test to assess whether genes overlapping between studies occurred significantly more often than expected by chance. we compared our study to three microarray studies and one qpcr study on honey bee gene expression after either mating or co2-narcosis. 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.  rna-sequencing the sequencing produced million reads per sample on average. 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 genes were significantly differentially expressed in at least one of three pairwise comparisons across treatments; this represents  % of the 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. this suggests that the mating treatment was the major driver of gene expression in the honey bee brain.fig. patterns of gene expression in individual queens. the heatmap was obtained after hierarchical clustering analysis of log2-transformed and normalized read counts that corresponded to genes significantly differentially expressed across treatments. the analysis shows that mated queens cluster separately as compared to virgin and co2-treated queens. genes in red are significantly up-regulated while genes in blue are significantly down-regulated. v = virgin, m = mated, c = co2-treated we also performed a principal component analysis on the same set of significantly differentially expressed genes. again, most of the changes in gene expression 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.fig. multivariate analysis of gene expression after mating and double narcosis with co the principal component analysis of significantly differentially expressed genes shows that the difference was primarily due to the mating process and to a minor extent the co treatment 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: genes at fdr < . of these, genes were up-regulated and were down-regulated in mated queens. our go analysis was based on the d. melanogaster orthologs for these genes. six metabolic pathways and go terms were significantly overrepresented. of particular interest were go terms associated with the following key biological processes: sensory perception, detection of stimulus, multiple metabolism-related go terms such as the fatty acid metabolic processes and the immune-related go terms defense response, melanisation defense response and innate immune response. see additional file for a list of the genes associated with these go terms.fig. pairwise comparisons of gene expression between treatments. all = genes that were significantly differentially expressed between treatment groups. up = total numbers of unique genes that were up-regulated across all comparisons for each treatment group. down = total numbers of unique genes that were down-regulated across all comparisons for each treatment group e-05 e-04 e-03 e-05 mated vs. co2 e-03 e-03 e-04 mated vs. co2-treated queens there were significantly differentially expressed genes between mated and co2-treated queens: genes were up-regulated and were down-regulated in mated queens. go analysis of the matching d. melanogaster orthologs revealed that go terms were significantly overrepresented, including sensory perception and response to organic substance. see additional file for a list of the genes associated with these go terms. co2-treated vs. virgin queens a much smaller set of genes were significantly differentially expressed between co2-treated and virgin queens: genes were up-regulated and were down-regulated in co2-treated queens. the d. melanogaster orthologs for these genes produced significantly overrepresented go terms, including the biological process cognition. see additional file for a list of the genes associated with this go term. comparative studies 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; 3) niño et al., where the brain transcriptomic profile was evaluated in virgin, co2-treated, and physically manipulated honey bee queens, 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. 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.. 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.fig. best candidate genes for mating process and reproductive activation. patterns of expression for genes that were significantly differentially expressed in our study and in kocher et al. and kocher et al.. the table shows the total read counts for each gene in virgin, mated and co2-treated queens from our study plus the direction of expression in mated vs. virgin queens in the other two studies: up = up-regulated in mated queens; down = down-regulated; bold red = same pattern of expression between studies; n/a = difference in gene expression occurs in a comparison which is not relevant for our study fig. response to carbon dioxide. patterns of expression of genes that were included in brito et al. and were also significantly differentially expressed in our study. lines above bars indicate the pairwise comparisons where the genes were significantly differentially expressed. arrows show direction of expression for brito et al.: green = same trend as in our study; red = opposite trend. see additional file 7 for more information on this comparative study discussion 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, neither inactivation nor afterpotential a, arrestin, g protein beta-subunit 76c, chaoptin and calhotin were all down-regulated in mated queens compared to virgins. 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 and odorant binding proteins. odour and pheromone perception are central to the social life of honey bees and play an important role in mate location and efficient mating. of particular interest is the presence of the pheromone-binding protein-related protein, 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. as previously reported, genes broadly related to metabolism were also differentially expressed between mated and virgin queens. 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, peptidoglycan recognition protein sa and gram-negative bacteria binding protein 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, serpin-27a and hemolectin are involved in the melanisation response 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 and other organisms. in drosophila, for example, 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. 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. 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 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 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. this gene was up-regulated in mated queens. the fact that pdk was up-regulated in the only group of queens with activated ovaries in our studies is in line with the hypothesis that pdk 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 pdk 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. 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, 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. 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 and no receptor potential a. 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 and obp19d were up-regulated in co2-treated queens, while three other obps 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, 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, 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. comparative studies 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. 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 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, 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 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 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. 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 gse http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=oludywiunlwblgd&acc=gse other data sets supporting the results of this article are included within the additional files. ethics for this study honey bees from commercial colonies were used and therefore no ethics statement is required. additional files additional file 1: summary statistics for illumina hiseq, read preparation and alignment to the honey bee genome. additional file 2: comparative studies. overlap analyses with previous microarray studies on the effect of mating and co2-narcosis on gene expression in honey bees. additional file 3: overlap analysis with kocher et al. [ 22]. genes that were significantly differentially expressed between mated and virgin queens in both studies. additional file 4: overlap analysis with kocher et al. [ 31]. genes that were significantly differentially expressed between mated and virgin queens in both studies. additional file 5: overlap analysis with niño et al. [ 50]. genes that were significantly differentially expressed between co2-treated and virgin queens in both studies. additional file 6: genes involved in eye development and visual perception that were differentially expressed between treatments in this study. additional file 7: results of comparative studies. additional file 8: pairwise comparisons. list of the genes that were differentially expressed between mated and virgin queens at fdr < . additional file 9: gene ontology analysis for the pairwise comparison mated vs. virgin. significantly enriched go terms and kegg pathways and details on associated genes for terms of interest. additional file 10: pairwise comparisons. list of the genes that were differentially expressed between co2-treated and mated queens at fdr < . additional file 11: gene ontology analysis for the pairwise comparison co 2 -treated vs. mated. significantly enriched go terms and details on associated genes for terms of interest. additional file 12: pairwise comparisons. list of the genes that were differentially expressed between co2-treated and virgin queens at fdr < . additional file 13: gene ontology analysis for the pairwise comparison co 2 -treated vs. virgin. significantly enriched go terms and details on associated genes for term of interests. abbreviations rnaribonucleic acid co2carbon dioxide rnaseqhigh-throughput rna sequencing gdnagenomic deoxyribonucleic acid ddh2odouble-distilled water rrnaribosomal rna fdrfalse discovery rate qpcrquantitative real-time polymerase chain reaction gogene ontology orodorant receptor obpodorant binding protein iisinsulin/insulin-like signaling pathway tortarget of rapamycin pathway keggkyoto encyclopedia of genes and genomes fabio manfredini and mark j f brown contributed equally to this work. competing interests the authors declare that they have no competing interests. authors’ contributions 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. this study was funded by a bbsrc isis grant bb/j019453/ 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 dp and dp fm was supported by a marie curie international incoming fellowship fp7-people-2013-iif- to mjfb. we would like to thank dave galbraight and alberto paccanaro for support with analysis of rnaseq data and four anonymous reviewers for providing thoughtful insights that helped to improve the manuscript.
in addition, we compared our results to previous studies that used microarray techniques across a range of experimental time-points. 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. our study is an important step in understanding the molecular mechanisms regulating post-mating behavioural transitions in a natural system. 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. apis melliferamatingcarbon dioxidedouble narcosistranscriptomicsrnaseqaustralian honey beesgene expressionbehaviourbrainbbsrc isisid0eepbg1329brown mark j f royal holloway research strategy foundid0esqbg1330brown mark j f leverhulme grantid0e5rbg1331brown mark j f australian research council discovery grantid0emtbg1332oldroyd benjamin p australian research council discovery grantid0e1ubg1333oldroyd benjamin p marie curie international incoming fellowshipid0eiwbg1334manfredini fabio issue-copyright-statement© the author 2015 mating is a complex process, which is frequently associated with behavioural and physiological changes. 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. 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. a subset of these changes in gene expression were also detected in co2-treated queens, as predicted from previous physiological studies.
 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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords apis melliferamatingcarbon dioxidedouble narcosistranscriptomicsrnaseqaustralian honey beesgene expressionbehaviourbrainbbsrc isisid0eepbg1329brown mark j f royal holloway research strategy foundid0esqbg1330brown mark j f leverhulme grantid0e5rbg1331brown mark j f australian research council discovery grantid0emtbg1332oldroyd benjamin p australian research council discovery grantid0e1ubg1333oldroyd benjamin p marie curie international incoming fellowshipid0eiwbg1334manfredini fabio issue-copyright-statement© the author 2015 