Regulation of behaviorally associated gene networks in worker honey bee ovaries

Ovaries play important roles in behavioral regulation in animals as diverse as mammals (Kemnitz et al., 1989) and insects (Cupp and Collins, 1979; Hancock and Foster, 2000; Klowden, 1990). Recently, it was suggested that ovaries also affect complex behavioral patterns in functionally sterile honey bee (Apis mellifera) workers (Amdam et al., 2006; Wang et al., 2009; Wang et al., 2010). The honey bee is the most-studied social insect and is a model organism for understanding the regulation of behavior (Amdam et al., 2006; Evans and Wheeler, 1999; Giray and Robinson, 1996; Hunt et al., 1995; Lindauer, 1953; Pennisi, 2006). Worker honey bees progress through an age-associated behavioral sequence. They primarily conduct within-nest tasks during their first 2–3 weeks of adult life and then transition to performing foraging tasks outside the colony. The timing (age) of foraging onset varies among individuals. Foraging workers also show different behaviors; some bees prefer nectar collection and others prefer pollen collection (Page and Fondrk, 1995).

This behavioral plasticity is influenced by ovary size (Amdam et al., 2006; Wang et al., 2009; Wang et al., 2010). Workers that receive an implant of additional ovaries, i.e. surgery to enhance total ovary mass, initiate foraging earlier than sham controls (Wang et al., 2010). Bees with naturally large ovaries also transition from nest tasks to foraging at younger ages, and they collect a greater amount of pollen than bees with small ovaries. These associations are found within wild-type (unselected) European honey bees (Amdam et al., 2006), in crosses between unselected European and Africanized honey bees (Siegel, 2011), and between strains of two European honey bee genotypes that were bidirectionally bred for differential foraging behavior (Hunt et al., 2007; Page et al., 1995). These artificially selected genotypes are referred to as High vs Low pollen hoarding strains and differ in several traits: the High strains have larger ovaries, forage earlier in life and collect more pollen than the Low strain bees. Thus, ovary size and foraging behavior are heritable traits in honey bees (Amdam et al., 2009), and variation in both of these phenotypes maps to overlapping genomic regions [quantitative trait loci (QTL)] (Graham et al., 2011; Linksvayer et al., 2009b; Wang et al., 2009). Linkage mapping of these traits has been performed repeatedly and independently in crosses derived from different sources, i.e. from pollen-hoarding strain bees as well as Africanized and European wild-type bees (Graham et al., 2011; Linksvayer et al., 2009b; Wang et al., 2009).

The mapping studies have provided positional candidate genes for phenotypic variation in worker ovary size and behavior. Expression levels of some of these candidate genes have been examined in worker brains and fat bodies (abdominal adipose tissue). Fat body expression of hormone receptor-like in 46 (HR46) and phosphoinositide-dependent kinase-1 (PDK1) correlates with ovary size as well as worker division of labor (Wang et al., 2009). HR46 is a possible nuclear hormone receptor for ecdysteroid hormones and PDK1 is an element in the nutrient sensing target of rapamycin-, epidermal growth factor (Egf)and insulin-signaling pathways (Kamakura, 2011). However, although previous studies have documented phenotypic, functional genomic and genetic evidence for the links between honey bee division of labor and reproductive physiology (Amdam et al., 2006; Wang et al., 2009; Wang et al., 2010), less is known about how ovarian tissue affects worker behavior per se.

Here, we studied shared features of worker ovary and brain transcriptomes in order to begin to tease apart how ovaries may function in honey bee division of labor. First, we analyzed gene expression globally, as well as locally in co-expressed gene clusters between ovaries obtained from High and Low pollen-hoarding strain bees. This analysis was contrasted with six published brain transcriptional studies. Two of these studies described associations between gene expression and division of labor in the honey bee (Whitfield et al., 2003) and in paper wasps (Toth et al., 2010). Four others examined the pheromonal and endocrine effects on worker behaviors: the delay of foraging onset by queen mandibular pheromone (QMP) (Grozinger et al., 2003; Kocher et al., 2010); the acceleration of foraging onset by methoprene, a juvenile hormone (JH) mimic (Whitfield et al., 2006); and honey bee aggression by alarm pheromone (Alaux et al., 2009b). Our comparative analysis generated a list of candidate biological processes for shared behavioral and ovarian gene regulation in honey bees. Finally, as variation in ovary size is associated with worker behavioral differences (Amdam et al., 2006), we independently tested whether a subset of genes was also differentially expressed in wild-type worker ovaries of different sizes. Our results show that gene expression patterns in worker ovaries mirror brain transcriptional responses triggered by socio-environmental factors, and identify several physiological mechanisms that may facilitate the link between ovarian activity and division of labor in honey bee workers, including differences in metabolic pathways and endocrine physiology.

Honey bees, Apis mellifera L., from the High and Low strains were maintained at Arizona State University, Phoenix, AZ, USA, and at the University of California, Davis, CA, USA. Two High strain colonies and two Low strain colonies were used as donors for experimental workers for the microarray study. Ovaries were collected from newly emerged worker bees (<24 h old) and used for transcriptional profiling. Newly emerged workers have encountered few adult stimuli that could affect their ovary physiology, thus potential differences in the ovary between strains should be primarily due to genotype or ovariole number. Both ovaries from each worker were carefully dissected in sterile Ringers solution (on ice). For both strains, ovaries of 10 individual bees (five bees from each of the two colonies) were pooled to form a biological sample. Four biological samples were prepared from each genotype.

Ovaries were flash-frozen in liquid nitrogen (N₂) and stored at –80°C. RNA was extracted using TRIzol Reagent (Invitrogen, cat. no. 15596-018, Carlsbad, CA, USA) (Kobayashi et al., 2006; Michaut et al., 2003). DNase (Qiagen, Valencia, CA, USA) treatment was used to remove residual genomic DNA. The integrity of extracted RNA was tested by capillary electrophoresis on an RNA6000 BioAnalyzer (Agilent Technologies, Palo Alto, CA, USA), and only high-quality samples were used for subsequent analyses.

RNA (200 ng) was amplified with the MessageAmp II aRNA Amplification Kit (Ambion, Austin, TX, USA). Amplified RNA (2μg) was directly labeled with Cy3 or Cy5 dye using a Kreatech labeling kit (BioMicro, Salt Lake City, UT, USA). Labeled probe (120 pmol) from two samples was then hybridized to the whole-genome oligonucleotide arrays supplied by the W. M. Keck Center for Comparative and Functional Genomics at the University of Illinois, Urbana-Champaign. The High and Low ovaries were hybridized with dye-swaps (N=4 for each group). Arrays were scanned using the Axon Genepix 4000B scanner (Molecular Devices, Sunnyvale, CA, USA) and GENEPIX software (Agilent Technologies, Santa Clara, CA, USA), and raw data were imported into SAS (SAS Institute, Inc., Cary, NC, USA) for analysis.

Features with an intensity below 300 (the average background intensity on the arrays) were removed from the analysis. Genes with less than six observations (out of a possible 16) were also removed. Data were log-transformed and normalized using a mixed-model ANOVA (SAS) with the following model: y = dye array + block + dye×array + ε, where y is expression, ε is error, dye and block are fixed effects, and array and its interactions are random effects. Detection of significance for differential expression of residuals was performed using a mixed-model ANOVA with the model: y μ + genotype + spot + dye + array + ε, where y is the residual from the previous model, genotype and dye are fixed effects, and array is a random effect. P-values were corrected for multiple testing using a false discovery rate (FDR) adjustment (proc MULTTEST, SAS). We chose an FDR threshold of ≤0.05, suggesting that less than 5% of the significant transcripts should be false positives.

Gene ontology (GO) enrichment analysis was performed only on differentially expressed genes in either High and Low strain ovaries. Only transcripts (N=1534) identified as having orthologs in Drosophila were included in the analysis to exclude multiple transcripts from the same gene and transcripts for which the genomic location has not been assigned. Genes were included as being significantly differentially expressed if the FDR-corrected P-values were ≤0.05. Enrichment was determined using GO-getter (http://chem.colorado.edu/knightgroup/) (Argast et al., 2009).

The differentially expressed genes were then mapped onto honey bee chromosomes for detecting genomic clusters. We used transcripts (1794) that have been assigned for genomic locations and tested whether these genes were clustered more than expected and whether such clusters were located within the previously mapped QTL for foraging behavior (Hunt et al., 2007; Hunt and Page, 1995), sucrose sensitivity (Rueppell et al., 2006) and ovary size (Linksvayer et al., 2009b; Wang et al., 2009). Clustering was assessed across each chromosome by calculating a chi-square value for the observed number of significant genes within a sliding window of 15 loci and a step size of three loci, relative to the expected number of significant genes within that window, given the number of significant differentially expressed genes for each chromosome. We established conservative 0.05 chromosome-specific significance thresholds for the chi-square statistic using 10,000 random permutations of the data set with respect to significance of differential expression.

We contrasted our list of 2151 differentially expressed ovarian genes to four previously published brain transcriptomes from wild-type worker bees, which were exposed to different behavioral regulators (Grozinger et al., 2003; Kocher et al., 2010; Whitfield et al., 2006; Whitfield et al., 2003). We also contrasted our ovarian transcriptomes to brain transcriptomes from paper wasps, which were related to foraging behavior and reproduction (Toth et al., 2010). Finally, we used the results of a study that identified brain transcriptional changes associated with exposure to alarm pheromone as a negative control (Alaux et al., 2009b). As there is no evidence that alarm pheromone affects reproductive physiology in honey bee workers, we did not expect overlap between our results and those of Alaux et al. (Alaux et al., 2009b). Significantly overlapping genes, i.e. those that were differentially expressed in both the ovarian and brain transcriptomes, were identified using a right-tailed Fisher's exact test. Furthermore, GO analyses were performed on the overlapping genes within each pollen-hoarding strain, using GOToolBox (http://genome.crg.es/GOToolBox/) (Alaux et al., 2009a).

Five genes [HR46, ftz transcription factor 1 (ftz-f1), tyramine receptor (TYR), major royal jelly protein 1 (MRJP1) and juvenile hormone-inducible protein 26 (JHi-26)] were chosen for verification of gene expression in samples that were used for the array. HR46 and ftz-f1 are nuclear receptors that may bind ecdysteroid hormones produced by the ovary of adult female insects, whereas TYR is a G protein-coupled receptor (GPCR) predicted to bind tyramine neurotransmitter. Both HR46 and TYR are also positional candidate genes in QTL regions that influence foraging behavior. MRJP1, a protein component in royal jelly, affects queen–worker differentiation (Kamakura, 2011) and potentially influences ovary physiology in workers. JHi-26 is a JH-inducible gene and can be triggered by methoprene, a JH analog. Ftz-f1, TYR, MRJP1 and JHi-26 were selected for profiling because they were differentially expressed between High and Low ovaries in the microarray assay (see Results). HR46 signals on the array were too low to pass the threshold for the analysis, but HR46 candidacy was inferred from prior testing of HR46 expression in ovaries from newly emerged bees of the High and Low strains [for which each biological sample consisted of 20 ovaries (pooled) from 10 individual bees from two colonies per strain, N=12, Student's t-test: t_1,22=–2.4310, P=0.0237; supplementary material Fig. S1].

Both ovary size and behavioral phenotype affect ovarian physiology and gene expression. Using qRT-PCR, we tested whether relative quantification (RQ) of the candidate genes for ovary and behavioral regulation (HR46, ftz-f1, TYR, MRJP1 and JHi-26) were associated with ovary size in wild-type, newly emerged bees. For these bees, we used five wild colony sources for testing the association between candidate genes and ovary size. After ovary dissections, ovariole number was counted using a dissecting microscope (Leica, MZ 125, Buffalo Grove, IL, USA). In order to maximize the difference in gene expression between big and small ovaries, we used five and 12 as cut-off thresholds for the ovariole number. Ovaries from bees in which both ovaries had either many ovarioles (13±0.11, i.e. summing over the two ovaries) or few ovarioles (4±0.06 in total) were flash-frozen in liquid N₂ and stored at –80°C. Biological samples were made by pooling 10 pairs of ovaries within the same size category, and these were assigned to either a large ovary (LAO) or a small ovary (SMO) group. We prepared 18 samples for each of the LAO and SMO groups.

Protocols for RNA extraction and DNase treatment were identical to those described above for the microarray study. cDNA was synthesized using TaqMan Reverse Transcription Reagents (Applied Biosystems, Foster City, CA, USA). Two-step qRT-PCR (real-time) was performed in triplicate using ABI Prism 7500 (Applied Biosystems), and the data were analyzed using the Delta-Delta CT method (Pfaffl, 2001) with actin (GenBank: XM_623378) as a reference gene. This gene is stably expressed across honey bee tissues, and is standardly used for gene expression studies in the bee (Lourenço et al., 2008). By monitoring negative control samples (without reverse transcriptase) and melting curves, we verified that the qRT-PCR assay was not confounded by DNA contamination or primer dimer (Vandesompele et al., 2002).

The primer sequences of the five target genes and actin were developed using Primer3 (http://frodo.wi.mit.edu/primer3/) (see supplementary material Table S1).

Expression data of candidate genes were log transformed to confer approximate normality (Wang et al., 2009). The resulting values conformed to assumptions of Student's t-test and ANOVA as assessed by normal probability plots of residuals, as well as by Bartlett's and Levene's tests for homogeneity of variance. Student's t-test was used to compare the candidate gene expressions in the ovary between High and Low strain bees. One-way ANOVA was used to compare the expression level of candidate genes between LAO and SMO groups. The log-transformed expression level of candidate genes was a dependent variable and ovary size (LAO or SMO) was an independent variable. Fisher's least significant difference (LSD) test was used for post hoc comparisons.

We compared ovarian transcript patterns from newly emerged High and Low pollen-hoarding strain bees using an established microarray approach (Grozinger et al., 2003; Whitfield et al., 2002). We found that 2151 transcripts, 20.3% of the 10,586 transcripts of the array, were differentially expressed between the High and Low strain sources at an FDR of 0.05 (supplementary material Table S2). Thus, the ovarian transcriptomes of the two pollen-hoarding strains were measurably different during the first hours of adult life.

GO analysis revealed the functional distribution pattern of differentially regulated genes in High and Low strain ovaries annotated against the Drosophila genome (supplementary material Fig. S2). A number of functions were ascribed, including the category ‘sensory perception’ with 10 olfactory receptors and two gustatory receptors expressed in ovaries. It is possible that the corresponding gene products incorporate into ovarian membranes and receive signals by binding molecules that circulate in the insects' hemolymph (blood). Such molecules can be internalized from the environment or secreted from tissues.

All differentially expressed genes were assigned to function groups within three functional domains of the GO analysis – ‘biological process’,‘molecular function’ and ‘cellular component’ (Fig. 1). This comparison between genotypes demonstrates that ovarian gene expression can differ dramatically between young worker bees. Transcripts putatively involved in defense responses (immune response) (P=0.033), regulation of apoptosis (P=0.041), detection of abiotic stimulus (P=0.010) and electron transport (P=0.041) were significantly enriched in Low strain ovaries, while kinase binding (P=0.016) was significantly higher in High strain ovaries (Fig. 1, see supplementary material Table S3 for details on statistics). These results indicate that Low strain worker ovaries may be more sensitive to stimuli and show higher levels of apoptosis compared with High strain worker ovaries. This interpretation matches the finding that Low strain workers are more sensitive to ovarian inhibition and less likely to complete egg development than High strain workers (Amdam et al., 2006). Thus, ovarian gene expression differences detected early in life (Fig. 1) may have physiological implications for honey bee workers.

Using chromosome-wide thresholds of significance (P≤0.05), we found clusters of differentially expressed genes in four genomic regions, one on each of chromosomes 5, 9, 11 and 14 (Table 1). None of these clusters overlapped with previously identified QTL for ovary size or behavior. Interestingly, on chromosome 11, the cluster contained six transcripts (GB14888, GB11786, GB11890, GB16246, GB11022 and GB14639) of the yellow protein family, which putatively encode the major royal jelly proteins 1, 4, 6, 2, 7 and 8 (MRJP1, MRJP4, MRJP6, MRJP2, MRJP7 and MRJP8, respectively). MRJPs are abundant in royal jelly, a food secreted from hypopharyngeal (head) glands of workers. Yet, transgenic fruit flies with MRJP1 monomer expressed internally (in fat) show activation of Egf signaling and phenotypes (Kamakura, 2011), indicating internal roles of MRJP transcripts – perhaps in a variety of tissues including ovaries.

The genes that were differently expressed between ovaries from the High and Low strains were compared with previously published wild-type brain transcriptomes regulated by QMP, nurse/foraging behavior and treatment with methoprene (a JH mimic) (Grozinger et al., 2003; Kocher et al., 2010; Whitfield et al., 2006; Whitfield et al., 2003). We found significant overlap between ovarian and brain data (right-tailed Fisher's exact test; Table 2). P-values less than 0.001 were significant after Bonferroni correction for multiple testing (P<0.004).

We selected genes for GO analysis that were differently regulated in ovaries of High or Low strain bees as well as regulated by QMP, nurse/foraging behavior or JH mimic in the brain of wild-type worker bees (Table 2, groups with light, medium or dark shading, respectively; Fig. 2). Similar differences in biological processes were found between High and Low strain ovaries in each of the three comparisons. The overlap between ovarian gene expression and brain transcriptomes was biased towards carbohydrate, lipid and amine metabolism in High strain bees, whereas the corresponding overlap for Low strain ovaries was neural network development and responses to different stimuli, including ecdysone hormone response (Fig. 2). Thus, we found that selection on pollen hoarding has changed ovarian expression of genes involved in metabolic control and signaling, sensory physiology and hormone responses, and the same genes respond in the brain during behavioral plasticity in wild-type bees.

There was no significant overlap between ovarian transcriptome in this study and the alarm-pheromone brain transcriptome selected as our negative control (P>0.004; Table 2). We also contrasted our ovarian gene list to the paper wasp brain transcriptomes related to foraging behavior and reproduction (Toth et al., 2010). We did not find significant overlap in gene expression pattern related to either foraging or reproduction (Table 2). These results are perhaps unsurprising because paper wasps (Polistes) and honey bees have fundamental difference in reproductive physiologies (Toth et al., 2010) that could differentially affect behavior. However, the common mechanisms between honey bees and wasps may be better revealed in the future by increasing the coverage of the wasp genome (Toth et al., 2010).

We selected five genes for expression verification that are involved in metabolic control and signaling, sensory physiology and hormone responses (i.e. HR46, Ftz-f1, JHi-26, MRJP1 and TYR; for a summary on function, see Materials and methods). The results of all five genes were consistent with the results from the array (Fig. 3) – all have significantly (P≤0.05) higher expression in Low strain ovaries than in High strain ovaries (one-tailed Student's t-test, TYR, P=0.039; ftz-f1, P=0.0458; HR46, P=0.0225; MRJP1, P=0.0004; JHi-26, P=0.0069).

Gene expression was quantified in wild-type worker ovaries of different size. These ovary size groups had significant effects on overall gene expression (one-way ANOVA, F_5,26=13.7682, P<0.0001, N=16; Fig. 4). TYR mRNA levels were significantly higher in the SMO group (Fisher's LSD, P=0.0016; Fig. 4A), whereas ftz-f1 and HR46 were more strongly expressed in the LAO group (Fisher's LSD, P=0.0102 and 0.0416, respectively; Fig. 4B,C). MRJP1 and JHi-26 expression, in contrast, did not vary by ovary size (Fisher's LSD, P=0.3844 and 0.9697, respectively; 3). These outcomes suggest that the ovarian expression levels of HR46, ftz-f1, MRJP1 and JHi-26 are best explained by selection (genotype) in this study, whereas the amount of TYR transcript can be robustly associated with ovary size in the wild-type as well as the pollen-hoarding strain bees.

Here, we describe transcriptional differences between the ovaries of newly emerged bees from honey bee strains bidirectionally selected for behavior. A number of genes were associated with the biological process ‘sensory perception’ or the functional category ‘reproduction’ (GO analysis, supplementary material Fig. S2). Differentially expressed genes included ecdysone-inducible nuclear hormone receptors and olfactory and gustatory receptors. Taken together, these results suggest that the honey bee worker ovary is a dynamic organ (Wang et al., 2010), although it is often presumed to be inactive under normal social conditions (Foster and Ratnieks, 2001). Furthermore, the significant overlap observed between the present study and previously published brain transcriptomes suggests that both ovary and brain physiologies are affected by systemic factors such as JH, ecdysteroids and QMP, and further supports the hypothesis that each of these factors are key participants in regulatory networks for social behavior in this species.

Honey bee worker ovaries produce ecdysteroids during the first day of adult life (Amdam et al., 2010), and the ecdysone receptor gene (EcR) is expressed in worker and queen ovaries (Takeuchi et al., 2007). Ecdysteroids signal oogenesis in Tribolium, Drosophila and mosquitoes (Culex pipiens) (Bernardi et al., 2009; Khater et al., 1994; Xu et al., 2010), but oogenesis is uncommon in worker bees because yolk uptake and oocyte maturation is suppressed pheromonally when a queen is present in the colony (Atkins et al., 1975).

We found that three genes associated with the ecdysteroid cascade [HR46, ftz-f1 and ecdysone-induced protein 75 (E75)] were more highly expressed in Low pollen-hoarding strain ovaries than in the High strain. At adult emergence, Low strain bees also have the higher ecdysteroid titers, and the plausible source of this hormone is the ovary (Amdam et al., 2010). Yet, elevated ecdysteroid signaling and/or sensing in Low strain ovaries does not correlate with oogenesis (Amdam et al., 2006). Rather, ecdysteroids can induce oocyte apoptosis and follicle atresia (Paul et al., 2005; Soller et al., 1999; Uchida et al., 2004), and can be supported by increased apoptosis regulation in Low strain ovaries (Fig. 1A). This hypothesis is consistent with Low strain workers being more readily inhibited from activating their ovaries for egg laying (Amdam et al., 2006), perhaps because of heightened sensitivity to social inhibition and increased rates of follicle apoptosis (Hunt et al., 2007). Future studies may answer whether ecdysteroids regulate apoptosis in honey bee ovaries.

Ecdysteroid titers are very low in the hemolymph of adult bees, suggesting that ecdysteroids have lost their function in adult workers (Hartfelder et al., 2002). Recently, many studies have indicated that ecdysteroid cascades in the fat body and brain may be involved in social behavior (Velarde et al., 2009; Wang et al., 2009; Wang et al., 2010). In this study, we found a significant overlap of ecdysteroid-related processes between the Low ovarian and brain transcriptomes associated with QMP exposure (Fig. 2A) and foraging behavior (Fig. 2C). This implies that ecdysteroids may be involved in both brain and ovary physiologies related to QMP and foraging behavior. However, further investigation is needed to test the role of ecdysteroids in honey bee brain regarding social behavior.

Our comparison between ovarian and brain transcriptomes relates biological processes of the ovary to brain mRNA expression profiles associated with exposure to QMP, methoprene and the nursing or foraging behavioral state. In every comparison, similar patterns of ovary–brain overlap were found (Fig. 2). Was this surprising? Worker behavior, including nursing and foraging, is modulated by QMP (Pankiw et al., 1998), JH (Huang et al., 1991) and ovarian physiology (Amdam et al., 2006; Nelson et al., 2007; Robinson, 1992). Therefore, the overlap we revealed (Fig. 2) may be a signature of (largely) systemically acting gene networks associated with the division of labor. Moreover, the clear separation between pollen-hoarding strains in their overlap of ovary and brain transcriptomes supports the hypothesis that the bidirectional selection on these strains acted on gene networks that connect behavior to physiology (Amdam et al., 2010; Page et al., 1998).

QMP emitted by honey bee queens inhibits ovary activation (Hoover et al., 2003), influences age-related division of labor (Pankiw et al., 1998) and affects brain development in worker bees (Morgan et al., 1998). The brain's transcriptional response to QMP is also negatively correlated with worker ovariole number (Kocher et al., 2010). This correlation points to unexplained connections between brain function and ovary size. We found that the functional category‘olfactory behavior’ was over-represented in the overlap between Low strain ovaries and brains regulated by QMP (Fig. 2A). Moreover, we identified 10 olfactory receptor (OR) genes differentially expressed between High and Low strain ovaries. ORs belong to the GPCR gene family, which is also enriched in ovaries of Low strain genotype bees (P=0.0673; supplementary material Table S3). GPCRs comprise a large protein family of trans-membrane receptors that bind light-sensitive compounds, odors, pheromones, hormones and neurotransmitters. ORs are generally expressed in cell membranes of olfactory receptor neurons (antennae) (Getz and Akers, 1994; Murphy et al., 2003; Vosshall et al., 1999). However, additional functions of OR proteins have been suggested (Feingold et al., 1999; Itakura et al., 2006; Xu et al., 2000) because of their localization to different ectopic tissues (non-olfactory tissues), such as germ cells and testis tissue in mammals (Fukuda et al., 2004; Marchand, 2003; Spehr et al., 2003; Vanderhaeghen et al., 1997; Ziegler et al., 2002). Recent experiments also show that one OR in honey bee antennal neurons can respond to a component of QMP (Wanner et al., 2007). It is possible that connections between brain function and ovarian traits in worker bees could be caused by correlated expression of ORs. This hypothesis can be tested in future studies.

JH is a systemic hormone with effects on insect metabolic biology, reproductive physiology and adult behavior. JH and ecdysteroid titers are often negatively correlated in insects, with opposing regulation demonstrated in Drosophila melanogaster (reviewed by Riddiford, 1996), mosquitoes (Aedes atropalpus) L (Birnbaum et al., 1984) and honey bees (Rembold, 1987). Using brain transcriptome data from exposure to methoprene (a JH mimic), we found that Low strain ovaries bias functional processes towards ecdysteroid and neural regulation in contrast to the High strain bees (Fig. 2C). It has previously been demonstrated that worker sensitivity towards JH responses is conditional on their Low vs High genotypes (Amdam et al., 2007; Amdam et al., 2010), i.e. Low strain worker bees show reduced JH sensitivity, perhaps because of the genetic bias towards ecdysteroid signaling that is suggested by our results.

Our study presents a genomic mapping approach that can be used on large-scale transcript data. This method can enhance functional genomic studies by adding spatial information as to where significantly expressed genes are located in the genome. Using this method, we found four genomic clusters of differentially expressed genes, suggesting that, in these regions of the genome, nearby genes are consistently and differentially expressed between High and Low genotypes (Table 1, supplementary material Fig. S3). Similar gene clusters are often co-regulated (Cohen et al., 2000; Spieth et al., 1991; Takadera et al., 1996). According to our conservative chromosome-wide thresholds, previously mapped QTL for worker behavior and ovary size do not contain putatively co-regulated gene clusters, as would be expected if the QTL led to changes in cis regulation of groups of genes. More quantitative measures of differential gene expression, i.e. with RNA sequencing and subsequent genome mapping, may better reveal genomic patterns of gene expression and the degree to which gene clusters are co-regulated. Alternatively, the identified QTL may affect ovariole number during ovary development in late larval stages. Because we used newly emerged bees in this study, we might have missed the developmental window when the QTL affect the ovariole number.

One of the clusters we detected contained MRJP-encoding genes. A few MRJP transcripts are generally expressed in honey bee reproductive tissues, suggesting roles in sex-specific reproductive maturity (Drapeau et al., 2006). We found the gene cluster to be upregulated in Low strain ovaries. Ecdysteroid signaling, which presumably is also upregulated in Low strain ovaries, enhances the royal-jelly-producing (and strongly MRJP-expressing) hypopharyngeal glands of workers in addition to influencing ovarian physiology (Wegener et al., 2009). Cross-fostering studies indicate that low strain nurse workers provide different nutritional environments to developing larvae than high strain nurse workers (Linksvayer et al. 2009b, 2011), perhaps as a result of differential MRJP1 expression. The elevated expression of MRJP genes in Low strain ovaries might, in other words, be an ecdysteroid-driven response. An MRPJ1 monomer influences honey bee development when secreted in larval food (Kamakura, 2011). MRJP1 was not differentially expressed in wild-type ovaries of different sizes (Fig. 3D), suggesting that it is not involved in associations between ovary size and behavior in worker bees. Yet, candidate gene expression is a relative quantity, an amount corrected towards housekeeper genes that are equally present per unit total RNA. Between individuals, however, the summed total RNA of variably sized tissues can correlate with organ size, e.g. if the larger-sized organs have more cells. Large ovaries have more cells that make up the ovariole count. This consideration could be relevant for translated proteins such as MRJP1, as larger ovaries might release more product than small ovaries.

Ovariole number plays a role in social behavior (Amdam et al., 2006; Wang et al., 2010) and ovarian physiology (Makert et al., 2006) in honey bees. However, High and Low pollen-hoarding genotypes can also be factors influencing the ovarian transcriptomes. Therefore, we tested the association between ovary size and the expression of five candidate genes in wild-type bees: HR46, ftz-f1, TYR, MRJP1 and JHi-26. For the two latter genes, no association was supported. For HR46 and ftz-f1, differences were significant but levels were opposite to those predicted by the average ovariole number of High and Low strain genotypes, suggesting that differential expression of these genes can be influenced by the genetic background. TYR, however, showed consistent relationships with ovary size between wild-type and the selected stocks.

It was suggested previously that selection for pollen-hoarding strains has influenced relationships between HR46, Ftz-f1 and ovary size in newly emerged bees (Wang et al., 2009). HR46 and ftz-f1 mRNA levels coincide with signaling by ecdysteroid hormones (Velarde et al., 2009), and ecdysteroid titers are shifted in late pupal stages and during early adult life between High and Low pollen-hoarding strain bees (Amdam et al., 2010). This shift may also change the induction of HR46 and ftz-f1 relative to ovary size when selected strains are compared with wild type. Similar phenomena are known from other systems (Toledo-Rodriguez et al., 2004).

TYR, which is activated by the neurotransmitter tyramine, was consistently negatively associated with ovariole number in the two selected strains and in wild-type bees. TYR is a positional candidate gene in the QTL Pln2 that is genetically linked to foraging behavior and ovary size (Hunt et al., 2007; Linksvayer et al., 2009b). Tyramine, furthermore, influences honey bee sucrose responsiveness, which diverges between pollen-hoarding strains, and predicts foraging behavioral traits in wild-type bees (Pankiw and Page, 2003; Scheiner et al., 2002). These connections might point to systemic patterns in TYR-dependent signaling. Yet, the modality of TYR expression differs between the brain and ovary of the pollen-hoarding strains (supplementary material Fig. S4), which suggests a more complex, tissue-specific function of TYR and tyramine that are supported by previous studies (e.g. Huang et al., 2004; Nykamp and Lange, 2000; Thompson et al., 2007). Taken together, our results support TYR as a candidate gene for associations between ovary size and behavior in worker bees.

Ovary size is affected by genetic and environmental factors and influences social behavior in honey bee workers. Social environmental signals (such as QMP) and physiological systems (such as the ecdysteroid and JH axes) regulate worker ovarian development and physiology. Yet, the worker ovary may also be an active endocrine signaling system and transcriptional feedback generator (Fig. 5). Therefore, our study supports the hypothesis that both ovaries and brains participate in the regulatory networks, which respond to the systemic stimuli such as hormones and pheromones, and modulate honey bee division of labor. These insights provide a platform for further functional studies that can determine how ovarian processes influence brain function and social behavior.

 
• Egf

  epidermal growth factor

 
• E75

  ecdysone-induced protein

 
• 75

  FDR false discovery rate

 
• ftz-f1

  ftz transcription factor-1

 
• GO

  gene ontology

 
• GPCR

  G protein-coupled receptor

 
• HR46

  hormone receptor-like in 46

 
• JH

  juvenile hormone

 
• JHi-26

  juvenile hormone-inducible protein 26

 
• LAO

  large ovary

 
• LSD

  least significant difference

 
• MRJP

  major royal jelly protein

 
• MRJP1

  major royal jelly protein 1

 
• OR

  olfactory receptor

 
• PDK1

  phosphoinositide-dependent kinase-1

 
• QMP

  queen mandibular pheromone

 
• QTL

  quantitative trait loci

 
• RQ

  relative quantification

 
• SMO

  small ovary

 
• TYR

  tyramine receptor gene

We thank N. Mutti and M. K. Fondrk for experimental assistance and we are also grateful to Z. Simões, C. S. Brent, A. Dolezal, K. Traynor and K. Dolezal for helpful comments on the manuscript.