Honeybee ( Apis mellifera ) Maternal Effect Causes Alternation of DNA Methylation Regulating Queen Development

Maternal effects have been recognized in animals and plants for a long time (Bernardo 1996; Roach & Wulff, 1987). Maternal effects are that the phenotype of offspring is influenced by maternal phenotype rather than its own genotype (Van Dooren et al., 2016; Schwabl & Groothuis, 2019). Vertebrate females can adjust their investment in eggs to get higher quality offspring (Cunningham & Russell, 2000). Insects can also adaptively vary investment in their eggs (Passera, 1980). Honeybees are typical eusocial insect species, and the queen is the only female individual amongsterile workers (Winston, 1991). Previous studies showed that female larvae fed with richer nutritional diet develop into queens, whereas lower nutritional diet results in the development of workers (Haydak, 1970). The space in Abstract Queen-worker caste dimorphism is a typical trait for honeybees (Apis mellifera). We previously showed a maternal effect on caste differentiation and queen development, where queens emerged from queen-cell eggs (QE) had higher quality than queens developed from worker cell eggs (WE). In this study, newly-emerged queens were reared from QE, WE, and 2-day worker larvae (2L). The thorax size and DNA methylation levels of queens were measured. We found that queens emerging from QE had significantly larger thorax length and width than WE and 2L. Epigenetic analysis showed that QE/2L comparison had the most different methylated genes (DMGs, 612) followed by WE/2L (473), and QE/WE (371). Interestingly, a great number of DMGs (42) were in genes belonging to mTOR, MAPK, Wnt, Notch, Hedgehog, FoxO, and Hippo signaling pathways that are involved in regulating caste differentiation, reproduction and longevity. This study proved that honeybee maternal effect causes epigenetic alteration regulating caste differentiation and queen development. Sociobiology An international journal on social insects


Introduction
Maternal effects have been recognized in animals and plants for a long time (Bernardo 1996;Roach & Wulff, 1987). Maternal effects are that the phenotype of offspring is influenced by maternal phenotype rather than its own genotype (Van Dooren et al., 2016;Schwabl & Groothuis, 2019). Vertebrate females can adjust their investment in eggs to get higher quality offspring (Cunningham & Russell, 2000). Insects can also adaptively vary investment in their eggs (Passera, 1980). Honeybees are typical eusocial insect species, and the queen is the only female individual amongsterile workers (Winston, 1991). Previous studies showed that female larvae fed with richer nutritional diet develop into queens, whereas lower nutritional diet results in the development of workers (Haydak, 1970). The space in which larvae develop within the cell also contributes to the queen-worker differentiation (Shi et al., 2011). However, our previous study showed a maternal effect on honeybee caste differentiation, resulting in high quality queens. Queens lay larger eggs into queen cells compared to worker cells, which results in queens with heavier body, more ovarioles and different gene expression patterns (Wei et al., 2019). This maternal effect likely depends on the nutritional content of the fertilized queen cell egg, and not on genomic differences between queen-cell eggs and worker cell eggs. However, the underlying epigenetic mechanism for this maternal effect remains unclear.
More importantly, the quality of the queen has strong effects on the fitness of a colony (Winston, 1991;Gilley et al., 2003;Amiri et al., 2017). In recent years, honeybee colony loss has been frequently reported in public media, and the farming industry has been significantly impacted by the death of honeybees (Michael et al., 2009;Steinhauer et al., 2013;Antúnez et al., 2016). For example, scientists showed a high rate of honeybee winter colony loss over 13 years from 2006-2019 (Milius, 2019). The poor quality of queens was considered as one of main factors for honeybee winter loss (vanEngelsdorp et al., 2010;Delaney et al., 2011). Rudolf Steiner argued in 1923 that honeybees would become extinct within 100 years due to the weakening of queens in the commercial queen rearing industry (Thomas, 1998). Commercial beekeepers rear new queens by transplanting young worker larvae into queen cells, altering the larva's diet and depriving the maternal effect from the queen (Doolittle, 1888;Büchler et al., 2013). There are many concerns about the adverse consequences of this queen rearing technology for queen quality and colony health. Woyke (1971) showed that queens reared from young worker larvae were smaller and had a smaller spermatheca and fewer ovarioles compared to queens from worker eggs. Rangel et al. (2013) reported that queens reared from old worker larvae produced fewer workers and drone combs and less honey than colonies with queens reared from young worker larvae. He et al. (2017) and Wei et al. (2019) found that queens developed from queen cells are better than queens reared from worker cell eggs and young larvae. A recent study clearly showed that honeybee queens have an ability to alter egg size in response to both genetic and environmental factors (Amiri et al., 2020). Therefore, honeybee maternal effects may potentially influence the quality of the queen.
DNA methylation plays an important role in regulating honeybee queen development (Kucharski et al., 2008;Shi et al., 2011;Maleszka, 2008). Honeybee queens have lower genome-wide methylation levels than workers (Shi et al., 2013). The DNA methyltransferase Dnmt3 profound shifts honeybee reproductive status (Kucharski et al., 2008). Li-Byarlay et al. (2013) showed that knocking down DNA methyltransferase 3 (Dnmt3) changes gene alternative splicing in honeybees. We previously found that transplanting younger larvae from worker cells result in better queens with lower DNA methylation levels (He et al., 2017). Cardoso-Júnior et al. (2018) reported that DNA methylation altered queen lifespan by regulating the expression of vitellogenin gene. Therefore, we used honeybee eggs laid in queen cells to rear queens and compare DNA methylation among queens reared from queen-cell eggs (QE), worker-cell eggs (WE) and 2-day worker larvae (2L). Our results showed that maternal effect caused DNA methylation alternation in honeybee queens.

Insects
Three honey bee colonies (Apis mellifera) were used throughout this study. Virgin queens were sisters reared from the same mother queen in order to minimize differences in genetic background. Each colony had nine frames with approximately 30,000 workers and a single drone inseminated queen (SDI). These colonies were kept at the Honeybee Research Institute, Jiangxi Agricultural University, Nanchang,China (28.46 uN,115.49 uE). Queens were collected from the same colonies as in a previous study (Wei et al., 2019).

Queen rearing
Queens were caged in a plastic queen-cell frame for 6 hrs, and then transferred into a plastic worker-cell frame to lay worker-cell eggs following methods in Wei et al., (2019). Plastic queen cells frames were arranged horizontally. Generally, 20-50 eggs were harvested from queen cells after 6 hrs laying. Some worker-cell eggs were placed into plastic queen cells by transferring the base of each cell (Pan et al., 2013). The rest of the worker-cell eggs were kept in their native colonies until they reached 2d worker larvae stage. These 2d worker larvae were removed and placed into plastic queen cells. Cells from QE, WE and 2L were mixed together and placed into the same colony to rear queens. Newly emerged queens were harvested immediately for morphological measurement and DNA methylation sequencing.

Morphological measurements
Thoraxes of queens were collected and placed under a zoom-stereo microscope system (Panasonic Co., Ltd., accuracy: 0.01 mm). The width and length were measured according to the manufacturer's instructions. Each queen group had seven biological replicates, therefore, totally 21 queens were used for morphological measurements.

DNA methylation sequencing
Newly emerged queens were placed into liquid nitrogen, and heads and thoraxes were collected for total DNA extraction. Each sample contained tissue from two queens, and three biological replicates were conducted for each rearing condition. The QE group had two biological replicates since the sequencing of the third sample failed. Total DNA was extracted using Universal Genomic DNA Extraction Kit (TaKaRa, DV811A) following the manufacturer's protocol. Concentrations of all samples were measured and adjusted to the same level, and DNA samples (300 ng DNA for each sample) were used for DNA methylation sequencing by Illumina HiSeq™ 2500 (Illumina Inc., CA, USA). The detailed methods are listed in He et al. (2017). Briefly, DNA was sheared into 200-300 bp insert size targets using Covaris ultrasonicator (Life Technology) and then purified using AMPure XP beads and end repaired. A single 'A' nucleotide was added to the 3'ends of the blunt fragments followed by ligation to methylated adapter with a T overhang. Insert size targets (200-300 bp) were purified by 2% agarose gel electrophoresis. A ZYMO EZ DNA Methylation-GoldTM Kit (ZYMO, Irvine, CA, USA) was used to conduct bisulfite conversion. Bisulfite libraries were generated by PCR amplification and quantified by qPCR (Agilent QPCR NGS Library Quantification Kit). In total, eight libraries were prepared and sequenced. These processes were performed in Beijing Biomarker Technology Co., Ltd (Beijing, China).
Low quality reads (reads containing adapter sequences, reads containing unknown nucleotide "N" over 10 % and reads with a quality value lower than 10 occupying more than 50% of the whole read) from eight bisulfite libraries were filtered. These filtered genomic fragments were then mapped to the honeybee genome (Apis mellifera. Amel 4.5), transformed into bisulfite-converted versions of the sequences (C-to-T and G-to-A) and then assigned to a digital index using bowtie2 according to Langmead and Salzberg (2012). Methylation sites were predicted using a bismark methylation extractor (Krueger & Andrews, 2011). Only uniquely mapped reads (clean reads) were retained. The ratio of C to CT was used to measure methylation levels.
The Bismark package was used to test 5mC (same to Krueger & Andrews, 2011) and different methylation regions (DMRs) were identified using Bisulfighter (Yutaka et al., 2014), according to our previous study (He et al., 2017). DMR related genes (DMGs) were annotated by mapping the target regions [gene body region (from transcription start sites to transcription end sites) or promoter region (upstream 2kb from the transcription start sites)] to the honeybee genome (Amel 4.5), since DNA methylation data was used to compare with our previous RNA-Seq data (Wei et al., 2019).

Data analysis
For queen thorax length and width, data were analyzed using an analysis of variance (ANOVA) followed by Fisher's PLSD test in Statview 5.01 package (SAS, USA), where p < 0.05 was considered as significantly different.

QE were larger than WE and 2L
The thorax length and width of QE were significantly higher than that of WE and 2L. In addition, the thorax size of WE was significantly higher than 2L (Fig 1, p < 0.05).

Quality of methylation sequencing
In total, 43.44 G clean bases were detected from eight libraries. Q30 values of samples were over 85 % and the average ratio of bisulfite conversion reached to 99 % (Table S1). These values indicate a considerably high sequencing quality.

DMRs and DMGs contribute to queen development
There were hundreds of DMRs among three queen types (Fig 2 and S1). DMRs (CG type) were mapped to all 16 chromosomes. Most chromosomes, such as Chr4, Chr5, Chr10 and Chr14, had more DMRs in QE/2L when comparing to WE/2L and QE/WE (Fig 2). The DMRs in CHH type were less than the CG type but showed a similar pattern (Fig S1). There were no DMRs of CHG type could map to the chromosomes.    QE/2L had more DMGs (612) than WE/2L (473) and QE/WE (371) comparisons (Fig 3, Table S2-4). These DMGs were enriched in 72 KEGG pathways which belong to four KEGG categories including metabolism, gene information processing, environmental information processing, cellular processes and organismal systems (Fig 4). Moreover, 42 DMGs enriched in environmental information processing were genes belonging to important signaling pathways, such as mTOR, MAPK, Notch, and Wnt signaling (Fig 5), which are involved in regulating queen development and reproduction (Chen et al., 2012;Chen et al., 2017). Similarly, QE/2L comparison had more DMGs (21) of these 42 genes than WE/2L (20) and QE/WE (17) comparisons (Fig 5).

DNA methylation had a low relationship with previous RNA-Seq data
By comparing the DNA methylation data to our previous RNA-Seq data from Wei et al., (2019), only 2 (GB41428 and GB46222), 1 (GB41428) and 0 DMGs were mapped to DEGs in QE/2L, WE/2L and QE/WE comparisons respectively (Fig S2). We also compared 35 immunity and development DEGs with our DNA methylation data, and the result (Fig 6) showed that the differences of gene expression and DNA methylation of these genes both increased as the age of transplant of the worker larvae increased. However, there wasn't a clear negative correlation between gene expression and DNA methylation in these genes (Fig 6).

Discussion
Environmental factors, such as larval diet and maternal effect, contribute towards honeybee caste differentiation and queen development (Haydak, 1970;Wei et al., 2019). In this study, queens developed from queen-cells had significantly larger thorax size than those from worker eggs and young worker larvae (Fig 1). Many studies have shown that queen weight and body size are strongly correlated with queen ovariole number which influences queen fecundity and quality (Borodacheva, 1973;Bilash et al., 1983;Huang & Zhi, 1985;He et al., 2017;Wei et al., 2019). Consequently, these results from the present study are consistent with our previous study (Wei et al., 2019) where we showed that QE had significantly higher weight and ovariole number than both WE and 2L. We conclude there is a strong maternal effect on honeybee queen development and potentially contributes to the queen quality.
Previous studies showed that queens have lower global DNA methylation level than workers, and DNA methylated genes such as genes involved in mTOR pathways can control queen-worker dimorphism (Kucharski et al., 2008;Chen et al., 2012;Shi et al., 2013). Here, we found that QE/2L had the most DMRs and DMGs, followed by WE/2L and QE/WE (Fig 2 and 3). In addition, 42 DMGs were involved in 11 signaling pathways such as mTOR, MAPK, Wnt and Notch pathways (Fig 4 and 5) that are known to regulate honeybee caste differentiation and queen development (Chen et al., 2012;Chen et al., 2017). Blue indicates significant downregulation in each comparison, red indicates upregulation and green indicates no difference. Left side is the gene ID and KEGG pathways. Some genes were involved in two or more pathways; therefore they were marked with different color lines. These DMGs induced by maternal effect are strongly related to queen-worker differentiation, reflecting that honeybee maternal effect could cause various epigenetic alteration in caste differentiation and queen development.
How is DNA methylation regulated by maternal effects? Different diet during larval stage alters DNA methylation (Wang et al., 2006;Kucharski et al., 2008;Shi et al., 2011;Foret et al., 2012;Maleszka, 2008). We previously showed that queens deposited larger eggs into queen cells (Wei et al., 2019). During development, eggs absorb nutrients such as vitellin from the ovary epidermal cells for 21 days (Torres, 1980;Li & Zhang, 2017). Larger eggs found in queen cells may contain more nutrients than eggs in worker cells, therefore, the nutrients available to the developing embryo may alter DNA methylation and promote queen development. However, the kind of nutrients in the eggs and how queens control egg size requires further investigations.
Moreover, DNA methylation plays an important role in soft inheritance and animal evolution (Bird, 2002;Dickins & Rahman, 2012;Klironomos, 2013). Maternal effect also dramatically contributes to the process of animal evolution (Galloway et al., 2009;Marshall & Uller, 2007). The commercial queen rearing technology started in the 19 th century (Doolittle, 1888) artificially transplants young worker larvae into queen cells to rear queens rather uses queen cell eggs. The widely use of commercial queen rearing technology for over 100 years continuously deprives the investment from mother queens (namely maternal effect) and may consistently weakens the queen and honeybee colony via altering DNA methylation. As a proximal remedy, rearing queens from queen-cell eggs may be a good strategy to yield high quality queens and healthy colonies.
Similar to our previous study (He et al., 2017), our DMGs had a low relationship with previous RNA-Seq results (Wei et al., 2019). Only 3 of 1456 DMGs were mapped to our previous DEGs data ( Fig S2). Evidence showed that honeybee DNA methylation has an association with alternative splicing and gene duplication (Elango et al., 2009;Dyson & Goodisman, 2020), but there isn't a completely direct link between DNA methylation and gene expression. The present study and our previous study showed a same pattern that QE/2L comparison had more DMGs and DEGs than QE/ WE and WE/2L comparisons (see Fig 3 from this study and Fig 3 from Wei et al., 2019). Many DMGs and DEGs were involved in queen-worker differentiation, body development and immunity etc., however they did not show a clear negative or positive relationship with our DNA methylation data (Fig 6). How DNA methylation regulating gene expression remains essentially unknown in honeybees. DNA methylation has been shown to associate with chromosome structure and histone modifications (Hunt et al., 2013;Dmitrijeva et al., 2018). Perhaps the alternative splicing, chromosome structure, histone modification and other undiscovered factors jointly contribute to the regulation on honeybee gene expression by DNA methylation, which needs further investigations.
In conclusion, this study firstly indicates that honeybee maternal effect causes DNA methylation changes in queen rearing and potentially contributes to the queen development and colony health, due to the function of DNA methylation in soft inheritance and animal evolution. Since the poor quality has been frequently reported as a main factor for colony losses (vanEngelsdorp et al., 2010;Delaney et al., 2011), therefore, the maternal effect should be deeply considered and used in commercial queen rearing.
We have declared that no conflict of interests exists.