Current Issue

The Western honeybee, Apis mellifera L., plays a key pollinator role in natural and agricultural ecosystems (Leonhardt et al., 2013), with an estimated commercial value of over $15 billion per year in the United States (Walsh, 2013). In the Republic of Korea, the economic value of bee pollinating fruit and vegetable agriculture was estimated to be approximately $5 billion (Jung, 2008). In addition to its economic value, the honeybee has been extensively studied as a model insect species due to its highly developed sociability, specialization in division of labor, and ability to manage colonies (Lee and Kim, 2017). Honeybee colonies consist of eggs, larvae, pupae, and adults, and the adults are divided into drones, female workers, and queens (Winston, 1991). The queen bee only consumes royal jelly during the larval stage, and its average development period is usually 16 days which is shorter than that of worker bees (18 days). The queen bee is responsible for oviposition; it can lay as many as 1,500 eggs per day for an average of three to four years of life (Winston, 1991). Unlike the queen bee, which receives royal jelly for the entire larval period, worker bees are developed when larvae are fed regular honey after feeding royal jelly for three days (Asencot and Lensky, 1988). Following emergence, age-related division of labor is based on a form of behavioral development by workers known as age polyethism (Robinson and Huang, 1998). During the first 2~3 weeks of adult life, the nurse bees are typically responsible for in-hive tasks such as brood care (nursing) and hive maintenance. Followed by a transition to tasks outside the hive, forager bees predominantly collect honey and pollen 2~3 weeks before death (Dukas, 2008).

Interestingly, the polyethism of honeybees is not rigid; instead, the honeybee can respond to changes in environment and colony conditions and flexibly alter their typical patterns of age polyethism (Robinson and Huang, 1998). For example, when few nurse bees have remained in the colony, foragers expand their duties to the tasks usually performed by nurses, while in the condition that the number of the forager is limited, nurse bees perform the task of foragers (Robinson and Huang, 1998). Considering the bee colony composing different developmental stages of the honeybee, development of queen/worker by royal jelly feeding, division of nurse/forager, and flexibility of age polyethism, determination of genes putatively involved in bee physiology and investigation of gene expression are essential to extend our knowledge to development and evolution of sociality in insects.

To date, quantitative real-time PCR (qRT-PCR) has been widely suggested as the primary technique for quantifying gene expression levels because of its relatively high speed, sensitivity, reproducibility, and accuracy (Ling and Salvaterra, 2011; Zhai et al., 2014). However, normalization should be performed with reliable reference genes stably expressing across different samples to determine target gene expression levels accurately. Therefore, suitable reference genes should be selected (Wong and Medrano, 2005; Ling and Salvaterra, 2011).

Therefore, in order to study honeybee physiology at different developmental stages, we investigated the expression stabilities of housekeeping genes in different developmental stages (egg, 1^st instar larva, 3^rd instar larva, 5^th instar larva, pupa, nurse, and forager) to suggest appropriate reference genes for gene expression study using qRT-PCR. In the present study, we selected five candidate reference genes, including ribosomal protein S5 (RPS5), ribosomal protein S18 (RPS18), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), ADP-ribosylation factor 1 (ARF1), and Ras-related protein Rab-1A (RAD1a), which have been previously suggested as reference genes in honeybee studies (Moon et al., 2018; Jeon et al., 2020; Kim et al., 2021). However, these genes have not been evaluated in the differenet developmental stages in honeybee, therefore, the stabilities of these five genes across developmental stages were evaluated with C_q distribution analysis and three programs (geNorm, NormFinder, and BestKeeper).

Considering that the honeybee is a good model insect for molecular physiological studies of social development (Consortium, 2006), determining the expression level of genes putatively involved in honeybee physiology is essential. In addition, appropriate reference genes stably expressing across the different sample conditions should be selected for accurate normalization of the gene of interest with qRT-PCR (Lourenço et al., 2008; Scharlaken et al., 2008; Reim et al., 2013). Therefore, several studies have been conducted to evaluate qRT-PCR reference genes in honeybees (Lourenço et al., 2008; Scharlaken et al., 2008; Reim et al., 2013; Moon et al., 2018; Jeon et al., 2020; Kim et al., 2021).

In this study, in order to select suitable reference genes in different developmental stages of the honeybee, we evaluated the expression stability of five candidate reference genes (RPS5, RPS18, GAPDH, ARF1, and RAD1a) across seven developmental stages (egg, 1^st instar larvae, 3^rd instar larvae, 5^th instar larvae, pupa, nurse, and forager) by analyzing C_q distribution of the genes and three analysis algorithms programs (geNorm, NormFinder, and BestKeeper). Although the selected five candidates have not been evaluated in different developmental stages, they were previously suggested as the appropriate reference genes in the honeybee. According to the previous studies, RPS5, RPS18, and GAPDH represented optimal reference genes in various tissue of honeybee collected across seasons (Jeon et al., 2020). In addition, RPS18 and GAPDH were suggested to be employed as suitable reference genes for qRT-PCR-based determination of seasonal and labor-specific gene expression profiles (Moon et al., 2018). Furthermore, RAD1a and ARF1 were suggested to be the reliable reference genes in different tissues and two adult stages (nurse and forager) of honeybees exposed to various pesticides (Kim et al., 2021). The selection of these five candidates was reliable because they mostly show stable expression regarding their cellular physiological function. The ribosomal proteins (RPS5 and RPS18) and metabolic enzyme (GAPDH) are essentail housekeeping genes involved in protein expression (Zhou et al., 2015) and glycolysis (Tunio et al., 2010), respectively. Two genes (ARF1 and RAD1a) belong to the Ras superfamily and are responsible for intracellular and membrane transport (Chavrier and Goud, 1999).

According to the results analyzed by geNorm, NormFinder, and BestKeeper, the three algorithms revealed slightly different stability values of the five candidate reference genes, as observed in previous studies (Moon et al., 2018; Wang et al., 2019; Jeon et al., 2020; Kim et al., 2020); therefore, the combined use of results analyzed by all programs would ensure more credible suggestion for optimal reference genes (Wang et al., 2019; Jeon et al., 2020). In the analysis of NormFinder, the NormFinder program does not provide an appropriate cutoff value, but recent studies suggested that a suitable cutoff value of gene expression stability is ≤0.15 (McMillan and Pereg, 2014; Julian et al., 2016). Based on these criteria, all five genes met the values, implying that RPS5, RPS18, GAPDH, ARF1, and RAD1a would be possibly suggested to be applied as the qRT-PCR reference genes (Fig. 4). However, considering M≤0.5 as a criterion for selection of an appropriate reference gene in the geNorm analysis (Hellemans et al., 2007; Liu et al., 2014), geNorm resulted in three genes (RPS5, RAD1a, and ARF1) as optimal reference genes (Fig. 3A). Although the rankings of the stability values of these three genes were different between geNorm and NormFinder analysis, RPS5, RAD1a, and ARF1 were also ranked in the high stable genes when compared with GAPDH and RPS18 in NormFinder (Fig. 4). Furthermore, BestKeeper revealed that RAD1a, ARF1, and RPS5 were the first to third stable genes, respectively, as judged by their lower SD values than RPS18 and GAPDH, although SD values of all five genes were ≤1.0, which is the cutoff line in BestKeeper analysis (Fig. 5) (Chechi et al., 2012). Likewise, C_q distribution analysis demonstrated that RAD1a, ARF1, and RPS5 exhibited relatively lower CV values than RPS18 and GAPDH (Fig. 2B), suggesting that RAD1a, ARF1 and RPS5 are suggested to be applicable as optimal reference genes in different developmental stages of the honeybee. However, when each C_q value of RAD1a, ARF1, and RPS5 in seven different developmental stages was compared, the highest C_q and lowest C_q of RPS5 and ARF1 were significantly different (p<0.05). In contrast, RAD1a showed statistically similar expression levels across developmental stages from egg to forager (p>0.05) (Fig. 2A), supporting that RAD1a is the most optimal reference gene for normalization of target gene using qRT-PCR in honeybee developmental stage.

In the present study, pairwise variations were also analyzed by geNorm to suggest an optimal number of reference genes for target gene normalization. According to the previous studies, V_n/V_n+1≤0.15 is the criterion (Vandesompele et al., 2002; Wan et al., 2010), but the value of pairwise variation for a combination of multiple reference genes were higher than 0.15 in this study (Fig. 3B). This indicates that the addition of the reference gene does not significantly contribute to the accurate normalization of the target gene in honeybee developmental stages, as suggested in the previous study (Silveira et al., 2009). A single gene, particularly RAD1a, showed the most negligible variation of C_q across different developmental stages (Fig. 2) and SD value in BestKeeper (Fig. 5). In addition, the M and stability value of RAD1a analyzed by geNorm (Fig. 3A) and NormFinder (Fig. 4) were under cutoff lines, suggesting that a single use of RAD1a as a reference gene can be appropriated in the gene expression study of the developmental stages of honeybee.

In conclusion, we evaluated the expression stabilities of five candidate reference genes in different developmental stages of the honeybee with C_q distribution analysis and three programs. Considering the cutoff line of CV≤0.5 for C_q distribution (Fig. 2B), stability value≤0.15 for NormFinder (Fig. 4), and SD≤1.0 for BestKeeper (Fig. 5), all five genes were applicable as the reference genes. However, in geNorm analysis, three genes (RPS5, RAD1a, and ARF1) were suggested as reference genes because their M values were ≤0.5 (Fig. 3A), but C_q values of RAD1a, rather than RPS5 and ARF1, were statistically similar across seven developmental stages (Fig. 2A), suggesting that RAD1a is suitable reference gene. Pairwise variation analysis further suggested a single use of the reference gene (Fig. 3B). Therefore, in this study, we finally suggest that RAD1a is the most appropriate reference gene for the accurate normalization of target gene expression at the developmental stage of the honeybee, but selected reference gene is remained to be validated by normalization of target gene.