Insights into the dynamics of hind leg development in honey bee (Apis mellifera L.) queen and worker larvae - A morphology/differential gene expression analysis.

Phenotypic plasticity is a hallmark of the caste systems of social insects, expressed in their life history and morphological traits. These are best studied in bees. In their co-evolution with angiosperm plants, the females of corbiculate bees have acquired a specialized structure on their hind legs for collecting pollen. In the highly eusocial bees (Apini and Meliponini), this structure is however only present in workers and absent in queens. By means of histological sections and cell proliferation analysis we followed the developmental dynamics of the hind legs of queens and workers in the fourth and fifth larval instars. In parallel, we generated subtractive cDNA libraries for hind leg discs of queen and worker larvae by means of a Representational Difference Analysis (RDA). From the total of 135 unique sequences we selected 19 for RT-qPCR analysis, where six of these were confirmed as differing significantly in their expression between the two castes in the larval spinning stage. The development of complex structures such as the bees' hind legs, requires diverse patterning mechanisms and signaling modules, as indicated by the set of differentially expressed genes related with cell adhesion and signaling pathways.

Introduction

Within the framework of a relatively fixed body plan, winged insects have become the most species-rich group in terrestrial ecosystems, especially so the holometalous orders, which have made their appearance about 350 mya (Misof ). The first order to branch off as a clade in the Holometabola is the Hymenoptera (wasps, ants and bees). While most of the hymenopteran species and families have a carnivorous or parasitic life style, bees have switched to a plant-based diet, in particularly pollen and nectar. This switch, which is thought to have occurred in the early-to-mid Cretaceous (Michener, 2007), represents not only a major lifestyle novelty within this group, but also prompted a very strong co-evolutionary relationship with angiosperm plants. Not only are bees the major pollinators of angiosperm plants and, by doing so, guarantee the bio-diversity of most terrestrial ecosystems, they nowadays also provide multibillion dollar ecosystem services to agriculture worldwide. Notably, most of the high value crops (fruits, nuts, legumes) are strongly dependent on pollination by bees, especially so the corbiculate bees (Klein ; Gallai ).

In the corbiculate bees (corbiculate Apidae, sensu Michener, 2007), the female’s hind legs present a special set of structures adapted for pollen collection and transport, this being the corbicula and the pollen comb on the outside of the tibia, and rows of bristles on the enlarged basitarsus forming the pollen brush. The corbicula is a smooth area on each hind tibia that is surrounded by a fringe of stiff bristles, some of which are thought to be mechanoreceptors that allow the bee to estimate the amount of pollen carried during foraging (Ford ; Proctor ).

While the presence of a corbicula is a synapomorphic character of this clade, its expression has undergone a considerable change in the highly eusocial honey bees and the stingless bees, where it makes its appearance in the worker caste only, but not in queens. Thus, although the presence of a corbicula is an ancestral character in the female sex of corbiculate bees (Cardinal and Danforth, 2011), the secondary loss of these structures in queens apparently reverted their hind leg architecture to a more ancestral state of the Hymenoptera.

This interesting phenomenon is a result of caste determination in the premetamorphic larval stages of the highly eusocial bees, and this process is best understood in the honey bee, Apis mellifera, where caste differentiation is trigged by a switch in larval diets during the early developmental stages (for reviews see Haydak, 1970; Leimar ; Hartfelder ). Such triggers are the elevated sugar levels in royal jelly fed to queen larvae (Asencot and Lensky, 1988), and a specific protein moiety present in royal jelly, the MRJP1 monomer named Royalactin (Kamakura, 2011). While the former may affect caste development via the insulin/insulin-like signaling (IIS) pathway (Wheeler ), the latter has been shown to drive queen differentiation via the Egfr pathway (Kamakura, 2011). These upstream nutritional factors and associated signaling pathways eventually converge onto a coordinated endocrine response resulting in caste-specifically modulated juvenile hormone (JH) and ecdysteroids hemolymph titers (Rachinsky ). JH controls caste-specific cellular developmental programs in honey bee larvae, including programmed cell death in the ovaries (Schmidt-Capella and Hartfelder, 1998; Leimar ), and it also regulates the development of other tissues and organs via a modulatory action on the hemolymph ecdysteroid titer (Rachinsky and Engels, 1995).

Throughout larval development, the locally acting signaling pathways [IIS (Wheeler , 2014; Azevedo and Hartfelder, 2008), Egfr (Kamakura, 2011), TOR (Patel ) and hypoxia (Azevedo )], and globally acting morphogenetic hormones (JH and ecdysteroids) are thought to directly or indirectly affect gene expression patterns underlying the formation of the distinct queen and worker phenotypes (Evans and Wheeler, 2000; Cristino ; Barchuk ). While the above-cited transcriptome analyses were all done on RNA extracted from whole body preparations, little is known on expression profiles associated with tissue-specific development. First approaches towards elucidating the molecular underpinnings of developmental plasticity in specific tissues were done on ovaries (Hepperle and Hartfelder, 2001; Humann and Hartfelder, 2011), due to their fundamental differences associated with reproductive division of labor.

With respect to hind leg and specifically so corbicula development, a microarray screen on differentially expressed genes in hind legs of third instar worker larvae (Bomtorin ) drew attention to four genes (grunge, dachshund, cryptocephal, ataxin-2), all previously known as regulators of leg development and bristle morphology in Drosophila melanogaster. This study also investigated the localization pattern for two homeobox proteins, the honey bee Ultrabithorax (Ubx) and Abdominal A (AbdA) homologs, in the developing hind legs of honey bee pupae, showing Ubx-free spots in the tibia epidermis of worker pupae where the characteristic corbicula bristles are expected to be formed (Bomtorin ).

So as to further our understanding of the processes underlying the complex structural difference in the hind legs of queens and workers, we addressed this question through a histological analysis on the development of their hind legs throughout the critical stages of postembryonic development. Furthermore, we generated differential gene expression libraries for hind leg imaginal discs of last instar queen and worker larvae. This was done in a developmental stage when the imaginal discs just start to evaginate and form the adult compartments with their respective caste-specific structures, as seen in the histological analysis. The libraries were generated by means of a Representational Difference Analysis (RDA) protocol (Hubank and Schatz, 2000; Pastorian ) that had already been successfully adapted for studies in bees (Judice ; Colonello-Frattini and Hartfelder, 2009; Humann and Hartfelder, 2011) and also in other social insects (Weil ; Oppelt ). Though considerably more labor intensive than microarray protocols, RDA has the advantage that it is not limited to a pre-existing platform. Furthermore, it specifically selects and enriches for differentially expressed transcripts, unlike RNAseq, which gives a general transcriptome result. Due to the sequential amplification steps, RDA does, however, not give strictly quantitative readings, thus requiring secondary confirmation of differential expression, generally done by quantitative PCR. Among the total of 135 contigs and singlets assembled for the two RDA libraries generated herein we selected 19 for RT-qPCR analysis, whereby six of these were confirmed as differing significantly in their expression levels among the two castes in the larval spinning stage. Two of these genes had Gene Ontology attributes related to cell adhesion and epithelial development, while one represented a putative microRNA.

Materials and Methods

Honey bee larva

Worker larvae were directly retrieved from brood combs of colonies kept in the Experimental Apiary of the University of São Paulo, Ribeirão Preto, Brazil. Queens were reared from first instar worker larvae transferred to queen cups to be raised in queenless hives, following standard apicultural procedures. The developmental stages were identified following the criteria of Rachinsky and Michelette and Soares (1993).

Histology

For the histological analyses, larvae were collected from the fourth instar until the prepupal stage of the last, fifth instar. Imaginal discs from queens and workers larvae were dissected in Ringer solution and immersed in Bouin’s fixative (35% formaldehyde/saturated picric acid/glacial acetic acid 5:15:1). The samples were processed according to standard procedures, following inclusion in paraffin or methacrylate resin (Historesin, Leica). Briefly, tissues were dehydrated in a graded alcohol series (50–100%). For paraffin embedding they were further dehydrated in benzene, infiltrated in two baths of liquid paraffin (62 °C) for 1 h each, and a final bath under vacuum conditions for 30 min, before embedding in paraffin blocks. Embedding in Historesin was done directly after alcohol dehydration, following the manufacturer’s instructions. Sections of 5–6 μm thickness were stained with hematoxylin and eosin, and analyzed and photographed with a Zeiss microscope equipped with an AXIOCam-HR3 system (Zeiss). Images were processed with Adobe Photoshop CS software for brightness correction.

Cell proliferation assay

After dissection, imaginal discs were fixed in 4% paraformaldehyde solution (pH 7.2) for 2–3 h at room temperature. Samples were next incubated in PBS/glycine (0.1 M/0.1 M; pH 7.2) for 1 h, following three washing steps in PBS (15 min each). The tissues were then washed in PBS containing 2% Triton X-100 (PBST) for 15 min before an overnight (4 °C) permeabilization in this same solution. Subsequently, samples were treated with PBST supplemented with 1% bovine serum albumin (BSA) for 1 h at room temperature before incubation overnight at 4 °C in the mitosis marker Anti-phospho-Histone H3 (Ser10) antibody generated in rabbit (Upstate/Millipore-Merck), diluted at 10 μg/mL in PBST. After washing twice in PBS and once in PBST, they were then incubated overnight at 4 °C in a sheep anti-rabbit IgG Cy3-conjugated secondary antibody (Sigma-Aldrich) diluted 1:100 in PBST.

For immunofluorescence analysis, at least, three pairs of leg discs of each developmental stage were whole-mounted on slides in a drop of glycerol and viewed in a laser confocal microscope system (Leica TCS-SP2 and TCS-SP5). Images were captured with Leica LAS AF Lite software. Adobe Photoshop CS5 software was used for brightness adjustments. Negative controls were done with leg imaginal discs that were not incubated in primary antibody.

Representational Difference Analysis (RDA)

RNA extraction and preparation of double-stranded cDNA

A total of 25 pairs of hind leg discs were used to construct each cDNA library. The pairs of hind leg imaginal discs were dissected from early spinning-stage larvae (L5S1) and transferred to TRIzol reagent (Invitrogen) for RNA extraction following the manufacturer’s standard protocol. After treatment with RNase-free DNaseI (Fermentas), RNA quality was checked by electrophoresis in an agarose gel (1.2%) run under denaturing conditions, and RNA quantity was determined spectrophotometrically (Nanovue, GE Healthcare) at 260/280 nm. cDNA libraries were generated from 2 μg of total RNA by reverse transcription and long distance PCR using the SMARTPCR cDNA Synthesis (Clontech) protocol. The reverse transcription step was carried out using the SMART IIA and SMART CDS IIA primers (Clontech) and Superscript II (Life Technologies) enzyme. These primers contain a recognition sequence for MboI restriction enzyme. Double-stranded cDNA was then produced using the long distance primer PCR SMART IIA and Platinum Taq DNA Polymerase High Fidelity (Life Technologies) in a PTC 200 thermocycler (MJ Research) with a protocol of 94 °C for 5 s, 65 °C for 5 s and 68 °C for 8 min. Aliquots of the reaction mixture were collected at two cycle intervals between cycles 14 and 28 to find the optimal cycle number (clear bands and no signs of product overcycling). This optimal cycle number was then used in a mass-production PCR step of double-stranded cDNA.

Representational difference analysis

A Representational Difference Analysis (RDA) protocol (Pastorian ) optimized for honey bee studies (Judice ; Colonello-Frattini and Hartfelder, 2009; Humann and Hartfelder, 2011) was employed. This suppression subtractive hybridization technique is based on two subsequent hybridization steps that enrich for differentially expressed transcripts. Briefly, double-stranded cDNA (1.5 μg) was digested with MboI restriction enzyme (Fermentas), and 1 μg of the product purified by means of a GFX Illustra PCR DNA and Gel band purification (GE Healthcare) protocol was ligated to R-adaptors (Judice ), followed by PCR amplification. This created the driver populations for L5S1 queen and worker imaginal discs. For production of the respective tester populations, the R-adaptors were removed from the PCR products and these were then ligated with two sets of new adaptors (J and N) in two successive rounds of PCR amplifications. In the first round, the tester cDNA was hybridized to a 1:100 excess of driver cDNA, in the second round the tester/driver ratio was set at a higher stringency of 1:800 (Hubank and Schatz, 2000). Before each round of hybridization, the adaptor-ligated cDNAs were purified using a GFX Illustra PCR DNA and Gel band purification kit (GE Healthcare). In the permutational hybridizations, tester cDNA from L5S1 worker hind legs was hybridized with driver cDNA from L5S1 queens and vice versa.

The differential products represented in each library were then ligated into pGEM-T Easy Vector (Promega) and used to transform E. coliDH5α chemocompetent cells. A total of 192 clones grown on LB agar plates in the presence of ampicillin were picked from the two libraries and the inserts purified in a miniprep protocol. Sequencing of the miniprep products amplified by a Big Dye Terminator Cycle Sequencing Ready Reaction (Applied Biosystems) protocol using M13 forward primer was done on an ABI-PRISM 3100 (Applied Biosystems) automated gene analyzer.

Bioinformatics analyses

A sequence analysis pipeline was customized using the E-Gene platform (Durham ). Quality of the sequencing reads was evaluated by Phred before primer and vector sequences were removed. The trimmed reads were then filtered for quality and read size and removal of mitochondrial DNA and ribosomal RNA, and the resulting expressed sequence tags (ESTs) were assembled using Phrap and CAP3 (Huang and Madan, 1999) software to produce the final list of contigs and singlets for each library. These unique sequences (US) were then fed into BLAST searches (Blastn and Blastx) against GenBank sequences. Sequences that did not match to any of these were entered into a new search for transcripts from other honey bee EST projects deposited in the NCBI sequence read archive (SRA) representing transcriptome shotgun assembly (TSA) reads. The US without matches to other organisms were also queried against the microRNA stem loop sequence database miRBase (Griffiths-Jones ). And finally, all US were mapped to the latest version of the A. mellifera whole genome assembly (Amel_4.5, available in BeeBase).

Quantitative RT-PCR

A subset of 19 genes represented in the two libraries was chosen for quantitative RT-PCR (RT-qPCR) analysis. Primers for these sequences were designed with Primer3 software and Primer Blast (NCBI), with predicted products ranging from 120 to 193 bp (Table 1). For calculation of relative expression ratios, two previously validated endogenous control genes were used, one encoding a cytoplasmic actin (act) and the other representing ribosomal protein 49 (rp49, also known as rpl32) (Lourenço ).

Table 1

Primers used in quantitative RT-PCR analyses. Gene names and gene IDs are those attributed in BeeBase Official Gene Set v3.2. Unpredicted genes are named by the genomic scaffold numbers where the respective transcripts were mapped.

Gene name	Gene ID	Primer sequence	Product size
tenectin	GB53632	5′ TGCGCCCGTTTACGACCAGG 3′5′ ACGGGCAAACCCTGAGGTAGCA 3′	176 bp
headcase	GB50117	5′ GCGCTTCGGTAGTACCACGACC 3′5′ ACCAGGAACGGGGAATGCGC 3′	164 bp
bowl	GB49901	5′ CGGTGGAGGACCAGTCAGATGAGG 3′5′ TGGTCGCTGGATTGGGCCAC 3′	125 bp
NDRG3-like	GB54474	5′ GGCCGCCGCCCTAGAAATGG 3′5′ TCGACGCTCGCTATGGAAACG 3′	169 bp
Maternal gene required for meiosis	GB50425	5′ GCTCACGGTACTTGAGCGAC 3′5′ GGAGTTGATGGGTGTGGATG 3′	125 bp
Glutaredoxin-like	GB52955	5′ TGGACGAGGTCACCACGATGA 3′5′ ACCGTTGCTCCGGCCCCTGT 3′	145 bp
Grainy head	GB46725	5′ CTCTCGACGAGCTCAAATCC 3′5′ AACCAACAGTCCCACAGTCC 3′	145 bp
14-3-3ε	GB42560	5′ CCTGATAGAGCATGTCGTCT 3′5′ CCTTGCATGTCTGACGTCC 3′	149 bp
Immunoglobulin-like and fibronectin	GB44117	5′ TACCCGACCCATCCTCATTA 3′5′ GATAGGGTAGACCGGACACG 3′	150 bp
C-type lectin 1	GB49260	5′ CGGACGTGGAAAAGGGGAACGAAGG 3′5′ ACACAAGTTGAACATTCGCGCCG 3′	193 bp
SCAR	GB47014	5′ CCTCCTCCTCCTGGTGAAAC 3′5′ GGTGGGGTGGGTGTGGAGTG 3′	128 bp
gustavus	GB55797	5′ CGATACAACACGCGACACAC 3′5′ GTACGAACGAGCGAAAAACG 3′	163 bp
troponin	GB55598	5′ TCCGAGAAGGCAGAGTGGCAG 3′5′ GCGCGGTTCGGTTGTTGCAC 3′	174 bp
Group 15.19	GB50099	5′ GACGAGACCAGCGACGTTTC 3′5′ GATCTTCTTCTTCTTCGTCTTGGTC 3′	120 bp
Group13.5	unpredicted	5′ AGCCGCGATAGAGAGAGGCG 3′5′ GGGAAGAAAGGGCGAGGGGG 3′	147 bp
Group3.1	unpredicted	5′ GAGACGTCGGGGTCGTTGTCG 3′5′ CCTTCTCGAGGATGTCTCTCTCCC 3′	127 bp
Group2.14	unpredicted	5′ AGACTCGTCCCCTGTCACC 3′5′ CTCGACGACTCTTTCTGCTCG 3′	154 bp
Group1.37	unpredicted	5′ TCATAAGCCACGCGAATACC 3′5′ AGGCAGACACGTGAATACGC 3′	146 bp
Group2.11	unpredicted	5′ TTGCCGACTTCCCTTACCTAC 3′5′ ACCTTGAAAATCCCGGAGAG 3′	124 bp

Total RNA was extracted from imaginal discs of queen and worker larvae with TRIzol reagent (Invitrogen), treated with RNase-free DNase I (Fermentas) to remove any contaminating genomic DNA, and then purified with the RNeasy kit (Qiagen). Aliquots of 2 μg of RNA were reverse transcribed using SuperScript II enzyme (Life Technologies). The RT-qPCR assays were done using a Maxima SYBR Green/RoxqPCR Master Mix (Fermentas) protocol in an ABI-Prism 7500 Real-Time PCR System (Applied Biosystems). Quantification assays to determine primer efficiency were prepared from standard curves of serial cDNA dilutions (1:5 to 1:500). Dissociation curves were analyzed to verify product specificity. The amplification protocol comprised one cycle at 50 °C for 2 min and 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Three biological replicates were prepared, each composed of 10 pairs of imaginal discs, and each biological replicate was assayed as technical triplicates. The cycle threshold values for each biological replicate were used as inputs to the Relative Expression Software Tool (REST) for statistical analyses (Pfaffl ). A p ≤ 0.05 was considered as statistically significant.

Results

Imaginal disc histology and cell proliferation analysis

The legs of honey bees develop from imaginal discs. Sitting in the ventral midline of the larval body they are contiguous with the body wall epidermis, representing pockets of thickened epidermis beneath the ventral thoracic cuticle (Figure 1B,C,F). We analyzed leg disc development in queens and workers from the fourth instar larvae (L4) until the prepupal stage (L5PP) and all subsequent descriptions equally suit both castes, as there were no obvious histological differences between the two.

Figure 1

Hind leg imaginal discs development in honey bee larvae. (A–B) In the fourth larval instar (L4,) the disc occupies the peripodial sac (PS and arrow in A). Mesenchymal cells (MC) are observed inside the disc, which is covered by the larval cuticle (C). (C–E) During the feeding stage of the fifth instar (L5F) the discs grow considerably and the leg compartments start to become apparent, but the epithelium (E) remains unchanged. Muscles cells (M) are adjacent to the epithelium where they will establish future connections. (F) In the legs of early spinning stage (L5S) larva, segmentation is more apparent, but the epithelium still remains the same. A, B - L4 queen; C - L5F1 worker; D - L5F2 worker; E - L5F3 worker; F, L5S1 queen.; C, cuticle; E, epidermis; FB, fat body; MC, mesenchymal cells; M, muscle; PS, peripodial sac; Tr, trachea; TG, thoracic ganglion. Scale bars: A, C, D, E, 100 μm; B, 50 μm; F, 200 μm.

As in other holometabolous insects, the honey bee leg discs develop within a peripodial space, limited by a peripodial membrane, which is an epithelial monolayer continuous with body wall epidermis on the one hand and with the thickened disc epithelium on the other (Figure 1A–E). The disc epithelium is extensively folded within the peripodial sac (PS), giving it a multilayered appearance (Figure 1). Inside the disc pouch, small mesenchymal cells are visible (Figure 1A–D).

During the feeding stage of the last larval instar (L5F) the discs grow and the leg compartments start to appear (Figure 1C–E), but the epithelial structure remains the same as that seen in the more cup-shaped discs of fourth instar larvae. This growth occurs through cell proliferation, detected in the entire structure, including the peripodial membrane, mesenchymal cells and disc epithelium (Figure 2). During this larval stage the most distal portion of the leg disc, the future tarsal region, folds and starts to project towards the posterior portion of the respective thoracic segment. Furthermore, myoblasts that originated from mesenchymal cell precursors start to aggregate close to the epithelium, where they form the future muscle fibers (Figure 1E,F). At this point, tracheae were observed to invade the internal disc space (Figure 1D).

Figure 2

Cell proliferation detection in the hind leg imaginal discs from queens and workers by anti-phospho-histone H3 immunofluorescence (yellow). Cell proliferation was evidenced during all stages analyzed, but it appears to be more prominent during the feeding stage. (A) L5F2 queen; (B) L5F3 queen; (C) L5F2 worker; (D) negative control. Scale bars: A, B, D, 100 μm; C, 50 μm.

During the larval-pupal transition, here represented as the L5S1 phase (Figure 1F), the leg discs evaginate out of the peripodial space and expand beneath the larval cuticle where they gradually acquire the pupal leg shape. The leg segments can now clearly be distinguished, but the cuticle that covers the leg is still the larval one and, thus, does not yet present any bristles. Being characteristics of the adult legs, these will only develop later, during the pupal-adult transition.

As far as caste differences during larval development are concerned, the leg imaginal discs of queens and workers did not show any overt differences in overall structure, nor were there any noteworthy differences with respect to cell proliferation. Nonetheless, the fourth and early fifth instars are the stages when caste fate of the legs and other caste-specific structures is decided (Dedej ).

Differential gene expression - RDA library characteristics

So as to gain insight into molecular underpinnings of these caste-specific differentiation processes we performed a Representational Difference Analysis. From the 192 sequencing reads representing the queen imaginal disc RDA library, 138 were considered of good quality (Phred quality ≥ 20). After E-Gene pipeline and CAP3 analysis, these were assembled into 11 contiguous sequences (contigs) and 60 single sequences (singlets). Blastn and Blastx analyses against the GenBank non-redundant database revealed significant similarities for 62 sequences, considering a cutoff E-value of 1.0 × e−5. Among these, 25 sequences had significant similarity to Drosophila genes (Table 2), and for most of these, Gene Ontology attributes could be retrieved from Flybase. Another five sequences were similar to stem loop microRNA sequences deposited in miRBase (Table 2).

Table 2

Unique sequences assembled from expressed sequence tags of the queen imaginal discs RDA library. Listed are genes with similarity to Drosophila melanogaster genes or stem-loop micro RNAs. The genomic scaffold (Group in Amel_4.5) and predicted CDS (GB number of Official Gene Set v3.2). Gene Ontology terms for Biological Process and Molecular Function are from Flybase and given as inferred from electronic annotation (IEA), inferred from sequence or structural similarity (ISS), inferred from physical interaction, (IPI), or inferred from direct assay (IDA).

Group	GB	Blast X (CG)	miR base	Molecular function	Biological process
Group1.3	GB50376	CG33197		metal ion bindingIEA	apoptotic process, embryo development, muscle organ development, regulation of gene expression, compound eye photoreceptor cell differentiation
Group1.40	GB49260	CG43164		unknown	unknown
Group3.4	GB49076	CG8671		unknown	dsRNA transport
Group3.8	GB55797	CG2944 gustavus		protein bindingIPI	dorsal appendage formation, germ cell development, oocyte anterior/posterior axis specification
Group3.9	GB53632	CG13648 tenectin		integrin bindingIDA	embryonic hindgut morphogenesis, imaginal disc-derived wing morphogenesis, epithelial tube morphogenesis
Group3.15	GB47014	CG4636 SCAR		protein bindingIPI, actin bindingIEA	actin cytoskeleton organization, cell adhesion mediated by integrin, compound eye morphogenesis
Group4.1	GB49650	CG30497		unknown	unknown
Group4.1	GB50425	CG45477 maternal gene required for meiosis		chromatin bindingIDA, sequence-specific DNA bindingIDA, metal ion bindingIEA	female meiotic division, regulation of chromatin organization, regulation of transcription from RNA polymerase II promoter
Group5.15	GB47845	CG2302 Nicotinic Acetylcholine Receptor alpha 7E		acetylcholine receptor activityIDA	cation transport
Group6.10	GB54474	CG15669 Misexpression suppressor of KsR2		unknown	unknown
Group6.14	GB52211	CG5215 Zinc-finger protein at 72D		mRNA bindingISS, zinc ion bindingIEA	phagocytosis, mRNA splicing, via spliceosome
Group7.9	GB42560	CG31196 14-3-3ε		protein binding IPI, protein heterodimerization activityIPI	determination of adult lifespan, DNA damage checkpoint, germarium-derived oocyte fate determination, regulation of mitosis; wing disc dorsal/ventral pattern formation
Group10.23	GB48399	CG6890 Tollo		transmembrane receptor activityISS	innate immune response in mucosa, peripheral nervous system neuron development, regulation of protein glycosylation
Group10.25	GB52624	CG17077 pointed		protein bindingIPI, repressing transcription factor bindingIPI	anatomical structure development; localization; post-embryonic organ morphogenesis; anterior/posterior axis specification
Group11.16	GB47245	CG33141 sticks and stones		protein bindingIPI	actin cytoskeleton organization, compound eye morphogenesis, myoblast fusion, nephrocyte filtration, spermatid development
Group11.18	GB44976	CG14029 vrille		RNA polymerase II core promoter proximal region sequence-specific DNA binding transcription factor activity involved in negative regulation of transcriptionIDA	circadian rhythm; bristle morphogenesis; imaginal disc-derived wing hair organization
Group11.18	GB44972	CG17484 Adherens junction protein p120		cadherin bindingIEA	adherens junction maintenance; cell adhesion; compound eye morphogenesis
Group13.12	GB49901	CG10021 brother of odd with entrails limited		metal ion binding IEA; nucleic acid binding IEA, sequence-specific DNA binding RNA polymerase II transcription factor activityISS	imaginal disc-derived leg morphogenesis, lateral inhibition, embryonic foregut morphogenesis, embryonic hindgut morphogenesis
Group13.12	GB49872	CG6831 rhea		integrin bindingIPI, protein bindingIPI	apposition of dorsal and ventral imaginal disc-derived wing surfaces, border follicle cell migration, cell adhesion, muscle attachment, regulation of cell shape
Group14.14	GB52743	CG10847 encore		nucleic acid binding IEA	germarium-derived oocyte fate determination; cystoblast division; oogenesis; germarium-derived egg chamber formation
Group14.15	GB41625	CG13676		chitin bindingIEA	chitin metabolic process
Group15.19	GB50099	CG33521		zinc ion bindingIEA	unknown
Group15.19	GB50117	CG15532 headcase		unknown	RNA interference; axon extension involved in development; neurogenesis
Group16.2	GB54355	CG5654 ypsilon schachtel		DNA bindingIEA	mRNA splicing, via spliceosome, oogenesis, regulation of transcription, DNA-template
Group16.8	GB45968	CG4145 Collagen type IV		extracellular matrix structural constituentIEA	cardiac muscle cell development, oviduct morphogenesis
Group1.20			osa-MIR5149
Group1.37			sbi-MIR396c
Group3.3			mmu-mir-466f-3
Group8.21			sbi-MIR396c
Group16.2			smo-MIR1082a

The same analysis performed on the sequencing reads for the worker RDA library resulted in 131 quality ESTs, which assembled into 19 contigs and 45 singlets, with 56 sequences presenting significant matches to GenBank sequences. Among these, 15 had best matches with Drosophila genes (Table 3) and six sequences to stem loop microRNAs (Table 3). As in the queen library, some of the worker RDA library sequences also did not present similarities with nr database sequences.

Table 3

Unique sequences assembled from expressed sequence tags of the worker imaginal discs RDA library with similarity to Drosophila sequences or stem-loops micro RNAs. Listed are the genomic scaffold (Group in Amel_4.5) and predicted CDS (GB number of Official Gene Set 3). Gene Ontology terms for Biological Process and Molecular Function are from Flybase and given as inferred from electronic annotation (IEA), inferred from sequence or structural similarity (ISS), inferred from physical interaction, (IPI), or inferred from direct assay (IDA).

Group	GB	Blast X (CG)	miR base	Molecular function	Biological process
Group1.40	GB49260	CG43164		unknown	unknown
Group2.19	GB55598	CG7178 wings up A		unknown	cardiac muscle tissue development, nervous system development, nuclear division
Group3.8	GB55828	CG9160 NADH dehydrogenase (ubiquinone) acyl carrier protein		ACP phosphopantetheine attachment site binding involved in fatty acid biosynthetic processISA	fatty acid biosynthetic process, mitochondrial electron transport, NADH to ubiquinone
Group4.13	GB44117	CG12676 echinoid		protein bindingIPI	dorsal appendage formation, imaginal disc-derived wing morphogenesis, sensory organ development, morphogenesis of an epithelium, negative regulation of epidermal growth factor receptor signaling pathway
Group5.9	GB46725	CG42311 grainy head		sequence-specific DNA binding transcription factor activityIDA, protein homodimerization activityIPI	chitin-based cuticle development, epithelial cell morphogenesis, nervous system development, open tracheal system development
Group5.11	GB42054	CG5670 Na pump α subunit		sodium:potassium-exchanging ATPase activityIMP, ATP bindingIEA	adult locomotory behavior, determination of adult lifespan, neuromuscular process, cation transport, sensory perception of sound
Group6.23	GB52955	CG31559		electron carrier activityIEA, protein disulfide oxidoreductase activityIEA	cell redox homeostasis
Group8.21	GB51653	CG17927 Myosin heavy chain		actin-dependent ATPase activityIDA, protein homodimerization activityIDA, structural constituent of muscleIMP	adult somatic muscle development, locomotion, muscle organ development, myofibril assembly
Group10.24	GB54239	CG11352 jim		metal ion bindingIEA, sequence-specific DNA binding transcription factor activityISS	dendrite morphogenesis, regulation of chromatin silencing
Group10.26	GB50974	CG30084 Z band alternatively spliced PDZ-motif protein 52		muscle alpha-actinin bindingIPI, protein bindingIPI	muscle structure development, myofibril assembly
Group11.18	GB44972	CG17484 Adherens junction protein p120		cadherin bindingIEA	adherens junction maintenance, cell adhesion, compound eye morphogenesis
Group12.13	GB40670	CG10573 knockout		unknown	motor neuron axon guidance
Group12.5	GB40240	CG2184 Myosin light chain 2		calcium ion bindingIEA	flight, muscle system process, myofibril assembly
Group14.13	GB54192	CG4651 Ribosomal protein L13		structural constituent of ribosomeIDA	centrosome duplication, mitotic spindle elongation
Group15.19	GB50117	CG15532 headcase		unknown	axon extension, imaginal disc-derived wing morphogenesis, RNA interference, neurogenesis
Group1.1			aly-MIR408
Group1.35			oan-mir-153-1:
Group2.19			mtr-MIR2608
Group8.16			mtr-MIR2676f:
Group15.19			aly-MIR169j
GroupUn98			mmu-mir-6240:	GroupUn98

For validation of origin, Blastn and Blastx analyses were done on all library sequences against the nucleotide and the consensual protein databases of predicted A. mellifera genes (Amel Official Gene Set v3.2). Sequences without BLAST matches to predicted honey bee genes were mapped to the assembled honey bee genome (Amel 4.5). Thus, all sequences could be grouped into six categories: (i) sequences with similarity in the nr database, represented herein as similar to Drosophila, since this is the best annotated insect genome; (ii) sequences corresponding to predicted honey bee genes, including those from the bee Transcriptome Shotgun Assembly (TSA) database; (iii) sequences with similarity to stem loop microRNAs; (iv) sequences mapping in intronic regions of predicted coding sequences in the honey bee genome; (v) sequences mapping close to predicted coding sequences in the honey bee genome, thus indicating possible location in an untranslated region (5′ or 3′ UTR); and (vi) sequences without information in public databases, including predicted honey bee genes, but mapping to the honey bee genome. Figure 3 shows the relative representation of sequences for the two libraries within these six categories.

Figure 3

Characteristics of the queen and worker RDA libraries prepared from queen and worker larval hind leg discs. The assembled sequences were grouped into (i) genes similar (E ≤ e−5)to Drosophila melanogaster, (ii) genes similar to predicted or exclusive transcripts for A. mellifera and other bees represented in the Transcriptome Shotgun Assembly (TSA) database, (iii) genes similar to stem-loop microRNAs (miRBase), (iv) genes with intronic location, (v) genes with a putative UTR location, and (vi) sequences representing unpredicted genes but mapping to the honey bee genome sequence.

Quantitative RT-PCR analyses

We designed primers for 19 sequences that were selected based on the following criteria: representation in the libraries (high number of ESTs per unique sequence) and similarity with Drosophila genes related with appendage development and/or body size. The genes selected for this analysis are listed in Table 1, and the REST analysis results for their differential expression in hind legs of worker vs. queen spinning phase larvae (L5S1) are shown in Figure 4. In this analysis, a statistically significant difference in expression was confirmed for six genes, four of which overexpressed in workers and two in queens. These were: Group2.14, Group1.37, gustavus, and C-type lectin 1, and tenectin and IGFn3-5, respectively.

Figure 4

Relative expression levels measured by RT-qPCR of genes selected from the RDA libraries. Each whisker plot represents worker gene expression relative to queens, the latter set at equal to 1. The plot was generated by REST analysis (Pfaffl ). Black asterisks represent genes validated as overexpressed in queens and white ones those overexpressed in workers (random test, REST, p ≤ 0.05).

Discussion

In the grub-like larvae of holometabolous insects, the adult legs, like other appendages, develop from ectodermal cells set aside during embryonic development to form larval imaginal discs. Under the endocrine milieu prevailing during larval development, these epidermal imaginal discs proliferate but do not secrete cuticle (Nijhout, 1994; Truman and Riddiford, 2007). Their high regenerative potential was fundamental for the transplantation studies that led to the concept of homeosis and discovery of homeotic genes in Drosophila melanogaster. Leg development is, thus, a complex process that involves several patterning genes coordinating proper development along the three axes. Most of these genes were functionally identified in Drosophila melanogaster (Campbell and Tomlinson, 1995; Diaz-Benjumea ; Couso and Bishop, 1998; Kubota ; Brook, 2010; Giorgianni and Mann, 2011).

In contrast, most studies on imaginal disc development in bees are rather old (Myser, 1954) or not very detailed. The most complete study is, to our knowledge, a description of imaginal disc development in the stingless bee Scaptotrigona postica (Cruz-Landim, 2009), but this one is more concerned with differences in developmental timing between the leg and other imaginal discs, than with caste-related processes. Information on molecular underpinnings of caste-specific leg development is also lacking, as differential gene expression in caste development has mostly been assessed through whole body analyses, which may mask mechanisms underlying the development of individual caste-specific structures. The current study thus aimed at analyzing hind leg development of honey bee larvae, both from a morphological, as well as a gene expression perspective.

Morphology of the developing hind leg of the honey bee

During leg disc development in A. mellifera, the imaginal disc epithelium grows within a peripodial sac, as in Drosophila(Ursprung, 1972), and although consisting of a monolayer it has a multilayered appearance due to extensive folding inside the peripodial space (Figure 1). The leg discs rapidly increase in cell number during the larval feeding stage (here shown by labeling of the mitosis marker phospho-histone H3, Figure 2). Subsequently, during the spinning stage, the leg compartments start to differentiate, while cell proliferation continues. Mesenchymal cells in the luminal side of leg discs were described in Drosophila (Ursprung, 1972) and S. postica (Cruz-Landim, 2009) as precursors of neural cells and muscles. Progressive compartimentalization of the leg imaginal discs became evident during the fifth larval instar, culminating with the evagination of the legs and other imaginal disc-derived appendages during the larval-pupal molt (Figure 1F).

As expected, in the larval stages studied here, the imaginal discs did not present obvious morphological differences related to caste-specific structures, such as a corbicula or bristle patterning. These structures become discernible in brown-eyed pupae only (Bomtorin ), i.e. at the beginning of the pharate-adult stage, and are concomitant with general bristle and trichome development in the bee’s epidermis. These events are synchronized by the hemolymph ecdysteroid titer which shapes the cuticle structure and pigmentation of pharate-adult honey bees (Elias-Neto ; Soares ). While these processes of pupal/adult cuticle structuring are now reasonably well understood, both in terms of morphology and gene expression, comprehension of the developmental processes underlying appendage formation and differentiation in the larval stages of the honey bee is still lacking. The current study is the first one to investigate the morphological dynamics of leg development during the larval/pupal transition, a stage that is critical for the caste-specific characteristics of these structures (Dedej ).

Differential gene expression in honey bee leg development

Although the leg imaginal discs of queens and workers do not show any notable differences in their morphologies in the larval stages, the caste-related endocrine millieu (Rembold, 1987; Hartfelder and Engels, 1998) and the gene expression environment differ considerably (Barchuk ; Bomtorin ). For our investigation on molecular events underlying the differentiation of caste-specific hind leg structure we, thus, chose a developmental stage right at the onset of metamorphosis. At this stage, growth of the discs ended and disc eversion begins to take place.

Overall, the queen and worker larval hind leg libraries were similar in terms of the total number of ESTs and the numbers of contigs and singlets. Furthermore, the hind leg libraries comprised only a relatively small number of unpredicted Apis mellifera genes, i.e. without similarity to non-redundant sequences in public databases, which is in contrast with RDA honey bee libraries generated earlier for other tissues (male accessory glands, Colonello-Frattini and Hartfelder, 2009; and larval ovaries, Humann and Hartfelder, 2011). These previous libraries had, however, been analyzed against an earlier gene prediction set (Amel Official Gene Set 2.0), which comprised an underestimated number of genes (10,157). This was subsequently corrected through Transcriptome Shotgun Assembly data and is now updated in a revised honey bee genome version which comprises over 15,000 protein coding genes (Elsik ). Some of the still unpredicted sequences probably represent regulatory long non-coding RNAs, which represent a considerable part of the mammalian genomes. For the honey bee, only few such regulatory RNAs are known, including two that have recently been revealed in an RDA study on ovary development in honey bee larvae (Humann ).

The RT-qPCR analyses revealed six sequences among the 19 selected for validation as significantly differentially expressed (Figure 4). The two sequences that were more expressed in queens were tenectin and immunoglobulin-like and fibronectin type III domain containing 5 (IGFn3-5). IGF n3-5, known as echinoid in D. melanogaster, is predicted as GB49260 in the honey bee genome. The protein is related with cell adhesion, playing a role in lateral inhibition of proneural clusters and the Notch signaling pathway. It has been described as an Egfr antagonist in Drosophila, acting during bristle patterning (Escudero ) and compound eye development (Fetting ). It also acts upon actin/myosin structures during epithelial development and modulates cell adhesion (Bogdan and Klämbt, 2001; Tepass and Harris, 2006; Laplante and Nilson, 2011).

tenectin is predicted in the honey bee Official gene Set 3.2 as GB53632, located in the genomic scaffold 3.9. It encodes an integrin binding protein related to developmental processes and is regulated by ecdysone during wing development. Transcripts for this gene were recently identified in Drosophila leg discs, eye and brain (Fraichard ).

Thus, interestingly, the genes that had their differential expression confirmed in queens were both cell adhesion molecules and related to certain signaling pathways. The possible responsiveness to ecdysone and involvement in the Egfr pathway furthermore indicate that these proteins not only maintain the adhesive properties of the disc epithelium, but may also coordinate or integrate upstream or downstream signals in leg differentiation. The role of ecdysone in caste development is well established (Rachinsky ; Hartfelder ; Mello ), and the Egfr pathway has recently been revealed as a major player in caste fate determination (Kamakura, 2011).

Among the four sequences that were confirmed by RT-qPCR as differentially expressed in workers, only two had a prediction in the public database. The first one encodes a C-type lectin I, predicted as GB49260 in the honey bee genome. It encodes a calcium-dependent lectin that belongs to a large family of endocytic receptors, proteoglycans, colectins and selectins. This lectin type presents high evolutionary divergence from the ones found in vertebrates, which suggests that it may actually have functions other than those typically assigned to lectins (Dodd and Drickamer, 2001).

GB55797 encodes the honey bee SOCS box 1 protein containing an SPRY domain, similar to the D. melanogaster gene gustavus(CG2944). According to its Gene Ontology attributes, the gustavus protein is related to dorsal appendage formation, oocyte anterior/posterior axis specification and intracellular signal transduction. SOCS-box proteins act as E3 ubiquitin ligases.

Also overexpressed in workers were the transcripts Group2.14 and Group1.37. These two sequences did not have a prediction in the A. melliferaOfficial Gene Set (v3.2), but all could precisely be mapped to the genome sequence. While no further information exists for the Group2.14 sequence, Group1.37 has similarity to stem-loop microRNAs. This microRNA, however, does not belong to any of the known miR families listed in miRBase for the honey bee. Nonetheless, not all honey bee microRNAs have been identified so far. Alternatively, stem loops are also frequently found as structural elements of long non-coding RNAs, as seen in the two honey bee long non-coding RNAs identified in ovary development (Humann ).

Genes previously associated with honey bee leg development (Barchuk, ) and the patterning genes analyzed by Bomtorin did not appear in the RDA libraries. The developmental stages analyzed by Barchuk , as well as the analysis method used, were however different from those of the current study. The tissue-specific quantitative analyses on the leg development candidate genes abdominal-A, cryptocephal and ultrabithorax on the other hand did not reveal expression differences in queens and workers during the last instar spinning stage (Bomtorin ), which was the developmental stage in focus in the present study.

Among the other honey bee genes with reads in the RDA libraries, several are similar to Drosophila melanogaster genes and with Gene Ontology (GO) attributes, related to muscle development, cell adhesion, wing imaginal disc development and epithelial morphogenesis. The gene encoding the protein SCAR (suppressor of cyclic AMP receptor) was represented in the queen library. SCAR is an actin nucleation factor that acts on myoblasts migration and fusion during D. melanogaster embryogenesis (Gildor ). The Drosophila gene bowl (brother of odd with entrails limited, also known as bowel) regulates joint formation, with mutants presenting fusion of legs segments. It belongs to the odd-skipped family and acts downstream of the Notch pathway during leg joint morphogenesis (Hao ; Kojima, 2004).

Another interesting finding was the identification of 14-3-3ε and ypsilon-schachtel transcripts in the queen library, as both had previously been identified using the RDA methodology on honey bee male accessory glands (Colonello-Frattini and Hartfelder, 2009). Furthermore, a proteomic approach had also identified a 14-3-3ε protein as overexpressed in 72 and 120 hours workers larvae (Li ). 14-3-3 proteins are a family of conserved regulatory molecules that participate in a wide range of cellular processes, such as cell proliferation, cancer progression, apoptosis, cell cycle control, actin cytoskeleton regulation. Since several isoforms are present in all eukaryotes, some of them could have overlapping functions (Freeman and Morrison, 2011; Gardino and Yaffe, 2011; Aitken, 2011). Furthermore, another 14-3-3ε function appears to be that of a modulator of signal integration by FoxO (Nielsen ), downstream in the insulin/TOR signaling pathway. Especially the latter is of importance in honey bee caste development (Patel ). The fact that 14-3-3 proteins represent a large family of relatively conserved proteins may have been a reason why we could not find statistical differences in the RT-qPCR assays.

ypsilon-schachtel (yps) is a Y-box protein component of an ovarian ribonucleoprotein complex that acts as translational regulator of oskar mRNA in Drosophila (Mansfield ). Again, the quantitative real time PCR did not validate it as differentially expressed in queens and workers. Nevertheless, the recurrent appearance of 14-3-3ε and yps in honey bee RDA libraries generated for different tissues and developmental processes indicates that theses genes may be redundantly involved in developmental pathways common to both castes and sexes.

Another gene related to epithelial development revealed in the worker RDA library was the transcription factor grainy head.In Drosophila development it plays a role as an ecdysone-response gene in cuticle differentiation (Gangishetti ). In bees it has previously been denoted as differentially expressed in caste development of the stingless bee Melipona quadrifasciata (Judice ) and it was also identified as a major factor promoting thelytokous reproduction in the Cape honey bee Apis mellifera capensis (Lattorff ).

In conclusion, our histological analysis of honey bee hind leg development, represents the first detailed description of general growth and differentiation processes that occur during the larval developmental stages critical for determining the caste fate of these structures (Dedej ). As comprehensively shown in the fruit fly, development of the body appendages depends on several patterning mechanisms and signaling modules. In honey bee development, an additional layer of complexity is added to these processes due to the caste-specificity in anatomical specializations, such as the corbicula and pollen brush on the hind legs of workers. In the differential gene expression analysis by means of RDA library sequencing we could see this reflected in the detection of a series of genes predicted as related to cell adhesion and signaling pathways, as well as a set of reads mapping to genome sequences for which no gene function has been predicted so far. Obviously, especially this set of novel genes should deserve attention in future studies.