1 Introduction

Royal jelly (RJ) is the food for all young larvae in the first 3 days after hatching, and is also the exclusive food for the queen throughout her lifetime (Winston 1987), playing a vital role in larval development. It has been proved to possess a wide range of pharmacological effects, including antimicrobial, antitumor, and anti-inflammatory activities (Nikokar and Shirzad 2014; Arzi et al. 2015). As RJ is helpful for human health, it has been widely used as supplementary food and medicine. To meet the increasingly demand for RJ, Chinese scientists and beekeepers have selected and bred high-RJ-producing strains of bees generating as much as ten times larger amounts of RJ than normal Italian bees for nearly four decades (Altaye et al. 2019).

The hypoharyngeal glands (HGs) are the major organ for the synthesis and secretion of RJ in honeybees (Painter and Biesele 1966; Cruz-Landim and Costa 1998); accordingly, extensive research has been carried out in terms of their morphology and physiology. Morphologically, HG is a paired long tuberous organ connected to a series of acini that reach peak size on 6–12 days and begin to shrink after day 15 during the summertime (Mao et al. 2009; Qi et al. 2015). Ultrastructure of the HG cells showed that the rough endoplasmic reticulum and secretory vesicles reach the maximum during nursing phase and decrease in forager bees (Knecht and Kaatz 1990; Jeroen and Johan 2005). Consisting of eight secretory cells, the size of acini is positively correlated with the RJ secretory activity (Albert et al. 2014). With regard to the physiology, HGs change along with the age-dependent role of workers. Specifically, major royal jelly proteins (MRJPs) are mainly synthesized by nurse bees, while carbohydrate-metabolizing enzymes (including α-glucosidase, α-amylase, and glucose oxidase) tend to be produced by foragers to process nectar to honey (Kubo et al. 1996; Ohashi et al. 1997; Ohashi et al. 1999). The physiological state of workers’ HGs is cooperatively regulated by ecdysone/juvenile hormone-signaling and diets (DeGrandi-Hoffman, Chen et al. 2010; Ueno et al. 2015). With the development of omics technologies, transcriptome and proteomics of HGs at different developmental stages in Apis mellifera were studied comprehensively (Mao et al. 2009; Liu et al. 2013; Ji et al. 2014; Qi et al. 2015). As high-RJ-producing bees possess extremely stronger ability of RJ production than Italian bees, it is a good strain for studying the molecular basis of RJ secretion. Proteomic analysis of HGs revealed that high-RJ-producing bees tend to have a higher level of proteins associated with protein synthesis and energy metabolism for enhancing RJ production compared with Italian bees (Li et al. 2010; Hu et al. 2019). In worker honeybees, the age-dependence of the role switch is flexible and largely depends on colony demands (Hrassnigg and Crailsheim 1998; Ohashi et al. 2000).

Previous studies have focused on the functions of HGs in honeybees with different ages; however, the influence of ages on the production of RJ is usually ignored. Thus, in the present study, to clarify the molecular basis of RJ production, acinus size and total protein content of HGs from age-matched nurses and foragers were initially determined, and then genes associated with RJ production were identified using RNA-seq. This study will provide valuable information for not only better understanding of the molecular mechanism for RJ production but also a potentially new avenue to breed a high-royal-jelly-producing bee strain.

2 Materials and methods

2.1 Establishment of experimental colonies and sample collection

High-RJ-producing honeybees, Fengqiang No. 1 Italian honey bees (A. mellifera lingustica), were maintained at the apiaries in the College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University. To collect age-matched samples of nurse and forager bees, experimental colonies were built in the apiaries as described previously (McQuillan et al. 2012) with minor modifications. In brief, capped frames were sourced from five queen-right colonies and housed in an incubator at 34 °C. Every morning, newly emerged bees (NEBs, < 24 h) were marked on their thoraxes with paint and put in a new hive. Markings continued for 5 days to make sure that the new hive could maintain its colony. NEBs were marked with different colors each day so that we could easily distinguish the age of each worker and collect age-matched samples (10-day-old and 21-day-old) of nurse and forager bees. A total of > 12,500 NEBs from five colonies (about 2500 NEBs each) were labeled and placed into the new hive containing a young mated queen, a honey comb, a pollen comb, a larvae comb, as well as an empty comb for egg laying. Three-day-old worker bees were collected from the new hive when marked bees were on the third day after emergence. Nurses and forager bees were captured according to a previous method (Liu et al. 2011). On the tenth day after emergence, 10-day-old nurses were collected by opening the nest and identifying individuals keeping their heads and thoraxes inside the brood cells for at least 10 s; 10-day-old foragers, identified by the presence of pollen loads on their hind legs, were collected at the entrance of new hive using soft forceps. Similarly, 21-day-old nurses and foragers were collected with the same method described above when marked bees were 21 days post-eclosion. The bees were anesthetized on ice, and HGs were dissected from bee heads using fine dissection scissors and a dissecting microscope.

2.2 Microscopic observation and calculation of HG acinus areas

After dissecting, HGs were stained with Giemsa solution following methods described previously (Corby-Harris and Snyder 2018). Briefly, each HG was cut into five pieces with the micro-scissors and put on the microscope slides. They were stained with Giemsa solution for 5 min and then immersed in phosphate-buffered saline to acquire an image using optical microscopy. The morphologies of HG acinus were observed and photographed at 10× magnification in a Leica DM2000 microscopes (Leica, Germany). A photo of the HG acinus was imported into Photoshop CS6, and each acinus area was calculated using pixel-ratio method according to previous study (Peng et al. 2018). To improve the accuracy of the measurement, five bees were measured in each group and ten acini were measured for each bee. The average size of acinus for each bee in each group was calculated on the basis of 50 acini from five bees.

2.3 Determination of total protein content of HGs

Ten HGs from each group were dissected and rinsed three times with phosphate buffered saline. The samples were homogenized in 100 μL of total protein extraction solution containing a mixture of proteinase inhibitors (TransGen Biotech, China) and kept for 10 min at 4 °C. The protein fraction was centrifuged at 12,000 ×g for 15 min at 4 °C to remove precipitate. Protein concentrations in the supernatant were measured using TaKaRa BCA Protein Assay Kit according to the instructions.

2.4 Library construction and transcriptome sequencing

Total RNA (3 μg) of HGs from 3-day-old workers, 10-day-old nurses, 10-day-old foragers, 21-day-old nurses, and 21-day-old foragers was extracted using TRIzol reagent (Invitrogen, Burlington, ON, Canada) according to the manufacturer’s protocol. Library construction and sequencing of 15 RNA samples (three biological replicates per stage) were performed as described previously with minor modification (Nie et al. 2018). Sequencing of the 15 HG samples was conducted using Illumina HiSeq X Ten (Illumina, Inc. San Diego, USA), and 150 bp paired end reads were generated. The raw data of RNA-seq presented in this article have been deposited to NCBI Short Read Archive and can be accessed from corresponding BioProject via the submission portal:

https://submit.ncbi.nlm.nih.gov/subs/bioproject/SUB7391751/overview.

2.5 RNA-Seq analysis

Clean reads were mapped to the A. mellifera genome using Hisat2 (V2.0.5) (Daehwan et al. 2015). FeatureCounts (V1.5.0-p3) was employed to count the number of reads mapped to each gene and gene expression levels were calculated using FPKM. For differential expression analysis, DESeq2 R package (1.16.1) was employed to determine differential expression of two groups using a model based on the negative binomial distribution (Love et al. 2014). The false discovery rate control was used to identify the adjusted p values using Benjamini and Hochberg’s approach. An adjusted p value < 0.05 was set as the criteria for screening differentially expressed genes (DEGs) by DEGSeq2. GO/KEGG enrichment analysis was carried out using the clusterProfiler R package with default parameters (Kanehisa and Goto 2000; Tarazona et al. 2011). A corrected p value was set as a threshold value to test the statistically significant enrichment. Weighted correlation network analysis (WGCNA) was performed using WGCNA package as described previously (Langfelder and Horvath 2008).

2.6 RNA-seq data validation using qRT-PCR

qRT-PCR was performed following our previous method (Nie et al. 2018). Briefly, cDNA of each sample was synthesized using 1 μg of HGs total RNA and PrimeScript RT reagent Kit (RR037A, Takara). The A. mellifera housekeeping gene actin (GenBank accession number NM_001185146.1) was used as the internal control for RNA normalization. The relative gene expression was calculated relative to the actin using the 2−ΔΔCt method (Livak and Schmittgen 2001).

2.7 Statistical analysis

Statistical differences were evaluated with an analysis of variance (ANOVA) with Tukey’s HSD using SPSS20.0 software. Asterisks indicate significant differences (*, p < 0.05; **, p < 0.01; ***, p < 0.001).

3 Results

3.1 Comparison of acinus area and total protein content in HGs of age-matched nurse and forager workers

Usually, the initial increase and subsequent decrease in the size of acini are associated with the change of age-dependent roles from nurse bees to foragers. To investigate the morphologic change of HGs from age-matched workers, HGs from 10-/21-day-old nurse and forager bees were stained with Giemsa solution (Figure 1A–D). Using pixel-ratio method, HG acinus areas of 10-day-old nurses possessed significantly larger acinus area (0.004705 ± 0.001162) than those of the 10-day-old foragers (0.002212 ± 0.000231) (p < 0.05). Similarly, 21-day-old nurses possessed significantly larger acinus area (0.007712 ± 0.001083) than the foragers (0.00344 ± 0.000883) at the same age. Meanwhile, we observed that the acinus size of 21-day-old nurses was also significantly larger than that of 10-day-old nurse bees (p < 0.05). As the RJ protein synthesis in HGs depends on the age-dependent role change of worker bees, total proteins in HGs from age-matched workers were determined using BCA Protein Assay Kit. It showed that total protein contents of each HG sample from both 10- and 21-day-old nurses were extremely significantly higher than those in corresponding foragers (Figure 2A). Since major royal jelly protein 1 (MRJP1) is the most abundant water-soluble protein in RJ, it can be used to determine the secretory activity of HGs (Hu et al. 2019). To investigate whether HGs of nurses has higher secretory activity than those of the same aged-matched foragers, the expression of MRJP1 gene was detected using qRT-PCR (Figure 2B). We found that the MRJP1 gene was highly expressed in the HGs of 10-day-old and 21-day-old nurses, which was significantly higher than the expression in corresponding age-matched foragers. The expression profile of MRJP1 was almost comparable to the pattern of total protein content from HGs, suggesting that HGs from nurses have a higher total protein content and secretory activity.

Figure 1.
figure 1

Stained hypopharyngeal gland (HGs) from age-matched nurse and forager workers. A, B HGs from 10-day-old nurses and forager bees, respectively. C, D HGs from 21-day-old nurses and forager bees, respectively. E) Size of HG acini age-matched nurse and forager bees. Asterisks indicate statistically significant differences between the mean acinus area of each group.

Figure 2.
figure 2

Total protein content and relative expression of MRJP1 in the HGs from age-matched nurses and forager workers. A Total protein content of each HG from 10-day-old/21-day-old nurses and forager bees. B Relative expression of MRJP1 gene in HGs of 10-/21-day-old nurses and forager bees using qRT-PCR. Asterisks indicate the statistically significant differences between each group.

3.2 Validation of data reliability

Ten-pooled total RNA samples were collected from the HGs of 3-day-old workers, 10-day-old nurses, 10-day-old foragers, 21-day-old nurses, and 21-day-old foragers for paired-end cDNA library construction and sequencing. The statistics for all samples are summarized in Table S2. Previous studies revealed that the expression of AmGR10 and Ambuffy was remarkably up-regulated in HGs of nurses than that of foragers (Paerhati et al. 2015; Ueno et al. 2015). The RNA-seq results showed that the two genes were significantly up-regulated in HGs of 10-day-old nurses compared with those of 10- and 21-day-old foragers, and also displayed significant upregulation in HGs of 21-day-old nurses relative to those of 21-day-old foragers (Table S3). This expression profile was similar to previous studies (Paerhati et al. 2015; Ueno et al. 2015). To further verify the RNA-seq results, six DEGs were selected and validated using qRT-PCR (Fig. S1). The results showed that the expression patterns of the six DEGs measured by qRT-PCR were consistent with those of the DEGs expression profile, indicating the reliability of RNA-seq results.

3.3 Profile of MRJPs in HGs from age-matched nurses and foragers

Approximately 90% of total RJ proteins are major royal jelly proteins (MRJPs) (Yu et al. 2010). In A. mellifera, MRJPs consist of ten members, including MRJP1-9 and non-transcribed pseudogene MRJP-ψ/10 (Drapeau et al. 2006; Buttstedt et al. 2017). As HGs are important exocrine glands for RJ production, the expression patterns of the MRJPs from age-matched nurses and forgers were investigated by RNA-seq (Figure 3 and Table S4). MRJPs, except for MRJP6, displayed similar expression patterns that they were much highly expressed in 10-day-old nurses than those in 3-day-old workers, 10-day-old foragers, and 21-day-old foragers; they were also expressed at higher levels in the HGs of 21-day-old nurses than those of foragers at the same age, suggesting that they might play a crucial role in the synthesis of royal jelly protein. Among MRJPs, MRJP1 exhibited the highest gene expression, with FPKM values of 9242.4, 480,014.2, 10,046.1, 659,679.4, and 65,100.2 in the HGs of 3-day-old workers, 10-day-old nurses, 10-day-old foragers, 21-day-old nurses, and 21-day-old foragers, respectively. MRJP8 and MRJP9 showed a relatively lower expression level (FPKM < 20), although they still had a similar gene expression profiles to others. In contrast, the expression of MRJP6 was higher in the HGs of 10-/21-day-old foragers than those of nurses with corresponding age, and it was also significantly increased in the HGs of 21-day-old foragers than those of 10- and 21-day-old nurses. These expression profiles were consistent with previous studies (Dobritzsch et al. 2019).

Figure 3.
figure 3

Expression profiles of MRJPs in the HGs of 3-day-old workers, 10-day-old nurses/foragers, 21-day-old nurses/foragers by FPKM value. MRJPs: major royal jelly proteins; HG: hypopharyngeal glands; N: nurses; F: foragers.

3.4 DEGs associated with synthesis and secretion of RJ proteins in HGs

To obtain DEGs associated with RJ production, rather than age, we screened DEGs according to the following criteria: the gene expression was up-regulated in the HGs of 10-day-old nurses compared with both 3-day-old workers and 21-day-old foragers; gene expression was up-regulated in HGs of nurses (10−/21-day-old) than foragers of the same age; there was no significant difference between 10-day-old nurses and 21-day-old nurses. A total of 807 common genes were identified in the HGs between 10-day-old nurses and 3-day-old workers, 10-day-old nurses and 21-day-old foragers, 10-day-old nurses and 10-day-old foragers, and 21-day-old nurses and 21-day-old foragers, without significant difference between 10-day-old nurses and 21-day-old nurses (Figure 4A and Table S3). Hierarchical clustering analysis of 807 DEGs showed five differential expression patterns (Figure 4B). Among these genes, only 510 genes met the above criteria for screening DEGs associated with RJ production rather than age (Table I and Table S5).

Figure 4.
figure 4

Venn diagram and hierarchical clustering analysis of differentially expressed genes. A Venn diagram of differentially expressed genes in the HGs of five honeybee samples of defined age. B) Hierarchical clustering analysis of 807 common differentially expressed genes in the HGs of 3-day-old workers, 10-day-old nurses, 10-day-foragers, 21-day-nurses, and 21-day-foragers. HG: hypopharyngeal glands; N: nurses; F: foragers; 10d: 10-day old; 21d: 21-day-old.

Table 1 Significance of 807 DEGs in HGs of five honeybee samples of defined age

3.5 GO and KEGG enrichment analysis of up-regulated DEGs in HGs of nurses

The 510 genes were subjected to GO analysis and mapped to reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG). Based on the category of GO term/KEGG, these genes were mainly related to translation (157 genes, Table S6), oxidoreductase activity (26 genes, Table S7), ribosome biogenesis (22 genes, Table S8), transmembrane transport (21 genes, Table S9), transcription (20 genes, Table S10), energy metabolism (13 genes, Table S11), amino acid metabolism (11 genes, Table S12), fatty acid metabolism (7 genes, Table S13), and DNA replication (7 genes, Table S14) (Figure 5). Under the term of translation, they were associated with ribosome (76 genes), amino acid activation (18 genes), translation factor activity (19 genes), protein modification process (26 genes), protein processing in endoplasmic reticulum (7 genes), as well as protein export (5 genes). A total of 85 different KEGG pathways were identified (Table S15). The top 20 KEGG pathways were showed in Figure 6, including ribosome, aminoacyl-tRNA biosynthesis, ribosome biogenesis in eukaryotes, citrate cycle (TCA cycle), fatty acid elongation, protein export, N-Glycan biosynthesis, carbon metabolism, folate biosynthesis, and cysteine and methionine metabolism. Only the top two KEGG pathways, ribosome (71 genes) and aminoacyl-tRNA biosynthesis (17 genes), were notably enriched in DEGs associated with RJ secretion. Many up-regulated genes in the HGs of nurses were involved in GO terms/KEGG related to protein biosynthesis, suggesting that the biological effects of these genes were activated on some occasions to enhance the activity of HGs, thereby promoting RJ secretion.

Figure 5.
figure 5

Enrichment analysis of upregulated DEGs in HGs of nurses based on GO terms and KEGG.

Figure 6.
figure 6

Information of the top 20 KEGG pathway.

3.6 Interaction network analysis of MRJP1, MRJP4 and MRJP5

Interaction network of 510 up-regulated DEGs was analyzed using WGCNA. Among these genes, as the major components of royal jelly, MRJP1, MRJP4, and MRJP5 also met the screening criteria mentioned above with extremely high expression level in HGs of nurses compared with that in age-matched foragers. To identify the interaction network of DEGs, the top 150 up-regulated genes interacting with MRJP1, MRJP4, and MRJP5, respectively, were selected based on the weight value. Gene network analysis showed that MRJP1, MRJP4, and MRJP5 interacted with many up-regulated DEGs associated with ribosome, 37 of which interacted with the three genes simultaneously (Figure 7).

Figure 7.
figure 7

Gene interaction network analysis of MRJP1, MRJP4, and MRJP5 with up-regulated differentially expressed genes. The gene interaction networks were constructed using Cytoscape. Cycle nodes represent genes. DEGs with yellow background belong to ribosome.

4 Discussion

This study compared the HG acinus area and total protein content, in parallel with a systematic analysis of HGs from age-matched nurse and forager workers using RNA-Seq. It was found that 510 genes, mainly enriched in translation, transcription, DNA replication, and energy metabolism, were associated with synthesis and secretion of RJ proteins in HGs. It is the first comprehensive analysis of HGs from nurses and foragers at the same ages (10- and 21-day-old), which can eliminate the interference of age-dependent nursing. The results provided an integrated perspective on the genes involved in the protein synthesis and secretion, facilitating better understanding of the molecular basis of RJ production in HGs of nurses.

4.1 Up-regulated DEGs were associated with DNA replication in HGs of nurses

The process of DNA replication is strictly modulated to ensure the fidelity of genetic information to progenies. In this study, seven genes associated with DNA replication were up-regulated in HGs of nurses compared with those of other groups. Particularly, DNA ligase-1 like gene, homologous gene of Saccharomyces cerevisiae and mammalian used as a key ligase for DNA replication (Johnston and NASMYTH 1978; Petrini et al. 1995), was highly expressed in the HGs of nurses. The expression was about six times higher in the HGs of 10-day-old nurses than those of same age foragers and around four times higher in the HGs of 21-day-old nurses than those of same age foragers (Table S3 and Table S14). It indicated that HGs might possess stronger ability of DNA replication to enhance DNA synthesis.

4.2 Up-regulated DEGs were associated with transcription process in HGs of nurses

Transcription was divided into three phases: initiation, elongation, and termination. During this process, RNA polymerases and transcriptional factors are necessary for the transcription of DNA into RNA (Griesenbeck et al. 2017). We found that 20 up-regulated genes related to the transcription process in the HGs of nurses included RNA polymerase, initiation factors, elongation factors, and the other four transcriptional factors (Table S10). These genes might be involved in RNA synthesis, improving the gene expression during RJ secretion.

4.3 Expression of RJ production–related genes was enhanced in the HGs of nurses

Of the 510 up-regulated genes, 157 genes were associated with translation (Figure 5). In terms of KEGG enrichment analysis, only ribosome (71 genes) and aminoacyl-tRNA biosynthesis (17 genes) were markedly enriched in the HGs of nurses. It was consistent with the previous analysis that the proteins in the HGs of nurses were dramatically enriched in ribosome, aminoacyl-tRNA biosynthesis, protein processing in endoplasmic reticulum, and protein export (Qi et al. 2015). Notably, 76 ribosome genes and 18 aminoacy-tRNA synthetase genes were up-regulated in HGs of nurses. Recently, it was reported that initiation factor eIF5B was vital in the transition from translation initiation to elongation (Wang et al. 2019). In this work, we found that 19 genes were enriched in translation factor activity, including initiation factor, elongation factor, and release factor, which might enhance the procedure of RJ protein synthesis—initiation, elongation, and termination. Additionally, 26 genes were associated with protein modification, such as glycosylation, phosphorylation, ubiquitination, deubiquitination, and protein folding (Table S6). Seven hundred and twenty phosphoproteins were identified in the HGs of nurses (Qi et al. 2015). As HGs are the major organ for RJ protein secretion, the identification of glycoproteins and phosphoproteins from Apis mellifera ligustica and Apis cerana cerana RJ implied the stronger RJ-secretory activity of the HGs of nurses boosted RJ secretion via a variety of protein modification (Han et al. 2014; Feng et al. 2015).

MRJPs account for ~ 90% of total RJ proteins secreted by HGs (Drapeau et al. 2006). In this study, the expression of MRJPs (except for MRJP6) in the HGs of nurses (10- and 21-day-old) was higher than those of age-matched foragers. MRJP1-MRJP7 (except for MRJP6) were highly expressed with FPKM > 20,000 in the HGs of nurses which constituted the main components of total RJ proteins. Among them, three genes (MRJP1, MRJP4, and MRJP5) were extremely up-regulated in HGs of nurses compared with those in age-matched foragers, and interaction network analysis showed they interacted with many up-regulated DEGs involved in ribosome. Surprisingly, major royal jelly protein 1 (551813) was detected and also met the criteria, with FPKM of 9.2, 38,593.7, 232.2, 28,490.6, and 1063 in the HGs of 3-day-old workers, 10-day-old nurses, 10-day-old foragers, 21-day-old nurses, and 21-day-old foragers, respectively. However, this gene was not in current annotation release (Amel_HAv3) of NCBI. To validate its authenticity, special primers were designed to investigate the gene expression in our samples by qRT-PCR (Fig. S1 and S2). Its expression profile was consistent with transcriptome sequencing, which was up-regulated significantly in the HGs of nurses relative to those of foragers at the same age, indicating that it might be a key component of RJ protein derived from HGs of nurses.

4.4 Up-regulated DEGs were associated with energy metabolism in HGs of nurses

HGs of high-RJ-producing nurse bees manifest stronger energy supplement than those of Italian bees (Hu et al. 2019). In cells and organs, citric acid cycle (TCA cycle) and glycolysis are the key process for energy generation, in which a series of enzymes are involved (Keighren et al. 2016; Hallan et al. 2017). In this study, 13 genes participating in TCA cycle and glycolysis were significantly up-regulated in HGs of nurses compared with those of age-matched foragers (Table S11). In TCA cycle, pyruvate dehydrogenase (551103), citrate synthase (410059), succinyl-CoA synthetase (102656009), succinate dehydrogenase (408734, 550667, and 100578223), fumarate hydratase (724321), and phosphoenolpyruvate carboxykinase (412843) were up-regulated in HGs of nurses to create a large number of ATP via aerobic oxidation of glucose. Under anaerobic conditions, glycolysis is the major route for the catabolism of glucose, fructose, and galactose to generate energy. Four genes, including glucose-6-phosphate isomerase (551154), triosephosphate isomerase (726117), phosphoglycerate mutase (552736), and enolase (552678) associated with glycolysis, were up-regulated in HGs of nurses, suggesting that it generated greater amount of energy to facilitate protein synthesis and other physiological processes.

As mentioned above, in combination with our previous study that odorant binding protein 17 gene was up-regulated in the antennae of nurses compared with those in foragers (Nie et al. 2018), it was speculated that nurses could perceive volatile larval pheromones (such as β-ocimene and allo-ocimene) with chemosensory organs, transmit these signals to brain, process it in the higher processing centers of the brain, induce stronger DNA replication and transcription in HGs, and ultimately drive the process of RJ protein synthesis to feed small larvae and queen.

5 Conclusions

As task performance of workers is correlated to age of bees in a normal colony, we constructed an artificial colony to obtain the even-aged 10-/21-day-old nurses and foragers to exclude the effects of age on the morphology, total protein content, as well as candidate genes associated with RJ synthesis and secretion in the HGs of nurses. As a consequence, we found that HGs of 10-/21-day-old nurses possessed a larger acinus area and higher total protein content than those of foragers with the same age. Furthermore, RNA-seq analysis revealed that 510 genes associated with synthesis and secretion of RJ were mainly enriched in translation, transcription, DNA replication, and energy metabolism, suggesting that HGs of nurses might employ diverse strategies described above to stimulate highly efficient RJ protein synthesis. This is the first identification of genes that are possibly involved in the RJ production in the HGs of age-matched nurses and foragers, providing a new perspective into the molecular basis underlying RJ secretion in honey bees.