1 Introduction

The honey bee presents one of the most impressive documented cases of developmental plasticity (Shell and Rehan 2018; Wheeler 1986). There are no genetic differences between the queen and worker castes. Both castes develop from divergent developmental pathways (Winston 1991). Recent and rapid progress in comparative functional genomics has revealed how the developmental and nutritional environment in the early larval stages establishes a cascade of epigenomic changes that alter patterns of both gene methylation and gene expression that ultimately result in different developmental systems and two radically different phenotypes (He et al. 2017; Kucharski et al. 2008; Yin et al. 2018).

In nature, the developmental environment is established by both the queen and the workers. Workers create just a few large queens cells, and many small hexagonal worker cells. The queen lays eggs in both cell types, but she lays slightly larger eggs in queen cells ensuring that even prior to hatching a queen-destined bee is provisioned with more food (Wei et al. 2019). Workers supply queen cells with an abundance of royal jelly, and provide worker cells with far less worker jelly. The two food types differ in nutrient content and provide very different nutrition for the two types of developing larvae (Wang et al. 2016).

Commercial queen-rearing practices in contemporary apiculture exploit this plastic developmental system by creating artificial queen cells into which eggs or young larvae from worker cells are transplanted (Doolittle 1888). If these artificial queen cells are placed in a recently queenless colony, they will be provisioned and reared by the workers as queens. Honey bee eggs are very fragile and transplanting eggs is technically difficult. For this reason, it is common for queen breeders to transfer larvae from worker cells. Transplanting larvae from worker to queen cells up to the 3rd instar stage will result in a functional queen phenotype (Doolittle 1888). While the queen developmental pathway is robust enough to tolerate this kind of manipulation, diverting larvae from a worker developmental pathway to a queen developmental pathway in the mid larval stage is not without consequence.

Several authors have now documented how queens reared from transplanted worker larvae are smaller, weighed less, have smaller spermathecae and fewer ovarioles than queens reared naturally or from transplanted eggs (He et al. 2017; Wei et al. 2019; Woyke 1971). The differences get more pronounced as the age of the transplanted larva increases. More recent analyses have documented differences in gene expression in the developing queens depending on whether they were produced from transplanted larvae or eggs (He et al. 2017; Yin et al. 2018), and He et al. reported even changes to the methylation of some genes (He et al. 2017).

The gene methylation changes are of particular interest since differential methylation is considered a key part of the mechanism that both enable the genome to be sensitive to the nutritional and developmental environment, and that establishes the different developmental pathways (Kucharski et al. 2008). Here, we examined the changes to the honey bee methylome resulting from rearing queens from eggs and larvae of different ages. To assess the impacts of different methods of queen rearing on the queen epigenome, our analysis considered changes across the whole genome, and also focused on genes that may be functionally related to caste differentiation processes that differ between workers and queens.

2 Materials and methods

2.1 Queen rearing and sampling

Colonies of the Western honeybee, Apis mellifera, were maintained at the Honeybee Research Institute, Jiangxi Agricultural University, Nanchang, China (28.46 uN, 115.49 uE). Experiments were performed under the supervision of the Animal Care and Use Committee of Jiangxi Agricultural University, China.

In this study, we compared four groups of queens: queens raised from worker eggs transplanted to queen cells on the 2nd day after laying, and queens raised from worker larvae transplanted to queen cells on the 4th, 5th, and 6th day after laying.

To produce these queens, a naturally mated laying queen was caged for 6 h on a frame of plastic worker cells to lay (Pan et al. 2013). This frame was moved to an incubator and 30–50 eggs or larvae were transferred to plastic queen cells on the 2nd day (egg queens raised from eggs, E), 4th day (1-day-old larvae queens, L1), 5th day (2-day-old larvae queens, L2), and 6th day after laying (3-day-old larvae queens, L3).

The queen cells were placed in racks in two queenless colonies to be reared. Half of each queen rearing group was assigned to each colony. Once sealed, queen cells were numbered, the length of each queen cell was measured and they were then relocated to a dark incubator (35 °C, 80%) to emerge. From the 15th day post-laying, the cells were checked every 2 h to see if any queens had emerged, and hourly after the first queen emerged. We measured the weight of royal jelly remaining in the queen cell after queen emerged.

The first four queens in each rearing group to emerge were flash frozen in liquid nitrogen and stored in a – 80 °C freezer for later genomic methylation analysis. The remaining queens were transferred to queen cages, which were placed in queenless colonies for 4–5 days and then flash frozen in liquid nitrogen and stored in a – 80 °C freezer for later ovariole scoring.

2.2 Scoring ovariole number

Sampled queens were thawed to room temperature and the ovaries dissected from the abdomen. For dehydration and fixing, we used an Automatic dehydrator (Leica, TP1020). First, ovaries were fixed in 4% paraformaldehyde fix solution (BBI Life Sciences) for 12 h. They were then dehydrated in a graded ethanol series (70–100%). Following dehydration, ovaries were placed in a 1:1 absolute ethanol/xylene mixture for 30 min, then 100% xylene for 10 min, followed by immersion in fresh xylene for 5 min. Finally, they were transferred to a 1:1 xylene-paraffin mixture for 30 min, and embedded in paraffin wax using a heated paraffin embedding station (Leica).

Once blocked in paraffin wax, 5–7 μm sections of the ovary were cut with a Leica RM 2245 microtome. Sections were placed on histological slides (Autostainer XL), stained with HE Staining Kit (BOSTER AR1180), mounted with neutral balsam mounting medium (BBI Life Sciences), and covered with a coverslip. Slides were imaged and photographed using a × 100 transmission light microscope (OLYMPUS, DP80) (Figure S1).

We counted the number of ovarioles in the left ovary by identifying sections in which the ovarioles were very clear and counting sections (fiji-win64 software, Figure S1) until we found at least two giving exactly the same number of ovarioles (Gan et al. 2012). We used a one-way ANOVA followed with Fisher’ s PLSD test to analyze differences between queen types.

2.3 Genome-wide methylation analysis

Four queens in each group were used for genome-wide methylation testing. Brains, ovaries and thorax of the queens were dissected over ice and pooled as one sample for genomic DNA extraction.

Genomic DNA was extracted (StarSpin Animal DNA Kit, GenStar). DNA samples were then sent for whole-genome bisulfite sequencing analysis (Novogene Bioinformatics Technology Co., Ltd./www.novogene.cn) using the method summarized below.

2.4 Library preparation and quantification

One hundred nanograms of genomic DNA were spiked with 0.5 ng lambda DNA and sonicated using a Covaris S220 DNA Sequencing/gene analyzer to break the genomic DNA into short 200–300 base pair fragments. Fragmented DNA was treated with bisulfite (EZ DNA Methylation–GoldTM Kit, Zymo Research), and the bisulfite-converted sequences were then processed using an Accel–NGS Methyl–Seq DNA Library Kit to add to each fragment a truncated adapter sequence, including an index sequence.

Library DNA concentration was quantified by a Qubit® 2.0 Flurometer (Life Technologies, CA, USA) and quantitative PCR. The insert size was assayed on an Agilent Bioanalyzer 2100 system.

The prepared genomic DNA library samples were sequenced (Illumina Hiseq XTen) and 125–150 bp paired end-reads were generated. Image analysis and base identification were performed with the Illumina CASAVA pipeline.

The quality of the raw reads was assessed with FastQC (fastqc_v0.11.5). Read sequences produced by the Illumina pipleline in FASTQ format were pre-processed through Trimmomatic (Trimmomatic-0.36) software use the parameter (SLIDINGWINDOW:4:15, LEADING:3, TRAILING:3, ILLUMINACLIP:adapter. fa:2, 30:10, MINLEN:36). Reads that passed these filtering steps were counted as clean reads and all subsequent analyses were performed on clean reads only. For all samples, two samples with low mapping values did not meet our requirements. Consequently, we lost one E sample and one L3 sample. The average number of clean reads was 41, 485,779, with 10.89 G clean base sequences (the two eliminated samples are not included in these statistics). The average Phred scores Q30% (Ewing and Green 1998) was 92.5% (Table SI).

2.5 Mapping reads to the Apis mellifera reference genome

Bismark software (version 0.16.3) was used to align bisulfite-treated library reads to the reference honey bee genome (Amel_HAv3.1 (GCF_003254395.2)) (Krueger and Andrews 2011). For alignment, both the reference genome and library reads were transformed into bisulfite-converted versions (C-to-T and G-to-A) and then given a digital index with bowtie2, which included data on the sequences, and sample (Langmead and Salzberg 2012).

Aligned sequence reads were then compared with the normal (nonbisulfite-converted) genomic sequence so that the methylation state of cytosines could be inferred. Reads that aligned to the same regions of the genome were regarded as duplicates. Sequencing depth and coverage were calculated assessing number of overlapping reads relative to number of duplicate reads.

The conversion rate was obtained by aligning to the lambda spike. The average bisulfite conversion rate was 99.7% (Table SI). The average site coverage rate was 26.74 (Table SII).

2.6 Methylation analysis

To identify the level of methylation at each site, we modeled the count of methylated cytosines (mC) at a site as a binomial (Bin) random variable with methylation rate r: mC~Bin (mC + umC*r) (Wang et al. 2019).

To calculate the methylation level of a sequence, we divided the sequence into multiple bins, of 10 kb and measured the sum of methylated and unmethylated read counts in each bin. Methylation level (ML) for each bin or C site shows the fraction of methylated Cs, and is defined as follows:

$$ \mathrm{ML}\ \left(\mathrm{C}\right)=\left(\mathrm{reads}\ \mathrm{mC}\right)/\left(\mathrm{reads}\ \left(\mathrm{mC}\right)+\mathrm{reads}\ \left(\mathrm{C}\right)\right). $$

Calculated ML was further corrected with the bisulfite non-conversion rate according to Listers’ method (Lister et al. 2013). Given the bisulfite non-conversion rate r, the corrected ML was estimated as: ML (corrected) = (ML - r) / (1 - r).

When comparing samples, differentially methylated regions (DMRs) were identified using DSS software (Feng et al. 2014; Park and Wu 2016; Wu et al. 2015) with the parameters: smoothing = TRUE, smoothing.span = 200, delta = 0, p.threshold = 1e-05, minlen = 50, minCG = 3, dis.merge = 100, pct.sig = 0.5. These parameters mean regions with the minimal length of 50 bp, more than 3 CG sites, and more than 50% of differentially methylated loci in the regions with P < 1 × 10−05 were considered as DMRs. DSS is an common used R library performing differential analysis for count-based sequencing data (MacKay et al. 2019; Mendizabal et al. 2019; Wang et al. 2019). The core of DSS is a new dispersion shrinkage method for estimating the dispersion parameter from Gamma-Poisson or Beta-Binomial distributions. Adjacent DMRs were combined when the distance between two DMRs was less than 100 bp. We analyzed the DMRs in all contrasts (L1 vs. E, L2 vs. E, L3 vs. E, L2 vs. L1, L3 vs. L2, L3 vs. L1).

We related DMRs to genes if the DMR overlapped at all with the gene body or promoter region. The KEGG database classifies genes according to high-level functions (http://www.genome.jp/kegg/). We used KOBAS software (Mao et al. 2005) to test the statistical enrichment of DMGs to different KEGG pathways. We analyzed the 20 most enriched KEGG pathways to explore the most differ pathway between the comparisons of L1–3 with E.

3 Results

3.1 Effect of queen rearing method on queen morphology

As the age of the larva transferred to the queen cell increased, queen cell length (Fig. 1a, ANOVA: F = 84.240, DF = 3, P < 0.0001), remaining royal jelly (Fig. 1b, ANOVA: F = 95.408, DF = 3, P < 0.0001), and queen weight (Fig. 1c, ANOVA: F = 14.772, DF = 3, P < 0.001) all decreased significantly. There was no significant difference in the number of ovarioles between experimental groups (Fig. 1d, ANOVA: F = 2.918, DF = 3, P = 0.0607).

Figure 1.
figure 1

Queen cell length (a), remaining royal jelly (b), weight (c), and number of ovarioles (d) of the four queen groups (E, L1, L2, and L3). Bars show mean ± SEM. Sample size of each queen group shown in each bar. One-way ANOVA was performed and letters above bars identify groups that did and did not differ significantly with each other (P < 0.05, ANOVA test followed with Fisher’ s PLSD test).

Full size image

3.2 Effect of queen rearing method on genome methylation

DMGs in all comparisons (L1 vs. E, L2 vs. E, L3 vs. E, L2 vs. L1, L3 vs. L2, L3 vs. L1) were shown in Table SIII. When considering the DMGs in different functional classes, we focused on genes identified as related to the function of reproduction or longevity (Corona and Robinson 2006; He et al. 2017; Yin et al. 2018), immunity (Barribeau et al. 2015; Boutin et al. 2015; He et al. 2017), metabolism (He et al. 2017), and DMGs identified in our own data related to these groups from KEGG pathway analyses. If genes were related to multiple functions, they were only counted once. Our classification criteria was reproduction or longevity > immunity > metabolism (Table SIV). Then, we mainly focused on the DMGs of L1–3 vs. E for further study.

The 20 most enriched KEGG pathways indicated that for all comparisons of queens reared from larvae with E queens most DMGs were involved in metabolic pathway, and as the age of the transferred larvae increased, the number of DMGs and degree of difference between treatment groups increased (Fig. 2). In these three comparisons, there were 79 DMGs (22, 34, and 52 DMGs in the comparisons of L1 with E, L2 with E, L3 with E respectively) related to the function of metabolism (Table SV). Of these, some genes were involved in multiple important pathways like caste differentiation, and immunity (Lyko et al. 2010).

Figure 2.
figure 2

Scatter plot of the 20 most enriched KEGG pathways enriched in DMGs in comparisons of L1 vs. E, L2 vs. E, and L3 vs. E. The size of the dots indicates the number of DMGs in that functional pathway. The color of the dots indicates q value (corrected p value). The rich factor indicates the relative enrichment of DMGs in the KEGG pathway (number of DMGs in the KEGG pathway/total number of annotated genes in the same KEGG pathway). The larger the rich factor, the greater the degree of enrichment.

Full size image

We analyzed the methylation differences at CG sites in the exons of all the DMGs (L1–L3 vs. E) functions related to reproduction or longevity, immunity, body development, or metabolism, and found 55 genes that showed a uniform increasing or decreasing trend with the age of larval transfer. This analysis revealed that then comparing queens reared from larvae (L1–L3) with E queens, the number of DMGs (Fig. 3), and the ratio of methylation differences between the compared sequences (Fig. 3) increased with the age of larval transfer to queen cells.

Figure 3.
figure 3

Fifty-five focal genes, related to functions in reproduction, longevity, immunity, learning and memory, and body development or metabolism, were selected for further analysis. Relative methylation level (given as the value in the box and indicated by the color of each cell) was calculated by comparing the proportions of methylated reads at each CG site in all exons of the gene between larval reared queen groups (L1–L3) compared with E. Green and blue indicates hypo-methylated genes in groups compared with E. For the scare bar, red, purple, and black indicates hyper-methylated genes, and yellow indicates no difference. Black borders indicate that there was at least one exon in this gene that was significantly differentially methylated in the comparisons. Boxes with no black border indicates that no exon in this gene was significantly differentially methylated in the comparisons. A differentially methylated exon is part of the differentially methylated region (DMR). More detailed information is given in Table SVI.

Full size image

4 Discussion

While it appears that functional queens can be produced by transplanting larvae from worker cells up to 3 days old, our study revealed significant consequences of this for queen development and the epigenome of the adult queen.

Queens raised from older worker larvae weighed less (Fig. 1c), had a tendency but not significantly fewer ovarioles than queens reared from eggs (Fig. 1d) which might be due to low sample size. This confirms reports from earlier studies (Woyke 1971), and is concerning since ovariole number will limit the queens fecundity (Bouletreau-Merle 1978; Woyke 1971), which can reduce colony growth and productivity.

It was also notable in this study that the queen cells were smaller for queens reared from larvae than queens reared from eggs (Fig. 1b) and on emergence contained less residual royal jelly. This demonstrates that the developmental and nutritional environment of queens reared from transplanted larvae is not the same as queens reared from transplanted eggs. We propose that since queens raised from larvae have less time in the queen cell before pupation than E queens there is less time for workers to provision them with an abundance of food, and they have less time prior to pupation to consume the rich royal jelly diet with its distinctive nutrition. This may explain both the drop in weight of the emerging adult queen and the smaller size of the sealed queen cell.

In this study, the parental queen was naturally mated, which is the closest scenario to natural conditions. Thus the offspring queens could carry a different set of chromosome from the parental drones. Therefore, the variance within each treatment group may be increased, but this could reduce any false positive loci. However, we have to acknowledge that allowing the introduction of multiple genotypes might have affected the results because the exactly intragenic DNA methylation remains unclear in social insects’ genomes.

In our data, the degree of difference of DMRs in some genes increased with the age at which the worker larva was transplanted, when compared to queens reared from transplanted eggs (Figs. 2 and 3).

Our analyses focused on the functions of the identified DMGs. Many of the DMGs have known functions in immunity, reproduction, and longevity. Some of these 55 genes may play an important role in caste differentiation or the development of queens and deserve further study. These are processes that might be compromised by immune stress during development (Butler and McGraw 2011; Searcy et al. 2004). Moreover, we found some DMGs identified in our comparisons of queens from the E group and the L groups were also identified in comparison of queens and workers (Lyko et al. 2010).

Lyko et al. (2010) found 561 DMGs between brains of queens and workers, out of which 36 genes were also found in the DMGs when comparing the L1 and E groups. Thirty-nine genes occurred the DMG comparison of L2 with E. Sixty-nine DMGs from Lyko et al. (2010) were also identified in the DMG comparison of L3 with E (Table SVI). Furthermore, some genes with functions related to reproduction, longevity, immunity, body development, metabolism, or others identified as differentially methylated in the comparisons in the brains of queens and workers by Lyko et al. 2010 were also found in our comparison of queens from different rearing treatments (Table SVI). In this case, it is especially concerning because these traits might compromise the health, fecundity, and longevity of the queen, and thereby the vigor of the whole colony. However, since DNA methylation might not always be associated to change in gene expression, the role of these differentially methylated genes remains to be further studied.

These characteristics also distinguish queens and workers, and our epigenomic analyses revealed that for these key traits queens reared from worker larvae might be more worker-like than queens reared from eggs. In effect, our epigenomic analyses suggest that the adult phenotype of queens reared from larvae is partially intercaste and more worker-like than naturally reared queens.

Further supporting this interpretation, KEGG analyses revealed that the most common functional class for the DMGs in our comparisons of queen of different rearing types was metabolic processes (Fig. 2). Of the 79 DMGs that were related to the function of metabolism, some genes were involved in other important functions like caste differentiation, and immunity (Lyko et al. 2010). For example, the gene Cat (LOC443552) was involved in peroxisome, glyoxylate, dicarboxylate metabolism, and the FoxO signaling pathway. The gene Plc (LOC408791) was involved in the metabolic pathways, phototransduction, inositol phosphate metabolism, and the Wnt signaling pathway. The FoxO signaling pathway and Wnt signaling pathway have been related to caste differentiation (Yin et al. 2018). Moreover, some DMGs we identified in the metabolic pathway functional group were also identified by Lyko et al. (2010) in their comparison of the metyhylation of brains of queens and workers (Table SVI). Ndufs1 (LOC726035) and Ndufs2 (LOC413891) were all involved in metabolic pathways, and these two genes have also been related to immune functions (Boutin et al. 2015).

Metabolic processes normally differ strongly between workers and queens (Corona et al. 1999; Evans and Wheeler 1999; Severson et al. 1989). In comparisons of gene expression differences between workers and queens, the greatest and most abundant differences in gene expression are also in genes related to metabolic processes (Corona et al. 1999; Evans and Wheeler 1999; Severson et al. 1989). The above suggests that queens reared from larvae are shifted away from the metabolic program of queens reared from transplanted eggs.

To summarize, our analyses revealed the epigenomic mechanisms that enable the differentiation of the worker and queen phenotypes are compromised if queens are reared from transplanted worker larvae, with the issue becoming more severe as the transplanted larvae become older. This is significant because the most common commercial practice in apiculture is to rear queens from larvae rather than eggs because larvae are far easier than eggs to transplant. Our data argue, however, that rearing queens from eggs gives a better developmental outcome and a more fecund and robust developmental outcome queen than rearing queens from larvae.