ORIGINAL RESEARCH ARTICLE published: 22 December 2014 doi: 10.3389/fgene.2014.00445 Developmental regulation of ecdysone receptor (EcR) and EcR-controlled gene expression during pharate-adult development of honeybees (Apis mellifera) Tathyana R. P. Mello1, Aline C. Aleixo1, Daniel G. Pinheiro2, Francis M. F. Nunes3, Márcia M. G. Bitondi4, Klaus Hartfelder5, Angel R. Barchuk6* and Zilá L. P. Simões4 1 Departamento de Genética, Faculdade de Medicina de Ribeirão Preto, Universidade de São Paulo, São Paulo, Brazil 2 Faculdade de Ciências Agrárias e Veterinárias, Universidade Estadual Paulista, São Paulo, Brazil 3 Departamento de Genética e Evolução, Centro de Ciências Biológicas e da Saúde, Universidade Federal de São Carlos, São Carlos, Brazil 4 Departamento de Biologia, Faculdade de Filosofia, Ciências e Letras de Ribeirão Preto, Universidade de São Paulo, São Paulo, Brazil 5 Departamento de Biologia Celular, Molecular e de Bioagentes Patogênicos, Faculdade de Medicina de Ribeirão Preto, Universidade de São Paulo, São Paulo, Brazil 6 Laboratório de Biologia Animal Integrativa, Departamento de Biologia Celular, Tecidual e do Desenvolvimento, Instituto de Ciências Biomédicas, Universidade Federal de Alfenas, Alfenas, Brazil Edited by: Greg J. Hunt, Purdue University, USA Reviewed by: Xavier Belles, Institute of Evolutionary Biology (CSIC-UPF), Spain Joanna Kelley, Stanford University, USA *Correspondence: Angel R. Barchuk, Departamento de Biologia Celular, Tecidual e do Desenvolvimento, Instituto de Ciências Biomédicas, Universidade Federal de Alfenas, Rua Gabriel Monteiro Silva, 700, CEP 37130-000 Campus Sede, Alfenas, Brazil e-mail: barchuk@unifal-mg.edu.br; arbarchuk@yahoo.com Major developmental transitions in multicellular organisms are driven by steroid hormones. In insects, these, together with juvenile hormone (JH), control development, metamorphosis, reproduction and aging, and are also suggested to play an important role in caste differentiation of social insects. Here, we aimed to determine how EcR transcription and ecdysteroid titers are related during honeybee postembryonic development and what may actually be the role of EcR in caste development of this social insect. In addition, we expected that knocking-down EcR gene expression would give us information on the participation of the respective protein in regulating downstream targets of EcR. We found that in Apis mellifera females, EcR-A is the predominantly expressed variant in postembryonic development, while EcR-B transcript levels are higher in embryos, indicating an early developmental switch in EcR function. During larval and pupal stages, EcR-B expression levels are very low, while EcR-A transcripts are more variable and abundant in workers compared to queens. Strikingly, these transcript levels are opposite to the ecdysteroid titer profile. 20-hydroxyecdysone (20E) application experiments revealed that low 20E levels induce EcR expression during development, whereas high ecdysteroid titers seem to be repressive. By means of RNAi-mediated knockdown (KD) of both EcR transcript variants we detected the differential expression of 234 poly-A+ transcripts encoding genes such as CYPs, MRJPs and certain hormone response genes (Kr-h1 and ftz-f1). EcR-KD also promoted the differential expression of 70 miRNAs, including highly conserved ones (e.g., miR-133 and miR-375), as well honeybee-specific ones (e.g., miR-3745 and miR-3761). Our results put in evidence a broad spectrum of EcR-controlled gene expression during postembryonic development of honeybees, revealing new facets of EcR biology in this social insect. Keywords: honey bee, adult development, 20E, ecdysteroid, juvenile hormone, JH, RNAi, miRNA INTRODUCTION Most multicellular organisms go through developmental tran- sitions that enable them to cope with environmental changes and/or broaden their niche possibilities. Such transitions are generally timed and synchronized by morphogenetic hor- mones in a broad range of species, including insects, amphib- ians, metamorphic fish, tunicates, echinoderms, and plants. In insects, developmental transitions, such as larval and meta- morphic molts, are driven by steroid hormones (ecdysteroids) acting in conjunction with juvenile hormone (JH). These hormones also control reproduction and aging (Flatt et al., 2005; Gáliková et al., 2011), and, in social insects, play important roles in caste polyphenism (Hartfelder and Emlen, 2012). The steroid hormone ecdysone is produced by the prothoracic glands. After secretion, it is transported via the hemolymph to its target organs. Due to its lipophilic nature it passes directly into the cytoplasm of target and/or modification center cells (Iga and Kataoka, 2012; Ono, 2014), where it can be modified to 20-hydroxyecdysone (20E) by a 20-monooxygenase encoded by the shade gene, a member of the cytochrome P450 family (CYP314a1) known as Halloween (Petryk et al., 2003). The mode of action of JH, which is a sesquiterpenoid morpho- genetic molecule, has only recently become clear (for review see www.frontiersin.org December 2014 | Volume 5 | Article 445 | 1 http://www.frontiersin.org/Genetics/editorialboard http://www.frontiersin.org/Genetics/editorialboard http://www.frontiersin.org/Genetics/editorialboard http://www.frontiersin.org/Genetics/about http://www.frontiersin.org/Genetics http://www.frontiersin.org/journal/10.3389/fgene.2014.00445/abstract http://community.frontiersin.org/people/u/186480 http://community.frontiersin.org/people/u/70980 http://community.frontiersin.org/people/u/139472 http://community.frontiersin.org/people/u/199147 mailto:barchuk@unifal-mg.edu.br mailto:arbarchuk@yahoo.com http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Bellés and Santos, 2014), both in terms of its receptor and downstream cascade, as well as its molecular interaction with ecdysteroids. Produced by the corpora allata in the retrocerebral complex, JH relies on binding proteins for its transport in the hemolymph to target cells. There, it first binds to its intracellular receptor, the Methoprene- tolerant (Met) protein, which then forms a complex with Taiman (Charles et al., 2011). This dimeric hormone- receptor complex then regulates the expression of target genes. Knowledge on the mechanism of action of insect ecdysteroids initiated with the early work of Clever and Carlson (for a his- torical review see Bellés and Santos, 2014), which eventually resulted in the so-called Ashburner model (Ashburner et al., 1974), which proposed a general model for the action of ecdysone, based on its participation in the regulation of gene expression (puffing) in the polytenic salivary gland chromosomes during Drosophila melanogaster molting and metamorphosis. Briefly, the model states that ecdysone associates with an intracellular recep- tor protein to activate early genes encoding transcription factors, which then activate late genes and, on the other, inhibit the transcription of previously activated early genes. The receptor protein and certain other members of this cascade belong to a large family of proteins, the nuclear hormone receptors (NR, see Fahrbach et al., 2012). NR proteins are generally comprised of four independent but functionally interacting domains. A/B is a highly variable domain that may contain a motif (AF-1) driving ligand-independent transcription. The second, the C domain, is a DNA-binding domain (DBD), the most conserved region of NRs. The D or hinge domain provides a link between DBD and the next domain, LBD, a multifunctional domain that mediates ligand binding, dimerization, and interaction with heat shock proteins, nuclear localization, and transactivation functions. Functional NRs form homodimers and/or heterodimers that recognize spe- cific DNA sequences. In the absence of a ligand molecule they act as repressors maintaining target genes inhibited by co-repressor complexes. In the presence of hormone they are activators of target genes by recruiting co-activator proteins and displacing co- repressors (Hill et al., 2013; Yamanaka et al., 2013; Evans and Mangelsdorf, 2014). The functional ecdysone receptor is a heterodimeric NR formed by the Ecdysone Receptor (EcR) and the ultraspiracle (usp) gene products (for a comprehensive review, see Hill et al., 2013). USP is an ortholog of the vertebrate retinoid-X receptor (RXR) (Yao et al., 1992) and is most commonly considered a kind of orphan NR. Though its ligand is not known, its participation as a mediator of JH action has been postulated, probably through a direct binding of JH (Barchuk et al., 2004). Furthermore, the EcR/USP complex can also bind to the let-7-C gene after a 20E pulse has triggered the larval-to-pupal metamorphic molt, thus inducing the transcription of a cluster of three microRNAs (miR-100, let-7 and miR-125). These then post-transcriptionally regulate the expression of genes involved in neuromuscular mor- phogenesis, leading to adult body characteristics (Chawla and Sokol, 2012; see also Rubio and Bellés, 2013). The EcR protein can also per se regulate the expression of target genes (Davis and Li, 2013), thus adding extra levels of complexity to the mechanisms and gene regulatory networks involving hormone/transcription factor activities. In insects showing caste polyphenism, there is evidence that ecdysteroids are important players in caste differentiation, not only during post embryonic, but possibly even during embryonic development (Schwander et al., 2008). The role of ecdysteroids in caste development and regulation of adult reproduction is cur- rently best understood in bees, especially so in the bumblebee Bombus terrestris (Geva et al., 2005) and in the honeybee Apis mellifera (Hartfelder and Engels, 1998), where they participate in the regulation of the differential morphogenesis programs by interacting with JH and possibly other mediating environmental modulators. Receptor proteins mediating ecdysteroids action in social insects have been studied mainly in the honeybee (The Honey Bee Genome Sequencing Consortium, 2006; Velarde et al., 2006), where USP and EcR cDNAs have been cloned (Barchuk et al., 2004; Takeuchi et al., 2007), and the expression profiles of the respective genes were determined in several organs, tissues, and conditions (Barchuk et al., 2004, 2008; Takeuchi et al., 2007; Velarde et al., 2009). However, and despite all these works, sev- eral responses to differential hormone signaling in honeybee caste development are still poorly understood (Barchuk et al., 2007). For instance, ecdysteroid titers in developing females are higher in queens during the second half of the last larval instar (Rachinsky et al., 1990) and differ in their profiles during pupal and pharate- adult development of queens and workers (Pinto et al., 2002). These hormone titer differences are associated with the differ- ential development of specific structures (e.g., brain and ovary, Barchuk et al., 2007) and also the onset of vitellogenin synthesis (Barchuk et al., 2002), but this is essentially correlative informa- tion lacking functional support. Herein we aimed at determining the extent to which EcR transcription follows ecdysteroids titers during honeybee postembryonic development and can actually mediate the action of molecular determinants of caste develop- ment in honeybees. Moreover, we expected that knocking-down EcR gene expression during pharate-adult development would bring to light new downstream targets of EcR. MATERIALS AND METHODS BEES Embryos and the successive developmental phases in the larval and pupal stages, as well as newly-emerged adults were obtained from A. mellifera colonies (Africanized hybrids) maintained at the Experimental Apiary of the University of São Paulo at Ribeirão Preto, Brazil. The developmental phases of workers and queens (Table 1) were identified according to Rembold (1987) and Michelette and Soares (1993). Immediately after sampling, the bees were immersed in TRIzol reagent (Life Technologies) and frozen at −80◦C until RNA extraction. NORTHERN BLOT ANALYSIS Approximately 15 μg of total RNA extracted from queens and workers at the PP1 and Pb developmental stages were subjected to electrophoresis in a denaturing 1.5% agarose/formaldehyde gel, and the RNA was then transferred to a PVF (Polyvinylidene Fluoride, GE) membrane using a VacuGene XL Vacuum Blotting Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 2 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 1 | Developmental stages (embryonic, larval, pupal, adult) of A. mellifera considered in this work. Abbreviation Developmental stage E Embryo L1 First instar larva L2 Second instar larva L3 Third instar larva L4 Fourth instar larva L5F1 Fifth instar larva, feeding phase 1 L5F2 Fifth instar larva, feeding phase 2 L5F3 Fifth instar larva, feeding phase 3 L5S1 Fifth instar larva, cocoon-spinning phase 1 L5S2 Fifth instar larva, cocoon-spinning phase 2 L5S3 Fifth instar larva, cocoon-spinning phase 3 PP1 Fifth instar larva, prepupa 1 PP2 Fifth instar larva, prepupa 2 PP3 Fifth instar larva, prepupa 3 Pw White-eyed pupa, unpigmented cuticle Pp Pink-eyed/pharate-adult transition, unpigmented cuticle Pdp Dark pink-eyed pharate-adult, unpigmented cuticle Pb Brown-eyed pharate-adult, unpigmented cuticle Pbl Brown-eyed pharate-adult, light pigmented cuticle Pbm Brown-eyed pharate-adult, intermediary pigmented cuticle Pbd Brown-eyed pharate-adult, dark pigmented cuticle NE Newly emerged adult system (GE Healthcare). An EcR cDNA fragment of 160 bp encoding the 3′ part of the DNA-binding domain was used for probe synthesis by means of the Random Primers DNA Labeling System (Life Technologies) and Redivue 32P-nucleotides (Amersham). After 3 h of hybridization at 42◦C, the membranes were washed during 20 min with 0.1 × SSC solution contain- ing 0.1% SDS and then exposed to a Super Sensitive ST film and the bands revealed with Cyclone™ Storage Phosphor System (PerkinElmer). HORMONE TREATMENTS For the analysis of the EcR expression response to artificially aug- mented levels of hormones, workers at the brown-eyed pupal phase (Pb) were removed from the brood frames and main- tained in an incubator at 34◦C and 80% relative humidity. For the ecdysone response, three groups of 3–7 workers were injected with 5 μg of 20-hydroxyecdysone (20E; Sigma) dissolved in 2 μL Ringer saline containing 12.5% ethanol. For the JH response, a similar number of Pb-phase workers received a topical applica- tion of 10 μg JH-III (Fluka) dissolved in 2 μL acetone. Controls received 2 μL of the respective solvents. The amounts of applied hormone were based on previous experiments in which we had examined their effects on inducing gene expression during pupal stage (Barchuk et al., 2002, 2004). RNA was isolated from fat bod- ies after 1, 12, and 24 h (independent experiments). Fat bodies were obtained via a longitudinal incision in isolated abdomens, which were then kept under gentle agitation in Petri dishes con- taining 0.9% NaCl. The resultant suspension of dispersed fat body cells was centrifuged during 1 min at 2500× g and the pellet was transferred into TRIzol reagent and frozen at −80◦C until RNA extraction. We used fat bodies because this allowed us to specif- ically assay this metabolically important organ, especially with regard to vitellogenin (vg) gene expression in honeybees. RNA EXTRACTION, REVERSE TRANSCRIPTION AND QUANTITATIVE PCR ASSAYS Total RNA was isolated using TRIzol (Life Technologies), fol- lowing the manufacturer’s protocol, and purified by column purification (RNeasy Mini Kit, QIAGEN), as described previously (Barchuk et al., 2004, 2007). For the quantification of mRNA lev- els (except those validating the RNA-Seq data), first strand cDNA was synthesized by reverse transcription from 2 μg of RNA with SuperScript II Reverse Transcriptase (Life Technologies) and an oligo(dT)12–18 primer (Life Technologies). For the validation of the RNA-Seq libraries, cDNA was synthesized using NCode™ miRNA First-Strand cDNA Synthesis and qRT-PCR (Invitrogen) kits and their instructions, adding a DNase (Promega) treatment step. Comparative analyses of transcript levels were performed by Real Time quantitative PCR (qPCR) using a 7500 Real-Time PCR System (Applied Biosystems) or a StepOne Plus system (Applied Biosystems). Amplifications were carried out in 20 μL reaction mixtures, each containing 10 μL of SYBR® Green Master Mix 2× (Applied Biosystems), 0.8 μL of a 10 mM stock solution of each of the gene-specific forward and reverse primers (Table S1), and 1 μL of first-strand cDNA diluted 1:4 (or 1:10, for cDNA samples used to validate RNA-Seq data) in ultrapure water. The sequences of forward primers were identical to the mature miRNA sequences available at miRBase, but replacing U by T, while the reverse Universal qPCR primer is supplied by NCode kit. Reaction conditions were 50◦C for 2 min, 95◦C for 10 min, followed by 40 cycles of 95◦C for 15 s and 60◦C for 1 min (or 33 s for miRNA amplification). Three biological replicates were run in three technical replicates each. Actin related protein 1 (Arp1, GenBank accession number NM_001185145.1), rpl32 (accession number NM_001011587.1), or a U5 snRNA gene were used as reference genes (for confirmation, cDNAs for all three reference genes were partially sub-cloned and sequenced in our labora- tory). Relative quantities of transcripts were calculated using the comparative Ct method (Applied Biosystems, User bulletin#2). Statistical analyses were carried out with Statistica version 7.0 (http://statistica.software.informer.com/). ASSESSING GENE TRANSCRIPTION PATTERNS ASSOCIATED TO EcR FUNCTION DURING HONEYBEE DEVELOPMENT USING RNAi dsRNA synthesis and treatment We employed a general protocol for dsRNA synthesis and injection in honeybees (Pbl phase) (Amdam et al., 2003). For EcR dsRNA synthesis, a 391 bp clone of EcR cDNA was amplified to serve as template, this comprising a fragment shared by the two transcript variants (A and B). The primers with the respective recognition site for T7 RNA polymerase (underlined) were: EcR-forward 5′-TAATACGACTCACTATAGGGCGAGAAT GGCGAGGAAGTACGAC and EcR-reverse 5′-TAATACGACT CACTATAGGGCGATTCTTGAACTTGAGGCTGAAG. A green fluorescent protein (GFP) gene clone was used as template to www.frontiersin.org December 2014 | Volume 5 | Article 445 | 3 http://statistica.software.informer.com/ http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development synthesize the respective dsRNA used as a non-target control (GFP-forward 5′-TAATACGACTCACTATAGGGCGAAGTGGA GAGGGTGAAGGTGA-3′ and GFP-reverse 5′-TAATACGACTC ACTATAGGGCGAGGTAAAAGGACAGGGCCATC-3′; see Nunes et al., 2013a). The amplification products were visualized and retrieved after agarose gel electrophoresis and purified using QIAquick™ (QIAGEN) columns. In vitro transcription reactions were performed by using the RiboMax™ T7 system (Promega) and the obtained dsRNA was isolated using TRIzol LS reagent (Invitrogen), subjected to a denaturation step at 98◦C for 5 min, followed by 30 min at room temperature, and diluted with nuclease free water to a final concentration of 2.5 μg/μL. The dsRNA quality was assessed by agarose gel electrophoresis. Pbl-phase workers (n = 30 for each experimental group) received an intra-abdominal injection of 2 μL of EcR dsRNA solution (2.5 μg/μL). Controls of the same developmental phase received the same volume of GFP dsRNA solution. dsRNA- injected bees were kept in an incubator at 34◦C and 80% relative humidity until adult eclosion (∼2 days), when they were trans- ferred to TRIzol reagent (Invitrogen) and frozen at −80◦C until RNA extraction. Total RNA extraction and first strand cDNA syn- thesis were carried out as described above. EcR knockdown effi- ciencies were assessed by RT-qPCR using variant-specific primers (EcRA-F, EcRB-F, and EcRA/B-R; see Table S1). Bees not used for gene expression analysis were used for evaluation of the adult phenotype. Analysis of gene expression patterns by RNA-Seq RNA pools of equal concentration from each group of EcR- and GFP-dsRNA treated bees were used for RNA-sequencing. Libraries were prepared using the TruSeq RNA™ Sample Preparation kit (Illumina) for poly-A+ RNA, and the TruSeq™ Small RNA Sample preparation kit (Illumina) for small RNAs (shorter than 200 nt). These libraries were shipped to the University of North Carolina (Chapel Hill, USA) facility where they were sequenced on an Illumina platform (Genome Analyzer II, Life Sciences). RNA-Seq reads for the poly-A+ RNA library were first submit- ted to adapter clipping using Scythe (Buffalo, 2011) (v.0.981— default parameters) for the 3′-end adapter and CutAdapt (Martin, 2011) (v.1.1—minimum overlap of 5 bp) for the 5′-end adapter. The next step was read trimming based on quality scores (mean Q ≥ 25), Ns (number of N bases lower than 10%) and poly-A tail prediction (minimum of 5 bp of A/T at both ends). This step was performed using PRINSEQ (v.0.19.5) (Schmieder and Edwards, 2011), which also filtered very small reads (length < 15 bp). An alignment against the A. mellifera genome (assembly version 4.5) was run using TopHat (Trapnell et al., 2009) (v.2.0.7), guided by the respective RefSeq (Release 55) transcript coordinates. The genomic alignments were then submitted to Cufflinks (Trapnell et al., 2010) (v.2.0.2) for transcript assembly, estimation of their abundances and testing for differential expression between EcR- KD and control samples. The Cufflinks procedures were also guided by the RefSeq transcript coordinates. The expression esti- mates were properly normalized considering ambiguous align- ments, and corrected for fragment bias (Roberts et al., 2011). The Poisson fragment dispersion model was used in the comparison analysis. Cufflinks calculates the FPKM (Fragments Per Kilobase of exon per Million fragments mapped), log2-fold-change and q- value (p-value adjusted by False Discovery Rate, FDR). However, the log2-fold-change was recalculated after adding an offset of 1 to FPKM values in order to enable comparison involving sam- ples without expression (zero) and to reduce the variability of the log ratios for low expression values (less than one). The func- tional annotation was done using Blast2GO (Conesa et al., 2005) (v.2.5), InterProScan (Mulder and Apweiler, 2007) (v.5-RC6), RefSeq transcript annotation and finding the Reciprocal-Best-Hit of A. mellifera RefSeq proteins against D. melanogaster proteins database (FlyBase r5.49) using blastp. The Blast2GO annotation pipeline was run based on blastp results of RefSeq proteins against nr database. Computational processing of the Small RNA-Seq reads com- prised the following steps: (i) initial sequence quality filtering based on unidentified bases; (ii) rRNA read filtering based on matches against SILVA database (Release 115); (iii) sequence adapter clipping using CutAdapt and Scythe; (iv) read trimming based on quality scores, Ns and poly-A+ tail prediction. All of these procedures were performed using PRINSEQ in the same way as described above. After each one of these preprocessing steps, an alignment against the A. mellifera genome (assembly ver- sion 4.5) was performed using the reads that had not already been aligned at each previous alignment step. Finally, all the alignment results were concatenated and transformed into a proper format to identify miRNAs. For this purpose, any splitted alignments were excluded. Genomic alignments were performed using TopHat and the other alignments were performed using Bowtie2 (Langmead and Salzberg, 2012) (v.2.0.6). miRNA digital expression (MDE) lev- els were obtained by analysis with miRDeep2 (Friedländer et al., 2012) (v.2.0.0.5∗), which provides the counts of reads mapped to the A. mellifera miRNA dataset in miRBase (Release 19). The original miRDeep2 code was modified to provide read counts for mature miRNAs instead of each precursor, and then the log2-fold-change was calculated and statistical significance was assessed using the method proposed by Audic and Claverie (1997) with adjustment by FDR. RESULTS EcR-A AND EcR-B TRANSCRIPT VARIANT IDENTIFICATION IN HONEYBEES Two transcript variants, EcR-A (Accession numbers NM_001098215.2) and EcR-B (NM_001159355.1) of 2635 and 2782 nucleotides, respectively, have been identified for the A. mellifera EcR gene (Takeuchi et al., 2007 and Watanabe et al., 2010). The difference in nucleotide length was shown to reside within the 5′ end, resulting in amino acid sequence variation in the N- modulator A/B domain. Conceptual translation of the nucleotide sequences resulted in a putative EcR-A protein consisting of 629 amino acid residues and an EcR-B protein of 557 amino acids, both sharing a 452 amino acid sequence in the carboxy terminal (Figure 1A). Northern blot analysis using a C-terminal EcR probe showed hybridization bands of approxi- mately 2.7 kb and 2.6 kb (Figure 1B), mainly in queen samples, but since we did not aim at quantifying, the respective band Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 4 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development FIGURE 1 | Gene and protein organization of honeybee EcR. (A) Genome: rectangles represent exons, and lines are introns, both are indicated with their respective numbers of nucleotides. Protein domains: Shown are the five domains with their respective amino acid numbers. The lenght of the coding sequence of each exon within the respective EcR domain is marked by oblique lines. 5′ and N termini are on the left. (B) Northern blot showing the expression of the two A. mellifera EcR transcript variants. Approximately 15 μg of total RNA was applied per lane. The radioactively labeled probe is a 160 bp fragment that included part of the DBD coding region. PP1, early prepupa; Pb: brown-eyed pharate-adult with unpigmented cuticle. density does not necessarily represent difference in transcript levels between the two castes. Nonetheless, this result reveals that the two transcripts indeed have small differences in length, this supporting the in silico evidence. DEVELOPMENTAL PROFILES OF THE EcR TRANSCRIPT VARIANTS A AND B Using variant-specific primers we quantified the transcript levels of EcR-A and EcR-B covering the entire postembryonic develop- ment for honeybee queens and workers (Figure 2). Three major findings are worthy of note: (i) transcripts representing the EcR- B variant are predominant in embryos (Mann–Whitney Test, P ≤ 0.05), but these transcript levels decline at the transition to the first larval instar, and it is the EcR-A variant which is then predominantly expressed during the end of the larval stage (fifth instar) and pupal stage; (ii) at several time-points, EcR expression is higher in workers than in queens (Mann–Whitney Test, P ≤ 0.05); and (iii) there is a clear discrepancy between circulating ecdysteroid levels and the developmental expression of EcR-A. Major caste differences in EcR-A expression were seen to accompany the larval/pupal metamorphic molt. As soon as the larvae were no longer fed by nurse bees and the brood cells were closed, the EcR-A levels were increased by two orders of magnitude in cocoon-spinning worker larvae (S1–S3 phases). A rise was also seen in EcR-A levels in cocoon spinning queen larvae, but this was significantly lower than in workers (Mann– Whitney Test, P ≤ 0.05). A similar pattern was also seen for the EcR-B variant, but at much lower modulation. Interestingly, the EcR expression levels were then decreased for both variants and in both castes at the onset of the prepupal development (PP1), marked by the appearance of apolysis fluid separating the fifth instar larval cuticle from the newly synthesized pupal cuti- cle in the head region. A new rise in the transcript levels of both variants was then seen at the end of the prepupal develop- ment (PP3), but this was primarily evident in workers (Mann– Whitney Test, P ≤ 0.05). EcR-A and EcR-B transcript variants remained at low levels during the pupal and early pharate-adult stages (Pw to Pbl phases) before they showed another steady increase, but again mainly so in workers (Mann–Whitney Test, P ≤ 0.05). TRANSCRIPTIONAL RESPONSE OF EcR TO ARTIFICIALLY AUGMENTED ECDYSTEROID AND JH TITERS So as to better understand the relationship between hemolymph hormone titers and hormone receptor expression, especially the remarkable divergence in the pupal stage, we treated Pb-phase workers and queens, as these are at the transition from pupal development per se to the pharate adult stage, with JH and 20E. At the Pb-phase the ecdysteroid titer is rapidly declining in both castes after having gone through the maximum peak at the pre- ceding Pp phase (Pinto et al., 2002), while JH levels are still basal (Rembold, 1987). The transcriptional responses for the two EcR variants assayed by RT-qPCR revealed a general repressive effect of both hormones at 24 h after application (Figure 3). In queens, 20E injection elicited a repressive effect on both EcR variants. www.frontiersin.org December 2014 | Volume 5 | Article 445 | 5 http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development FIGURE 2 | Developmental expression profile of the EcR gene in A. mellifera queens and workers. EcR-A and EcR-B transcript levels were measured by RT-qPCR. Bars represent means ± S.E.M. of three biological replicates, each run as three technical replicates. Asterisks indicate a statistical difference (Mann–Whitney Test, P ≤ 0.05) between queens and workers for the respective developmental phase. See Table 1 for a description of the developmental phases. Hormone titers are based on Rachinsky et al. (1990) and Pinto et al. (2002). Mean transcript levels were diminished at 12 h after 20E injection and were significantly lower at 24 h (Mann–Whitney Test, P ≤ 0.05). In workers this was the case only for the EcR-B transcript and only at 24 h (Figure 3A). The effect of exogenous JH on EcR expression was not as clear- cut as that elicited by 20E treatment. While there was no apparent effect on EcR-A transcripts in workers, the EcR-B levels showed slightly elevated means at all time points (Figure 3B), and these were significantly higher at 12 h following hormone treatment (Mann–Whitney Test, P ≤ 0.05). Interestingly, in the queen caste the response to JH treatment appeared to be opposite to that seen in workers, with mean EcR-A and EcR-B transcript levels dimin- ished already at 1 h after treatment and significant differences apparent at 1 h in the case of EcR-B and at 24 h for EcR-A (Mann– Whitney Test, P ≤ 0.05). These results indicate a repressor effect of high circulating ecdysteroid levels on EcR expression in both castes and a differential response to JH, with workers responding positively and queens negatively to elevated JH levels. EcR KNOCKDOWN IN PHARATE-ADULT HONEYBEE WORKERS SIGNIFICANTLY DOWNREGULATES THE EXPRESSION OF CANDIDATE TARGET GENES So as to understand the role of the EcR gene in honeybee devel- opment, beyond the correlation analysis between transcript levels of the two EcR variants and hormone levels, we experimentally decreased the EcR gene functionality by an RNA interference approach. We herein focused on the EcR response in workers during the pharate-adult to adult transition because only one of the two transcript variants, viz. EcR-A, undergoes a grad- ual increase at this developmental interval, and only so in the worker caste (Figure 2). We expected this to give not only more clear-cut results and insights into the role of the predominant EcR variant, but also into still very little understood aspects of morphogenetic processes taking place in developing adult honeybees. The dsRNA fragment used in this experiment represented an EcR region shared by the two transcript variants and its injection resulted in a reduction of 79.8 and 74.9% for EcR-A and EcR-B mRNA levels, respectively (P < 0.001, Student’s t- test; see Figure 4). A mortality of 10% was observed in both EcR- (KD) and GFP-dsRNA treated (control) bees. A proportion of dsRNA-injected bees showed alterations in cuticle pigmen- tation and wing development, similar to previously reported observations by Barchuk et al. (2008) when studying ultraspir- acle function. Based on the strong knockdown response we next assayed the transcriptional response of four candidate tar- get genes, these being a homolog of the D. melanogaster ftz-f1 gene, the vg gene, and two genes involved in adult cuticle for- mation (AmelCPR14 and BursA). The ftz-f1 gene was included in this analysis because in D. melanogaster it acts as a com- petence factor for the response to 20E; furthermore, EcR also inhibits ftz-f1 expression in D. melanogaster mid-prepupa, thus temporarily impairing the larval-to-pupal transition in response to the second 20E peak (King-Jones and Thummel, 2005). In pharate-adult honeybees, the levels of ftz-f1 transcripts were seen to increase (data not shown) concomitantly with the levels of EcR-A, suggesting a synergistic action of the two genes. In addi- tion, the increase in the levels of the two genes coincides with the increase in the expression of genes encoding enzymes and proteins needed for the complete differentiation of the adult cuticle (Soares et al., 2007, 2011, 2013; Elias-Neto et al., 2010). Similarly, in D. melanogaster the expression of ftz-f1 has recently been related to adult cuticle formation and eclosion (Sultan et al., 2014). The analysis of ftz-f1 transcript levels in newly emerged workers (N = 12), i.e., approximately 2 days after injecting dsR- NAs, showed that ftz-f1 expression was significantly decreased in EcR-KD bees (P ≤ 0.05, Student’s t-test) (Figure 4). A significant effect of the EcR knockdown was also seen for the cuticular pro- tein gene AmelCPR14, but not for the Burs A gene that encodes a subunit of the neurohormone Bursicon. The significant reduction in the expression of a cuticular protein gene following EcR-RNAi is consistent with the ecdysteroid-related expression of these genes in developing honeybees (Soares et al., 2007, 2011, 2013). An interesting though not easily explained finding was that vg gene expression was not significantly affected by reducing EcR func- tion, although the mean vg transcript levels were slightly reduced Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 6 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development FIGURE 3 | EcR expression response to experimentally augmented levels of (A) 20-hydroxyedysone (20E) and (B) juvenile hormone (JH) in honeybee castes. Aliquots of 5 μg of 20E or 10 μg of JH-III were applied to brown-eyed pharate-adults with unpigmented cuticle (Pb phase). RNA samples from fat bodies were obtained after 1, 12, and 24 h after hormone applications. Bars represent means ± S.E.M. of three samples. Each biological replicate consisted of 3–7 Pb-phase workers and each was run as three technical replicates. Statistical differences (∗P = 0.05) in gene expression between treated and control groups were assessed by Mann–Whitney Test. compared to the non-target dsRNA control. This was surprising as vg gene expression has been shown to gradually increase in pharate-adult honeybee females, and this increase was thought to be related to ecdysone levels (Barchuk et al., 2002; Piulachs et al., 2003). EcR KNOCKDOWN AFFECTS THE POLY-A+ PROFILE OF NEWLY EMERGED WORKERS So as to understand EcR functions during the pharate-adult to adult transition of honeybee workers on a more global scale we compared the poly-A+ transcriptomes of EcR-KD and GFP-injected (control) bees. After filtering of the raw data we obtained 112,659,148 reads for the KD and 71,050,536 reads for the control samples. Most of these reads were 50 nt long. This data has been submitted to the Sequence Read Archive (SRA, NCBI, http://www.ncbi.nlm.nih.gov/sra) under the Accession Number SRX700299. As we had only one RNA sample set per group (two libraries, no replicates), the estimate obtained by Cuffdiff analysis was that 234 loci were differentially expressed [absolute log2 (fold change) >1; q-value = 0.05; FPKM > 5 in at least one library] (Table S2). Among these, 121 code for known protein products, and 100 of these were upregulated in KD pharate-adults and 21 www.frontiersin.org December 2014 | Volume 5 | Article 445 | 7 http://www.ncbi.nlm.nih.gov/sra http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development FIGURE 4 | Relative transcript levels of four candidate EcR-target genes following EcR knockdown in honeybee workers. Twelve pharate-adult workers (Pbl phase) were injected with 5 μg of EcR-dsRNA (KD) or GFP-dsRNA (C, control) and sampled just after adult eclosion. Bars represent means and S.E.M of 12 biological replicates, each run as three technical replicates. ∗Indicates a statistically significant difference (Mann–Whitney Test, P ≤ 0.05). were downregulated (overexpressed in control bees; Table 2). The five times higher number of differentially expressed genes in the EcR-KD group indicates that during the pharate-adult to adult transition more genes may be repressed by ecdysone than are induced. In terms of functional assignments the following conclu- sions can be drawn. Seven genes among the ones upregulated in the KD group code for cytochrome P450 proteins (Table 2 and Figure S1A), six of these belonging to CYP clade 3 (CYP6AS2, CYP6AS3, CYP6AS4, CYP6AS5, CYP6AS12, and CYP6BD1) and one (CYP305D1) to CYP clade 2 [for clade assignments of hon- eybee cytochrome P450 genes see (Claudianos et al., 2006)]. A second protein family that was well-represented among the upregulated genes in the KD group is that encoding Major Royal Jelly Proteins (MRJP1 and MRJP9) and an MRJP-associated pro- tein, apisimin. A third class is represented by hormone response- related genes: a gene encoding a JH-inducible protein, a gene encoding a honeybee eclosion hormone (EH) homolog, and krüppel-homolog 1, an immediate response gene regulated by the JH receptor (Bellés and Santos, 2014). Nonetheless, the genes with the highest differential expression index are three genes encoding transcripts of unknown function and without conserved domain evidence (LOC100576540, LOC727013, and LOC727546). The fourth highest upregulated gene in the KD group encodes an α-glucosidase, an enzyme that converts the disaccharide sucrose into glucose and fructose and is, thus, critically involved in car- bohydrate metabolism. Another three genes in the top gene list are also related to metabolic functions, these being transcripts for a glycine N-methyltransferase-like, a glycine-methanol-choline (GMC) oxidoreductase 3 and a lipase 3-like protein. Furthermore, three genes upregulated in KD bees, the GMC oxidoreductase 3, a UDP-glycosyltransferase (LOC 413043) and a glucuronosyltrans- ferase (LOC 725997), could be related to ecdysteroid metabolism and function. The genes downregulated in EcR-KD bees are listed at the bot- tom of Table 2. They are represented with positive fold change values, as these were calculated as relative to the control group. In contrast to the upregulated genes, those that were downreg- ulated are not as clearly associated to putative functions during the pharate-adult to adult transition, except for the LOC724735 and Grp genes that encode structural cuticle proteins needed for the construction of the adult cuticle at this stage. The gene with the highest overexpression index in the control group codes for a Niemann-Pick type protein (NPC2), that is, genes involved in cholesterol metabolism-related syndromes and diseases (another npc2-type gene was found slightly overexpressed in the KD bees). Next are three transcripts possibly related to venom gland func- tion, encoding a phospholipase, secapin and a putative mast cell degranulating peptide (Table 2 and Figure S1A). Also down- regulated was the brown gene, which encodes an ABC-2 type transporter protein, and a gene coding for a Major Royal Jelly Protein (MRJP3). A more global analysis on the entire set of differentially expressed genes was done based on Gene Ontology (Blast2GO and InterProScan) using Fisher and Kolmogorov–Smirnov statis- tics. This confirmed that the poly-A+ RNAs representing genes upregulated in the KD group are enriched in proteins participat- ing in metabolic pathways, particularly ones with catalytic and oxidoreductase activities (Table S3). So as to validate the poly-A+ RNA-sequencing results we then chose two genes revealed as upregulated in the KD group (Cyp6as5, a P450 protein coding gene, and kr-h1, a gene encoding the JH response factor Krüppel homolog-1) and two down- regulated genes (LOC406145, secp and LOC724386, npc2). For these we designed or selected from the literature gene-specific primers and ran RT-qPCR assays. The expression pattern was confirmed for all four genes (Figure S1A), thus providing further evidence that the 234 poly-A+ RNA coding genes found as differ- entially expressed between treated and control bees are under EcR control. EcR KNOCKDOWN AFFECTS THE miRNA PROFILE OF NEWLY EMERGED WORKERS We obtained a total of 31,171,886 and 33,683,147 reads of small RNAs from the KD and control sequence libraries, respectively. This data has been submitted to the Sequence Read Archive (SRA, NCBI, http://www.ncbi.nlm.nih.gov/sra) under the Accession Number SRX700299. After filtering the raw data, we focused on the discovery of miRNAs linked to the EcR network. A total of 4,436,511 reads of the KD samples (∼13.2%) and 10,557,117 reads of the control sample (∼33.9%) mapped to known honey- bee mature miRNAs (available in miRBase version 19), suggesting that EcR disruption causes a general downregulation of miRNA families. We considered as “expressed” those miRNAs with more than 10 reads represented in at least one library. By doing so we retrieved a total of 132 known miRNAs expressed in newly emerged workers, most of them (124) in both conditions (Table S4). In order to find a set of miRNAs whose transcription is signif- icantly affected by the EcR pathway, we filtered the Cuffdiff results by selecting miRNAs with expression differences higher than 1.2- fold and a q-value < 0.05 between KD and control bees. We found 60 downregulated and 10 upregulated miRNAs in KD samples compared to controls (Table 3). These data were then further val- idated by RT-qPCR assays for the following miRNAs: miR-14, miR-100, miR-125, miR-133, miR-375, miR-3728, and miR-3771 (Figure S1B). Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 8 http://www.ncbi.nlm.nih.gov/sra http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 2 | Protein-coding genes (121) that were differentially expressed in EcR knockdown (EcR-KD) bees (FC ≥ 2; q-value ≤ 0.0015). Gene ID Symbol Annotation EcR-KD Control log2 FC FC 100576540 LOC100576540 Uncharacterized protein 102,312 14,520 −2.817 7.935 727013 LOC727013 Uncharacterized protein 1184,780 181,732 −2.705 7.316 727546 LOC727546 Uncharacterized protein 836,105 128,550 −2.701 7.297 409889 AGLU2 Alpha glucosidase 2 447,420 74,553 −2.585 6.684 724275 LOC724275 FBgn0051324 9486 1630 −2.541 6.455 409677 Cyp6as5 Cytochrome P450 6AS5 159,673 33,817 −2.239 5.014 551405 LOC551405 4-nitrophenylphosphatase-like 12,995 2764 −2.233 4.987 406093 LOC406093 Apisimin 943,098 205,225 −2.200 4.841 727213 LOC727213 FBgn0036592 11,369 2517 −2.175 4.732 552832 LOC552832 Glycine N-methyltransferase-like 97,380 23,780 −2.034 4.137 726965 LOC726965 CHK kinase-like; Protein kinase-like domain; Protein of unknown function DUF227; juvenile hormone-inducible protein 13,029 3194 −2.028 4.115 677664 Obp15 Odorant binding protein 15 44,643 11,157 −2.001 4.002 410747 GMCOX3 GMC oxidoreductase 3 14,032 3601 −1.962 3.850 409628 CDase Neutral ceramidase 90,558 23,540 −1.944 3.778 411353 LOC411353 Lipase 3-like 358,486 93,415 −1.940 3.764 406069 Kr-h1 Krüppel homolog 1 19,403 5207 −1.898 3.601 725165 LOC725165 Aquaporin-4-like 73,542 19,826 −1.891 3.577 412209 CYP6AS4 Cytochrome P450 6AS4, transcript variant 1 154,962 42,262 −1.874 3.514 552388 LOC552388 Major royal jelly protein 1-like 34,026 9661 −1.816 3.299 413908 CYP6AS12 Cytochrome P450 6AS12 71,569 20,875 −1.778 3.160 724312 LOC724312 Vanin-like protein 1-like 16,428 4810 −1.772 3.141 100578616 LOC100578616 Uncharacterized protein 65,980 19,436 −1.763 3.109 726377 LOC726377 Eclosion hormone-like 7230 2141 −1.756 3.083 100576979 LOC100576979 Apidaecins type 73-like 41,592 12,496 −1.735 3.010 411257 Hbg2 Alpha-glucosidase 213,375 66,083 −1.691 2.860 409709 LOC409709 Glucocerebrosidase 59,543 18,587 −1.680 2.821 724367 LOC724367 Protein lethal(2)essential for life-like 61,621 19,247 −1.679 2.818 726372 LOC726372 Trypsin-1-like 267,513 83,558 −1.679 2.818 411330 LOC411330 EF-hand domain; EF-hand-like domain; EF-Hand 1, calcium-binding site 16,958 5405 −1.649 2.721 726362 LOC726362 Luciferin 4-monooxygenase-like 18,728 6023 −1.637 2.679 408379 TpnCIIIa Troponin C type IIIa 384,660 126,632 −1.603 2.569 411188 LOC411188 L-lactate dehydrogenase 12,273 4046 −1.601 2.563 100578106 LOC100578106 Leucine-rich repeat, cysteine-containing subtype 198,220 65,924 −1.588 2.522 100576902 LOC100576902 Uncharacterized protein 7740 2611 −1.567 2.457 408788 LOC408788 Glucuronosyltransferase 13,777 4672 −1.560 2.434 724126 LOC724126 GNAT domain; Acyl-CoA N-acyltransferase 38,134 12,950 −1.558 2.428 551223 CYP305D1 Cytochrome P450 305D1 17,548 5976 −1.554 2.415 100576134 LOC100576134 FBgn0036665 37,976 13,190 −1.526 2.327 409873 Mrjp9 Major royal jelly protein 9 47,842 16,644 −1.523 2.320 552320 LOC552320 FBgn0263216 38,886 13,772 −1.497 2.242 409626 SP40 Serine protease 40 251,004 89,994 −1.480 2.190 408395 LOC408395 Venom carboxylesterase-6-like 39,507 14,192 −1.477 2.182 406142 LOC406142 Hymenoptaecin 1436,990 516,557 −1.476 2.179 (Continued) www.frontiersin.org December 2014 | Volume 5 | Article 445 | 9 http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 2 | Continued Gene ID Symbol Annotation EcR-KD Control log2 FC FC 726412 LOC726412 Solute carrier family 22 member 8-like, transcript variant 1; solute carrier family 22 member 8-like, transcript variant 2 9953 3585 −1.473 2.170 726737 LOC726737 Venom acid phosphatase Acph-1-like 7268 2629 −1.467 2.153 551560 CYP6BD1 Cytochrome P450 6BD1 106,633 39,084 −1.448 2.097 100577477 LOC100577477 Uncharacterized protein 5897 2167 −1.444 2.085 726611 LOC726611 Uncharacterized protein 75,931 28,272 −1.425 2.031 100576555 LOC100576555 Cytochrome b561, eukaryote; Cytochrome b561/ferric reductase transmembrane 36,621 13,730 −1.415 2.003 725381 LOC725381 Uncharacterized protein 163,884 62,397 −1.393 1.941 409716 apd-3 Apidermin 3 70,996 27,036 −1.393 1.940 724239 LOC724239 Kynurenine/alpha-aminoadipate aminotransferase, mitochondrial-like 8440 3274 −1.366 1.867 100576895 LOC100576895 Fatty acyl-CoA reductase 10,988 4279 −1.361 1.851 100578995 LOC100578995 Vanin-like protein 1-like 5902 2301 −1.359 1.846 413043 LOC413043 UDP-glycosyltransferase 293,945 116,620 −1.334 1.779 413575 LOC413575 Facilitated trehalose transporter Tret1-like 128,316 51,327 −1.322 1.747 725114 LOC725114 Trypsin Inhibitor-like, cysteine rich domain 6914 2772 −1.319 1.739 725273 CTL1 C-type lectin 1 43,402 17,672 −1.296 1.680 725204 LOC725204 Tyrosine aminotransferase-like 8411 3454 −1.284 1.648 726934 SPH50 Serine protease homolog 50 13,552 5672 −1.257 1.579 724993 LOC724993 FBgn0036202 213,583 89,752 −1.251 1.564 725997 LOC725997 Glucuronosyltransferase 26,069 10,989 −1.246 1.553 551458 LOC551458 Leucine-rich repeat-containing protein 20-like, transcript variant 2; leucine-rich repeat-containing protein 20-like, transcript variant 1 107,500 46,024 −1.224 1.498 724756 LOC724756 Gadd45 92,117 39,653 −1.216 1.479 100578635 LOC100578635 Uncharacterized protein 822,735 356,743 −1.206 1.453 413117 LOC413117 Proton-coupled amino acid transporter 4-like 41,718 18,176 −1.199 1.437 552366 LOC552366 Hypothetical LOC552366 15,048 6559 −1.198 1.435 724721 LOC724721 Dehydrogenase/reductase SDR family member 11-like 356,294 158,024 −1.173 1.376 100578329 LOC100578329 Putative fatty acyl-CoA reductase CG5065-like 44,956 19,991 −1.169 1.367 412458 LOC412458 Dehydrogenase/reductase SDR family member 11-like 337,026 151,686 −1.152 1.327 726803 LOC726803 FBgn0037126 305,241 139,173 −1.133 1.284 406115 Apid73 Apidaecin 38,277 17,505 −1.129 1.274 552301 SP35 Serine protease 35 156,146 71,778 −1.121 1.257 724654 LOC724654 Cytochrome b5 type B-like 285,205 131,483 −1.117 1.248 409740 LOC409740 Clavesin-1-like 18,383 8491 −1.114 1.242 724899 Lys-2 Lysozyme 2 15,635 7245 −1.110 1.232 552193 LOC552193 Proton-coupled amino acid transporter 4-like 12,681 5892 −1.106 1.223 406065 Wat Worker-enriched antennal transcript 187,038 87,451 −1.097 1.203 (Continued) Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 10 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 2 | Continued Gene ID Symbol Annotation EcR-KD Control log2 FC FC 408622 LOC408622 Phenylalanine hydroxylase, transcript variant 2 19,241 9009 −1.095 1.199 725038 LOC725038 Protein npc2 homolog, similar to that of GeneID724386 619,928 290,849 −1.092 1.192 551094 LOC551094 Fatty-acid amide hydrolase 2-A-like, transcript variant 2 6524 3061 −1.092 1.191 408807 LOC408807 Hypothetical LOC408807 22,855 10,746 −1.089 1.185 409638 LOC409638 Elongation of very long chain fatty acids protein AAEL008004-like 112,535 52,938 −1.088 1.184 406109 Obp6 Odorant binding protein 6 12,054 5672 −1.088 1.183 727522 LOC727522 Insect allergen-related 16,753 7972 −1.071 1.148 100577337 LOC100577337 Glucose dehydrogenase 6428 3072 −1.065 1.134 100577132 LOC100577132 Calcium calmodulin-dependent protein kinase kinase 2 7270 3491 −1.058 1.120 408299 LOC408299 Purine nucleoside phosphorylase-like 159,071 76,456 −1.057 1.117 408733 LOC408733 Pinocchio 210,857 101,474 −1.055 1.113 727367 LOC727367 Glucose dehydrogenase [acceptor]-like 29,151 14,038 −1.054 1.111 410087 LOC410087 Protein lethal(2)essential for life-like, transcript variant 1 101,566 48,915 −1.054 1.111 100576662 LOC100576662 Uncharacterized protein 60,105 29,083 −1.047 1.097 411615 CYP6AS2 Cytochrome P450 6AS2 76,760 37,313 −1.041 1.083 551044 Gld2 Glucose dehydrogenase 2, transcript variant 1 5627 2756 −1.030 1.060 726871 LOC726871 Synaptic vesicle glycoprotein 2C-like 77,575 38,301 −1.018 1.037 406144 LOC406144 Abaecin 979,123 484,043 −1.016 1.033 408868 LOC408868 Inositol-1-(or 4-)monophosphatase 87,544 43,471 −1.010 1.020 408421 LOC408421 CHK kinase-like; Protein kinase-like domain; Protein of unknown function DUF227 391,146 194,281 −1.010 1.019 725344 LOC725344 Histone H2B.3-like 13,796 6863 −1.007 1.015 726690 CYP6AS3 Cytochrome P450 6AS3 34,408 17,138 −1.006 1.011 725217 LOC725217 Armadillo-type fold 9786 19,784 1.016 1.031 552421 LOC552421 Glycogenin-1-like 37,050 75,154 1.020 1.041 100577901 LOC100577901 FBgn0030763 18,696 38,195 1.031 1.062 100578838 LOC100578838 Chymotrypsin inhibitor-like 48,355 99,598 1.042 1.087 100579045 LOC100579045 Uncharacterized protein 3038 6320 1.057 1.116 412399 LOC412399 Organic cation transporter protein-like 8913 19,674 1.142 1.305 412797 LOC412797 Facilitated trehalose transporter Tret1-like, transcript variant 1 7785 17,708 1.186 1.406 100578811 LOC100578811 Transmembrane protein 223-like 7839 18,013 1.200 1.441 724735 LOC724735 Endocuticle structural glycoprotein SgAbd-2-like 2941 7068 1.265 1.600 678674 LOC678674 Venom allergen Api m 6 1994 5014 1.330 1.769 552281 Grp Glycine-rich cuticle protein 7302 18,853 1.369 1.873 100578552 LOC100578552 Chitin binding domain 4141 10,943 1.402 1.965 406145 LOC406145 Secapin 32,484 92,219 1.505 2.266 100578873 LOC100578873 Allergen Api m 6-like 1799 5256 1.547 2.393 100576769 LOC100576769 Mast cell degranulating peptide-like 31,784 95,466 1.587 2.518 406141 Pla2 Phospholipase A2 5068 15,586 1.621 2.627 (Continued) www.frontiersin.org December 2014 | Volume 5 | Article 445 | 11 http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 2 | Continued Gene ID Symbol Annotation EcR-KD Control log2 FC FC 100578863 LOC100578863 Uncharacterized protein 1641 5082 1.631 2.660 725163 LOC725163 Trypsin Inhibitor-like, cysteine rich domain 5520 22,669 2.038 4.153 406121 Mrjp3 Major royal jelly protein 3 1102 5592 2.343 5.491 412203 bw Brown 1347 7074 2.392 5.724 724386 LOC724386 Protein npc2 homolog 757,358 4202,210 2.472 6.111 Negative values of log2 FC indicate genes that are upregulated in EcR-KD bees (highlighted in red); positive values indicate downregulation in EcR-KD bees (highlighted in green). Genes shown in bold had their transcription patterns validated by qPCR. Expression values measured as FPKM (Fragments Per Kilobase of exon per Million fragments mapped). This list contains only the genes with FPKM values of 5 for any of the two samples. DISCUSSION THE HONEYBEE EcR TRANSCRIPT VARIANTS AND THEIR DEVELOPMENTAL REGULATION The existence of more than one EcR isoform is commonplace in insects, including the honeybee, for which two transcript variants, EcR-A and EcR-B had been found (Takeuchi et al., 2007). First shown for D. melanogaster (Talbot et al., 1993) and then for the red flour beetle Tenebrio molitor (Mouillet et al., 1997), the extensive review of insect EcR isoforms by Watanabe et al. (2010) showed a high similarity in their nucleotide and amino acid sequences in most of their functional domains, except for the N-terminal region including the variable A/B modula- tor domain, which might allow for the recruitment of different co-activators/co-repressors (Tora et al., 1988; Kato et al., 1995; Watanabe et al., 2010). First, we confirmed by northern blotting the expression of the two EcR variants in honeybee queens and workers. Then, we compared their temporal expression profiles to the hemolymph ecdysteroid titers of fifth instar queen and worker larvae (F1-PP3 phases) (Rachinsky et al., 1990). The results for the developmen- tal expression profiles of the two ecdysteroid receptor variants are surprising in two aspects. First, contrasting with the hormone titers, which are higher in queens than in workers, the EcR tran- script levels were found to be higher in workers, especially so for EcR-A. Second, there was a marked drop in EcR expression at the beginning of the prepupal phase (PP1), i.e., exactly when the hemolymph ecdysteroid levels increase to reach a develop- mental peak at the subsequent PP2 phase. Strikingly as well, the transcript levels for both EcR variants remained at low or basal levels during the pupal and early pharate-adult stages (Pw to Pbl phases), even though the ecdysteroid hemolymph titers are at a maximum during this period (Feldlaufer et al., 1985; Pinto et al., 2002). The switch from EcR-B expression in the embryonic stage to EcR-A as the predominant isoform in the fifth larval instar and pupal stage reflects a change in the processing of an eventual long pre-mRNA, or a shift in transcription start site utilization (our RNA-Seq data are in support of the latter possibility and even suggest the existence of a third EcR transcript variant). Since the EcR gene is known to be induced after an ecdysteroid pulse (Karim and Thummel, 1992; Davis and Li, 2013), the production of EcR-B mRNA in honeybee embryos would require the pres- ence of steroid hormones, which is indeed the case. Makisterone A, the predominant ecdysteroid in A. mellifera (Feldlaufer et al., 1985), has been shown to be present in ovaries in quite large amounts (Feldlaufer et al., 1986a), and unpublished data from our laboratory also confirm the presence of ecdysteroids in devel- oping embryos. High levels of ecdysteroids in ovaries have also been shown for bumblebee queens (Geva et al., 2005) and queens of a swarm-founding neotropical wasp, Polybia micans (Kelstrup et al., 2014). Embryonic ecdysteroids can be synthesized by enzy- matic conversion from inactive conjugates stored during oogen- esis (Dorn, 2000) or, as seen in mosquitoes, transferred by males during copulation (Baldini et al., 2013). Since makisterone A is the predominant ecdysteroid com- pound in queen ovaries (Feldlaufer et al., 1986a) and also in pupal-stage hemolymph (Feldlaufer et al., 1986b) and 20E is not negligible in prepupal hemolymph (Rachinsky et al., 1990), the observed embryonic-to-larval EcR isoform switch may be linked to variation in the ecdysteroid composition circulating in the hemolymph, throughout a bee’s life cycle. This could not only be responsible for the observed differential EcR transcription, but also for the formation of different hormone/receptor complexes with potentially different target genes thus, governing separate physiological processes. 20E, for example, might have retained a role in reproductive physiology, as suggested by Takeuchi et al. (2007), whereas makisterone A may have been co-opted for gov- erning postembryonic development, as suggested for Dysdercus fasciatus (Feldlaufer et al., 1991). Nonetheless, for honeybees such “division of labor” in ecdysteroid compounds is still highly speculative, especially since the ecdysteroid hemolymph levels in adult honeybee queens and workers are continuously low, this making it rather unlikely that these steroid hormones may play a major role in the reproductive female physiology (Hartfelder et al., 2002). Instead, they seem to be preferentially stored in the developing follicles. The second and third major findings mentioned above are that the EcR-A transcript levels are higher in worker than in queen development, and that there is no positive, but rather an appar- ently negative correlation between hormone levels and hormone receptor transcript levels. This stands in stark contrast to the developmental pattern of the hemolymph ecdysteroid titers in the two castes, which are higher in queens than in workers, partic- ularly so during larval-pupal metamorphosis (Rachinsky et al., 1990). The ecdysteroid titer in last instar queen larva rises as soon as the brood cells are closed and the larvae start to spin their Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 12 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 3 | miRNAs that were differentially expressed in EcR-knockdown bees. Effect miRNA Read counts Normalized number Fold change of reads Control KD Control KD D ow nr eg ul at ed ame-miR-3771-3p 12 1 0.79 0.13 6.06 ame-miR-6067-5p 19 2 1.25 0.25 4.99 ame-miR-6046-3p 31 6 2.04 0.76 2.68 ame-miR-6057-5p 81 17 5.34 2.16 2.48 ame-miR-100-5p 1,117,633 257,905 73,616.85 32,712.83 2.25 ame-miR-3770-5p 57 14 3.75 1.78 2.11 ame-miR-133-3p 7726 2092 508.90 265.35 1.92 ame-miR-375-3p 64,153 17,597 4225.66 2232.01 1.89 ame-miR-1-3p 78,433 21,559 5166.27 2734.56 1.89 ame-miR-6043-3p 352 98 23.19 12.43 1.87 ame-miR-2765-5p 1903 548 125.35 69.51 1.80 ame-miR-125-5p 118,218 35,904 7786.85 4554.09 1.71 ame-miR-6047a-3p 298 92 19.63 11.67 1.68 ame-miR-927a-5p 129,279 40,144 8515.42 5091.89 1.67 ame-miR-971-3p 110 35 7.25 4.44 1.64 ame-miR-87-3p 14,514 4783 956.02 606.68 1.58 ame-miR-263a-5p 857,064 283,563 56,453.55 35,967.31 1.57 ame-miR-137-3p 3351 1119 220.73 141.93 1.56 ame-miR-1175-3p 2179 731 143.53 92.72 1.55 ame-miR-210-3p 3330 1121 219.34 142.19 1.55 ame-miR-316-5p 21,370 7287 1407.61 924.29 1.53 ame-miR-6038-5p 445 158 29.31 20.04 1.46 ame-miR-219-5p 70 25 4.61 3.17 1.45 ame-miR-980-3p 235 86 15.48 10.91 1.41 ame-miR-92a-3p 190 70 12.52 8.88 1.41 ame-miR-92b-3p 42,256 15,578 2783.34 1975.92 1.40 ame-miR-279a-3p 55,440 20,620 3651.75 2615.45 1.39 ame-miR-279c-3p 6610 2482 435.39 314.82 1.39 ame-let-7-5p 119,020 44,921 7839.67 5697.81 1.38 ame-miR-989-3p 528 201 34.78 25.49 1.37 ame-miR-252a-5p 19,053 7232 1254.99 917.31 1.37 ame-miR-184-3p 1,460,229 568,579 96,183.15 72,118.91 1.34 ame-miR-6041-3p 95 37 6.26 4.69 1.34 ame-miR-252b-5p 420 164 27.66 20.80 1.33 ame-miR-3747b-5p 793 310 52.23 39.32 1.33 ame-miR-2788-3p 1137 445 74.89 56.44 1.33 ame-miR-317-3p 6091 2391 401.21 303.28 1.32 ame-miR-6012-3p 122 48 8.04 6.09 1.32 ame-miR-3791-3p 697 277 45.91 35.13 1.31 ame-miR-71-5p 1667 660 109.80 83.71 1.31 ame-miR-34-5p 3245 1289 213.74 163.50 1.31 ame-miR-3718a-3p 3817 1517 251.42 192.42 1.31 ame-miR-927b-5p 3108 1250 204.72 158.55 1.29 ame-miR-929-5p 3423 1379 225.47 174.91 1.29 ame-miR-305-5p 62,480 25,699 4115.47 3259.68 1.27 ame-miR-124-3p 2142 882 141.09 111.87 1.26 ame-miR-996-3p 44,757 18,495 2948.08 2345.92 1.26 ame-miR-10-5p 1,146,362 477,869 75,509.19 60,613.20 1.25 ame-miR-276-3p 808,220 339,646 53,236.27 43,080.91 1.24 (Continued) www.frontiersin.org December 2014 | Volume 5 | Article 445 | 13 http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Table 3 | Continued Effect miRNA Read counts Normalized number Fold change of reads Control KD Control KD ame-miR-3477-5p 46,537 19,658 3065.32 2493.43 1.23 ame-miR-6001-5p 173 73 11.40 9.26 1.23 ame-miR-283-5p 42,234 17,884 2781.89 2268.42 1.22 ame-miR-263b-5p 41,279 17,642 2718.99 2237.72 1.21 ame-miR-14-3p 94,785 40,489 6243.35 5135.65 1.21 ame-miR-79-3p 3701 1589 243.78 201.55 1.21 ame-miR-3719-3p 2799 1217 184.37 154.37 1.20 ame-miR-12-5p 89,841 38,948 5917.70 4940.19 1.20 ame-miR-279b-3p 650 281 42.81 35.64 1.20 ame-miR-6039-5p 494 214 32.54 27.14 1.20 ame-miR-6051-3p 343 149 22.59 18.90 1.20 U pr eg ul at ed ame-miR-3728-5p 105 124 6.92 15.73 2.27 ame-miR-965-5p 21 20 1.38 2.54 1.84 ame-miR-6005-3p 76 64 5.01 8.12 1.62 ame-miR-3798-3p 342 259 22.53 32.85 1.45 ame-miR-6052-5p 43 31 2.83 3.93 1.39 ame-miR-3761-5p 342 232 22.53 29.43 1.31 ame-miR-3720-5p 1586 1004 104.47 127.35 1.22 ame-miR-6001-3p 27,917 17,745 1,838.85 2250.79 1.22 ame-miR-306-5p 186,621 117,263 12,292.45 14,873.71 1.21 ame-miR-9a-5p 88,013 55,087 5797.29 6987.27 1.21 Only miRNAs with Fold-change >1.2 are listed. cocoons (S1-stage), while in workers this was only seen in the late spinning (S3) to early prepupal (PP1) phases. Furthermore, the peak in edysteroid levels reached during the prepupal phase (PP2) is twice as high in queens compared to workers (Rachinsky et al., 1990). The negative correlation between ecdysteroid lev- els and EcR expression is particularly evident at two time points: in prepupae, when ecdysteroid levels are high in the PP1–PP2 phases, just as the EcR-A expression pattern undergoes a valley, and in the pupal stage, when the ecdysteroid levels are high in both castes at the Pp phase, before dropping in the Pb-Pbl phases (Pinto et al., 2002). It is only after this drop in circulating hor- mone levels that EcR-A transcription is resumed, particularly so in the worker caste, and strikingly, it is during the subsequent Pbm and Pbd phases that the ecdysteroid levels are again lower in workers than in queens (Pinto et al., 2002). This apparent negative correlation between hormone and hor- mone receptor levels was suggestive of a repressive action of high concentrations of circulating ecdysteroids on the expression of their receptor gene. To test this we manipulated the endogenous hormone levels by treating Pb pharate-adults with either 20E or JH. The results of the 20E injection experiments showed that a prolonged and excessive presence of ecdysteroids had a repressive effect on EcR-A and EcR-B expression, especially so in queens (Figure 3A). Interestingly, workers seem to be more resilient to this repressor effect, as there was no significant reduction in EcR- A transcript levels, comparable to that seen in queens, or for EcR-B in both castes. Such resilience was also denoted in the JH application experiment, where EcR-A transcript levels in work- ers remained little affected compared to those in queens and for EcR-B in both castes 24 h after the 20E injection. Strikingly, JH appeared to have opposite effects on EcR-B expression in the two female castes, showing a positive effect in workers and a negative one in queens. These differences of hormone effects on EcR- expression related to caste certainly deserve a closer look in future experiments. Repressive effects of high concentrations of ecdysteroids on EcR-expression are, however, not new and are likely to be a gen- eral feature of hormone systems that underly cyclical events in morphogenesis and physiology. For instance, similar results were described for Manduca sexta, where low concentrations of 20E induced EcR expression while high concentrations repressed the expression of this gene (Jindra et al., 1996). Like ours, these results suggest that the EcR gene responds positively to a slight increase in ecdysteroids, whereas high hormone levels are repressive. In fact, as we could see, EcR expression actually appears to precede the rise in hormone levels, for instance in the S1–S3 phases, when the circulating ecdysteroid levels start increasing (Rachinsky et al., 1990), but EcR-A and also EcR-B transcript levels have already undergone a steep rise. A repetition of this pattern can be inferred for the pupal ecdysis event, occurring between PP3 and Pw, when the ecdysteroid titers undergo a sharp drop, but EcR-A and EcR-B are on the rise (mainly in workers), and drop once the ecdysteroid Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 14 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development titers build up to maximal values in the Pp phase. It is such cyclical events, the molts, that are synchronized by the ecdysone/ecdysone receptor complex action, and this is primarily seen in the epider- mis, the main organ of cuticle synthesis. In the honeybee, several cuticle protein genes were shown to be regulated by ecdysteroids (Soares et al., 2007, 2011, 2013; Elias-Neto et al., 2010). RNAi-MEDIATED KNOCKDOWN REVEALS EcR REGULATED GENES IN DEVELOPING ADULTS Upon comparing the sequencing results of the poly-A+ libraries for EcR knockdown (EcR-KD) and control groups, the Cutdiff analysis classified 234 loci as differentially expressed. Among these, 121 were annotated as coding for known protein products or, from another point of view, 113, i.e., one half, represent loci for unknown, not annotated products, which could be either pro- teins or long non-coding RNAs. Especially the latter are still “dark matter” in the honeybee genome, represented by many ESTs in the databases, but only four long non-coding RNAs are so far characterized to some detail (Sawata et al., 2002; Humann et al., 2013). Among the genes with known orthologs or sequence similar- ity in functional domains, 100 were overexpressed (fold change > 1) in the EcR-KD group and 21 in the control group, this indicating that apparently more genes are repressed by the ecdysone/EcR receptor complex than are activated. Furthermore, a Gene Ontology and KEGG pathway analysis showed that there is little overlap in gene functions between the two sets of differen- tially expressed genes (DEGs). As mentioned above, cytochrome P450 genes are strongly represented among the DEGs. While cytochrome P450 genes are a large gene family, strongly related to detoxification processes, this family has undergone consider- able reduction in honeybee genome evolution (Claudianos et al., 2006). This reduction is, however, denoted only in certain clades of the P450 enzymes, but not in the clades comprising the genes found in our EcR-RNAi experiment. Unfortunately, there is no further functional or tissue/cell type information available for the five cytochrome P450 genes, especially whether or not they may be related to steroid synthesis or metabolism. Nonetheless, sim- ilar findings as the ones we report here were also denoted by Davis and Li (2013) in their genomic screen for ecdysone and EcR-dependent gene expression in D. melanogaster. A second group of overrepresented genes that called attention was the hormone response-related genes, as these may provide a link between JH and ecdysteroid action during the pharate-adult to adult transition in honeybees. For this group we found three genes as overexpressed in the EcR-KD group, viz. a JH-induced protein, kr-h1, and an Eclosion hormone-like (EH-like) gene. kr- h1 is certainly the most interesting gene in this set, as it represents a direct readout of the activity of the JH response in target tissues (Lozano and Belles, 2011; Bellés and Santos, 2014). As kr-h1 has previously been identified in a screen for ecdysone-response genes in D. melanogaster (Beckstead et al., 2005), the current identifi- cation of this gene in the EcR-KD group provides experimental evidence toward a mechanistic explanation for the modulation of vitellogenin induction in honeybee pharate-adults, where vg expression is caste-specifically induced by JH and counteracted by ecdysteroids (Barchuk et al., 2002). Overexpression of an EH-like gene in the EcR-KD group was not unexpected, as EH is syn- thesized in response to declining ecdysteroid titers and is part of the ecdysis triggering signaling cascade (Zitnan and Adams, 2012). Interestingly, other three upregulated genes in EcR-KD bees may have roles in ecdysteroid metabolism and function. Several GMC oxidoreductase genes in diverse insects, including A. mellifera, are clustered in an evolutionary conserved tandem array with potential to be co-regulated for a common function related to ecdysteroid metabolism (Iida et al., 2007). The prod- ucts of LOC413043 and LOC725997 may regulate ecdysteroid titer and function since the enzymes encoded by these genes cat- alyze the transfer of glucose from UDP-glucose to ecdysteroids, and thus are possibly related to ecdysteroid inactivation (O’Reilly and Miller, 1989). Overexpression of members of the Major Royal Jelly Protein (MRJP) family can be interpreted in the context of a repressive action of the ecdysone/EcR receptor complex on genes of the adult honeybee life cycle (the only exception being the gene cod- ing for the MRJP3, which was overexpressed in control bees). The mrjp gene family with its nine members is a lineage-specific exten- sion in the genus Apis, from a single mrjp-like gene within the yellow genes complex (Drapeau et al., 2006). Even though these proteins are highly expressed in the hypopharyngeal glands of nurse worker bees, constituting the major protein fraction of the glandular secretions fed to larvae (royal jelly and worker jelly), expression of the mrjp genes is neither exclusive to this tissue nor is it restricted to the worker caste. Especially mrjp9 has been shown to be broadly expressed, in different tissues of adult work- ers and also in queens and even drones (Buttstedt et al., 2013). In contrast to mrjp9, mrjp1 expression is more tissue-specific, being highest in heads (viz. hypopharyngeal glands) of nurse bees, with expression levels being considerably lower in other body parts, castes and sexes (Buttstedt et al., 2013). MRJP1 is the predomi- nant MRJP moiety in royal jelly, present as oligomers of MRJP1 subunits, which are held together by apisimin, a small 5 kDa pro- tein (Tamura et al., 2009). ESTs corresponding to apisimin were found as overrepresented in the EcR-KDS group, indicating that its expression is co-regulated with that of mrjp1. But this co- regulation is not due to genomic proximity, as the mrjp/yellow gene cluster maps to chromosome 11, whereas apisimin is located in chromosome 6. Interestingly, an MRJP1 monomer, royalactin, was found to be an important factor in caste development, acting through the Egfr signaling pathway (Kamakura, 2011). Among the EcR-KD group genes we also identified obp15, which encodes a putative odorant binding protein. Some obp genes were also found to be under negative EcR control in D. melanogaster, including the obp15 and obp6 genes (Davis and Li, 2013). Forêt and Maleszka (2006) had previously shown that obp15 is expressed in the antennae of adult bees and also in young larvae but not in pupae. The high ecdysteroid levels in honey- bee hemolymph during the pupal to pharate-adult transition, thus, appear to repress obp15 expression, and possibly also other members of this complex gene family. Among the genes overrepresented in the transcriptome of the control group (downregulated in EcR-KD bees), the first in the top ten list is annotated as npc2. Genes of this family are associated with Niemann-Pick syndromes and diseases affecting www.frontiersin.org December 2014 | Volume 5 | Article 445 | 15 http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development cholesterol metabolism (Carstea et al., 1997). In D. melanogaster, NPC mutations cause intracellular enrichment of cholesterol, reduced ecdysteroidogenesis and death in the first larval instar. The fact that this condition could be fully rescued when an excess of dietary cholesterol was given to these mutants indicated that the ecdysone biosynthesis pathway is intact, but precursor processing is not (Huang et al., 2007). Interestingly, in honey- bees, as in other insects, the major ecdysteroid is not ecdysone or its derivative 20E, but makisterone A, an ecdysteroid methy- lated at C24 (Feldlaufer et al., 1985, 1986a,b; Rachinsky et al., 1990), possibly due to a lack or restriction in C24-demethylation of a phytosterol precursor. The expression of two other genes overrepresented in the transcriptome of the control group has been shown to be dependent on the ecdysteroid titer. The pro- tein encoded by LOC724735, an endocuticle structural protein (Márcia M. G. Bitondi, unpublished results) and also the Grp gene, renamed as tweedle1 (AmelTwdl1) (Soares et al., 2011), were induced in the integument by the ecdysteroid pulse that promotes the pupal to pharate-adult transition. Thus, the functionality of the ecdysone/EcR complex is necessary for the activation of these genes. Interestingly, mrjp3, the third among the genes overrepresented in control bees (and thus induced by the EcR pathway), encodes one of the main MRJPs produced by nurse bees (Buttstedt et al., 2013). The mrjp3 gene thus seems to be highly expressed by the time of adult emergence and the first days of adult life. Unlike the mrjp1 and mrjp2 genes, mrjp3 reaches negligible expression levels in foragers, which, together with its distinctive amino acid sequence (Drapeau et al., 2006), supports the notion of its main function as food protein supplier to lar- vae by nurse bees. However, the fact that mrjp1, another MRJP gene highly expressed in nurse bees, was found to be repressed by the ecdysteroid pathway (see above), suggests the mrjp3 gene is regulated in a distinct mode from the other mrjp genes. miRNAs AS ACTORS IN THE EcR REGULATORY NETWORK Here we demonstrate that the RNAi-mediated knockdown of EcR function perturbs the expression of 70 miRNAs (∼1/3 of the hon- eybee miRNAs known to date). Most of these (60) were downreg- ulated and 10 were upregulated and we assume that these down and upregulated miRNAs are “induced” or “repressed,” respec- tively, by the EcR pathway as bees undergo the pharate-adult to adult transition. Among the miRNAs that showed significant changes in abun- dance following EcR knockdown, most had already been identi- fied in a large-scale sequencing project (Chen et al., 2010), but had no function(s) associated. Our data now lead to infer that these miRNAs are, at least, closely associated with EcR action and, consequently, connected to pupal-adult metamorphosis. In addi- tion to these miRNAs of yet unclear functions, we also found conserved and functionally well-defined miRNAs, such as let-7, miR-1, miR-133, miR-375, miR-184, and miR-34. For example, miR-133 and miR-1 are both clustered in the mouse and fly genomes, and they play important roles in muscle development and differentiation in vertebrates and invertebrates (Sokol and Ambros, 2005; Chen et al., 2006; Boutz et al., 2007). In the honey- bee, however, we found these two miRNAs to be located far apart from one another on chromosome 16. Nonetheless, they still seem to be linked in their cooperative functions, such as formation and physiology of flight muscle tissue. miR-133 has also been impli- cated in dopamine production (Yang et al., 2014), and high levels of dopamine were shown to coincide with rapid growth and com- partmentalization of the antennal lobe neuropil, suggesting a role in the developing brain of the honeybee (Kirchhof et al., 1999). Furthermore, dopamine-derivatives are substrates for oxidation by laccases (Andersen, 2010) that are involved in tanning of the developing adult cuticle (Elias-Neto et al., 2010). Members of the D. melanogaster let-7-C locus (a cluster containing the let- 7, miR-100, and miR-125 genes) are also found in the honeybee genome. In D. melanogaster they are expressed in neuromuscula- ture development of pupae and adults, and knockout flies showed disturbances in flight, reproduction and locomotion (Sokol et al., 2008). Moreover, ecdysteroid signaling was shown to be linked to the expression levels of the let-7-C cluster genes, as well as of miR-14 and miR-34 during insect development (for review see Kucherenko and Shcherbata, 2013). Many of the miRNAs affected by EcR knockdown in honeybees (let-7, miR-1, miR-9a, miR-12, miR-14, miR-34, miR-79, miR- 92b, miR-124, miR-184, miR-210, miR-219, miR-263a, miR-276, miR-279, miR-283, miR-305, miR-306, miR-316, miR-317) have previously been reported as putatively involved in the regulation of D. melanogaster immune genes, particularly those belonging to the JNK, Imd and Toll signaling pathways (Fullaondo and Lee, 2012). Accordingly, ecdysone and the ecdysone receptor complex (EcR/USP) are considered critical for innate cellular immunity (Flatt et al., 2008; Regan et al., 2013). Among these miRNAs, miR- 184 is highly and/or broadly expressed in a number of tissues and developmental stages of vertebrates (Wienholds and Plasterk, 2005) and invertebrates (Jagadeeswaran et al., 2010), including A. mellifera (Chen et al., 2010; Nunes et al., 2013b). Moreover, several studies reported a wide spectrum of roles for miR-184, such as germline differentiation, axis formation of the egg cham- ber, anteroposterior patterning and cellularization of the embryo, gastrulation and neuroectoderm formation, apoptosis, and pro- cesses involved in the development and differentiation of imaginal discs (head, wing, and eyes) (see Iovino et al., 2009; Li et al., 2011, and references therein). The ecdysone response of miR-184 seen here in pharate-adult honeybees is associated with a period of extensive tissue remodeling, suggesting that miR-184 may play a role in the differentiation of honeybee imaginal disc-derived structures and maintenance of their tissue identities. Interestingly, the EcR mRNA has predicted binding sites for miR-14 (data not shown), and our global gene expression assays revealed a down- regulation of this miRNA in bees silenced for EcR gene function. These results suggest that in A. mellifera, EcR gene expression is regulated in a loop-type mechanism involving miR-14, as already demonstrated for D. melanogaster (Varghese and Cohen, 2007; for a comprehensive review see Yamanaka et al., 2013). CONCLUDING REMARKS Our results suggest a differential use of EcR isoforms during the honeybee life-cycle stages. We could show that there is a generally positive EcR gene response to slight increases in ecdys- teroids, whereas high levels of these hormones are repressive. The EcR knockdown experiments revealed that the expression of Frontiers in Genetics | Evolutionary and Population Genetics December 2014 | Volume 5 | Article 445 | 16 http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development several hormone response-related genes (e.g., kr-h1) is contingent on a functional ecdysone/EcR receptor complex, thus provid- ing a possible link between JH and ecdysteroid action during preimaginal honeybee development. These knockdown experi- ments also hightlighted the relevance of a set of miRNAs involved in the regulation of immune response genes and in the gen- eral morphogenesis processes during pharate-adult development (e.g., miR-184 and let-7 locus genes). Within this framework and on the background of current knowledge on honeybee biology, our results highlight the relevance of the drop in the ecdysteroid pathway function for the appropriate timing in the expression of adult-specific genes, such as the Major Royal Jelly Protein (MRJP) family members. AUTHOR CONTRIBUTIONS Tathyana R. P. Mello, Aline C. Aleixo, Angel R. Barchuk, and Zilá L. P. Simões conceived the project; Tathyana R. P. Mello, Aline C. Aleixo, Daniel G. Pinheiro, Francis M. F. Nunes, Klaus Hartfelder, and Angel R. Barchuk performed the experiments; Tathyana R. P. Mello, Aline C. Aleixo, Daniel G. Pinheiro, Francis M. F. Nunes, Márcia M. G. Bitondi, Klaus Hartfelder, Angel R. Barchuk, and Zilá L. P. Simões analyzed and interpreted the data and drafted the MS. All authors approved the final version of the MS. FUNDING This study was funded by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), grants 2011/03171-5; 2008/1446- 4; 2008/10757-3. ACKNOWLEDGMENTS We thank Dr. Nilce M. M. Rossi, Roseli de Aquino P. Ferreira and Diana Gras for lab facility and support in Northern blotting assays, and Luiz Aguiar for technical assistance in the apiary. We also thank Felipe Martelli, Denyse C. Lago, Fabiano Abreu, and Camilla Valente Pires for help with the qPCR assays. SUPPLEMENTARY MATERIAL The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/journal/10.3389/fgene. 2014.00445/abstract REFERENCES Amdam, G. V., Simões, Z. L. P., Guidugli, K. R., Norberg, K., and Omholt, S. W. (2003). Disruption of vitellogenin gene function in adult honeybees by intra-abdominal injection of double-stranded RNA. BMC Biotechnol. 3:e1. doi: 10.1186/1472-6750-3-1 Andersen, S. O. (2010). Insect cuticular sclerotization: a review. Insect Biochem. Mol. Biol. 40, 166–178. doi: 10.1016/j.ibmb.2009.10.007 Ashburner, M., Chihara, C., Meltzer, P., and Richards, G. (1974). Temporal control of puffing activity in polytene chromosomes. Cold Spring Harb. Symp. Quant. Biol. 38, 655–662. doi: 10.1101/SQB.1974.038.01.070 Audic, S., and Claverie, J. M. (1997). The significance of digital gene expression profiles. Genome Res. 7, 986–995. Baldini, F., Gabrieli, P., South, A., Valim, C., Mancini, F., and Catteruccia, F. (2013). The interaction between a sexually transferred steroid hormone and a female protein regulates oogenesis in the malaria mosquito Anopheles gambiae. PLoS Biol. 11:e1001695. doi: 10.1371/journal.pbio.1001695 Barchuk, A. R., Bitondi, M. M. G., and Simões, Z. L. P. (2002). Effects of juvenile hormone and ecdysone on the timing of vitellogenin appearance in hemolymph of queen and worker pupae of Apis mellifera. J. Insect Sci. 2, 1–8. doi: 10.1673/031.002.0101 Barchuk, A. R., Cristino, A. S., Kucharski, R., Costa, L. F., Simões, Z. L. P., and Maleszka, R. (2007). Molecular determinants of caste differentiation in the highly eusocial honeybee Apis mellifera. BMC Dev. Biol. 7:70. doi: 10.1186/1471- 213X-7-70 Barchuk, A. R., Figueiredo, V. L. C., and Simões, Z. L. P. (2008). Downregulation of ultraspiracle gene expression delays pupal development in honeybees. J. Insect Physiol. 54, 1035–1040. doi: 10.1016/j.jinsphys.2008.04.006 Barchuk, A. R., Maleszka, R., and Simões, Z. L. P. (2004). Apis mellifera ultraspiracle: cDNA sequence and rapid up-regulation by juvenile hormone. Insect Mol. Biol. 13, 459–467. doi: 10.1111/j.0962-1075.2004.00506.x Beckstead, R. B., Lam, G., and Thummel, C. S. (2005). The genomic response to 20-hydroxyecdysone at the onset of Drosophila metamorphosis. Genome Biol. 6:R99. doi: 10.1186/gb-2005-6-12-r99 Bellés, X., and Santos, C. G. (2014). The MEKRE93 (Methoprene tolerant-Krüppel homolog 1-E93) pathway in the regulation of insect metamorphosis, and the homology of the pupal stage. Insect Biochem. Mol. Biol. 52, 60–68. doi: 10.1016/j.ibmb.2014.06.009 Boutz, P. L., Chawla, G., Stoilov, P., and Black, D. L. (2007). MicroRNAs regulate the expression of the alternative splicing factor nPTB during muscle development. Genes Dev. 21, 71–84. doi: 10.1101/gad.1500707 Buffalo, V. (2011). Scythe. Available online at: https://github.com/vsbuffalo/scythe Buttstedt, A., Moritz, R. F. A., and Erler, S. (2013). More than royal food - Major royal jelly protein genes in sexuals and workers of the honeybee Apis mellifera. Front. Zool. 10:e72. doi: 10.1186/1742-9994-10-72 Carstea, E. D., Morris, J. A., Coleman, K. G., Loftus, S. K., Zhang, D., Cummings, C., et al. (1997). Niemann-Pick C1 disease gene: homology to mediators of cholesterol homeostasis. Science 277, 228–231. doi: 10.1126/science.277.53 23.228 Charles, J.-P., Iwema, T., Epa, V. C., Takaki, K., Rynes, J., and Jindra, M. (2011). Ligand-binding properties of a juvenile hormone receptor, Methoprene-tolerant. Proc. Natl. Acad. Sci. U.S.A. 108, 21128–21133. doi: 10.1073/pnas.1116123109 Chawla, G., and Sokol, N. S. (2012). Hormonal activation of let-7-C microRNAs via EcR is required for adult Drosophila melanogaster morphology and function. Development 139, 1788–1797. doi: 10.1242/dev.077743 Chen, J.-F., Mandel, E. M., Thomson, J. M., Wu, Q., Callis, T. E., Hammond, S. M., et al. (2006). The role of microRNA-1 and microRNA-133 in skeletal muscle proliferation and differentiation. Nat. Genet. 38, 228–233. doi: 10.1038/ng1725 Chen, X., Yu, X., Cai, Y., Zheng, H., Yu, D., Liu, G., et al. (2010). Next- generation small RNA sequencing for microRNAs profiling in the honey bee Apis mellifera. Insect Mol. Biol. 19, 799–805. doi: 10.1111/j.1365-2583.2010. 01039.x Claudianos, C., Ranson, H., Johnson, R. M., Biswas, S., Schuler, M. A., Berenbaum, M. R., et al. (2006). A deficit of detoxification enzymes: pesticide sensitivity and environmental response in the honeybee. Insect Mol. Biol. 15, 615–636. doi: 10.1111/j.1365-2583.2006.00672.x Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610 Davis, M. B., and Li, T. (2013). Genomic analysis of the ecdysone steroid signal at metamorphosis onset using ecdysoneless and EcRnull Drosophila melanogaster mutants. Genes Genomics 35, 21–46. doi: 10.1007/s13258-013-0061-0 Dorn, A. (2000). “Arthropoda—insecta: embryology,” in Reproductive Biology of Invertebrates – Progress in Developmental Endocrinology, Vol. 10, Part B, ed A. Dorn (Chichester: Wiley & Sons), 72–116. Drapeau, M. D., Albert, S., Kucharski, R., Prusko, C., and Maleszka, R. (2006). Evolution of the yellow/major royal jelly protein family and the emer- gence of social behavior in honey bees. Genome Res. 16, 1385–1394. doi: 10.1101/gr.5012006 Elias-Neto, M., Soares, M. P. M., Simões, Z. L. P., Hartfelder, K., and Bitondi, M. M. G. (2010). Developmental characterization, function and regulation of a Laccase2 encoding gene in the honey bee, Apis mellifera (Hymenoptera, Apinae). Insect Biochem. Mol. Biol. 40, 241–251. doi: 10.1016/j.ibmb.2010. 02.004 Evans, R. M., and Mangelsdorf, D. J. (2014). Nuclear receptors, RXR, and the Big Bang. Cell 157, 255–266. doi: 10.1016/j.cell.2014.03.012 Fahrbach, S. E., Smagghe, G., and Velarde, R. A. (2012). Insect nuclear receptors. Annu. Rev. Entomol. 57, 83–106. doi: 10.1146/annurev-ento-120710-100607 www.frontiersin.org December 2014 | Volume 5 | Article 445 | 17 http://journal.frontiersin.org/journal/10.3389/fgene.2014.00445/abstract http://journal.frontiersin.org/journal/10.3389/fgene.2014.00445/abstract https://github.com/vsbuffalo/scythe http://www.frontiersin.org http://www.frontiersin.org/Evolutionary_and_Population_Genetics/archive Mello et al. Ecdysone receptor in honeybee development Feldlaufer, M. F., Herbert, E. W. J., and Svoboda, J. A. (1985). Makisterone A: the major ecdysteroid from the pupae of the honey bee, Apis mellifera. Insect Biochem. 15, 597–600. doi: 10.1016/0020-1790(85)90120-9 Feldlaufer, M. F., Herbert, E. W. Jr., Svoboda, J. A., and Thompson, M. J. (1986b). Biosynthesis of makisterone A and 20-hydroxyecdysone from labeled sterols by the honey bee, Apis mellifera. Arch. Insect Biochem. Physiol. 3, 415–421. doi: 10.1002/arch.940030502 Feldlaufer, M. F., Svoboda, J. A., and Herbert, E.W. Jr. (1986a). Makisterone A and 24-methylenecholesterol from the ovaries of the honey bee, Apis mellifera L. Experientia 42, 200–201. doi: 10.1007/BF01952468 Feldlaufer, M. F., Weirich, G. F., Lusby, W. R., and Svoboda, J. A. (1991). Makisterone C a 29-carbon ecdysteroid from developing embryos of the cot- ton stainer bug, Dysdercus fasciatus. Arch. Insect Biochem. Physiol. 18, 71–79. doi: 10.1002/arch.940180202 Flatt, T., Heyland, A., Rus, F., Porpiglia, E., Sherlock, C., Yamamoto, R., et al. (2008). Hormonal regulation of the humoral innate immune response in Drosophila melanogaster. J. Exp. Biol. 211, 2712–2724. doi: 10.1242/jeb.014878 Flatt, T., Tu, M. P., and Tatar, M. (2005). Hormonal pleiotropy and the juvenile hormone regulation of Drosophila development and life history. Bioessays 27, 999–1010. doi: 10.1002/bies.20290 Forêt, S., and Maleszka, R. (2006). Function and evolution of a gene family encoding odorant binding-like proteins in a social insect, the honey bee (Apis mellifera). Genome Res. 16, 1404–1413. doi: 10.1101/gr.5075706 Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688 Fullaondo, A., and Lee, S. Y. (2012). Identification of putative miRNA involved in Drosophila melanogaster immune response. Dev. Comp. Immunol. 36, 267–273. doi: 10.1016/j.dci.2011.03.034 Gáliková, M., Klepsatel, P., Senti, G., and Flatt, T. (2011). Steroid hormone reg- ulation of C. elegans and Drosophila aging and life history. Exp. Gerontol. 46, 141–147. doi: 10.1016/j.exger.2010.08.021 Geva, S., Hartfelder, K., and Bloch, G. (2005). Reproductive division of labor, dom- inance, and ecdysteroid levels in hemolymph and ovary of the bumble bee Bombus terrestris. J. Insect Physiol. 51, 811–823. doi: 10.1016/j.jinsphys.2005. 03.009 Hartfelder, K., Bitondi, M. M. G., Santana, W. C., and Simões, Z. L. P. (2002). Ecdysteroid titers and reproduction in queens and workers of the honey bee and of a stingless bee: loss of ecdysteroid function at increasing levels of sociality? Insect. Mol. Biol. 32, 211–216. doi: 10.1016/S0965-1748(01)00100-X Hartfelder, K., and Emlen, D. (2012). “Endocrine control of insect polyphenism,” in Insect Endocrinology, ed L. I. Gilbert (San Diego, CA: Academic Press), 464–522. Hartfelder, K., and Engels, W. (1998). Social insect polymorphism: hormonal regu- lation of plasticity in development and reproduction in the honeybee. Curr. Top. Dev. Biol. 40, 45–77. doi: 10.1016/S0070-2153(08)60364-6 Hill, R. J., Billas, I. M. L., Bonneton, F., Graham, L. D., and Lawrence, M. C. (2013). Ecdysone receptors: from the Ashburner model to structural biology. Annu. Rev. Entomol. 58, 251–271. doi: 10.1146/annurev-ento-120811-153610 Huang, X., Suyama, K., Buchnan, J., Zhu, A. J., and Scott, M. P. (2007). A Drosophila model of the Niemann-Pick type C lysosome storage disease: dnpc1a is required for molting and sterol homeostasis. Development 132, 5115–5124. doi: 10.1242/dev.02079 Humann, F. C., Tiberio, G. J., and Hartfelder, K. (2013). Sequence and expres- sion characteristics of long noncoding RNAs in honey bee caste development – potential novel regulators for transgressive ovary size. PLoS ONE 8:e78915. doi: 10.1371/journal.pone.0078915 Iga, M., and Kataoka, H. (2012). Recent studies on insect hormone metabolic path- ways mediated by Cytochrome P450 enzymes. Biol. Pharm. Bull. 35, 838–843. doi: 10.1248/bpb.35.838 Iida, K., Cox-Foster, D. L., Yang, X., Ko, W.-Y., and Cavener, D. R. (2007). Expansion and evolution of insect GMC oxidoreductases. BMC Evol. Biol. 7:75. doi: 10.1186/1471-2148-7-75 Iovino, N., Pane, A., and Gaul, U. (2009). miR-184 has multiple roles in Drosophila female germline development. Dev. Cell 17, 123–133. doi: 10.1016/j.devcel.2009. 06.008 Jagadeeswaran, G., Zheng, Y., Sumathipala, N., Jiang, H., Arrese, E. L., Soulages, J. L., et al. (2010). Deep sequencing of small RNA libraries reveals dynamic regula- tion of conserved and novel microRNAs and microRNA-stars during silkworm development. BMC Genomics 11:52. doi: 10.1186/1471-2164-11-52 Jindra, M., Malone, F., Hiruma, K., and Riddiford, L. M. (1996). Developmentl pro- files and edysteroid regulation of the mRNAs for two edysone receptor isofroms in the epidermis and wings of the tobacco hornworm, Manduca sexta. Dev. Biol. 180, 258–272. doi: 10.1006/dbio.1996.0299 Kamakura, M. (2011). Royalactin induces queen differentiation in honeybees. Nature 473, 478–483. doi: 10.1038/nature10093 Karim, F. D., and Thummel, C. S. (1992). Temporal coordination of regulatory gene expression by the steroid hormone ecdysone. EMBO J. 11, 4083–4093. Kato, S., Endoh, H., Masuhiro, Y., Kitamoto, T., and Uchiyama, S., Sasaki H., et al. (1995). Activation of the estrogen receptor through phosphorylation by mitogen-activated protein kinase. Science 270, 1491–1494. doi: 10.1126/sci- ence.270.5241.1491 Kelstrup, H., Hartfelder, K., Nascimento, F. S., and Riddiford, L. M. (2014). Reproductive status, endocrine physiology and chemical signaling in the Neotropical, swarm-founding eusocial wasp Polybia micans Ducke (Vespidae: Epiponini). J. Exp. Biol. 217, 2399–2410. doi: 10.1242/jeb.096750 King-Jones, K., and Thummel, C. S. (2005). Nuclear receptors, a perspective e from Drosophila. Nat. Rev. Genet. 6, 311–323. doi: 10.1038/nrg1581 Kirchhof, B. S., Homberg, U., and Mercer, A. R. (1999). Development of dopamine- immunoreactive neurons associated with the antennal lobes of the honey bee, Apis mellifera. J. Comp. Neurol. 411, 643–653. Kucherenko, M. M., and Shcherbata, H. R. (2013). Steroids as external tem- poral codes act via microRNAs and cooperate with cytokines in differential neurogenesis. Fly 7, 173–183. doi: 10.4161/fly.25241 Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357–359. doi: 10.1038/nmeth.1923 Li, P., Peng, J., Hu, J., Xu, Z., Xie, W., and Yuan, L. (2011). Localized expres- sion pattern of miR-184 in Drosophila. Mol. Biol. Rep. 38, 355–358. doi: 10.1007/s11033-010-0115-1 Lozano, J., and Belles, X. (2011). Conserved repressive function of Krüppel homolog 1 on insect metamorphosis in hemimetabolous and holometabolous species. Sci. Rep. 1:e163. doi: 10.1038/srep00163 Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. doi: 10.14806/ej.17.1.200 Michelette, E. R. F., and Soares, A. E. E. (1993). Characterization of preimagi- nal developmental stages of Africanized honey bee workers (Apis mellifera L.). Apidologie 24, 431–440. doi: 10.1051/apido:19930410 Mouillet, J. F., Delbecque, J. P., Quennedey, B., and Delachambre, J. (1997). Cloning of two putative ecdysteroid receptor isoforms from Tenebrio molitor and their developmental expression in the epidermis during metamorphosis. Eur. J. Biochem. 248, 856–863. doi: 10.1111/j.1432-1033.1997.00856.x Mulder, N., and Apweiler, R. (2007). InterPro and InterProScan: tools for protein sequence classification and comparison. Methods Mol. Biol. 396, 59–70. doi: 10.1007/978-1-59745-515-2_5 Nunes, F. M. F., Aleixo, A., Barchuk, A. R., Bomtorin, A. D., Grozinger, C., and Simões, Z. L. P. (2013a). Non-target effects of Green Fluorescent Protein (GFP)-derived double-stranded RNA (dsRNA-GFP) used in honey bee RNA interference (RNAi) assays. Insects 4, 90–103. doi: 10.3390/insects4010090 Nunes, F. M. F., Ihle, K. E., Mutti, N. S., Simões, Z. L. P., and Amdam, G. V. (2013b). The gene vitellogenin affects microRNA regulation in honey bee (Apis mellifera) fat body and brain. J. Exp. Biol. 216, 3724–3732. doi: 10.1242/jeb. 089243 O’Reilly, D. R., and Miller, L. R. (1989). A baculovirus blocks insect molting by producing ecdysteroid UDP-glucosyl transferase. Science 245, 1110–1112. Ono, H. (2014). Ecdysone differentially regulates metamorphic timing rela- tive to20-hydroxyecdysone by antagonizing juvenile hormone in Drosophila melanogaster. Dev. Biol. 391, 32–42. doi: 10.1016/j.ydbio.2014.04.004 Petryk, A., Warren, J. T., Marques, G., Jarcho, M. P., Gilbert, L. I., Kahler, J., et al. (2003). Shade is the p450 enzyme that mediates the hydroxylation of ecdysone to the steroid insect molting hormone 20-hydroxyecdysone. Proc. Natl. Acad. Sci. U.S.A. 100, 13773–13778. doi: 10.1073/pnas.2336088100 Pinto, L. Z., Hartfelder, K., Bitondi, M. M. G., and Simões, Z. L. P. (2002). Ecdysteroid titers in pupae of highly social bees relate to distinct mode