Reproductive workers show queenlike gene expression in an intermediately eusocial insect, the buff‐tailed bumble bee Bombus terrestris

Bumble bees represent a taxon with an intermediate level of eusociality within Hymenoptera. The clear division of reproduction between a single founding queen and the largely sterile workers is characteristic for highly eusocial species, whereas the morphological similarity between the bumble bee queen and the workers is typical for more primitively eusocial hymenopterans. Also, unlike other highly eusocial hymenopterans, division of labour among worker subcastes is plastic and not predetermined by morphology or age. We conducted a differential expression analysis based on RNA‐seq data from 11 combinations of developmental stage and caste to investigate how a single genome can produce the distinct castes of queens, workers and males in the buff‐tailed bumble bee Bombus terrestris. Based on expression patterns, we found males to be the most distinct of all adult castes (2411 transcripts differentially expressed compared to nonreproductive workers). However, only relatively few transcripts were differentially expressed between males and workers during development (larvae: 71 and pupae: 162). This indicates the need for more distinct expression patterns to control behaviour and physiology in adults compared to those required to create different morphologies. Among female castes, reproductive workers and their nonreproductive sisters displayed differential expression in over ten times more transcripts compared to the differential expression found between reproductive workers and their mother queen. This suggests a strong shift towards a more queenlike behaviour and physiology when a worker becomes fertile. This contrasts with eusocial species where reproductive workers are more similar to nonreproductive workers than the queen.


Introduction
Eusociality, the division of adult females in reproductive queens and mainly sterile workers that care for the brood, has evolved multiple times independently within the Hymenoptera (bees, ants and wasps; Andersson 1984). The level of sociality varies within the Hymenoptera, ranging from nonsocial solitary species through primitively eusocial to highly eusocial taxa. Among highly eusocial hymenopterans, beside the clear division of reproduction between morphologically distinct workers and queens, further worker subcastes exist. These worker subcastes specialize in a particular set of tasks for a certain amount of time. Members of the subcastes may be responsible for, among others things, brood care, foraging or nest-guarding. In some ant groups, worker subcastes are morphologically distinct and display, at the least, a clear size polymorphism (Buckingham 1911;Detrain & Pasteels 1992). In other highly eusocial taxa, worker subcastes are monomorphic and task specialization is determined by age (Cameron 1989). In primitively eusocial taxa, such as the paper wasp Polistes, female adult castes are behaviourally distinct but monomorphic and behaviourally plastic, meaning an adult worker can potentially become the dominant, reproducing queen at any time by replacing the current queen or founding a new colony (Reeve et al. 2000;Sumner et al. 2006).
These distinct morphological and behavioural castes, which exist among adult females of a eusocial colony, are based on alternative expression of the same genome. The plasticity of the behavioural castes in the primitively eusocial paper wasp, Polistes canadensis, was demonstrated by the existence of overlapping gene expression patterns along a continuum from newly emerged females, through intermediate workers to the dominant queens (Sumner et al. 2006). Most gene expression studies in this area have, however, concentrated on highly eusocial taxa. Large differences in gene expression have been recorded both between the morphologically distinct queens and workers (Temnothorax longispinosus: Feldmeyer et al. 2014;Vespula squamosa: Hoffman & Goodisman 2007; Solenopsis invicta and S. richteri: Ometto et al. 2011; Apis mellifera: Grozinger et al. 2007) and between monomorphic, behavioural worker subcastes (T. longispinosus: Feldmeyer et al. 2014). The expression patterns of reproductive workers that lay unfertilized eggs later in a colony cycle become more 'queenlike' but they still remain more similar to nonreproductive workers than queens (Grozinger et al. 2007;Feldmeyer et al. 2014). Of the many genes found to be involved in caste differentiation, vitellogenin has perhaps received most attention and has been shown to be differentially expressed among female castes of the honeybee and several ant species (Amdam et al. 2003;Corona et al. 2013;Feldmeyer et al. 2014;Morandin et al. 2014). Often in such studies, a heavy focus has been placed on adult female castes; however, little work has been carried out to elucidate expression differences of males, but see Nipitwattanaphon et al. (2014). The haploid males are both morphologically and behaviourally distinct from their sisters and mother, but, although they differ in their ploidy level, they otherwise share the same genes as other colony members and are therefore also alternative expressions of the same genome.
Bumble bees represent an interesting taxon to study the phenomenon of eusociality as they possess both highly eusocial characteristics and more primitive features. For instance, whether a female will become a queen or a worker is irreversibly determined during development, as is the case for highly eusocial taxa. However, although a clear size dimorphism exists between queens and workers, generally both female adult castes are morphologically similar as in primitively eusocial species. Workers take on distinct tasks within a colony, but the division of labour is more plastic than is the case for higher eusocial bees and is generally not temporally fixed (Cameron 1989). Furthermore, towards the end of the colony cycle, the division of labour between workers and reproductive queens breaks down and queens and workers come into direct conflict over the parentage of males. At this stage, some workers activate their ovaries and begin to lay eggs and in the process become highly aggressive towards each other and also the queen (Bloch 1999;Alaux et al. 2004).
So far, no broad-scale studies have been conducted, which focus on the expression patterns involved in caste determination within bumble bees, although two previous studies did present some caste-specific genes (Pereboom et al. 2005;Colgan et al. 2011). Pereboom et al. (2005) investigated how and when females developed into queens or workers. They identified, using suppression subtractive hybridization, 12 genes whose expression differed in the comparisons: (i) worker and queen first-instar larvae; (ii) worker and queen fourthinstar larvae; (iii) adult queens and workers; (iv) reproductive and nonreproductive workers. Colgan et al. (2011), within their analysis of the bumble bee transcriptome, found a high number of transcripts (2185) that differ in their expression between adult castes, genders and developmental stages but considered their results as preliminary due to a lack of replication (1 larva, 1 pupa, 2 adult workers, 1 adult male and 1 virgin queen).
Here, using RNA-seq, we investigate genes involved in caste determination within the buff-tailed bumble bee, Bombus terrestris. We compare expression patterns of reproductive workers with those of nonreproductive workers and queens to isolate genes which are important for the acquisition of fertility as well as genes which may control behaviour differences compared to nonreproductive workers. Because of the flexible, plastic nature of bumble bee worker subcastes (Cameron 1989), reproductive workers are capable of becoming more 'queenlike' not only in their fertility but also in their behaviour. We therefore test the hypothesis that there is a greater similarity in gene expression patterns between queens and reproductive workers compared to those found in less plastic highly eusocial species.
Furthermore, we explore genes that control the specific behaviour and morphology of males. We investigate the question when, during the ontogeny of a male bumble bee, is the difference in gene expression to workers the greatest? Is the male gene expression pattern more distinct during larval development when the gonads and imaginal discs are generated? Are more genes involved in the development of the adult morphology during the pupal phase? Or does indeed the development and control of distinct behaviours among adults require the most distinct gene expression pattern? To address these questions, we compare gene expression patterns of males and workers both within larvae and pupae. In adults, we analyse differences in expression patterns between males, queens, reproductive workers and nonreproductive workers.

Colonies
Six young, commercially available Bombus terrestris audax colonies were obtained from Agralan Ltd. Initially, the colonies consisted of a mother queen and 8-20 workers. All colonies were kept in wooden nest boxes with the inner dimensions of 24 9 16 9 13.5 cm. The bees were supplied with pollen (mixed polyfloral pollen, www.naturallygreen.co.uk) and a sugar solution (BIOGLUC â ; Biobest) ab libitum. All colonies were kept in identical conditions within the same room at 26 C and 60% humidity in constant darkness.

Sampling
We aimed to collect samples from 11 different combinations of caste and developmental stage, each from three independent colonies. Within larvae and pupae, these were workers, males and queens, while in adults, we intended to collect males, reproductive workers, nonreproductive workers, mother queens and virgin queens.
Sampling was carried out under red light conditions. The gender of adults was determined by counting antennal segments (males: 13; females: 12) and checking for the presence or absence of a sting (Prys-Jones & Corbet 1987), while queens were identified via their superior mass (adult workers ranged from 26 to 325 mg, male adults from 127 to 347 mg and adult queens from 616 to 1191 mg). To identify reproductive adult workers, samples from each colony were anesthetized by cooling for approximately 10 minutes, and their abdomens were dissected to observe ovary development. To avoid loss of RNA, dissections lasted only a few seconds and samples were immediately snap-frozen in liquid nitrogen. Workers with developed ovaries were labelled 'reproductive'. The workers, in which ovaries were not visible, were categorized to be of 'undetermined reproductive status', because of the potential time lag between the expression of reproductive genes and subsequent changes in ovary morphology.
For the sampling of workers, queens and males during larval and pupal stages, the following protocol was followed. The colonies were photographed at regular intervals of one to two days to monitor the emergence of new batches and their development. With the term 'batch', we refer to a single cohort of offspring laid together. At intervals of at least three days, larvae and pupae were sampled from each batch while ensuring at least half of each batch was allowed to develop to adulthood. We collected larvae from each of the four larval instar stages based on their weight according to Cnaani et al. (1997) and assuming male instar masses were similar to worker instars. Pupae were collected both shortly after pupation (prepupae) and later in pupal development when appendages were developed.
Importantly, gender and caste of all sampled larvae and pupae were confirmed by isolating batches after pupation and sexing all emerging adults. Only if 100% of the unsampled adults emerging from a batch belonged to the same gender and caste would the samples from that batch be considered for analysis. All samples were collected between the hours of 9 am and 5 pm as soon as they became available. They were immediately weighed, snap-frozen in liquid nitrogen and then stored at À80 C.
Worker larvae and pupae were obtained from batches laid and reared in young colonies in the presence of the queen. After sufficient worker batches were available, the mother queen was removed from each colony for sampling. All batches laid in the presence of the queen but hatched shortly before or after the removal of the queen were considered potential queen batches (Pereboom et al. 2005). Any batches which were laid after queen removal were considered male batches. Additional male larvae and pupae were reared by isolating two to three groups of five workers from each colony in separate, small Perspex boxes containing pollen, sugar water and cat litter. The majority of male larvae and pupae (mean 67.2% AE 8.4% SEM) and all male adults were sampled from the main colonies.
As samples of the first larval stage were not obtained for workers from three separate colonies, L1 samples were excluded from all libraries. Adult virgin queens were only obtained from one colony, and the batches from which they emerged also produced adult workers. Therefore, larvae and pupae were only confirmed as queens if (i) they were sampled from batches from which adult queens emerged and (ii) if they exceeded 500mg (no sampled male or worker larva, pupa or adult exceeded 420 mg).

RNA extractions
Whole bodies were used for sampling for two reasons. First, we had no prior assumptions regarding the tissues, within which genes would be differentially expressed between castes, genders and developmental stages. Second, this would allow us to detect as many differentially expressed genes as possible across all comparisons. All samples were homogenized directly from À80 C. This was performed within the Eppendorf tubes with a plastic pestle for most larvae and in a ceramic mortar and pestle for large larvae, all pupae and all adults. The mortar was filled with liquid nitrogen to keep the samples frozen during homogenization. This was not necessary for the homogenizations which took place in Eppendorf tubes as the process was completed quickly. Total RNA was extracted from all samples using a GenElute Mammalian Total RNA Miniprep kit (Sigma-Aldrich) following the manufacturers' protocol. The quality and concentration of RNA were estimated with an Agilent 2100 Bioanalyzer.

RNA library construction
A total of 27 RNA libraries were constructed that covered all 11 combinations of caste and developmental stage from 1 or 3 colonies (Table 1). Based on the concentrations estimated with the Bioanalyzer, the larval libraries were prepared so as to contain equal quantities of RNA from each of the three larval stages 2-4 and equal quantities per individual within each larval stage. The same was also true for prepupae and pupae within the pupal libraries.

Sequencing and assembly
The 27 libraries were sequenced on three lanes of an Illumina HiSeq 2500 system in rapid mode at the Edinburgh Genomics facility of the University of Edinburgh. After quality control and raw read processing, the reads were mapped to the B. terrestris transcriptome, (BT_transcriptome_v2; Colgan et al. 2011), using bwa_0.6.1. Only reads which mapped uniquely were considered for further analysis. Counts per transcript were subsequently calculated for each library using custom scripts.

Differential expression analysis
The Blast2GO java program (Conesa et al. 2005) was used to annotate the transcriptome with gene descriptions and Gene Ontology (GO) terms (blastx against the NR database with e < 0.001). Differential expression analyses were carried out with the DESEQ package (1.16.0; Anders & Huber 2010) in R (3.1.1; Team 2012).
A neighbour-joining tree was created based on the expression differences between each of the 27 libraries. The distance matrix for the tree was calculated with the DESEQ package and contained euclidean distances between each library based on variance stabilization transformed counts. The tree was created with Phylip (3.695, Felsenstein 2005). A principal component analysis was performed on all adult libraries within the DESEQ package on variance stabilization transformed data. Euler diagrams were created with the R package venneuler (Wilkinson & Urbanek 2011).
Transcripts with a total of 50 reads or less across all 27 libraries were removed before performing the differential expression (DE) analyses. All remaining transcripts were tested for differential expression in each comparison. For each DE analysis, standard comparisons were performed between two conditions on normalized count data and with dispersion accounted for. Only transcripts with a Benjamini-Hochberg corrected P value (FDR) < 0.05 were considered as significantly, differentially expressed. For comparisons between castes within developmental stages, colonies were considered as replicates. No comparisons were made against queen pupae, queen larvae or adult virgin queens, as in each case, only one replicate existed. These libraries were, however, included in comparisons of expression between developmental stages.
Gene function enrichment analyses (Fisher's exact test) were carried out on DE transcripts with the R package topGO (2.16.0; Alexa & Rahnenfuhrer 2010). Enriched GO terms (FDR < 0.01) were subsequently summarized to meaningful clusters using Revigo (Supek et al. 2011). This method reduces redundancy of GO terms.

Assembly
A total of 469.3 million 50 base pair, single-end reads were generated, ranging from 13.9 to 23.7 million reads per library. The reads mapped to the Bombus terrestris transcriptome at an average of 85.27% (75.25-92.47%) per library. The transcripts ranged in length from 101 to 26 110 bases (mean: 1102; median: 721; Fig. 1). Average read depth across the 27 libraries ranged from 0 to 47 420

Overview of gene expression patterns
All replicates, that is libraries from the same caste and developmental stage but from different colonies, showed low variation in their gene expression patterns and thus grouped together well in a neighbour-joining tree (Fig. 2). The main clusters in the tree were formed by developmental stage (larvae, pupae and adults) rather than by caste. A differentiation in expression pattern between genders becomes more apparent in adults, where males form a distinct cluster. Within female adult castes, a further clustering seems to have occurred. All reproductive workers and mother queens clustered together, and two of the workers with undetermined reproductive status (W Au 8 and W Au 11) formed a separate branch, while the adult virgin queen remained more distant to all other female adult groups. The adult worker with undetermined reproductive status from colony 9 (W Au 9), on the other hand, grouped together with reproductive workers and mother queens. A principal component analysis (PCA) performed on all adult libraries indicated that W Au 9 was indeed reproductive although ovaries had not been visible (Fig. 3).
In the analysis, W Au 9 clusters strongly with all reproductive workers and mother queens. W Au 8 and W Au 11 form a distinct group, well separated from the reproductive workers and mother queens. This clustering pattern could mean that W Au 9 was reproductive and ovaries had not yet been visible in dissection. W Au 8 and W Au 11 were most probably nonreproductive. These conclusions were further supported by an overrepresentation of Apis mellifera reproductive genes (Cardoen et al. 2011) within W Au 9 but not in W Au 8 or W Au 11. The expression of the 299 genes, which were overexpressed in reproductive honeybee workers compared to nonreproductive workers, was lower in W Au 8 (median 122.8; mean 363.2) and W Au 11 (median 117.2; mean 389.4) than all three of our reproductive workers (W Ar 7: median 194.0, mean 638.9; W Ar 8: median 212.3, mean 627.6; W Ar 9: median 138.9, mean 457.0). These differences were significant compared to W Ar 7 (compared to W Au 8: P = 0.0045; compared to W Au 11: P = 0.0083) and W Ar 8 (compared to W Au 8: P = 0.0080; compared to W Au 11: P = 0.0018; Mann-Whitney U-test; Fig. 4). Expression of the Cardoen reproductive genes was  significantly higher in W Au 9 (median: 196.5; mean: 542.5) than in both W Au 8 and W Au 11 (P = 0.0238 and 0.0376, respectively; Mann-Whitney U-test; Fig. 4). For these reasons, W Au 9 was considered reproductive, and W Au 8 and W Au 11 were classed as nonreproductive for all further analyses.
The patterns shown in the neighbour-joining tree ( Fig. 2) and PCA of adult castes (Fig. 3) were reflected in the number of DE transcripts found between developmental stages and castes. From 6289 to 7483 (mean 7019) transcripts were differentially expressed between developmental stages. Only 71 and 162 DE transcripts were found between males and workers within larvae and pupae, respectively, while a mean of 4114 DE transcripts was found within adult comparisons ranging from 111 between reproductive workers and mother queens to 8706 between adult males and mother queens (Fig. 5).
For some analyses of differential expression, colonies were not uniformly distributed, for example adult males (colonies 7, 9 and 11) versus nonreproductive workers (colonies 8 and 11). For these, ANCOVAs were performed to test for significant colony effects (see Supporting information). In only one of nine cases was a significant main effect of colony found. There were five cases where a significant interaction between colony and expression was found. In each of these cases, a separate significant effect was still found for the important group difference: caste, developmental stage or gender. This means that any significant effects listed in the comparisons below can be attributed to differences in caste, gender or developmental stage rather than a colony effect.

Developmental stages
A total of 12 218 DE transcripts were recorded in the three comparisons between larvae, pupae and adults ( Fig. 6). As already suggested by the neighbour-joining tree ( Fig. 2), adults differed most greatly from the other two developmental stages, confirmed by 3237 transcripts which were differentially expressed compared to both pupae and larvae. A Gene Ontology (GO) term enrichment analysis showed that a heightened cellular metabolism distinguishes larvae from pupae and adults. The three main clusters of significantly overrepresented GO terms (Fisher's exact test, FDR < 0.01) in a Revigo treemap were 'translation', 'oxidative phosphorylation' and 'ribosomal biogenesis' (Fig. 7). Most overrepresented GO terms among transcripts upregulated in pupae are either related to cell communication and movement, 'signal transduction' and 'cellular component organization', or the development of morphological features, 'anatomical structure morphogenesis' (Fig. 8). Most enriched adult GO terms belonged to the supercluster 'G-protein-coupled receptor signalling pathway' (Fig. 9). This cluster included subclusters such as, 'phototransduction', 'detection of stimulus' and 'cell surface receptor signalling pathway', highlighting the higher sensory capabilities and requirements of adults. 42.6% of larval, 58.7% of pupal and 48.3% of adult DE transcripts either received no significant blast hit or were linked to genes of unknown function. Detailed results of the Fisher's tests can be found in the Supporting information. This test was repeated using lists of unique genes rather than transcript lists. All tendencies and the largest GO clusters remained unchanged. However, the number of significantly enriched GO terms was reduced. This was most probably due to the reduced number of genes in the test as a consequence of a high number of transcripts without an annotated gene match.

Male versus worker larvae
Within larvae, only a relatively small group of transcripts proved to be differentially expressed between males and workers (32 and 39 upregulated transcripts, respectively). Within the list of male larvae, DE transcripts nose resistant to fluoxetine protein 6-like, nrf-6, appeared six times with a fold change (FC) ranging from 3.86 to 25.34 and expression of 48-4576 mean normalized counts (mnc; Appendix S1, Supporting information). Nrf-6 is a transmembrane protein present in the intestine of various invertebrates (Choy & Thomas 1999;Yao et al. 2014) and has been reported as upregulated in the gut of Ostrinia nubilalis larvae (Lepidoptera) in response to a bacterial toxin (Yao et al. 2014). The presence of a further transcript within this list which encodes cytochrome p450 6k1-like (BTT39618_1; 2.07 FC; 1,682 mnc; Appendix S1, Supporting information) provides possible further evidence for an infection within the male larvae. Riddell et al. (2014) found in B. terrestris that the expression of 16 different cytochrome p450 transcripts was altered postinfection.
Takeout-like (XP_003397291.1; transcript BTT15842_1) was also strongly upregulated in male larvae compared to worker larvae (5.75 FC; 2693 mnc; Appendix S1, Supporting information). A close homolog to this transcript (blastp: 68% identity, e-value 3e À126 ) has been characterized for A. mellifera (Hagai et al. 2007). Takeout (to) was 0 2500 5000 7500 10 000 Differentially expressed transcripts within and between developmental stages. Darker colours: upregulated in first named caste; lighter colours: upregulated in second named caste. M = male; W = worker; MQ = mother queen; L = larva; P = pupa; A = adult; R/ N = reproductive/nonreproductive.  reported to be involved in the regulation of maturation in worker honeybees. In that study, only adult workers were investigated so that any gender or developmental effects are as yet unknown for Hymenoptera. However, the to gene family is known to be overexpressed in adult Drosophila males, affecting courtship behaviour (Dauwalder et al. 2002). The majority of the worker larvae DE genes (24 of 39; 8 of the top 10 in terms of FC) were either of unknown function or received no significant blast hits (Appendix S1, Supporting information). One vitellogenin transcript (BTT24408_1; 4.35 FC; 62 mnc; Appendix S1, Supporting information) was overexpressed in worker larvae compared to male larvae.

Male versus worker pupae
Differentiation was somewhat greater between males and workers during the pupal phase compared to the larval phase. A total of 128 transcripts were significantly upregulated in male pupae and 34 in worker pupae. The pupal list contained a high number of uncharacterized transcripts: 84 (66%) male and 24 (71%) worker pupae transcripts (Appendix S1, Supporting information).
Six male DE transcripts coded for tubulin-related genes (3 a-tubulin transcripts, 1 b-tubulin transcript and 2 tubulin-tyrosine ligases; 6.28-490.54 FC; 56-1200 mnc; Appendix S1, Supporting information). The tubulin-tyrosine ligase is involved in the post-transcriptional modification of a-tubulin (Ersfeld et al. 1993), so it appears tubulin transcripts, especially a, may be important for male pupal development. The same vitellogenin transcript upregulated in worker larvae (BTT24408_1) was also upregulated in worker pupae compared to male pupae (5.83 FC; 29 mnc; Appendix S1, Supporting information).

Fertility genes
Within the comparisons between adult castes (males, reproductive workers, nonreproductive workers and mother queens), reproductive workers and mother queens were most similar with only 111 DE transcripts (64 upregulated in reproductive workers and 47 in mother queens; Fig. 5). Nonreproductive workers, on the other hand, were distinct from both mother queens (2499 upregulated in nonreproductive workers, 2817 in mother queens) and, to a lesser extent, reproductive workers (844 upregulated in reproductive, 810 in nonreproductive workers). The majority (791, 93.7%) of the transcripts upregulated in reproductive workers compared to nonreproductive workers were also upregulated in mother queens compared to nonreproductive workers. As the common difference between nonreproductive and reproductive workers and between nonreproductive workers and mother queens is their fertility status, we have named these 791 transcripts 'fertility genes' (Fig. 10).
All differential expression values in this section are based on the comparison of reproductive and nonreproductive workers, although all transcripts were also upregulated in mother queens versus nonreproductive workers. A total of 267 (33.8%) of the fertility transcripts were of unknown function (1.78-336.18 FC; 18-39 927 mnc; Appendix S1, Supporting information). A large number of transcripts were involved in protein   144-491 mnc). Sixty-one of the fertility transcripts (1.77-9.53 FC; 12-9973 mnc) were direct homologs of genes upregulated in honeybee reproductive workers in a similar comparison (Cardoen et al. 2011). These transcripts encoded genes with functions such as oocyte meiosis, oocyte axis specification, oogenesis and female gonad development (Appendix S1, Supporting information).
The list also contained two vitellogenin (4.95 and 6.03 FC; 4103 and 111 595 mnc) and four vitellogenin   Fig. 11a,b). The vitellogenin transcripts (BTT24408_1 and BTT40935_1) are closely related to the 1772 amino acid vitellogenin genes ACQ91623 and ACU00433 of Bombus ignitus and Bombus hypocrita, respectively (Table 2). These genes correspond to the conventional Vg1 gene described by Morandin et al. (2014; blastp: E = 0.0, Id = 33%). The four receptor transcripts corresponded to the two B. terrestris genes vitellogenin receptor-like isoform 1 and isoform 2 (XP_003402703 and XP_003402704). Vitellogenin was, however, not restricted to female reproductive castes. The second highest expressed vitellogenin transcript across all libraries, BTT07410_1, constituted on average 28.1% of vitellogenin transcripts. This transcript together with three further transcripts (BTT35710_1, BTT41989_1 and BTT37349_1) is associated with the B. terrestris gene XP_003400264 (vitellogenin-6-like), which is 1,514 amino acids in length and corresponds with the Vg-like-A gene described by Morandin et al. (2014; blastp: E = 0.0, Id = 44%; Table 2). These four transcripts appear to be involved in development and independent of gender as they were upregulated in all larvae and pupae samples compared to adults irrespective of caste and gender (Fig. 11c). One vitellogenin transcript (BTT00708_1) was significantly upregulated in adults compared to pupae and larvae but was downregulated by reproductive workers (significantly compared to male adults) and mother queens (significantly compared to nonreproductive workers and male adults; Fig. 11d). This transcript is coded by the B. terrestris vitellogenin-like gene XP_003393940, which is much shorter than the two previously discussed vitellogenin genes (319 amino acids) and is similar to Vg-2 of A. mellifera (blastp: 66% identity, e-value 1e À142 ) and the Vg-C-like homolog described in Morandin et al. (2014; blastp: E = 1e À134 , Id = 57%; Table 2).
Seven a-glucosidase transcripts were differentially expressed within the fertility genes (6.93-9.14 FC; 16-35 260 mnc). An analysis of all 10 a-glucosidase transcripts within the B. terrestris transcriptome across all libraries showed raised expression levels for reproductive workers and mother queens compared to all other castes and developmental stages. Nonreproductive workers had the third highest levels of the 11 combinations of caste and developmental stage but a-glucosidase transcripts were eight times more abundant in reproductive workers and mother queens (Fig. 12). Four glucose dehydrogenase transcripts (BTT01220_1, BTT08099_1, BTT18258_1 and BTT20465_1), on the other hand, were downregulated in mother queens and reproductive workers, although upregulated in all adults compared to larvae and pupae (Fig. 13). These transcripts all related to the B. terrestris glucose dehydrogenase gene XP_003395668.1.
Interestingly, mean expression of the 10 a-glucosidase transcripts correlated significantly with mean expression of the two vitellogenin transcripts (BTT24408_1 and BTT40935_1), which were also upregulated in the fertility genes (q = 0.7247; P = 1.91 9 10 À5 ; Spearman's rho). Similarly, mean expression of the four glucose dehydrogenase transcripts, downregulated in fertility genes, significantly correlated with the downregulated vitellogenin transcript (BTT00708_1; q = 0.7888; P = 1.02 9 10 À5 ; Spearman's rho). Two transcripts (BTT20241_1 and BTT33633_1; 67.16 and Inf FC; 5528 and 37 mnc), which encode laccase-2like, were upregulated in reproductive versus nonreproductive workers but not in mother queens versus nonreproductive workers. Laccase 2 is a protein involved in the sclerotization of extracellular structures in invertebrates (Arakane et al. 2005).

Nonreproductive workers
For the majority (465 of 810; 57.4%) of the transcripts upregulated in nonreproductive workers compared to reproductive workers, the function was unknown (Appendix S1, Supporting information). Nineteen of the nonreproductive worker genes were direct homologs of genes upregulated in nonreproductive A. mellifera workers (Cardoen et al. 2011). Eight of those (1.77-3.74 FC; 72-653 mnc) had been attributed to the effect of the queen mandibular pheromone (QMP) in a previous study (Grozinger et al. 2003; Appendix S1, Supporting information).

Adult queens
Transcripts, which were upregulated in mother queens compared to both reproductive and nonreproductive workers, were considered 'queen genes' (Fig. 10). The 40 queen transcripts ranged in fold change compared to reproductive workers from 1.68 to 8.87 (29-245 472 mnc; Appendix S1, Supporting information). Eleven of the transcripts (27.5%) were of unknown function. Most notable among the queen genes were 5 transcripts relating to serine protease inhibitors, SPI (2.92-8.87 FC; 2145-10 596 mnc). These five SPIs were expressed together at a mean of 27 758 mnc AE 1247 SEM in mother queens compared to only 5026 mnc in the virgin queen ( Fig. 14; Appendix S1, Supporting information). The second highest levels were found in nonreproductive workers (7555 mnc AE 1527 SEM) followed by reproductive workers (6435 mnc AE 699 SEM).

Comparison with previous studies on Bombus terrestris
The top 10 transcripts upregulated in larvae in the study carried out by Colgan et al. (2011) related to cuticle proteins, the storage protein hexamerin and the metabolic proteins carbonic anhydrase and cytochrome p450. In this study, 5 cuticle, 2 hexamerin (70c and 70b), 10 carbonic anhydrase and 12 cytochrome p450related transcripts were also upregulated in larvae compared to pupae and adults. The 10 transcripts listed in Colgan et al. could be linked to one GO term (GO:0042302: 'structural constituent of cuticle'), which was also attributed to 17 of the larvae transcripts (upregulated relative to pupae and adults). In a further study, a cuticle protein and hexamerin were also present in larvae but absent in adults; pupae were not included in the analysis (Pereboom et al. 2005). The vitellogenin transcript BTT07410_1, which we found to be upregulated in larvae and pupae, was also overexpressed in pupae in the Colgan et al. study (2011), however, not detected in larvae. All 7 of the GO terms, which were associated with the top 10 pupal genes in the Colgan et al. study (2011), were also present in our list of upregulated pupae transcripts (P = 1.3 9 10 À4 , hypergeometric test).
In workers, Colgan et al. (2011) found overexpressed genes associated with flight, defence and metabolism (cytochrome p450, lipase and a-glucosidase). In this study, flight muscles were also overrepresented in nonreproductive workers and the metabolism genes lipase, cytochrome p450 and a-glucosidase were more highly expressed in workers than in males. Twenty of the 36 GO terms associated with the worker transcripts in the Colgan et al. study were also found in the transcripts upregulated in nonreproductive workers relative to adult males in the current study (P = 5.8 9 10 À9 , hypergeometric test). The genes differentially expressed between adult female castes and subcastes in the Pereboom et al. study (2005), 60-S ribosomal protein, chymotrypsin, cytochrome oxidase, peroxiredoxin, fatty acyl CoA-desaturase and ATP synthase beta subunit, could not be confirmed with our data. Colgan et al. (2011) found transcripts of the flight muscle gene titin to be overrepresented in male adults, as well as several immunity genes. Many flight muscle proteins were also upregulated in our study; however, we could not confirm the overrepresentation of immunity genes among the transcripts with known function. All 17 GO terms present in the top 10 male transcripts of the Colgan et al. study (2011) were also present in our list of male transcripts (upregulated in adult males relative to nonreproductive workers; P = 2.3 9 10 À14 , hypergeometric test).

Discussion
We compared gene expression patterns both between developmental stages and between castes within each developmental stage for the buff-tailed bumble bee Bombus terrestris. The number of differentially expressed transcripts ranged from 71 between male and worker larvae to 8706 between adult males and mother queens. We found gene expression patterns to differ more between developmental stages than between caste or gender. Genes upregulated in larvae were associated with a high cellular metabolism, whereas in pupae overexpressed genes were associated with cell communication and the development of morphological features. Most of the overrepresented GO terms in adults were related to the G-protein-coupled receptor signalling pathway. G-proteins are cell surface receptors, which respond to extracellular stimulants with an intracellular signal cascade (Strader et al. 1994;Dohlman 2002).
The number of genes differentially expressed became progressively larger through the three developmental stages as each caste became more distinct. These findings suggest a comparatively low number of genes are required to create distinct morphological castes compared to the high number involved in distinct behaviours between adult castes. Gender grouped more strongly than caste as expression was less variable between adult males than within each of the female castes. Similar findings have been presented for the social wasp Vespula squamosa, for which workers, queens and males clustered clearly into developmental stages (Hoffman & Goodisman 2007). A study on the two fire ant species Solenopsis invicta and S. richteri also found expression patterns between developmental stages to differ more than between gender followed by caste and species (Ometto et al. 2011).
Our data confirmed, to some extent, previous findings for B. terrestris (Pereboom et al. 2005;Colgan et al. 2011). Several associations of gene functions with specific castes or developmental stages detected by Colgan et al. (2011) were also found in the present study. Discrepancies can be explained by, in contrast to our study, a lack of replication in the 2011 study or a difference in analysis structure; Colgan et al. (2011) implemented R-STAT (Stekel et al. 2000) to calculate differential expression of a contig within all libraries, whereas we performed specific pairwise comparisons.
Little overlap could be found with an older study on caste determination in B. terrestris (Pereboom et al. 2005). However, due to the method implemented in that study, suppression subtractive hybridization, only a few differentially expressed genes could be isolated, and also, due to a different focus, fewer comparisons were performed than in our study (Pereboom et al. 2005).

Reproductive workers closely resemble queens
Towards the end of a bumble bee colony cycle a queen-worker conflict develops, in which reproductive workers compete with the mother queen for male parentage (Bloch 1999;Alaux et al. 2004). The expression patterns observed in this study support our hypothesis that when bumble bee workers become reproductive they would, in comparison with highly eusocial species, more strongly resemble queens in their behaviour and physiology due to the more plastic nature of worker castes in bumble bees. Of all adult expression patterns, those of reproductive workers and mother queens were most similar, in fact more similar than between reproductive and nonreproductive workers. Only 111 transcripts differed significantly between reproductive workers and mother queens compared to 1654 between reproductive and nonreproductive workers. Nonreproductive workers differed from mother queens even more strongly (5316 DE transcripts). These findings are in strong contrast to patterns found in two highly eusocial hymenopteran species. In Apis mellifera, over 2000 genes differed significantly in both comparisons between queens and either reproductive or nonreproductive workers; the expression of only 221 genes differed significantly between the two worker castes (Grozinger et al. 2007). Similarly, 2785 genes were significantly up-or downregulated between queens and reproductive workers in the myrmicine ant Temnothorax longispinosus compared to only 571 between reproductive and nonreproductive workers (Feldmeyer et al. 2014). Feldmeyer et al. (2014) suggested the high similarity between reproductive and nonreproductive workers in these two hymenopteran taxa indicates that a relatively low number of genes are required for ovary activation and egg laying compared to the high number involved in further physiological or behavioural differences which exist between queens and workers. Based on this assumption, our data indicate a greater similarity in behaviour and general physiology between bumble bee queens and reproductive workers than is the case for honeybees or myrmicine ants. The division of labour among bumble bees is not as clearly temporally or morphologically fixed as in the highly eusocial honeybees and most ants, indicating the capability of individual bumble bee workers to flexibly adapt their current role (e.g. from forager to nurse) to changing conditions within a colony at any given time (Cameron 1989). In honeybees, a shift towards a more 'queenlike' expression pattern was recorded in reproductive workers (Grozinger et al. 2007), but it is possible that the more flexible nature of the bumble bee worker roles in our study allowed a much stronger shift in behaviour and physiology, allowing the reproductive workers to more strongly resemble a queen.

Male expression patterns are most distinct among adults
Males, in contrast to both queens and all workers, do not possess a sting and their antennae contain an additional segment. Their sexual organs naturally also differ. It was therefore surprising that the expression of comparatively few transcripts significantly differed during development. Within the larval stage, no clear clusters could be formed based on expression patterns, and only 71 transcripts differed significantly in their expression levels between males and workers. During the pupal stage, when morphological features are being generated, expression patterns became more distinct with 162 transcripts differentially expressed. However, it was only in adulthood that the expression pattern of males became truly distinct from all other castes. In male adults, between 2411 and 8706 transcripts were either up-or downregulated compared to the three adult female castes mother queen, reproductive worker and nonreproductive worker. This indicates that a much greater number of genes may be required to control behaviour and the physiology of reproduction than to develop morphologies.
A high number (69; 59.5%) of the pupal transcripts upregulated in male pupae, and therefore likely to contain some genes linked to the development of the male morphology, were of unknown function. The six aand b-tubulin transcripts, which were overrepresented in male pupae, are possibly linked to spermatogenesis as both a2and b-tubulin are known to be testis specific in Drosophila (Kemphues et al. 1979;Theurkauf et al. 1986). One hundred and ninety transcripts were involved in mitochondrial processes, and a further 37 were associated with genes linked to muscle development. These 37 transcripts related to the proteins myosin, troponin, twitchin and titin, which are all integral parts of insect muscles (Hooper & Thuma 2005). In their mating flights, males have been recorded as covering significantly larger distances than workers from the same colony (Kraus et al. 2009). The apparent greater need for muscle development and higher energy levels in males compared to workers are possibly linked to their greater flight distances.

Vitellogenin
Vitellogenin was originally thought to be limited to reproductive egg-laying females due to its function as a yolk precursor in all oviparous animals, although it is now known to fulfil various functions in hymenopterans (Amdam et al. 2003). The reproductive ground plan model proposed by Amdam et al. (2004) describes how pleiotropic associations of reproductive genes, above all vitellogenin, with genes that control sensory perception, longevity and foraging behaviour have been utilized to control behaviour patterns in honeybee worker subcastes.
Previously, only one vitellogenin gene had been described for honeybees, which is differentially expressed in female castes (Amdam et al. 2012). However, in a more recent study on Formica ants, four vitellogenin homologs were found within the genome of all ant and bee species included in the study (Morandin et al. 2014). These vitellogenin homologs were classed as conventional vitellogenin (Vg-1), Vg-like-A, Vg-like-B and Vg-like-C, which were expressed at different levels and differently between queens and workers. Four copies of Vg-1 have been found in S. invicta (Wurm et al. 2011) and T. longispinosus (Feldmeyer et al. 2014) and two in Pogonomyrmex barbatus (Corona et al. 2013). In each of these cases, the gene copies showed differential expression between adult female castes.
Here, we have found only one copy of Vg-1 and two further vitellogenin genes which are closely related to Vg-like-A and Vg-like-C. Vg-1, as in Formica adults (Morandin et al. 2014), was the highest expressed of the three vitellogenin genes discovered in this study. We found Vg-1 to be highly upregulated in mother queens and reproductive workers compared to all other castes and developmental stages, which suggests it has maintained its conventional function in reproductive egg-laying females for B. terrestris. This also appeared to be the case for 3 of 7 Formica species, in which Vg-1 was upregulated in queens compared to workers (Morandin et al. 2014). Workers were not grouped according to reproductive status in the Morandin et al. study (2014), which could explain the lack of significant differences between castes in more than three species. The comparison of expression between queens and workers for Vg-A-like differed among Formica species (upregulated in queens for 3 and in workers for 1 species; Morandin et al. 2014). In B. terrestris, the homolog of Vg-A-like, XP_003400264, appears to play a lesser role in adults, as it was upregulated in larvae and pupae of both genders compared to adults. Expression of Vg-C-like was significantly higher in workers than queens in all 7 Formica species (Morandin et al. 2014). In the current study, the homolog of Vg-C-like, XP_003393940, was also downregulated in mother queens but also in reproductive workers compared to higher levels in nonreproductive workers and adult males.
Here, we have shown that three copies of vitellogenin genes are not only differentially expressed between adult females castes as shown for other hymenopteran taxa (Amdam et al. 2004;Wurm et al. 2011;Corona et al. 2013;Feldmeyer et al. 2014;Morandin et al. 2014), but that they are differentially expressed across all adult castes and between developmental stages.

Carbohydrate processing enzymes
We found the expression of two carbohydrate processing enzymes to be differentially expressed among adult castes. Expression of a-glucosidase was almost exclusively restricted to female adults but with levels eight times higher in mother queens and reproductive workers than in nonreproductive workers. This is in contrast to honeybees for which a-glucosidase is downregulated in reproductive compared to nonreproductive honeybee workers (Cardoen et al. 2011). In honeybees, a-glucosidase catalyses the splitting of the sucrose present in nectar in the production of honey (Ohashi et al. 1999;Kubota et al. 2004). The apparent restriction of this protein to reproductive workers and mother queens may indicate a different role for this protein in B. terrestris compared to honeybees. Glucose dehydrogenase, on the other hand, was present in all B. terrestris adults but was downregulated in reproductive workers and mother queens. The similar protein glucose oxidase is specifically found in the hypopharyngeal gland of forager honeybees and converts the glucose of nectar to gluconic acid and hydrogen peroxide in honey production (Ohashi et al. 1999). Glucose dehydrogenase may perform a similar function in B. terrestris as it also catalyses the oxidation of glucose to gluconic acid but without the by-product hydrogen peroxide (Bak 1967). Expression of a-glucosidase significantly correlated positively with Vg-1, while expression patterns of glucose dehydrogenase significantly correlated positively with Vg-C-like. These correlations indicate interactions between vitellogenin and the two carbohydrate enzymes, which may be associated with distinct foraging preferences among adult castes.

Further caste-specific genes
One highly represented gene in the list of transcripts overexpressed in mother queens compared to reproductive workers was serine protease inhibitor. Serine prote-ases have been detected in the venom of a variety of Hymenoptera species (Hoffman & Jacobson 1996;Winningham et al. 2004). One possibility is that serine protease inhibitor was produced to counteract the effect of stings, either as a reaction to sting attacks or as a preventative measure. This could be linked to the high aggression shown towards a bumble bee queen by workers late in a colony cycle often resulting in her death (Bourke & Ratnieks 2001).
Workers can become reproductive in queenright conditions, but whether workers or queens control worker reproduction is unresolved (Alaux et al. 2007). Intriguingly, we found eight transcripts upregulated in nonreproductive individuals (BTT06229 1, BTT09963 1, BTT20486 1, BTT15870 1, BTT22989 1, BTT27276 1, BTT17949 1 and BTT09790 1) whose expression is believed to be regulated by queen mandibular pheromone in A. mellifera and where expression shows similar patterns (Grozinger et al. 2003;Cardoen et al. 2011). It is clear that further research is needed to understand the relationship between pheromonal signalling and ovary development (Amsalem et al. 2009).
In each of the caste comparisons performed in this study, large numbers of differentially expressed transcripts either could not be associated with any known gene or were related to genes with so far unknown function. These range from 1636 to 2609 (32.0-54.4%) upregulated transcripts when comparing between developmental stages. The number of differentially expressed transcripts was much lower between male and worker larvae (34 and 39) and pupae (128 and 34), but still, the majority of these transcripts (58.7%) were of unknown function. A total of 267 of the 791 fertility transcripts, that is upregulated in reproductive workers and mother queens compared to nonreproductive workers, belonged to uncharacterized genes, while 465 and 526 transcripts in the comparison between nonreproductive workers and adult males were of unknown function. Clearly, further research is required in these areas.

Conclusions
We conducted the first large-scale RNA-seq analysis into caste differentiation within the genus Bombus, for which eusociality can be considered intermediate between that found in primitively eusocial taxa such as the paper wasp and highly eusocial species such as the honeybee or most ants. As in other similar studies on eusocial hymenopterans, a high number of genes were differentially expressed in all comparisons between castes, genders and developmental stages. Significant overlaps with analyses on higher eusocial taxa exist in terms of overall expression patterns as well as specific genes. One striking difference between Bombus terrestris and higher eusocial hymenopterans is how much more closely a bumble bee reproductive worker resembles the queen regarding its gene expression. Further research may be able to determine whether this finding is restricted to B. terrestris or if it is linked to the more plastic nature of worker subcastes in bumble bee taxa in general. The annotation of many unknown genes, which were differentially expressed in our analysis, and further research on B. terrestris following the imminent release of the genome will help us to better understand how distinct castes are created, maintained or altered within this important species.