BACKGROUND
caste development in honeybees  is a remarkable example of, and model for, developmental plasticity.
in this species, two adult female phenotypes are produced by one genome in response to diet  <cit> .
larvae destined to become queens are fed royal jelly  which is rich in carbohydrates and lipids  <cit> .
rj also contains major royal jelly proteins, one of which, royalactin, is vital for queen development  <cit> .
as a result of a higher food intake  <cit>  and a more nutritious food, queen bees are larger, longer-lived, have fully developed ovaries and differ in their behaviour and morphology, when compared to worker bees.
transferring larvae from worker cells to queen cells  during early larval development causes bees to develop some queen-like phenotypes.
based on these experiments it has been proposed that differentiation between the two castes begins on the first day of larval development and is progressive  <cit> .
manipulating larval diet in this way can induce queen-like properties in worker bees up to 60 hours of larval development  <cit>  demonstrating, that at least during early larval development, there is plasticity in the developmental pathways that give rise to queen and workers.
the molecular mechanisms by which rj triggers queen development are only beginning to be understood, but rj intake causes differences in juvenile hormone  titre between worker and queen larvae  <cit> .
these jh levels regulate gene expression  <cit> .
it seems likely, therefore, that a complex set of molecular mechanisms link rj exposure to changes in jh titre.
such mechanisms should be reflected in differences in gene expression between developing larvae destined for different castes.
differences in gene expression underlying queen and worker developmental trajectories have been investigated in a number of previous studies  <cit> .
these studies used a variety of techniques to profile gene expression, but have focused on mid-to-late larval development as few differences in gene expression have been identified before this stage  <cit> .
mid to late larval development includes the period , when jh titres differ between castes  <cit> , and so if is difficult to determine if these gene expression changes are due to rj exposure, or the action of jh.
while a number of mechanisms that regulate gene expression during caste development have been proposed, dna methylation is thought to be critical, as it can link environmental exposure to changes in gene expression  <cit> .
rna interference targeting the de-novo dna methyltransferase  <dig>  biases development towards the queen trajectory  <cit> .
honeybees, similar to other invertebrates, have relatively low overall levels of cpg methylation restricted to gene bodies  <cit> .
methylation within the gene body appears to influence alternative splicing in the honeybee  <cit>  and the ant  <cit> .
it is possible that the role of dna methylation in caste development may be via influencing alternative splicing of transcripts.
evidence from mammalian models indicates that gene body methylation may have other molecular functions, such as promoting rna polymerase pausing  <cit> , modulating other epigenetic marks  <cit> , or affecting altering promoter usage  <cit> .
it has been proposed that plasticity, such as caste development in eusocial insects and neuronal plasticity in invertebrates, is associated with dynamic changes in dna methylation  <cit>  and, by inference, alternative splicing.
dna methylation is likely to be one of many molecular mechanisms involved in caste development, but these mechanisms as a whole will lead to changes in gene expression, detectable using microarray techniques, or alternative splice isoforms detectable using transcriptomics.
we have examined gene expression throughout larval development, both before and after the wave of jh, in queen and worker castes, carrying out the widest survey of gene expression during caste development in this remarkable animal.
we used this to build a more complete understanding of how a honeybee larva responds to rj, and how this biases them into an alternative developmental trajectory, culminating in a queen bee.
RESULTS
experiment design
there are five stages of larval development  in the honeybee.
queen and worker samples were taken at the following time points : 6 hours , 12 hours , 36 hours , 60 hours , 84 hours , 108 hours  and 132 hours .
at each time point four replicate queen and worker samples were collected.
each replicate sample collected at 6 hours contained  <dig> larvae and samples collected at all other time points contained  <dig> larvae.
queen development was induced by grafting larvae  into artificial queen cells before returning them to the hive.
grafting, which is a standard apicultural technique, is carried out as soon as practical after the larvae hatch .
larvae in these queen cells are fed rj and this triggers queen development.
when collecting 6-hour larvae it was noted that rj was present in these cells.
worker larvae were not grafted.
gene expression data was generated from these samples and using custom two colour long-oligonucleotide microarrays .
at 60 hours, two additional biological replicates, each consisting of either  <dig> queen or worker larvae, were taken and transcriptome data was generated using high throughput sequencing .
the 60 hours time point was selected for hts as this time point has been highlighted as the earliest difference in gene expression seen a previous study  <cit> .
hts sequencing also gave us the opportunity to assess gene expression in greater depth, as well as analyse the effect of caste development on alternative splicing.
differential gene expression during caste development
microarray analysis indicates that queens and workers have relatively equal numbers of differentially expressed genes  throughout larval development .
hts of rna, however, found that 83 % of degs had higher expression in queens, similar to previous results  <cit> .
many studies have reported finding more differentially expressed genes using hts than using microarrays  <cit> .
this is at least partly because hts analysis is not hindered by the acquisition bias that can affect microarray analysis.
microarrays are sensitive to saturation effects and are also not able to detect transcripts that are expressed at a relatively low level due to high amounts of background fluorescence obscuring low hybridization signal to the probes.
we believe the latter phenomena accounts for the relatively high number of degs detected in queens by hts compared to microarray, as these genes had a median rpkm of  <dig>  compared with a median rpkm of  <dig>  for genes that are more highly expressed in worker larvae.
we believe the relatively low expression of these genes precluded them from being detected as differentially expressed by microarray analysis.
we have carried out extensive validation of these data sets using rt-qpcr  which indicates that both the microarray and hts experiments are producing high-quality, biologically relevant information.
genes we identified as being differentially expressed during caste development were assigned drosophila orthologs using blastx  <cit> .
approximately 90 % of the genes with higher expression in queen larvae have orthologs in drosophila as compared with 77 % of the genes with higher expression in workers, a statistically significant difference  .
barchuk et al.
also found a tendency for worker larvae to have higher expression of genes that lack orthologs in drosophila <cit> .
these data imply that genes more highly expressed by queen larvae have been conserved across ~ <dig> million years of evolution.
in contrast, genes more highly expressed in worker larvae are not as well conserved in drosophila and may be rapidly evolving.
differences in gene expression between queen and worker larvae are observed as early as 6 hours  post-grafting .
indeed, the largest numbers of degs identified by microarray analysis are detected at this time point .
while it is possible that the large number of differentially expressed genes at this time point are due to perturbation of the larvae during grafting, the appearance of genes that are differentially expressed at 6 hours and that retain higher expression in one caste during larval development implies that at least some of this differential gene expression represents an early response to rj.
we also do not observe gene ontology categories associated with stress response or physiological perturbation in the genes that are differentially expressed at 6 hours, indicating that differences we observe in expression at this early time-point are likely an early response to rj.
in queen larvae, several genes with increased expression early in larval development are maintained at high levels throughout larval development.
these genes are mitochondrial cytochrome c, phosphoenolpyruvate carboxykinase , phytanoyl-coa dioxygenase domain-containing protein  <dig> homolog and glycine n-methyltransferase-like .
worker larvae, in contrast, have higher levels of the apoptosis regulator bcl- <dig> throughout development.
these genes have functions that are difficult to interpret in terms of caste development, but are consistent markers of either worker or queen development.
given that queen development can be triggered in larvae by grafting until larval day three, well after the 6 hour time point by which we see substantial differences in gene expression, these early differences in gene expression must not commit the larvae to either a queen or worker fate.
we therefore propose these early changes in gene expression, in response to rj, may bias development toward one caste or the other, but that later gene expression can reverse that bias.
these early differences in gene expression may be attributable to the larval diets which differ, with the queen diet initially containing more sugar  <cit>  and less protein  <cit>  than worker jelly  though it is not clear if these differences occur at or before 6 hours after grafting.
these dietary differences have also been linked to morphological differences between the castes that are apparent as early as 2 days into larval development  <cit> .
the idea of an early, biased, phase of caste development is consistent with early larval grafting experiments which showed that worker larvae shifted into queen cells at ~12 hours of larval development result in bees that are queen-like, but retain some worker characteristics  <cit> .
a two-stage process of caste development has been previously proposed based on grafting, nutrition and physiological data  <cit> , but our data provides the first molecular evidence to support an early biased phase of gene expression and development followed by a commitment phase.
differential expression of hormone biosynthesis genes
one key event during larval development is the wave of jh synthesis in queen larvae, during l <dig> and l <dig> stages  <cit> .
juvenile hormone  is implicated in caste development in honeybees  <cit> .
our data indicates that several genes involved in jh and ecdysteroid biosynthesis are differentially expressed during larval development between castes  at, or around this jh wave.
queen destined larvae, have higher expression of eight cyp genes, many with roles in hormone biosynthesis.
cyp315a <dig>  acts in ecdysteroid biosynthesis  <cit> , and has higher expression at six hours of larval development.
cyp15a <dig>  which has higher expression in queens throughout larval development but peaking at 84 hours, acts as a jh epoxidase, catalysing the final step of jh biosynthesis  <cit> .
queens also have higher expression of jh methyltransferase, involved in the final steps of jh synthesis, with differential expression peaking at 84 hours after grafting .
in contrast, during mid larval development, workers have higher expression of jh esterase, which degrades jh in honeybees  <cit> .
seven genes induced in drosophila by jh titre  <cit>  are dynamically regulated, with three having higher expression in workers early in development and five having higher expression in queens during mid to late development, temporally correlating with the peak of jh titre during development  <cit> .
the trends in hormone-related gene expression in queen larvae, peaking as they do during mid-late larval development, are consistent with these genes either being involved in the synthesis of jh associated with the wave between the third and fourth larval instars, or with responding to the jh wave.
at this point we also observe the appearance of differential expression of a number of genes in both queen and worker larvae that remain stably differentially expressed throughout the remainder of larval development .
this change to a committed developmental trajectory, however, does not alter or reprogram the expression of genes we find differentially expressed at early stages.
our data is consistent with an initial bias towards one caste or the other induced quite early  by rj consumption, followed by a consolidation of that fate and final specification through circulating hormones.
differential expression of cytochrome p <dig> genes
cytochrome p <dig>  genes  which act in hormone and pheromone biosynthesis  <cit>  and the detoxification of foreign compounds  <cit>  show differential regulation between castes.
insect cyp genes are divided into four major clans based on sequence similarity  <cit> .
worker larvae have higher expression of  <dig> cyp genes at time points throughout larval development, the majority belonging to clan three, in particular sub-families six and nine.
there is evidence for rapid expansion of these two sub-families in the honeybee through recent gene duplication events  <cit> .
genes in cyp clan three are loosely classified as environmental response genes  <cit> .
the differential expression of cyp clan three genes is most marked at 60 hours of larval development, where worker larvae have higher expression of  <dig> genes of this clan.
this period coincides with a change in the worker larval diet which results in the introduction of pollen grains, whereas rj only contains trace levels of pollen  <cit> .
the introduction of pollen exposes worker larvae to a number of foreign compounds.
one example of this is the addition of quercetin, a flavinoid compound found in honey and pollen  <cit> .
four cyp genes are known to metabolise quercetin in the honeybee   <cit>  and all four have higher expression in workers during mid-larval development.
some cyp genes are involved in pheromone production  <cit>  and possibly caste specific synthesis of mandibular acids  <cit> .
cyp genes have also been linked to termite caste development  <cit> .
functional classes of genes differentially regulated between queen and worker larvae
two other families of genes were differentially expressed during caste development .
queens have higher expression of five genes involved in the response to reactive oxygen species , especially during mid larval development.
these genes include catalase, several thioredoxins and a peroxiredoxin.
the generation of ros by the mitochondria has been proposed to be a key mechanism in the aging process  <cit> .
this may reflect an increased respiration rate previously observed in queen larvae  <cit> .
the differential regulation of these genes may influence expression of vitellogenin, important in protecting honeybees from oxidative stress  <cit> .
in contrast workers have higher expression of three glutathione-s-transferase genes.
glutathione-s-transferases are important for detoxification of both endogenous  and exogenous  compounds  <cit> .
programmed cell death has a role in reducing ovariole numbers in larval worker ovaries during fourth and fifth larval instars  <cit> .
consistent with this we see higher expression of eleven genes involved in programmed cell death in worker larvae.
these belong to two major programmed cell death pathways - apoptosis and autophagic cell death.
hexamerin70b and caste development
hexamerin genes were differentially regulated between castes at several larval time-points.
rt-qpcr  <cit> , which is generally considered to be the gold-standard method to measure gene expression  <cit> , was used to confirm that, at 60 hours of larval development, workers have between  <dig> and  <dig> fold higher expression of the four hexamerin genes when compared to queens .
hexamerins are a family of insect amino acid storage proteins that have evolved from the copper-containing hemocyanins of ancestral aquatic insects and crustaceans  <cit> .
a role for the hexamerins has been suggested in caste development in termites, where two hexamerin genes have been shown to facilitate the jh dependent worker to soldier differentiation  <cit> .
several previous studies have found differences in expression of hexamerin genes between queens and workers implying a possible role for these proteins in caste development  <cit> .
hexamerin 70b  is expressed in the fat body and gonads during larval development and its expression, which is regulated by juvenile hormone  <cit> , peaks prior to the expression of other hexamerins during larval development  <cit> .
thus far no biological function has been attributed to hex 70b; it is found in the cytoplasm, nucleus, and in punctate dots in the cytoplasm of cells within the fat body  <cit> , implying that this protein has multiple biological roles, as reported for hex 70a  <cit> .
hexamerin70b knock-down
of the four hexamerins assayed by rt-qpcr  hex 70b showed the smallest fold change between queen and worker larvae .
this comparatively small fold-change provides us with the best opportunity, amongst the hexamerins, to knock-down expression of this gene to queen larvae-like levels, even if rnai treatment is inefficient and knock-down is weak.
given the technical challenge of rnai in honeybee larvae, and the high expression levels of all the hexamerin genes, we targeted hex 70b because moderate knockdown of this gene was more likely to produce queen-like expression levels, and give a clear, measurable, phenotype.
rnai knockdown of hex 70b resulted in significant changes  in the proportion of queen-like individuals produced when compared to control injected larvae .
emerging individuals were classified as queen-like based on the morphology of several key structures .
while hex 70b rnai produced a statistically significant change in phenotype, we did not see a significant difference in the expression of hex 70b between target and control larvae as determined by rt-qpcr , although we did observe a trend towards lower hex 70b expression in the hex 70b rnai injected group relative to controls.
the lack of statistical significance likely reflects the small sample sizes used for rt-qpcr, and large variation in expression of this gene in injected individuals.
samples were taken  <dig> – 18 hours post-injection for this analysis and it is possible that sampling too early, or too late, may explain why the difference in expression between the target and control groups was not significant.
knockdown of hex 70b biases individuals towards a queen fate , consistent with the high expression of these genes in worker larvae .
we propose that hex 70b has a key role in caste development at a time when we propose a commitment to worker or queen fate is set .
this implies that early biases in gene expression between worker and queen larvae can be overcome to change the developmental trajectory of a larva.
the possibility that hexamerins are acting to regulate jh activity  <cit>  raises the possibility that we are disrupting the final commitment phase of caste development by reducing hex 70b expression, rather than changing the bias in commitment.
gene ontology and pathway analysis
degs involved in caste development were categorized by proposed gene function by examining enriched gene ontology  terms,  and pathways  taking advantage of well annotated drosophila orthologs of these genes.
this analysis further reinforces the idea of unique caste-specific developmental trajectories in female honeybees; queen and worker larvae show different trends in terms of the pathways and processes that are more highly expressed during development.
differential activation of these pathways and processes are likely to reflect the different requirements of the adult phenotypes.
during early development queen larvae have higher expression of genes in pathways involved in cellular maintenance and growth.
genes involved in dna replication and amino acid metabolism pathways are enriched at 6 hours, and go terms include nucleotide binding, protein folding and atpase activity.
during mid larval development queen larvae have higher expression of genes involved in energy generation pathways, including the tricarboxylic acid cycle  and oxidative phosphorylation.
this switch to energy generating processes at  <dig> and 108 hours may reflect in the increased growth rate seen in queen larvae.
worker larvae are heavier than queen larvae up until 84 hours, when queens begin to gain weight faster than workers  <cit> .
in late larval development, queens have higher expression of genes with go terms associated with proteolysis.
workers have higher expression of genes involved in amino acid, xenobiotic and general metabolism from early through to mid development.
workers also have higher expression of genes involved in muscle development early in larval development.
the higher expression in workers of genes involved in muscle development was also noted by barchuk et al.
<cit> .
barchuk et al.
suggest that this may represent worker larvae preparing for their adult life as foragers.
as the thoracic muscles disappear entirely during metamorphosis  <cit> , this investment must be able to be carried through to adult muscle development to be beneficial to flight muscle strength.
similar to queens, workers have higher expression of genes involved in proteolysis during late larval development; however, in workers these genes have higher expression at 108 hours in contrast to 132 hours in queens.
this temporal shift in gene expression was not a common phenomenon in our data sets .
the difference in expression of these proteolysis genes may reflect preparation for pupation in queen larvae.
in workers, at 108 hours, the expression of these genes is unlikely to be related to pupation but instead may have a role in mediating pcd in the worker ovary  <cit> .
consistent with this hypothesis we see higher expression of genes involved in apoptotic  and autophagic cell death  in workers during mid to late larval development .
alternative splicing of mrna transcripts and caste development
one set of genes with consistently higher expression in queens were those associated with the spliceosome .
this is particularly interesting as spliceosome encoding genes are differentially methylated between castes  <cit> , and a link has been proposed between gene-body methylation and control of alternative splicing in the honeybee  <cit>  and the ant  <cit> .
dna methylation is believed to be crucial to caste development  <cit>  as rna interference targeting the de-novo dna methyltransferase  <dig>  biases development towards the queen trajectory  <cit> .
honeybees, similar to other invertebrates, have relatively low overall levels of cpg methylation restricted to gene bodies  <cit> .
over evolutionary time methylated cytosines are more often deaminated to thymines than non-methylated cytosines  <cit> , leaving a signature of low cpg content in methylated gene bodies .
the cpg content of a gene  can be used to infer the ancestral germ-line methylation state of a gene  <cit> .
while cpg ratio does not indicate the current methylation state of a gene, several studies have found good correlation between low cpg values and high dna methylation in the honeybee  <cit> .
analysis of the cpg ratios for the degs implies a dynamic pattern of switching between historically methylated and historically un-methylated transcripts as development proceeds .
if current dna methylation is related to cpg measurements, as suggested by previous studies, then our data implies a complex and changing relationship between methylation and gene expression during caste development.
dna methylation of introns and exons has been associated with alternative splicing of mrna transcripts in diverse species  <cit> .
differences in dna methylation have been detected during larval development  <cit>  and knockdown of dnmt <dig> has been shown to cause aberrant splicing in honeybee larvae  <cit> .
given the proposed link between alternative splicing and methylation in the honeybee  <cit> , and our evidence implying queens have higher expression of spliceosome related genes  we used our hts dataset to investigate alternative splicing.
discovery of alternative transcripts in our datasets indicates that ~50 % of genes expressed in larvae have alternative transcripts.
for degs both queens and workers have higher expression of more genes than expected that encode multiple splice variants in larvae  .
however, both queens and workers tend to show higher expression of a single variant only rather than multiple variants of the same gene .
thus in the majority of cases a single transcript is differentially regulated, while the other transcript for these genes remain constitutively expressed.
there are, of course, exceptions to this; in queens 11 % of degs and 16 % in workers have multiple transcripts that are differentially regulated.
in all cases, bar three, these genes are regulated co-ordinately.
for three degs, splice variants are regulated independently, i.e.
one transcript of the gene is more highly expressed in queens, another in workers.
these genes encode a dnaj-type chaperone , a rna exonuclease  and a tld-domain containing protein .
while queens and workers both have higher expression of genes that have the capacity to encode multiple splice variants, usually the transcription of only a single variant is affected.
our analysis indicates that while alternative splicing is abundant in both queen-destined and worker-destined larvae, we find little evidence for alternative splicing contributing to the queen and worker developmental trajectories.
that we find specific transcripts of genes that are differentially expressed, against a background of invariant expression of other splice variants, implies that alternative splicing may play a role in caste development, but that this role is neither a simple one nor limited to a specific caste.
the increase in spliceosome gene expression may reflect an increase in general transcription, with queens expressing more genes than workers, as shown in our hts data.
alternative splicing has recently been linked to dna methylation in the honeybee  <cit>  and the ant  <cit> .
until recently most of the evidence for a link between alternative splicing and methylation in the honeybee implies that methylation promotes exon skipping  <cit> , although evidence exists for promotion of exon inclusion  <cit> .
in general, however, alternative transcripts occur more often in methylated genes as compared with unmethylated genes  <cit> .
we know that dna methylation is integral for worker development  <cit>  and our data, that alternative splicing does not correlate with a particular developmental trajectory, hints at a role for dna methylation in caste development outside the regulation of alternative splicing.
CONCLUSIONS
caste development in honeybees is a complex process involving coordinated changes in gene expression triggered by larval diet.
our study indicates that these changes occur earlier than previously thought, implying an early biasing of gene expression in each developing caste, followed by a commitment phase which may be triggered by circulating jh .
both castes are likely to have changed from the non-social ancestor of honeybees, however we propose that the gene expression profile accompanying the queen developmental trajectory, based on the evolutionary history of the genes expressed, is more similar to that seen in the ancestral non-social bee.
the gene expression profile accompanying the worker developmental trajectory represents genes, and gene networks, that are more diverged from this ancestor implying that novel genes and pathways have been co-opted in worker development into regulating this complex polyphenism.

 caste development in honeybees is a remarkable example of, and model for, developmental plasticity. in this species, two adult female phenotypes are produced by one genome in response to diet. larvae destined to become queens are fed royal jelly which is rich in carbohydrates and lipids. rj also contains major royal jelly proteins, one of which, royalactin, is vital for queen development. as a result of a higher food intake and a more nutritious food, queen bees are larger, longer-lived, have fully developed ovaries and differ in their behaviour and morphology, when compared to worker bees. transferring larvae from worker cells to queen cells during early larval development causes bees to develop some queen-like phenotypes. based on these experiments it has been proposed that differentiation between the two castes begins on the first day of larval development and is progressive. manipulating larval diet in this way can induce queen-like properties in worker bees up to 60 hours of larval development demonstrating, that at least during early larval development, there is plasticity in the developmental pathways that give rise to queen and workers. the molecular mechanisms by which rj triggers queen development are only beginning to be understood, but rj intake causes differences in juvenile hormone titre between worker and queen larvae. these jh levels regulate gene expression. it seems likely, therefore, that a complex set of molecular mechanisms link rj exposure to changes in jh titre. such mechanisms should be reflected in differences in gene expression between developing larvae destined for different castes. differences in gene expression underlying queen and worker developmental trajectories have been investigated in a number of previous studies. these studies used a variety of techniques to profile gene expression, but have focused on mid-to-late larval development as few differences in gene expression have been identified before this stage. mid to late larval development includes the period, when jh titres differ between castes, and so if is difficult to determine if these gene expression changes are due to rj exposure, or the action of jh. while a number of mechanisms that regulate gene expression during caste development have been proposed, dna methylation is thought to be critical, as it can link environmental exposure to changes in gene expression. rna interference targeting the de-novo dna methyltransferase biases development towards the queen trajectory. honeybees, similar to other invertebrates, have relatively low overall levels of cpg methylation restricted to gene bodies. methylation within the gene body appears to influence alternative splicing in the honeybee and the ant. it is possible that the role of dna methylation in caste development may be via influencing alternative splicing of transcripts. evidence from mammalian models indicates that gene body methylation may have other molecular functions, such as promoting rna polymerase pausing, modulating other epigenetic marks, or affecting altering promoter usage. it has been proposed that plasticity, such as caste development in eusocial insects and neuronal plasticity in invertebrates, is associated with dynamic changes in dna methylation and, by inference, alternative splicing. dna methylation is likely to be one of many molecular mechanisms involved in caste development, but these mechanisms as a whole will lead to changes in gene expression, detectable using microarray techniques, or alternative splice isoforms detectable using transcriptomics. we have examined gene expression throughout larval development, both before and after the wave of jh, in queen and worker castes, carrying out the widest survey of gene expression during caste development in this remarkable animal. we used this to build a more complete understanding of how a honeybee larva responds to rj, and how this biases them into an alternative developmental trajectory, culminating in a queen bee.  experiment design there are five stages of larval development in the honeybee. queen and worker samples were taken at the following time points: 6 hours, 12 hours, 36 hours, 60 hours, 84 hours, 108 hours and 132 hours. at each time point four replicate queen and worker samples were collected. each replicate sample collected at 6 hours contained larvae and samples collected at all other time points contained larvae. queen development was induced by grafting larvae into artificial queen cells before returning them to the hive. grafting, which is a standard apicultural technique, is carried out as soon as practical after the larvae hatch. larvae in these queen cells are fed rj and this triggers queen development. when collecting 6-hour larvae it was noted that rj was present in these cells. worker larvae were not grafted. gene expression data was generated from these samples and using custom two colour long-oligonucleotide microarrays. at 60 hours, two additional biological replicates, each consisting of either queen or worker larvae, were taken and transcriptome data was generated using high throughput sequencing. the 60 hours time point was selected for hts as this time point has been highlighted as the earliest difference in gene expression seen a previous study. hts sequencing also gave us the opportunity to assess gene expression in greater depth, as well as analyse the effect of caste development on alternative splicing. differential gene expression during caste development microarray analysis indicates that queens and workers have relatively equal numbers of differentially expressed genes throughout larval development. hts of rna, however, found that 83 % of degs had higher expression in queens, similar to previous results. many studies have reported finding more differentially expressed genes using hts than using microarrays. this is at least partly because hts analysis is not hindered by the acquisition bias that can affect microarray analysis. microarrays are sensitive to saturation effects and are also not able to detect transcripts that are expressed at a relatively low level due to high amounts of background fluorescence obscuring low hybridization signal to the probes. we believe the latter phenomena accounts for the relatively high number of degs detected in queens by hts compared to microarray, as these genes had a median rpkm of compared with a median rpkm of for genes that are more highly expressed in worker larvae. we believe the relatively low expression of these genes precluded them from being detected as differentially expressed by microarray analysis. we have carried out extensive validation of these data sets using rt-qpcr which indicates that both the microarray and hts experiments are producing high-quality, biologically relevant information. genes we identified as being differentially expressed during caste development were assigned drosophila orthologs using blastx. approximately 90 % of the genes with higher expression in queen larvae have orthologs in drosophila as compared with 77 % of the genes with higher expression in workers, a statistically significant difference. barchuk et al. also found a tendency for worker larvae to have higher expression of genes that lack orthologs in drosophila. these data imply that genes more highly expressed by queen larvae have been conserved across ~ million years of evolution. in contrast, genes more highly expressed in worker larvae are not as well conserved in drosophila and may be rapidly evolving. differences in gene expression between queen and worker larvae are observed as early as 6 hours post-grafting. indeed, the largest numbers of degs identified by microarray analysis are detected at this time point. while it is possible that the large number of differentially expressed genes at this time point are due to perturbation of the larvae during grafting, the appearance of genes that are differentially expressed at 6 hours and that retain higher expression in one caste during larval development implies that at least some of this differential gene expression represents an early response to rj. we also do not observe gene ontology categories associated with stress response or physiological perturbation in the genes that are differentially expressed at 6 hours, indicating that differences we observe in expression at this early time-point are likely an early response to rj. in queen larvae, several genes with increased expression early in larval development are maintained at high levels throughout larval development. these genes are mitochondrial cytochrome c, phosphoenolpyruvate carboxykinase, phytanoyl-coa dioxygenase domain-containing protein homolog and glycine n-methyltransferase-like. worker larvae, in contrast, have higher levels of the apoptosis regulator bcl- throughout development. these genes have functions that are difficult to interpret in terms of caste development, but are consistent markers of either worker or queen development. given that queen development can be triggered in larvae by grafting until larval day three, well after the 6 hour time point by which we see substantial differences in gene expression, these early differences in gene expression must not commit the larvae to either a queen or worker fate. we therefore propose these early changes in gene expression, in response to rj, may bias development toward one caste or the other, but that later gene expression can reverse that bias. these early differences in gene expression may be attributable to the larval diets which differ, with the queen diet initially containing more sugar and less protein than worker jelly though it is not clear if these differences occur at or before 6 hours after grafting. these dietary differences have also been linked to morphological differences between the castes that are apparent as early as 2 days into larval development. the idea of an early, biased, phase of caste development is consistent with early larval grafting experiments which showed that worker larvae shifted into queen cells at ~12 hours of larval development result in bees that are queen-like, but retain some worker characteristics. a two-stage process of caste development has been previously proposed based on grafting, nutrition and physiological data, but our data provides the first molecular evidence to support an early biased phase of gene expression and development followed by a commitment phase. differential expression of hormone biosynthesis genes one key event during larval development is the wave of jh synthesis in queen larvae, during l and l stages. juvenile hormone is implicated in caste development in honeybees. our data indicates that several genes involved in jh and ecdysteroid biosynthesis are differentially expressed during larval development between castes at, or around this jh wave. queen destined larvae, have higher expression of eight cyp genes, many with roles in hormone biosynthesis. cyp315a acts in ecdysteroid biosynthesis, and has higher expression at six hours of larval development. cyp15a which has higher expression in queens throughout larval development but peaking at 84 hours, acts as a jh epoxidase, catalysing the final step of jh biosynthesis. queens also have higher expression of jh methyltransferase, involved in the final steps of jh synthesis, with differential expression peaking at 84 hours after grafting. in contrast, during mid larval development, workers have higher expression of jh esterase, which degrades jh in honeybees. seven genes induced in drosophila by jh titre are dynamically regulated, with three having higher expression in workers early in development and five having higher expression in queens during mid to late development, temporally correlating with the peak of jh titre during development. the trends in hormone-related gene expression in queen larvae, peaking as they do during mid-late larval development, are consistent with these genes either being involved in the synthesis of jh associated with the wave between the third and fourth larval instars, or with responding to the jh wave. at this point we also observe the appearance of differential expression of a number of genes in both queen and worker larvae that remain stably differentially expressed throughout the remainder of larval development. this change to a committed developmental trajectory, however, does not alter or reprogram the expression of genes we find differentially expressed at early stages. our data is consistent with an initial bias towards one caste or the other induced quite early by rj consumption, followed by a consolidation of that fate and final specification through circulating hormones. differential expression of cytochrome p genes cytochrome p genes which act in hormone and pheromone biosynthesis and the detoxification of foreign compounds show differential regulation between castes. insect cyp genes are divided into four major clans based on sequence similarity. worker larvae have higher expression of cyp genes at time points throughout larval development, the majority belonging to clan three, in particular sub-families six and nine. there is evidence for rapid expansion of these two sub-families in the honeybee through recent gene duplication events. genes in cyp clan three are loosely classified as environmental response genes. the differential expression of cyp clan three genes is most marked at 60 hours of larval development, where worker larvae have higher expression of genes of this clan. this period coincides with a change in the worker larval diet which results in the introduction of pollen grains, whereas rj only contains trace levels of pollen. the introduction of pollen exposes worker larvae to a number of foreign compounds. one example of this is the addition of quercetin, a flavinoid compound found in honey and pollen. four cyp genes are known to metabolise quercetin in the honeybee and all four have higher expression in workers during mid-larval development. some cyp genes are involved in pheromone production and possibly caste specific synthesis of mandibular acids. cyp genes have also been linked to termite caste development. functional classes of genes differentially regulated between queen and worker larvae two other families of genes were differentially expressed during caste development. queens have higher expression of five genes involved in the response to reactive oxygen species, especially during mid larval development. these genes include catalase, several thioredoxins and a peroxiredoxin. the generation of ros by the mitochondria has been proposed to be a key mechanism in the aging process. this may reflect an increased respiration rate previously observed in queen larvae. the differential regulation of these genes may influence expression of vitellogenin, important in protecting honeybees from oxidative stress. in contrast workers have higher expression of three glutathione-s-transferase genes. glutathione-s-transferases are important for detoxification of both endogenous and exogenous compounds. programmed cell death has a role in reducing ovariole numbers in larval worker ovaries during fourth and fifth larval instars. consistent with this we see higher expression of eleven genes involved in programmed cell death in worker larvae. these belong to two major programmed cell death pathways - apoptosis and autophagic cell death. hexamerin70b and caste development hexamerin genes were differentially regulated between castes at several larval time-points. rt-qpcr, which is generally considered to be the gold-standard method to measure gene expression, was used to confirm that, at 60 hours of larval development, workers have between and fold higher expression of the four hexamerin genes when compared to queens. hexamerins are a family of insect amino acid storage proteins that have evolved from the copper-containing hemocyanins of ancestral aquatic insects and crustaceans. a role for the hexamerins has been suggested in caste development in termites, where two hexamerin genes have been shown to facilitate the jh dependent worker to soldier differentiation. several previous studies have found differences in expression of hexamerin genes between queens and workers implying a possible role for these proteins in caste development. hexamerin 70b is expressed in the fat body and gonads during larval development and its expression, which is regulated by juvenile hormone, peaks prior to the expression of other hexamerins during larval development. thus far no biological function has been attributed to hex 70b; it is found in the cytoplasm, nucleus, and in punctate dots in the cytoplasm of cells within the fat body, implying that this protein has multiple biological roles, as reported for hex 70a. hexamerin70b knock-down of the four hexamerins assayed by rt-qpcr hex 70b showed the smallest fold change between queen and worker larvae. this comparatively small fold-change provides us with the best opportunity, amongst the hexamerins, to knock-down expression of this gene to queen larvae-like levels, even if rnai treatment is inefficient and knock-down is weak. given the technical challenge of rnai in honeybee larvae, and the high expression levels of all the hexamerin genes, we targeted hex 70b because moderate knockdown of this gene was more likely to produce queen-like expression levels, and give a clear, measurable, phenotype. rnai knockdown of hex 70b resulted in significant changes in the proportion of queen-like individuals produced when compared to control injected larvae. emerging individuals were classified as queen-like based on the morphology of several key structures. while hex 70b rnai produced a statistically significant change in phenotype, we did not see a significant difference in the expression of hex 70b between target and control larvae as determined by rt-qpcr, although we did observe a trend towards lower hex 70b expression in the hex 70b rnai injected group relative to controls. the lack of statistical significance likely reflects the small sample sizes used for rt-qpcr, and large variation in expression of this gene in injected individuals. samples were taken – 18 hours post-injection for this analysis and it is possible that sampling too early, or too late, may explain why the difference in expression between the target and control groups was not significant. knockdown of hex 70b biases individuals towards a queen fate, consistent with the high expression of these genes in worker larvae. we propose that hex 70b has a key role in caste development at a time when we propose a commitment to worker or queen fate is set. this implies that early biases in gene expression between worker and queen larvae can be overcome to change the developmental trajectory of a larva. the possibility that hexamerins are acting to regulate jh activity raises the possibility that we are disrupting the final commitment phase of caste development by reducing hex 70b expression, rather than changing the bias in commitment. gene ontology and pathway analysis degs involved in caste development were categorized by proposed gene function by examining enriched gene ontology terms, and pathways taking advantage of well annotated drosophila orthologs of these genes. this analysis further reinforces the idea of unique caste-specific developmental trajectories in female honeybees; queen and worker larvae show different trends in terms of the pathways and processes that are more highly expressed during development. differential activation of these pathways and processes are likely to reflect the different requirements of the adult phenotypes. during early development queen larvae have higher expression of genes in pathways involved in cellular maintenance and growth. genes involved in dna replication and amino acid metabolism pathways are enriched at 6 hours, and go terms include nucleotide binding, protein folding and atpase activity. during mid larval development queen larvae have higher expression of genes involved in energy generation pathways, including the tricarboxylic acid cycle and oxidative phosphorylation. this switch to energy generating processes at and 108 hours may reflect in the increased growth rate seen in queen larvae. worker larvae are heavier than queen larvae up until 84 hours, when queens begin to gain weight faster than workers. in late larval development, queens have higher expression of genes with go terms associated with proteolysis. workers have higher expression of genes involved in amino acid, xenobiotic and general metabolism from early through to mid development. workers also have higher expression of genes involved in muscle development early in larval development. the higher expression in workers of genes involved in muscle development was also noted by barchuk et al.. barchuk et al. suggest that this may represent worker larvae preparing for their adult life as foragers. as the thoracic muscles disappear entirely during metamorphosis, this investment must be able to be carried through to adult muscle development to be beneficial to flight muscle strength. similar to queens, workers have higher expression of genes involved in proteolysis during late larval development; however, in workers these genes have higher expression at 108 hours in contrast to 132 hours in queens. this temporal shift in gene expression was not a common phenomenon in our data sets. the difference in expression of these proteolysis genes may reflect preparation for pupation in queen larvae. in workers, at 108 hours, the expression of these genes is unlikely to be related to pupation but instead may have a role in mediating pcd in the worker ovary. consistent with this hypothesis we see higher expression of genes involved in apoptotic and autophagic cell death in workers during mid to late larval development. alternative splicing of mrna transcripts and caste development one set of genes with consistently higher expression in queens were those associated with the spliceosome. this is particularly interesting as spliceosome encoding genes are differentially methylated between castes, and a link has been proposed between gene-body methylation and control of alternative splicing in the honeybee and the ant. dna methylation is believed to be crucial to caste development as rna interference targeting the de-novo dna methyltransferase biases development towards the queen trajectory. honeybees, similar to other invertebrates, have relatively low overall levels of cpg methylation restricted to gene bodies. over evolutionary time methylated cytosines are more often deaminated to thymines than non-methylated cytosines, leaving a signature of low cpg content in methylated gene bodies. the cpg content of a gene can be used to infer the ancestral germ-line methylation state of a gene. while cpg ratio does not indicate the current methylation state of a gene, several studies have found good correlation between low cpg values and high dna methylation in the honeybee. analysis of the cpg ratios for the degs implies a dynamic pattern of switching between historically methylated and historically un-methylated transcripts as development proceeds. if current dna methylation is related to cpg measurements, as suggested by previous studies, then our data implies a complex and changing relationship between methylation and gene expression during caste development. dna methylation of introns and exons has been associated with alternative splicing of mrna transcripts in diverse species. differences in dna methylation have been detected during larval development and knockdown of dnmt has been shown to cause aberrant splicing in honeybee larvae. given the proposed link between alternative splicing and methylation in the honeybee, and our evidence implying queens have higher expression of spliceosome related genes we used our hts dataset to investigate alternative splicing. discovery of alternative transcripts in our datasets indicates that ~50 % of genes expressed in larvae have alternative transcripts. for degs both queens and workers have higher expression of more genes than expected that encode multiple splice variants in larvae. however, both queens and workers tend to show higher expression of a single variant only rather than multiple variants of the same gene. thus in the majority of cases a single transcript is differentially regulated, while the other transcript for these genes remain constitutively expressed. there are, of course, exceptions to this; in queens 11 % of degs and 16 % in workers have multiple transcripts that are differentially regulated. in all cases, bar three, these genes are regulated co-ordinately. for three degs, splice variants are regulated independently, i.e. one transcript of the gene is more highly expressed in queens, another in workers. these genes encode a dnaj-type chaperone, a rna exonuclease and a tld-domain containing protein. while queens and workers both have higher expression of genes that have the capacity to encode multiple splice variants, usually the transcription of only a single variant is affected. our analysis indicates that while alternative splicing is abundant in both queen-destined and worker-destined larvae, we find little evidence for alternative splicing contributing to the queen and worker developmental trajectories. that we find specific transcripts of genes that are differentially expressed, against a background of invariant expression of other splice variants, implies that alternative splicing may play a role in caste development, but that this role is neither a simple one nor limited to a specific caste. the increase in spliceosome gene expression may reflect an increase in general transcription, with queens expressing more genes than workers, as shown in our hts data. alternative splicing has recently been linked to dna methylation in the honeybee and the ant. until recently most of the evidence for a link between alternative splicing and methylation in the honeybee implies that methylation promotes exon skipping, although evidence exists for promotion of exon inclusion. in general, however, alternative transcripts occur more often in methylated genes as compared with unmethylated genes. we know that dna methylation is integral for worker development and our data, that alternative splicing does not correlate with a particular developmental trajectory, hints at a role for dna methylation in caste development outside the regulation of alternative splicing.  caste development in honeybees is a complex process involving coordinated changes in gene expression triggered by larval diet. our study indicates that these changes occur earlier than previously thought, implying an early biasing of gene expression in each developing caste, followed by a commitment phase which may be triggered by circulating jh. both castes are likely to have changed from the non-social ancestor of honeybees, however we propose that the gene expression profile accompanying the queen developmental trajectory, based on the evolutionary history of the genes expressed, is more similar to that seen in the ancestral non-social bee. the gene expression profile accompanying the worker developmental trajectory represents genes, and gene networks, that are more diverged from this ancestor implying that novel genes and pathways have been co-opted in worker development into regulating this complex polyphenism. 
our data, based on the biases in gene expression early in development together with published data, supports the idea that caste development in the honeybee consists of two phases; an initial biased phase of development, where larvae can still switch to the other caste by differential feeding, followed by commitment to a particular developmental trajectory. within our gene expression data we identified processes that potentially underlie caste differentiation. both queen and worker developmental trajectories are associated with the expression of genes that have alternative splice variants, although only a single variant of a gene tends to be differentially expressed in a given caste. to understand how these trajectories are established we have generated a comprehensive atlas of gene expression throughout larval development. some of these early changes in gene expression are maintained throughout larval development, indicating that caste-specific developmental trajectories are established much earlier than previously thought. we found substantial differences in gene expression between worker and queen-destined larvae at 6 hours after hatching. using rna interference, we were able to demonstrate that one of these genes, hexamerin 70b, has a role in caste differentiation.
 female larvae of the honeybee develop into either queens or workers depending on nutrition. this nutritional stimulus triggers different developmental trajectories, resulting in adults that differ from each other in physiology, behaviour and life span.  to understand how these trajectories are established we have generated a comprehensive atlas of gene expression throughout larval development. we found substantial differences in gene expression between worker and queen-destined larvae at 6 hours after hatching. some of these early changes in gene expression are maintained throughout larval development, indicating that caste-specific developmental trajectories are established much earlier than previously thought. within our gene expression data we identified processes that potentially underlie caste differentiation. queen-destined larvae have higher expression of genes involved in transcription, translation and protein folding early in development with a later switch to genes involved in energy generation. using rna interference, we were able to demonstrate that one of these genes, hexamerin 70b, has a role in caste differentiation. both queen and worker developmental trajectories are associated with the expression of genes that have alternative splice variants, although only a single variant of a gene tends to be differentially expressed in a given caste.  our data, based on the biases in gene expression early in development together with published data, supports the idea that caste development in the honeybee consists of two phases; an initial biased phase of development, where larvae can still switch to the other caste by differential feeding, followed by commitment to a particular developmental trajectory. caste developmentlarval developmentgene expressionapis mellifera 