Epigenetic Modification of Gene Expression in Honey Bees by Heterospecific Gland Secretions

Background In the honey bee (Apis mellifera), queen and workers have different behavior and reproductive capacity despite possessing the same genome. The primary substance that leads to this differentiation is royal jelly (RJ), which contains a range of proteins, amino acids, vitamins and nucleic acids. MicroRNA (miRNA) has been found to play an important role in regulating the expression of protein-coding genes and cell biology. In this study, we characterized the miRNAs in RJ from two honey bee sister species and determined their possible effect on transcriptome in one species. Methodology/Principal Findings We sequenced the miRNAs in RJ either from A. mellifera (RJM) or A. cerana (RJC). We then determined the global transcriptomes of adult A. mellifera developed from larvae fed either with RJM (mRJM) or RJC (mRJC). Finally we analyzed the target genes of those miRNA that are species specific or differentially expressed in the two honey bee species. We show that there were differences in miRNA between RJM and RJC, and that transcriptomes of adult A. mellifera were affected by the two types of RJ. A high proportion (23.3%) of the affected genes were target genes of differential miRNAs. Conclusion We show for the first time that there are differences in miRNAs in RJ between A. mellifera and A. cerana. Further, the differences in transcriptomes of bees reared from these two RJs might be related to miRNA differences of the two species. This study provides the first evidence that heterospecific royal jelly can modify gene expression in honey bees through an epigenetic mechanism.


Introduction
The Western honey bee (Apis mellifera) is one of the most important economical insects because of its crucial role in pollination [1]. A honey bee colony is composed of three castes, a fertile queen, hundreds of haploid drones, and thousands of nearly sterile workers [2,3]. Despite their identical genome, the queen and her workers exhibit vast differences in morphology, behavior, physiology, reproduction and longevity [4][5][6]. The primary substance that leads to this is royal jelly (RJ), which is a yellow milky substance from ''nurses'' with developed hypopharyngeal and mandibular glands [7][8][9].
RJ contains a range of proteins, carbohydrates, lipids, minerals, vitamins, and a large number of bioactive substances, especially immunological peptides and antibacterial proteins [10][11][12][13]. Major Royal Jelly Proteins (MRJPs) are the prime RJ ingredient, which are crucial in regulating reproductive maturation [11]. Royalactin is a 57-kDa protein, which can induce larvae developing into queens [14]. Royalactin helps increase body size, promote ovary development and shorten the developmental time. In addition, RJ also contains small amounts of nucleic acids. One study found that RJ contains both DNA and RNA, and there are quantitative differences in nucleic acids in fresh RJ between A. mellifera and A. cerana [13]. The most recent discovery is that RJ contains microRNAs which may play a role in caste differentiation [15].
While the mechanism by which miRNA modulate gene expression has been well studied [17,21,23], it is not clear whether there are differences in RJ originated miRNAs between two honey bee species (A. mellifera and A. cerana), and whether feeding A. mellifera with different RJ might cause differences in transcriptome. A. cerana is considered to the most closely related species to A.mellifera [33]. The two species share common features (such as open nesting) [34] but show differences in behaviors and physiology [35]. It is likely that the two species diverged from a common ancestor around three million years ago [36]. In this study we report the differences in miRNA of RJ either from A. mellifera (RJM) or A. cerana (RJC). We also tested the hypothesis that differentially expressed miRNAs in the RJ of the two species affect A. mellifera transcriptome.

Differences in miRNA between Royal Jelly of Two Different Species of Honey Bees
A total of 11,828,863 reads from the RJM library and 10,289,838 reads from the RJC library were obtained after discarding the empty adapters (Table 1). After discarding those sequences that were of low-quality, shorter than 18 nucleotides and single-read sequences, 10,318,386 and 9,493,118 reads for the RJM and RJC remained for analysis respectively. RNAs sequenced by Solexa were in the length of 10-44 nucleotides (nt), and length distributions of small RNAs in the two libraries were significantly different (Contingency Table Analysis, X 2 .1000, P,0.001, Fig. 1). miRNAs, those with 20-22 nt, in the two types of RJ also showed significantly different length distributions (Contingency Table Analysis, X 2 .1000, P,0.001).
Subsequently, small RNAs were classified into different categories according to their biogenesis and annotations ( Table 2). Both RJM and RJC contained several heterogeneous small RNA species which included miRNAs, degraded rRNA fragments and mRNA fragments. In royal jelly, miRNAs were the major fraction of small RNA species. As shown in Table 3, the total reads of miRNAs in RJC (1,872,895 reads) were higher compared to RJM (1,735,052 reads).
Solexa sequencing and RNA classification indicated that expression profiles of miRNAs in RJM and RJC are significantly different. By referencing to the mirBase release 13.0 [37], we identified 69 and 48 known miRNAs in RJM and RJC, respectively. There were 23 miRNAs specific to RJM, 2 miRNAs specific to RJC, and 46 shared in both RJ ( Table 4). The average expression level of all miRNA in RJM was about 8-fold higher than that of RJC (Table S1). Among these, RJM contained 31 upexpressed, 6 equally expressed, and 2 down-expressed miRNAs compared to RJC (Fig. 2). According to sequence homology, we noticed a high-percentage of miRNA from categories of metabolic process, cell part and catalytic activity (Fig. 3). Cellular process, cell and binding terms were dominant.

Transcriptome Modifications in A. mellifera by Two Different RJs
To test the hypothesis that miRNAs in RJM and RJC affect A. mellifera transcriptome, we determined the global transcriptomes of adult A. mellifera developed from larvae fed either with RJM (mRJM) or RJC (mRJC). The total number of reads for mRJM and mRJC were 48,971,186 and 49,358,642, respectively ( Table 5). The distributions of perfectly matched reads to the honey bee genome in mRJM and mRJC were not significantly different (Contingency Table Analysis, X 2 = 1.14, P.0.2), however, those uniquely matched to the genome in the two types of bees were significantly different (X 2 = 150, P,0.0001). When the mapping was to the honey bee genes (which do not include intron and UTRs), both perfectly matched and uniquely matched reads were significantly differently distributed in mRJM and mRJC (X 2 .150, P,0.0001 in both cases).
By use of the Reads Per kb per Million reads (RPKM) method [38], we explored gene expression levels in mRJM and mRJC. This method was adopted to eliminate the influence of variation in gene length and the total reads number. The number of downregulated genes was more than two times that of up-regulated genes, with 439 down-regulated and 179 up-regulated genes ( Fig. 4, Table S2). We systematically examined every differentially expressed gene (DEG) in order to identify genes involved in important pathways. As shown in Table S3, these DEGs were mainly located in endocytosis, focal adhesion, metabolic pathways, regulation of actin cytoskeleton, and RNA transport. Some DEGs were related to pathways on caste differentiation, such as insulin signaling, mTOR, and MAPK. According to sequence homology, we obtained DEGs from categories in metabolic process, cell part and catalytic activity (Fig. 5). Cellular process, cell and binding terms were dominant.

Analysis of miRNA Targeted Genes
We identified the target genes of the following miRNAs: 23 RJM-specific, 2 RJC-specific and 33 differentially expressed miRNA (Table S4). Among the 618 differentially expressed genes between mRJM and mRJC, 144 (23.3%) genes were identified as target genes of miRNAs (Table S4).

miRNA Differences in RJ of A. mellifera and A. cerana
miRNA in honey bees have been shown to correlate with behavioral plasticity [27][28][29]. In our paper, we were more concerned with miRNAs in royal jelly that may play roles in affecting transcriptome. High-throughput sequencing of small RNAs indicated that there were many small RNAs in the two types of royal jelly (Fig. 1). Small RNAs are of 18-35 nt, which includes three major types: miRNA (20-22 nt), siRNA (24-26 nt), and piRNA (32-34 nt). RJC has a higher percentage of piRNA, especially those with 33 nt (Fig. 1). We detected 23 unique miRNA in RJM and 2 in RJC. In addition, there were 33 miRNAs differentially expressed in the two types of RJ (Fig. 2). In the up-regulated miRNAs (Table S1), amebantam, ame-mir-184, ame-mir-14, ame-mir-252 were the most abundant in RJM. The four miRNAs (ame-let-7, ame-mir-34, ame-mir-100, ame-mir-375) commonly found in other animal bodies or products (such as milk [39] or humans [40] and mouse [41]) were also present in RJ. Only two miRNAs (ame-mir-10, ame-mir-2944) showed higher expression in RJC. Consistent with Guo [15], we also found ame-mir-263, ame-mir-277, and amemir-283 in the two types of RJ. However, RJM in our study contained ame-mir-263b, which was absent in their study. This might be due to the fact that our RJ was obtained from 3 day old larvae (largely 4th instar larvae [42]) while they obtained RJ from 4-6th instar larvae.

Transcriptome Modification in A. mellifera due to Two Different RJs
Though A. mellifera and A. cerana might diverge from a common ancestor, they show differences in morphology, physiology and disease resistance. After fed with heterospecific royal jelly, A. mellifera showed many DGEs. Kucharski et al. [43] proposed that important elements of glutamatergic synapses are G-protein coupled metabotropic glutamate receptors (GPC mGluRs), which contribute to synaptic plasticity and development. According to their sequence similarity, transduction mechanism and pharma-cological profile, mGluRs are divided into three groups: group I (mGluR1 and mGluR5 receptors), group II (AmGluRA), and group III. The mGluR1 receptor links to phospholipase C, which causes phosphoinositide hydrolysis and release of calcium from intracellular stores. In mRJC, the expression level of mGluR1 (GB406151) was lower than that in mRJM, which might affect synaptic plasticity and development in the queen. Myosins [44] are one of three superfamilies of transporting motor proteins, which is involved in organelle formation, vesicle transportation, and cytoskeleton organization. In mRJC, Mhc1 (GB409843, a member of myosins) also decreased. The expression levels of InR-2 (GB725827), NLG-1 (GB724358), and trpgamma (GB410823) also decreased in mRJC. InR-2 is a member of insulin and insulinlike growth factor [45], which is linked to reproductive division of labor and foraging behavior [28]. NLG-1 and trpgamma are closely related to sensory input arising from environmental stimuli [46][47][48]. Four DEGs (GB552209, GB550937, GB410013, and GB409278) were involved in melanogenesis. They showed higher expression in mRJC and could explain prior studies showing darker coloration in mRJC [49], as was also the case in this study. A. cerana has been shown to have a higher sensitivity to odor than A. mellifera [50]. Seven genes (GB552209, GB406100, GB724316, GB551935, GB410013, and GB725569) involved in olfactory transduction were up-regulated in mRJC.

General Conclusions
We show for the first time that miRNAs in royal jelly are different between A. mellifera and A. cerana. Further, transcriptomes are modified as a result of bees being fed royal jelly of different species. Because a high proportion of the differentially expressed genes were target genes from miRNA, we speculate that the transcriptome modifications are partly caused by miRNA differences of the two species. This study provides the first evidence that miRNA in heterospecific royal jelly can modify gene expression in honey bees. Our results suggested that royal jelly from A. mellifera and A. cerana have different epigenetic effect on gene expression, although these two species are evolutionarily closely related.

Materials and Methods
Honey bee colonies (Apis mellifera and Apis cerana) were raised at the Honeybee Research Institute, Jiangxi Agricultural University, Nanchang, China (28.46uN, 115.49uE) by standard beekeeping techniques.

Differences in miRNA between Royal Jelly of Two Different Species of Honey Bees
Harvest of RJM and RJC. RJM and RJC were produced according to standard practices in China [51]. Briefly, the queen was confined inside a queen excluding cage. Queen cups with young larvae (one day old) were introduced into the colony and allowed to be fed by workers for 2 days. We first carefully removed 3 day old larvae by using either a grafting tool or a pair of forceps, then removed the royal jelly by using a spatula.  Table 4. Summary of known miRNA in RJM and RJC.
To further analyze the RNA secondary structure comprised with the matched Solexa reads, digital-quality sequences with perfect match or one mismatch were maintained. Genomic sequences of 100 nt were taken from these sequences, then the secondary structure was predicted and analyzed by RNAfold and MIREAP [39] under default settings. There are three criteria for candidate miRNA genes: (a) mature miRNAs are present in one arm of the hairpin precursors lacking large internal loops or bulges; (b) the secondary structures of the hairpins are steady, with the free energy of hybridization being lower than 220 kcal/mol; (c) hairpins are located in intergenic regions or introns [39]. Finally, these candidate miRNA reads were analyzed by miRBase database 13.0.
To compare miRNA expression levels between RJM and RJC, the reads of every miRNAs were subjected to the following   analysis: a miRNA was considered ''altered'' only if it had both: (a) 10 copies by Solexa sequencing in both RJM and RJC, and (b) a two-fold difference in copy numbers between RJM and RJC. GO assignments usually has three ontologies: biological process, cellular component and molecular function, which is used to classify the functions of almost all miRNA in our paper.

Transcriptome Modification in A. mellifera Fed Two Different RJs
Honey bee (A. mellifera) larvae were reared inside 24-cell tissue culture plates (Costar, NY, USA) inside an incubator (35uC, 7563% RH). Each cell was primed with 200 ml of freshly collected royal jelly, either from RJM or from RJC before a 1 day old larva was transferred into it. Larvae were transferred every 8 hrs to another plate with new food. For pupation, 6 day old larvae were transferred to 6-cell tissue culture plates lined with a piece of Kimwipe and kept in an incubator (35uC and 7863% RH) [52]. After adult emergence, we obtained one sample per treatment (mRJM or mRJC), each consisted of 10 adult bees (5 from each of the two colonies) and used their heads for global transcriptome analysis. The two samples were kept at 280uC until use. A total of 1,400 larvae were reared for this experiment.

Measurement of transcriptomes in mRJM and
mRJC. Total RNA was extracted with TRIzol regent (Invitrogen, USA) and treated with RNase-free DNase I (Takara Biotechnology, China). Poly(A) mRNA was separated by oligo-dT beads and then treated with the fragmentation buffer. By use of reverse transcriptase and random hexamer primers, the RNA fragments were transcribed into first-strand cDNA. Second strand cDNA synthesis was performed with DNA polymerase I and RnaseH. End-repair was done with T4 DNA polymerase, Klenow fragment, T4 Polynucleotide kinase. Ligation was accomplished with adapter or index adapter using T4 quick DNA ligase. Adaptor ligated fragments were selected according to size. Desired range of cDNA fragments were then excised from the gel. Finally, after validation of Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System, the cDNA library was sequenced by Illumina HiSeq2000.
By use of SOAP [53], a specific transcript with uniquely mapped reads was counted and their sequences assembled. Mapped reads was evaluated according to RPKM (Reads Per kb per Million reads) value of each transcript [38]. The transcript fold change was then calculated by the formula of log 2 (mRJM)_RPKM/(mRJC)_RPKM. The formula to calculate the probability of a specific gene being expressed equally between the two samples was defined as P(yDx)~( N2 N1 ) y (xzy)!
x!y!(1z N2 N1 ) xzyz1 ð Þ Where N1 and N2 indicate the total number of clean reads in mRJC and mRJM, respectively, and x and y indicate the mapped clean read counts of the transcript in each sample respectively. Then, the FDR (False Discovery Rate) method was applied to determine the threshold of the p-value in multiple tests. In this study, 'FDR,0.001' and the absolute value of log 2 Ratio .1 were used as the threshold to judge the significance of differentiated gene expression. We used the Blastall program to annotate the pathways of DEGs against the KEGG database.

Analysis of miRNA Targeted Genes
To identify possible target sequences of RJ, we used the RNA hybrid software and ftp.ncbi.nih.gov/genomes/Apis_mellifera/ RNA/rna.fa.gz, which provided us with a single predicted site of interaction with a minimum free energy. Table S1 The miRNAs expression analysis in RJM and RJC.