A dietary phytochemical alters caste-associated gene expression in honey bees.

In the eusocial honey bee Apis mellifera, with reproductive queens and sterile workers, a female larva's developmental fate depends on its diet; nurse bees feed queen-destined larvae exclusively royal jelly, a glandular secretion, but worker-destined larvae receive royal jelly for 3 days and subsequently jelly to which honey and beebread are added. RNA-Seq analysis demonstrated that p-coumaric acid, which is ubiquitous in honey and beebread, differentially regulates genes involved in caste determination. Rearing larvae in vitro on a royal jelly diet to which p-coumaric acid has been added produces adults with reduced ovary development. Thus, consuming royal jelly exclusively not only enriches the diet of queen-destined larvae but also may protect them from inhibitory effects of phytochemicals present in the honey and beebread fed to worker-destined larvae.

INTRODUCTION

In the eusocial honey bee Apis mellifera, with reproductive queens and sterile workers, a female larva’s developmental fate largely depends on the diet it receives. Larvae fed exclusively royal jelly, a glandular secretion of nurse bees, become queens, whereas those fed royal jelly for 3 days and subsequently worker jelly containing honey and beebread become sterile workers. The ensuing differences in phenotype are thought to result from either worker castration via nutritional deprivation or queen enrichment by nutritional supplementation. Suboptimal nutrition demonstrably contributes to “trophic castration” and corresponding worker sterility (, ); withholding sugars during the first 3 days of larval life promotes development of sterile workers (–). By contrast, queen-destined larvae consuming only royal jelly develop more rapidly and attain a larger adult body size than do worker-destined larvae (). Attention has historically focused on royal jelly to identify constituents determining developmental fate. The major royal jelly protein royalactin (MRJP1), for example, acts via the epidermal growth factor receptor (Egfr) pathway, which determines cell proliferation, and insulin/insulin-like signaling (IIS)/target of rapamycin (TOR), which contributes to metabolism and growth (–). In addition, DNA methylation status of the genome affects caste differentiation; knockdown of DNA methyltransferase Dnmt3 expression in neonates resulted in queen-destined larval development ().

A key dietary difference between queen and worker larvae that has received little attention is the presence of phytochemicals in the honey and beebread deliberately added by nurses to worker jelly. As described by Haydak () nearly 50 years ago, “Nurse bees apparently recognize the sex and the caste of larvae and may exercise some choice in feeding pollen to larvae of different castes.” Adding pollen (in the form of beebread, or processed pollen) does not substantially increase the protein content of modified worker jelly, so the function of this behavior has never been clear. Both honey and beebread contain plant-derived phenolics (), which are known to function as modulators of DNA methylation, intracellular signaling networks, chromatin structure, and microRNA (, ). By feeding queen-destined larvae only royal jelly, nurse bees may also protect these larvae from biologically active phytochemicals in honey and beebread. Mao et al. () demonstrated that p-coumaric acid, a virtually universal constituent of honey and beebread, up-regulates multiple detoxification genes in worker adults. Here, we examined the effects of p-coumaric acid on larvae to test the hypothesis that a food-based phenolic constituent of beebread and honey can influence female larval developmental pathways in A. mellifera.

RESULTS

An analysis aimed at determining the amount of p-coumaric acid ingested during larval development (Materials and Methods) revealed that, whereas royal jelly is entirely devoid of detectable phenolic acids, honey and beebread in hives at the University of Illinois at Urbana-Champaign (UIUC) contained 5.88 ± 0.12 μg/g and 31.92 ± 5.74 μg/g (mean ± SEM, with a high value of 38.5 μg/g), respectively. For bioassay, we chose a concentration of 82 μg/g, because at this level, it is nontoxic to adults () and within the general magnitude of the concentrations we found in local honey and beebread (likely derived from soybean, Glycine max), particularly in view of the fact that p-coumaric acid is one of multiple phenolic constituents found in honey and beebread that may similarly influence gene expression ().

We conducted an RNA-Seq analysis to identify differences in gene expression in larvae reared on a semi-artificial diet with and without p-coumaric acid, as well as an in vitro larval rearing experiment to determine the effects of ingesting this compound on adult female reproductive status (Materials and Methods). To obtain sufficient sequencing depth, and thereby increase the number of genes discovered, and assess with accuracy genes with less abundant transcripts, three samples per lane (with 10 to 12 whole larvae per sample) were sequenced using Illumina HiSeq 2000. After a stringent quality filtering process, more than 50 million high-quality reads for each sample were obtained; about 40 million of these reads for each sample (75.5 to 76.0%) were uniquely aligned to the reference genome (table S1) using Bowtie 2 ().

Before differential expression analysis, the variation between the three biological replicates and between the two treatments (rearing on diets with and without 0.5 mM p-coumaric acid) was assessed using a multidimensional scaling (MDS) plot. Biological replicates of control treatments clustered closely, whereas those of the 0.5 mM p-coumaric acid treatments did not, indicating variation among replicates (in the first dimension) and among library preparations (in the second dimension) (Fig. 1). Despite this variation, the samples of both treatments were well separated in the first dimension, indicating that the sequencing data were qualified for differentially expressed gene (DEG) identification. To investigate p-coumaric acid effects on larval development, we used edgeR () to identify DEGs between control and 0.5 mM p-coumaric acid treatments (Fig. 2). Among the 5344 DEGs [false discovery rate (FDR), <0.05], 2663 were up-regulated (representing about one-third of the genome) and 2681 were down-regulated (table S2).

Fig. 1

edgeR’ multidimensional scaling plot (MDS) showing relationships among all pairs of samples (two treatments) and replicates (3 replicates per treatment).

Treatments comprise larval honey bees consuming artificial diet and larval honey bees consuming artificial diet with p-coumaric acid. CK denotes control treatment; CA denotes 0.5 mM p-coumaric acid treatment. Control and treatment samples group separately when plotted.

Fig. 2

edgeR’ smearplot showing differences in gene expression in larval honey bees consuming artificial diet with and without p-coumaric acid; red dot indicates gene that is differentially expressed (FC = fold change in expression).

To gain insight into the biological processes associated with the DEGs, we performed DAVID (Database for Annotation, Visualization and Integrated Discovery) functional annotation clustering analysis (, ) on 2077 DEGs (logFC < −0.6 and logFC > 0.6) of the 3367 DEGs that have Drosophila orthologs with FlyBase annotations. Among 276 clusters, there were three clusters with enrichment scores higher than 5 (Table 1 and table S3). In cluster 1 (enrichment score: 11.08), all of the 17 enriched terms are related to larval development. In cluster 2 (enrichment score: 9.73), four enriched terms are associated with gene transcription in the nucleus and mRNA metabolism. In cluster 3 (6.98), two enriched terms concern noncoding RNA (ncRNA) processing and metabolism, and two others concern transfer RNA (tRNA) processing and metabolism.

Table 1

DAVID functional annotation clustering analysis on 2077 of the 3367 DEGs in honey bee larvae consuming p-coumaric acid.

Annotation cluster	Term	Count	P	Benjamini
Cluster 1 (enrichment score: 11.08)	Imaginal disc development	131	8.90 × 10−16	2.40 × 10−12
	Post-embryonic organ development	100	7.10 × 10−14	9.40 × 10−11
	Imaginal disc morphogenesis	96	1.80 × 10−13	1.60 × 10−10
	Post-embryonic organ morphogenesis	96	1.80 × 10−13	1.60 × 10−10
	Post-embryonic morphogenesis	109	2.70 × 10−12	1.80 × 10−9
	Metamorphosis	111	2.80 × 10−12	1.50 × 10−9
	Wing disc development	92	4.90 × 10−12	2.20 × 10−9
	Instar larval or pupal morphogenesis	107	4.90 × 10−12	1.90 × 10−9
	Post-embryonic development	126	5.30 × 10−12	1.80 × 10−9
	Instar larval or pupal development	120	4.10 × 10−11	1.10 × 10−8
	Imaginal disc–derived appendage morphogenesis	79	1.10 × 10−10	2.80 × 10−8
	Wing disc morphogenesis	74	2.00 × 10−10	4.30 × 10−8
	Post-embryonic appendage morphogenesis	76	2.20 × 10−10	4.50 × 10−8
	Appendage morphogenesis	79	2.20 × 10−10	4.20 × 10−8
	Imaginal disc–derived appendage development	79	2.80 × 10−10	4.90 × 10−8
	Imaginal disc–derived wing development	73	3.30 × 10−10	5.50 × 10−8
	Appendage development	79	5.40 × 10−10	8.40 × 10−8
Cluster 2 (enrichment score: 9.73)	Transcription	131	7.10 × 10−12	2.10 × 10−9
	Regulation of transcription	192	1.10 × 10−9	1.60 × 10−7
	Regulation of transcription, DNA-dependent	139	9.40 × 10−6	5.80 × 10−4
	Regulation of RNA metabolic process	147	6.60 × 10−5	3.10 × 10−3
Cluster 3 (enrichment score: 6.98)	ncRNA metabolic process	49	1.30 × 10−8	1.80 × 10−6
	ncRNA process	38	2.30 × 10−8	3.00 × 10−6
	tRNA metabolic process	36	5.10 × 10−7	4.70 × 10−5
	tRNA process	23	7.60 × 10−7	6.30 × 10−5

To characterize larval response profiles induced by p-coumaric acid at the level of a single pathway, we performed an RNA-Seq data pathway and gene-set analysis by using the two R packages GAGE and Pathview (, ). Among the two gene sets of signal transduction and metabolism extracted from the KEGG pathway database (www.genome.jpkeggpathway.html), only the Hippo signaling pathway, which contributes to regulating organ size in animals, was enriched (P < 0.1) (Fig. 3), suggesting that p-coumaric acid influences organ size in developing larvae. In this pathway, p-coumaric acid altered the expression of 26 genes among the 46 total genes (54.2%) (table S4). Among the five down-regulated genes, the only fold change less than −2.00 was that of Mats (−2.04). In contrast, 7 of the 21 up-regulated genes had fold changes greater than 2.00; these included Ds, Ft, Atrophin, Zyx, Ex, dJub, and Km.

Fig. 3

Hippo signaling pathway.

Up-regulated genes appear in red, and down-regulated genes appear in green. For specific fold change and P values, see table S3.

As in adult workers (), larval consumption of p-coumaric acid results in the up-regulation of multiple xenobiotic-metabolizing cytochrome P450 genes and immunity genes (Table 2), generally congruent with those up-regulated by p-coumaric acid in adults (). Of the nine up-regulated xenobiotic-metabolizing genes, eight are P450 genes and one is the γ-glutamyltranspeptidase-1–like gene. Among the P450 genes, CYP9Q1 was up-regulated 1.85-fold, and the six CYP6AS genes were up-regulated from 1.88- to 47.11-fold. Also in the context of xenobiotic metabolism, the nuclear xenobiotic sensor Hr96 gene () was down-regulated 2.11-fold in the p-coumaric acid treatment. In the context of immunity, six genes were up-regulated, including the antimicrobial peptide genes apidaecin type 73–like and apid73 apidaecin, induced ~25- and 11-fold, respectively. Other genes related to immunity were up-regulated 1.68- to 1.87-fold; no immunity genes were down-regulated in response to p-coumaric acid.

Table 2

Genes differentially regulated by at least 1.4-fold by p-coumaric acid in honey bee larvae.

Group	Gene ID	Gene name	Fold change	Adjusted P
Xenobiotic metabolism genes	GB40288	CYP6AS1	1.88	1.21 × 10−6
	GB49886	CYP6AS2	2.55	1.27 × 10−12
	GB49887	CYP6AS3	4.83	7.67 × 10−20
	GB49885	CYP6AS4	5.33	8.44 × 10−23
	GB49890	CYP6AS5	47.11	1.98 × 10−40
	GB49888	CYP6AS12	15.32	4.69 × 10−25
	GB43713	CYP9Q1	1.85	1.10 × 10−7
	GB43727	CYP9Q2	1.84	2.67 × 10−8
	GB43728	CYP9Q3	2.25	5.52 × 10−8
	GB43714	CYP9P1	2.57	1.70 × 10−5
	GB43709	CYP9P2	1.46	5.33 × 10−4
	GB43716	CYP9R1	1.79	6.20 × 10−6
	GB43715	CYP9S1	3.12	1.04 × 10−13
	GB51356	CYP4G11	1.68	4.09 × 10−7
	GB45654	γ-Glutamyltranspeptidase 1–like	1.29	0.020213
Immunity genes	GB46236	Apidaecins type 73–like	25.41	5.92 × 10−17
	GB51306	Apid73 apidaecin	10.90	2.79 × 10−15
	GB40163	LYZ1	1.75	1.31 × 10−5
	GB42685	β-1,3-Glucan recognition protein 1	1.68	0.005207
	GB47805	Peptidoglycan recognition protein S2	1.87	5.78 × 10−6
	GB45648	Immune deficiency	1.74	4.76 × 10−7
Caste differentiation genes	GB48059	Ecdysone receptor	2.31	2.69 × 10−13
	GB47037	Hormone receptor–like in 4	2.00	9.95 × 10−5
	GB49105	Ecdysteroid-regulated gene E74	2.28	3.59 × 10−15
	GB42692	Ultraspiracle	2.11	3.60 × 10−9
	GB40074	Probable nuclear hormone receptor HR38–like	2.00	0.00236
	GB45414	Hormone receptor–like in 39	1.29	0.01531
	GB54477	EGFR	1.64	1.28 × 10−6
	GB46903	S6K	1.51	7.88 × 10−5
	GB55485	Dnmt3	1.78	4.20 × 10−5
	GB55425	InR2	2.07	7.65 × 10−9
	GB53353	InR	2.06	0.000109
	GB49911	IRS	2.05	4.84 × 10−11
	GB43560	ILP2	−1.91	7.05 × 10−6
	GB48301	Forkhead box protein O	2.38	2.69 × 10−17
	GB44043	Juvenile hormone acid methyltransferase	2.73	1.42 × 10−13
	GB55205	Major royal jelly protein 1	−6.61	7.32 × 10−8
	GB55207	Major royal jelly protein 6	−3.80	0.003876
	GB55206	Major royal jelly protein 4	−3.76	0.007108
	GB55208	Major royal jelly protein 5	−2.95	0.020384
	GB55204	Major royal jelly protein 3	−2.56	0.00404
	GB55213	Major royal jelly protein 7	−1.89	0.005034
	GB55212	Major royal jelly protein 2	−1.95	0.002019
	GB55211	Major royal jelly protein 2–like (pseudo)	−2.56	2.36 × 10−5

Because both DEG functional annotation and pathway analysis showed that p-coumaric acid influenced gene expression linked with organ development and maturation and transcription in nuclei, we analyzed DEGs associated with caste differentiation to determine whether it was similarly affected (Table 2). Other than ILP-2, which was down-regulated 1.91-fold, 14 genes known to influence worker-queen caste differentiation were up-regulated from 1.29- to 2.73-fold. Surprisingly, major royal jelly protein 1 to 7 genes and major royal jelly protein 2–like gene (a pseudogene) were highly expressed in larvae consuming control diet, whereas all of these genes were down-regulated in larvae reared on the p-coumaric acid–containing diet. Particularly striking was the 6.61-fold down-regulation of the royalactin (MRJP1) gene.

Because the results of our RNA-Seq analysis collectively suggested the potential for p-coumaric acid to alter larval developmental fate, we reared the larvae in vitro on nutrient-rich diet containing royalactin, a queen-inducing feeding regimen (), with and without p-coumaric acid. Within 12 hours of eclosion (), adults were dissected and the level of ovary development was visually scored on a standard five-point scale (). Reduced ovary development was exhibited in bees reared on a diet containing p-coumaric acid (Fisher’s exact test, P < 0.0001), with maximum ovary development observed in bees reared on a diet lacking p-coumaric acid (Fig. 4).

Fig. 4

Effect of p-coumaric acid ingestion during larval growth on ovary development in A. mellifera.

(A) Ovary development scale of 1 to 5 (undeveloped to fully developed) modified from Pernal and Currie (). (B) Distribution of ovaries according to level of ovary development in newly eclosed adult bees reared in vitro on diets with and without 0.5 mM p-coumaric acid (two-tailed Fisher’s exact test, P < 0.0001; n = 130).

DISCUSSION

By consuming pollen, nectar, honey, and beebread, adult workers consume substantial quantities of phenolics, which can interfere with multiple signaling pathways and modulate epigenetic changes (, ). Adults use multiple xenobiotic-detoxifying enzymes, including midgut cytochrome P450s, to process these phenolics. Quercetin, for example, a flavonoid ubiquitous in honey and beebread, is detoxified by P450s in both CYP6AS and CYP9Q subfamilies (, ). Multiple flavonoids and phenolic acids in honey and beebread up-regulate the P450 genes encoding key detoxification enzymes in adults (, ). Of these, CYP6AS2, CYP6AS3, CYP6AS4, CYP6AS5, and CYP6BD1 in nurse midguts were up-regulated 1.90- to 3.11-fold by 1 mM p-coumaric acid, and CYP9Q3 was up-regulated 2.55-fold (). The RNA-Seq data presented here indicate that xenobiotic-metabolizing P450 genes are up-regulated even more markedly by 0.5 mM p-coumaric acid in 3-day-old larvae (Table 2). Six CYP6AS genes, six CYP9 genes, and the CYP4G11 gene are up-regulated 1.88- to 47.11-fold, 1.46- to 3.12-fold, and 1.68-fold, respectively. This up-regulation is likely required to enable larvae to process the multiple phenolics they encounter in honey and beebread in worker jelly for the remainder of their development. That this up-regulation is associated with the switch to worker jelly with its phytochemical supplements is indicated by Cameron et al. (), who examined a cohort of neonates either grafted into queen-rearing cells and placed in a queen-less hive or raised as workers in a queen-right hive. At 84 hours, comparable to the developmental stage of larvae examined in our study, P450 genes up-regulated in workers reared in the hive and fed by nurse bees included CYP6AS2, CYP6AS3, CYP6AS4, CYP6AS10, CYP6AR1, and CYP9Q2. By comparison, in our study, four of the same six P450 genes were up-regulated in larvae consuming p-coumaric acid (CYP6AS2, CYP6AS3, CYP6AS4, and CYP9Q2), as were five other related xenobiotic-metabolizing P450s (CYP6AS1, CYP6AS5, CYP6AS12, CYP9Q1, and CYP9Q3). Notably, Cameron et al. () found only a single P450 gene up-regulated in larvae consuming only royal jelly—CYP15A1, a Clan 2 gene not known to be involved in detoxification.

In addition to detoxification genes, feeding larvae a diet containing 0.5 mM p-coumaric acid resulted in the discovery of 5344 DEGs. Functional analyses of these DEG functions indicate that p-coumaric acid modulates expression of genes related to larval development, transcription, and Hippo signaling (Table 2, Fig. 3, and tables S3 and S4). This compound may drive worker development via modulating the two essential pathways of IIS/TOR and EGFR, as well as the DNA methylation status for caste determination. The IIS network is implicated in determining female caste fate in honey bees; of the two insulin-like peptides AmILP1 and AmILP2, AmILP1 regulates juvenile hormone production, whereas AmILP2 is necessary for ovary development, because its knockdown causes a reduction in ovary size (). In our study, consuming p-coumaric acid down-regulates AmILP-2 in larvae (Table 2) and reduces ovary development in adults (Fig. 4). Royalactin in royal jelly is thought to trigger queen development via the Egfr pathway (); in our study, Egfr in larvae, which requires royalactin to function, was up-regulated by p-coumaric acid. At the same time, the royalactin gene was highly expressed in larvae consuming the control diet but down-regulated in larvae consuming p-coumaric acid (Table 2).

The production of royalactin in worker larvae (–) suggests that impacts of this protein on developmental fate may involve endogenous processes in larvae beyond royal jelly consumption; these endogenous processes may be regulated by the presence or absence of beebread- or honey-derived phytochemicals such as p-coumaric acid.

Epigenetic processes bridge intrinsic and environmental signals through DNA cytosine-5-methyltransferases (Dnmts), histone deacetylases (HDAs), and methyl-binding proteins (MBPs) (). A. mellifera has a functional CpG methylation system consisting of two orthologs of Dnmt1, one ortholog of Dnmt2, and one ortholog of Dnmt3 (). Larval knockdown of DNA methyltransferase 3 (Dnmt3) causes development of queen-like traits in worker-destined larvae (); in our study, Dnmt3 expression was altered (up-regulated) by p-coumaric acid (Table 2). In bee larvae, nutrition modulates juvenile hormone III (JHIII) and 20-ecdysone titers via IIS/TOR and EGFR pathways to control caste differentiation. The fact that the Egfr pathway is responsible for 20-ecdysone synthesis suggests that its activation decreases development time for queen-destined larvae (), whereas AmILP1 regulates JHIII production via IIS/TOR to prevent autophagic programmed cell death in the ovary, thus increasing ovary size (, ). Nine receptors, S6K, IRS, FOX, and the JHA methyltransferase genes related to the two pathways and 20EJH III syntheses were up-regulated by p-coumaric acid (Table 2), a pattern of expression identical to that found in queen-destined larvae; p-coumaric acid appears to influence the developmental fate of female honey bee larvae via directly regulating expression of genes in caste determination pathways.

The fact that so many genes known to influence developmental fate of female larvae in A. mellifera are differentially expressed upon consumption of p-coumaric acid suggests that, for the honey bee, consuming biologically active phytochemicals may influence reproductive development in worker-destined larvae. Reliance upon royal jelly as the exclusive food for queen-destined larvae and for the long-lived adult queen may have evolved in concert with the exploitation of the biological activity of phenolics for caste regulation via a form of chemical castration. Reproductive development of solitary hymenopterans that consume pollen and nectar is presumably unaffected by phytochemicals found in these foods. Whether the pronounced regulatory effects of nectar and pollen phytochemicals on honey bees are widespread among social hymenopteran pollinators in general or are unique to the highly eusocial honey bee remains an open question, the answer to which may shed light on the impact of diet and food processing on the evolution of eusociality.

MATERIALS AND METHODS

Overall experimental design

Mao et al. () demonstrated that p-coumaric acid, a virtually universal constituent of bee foods (honey and beebread), up-regulates multiple detoxification and immunity genes. In view of its influence on gene regulation in adults, our overall objective was to examine the effects of p-coumaric acid on larvae, to test the hypothesis that ingestion of plant-derived phenolics plays an essential role in determining female larval developmental pathways in A. mellifera. We conducted RNA-Seq to identify differences in gene expression in larvae reared on a royal jelly diet with and without p-coumaric acid as well as an in vitro rearing experiment to determine the impact of p-coumaric acid ingestion on reproductive status in adult females.

Quantification of p-coumaric acid present in larval foods

To determine the appropriate concentrations of p-coumaric acid to use in rearing, we collected pollen, honey, and beebread from hives at the University of Illinois Bee Research Facility and subjected them to analysis by high-pressure liquid chromatography (HPLC). Royal jelly samples were obtained from a commercial supplier (vide infra). For each sample, 1 g of each product was placed in a test tube, to which 1 ml of water and 5 ml of methanol were added; the solution was then mixed on a tube rotator for 4 hours. Pollen, beebread, and royal jelly samples were centrifuged at 10,000g for 10 min at room temperature. To achieve the same volume as the honey samples, methanol was then added to extract precipitates in the pollen, beebread, and royal jelly samples for an additional extraction, and the two supernatants were combined. All samples for HPLC analysis were filtered with 0.22-μm filters. For each bee product, three independent samples were prepared. The HPLC analyses were performed on a Phenomenex Luna phenyl-hexyl column (250 × 3–mm inside diameter; 5-μm particle size; 100 Å pore size) using the following gradient run (1 ml/min) with water containing 0.1% formic acid (A) and acetonitrile containing 0.1% formic acid (B): 10% B from 0 to 10 min, and 10 to 57% B from 10 to 40 min. The content of p-coumaric acid was quantified at 310 nm by comparison with a standard curve.

We conducted a preliminary rearing experiment subjecting larvae to diets containing a range of concentrations of p-coumaric acid approximating those potentially encountered by larvae consuming honey and beebread () to identify levels that do not suppress growth and thus are unlikely to have outright toxic effects. Over this range of concentrations (0.1, 0.25, 0.5, and 1.0 mM), body masses of larvae reared through fifth instar (n = 15 to 19 per treatment) were unaffected (mean masses ranged from 173.2 to 174.4 mg; P = 0.92, 0.37, 0.18, and 0.55, respectively). Accordingly, we chose an intermediate concentration of 0.5 mM p-coumaric acid to rear larvae for RNA-Seq experiments to determine the impact of this ubiquitous natural diet constituent on gene expression in honey bee larvae.

Chemicals

Fresh-frozen royal jelly (organic) was purchased from GloryBee Foods; p-coumaric acid, d-glucose, and d-fructose were obtained from Sigma-Aldrich. Bacto yeast extract was obtained from Becton, Dickinson and Company.

Larval honey bee diet preparation

A semi-artificial diet was prepared with 50% royal jelly and 50% aqueous solution containing 12% d-glucose, 12% fructose, and 2% yeast extract. For preparation of diets with varying concentrations of p-coumaric acid, ecologically appropriate amounts of this compound () were added into the semi-artificial diets to yield concentrations of 0, 0.1, 0.25, 0.5, and 1 mM (w/v). The prepared diets were stored at −20°C before use in all experiments.

Larval grafting and rearing

All first instar honey bee larvae in these experiments were obtained from the University of Illinois Bee Research Facility in Urbana, IL. To collect eggs, a previously occupied empty frame with an excluder cage was placed into a hive with a healthy colony in the afternoon for cleaning by workers. A queen was transferred into the excluder cage in the evening of the next day, confined overnight, and released. In the morning of the fourth day after egg collection, the frame with newly hatched larvae was removed from the hive and brought to the laboratory for larval grafting. To evaluate p-coumaric acid effects on larval growth, 24 newly hatched larvae for each treatment (0, 0.1, 0.25, 0.5, and 1 mM) were grafted into the wells of a 24-well cell culture plate with 30 μl of the corresponding diet in each well. Our goal was to maximize the likelihood of detecting effects of beebread and honey ingestion, and the available literature is conflicting as to the age at which larvae are developmentally “bipotential,” with some literature indicating that even 3-day-old workers may begin to display developmental differences [reviewed in ()]. Because there are few records of the exact point at which beebread/honey feeding begins [with Haydak () the principal source of information on the timing], we began our treatments on the first day to be maximally inclusive. For RNA-Seq experiments, 12 newly hatched larvae for each sample for each of two treatments, 0 and 0.5 mM, were transferred as before, and each treatment was independently replicated three times (for a total of six libraries). The larval rearing plates were placed into a desiccator maintained at a relative humidity of 96% (saturated K2SO4) in a 34° ± 1°C dark incubator. Each day, excess diet in the wells was removed using a vacuum and replaced with fresh diet prewarmed to 34°C. On the third day, the larvae for the RNA-Seq experiments were collected; excess surface fluid was removed with soft tissue paper, frozen in liquid nitrogen, and stored at −80°C for total RNA extractions. On the fifth day of the experiment assessing the effects of p-coumaric acid on growth, each of the larvae from which excess fluid had been removed was weighed.

RNA extraction, complementary DNA library construction, and sequencing

After 10 to 12 larvae for each sample were ground in liquid nitrogen with mortar and pestle, half of the powder was added into 1 ml of TRIzol and total RNA was extracted following the TRIzol protocol. Six hundred microliters of the supernatant containing total RNA was mixed with an equal volume of 70% ethanol and transferred to a Qiagen RNeasy Mini kit column. The extraction was completed following the Qiagen protocol, which included an on-column DNase digest (Qiagen). Complementary DNA library construction and RNA sequencing were carried out at the W.M. Keck Center for Comparative and Functional Genomics at the Roy J. Carver Biotechnology Center at the UIUC. The RNA-Seq libraries were prepared using the TruSeq RNA-Seq Sample Prep kit according to the manufacturer’s instructions (Illumina) and sequenced using Illumina HiSeq 2000 for 100–base pair single-end reads (three samples per lane).

RNA-Seq analyses

To serve as the reference genome, 17 FASTA files were downloaded from BeeBase (ftp://ftp.ncbi.nih.gov/genomes/Apis_mellifera/Assembled_chromosomes/seq/). After the reads for each sample were mapped to the reference genome using Bowtie 2 with a “sensitive” option (), the resulting SAM files were converted to BAM files using SAMTools (). Read counts per gene were calculated for each sample with HTSeq-count (www-huber.embl.de/users/anders/HTSeq/doc/overview.html) using the BAM file and the genome annotation file (amel_OGSv3.2.gff3) as the input. Finally, edgeR using qCML tagwise dispersion model was used to identify DEGs between larvae fed a p-coumaric acid–free diet and those fed a diet containing 0.5 mM p-coumaric acid (), and significant DEGs were selected at FDR <0.05.

Analysis of the biological processes affected by p-coumaric acid

Biological processes that were statistically overrepresented with these DEGs were identified using the DAVID functional annotation clustering tool (, ). After BeeBase IDs of these DEGs were converted into FLYBASE_GENE_IDs of their fruit fly orthologs, the resulting DEGs (logFC < −0.6 and logFC > 0.6) were uploaded and analyzed with GOTERM_BP_FAT selected.

Pathway analysis of DEGs

Two R packages, GAGE and Pathview, were used for pathway analyses of these DEGs (, ). After two honey bee gene sets of signal transduction and metabolism were extracted from the KEGG Pathway database (www.genome.jp/kegg/pathway.html), the analyses were carried out according to the workflow jointed with edgeR on RNA-Seq Data Pathway and Gene-set Analysis Workflows. (www.bioconductor.org/packages/2.13/bioc/vignettes/gage/inst/doc/RNA-seqWorkflow.pdf).

Queen rearing

In a separate bioassay, newly hatched larvae were collected and reared using in vitro queen-rearing techniques as described by Mutti et al. (). For the control treatment, 96 larvae were fed the standard diets; for the treatment group, 96 larvae were reared on an identical diet, supplemented with 0.5 mM p-coumaric acid. Within 12 hours of eclosion, adult bees were dissected and the level of ovary development was visually scored on a scale of 1 to 5, modified from Pernal and Currie (), corresponding to the following: 1, underdevelopment (the diameter of the ovary posterior end is substantially smaller than that of the oviduct anterior end); 2, slight development (the diameter of the ovary posterior end is slightly expanded but still smaller than that of the oviduct anterior end); 3, intermediate development (the diameter of the ovary is uniform throughout its length); 4, advanced development (diameters of ovary posterior ends are significantly larger than those of oviduct anterior ends); 5, full development (ovaries are typical of virgin queens) (Fig. 4).

We used the standard methodology for in vitro queen rearing, which in reality does not invariably produce queens; generally, in vitro rearing produces a range of adults exhibiting a continuum of reproductive characteristics. Our goal was not to demonstrate that p-coumaric acid in the diet is in and of itself sufficient to produce fertile queens; rather, it was to determine whether its presence in the diet can influence caste-associated traits. Thus, we used a measure of ovary size at emergence as an indicator of developmental fate because caste is already determined before emergence (with ovary activation occurring post-emergence in response to other environmental conditions).