|Home | About | Journals | Submit | Contact Us | Français|
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
In honeybees, differential feeding of female larvae promotes the occurrence of two different phenotypes, a queen and a worker, from identical genotypes, through incremental alterations, which affect general growth, and character state alterations that result in the presence or absence of specific structures. Although previous studies revealed a link between incremental alterations and differential expression of physiometabolic genes, the molecular changes accompanying character state alterations remain unknown.
By using cDNA microarray analyses of >6,000 Apis mellifera ESTs, we found 240 differentially expressed genes (DEGs) between developing queens and workers. Many genes recorded as up-regulated in prospective workers appear to be unique to A. mellifera, suggesting that the workers' developmental pathway involves the participation of novel genes. Workers up-regulate more developmental genes than queens, whereas queens up-regulate a greater proportion of physiometabolic genes, including genes coding for metabolic enzymes and genes whose products are known to regulate the rate of mass-transforming processes and the general growth of the organism (e.g., tor). Many DEGs are likely to be involved in processes favoring the development of caste-biased structures, like brain, legs and ovaries, as well as genes that code for cytoskeleton constituents. Treatment of developing worker larvae with juvenile hormone (JH) revealed 52 JH responsive genes, specifically during the critical period of caste development. Using Gibbs sampling and Expectation Maximization algorithms, we discovered eight overrepresented cis-elements from four gene groups. Graph theory and complex networks concepts were adopted to attain powerful graphical representations of the interrelation between cis-elements and genes and objectively quantify the degree of relationship between these entities.
We suggest that clusters of functionally related DEGs are co-regulated during caste development in honeybees. This network of interactions is activated by nutrition-driven stimuli in early larval stages. Our data are consistent with the hypothesis that JH is a key component of the developmental determination of queen-like characters. Finally, we propose a conceptual model of caste differentiation in A. mellifera based on gene-regulatory networks.
Phenotypic variation among individuals of the same species triggered by environmental action is an intriguing biological phenomenon that can be found in quite striking manifestations in members of different insect orders . In highly eusocial bees (Hymenoptera) one or a few females (queens) specialize in reproductive tasks, whereas a large number of quasi-sterile individuals (workers) engage in colony maintaining activities [2,3]. This polyphenism is generally determined by discrete switches during postembryonic development, and commences with the differential feeding of female larvae . The nutritional stimuli trigger an endocrine response that is manifested by an elevated juvenile hormone (JH) titer in queen larvae when compared to workers [for review see ]. The queen-inducing properties of JH were first demonstrated by Wirtz and Beetsma , who topically applied JH on fourth and early fifth instar worker larvae [for similar results in stingless bees see [7,8]]. However, the molecular mechanisms underlying this phenomenon are not yet understood. In particular, we are largely ignorant of how nutritional factors affect the endocrine system and alter JH synthesis rates of queens and workers, and how these changes drive caste-specific developmental pathways during metamorphosis.
In Apis mellifera, a model system for caste development and division of labor in social Hymenoptera, young larvae of both castes are fed with royal jelly, a secretion produced by glands in the head of adult workers. Whereas these nurse bees feed copious amounts of royal jelly to queen larvae until they enter metamorphosis, they switch the diet for late instar worker larvae from pure royal jelly to a mixture of glandular secretions with honey and pollen (worker jelly). In addition, prospective queen larvae receive 10 times more food than worker larvae . As a consequence of this differential feeding regime the two types of larvae follow two very different developmental trajectories, in spite of having exactly the same genetic background. Conceptually, this process of caste differentiation involves two kinds of alterations in the original developmental pattern (or ground plan present in ancestral solitary bees): one type, which we can call incremental alterations, affects the general growth of the body or specific organs, especially the ovaries. The other type can be considered as character state alterations that result in the presence or absence of entire specific structures, such as the pollen-collecting apparatus on the hind legs, wax glands, etc. Both types of alterations can be envisaged as JH threshold responses controlling the expression of genes involved in the development of specific organs and in specifying the general body plan (Figure (Figure11).
The first large-scale study on the molecular biology of caste differentiation was done in A. mellifera by Severson et al. . These authors demonstrated by in vitro translation analyses that queens and workers differ in their mRNA profiles during larval and prepupal stages. Later studies by Corona et al.  and Evans and Wheeler [12,13] found that most of the differentially expressed genes between prospective queens and workers were related to metabolic processes, and specifically, that queens up-regulate metabolic enzymes. Conversely, workers were shown to up-regulate a member of the cytochrome P450 family, hexamerin 2, dihydrodiol dehydrogenase and a fatty-acid binding protein. In addition, these studies revealed that several regulatory genes such as the mitochondrial translation initiation factor (AmIF-2mt), a member of the Ets family of transcription factors with a DNA binding domain, were also up-regulated in worker larvae .
In a recent study, Cristino et al.  examined the up-stream regulatory elements associated with all transcripts previously found to be differentially expressed in worker and queen larvae. They confirmed that the majority of the annotated differentially expressed genes (DEGs) are related to metabolic processes, with an interesting dichotomy for enzymes with hydrolase and oxidoreductase activities, which were found to be up-regulated in workers and queens, respectively. Genes up-regulated in workers were also shown to share more common (or conserved) overrepresented cis-elements when compared to genes up-regulated in queens.
While the aforementioned studies support the notion that incremental alterations are associated with the differential expression of physiometabolic genes, the nature of genetic mechanisms underlying the development of worker distinctiveness or character state alterations remains to be understood.
Here we used cDNA microarrays to monitor differential gene expression in honeybee queen and worker larvae and to identify cis-acting elements associated with these two developmental trajectories. Graph theory and complex networks concepts were adopted to attain an objective visual representation of the connectivity between motifs and genes in both castes. We identify groups of genes responsible for the development of queen and worker singularities and describe how their expression is co-regulated during critical stages of larval development. We also discuss the role of morphogenetic hormones in the developmental process of queen-like character and finally, propose a model of caste differentiation in A. mellifera.
We performed cDNA microarray analyses on a set of more than 6,000 unique genes throughout larval development of honeybees and identified a total of 240 genes as differentially expressed between queens and workers in the analyzed stages, namely L3, L4 and L5S2 (see Additional file 1; Methods and Gene Expression Omnibus database (NCBI), accession number GSE6452). Only genes that met stringent statistical criteria (α < 0.05; B > 0; see Methods) were selected as primary candidates.
Out of the 240 differentially expressed genes (DEGs), the majority (167) were up-regulated in fourth instar larvae (L4), 105 were up-regulated in prospective queens and 62 in worker-destined larvae [see Additional file 1]. Third instar larvae (L3) showed 37 DEGs and the fifth instar spinning stage larvae (L5S2) contributed with 36 DEGs [see Additional file 1]. This indicates that major changes in gene expression take place during the period of nutritional switch, in an environment with relatively high levels of JH for both queens and workers. Interestingly, out of 37 DEGs in L3 larvae the overwhelming majority (34) were found in workers, possibly reflecting the higher growth rate of young worker larvae . Additional 52 genes were found to be up-regulated in JH-treated larvae. These ESTs most likely represent JH-responsive genes [see Additional file 1] because they reflect transcriptional changes occurring between 1 h and 24 hrs after hormone application.
Our functional classification of the honeybee DEGs based on Drosophila melanogaster Gene Ontology (GO) annotation (only reciprocal orthologs were considered) reveals that 56% of the queen up-regulated genes and 69% of the worker up-regulated genes do not have known counterparts in Drosophila. This predominance of novel (Drosophila unrelated) genes up-regulated in worker larvae has not been observed by Cristino et al.  who analyzed a small sample of 51 genes empirically selected as differentially expressed in queen and worker larvae [12,13,11,15].
Our finding that many of the up-regulated unique genes appear to be involved in the development of worker morphological and behavioral characteristics is compatible with the complexity of workers' body plan relatively to other insects such as Drosophila. When compared to honeybee queens, whose main purpose in life is to lay thousands of eggs, worker honeybees perform a myriad of activities inside and outside the colony and can even engage in reproductive tasks when released from the queen's inhibitory influence.
Based on GO terms for Biological Processes, we found that in all examined developmental stages workers up-regulate more developmental genes than queens (Figure (Figure2A).2A). The proportion of developmental genes is always around 50% in workers, whereas in queens this proportion is always very low, with a strong bias towards physiometabolic genes. For example, in L3 larvae nine developmental genes were found to be up-regulated in workers and none in queens and in L4 larvae there were eight such genes in workers and four in queens. Considering all developmental stages, workers up-regulated 17 genes classified as participating in developmental processes, whereas only five genes were up-regulated in queens. Interestingly, all these genes are classified as participating in developmental processes related to morphogenetic differentiation of specific organs, like pollen-collecting and reproductive apparati, nervous system, wax and Nasanov's glands, etc.
In agreement with previous work done in Apis mellifera [11-14] and in Bombus terrestris  we show that most of the known DEGs are related to physiometabolic processes (57 up-regulated in queen and 29 up-regulated in worker larvae). Among these, L4 queen larvae up-regulate five genes related to the metabolism of nitrogenous compounds (GB12123, GB13298, GB10789, GB10196, GB18599; see Additional file 1), while none of these genes appear to be up-regulated in workers. This is consistent with the fact that royal jelly is richer in nitrogen compounds (amino acids and nucleotides) than worker jelly [17,18]. On the other hand, in L3 stage there is only one nitrogen metabolism-associated DEG that is up-regulated in both castes, thus suggesting the existence of a maximization of feeding differences between developing queens and workers during L3–L4 stages. Furthermore, specifically during the L4 stage, queens up-regulate more genes associated with cellular localization, protein binding, nucleotide binding, nucleic acid binding, hydrolase and oxidoreductase activities than workers (Figure (Figure2B).2B). Proteins encoded by these genes are expected to participate in the physiometabolic processes leading to the differential growth of the queen's body.
Taken together, the relative proportion of differentially expressed physiometabolic genes is not unexpected, since most genes expressed during the honeybee life cycle are classified as belonging to three main categories; metabolic, cell growth and/or maintenance processes [see ]. This phenomenon is additionally aggravated in queen larvae, whose development is shifted towards a general growth.
Within the physiometabolic category there are some DEGs encoding metabolic enzymes and also genes whose products are known to regulate the rate of mass-transforming processes and the general growth of the organism (Table (Table1).1). Conversely, many genes that are well-characterized in D. melanogaster and in other model organisms may underlie processes leading to the development of caste-biased structures. For example, genes participating in neurogenesis, leg development, apoptosis, and genes coding for components of the cellular matrix (Table (Table1).1). The protein sequences of 81 DEGs clustered by functional groups were searched against a protein domains database (Pfam), as an additional support for the putative biological roles assigned by homology to the Drosophila's counterparts (Table (Table1).1). The first three processes mentioned above are the basis for the respective morphological differences favoring the worker caste, thus defining the adult skills early in development: learning and memory, pollen and propolis collection, and a reduced reproductive capacity. The up-regulation in worker larvae of components of the cytoskeleton may reflect an early production of muscle elements, fundamental for the adult flight activities.
With regard to hormonal control of larval development, the experiment with JH treatment of developing worker larvae allowed us to identify 52 'JH-responsive' genes involved either directly or indirectly in JH signaling, specifically during the critical period of caste development. Ten DEGs that were identified in both JH and caste assays were clustered as "JH responsive + caste" genes (Table (Table1).1). Six genes were found in both JH treated larvae and L4 queens, and four genes were held in common between JH controls and L4 workers.
Queens up-regulate more physiometabolic genes than workers. Moreover, a vast majority of these genes code for products with protein-binding activities like chaperonins or chaperon binding proteins and ribosomal and proteasome related proteins with endopeptidase activities (Table (Table11).
A fundamental key controller of physiometabolic processes, also up-regulated in L4 queens, is a member of the insulin-signaling pathway tor. This conserved gene regulates translation initiation in response to nutrients from yeast to mammals. In the insect fat body, Tor acts as a sensor of amino acid levels in the hemolymph . When there are sufficient amino acids available, the activation of the Tor system triggers the expression of ALS protein (mammalian acid labile subunit ortholog), a member of a systemic communicational pathway, which signals to other larval tissues, in a PI3K-dependent manner, the nutritionally favorable situation of the organism . In corpora allata, the activation of Tor pathway could trigger the synthesis of juvenile hormone, one of the central regulators of insect development. Other responsive organs are the ovaries, where cell proliferation has been demonstrated to be responsive to the insect's nutritional environment .
Acting together with Tor in determining differential growth between queens and workers are CREG and CRC (Table (Table1).1). CREG (cellular repressor of E1A-stimulated genes), that is up-regulated in L4 workers, is a secreted glycoprotein. In humans, it antagonizes cellular transformation by E1A and ras . CRC encodes a Drosophila homolog of vertebrate ATF4, a member of the CREB/ATF family of basic-leucine zipper (bZIP) transcription factors. This protein (CRC) participates in a conserved mechanism of sensing amino acid deprivation and stress induction . Thus, the up-regulation of tor in prospective queen larvae and of creg and crc in prospective worker larvae may constitute a dual system of growth determination in response to differential feeding in honeybees.
Learning and memory-related skills that honeybee workers use for navigation, foraging, nestmate recognition, and other activities are believed to be associated with the prominent regions of insect brains, the mushroom bodies . As expected, the ratios typically used to evaluate the relative size of specific brain areas versus the body size support the notion that workers have bigger and more developed mushroom bodies than queens [26,2]. We found that five genes encoding proteins that participate in neural system development in Drosophila and in vertebrates (dac, atx2, shot, ephR and fax) are up-regulated in developing worker larvae (Table (Table1).1). Since the nervous system in queens and workers begins to differentiate during post-embryonic stages  these five genes are candidates for molecular determinants of the observed morphological differences in developing brains.
The making of an insect leg is a complex process that is even more complicated in worker honeybees because of the presence of unique structures for pollen and propolis collection . We found three genes to be up-regulated in L3 workers that have been shown to act as regulators of leg development (gug, dac and crc) and one gene regulating bristle morphology (atx2; Table Table1)1) in Drosophila . This finding is consistent with the importance of bristles morphology and density for the pollen-collecting apparatus . The temporal expression of these genes during leg imaginal disc development in the critical period of caste differentiation suggests that they have been recruited, together with other (unknown) genes, to regulate the alternative leg structures in A. mellifera castes.
Programmed cell death (PCD) is a process that, in concert with cell proliferation, modulates the development of specific organs during metazoan ontogenesis. In insects, this process is commonplace during larval nervous system development and during metamorphosis, where major systems reorganization occurs. In the honeybee, PCD in the worker ovary reduces the number of ovarioles during metamorphosis from 150–200 primordia to less than 10. It has been shown that this phenomenon is prevented from occurring in queen larvae by high titres of JH  that may inhibit the ecdysteroid-triggered PCD, as suggested for Drosophila . As shown in Table Table1,1, several genes preferentially expressed in worker larvae are associated with PCD. Among these are the cathepsin D gene (cathD) , and a gene coding for a lysozyme that participates in autophagic cell death of salivary glands in Drosophila, CG11159 . We also found atx2 (Ataxin-2) gene up-regulated in L3 workers. The transgenic over-expression of atx2 in Drosophila results in female sterility, possibly by impairing the adhesion between oocytes and follicular cells . Another example of genes up-regulated in worker larvae is traf1 (TNF-receptor-associated factor 1). Its product participates in both autophagic cell death and induction of apoptosis in Drosophila. Several genes up-regulated in L4 queens are similar to Drosophila and human genes coding for proteins with anti-apoptotic activity. One of them is lethal(2) tumorous imaginal discs (l(2)tid; 34). Another one is Trap1, an unfolded protein binding protein and a member of the heat-shock family of proteins that may play an important role in the suppression of apoptosis caused by formation of ROS . In addition, tor, which is also up-regulated in L4 queens, is a negative regulator of autophagic cell death . Thus, when activated by factors acting downstream of the insulin receptor, Tor could provide a link between nutrient availability and the regulation of tissue disintegration, such as in the ovarian tissue of honeybees. Taken together, all these genes are likely to participate in the regulatory processes underlying the differential reproductive capacity of honeybee queens and workers.
Several DEGs that are up-regulated in worker larvae code for cytoskeleton constituents (Table (Table1).1). Among these are genes encoding myosin II heavy chain (mhc), actin, (act5C), troponin T, and the genes upheld, shot, lva and ank2. In honeybees, muscle development continues until the early fifth instar allowing for the replacement of larval muscles by imaginal muscles during metamorphosis . One explanation for this increased expression of several cytoskeletal genes associated with muscle development may be the differential flight ability of adult honeybee workers. The differential expression of muscle related genes has also been observed in a closely related stingless bee, Melipona quadrifasciata . Finally, L3 workers differentially express the usp (ultraspiracle) gene that codes for the ecdysone receptor partner and a strong candidate for a JH receptor in insects , including honeybees .
Inspired by the widely known effect of JH on caste differentiation in honeybees  we examined gene expression in worker larvae following the application of exogenous JH-III. We found four genes (ccp84Ad and three unknown ones) to be up-regulated in L4 larvae 1 h after JH treatment against three up-regulated genes (tsp5D and two unknown ones) in control L4 workers. Ccp84Ad that is highly expressed in L4 queens is annotated as a structural component of larval cuticle . Since L4 queen larvae have high titres of endogenous JH , this rapid and positive regulation of ccp84Ad by JH suggests that this is an early responsive gene in the expression cascade promoted by JH. Tsp5D, a tetraspanin protein gene, is up-regulated in control L4 workers. A member of the tetraspanin family of proteins resides in cell surfaces of growth cones facilitating synapse formation during Drosophila neurogenesis . Thus, since queen larvae have higher titres of JH, the repression by exogenous JH of the expression of a gene (tsp5D) that participates in cone growing, makes this protein a candidate player in worker neurogenesis.
Among the 45 DEGs identified as up-regulated at 24 h after JH treatment, 28 were found to be induced and 17 repressed by this hormone. Interestingly, the majority of JH-responsive genes with known orthologs in Drosophila are physiometabolic genes (Figure (Figure3).3). The exceptions are unzipped (zip), which participates in axonogenesis processes , cathepsin L (Cp1), implicated in cell death processes  and embryonic aldolase . Out of the 45 JH responsive genes, 26 have known Drosophila orthologs.
Previous studies have shown that JH treatment accelerates behavioral maturation in young bees  and plays a role as an organizer of the mushroom bodies . Recently, a microarray analysis was used to understand the effects of a JH analog (Methoprene) on brain gene expression profiles during behavioral maturation of honeybee workers . Methoprene induces significant forager-like changes in gene expression even in workers with no foraging experience suggesting that the increase in JH titres may be related to expressional changes occurring during the natural transition from hive to forager behavior [47,46]. Interestingly, 17 out of 52 JH responsive DEGs identified in our work (JH 1 h and JH 24 h) in worker L4 larvae overlap with DEGs in the brain of adult workers treated with Methoprene in the study reported by Whitfield et al. . The effects of JH-like compounds caused the same shift in gene expression for 10 DEGs in both assays. Of particular interest are two DEGs considered as markers for behavior, namely Hsc70Cb (GB10836; a forager marker) and RfaBp (GB11059; a nurse marker) [48,46]. In our study, Hsc70Cb and RfaBp were up-regulated in queen larvae and in worker larvae, respectively. Furthermore, these two DEGs change their expression profiles in a hormone-dependent manner as observed in worker L4 larvae (our work) and young adult worker bees . RfaBp transcription was down-regulated by JH treatment in two very different stages of bee development (caste determination in larvae and 'hive-bee-to-forager transition' in adults). On the other hand, Hsc70Cb was up-regulated by JH treatment in worker L4 larvae as observed in queen L4 larvae (high titres of JH) but was down-regulated by Methoprene treatment in young adult bees . The contrasting effects of JH on Hsc70Cb expression during different life stages are interesting examples of how hormones can exert their regulatory effects on genes expression in a context-dependent manner. Until recently, only a small number of insect genes had been described as being either directly or indirectly regulated by JH [39,49]. Consequently, this study and the reported by Whitfield  represent a significant expansion of this functional category.
Finally, and as shown in Figure 2(B), JH treatment induces a queen-like expressional pattern. Physiometabolic and Localization genes are up-regulated in JH-treated workers, whereas Developmental genes are up-regulated in control workers. These results support the notion that JH is a potent activator of physiometabolic processes and exerts its role by inducing incremental alterations.
Most of the genes identified as differentially expressed between the two castes have been assigned to functional groups according to the literature and functional databases (GO and Pfam). We hypothesized that the groups of DEGs sharing common expression patterns might be used to infer putative clusters of co-regulated genes by means of computational analyses of DNA-sequence motifs [50-52].
A motif discovery pipeline was run on 19 sets of functional groups (Table (Table1)1) and on two sets of top10 DEGs observed in the L4 stage (Table (Table2).2). Nine sets of genes show significant differences in the distribution of motif scores when compared to control (random) sets [see Additional file 2]. The datasets neurogenesis (2), leg development (3), apoptosis/other proteins (4.2), cytoskeleton/cytoskeleton proteins (5.1), cytoskeleton/transcription associated proteins (5.2), cytoskeleton/other proteins (5.3), JH responsive/protein binding (6.3), hormone+caste (7) and top10-WL4 (9) show evidences for the occurrence of more conserved motifs than expected by chance, but only 8 motifs were considered as overrepresented (Church <= 1e-10 and ROC-AUC >= 0.7) in 4 functional groups [Table [Table3;3; ;11 motif in (4.2), 3 motifs in (6.3), 1 motif in (7) and 3 motifs in (9)]. No similarity was found by aligning these eight motifs with D. melanogaster binding site sequences described in the TRANSFAC database .
The upstream control regions (UCRs) of worker top10 genes have significantly more overrepresented DNA-sequence motifs than queen top10 genes (Table (Table3;3; three motifs in workers and none in queens). This result is in agreement with Cristino et al.  who reported that the number of overrepresented motifs discovered in worker genes was much higher than in queen genes. Interestingly, at least two of the discovered motifs were found in UCRs of 13 genes (six in queens and seven in workers, Figure Figure3)3) that were not part of the dataset used in the motif discovery pipeline. Two queen DEGs showed putative binding sites for CF1-USP (GB12001, GB13929; 54). In workers, a group of four DEGs has binding sites for CF1-USP and three DEGs for EcR-USP (55; Figure Figure33).
The main interactions among the queen DEGs are around the motifs discovered in JH responsive datasets (Table (Table1;1; Figure Figure3).3). Interestingly, Hsc70Cb gene has all four hormone-related motifs (M6-3-1, -2, -3 and M7; Figure Figure3)3) in its UCR. On the other hand, worker DEGs are largely organized around the motifs found in the top10-WL4 dataset (WL4-1, -2, -3; Table Table3,3, Figure Figure3).3). The localization of three up-regulated genes, GB16648 (usp), GB19338 (crc) and GB11059 (RfaBp, Retinoic and Fatty acid Binding Protein), in the worker's network (Figure (Figure3)3) is particularly informative. The motifs for the binding of EcR-USP and the best motif of the top10-WL4 dataset (WL4-1; Table Table3)3) were found in the UCR of usp itself, indicating that this gene can be regulated by a protein that is important for transcriptional control in workers and also by its own protein when heterodimerized with EcR, as suggested by Barchuk et al. . In crc gene, up-regulated in workers, we found three motifs, CF1-USP and two top10-WL4 motifs (WL4-1 and WL4-2; Table Table33 and Figure Figure3),3), suggesting that crc can be regulated by USP and other regulatory proteins important for worker pattern development. RfaBp gene was also found to be up-regulated in workers by Evans and Wheeler [12,13] and we show that it is up-regulated in worker larvae and responds negatively (it is down-regulated) to JH treatment. In adult workers, RfaBp expression is repressed by a JH analog  and has been considered a nurse marker gene . Two motifs were found in the RfaBp UCR, CF1-USP and the best motif of the "Hormone+caste" dataset (M7, Table Table33 and Figure Figure3),3), suggesting that RfaBp can be regulated by USP and also by another regulatory protein involved in JH response. Taken together, the available data indicate that these caste determining genes (RfaBp, crc and usp) seem to play an additional important role in regulating major phenotypic changes in honeybee development.
After searching for matches of at least 60% similarity to the discovered motifs and the two canonical patterns, CF1-USP and EcR-USP, for the UCRs of 183 DEGs (105 for queens, 78 for workers), we designed putative transcriptional networks based on the occurrence of overrepresented motifs in the UCR of DEGs between A. mellifera castes (Figure (Figure3).3). The first type of network considered here, illustrated in Figure 3(A) is a bipartite network involving both motifs and genes. The second type of network, illustrated in Figure 3(B) depicts only genes and their interrelationship. Figure 3(C) shows that most motifs, except for "hormone" and "apoptosis" motifs, are plotted above the main diagonal line (dashed line) and regulate more genes in workers.
The similar clustering coefficients obtained for workers (cc = 0.37 ± 0.23) and queens (cc = 0.36 ± 0.25) indicates similar dense connectivity among the immediate neighbors of each gene. Nevertheless, the degree measure shows that the worker network is substantially denser (d = 62.21 ± 28) showing more interconnections than that obtained for queens (d = 31.23 ± 15.67), which has also a more modularized system (Figure (Figure3B).3B). In agreement with Cristino et al. , this result indicates that workers genes are more strongly interrelated (Figure (Figure3).3). Moreover, the obtained networks show that worker DEGs share significantly more conserved motifs than expected by chance.
All the results presented here suggest the existence of groups of co-regulated genes and highlight potential key genes (tor, usp, crc, RfaBp) which might determine developmental processes leading to the formation of caste-specific morphogenetic fields.
In an attempt to consolidate our results with published data, we propose the following model for caste differentiation in A. mellifera (Figure (Figure4).4). The type of food eaten by the larva must first be recognised by the receptor system in the gut epithelial cells followed by complex signalling via the stomatogastric nervous system [; Maleszka et al., in preparation] that sends the information to the brain and the retrocerebral endocrine system. Corpora allata (CA) activity and the behavior of target tissues may be even under the upstream regulation of insulin/IGF molecules eventually secreted by the neurosecretory cells of the larval brain , as suggested for Drosophila [see ]. The other sensing organ, the fat body (FB), receives the information directly from the hemolymph. Thus, high level inputs in these organs in prospective queens result in the activation of the insulin/IGF pathway and the Tor system, which in turn increases the levels of JH synthesis in the CA, and may trigger the ALS-mediated systemic communicational pathway in FB. "Worker jelly", on the other hand, affects insulin/IGF pathway in a less pronounced manner, and may not be able to increase the levels of JH above a specific threshold.
As a result of the activated signaling pathways, the high titres of JH in prospective queens regulate the expression of physiometabolic genes that together with the available nutrients from royal jelly determine the general body growth pattern. In this model, the up-regulation of tor in prospective queen larvae can be seen as a determinant of the observed differential growth rates. On the other hand, low levels of JH combined with limited nutrient availability in prospective workers lead to the development of smaller adults. In this case, the up-regulation of creg and crc (negative growth regulators) in workers may constitute the second determinant of a dual system of growth regulation in response to differential feeding in honeybees.
Furthermore, reaching a JH threshold in prospective queens (Figure (Figure1)1) not only permits general body growth but also acts by negatively regulating the development of some organismal systems that are characteristic of adult workers and are also present in the original developmental pattern. High levels of JH may, for instance, inhibit the development of worker-specific leg structures , as well as prevent cell death in the ovaries, leading to a higher number of ovarioles in adult queens. Moreover, JH is likely to control the extended brain development specifically that of MB, resulting in bigger worker brains when compared to queen brains [26,2]. Thus, JH has contrasting effects on growth and development, at least on certain structures and might be regulating trade-off processes during caste differentiation in A. mellifera.
Supporting our general model, the genes in the queen putative regulatory network are mainly associated with motifs discovered in the hormone responsive dataset (Figure (Figure3).3). In contrast, the genes (usp, crc, RfaBp and actin) in the worker network are mainly connected via worker- and apoptosis-biased motifs.
Several steps of the proposed pathway leading to caste development have better experimental support than others. Without doubt, the proximate mechanism(s) linking Tor activity and JH synthesis (a general biological issue) are still unclear. We also have to investigate whether an ALS-like mechanism of systemic communication of the organism's nutritional condition is functional in honeybees (the A. mellifera ortholog of als gene is GB20133). We further need to deepen our understanding of a number of molecular pathways that are critical for the establishment of caste identities, including the development of the nervous system, ovaries and legs. Finally, silencing key genes like usp, crc and RfaBp in developing workers and tor in queens would be particularly advantageous in gaining novel insights into the behavior of the proposed regulatory networks.
Two types of spotted microarrays were used in this study. A custom-made small array contained 768 ESTs implicated in processes believed to be important for caste determination; 723 of these ESTs were manually selected from the ORESTES project  and 45 were arbitrarily chosen from the existing ESTs encoding transcription factors and microRNAs. The small array was produced in the Adelaide Microarray Facility. The second array representing brains of mixed age workers was constructed by Robinson and colleagues (bEST_9000; University of Illinois at Urbana-Champaign, USA; ). Together, these slides account for more than 60% of the honeybee genes.
The differential feeding between the prospective queens and workers begins after larval stage 3, when the nurse bees continue to feed queen larvae with large quantities of royal jelly, whereas they include honey and pollen in the worker larvae's diet, reducing the amount of royal jelly. This switch in feeding is somehow linked to a differential JH synthesis determining higher titres of this hormone in queens than in workers's hemolymph . To identify genes that are responsive to these hormonal and nutritional changes we tested RNAs from third instar (L3), fourth instar (L4) and fifth instar spinning stage larvae (L5S2), characterized by high (L3–L4) and low levels of JH (L5S2).
As JH is known to govern the induction of queen development in highly eusocial bees, we also tested RNAs coming from worker larvae treated with JH-III in L4 (samples obtained 1 h and 24 h after hormone treatment). These hybridizations could give us information about the genes responsive to exogenous JH, whose profiles could be compared to those obtained in normal L3–L4 queens and thus, highlight those genes regulated by JH and responsible for the development of prospective queens.
Honeybee larvae were collected from A. mellifera colonies (Africanized hybrids) maintained at the Experimental Apiary of the University of São Paulo at Ribeirão Preto, Brazil. Larvae of the same and determined age were obtained as in Barchuk et al. . The developmental stages were classified according to the criteria proposed by Michelette and Soares .
To test the effects of exogenous JH on worker larvae gene expression profile, fourth instar larvae received a topical application of 10 μg JH III (C16-juvenile hormone, Fluka Biochemika, 59992; diluted in acetone to a stock solution of 2.5 μg/μL). The amount of applied hormone was based on our previous experiments in which we examined the induction of vg and usp expression by JH during post-embryonic development in honeybees [61,39] and on a pioneer work of Goewie and Beetsma  of "artificial" queens development by JH application. JH-applied larvae (n = 80–100) in brood frames were left in an incubator (34°C and 80% relative humidity) for 1 h previous putting them back into the colony where they were maintained for 24 h. After the appropriate time, larvae (JH treated and control) were included in TRIzol reagent (Invitrogen) and frozen at -80°C until RNA extraction.
Microarray experiments are described according to the MIAME specifications  and the resulted data have been deposited in the Gene Expression Omnibus (GEO, at NCBI database) under the accession number GSE6452 . The 768_EST microarrays were prepared from PCR amplified cDNA clones using either M13 or specific primers. The amplicons were washed and EtOH precipitated before preparing the final samples for spotting. cDNA fixation on glass supports was done by the Adelaide Microarray Facility . Prior to hybridization each slide was incubated in pre-hybridization solution (10 mg/mL BSA, 25% formamide, 5× SSC and 0.1% SDS) for 60 min at 42°C then rinsed in Milli-Q water and dried by centrifugation at 750 rpm for 5 min.
We used 30 arrays (768_EST and bEST_9000) for both kinds of experiments: developmental gene expression and JH-responsive genes evaluation. For internal controls we used heterologs or reference genes like phosphoglycerate kinase 1 and 2-microglobulin (cattle), rubisco small chain 1 and chlorophyll ab binding protein (soy; see ) and ribosomal protein S8 (honeybee). Each slide included two (768_EST) to four (bEST_9000) replicates. Pairs of RNA samples (labeled with Cy3 or Cy5 fluorophores) from prospective queens and workers from each developmental stage (L3, L4 and L5S2) and JH-treated and control workers were hybridized to the same slide. Dye swaps were done for each comparison and more than two slides were used to evaluate some developmental stages.
Total RNA was extracted following Invitrogen's protocol combined with column purification (RNeasy Mini Kit, QIAGEN, Cat. 74104). RNA quality was determined by electrophoresis in agarose gels as in Kucharski and Maleszka . cDNA synthesis, RNA amplification and probe synthesis (from 1 μg of amplified RNA) were done with and according to Low RNA Input Fluorescent Linear Amplification Kit (Agilent, Cat. 5184-3523). Probes purification was done using MinElute PCR Purification Kit (QIAGEN, Cat. 28004).
For hybridization, probes (in 90 μL 2× SSC) were pre-heated at 50°C during 3 min and placed on microarrays under lifter-slip cover glasses (22 × 60, 31.25 μL). cDNA microarrays were placed in single slide hybridization chambers which were incubated in a water bath for hybridization for 17h, at 50°C. Washing procedure included the following steps: 2× SSC and 0.1% SDS; 2× SSC; 0.1× SSC and bidistilled water; 3 min per wash, all at room temperature. Slides were dried by centrifugation at 750 rpm for 5 min previous scanning.
Slides were scanned using Affymetrix 428 Array Scanner and Jaguar 2.0 software; gain 40–68 dB, 10 micron resolution, Cy3 with Green Laser (532 nm), Filter FM570-10 and Cy5 with Red Laser (635 nm), Filter FM665-12. For images quantitation we used ScanAlyze 2.35  with default parameters.
Background correction was performed to avoid the correction of negative or zero intensities (offset correction) by adding to the background-corrected intensities a positive constant (= 50). It damps spurious variation in log-ratios, particularly at low intensity spots. The "print-tip loess" normalization was used to correct for within-array dye and spatial effects and single channel normalization was used to facilitate comparison between arrays . Since our slides had 4 (768Br) and 2 (bEST_9000) within array replicate spots, the normalization procedure also considered the variation among their values . After normalization we determined a log2 ratio (Queen's sample intensity/Worker's sample intensity), for each probe on each array. The fold-change in expression and its standard error for each gene were calculated by fitting a linear regression to the normalized expression data. A Bayesian smoothing procedure was used to shrink the estimated standard errors, with which we calculated the moderated t-statistic for each gene .
Differentially expressed genes (ESTs) from 768_EST and bEST_9000 slides were assigned to GB IDs of the GLEAN3-predicted protein sequences . All DEGs were annotated according to Gene Ontology terms for Biological Process and Molecular Function . FatiGO web tool  was used to annotate Biological Process and Molecular Function terms based in D. melanogaster sequence similarity for the target genes matched by the overrepresented cis elements.
In order to detect biologically relevant motifs (cis-elements) in the upstream control regions (UCR) of the sets of JH- and caste-related genes we selected the input data set based on two different criteria: (1) those genes that have shown significant differential expression between castes and in hormone assays observed in the microarray analyses, and (2) those genes differentially expressed in castes that shared functional similarity with empirical evidence already described in the literature (Tables (Tables1,1, ,2,2, and Additional file 1). A motif discovery script was designed based on reliable strategies proposed by MacIsaac and Fraenkel  and Cristino et al. . The pipeline run separately on the two main sets of UCR sequences combining the output of three programs: AlignAce , MEME  and MDscan . Default parameters values were used in all searches, except that GC-content background in intergenic regions was set to 25% for AlignAce and a background distribution file was computed for MDscan.
A database containing 10,123 UCR sequences was generated by parsing the Official Gene Set annotation file assembly 3.0 (downloaded in GFF format from 75) to extract upstream regions starting from the terminal 5' genomic coordinate of each predicted CDS. The UCRs were arbitrarily set to a size frame of 1000 nucleotides, but were trimmed whenever another predicted ORF was detected in any of these regions .
MAP (maximum a priori log likelihood), group specificity score (Church score) , ROC-AUC (area under the curve for a receiver operator characteristic plot) metric and MNCP (mean normalized conditional probability) metric  were used to score the discovered motifs. Motifs databases were generated for each subset of genes with MAP score equal to and greater than 5.0. A non-parametric test (Kolmogorov-Smirnov two-sample test) was conducted to infer significance levels for the 19 set of discovered motifs. The distribution of the four motif's score metrics for all 19 set of (1) hormone- and caste-differentially expressed genes and in (2) functional groups up-regulated in castes [see Additional file 2] was tested against motifs discovered in control random dataset. Control UCR dataset were 100 times randomly selected from the genomic background (10,132 UCRs) equal in size to each 19 input set and control motifs were generated by running the motif discovery script in randomly selected sequences. Only motifs with very stringent scores (MAP ≥ 5; Church ≤ 1e-10; ROC-AUC ≥ 0.7) were considered in this manuscript.
Functional ecdysone response elements have been identified and it is now well established that the EcR-USP complex binds to direct or inverted (palindromic) repeats (rGkTCAATGaMCy) [78,55,79]. Another binding site pattern involving USP (CF1-USP) was already identified in D. melanogaster s15 chorion gene (GGGGTCAcs) . Both known motifs were searched in all 183 DEG.
The discovered and known motifs were represented by position weight matrices (PWM) [80,81]. Each bee motif was aligned against the D. melanogaster sequences in the TRANSFAC database (release 4.0) . Only the alignments with a threshold of 80% identity for each PWM were considered as significant matches.
The interrelation between overrepresented motifs and genes was modeled as a bipartite network based on concepts from graph theory  and complex networks . The clustering coefficient (cc) and degree (d; 83) were calculated from subsumed networks, which were obtained from the respective bipartite representation by linking all pairs of genes connected to a same motif, thus resulting in networks with only one type of node.
An Ubuntu Linux (version 6.06; 84) operating system was used to implement all scripts and pipelines designed for annotation procedures and motif discovery. The Python programming language , Biopython , and TAMO (Tools for Analysis of Motifs) packages  were used in program design. For the microarray analyses, all normalizations and fold change calculations were performed using the functions in the library Limma of the R/Bioconductor package . For the detection of conserved domains, the 84 protein sequences used to motif search were screened against the Pfam database  using the HMMER platform (current release 2.3.2; 90), with a cutoff value set at 1e-5.
ARB, ZLPS and RM designed research; ARB and RK performed research; ARB coordinated the collaboration and drafted the manuscript; AdSC performed the bioinformatics analyses and participated in drafting the manuscript. LdFC made the statistical analyses of the networks and supervised AdSC. ZLPS supervised ARB and AdSC and helped drafting the manuscript. RM supervised ARB and RK and participated in writing the manuscript. All authors read and approved the final manuscript.
We thank Luiz Roberto Aguiar for technical assistance in the apiary of Ribeirão Preto and Francis Nunes for his help with the identification of the Brazilian EST clones. We also thank Sylvain Forêt for his suggestions on the use of Limma and Inês C. Schmidt Capella for help in JH application experiments. Klaus Hartfelder's reading and commenting on a previous version of the manuscript is greatly appreciated.
This work was funded by grants from Fundação de Amparo à Pesquisa do Estado de São Paulo, Brazil, Fapesp (Proc. 99/00719-6; 2005/03926-5; 04/08928-3) to ZLPS. LdFC was supported by CNPq (308231/03-1) and FAPESP (05/00587-5). ARB was supported by postdoctoral fellowships from Fapesp (Proc. 04/03408-1 and 05/53083-4) and AdSC by doctoral fellowship from MEC-Capes (Brazil).