The rise and fall of major royal jelly proteins during a honeybee (Apis mellifera) workers' life

Abstract The genome of the western honeybee (Apis mellifera) harbors nine transcribed major royal jelly protein genes (mrjp1‐9) which originate from a single‐copy precursor via gene duplication. The first MRJP was identified in royal jelly, a secretion of the bees' hypopharyngeal glands that is used by young worker bees, called nurses, to feed developing larvae. Thus, MRJPs are frequently assumed to mainly have functions for developing bee larvae and to be expressed in the food glands of nurse bees. In‐depth knowledge on caste‐ and age‐specific role and abundance of MRJPs is missing. We here show, using combined quantitative real‐time PCR with quantitative mass spectrometry, that expression and protein amount of mrjp1‐5 and mrjp7 show an age‐dependent pattern in worker's hypopharyngeal glands as well as in brains, albeit lower relative abundance in brains than in glands. Expression increases after hatching until the nurse bee period and is followed by a decrease in older workers that forage for plant products. Mrjp6 expression deviates considerably from the expression profiles of the other mrjps, does not significantly vary in the brain, and shows its highest expression in the hypopharyngeal glands during the forager period. Furthermore, it is the only mrjp of which transcript abundance does not correlate with protein amount. Mrjp8 and mrjp9 show, compared to the other mrjps, a very low expression in both tissues. Albeit mrjp8 mRNA was detected via qPCR, the protein was not quantified in any of the tissues. Due to the occurrence of MRJP8 and MRJP9 in other body parts of the bees, for example, the venom gland, they might not have a hypopharyngeal gland‐ or brain‐specific function but rather functions in other tissues. Thus, mrjp1‐7 but not mrjp8 and mrjp9 might be involved in the regulation of phenotypic plasticity and age polyethism in worker honeybees.

All MRJPs can be detected in food jelly, with MRJP1-3 and MRJP5 accounting for 82%-90% of total food jelly proteins (Schmitzová et al., 1998;Schönleben, Sickmann, Mueller, & Reinders, 2007;Zhang et al., 2014), and have undoubtedly a nutritional function. However, functions go far beyond that: An oligomeric form of MRJP1 (oligoMRJP1) builds a pH-dependent fibrillary network  in complex with apisimin, a serine-valine-rich small protein that is another proteinaceous component of royal jelly (Bíliková et al., 2002). This fibrillary network confers the needed viscosity to royal jelly to prevent queen larvae falling out of their vertically oriented queen cells Kurth, Kretschmar, & Buttstedt, 2019). In addition, the complex of oligoMRJP1/apisimin binds 24-methylenecholesterol and provides the developing larvae with essential sterols (Tian et al., 2018). MRJP3 binds and stabilizes RNA in royal jelly and is thought to share RNA among individuals (Maori et al., 2019).
When worker bees were fed with labeled RNA, this RNA was found bound to MRJP3 in the food jelly produced by these worker bees (Maori et al., 2019). It is thought that this transmission of RNA from workers to larvae could drive social immunity against pathogens (Maori et al., 2019). In addition, oligoMRJP1/apisimin, MRJP2, and MRJP4 have antibacterial activity in vitro (Bíliková, Wu, & Šimúth, 2001;Kim et al., 2019;Vezeteu, Bobiş, Moritz, & Buttstedt, 2017).
Besides their expression in food-producing HGs, mrjp mRNA and the resulting proteins have been detected in a variety of tissues, for example, antennae, brain, nerve chord, hemolymph, and the Malpighian tubule system, not only in worker bees but also in drones and queens (Buttstedt, Moritz, & Erler, 2014;Chan et al., 2013;Whitfield et al., 2002). However, the focus of the expression of mrjp1-7 was clearly assigned to the heads of worker bees (Buttstedt, Moritz, & Erler, 2013). Expression of these genes is not only upregulated in food jelly-producing nurse bees but also in forager bees when compared to caged worker bees outside of the hive context (Buttstedt et al., 2013). Surprisingly, for mrjp1, mrjp2, mrjp5, and mrjp7 expression in nurse bees was not significantly higher than in forager bees (Buttstedt et al., 2013) albeit nurse bees are feeding larvae whereas foragers do not. This partially contrasts earlier studies reporting on higher expression of mrjp1, mrjp3, and mrjp4 in nurse bee heads (Klaudiny, Kulifajová, Crailsheim, & Šimúth, 1994;Ohashi, Natori, & Kubo, 1997) and MRJP1-3 to only be detectable in HGs of nurse but not forager bees (Kubo et al., 1996). Apart from their occurrence in the HGs, transcripts of all mrjps were found in the F I G U R E 1 Details of worker honeybees (Apis mellifera) (a) Worker bees storing pollen and honey (upper left) on a frame next to nurse bees taking care of the brood (lower right). (b) Upon opening the head capsule of a bee, the hypopharyngeal glands (arrow) are visible as oval acini attached to a collecting duct. (c) After removal of the hypopharyngeal glands, the brain (arrow) can be seen below, surrounded by air sacs (c) (b) (a) honeybee brain in a study analyzing brain expressed sequence tag libraries (Whitfield et al., 2002) and the question arises whether the observed expression in forager heads might be caused by a shift of the expression from the HGs to the brain while the bees age.
To obtain a deeper understanding of the potential involvement of MRJPs in the appearance of phenotypic plasticity and age polyethism in honeybees, we combined quantitative real-time PCR (qPCR) with quantitative mass spectrometry, to elucidate the intensity and timing of mrjp transcription and translation in both HGs ( Figure 1b) and brains ( Figure 1c) throughout a worker honeybee's life from hatching to the forager stage using a fully active bee hive.

| Honeybee samples
Honeybees (Apis mellifera) were sampled in May and June 2016 from a queen-right brood-rearing colony located in Halle (Saale), Germany (latitude: 51.5046; longitude: 11.9493). To raise age-matched worker bees, a brood frame containing pupae with dark eyes was removed from the hive and incubated at 34°C and 60% relative humidity. A total of 600 freshly hatched bees were paint-marked (Edding 751 gloss paint markers, Edding) on their thoraces and returned to the hive. After 0 (directly after hatching), 4, 8, 12 (nurse bee period), 16, 20 (transition phase), and 24 days (foragers), ten bees per day were freeze-killed in liquid nitrogen and stored at −80°C until further processing. Bees were only sampled according to age, and it is not known whether, for example, 24-day-old bees foraged for pollen, nectar, water, or resin.

| Gene expression
Honeybee HGs (Figure 1b) and brains ( Figure 1c) of ten bees per time point were dissected, washed in insect saline (Carreck et al., 2013), and immediately placed into 200 µl buffer RA1 supplemented with β-mercaptoethanol (NucleoSpin ® RNA Kit, Macherey-Nagel). RNA was further extracted according to the manufacturer's protocol that included a DNase digestion step. The flow-through after binding of the RNA to the NucleoSpin ® RNA columns was retained for subsequent protein isolation (see Section 2.32.3). Quantity of total RNA was photometrically determined with a NanoDrop 1,000 (Thermo Fisher Scientific), and total RNA per HG pair or brain per bee was calculated ( Figure S1).
Consequently, only pros26 (SD: 0.77) was satisfactory as reference gene from the initial testing group. In addition, mrjp8 (SD: 0.97) showed a very stable expression within the analyzed samples and was used as second reference gene whereas all other target genes (mrjp1-7, mrjp9, and apisimin) were regulated over a wide range (SD: 1.95-5.47).
C q values were determined with the Bio-Rad CFX Manager 3.1 (Bio-Rad) using linear regression for each sample with C q determination mode. Specificity of qPCR products was analyzed with the capillary electrophoresis system QIAxcel (Qiagen) ( Figure S2). PCR efficiency was estimated by serial dilution qPCR (Table S1), and relative target gene expression was determined according to Pfaffl (2001) using the geometric mean of the reference genes mrjp8 and pros26. To determine relative transcript abundance (Table S2), relative gene expression was normalized to total RNA amount by multiplication of both values.

| Protein isolation and SDS-polyacrylamide (PA) gel electrophoresis (GE)
Proteins were precipitated from the flow-through after binding of the RNA to the NucleoSpin ® RNA columns (NucleoSpin ® RNA Kit, Macherey-Nagel) by sodium deoxycholate/trichloroacetic acid according to Arnold and Ulbrich-Hofmann (1999). Of the ten samples isolated per tissue on days zero and eight, five were retained for subsequent quantitative mass spectrometry (see Section 2.42.4).

| Quantitative mass spectrometric analyses
Nano-LC-HD-MSE data were finally acquired for three randomly selected biological replicates per group (HGs and brains, day 0 and day 8) and three technical replicates for each biological replicate.
If a protein in the database was not quantified in a sample, the amount was set to zero for statistical analyses. It should be noted that this does not necessarily mean that the protein was not present in the sample; its amount might be just too low for quantification.
For statistical analyses, the final protein list (Table S4) was further adapted as described in the Supplementary Material and Methods section. The final database included 1,552 different proteins/protein isoforms (Table S3, Tab "New database for mass spec"). All mass spectrometry data have been deposited to the ProteomeXchange Consortium (http://prote omece ntral.prote omexc hange.org) via the PRIDE partner repository (Vizcaino et al., 2013) with the dataset identifier PXD012618.

| Statistics
All statistical analyses were performed with Statistica 8.0 (StatSoft).
Total RNA amount data were log-transformed to achieve normal distribution (Kolmogorov-Smirnov test, p > 0.05) and subsequently ana- Mass spectrometry data did not meet criteria for normal distribution and were thus analyzed with a generalized linear model (GZLM) to reveal major effects (protein, age, tissue, and interactions between these effects). A Venn diagram was generated using VENNY 2.1 to illustrate tissue and time dependent changes (Oliveros, 2007).
Heat maps were built using the MultiExperiment Viewer (MeV, mev. tm4.org) version 4.9, and hierarchical clustering was performed using optimized gene and sample leaf order, Euclidean distance, and average linkage clustering (bootstrap resampling, 1,000 replications) (Eisen, Spellman, Brown, & Botstein, 1998). F I G U R E 2 Heat map of relative transcript abundances in the hypopharyngeal glands and brains of worker honeybees. Transcript abundance is represented for all genes as a color gradient across all samples from light yellow (highest) to deep blue (lowest). Values were log-transformed (nontransformed in brackets) and visualized using the MultiExperiment Viewer (MeV, mev.tm4.org) version 4.9. Mrjp8 is missing in the transcriptional analysis as it was so evenly expressed and neither influenced by tissue nor by age of the bees that mrjp8 was used as second reference gene in addition to pros26. Note: Box-Cox-transformed data were analyzed with one-way ANOVA followed by post hoc Bonferroni test. a-d Transcript abundances of days in the same row with different superscripts are significantly different (p < 0.05). The day with the highest transcript abundance for a specific gene is highlighted in gray, and days that do not differ from the day with the highest transcript abundance are depicted in bold.
In the brain, expression of mrjp1-5, mrjp7, and apisimin followed the same pattern as already observed in the HGs with a relative transcript abundance increase from days 0 to 12 and a further decrease to day 24 ( Figure 2b, Table 1, Table S2). All of these genes were highest expressed at day 12 ( mrjp7-1,589 ± 1,020; mrjp4-944 ± 554; mrjp5-689 ± 360; mrjp6-87 ± 56; and mrjp9-8 ± 4). Albeit highest expressed at day 12, transcript abundance of mrjp9 did not significantly increase from days 0 to 12 but showed a decrease at day 24 ( Figure 2b, Table 1). Mrjp6 did not show any significant difference in transcript abundance in the brain.
Because of the known complex formation of MRJP1 and apisimin in RJ with a stoichiometry of 4:4 (Mandacaru et al., 2017;Tian et al., 2018), we explicitly compared transcript abundance of these two genes which did not differ within the same day between days 0 and 20, but apisimin had significantly more transcripts than mrjp1 on day 24 in both brains and HGs (HG: 44-fold (p < 0.001); brain: 32-fold (p = 0.034)) (one-way ANOVA, p < 0.001, df = 27, F = 43.046; see Figure S3).
When calculating the relative transcript abundances, conspicuous high differences between some of the pools of the same day and gene attracted our attention. As a measure of relative variability, the coefficient of variation was calculated (Table 2). In the HGs, high (≥0.8, highlighted in bold) and very high (≥1.0, highlighted in italics) coefficients of variation were observed for mrjp1-4 and mrjp7 on day 20, mrjp2-4 and mrjp7 on days 16 and 24, and mrjp9 on day 9. In the brain, on day 20 mrjp1-4 and mrjp7, as well as on day 24 mrjp2, mrjp4, and mrjp7, and on day 8 mrjp6 showed very high coefficients of variation (≥1.0, highlighted in bold) ( Table 2).

| Protein amounts
The protein amounts isolated per single brain or per pair of HGs were more than sufficient to be analyzed by quantitative mass spectrometry ( Figure S4). Within the HGs, the high abundance of MRJP1 at day 8 is already visible on SDS-PA gels. Interestingly, the band corresponding to MRJP1, directly isolated from the HGs, migrates at a lower apparent molecular weight (~50 kDa) than when isolated from RJ (~55 kDa) ( Figure S4, MRJP1 marked with asterisk). This discrepancy was already observed upon the first isolation of MRJP1 (Hanes & Šimúth, 1992) but has not yet been clarified.
In total, 1,003 proteins were quantified (Table S4, Figure S5). In the HGs, less proteins were quantified at day 0 (253) than at day 8 (679) whereas in the brain the number of quantified proteins was TA B L E 2 Coefficients of variation (ratio of standard deviation to mean) for mrjp1-7, mrjp9, and apisimin transcript abundances almost equal on both days (767-brain at day 0; 758-brain at day 8). Regarding tissue specificity, 309 proteins were brain-specific but only 114 were HG-specific. Furthermore, 62 proteins were only quantified at day 0 whereas 218 proteins were specifically found at day 8 (Table S4, Figure S5). The dendrogram ( Figure S6)  p < 0.001) as all MRJPs either were found in higher amounts at day 8 compared to day 0 (6.5-to 209.6-fold) or were only detected at day 8 but not at day 0 ( Figure 3). An interaction was also found between age and tissue based on the independent increase in protein amounts from days 0 to 8 ( Figure S4 (Table S6). Within the brain, only MRJP1 was found within the 20 most abundant proteins at day 8 (37.8 ± 9.8 fmol; 5th place) (Table S6).

| D ISCUSS I ON
We here studied, along the time gradient of individual worker bee development, the potential involvement of MRJPs in caste-specific phenotypic plasticity in combination with developmental variance.
Whereas the expression of some mrjps changed with age, others do not seem to be influenced by age or age polyethism: Mrjp8 showed a very low and even expression (C q = 26.75 ± 0.97, mean ± SD), and the protein could not be quantified neither in the brain nor in the HGs. Among the other mrjps, mrjp9 is the lowest expressed mrjp at any time point in both tissues (Figure 2a, Table   S2) and only low amounts of the protein (3.1 ± 2.9 fmol) were quantified solely at day 8 in the HGs. Concordant with that, MRJP8 and MRJP9 were not detected in a comparative proteome study of nurse and forager bee brains (Hernández et al., 2012). Albeit detectable in royal jelly (Zhang et al., 2014), MRJP8 and MRJP9 only represent a minor portion as MRJP1-3 and MRJP5 account for 82%-90% of total food jelly proteins (Schmitzová et al., 1998). Taken together with the fact that mrjp8 and mrjp9 were defined as the most ancestral mrjp gene pair (Buttstedt et al., 2013(Buttstedt et al., , 2014Helbing et al., 2017) and that they were also identified as components of honeybee venom (Peiren et al., 2005(Peiren et al., , 2008, a tissue-specific function in the HGs or the brain and thus an involvement in phenotypic plasticity seem unlikely.

Mrjp1-7 show in the HGs an age-dependent expression pattern
with an increase in transcript abundance from day zero to days 4-12 and a subsequent decrease to day 24 ( Figure 2, Table 1). This is in accord with the general notion that mrjps are higher expressed in brood-feeding nurse bees than in foragers (Klaudiny et al., 1994;Kubo et al., 1996;Ohashi et al., 1997). However, whereas all other mrjps do have their highest transcript abundance during the nurse bee period, mrjp6 shows highest abundance at day 20 (Table 1) in accordance with previous studies reporting on a significantly higher expression of mrjp6 in forager compared to nurse bee heads and HGs (Buttstedt et al., 2013;Liu et al., 2014). In the brain, mrjp6 expression does not show any caste-related modulation supported by Hernández et al. (2012) who did not detect differences in MRJP6 between brains of nurses and foragers. Mrjp6 differs in its time-resolved expression pattern clearly from the other head-expressed mrjp1-5 and mrjp7, and thus, a coregulation of all head-expressed mrjps by the very same transcription factors is unlikely. This is supported by Winkler et al. (2018) who were able to downregulate the expression of mrjp1-3 with 20-hydroxyecdysone, a molting hormone in insects, whereas mrjp4-9 were not affected (Winkler et al., 2018).
However, mrjp6 expression in heads of foragers and nurses is higher than in heads of drones and queens (Buttstedt et al., 2013). Thus, although a nurse-specific function is unlikely, the focus of expression lays in the heads of workers whereas the reproductive castes express 200-to 6,500-fold less mrjp6 in their heads (Buttstedt et al., 2013).
At day 8, MRJPs represent 8.4% of all proteins in the HGs. In royal jelly, MRJPs account for 82%-90% of all proteins (Schmitzová et al., 1998). This difference is due to the fact that royal jelly only contains the proteins that are secreted by the HGs, whereas in our study whole HGs were used as sample. As expected, at day eight MRJP1-3 (66.4-29.7 fmol) and MRJP5 (20.9 fmol) are among the 20 most abundant HG proteins (Table S6). Interestingly, also MRJP4 and MRJP7 were found in medium quantities (13.5 and 13.3 fmol, respectively) albeit both proteins were not detected in early proteome studies on royal jelly (Li, Feng, Zhang, & Pan, 2008;Scarselli et al., 2005;Schmitzová et al., 1998). However, studies that are more recent confirm the presence of MRJP4 and MRJP7 in royal jelly (Feng et al., 2015;Zhang et al., 2014). MRJP1-7 were also quantified in the honeybee brain at day eight, and albeit detected in lesser amounts than in the HGs (1.8-fold), MRJP1 is the fifth most abundant protein at day eight in the brain. Its amount increased from hatching to the nurse bee period 23.4-fold (Table S4) and has been shown to decrease again 9.2-fold from the nurse bee to the forager state . Within the brain, immunolocalization revealed MRJP1 to be present in the antennal and the optical lobes, and in the intercellular spaces in mushroom bodies Kucharski, Maleszka, Hayward, & Ball, 1998;Meng et al., 2018), brain structures involved in associative learning (Menzel, 1993). In the buff-tailed bumblebee Bombus terrestris, the single-copy mrjp-like (mrjpl) gene is also transcribed in the brain (NCBI database BioProject PRJNA383917) and the protein can be immunohistochemically detected in the mushroom bodies (Albert, Spaethe, Grübel, & Rössler, 2014). Albert et al. (2014) suggest therefore that expression of mrjpls in the brain corresponds to the ancestral function rather than to a derived one. However, in A. mellifera the ancestral mrjp9 is not upregulated in heads compared to thoraces and abdomen of workers (Buttstedt et al., 2013), and based on our data, we excluded a tissuespecific function of the ancestral MRJP9 in the brain. However, this does not exclude that the single-copy ancestral mrjpl indeed fulfills, among others, a brain-specific function. This function might have been lost in the honeybee, as new mrjp copies, for example, MRJP1, adopted these functions. MRJP1 can regulate and affect the growth of cells across species (Wan et al., 2018;Watanabe et al., 1996), and thus, the protein might be involved in growth regulation of specific neurons (Kenyon cells) in the mushroom bodies, an idea already raised by Albert et al. (2014).
In royal jelly, the crucial viscosity-determining function of MRJP1  is only accomplished by oligomer and subsequent fibril formation of MRJP1 together with apisimin Mandacaru et al., 2017). The gene encoding apisimin has up to date only been found in the genus Apis and thus seems to be an orphan within the genus with the only so far described function in complex with MRJP1. Both genes show high and similar expression values until day 20 not only in the HGs but also in the brain (Figures   2, S3). However, we were not able to quantify the apisimin protein via mass spectrometry neither in the HGs nor in the brain. This is due to the fact that the 54 amino acids comprising apisimin possess only two unfavorably situated trypsin cleavage sites, resulting after trypsin cleavage in a single lysine, a hexapeptide of only 633.74 Da, and the residual 47 amino acids of the protein (4,796.49 Da). Thus, identification, for which at least two peptides are needed, and subsequent quantification were not possible with a tryptic digest. Furthermore, the plain identification of apisimin via mass spectrometry has been described to be difficult before Tamura et al., 2009).

Mrjp1-5 and mrjp7
show not only in the HGs but also in the brain an expression increase from hatching until the nurse bee period and then a subsequent decrease to day 24. Thus, the initial question, whether the observed elevated expression of mrjp1, mrjp2, mrjp5, and mrjp7 in forager heads is caused by a shift of expression from the HGs to the brain while the bees age, can be clearly answered with no.
But why are there these discrepancies between studies? Indeed, also in Buttstedt et al. (2013) expression of mrjp1, mrjp2, and mrjp7 was found to be fourfold to 36-fold lower in forager compared to nurse bee heads; however, these differences were deemed not significant due to high standard errors, particularly in forager heads (Buttstedt et al., 2013). Also, in this study, very high (≥1.0) coefficients of variation were, except for mrjp6 on day 8 in the brain, exclusively observed for mrjp1-4 and mrjp7 within the forager period at days 20 and 24 (Table 2). In addition, on day 16 in the HGs coefficients of variations for mrjp2-4 and mrjp7 are high (≥0.8). This indicates that individual variation especially for mrjp1-4 and mrjp7 is high during the transition phase from nursing to foraging. The cause for this is that the transition is not a sudden shift but rather characterized by a transition phase where nursing slowly fades out and foraging gradually begins (Seeley, 1982). The physiological changes that workers undergo during age-related polyethism seem to be uniquely timed for each individual worker. Indeed, although the majority of bees start to forage at an age around 20 days, individual bees are known to fly out for their first collecting flights as early as day 10 (Rösch, 1925;Seeley, 1982;zu Oettingen-Spielberg, 1949). In addition, whereas some bees do not continue to feed brood after starting to forage, others do so (Rösch, 1925 lies within the range of the transcript abundances detected at day 8 directly within the nurse bee period (14,649). This explains discrepancies that, for example, MRJP1 is described as "detected in the nurse-bee gland, but not in the forager-bee gland" (Kubo et al., 1996) or as "detected in both the nurse-bee and forager-bee glands, although its density was stronger in the nurse-bee gland" (Ohashi et al., 1997). Thus, conclusions on expression of mrjps should be made with caution, especially when working with forager bees.
Taken together, our results in combination with previous studies suggest the following: Time-resolved expression in the brain follows for each mrjp the expression in the HGs, that is, when expression increases in the HGs, an expression increase is also seen in the brain.
However, expression in the brain is always lower than in the HGs this study the scenario mentioned afore but without showing high (≥0.80) coefficients of variation. However, in previous studies it was described as higher in nurses than foragers, being similar between nurses and foragers, and being higher in foragers than in nurses (Buttstedt et al., 2013;Drapeau et al., 2006;Hernández et al., 2012;Ji et al., 2014). The cause of these differences is still unclear. (III) Mrjp6 expression does not vary significantly in the brain and shows its highest expression in the HGs during the forager bee period. (IV) Mrjp8 and mrjp9 show, compared to the other mrjps, low expression in the HGs and in the brain. Due to their higher occurrence in other body parts of the bees (Buttstedt et al., 2013;Peiren et al., 2005Peiren et al., , 2008, they might not have a HG-or brain-specific function.
Thus, whereas an involvement of mrjp1-7 in caste-specific phenotypic plasticity is possible, especially mrjp8 and mrjp9 do not seem to be influenced by age polyethism.

CO N FLI C T O F I NTE R E S T
None decalred.

AUTH O R CO NTR I B UTI O N S
DD and AB designed the research; DD, DA, and AB performed the research; DD, SE, MF, and AB analyzed the data; AB prepared the original draft; and DD, DA, SE, and AB reviewed and edited the first manuscript version. All authors reviewed and edited the revised version.

DATA AVA I L A B I L I T Y
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (http://prote omece ntral.prote omexc hange.org) via the PRIDE partner repository (Vizcaino et al., 2013) with the dataset identifier PXD012618. Transcriptomic data are provided in the Table S2.