Screening and Validation of Reference Genes for RT-qPCR Under Different Honey Bee Viral Infections and dsRNA Treatment

Honey bee viruses are one of the most important pathogens that have contributed to the decrease in honey bee colony health. To analyze the infection dynamics of honey bee viruses, quantification of viral gene expression by RT-qPCR is necessary. However, suitable reference genes have not been reported from viral and RNAi studies of honey bee. Here, we evaluated the expression of 11 common reference genes (ache2, rps18, β-actin, tbp, tif, rpl32, gadph, ubc, α-tubulin, rpl14, and rpsa) from Apis mellifera (Am) and Apis cerana (Ac) under Israeli acute paralysis virus (IAPV), chronic bee paralysis virus (CBPV), and Chinese sacbrood virus (CSBV) infection as well as dsRNA-PGRP-SA treatment, and we confirmed their validation by evaluating the levels of the defensin 1 and prophenoloxidase (ppo) genes during viral infection. Our results showed that the expression of selected genes varied under different viral infections. ache2, rps18, β-actin, tbp, and tif can be used to normalize expression levels in Apis mellifera under IAPV infection, while the combination of actin and tif is suitable for CBPV-infected experiments. The combination of rpl14, tif, rpsa, ubc, and ache2 as well as more reference genes is suitable for CSBV treatment in Apis cerana. Rpl14, tif, rps18, ubc, and α-tubulin were the most stable reference genes under dsRNA treatment in Apis mellifera. Furthermore, the geNorm and NormFinder algorithms showed that tif was the best suitable reference gene for these four treatments. This study screened and validated suitable reference genes for the quantification of viral levels in honey bee, as well as for RNAi experiments.


INTRODUCTION
Honey bees provide a significant pollination service for many agricultural crops and wild plants. However, a decline in bee populations has been recently observed in the United States and parts of European countries (Brutscher et al., 2015). Although multiple biotic and abiotic factors contribute to colony decline, several epidemiological and temporal monitoring studies have indicated that pathogens play an important role in colony loss (Brutscher et al., 2015). Among these pathogens, RNA viruses are the major impacting factors, including two Apis mellifera viruses, Israeli acute bee paralysis virus (IAPV) (Maori et al., 2007) and chronic bee paralysis virus (CBPV) (Genersch and Aubert, 2010), as well as one Apis cerana virus, namely, Chinese sacbrood virus (CSBV) .
Israeli acute paralysis virus, a positive-sense RNA virus in the family Dicistroviridae, has a widespread impact on honey bee health and has been linked with colony losses (Chen et al., 2014;Galbraith et al., 2015). RNA-seq data indicated that IAPV attacks several core genes in insect immune pathways (Toll, JAK-STAT and RNA interference pathways) (Galbraith et al., 2015). Although CBPV has not been classified, its genome, which is made up of two segments of single-stranded RNA in a non-enveloped anisometric capsid, shows homology to the Nodaviridae and Tombusviridae families (Coulon et al., 2017). Exposure of honey bees to CBPV alone promoted downregulation of immune-related genes, such as dorsal-1a and ppo (Coulon et al., 2017). CSBV, a geographic strain of SBV isolated from A. cerana in China, can cause fatal infections to A. cerana larvae by decreasing antimicrobial peptide levels (Shan et al., 2017). The quantification of viral and immune gene expression levels is one of the most important aspects when studying host and virus interactions (Marcial-Quino et al., 2016). RNA interference (RNAi) technology is frequently used as a research tool to study the function, transcription and regulation of gene including immune genes (Vyas et al., 2017;Zhu et al., 2019). Introducing exogenous double stranded RNA (dsRNA) into the insect cells activates RNAi pathways that normally function to induce antiviral responses (Dietrich et al., 2017). To analyze the expression differences in immune genes during viral infection and RNAi process, RT-qPCR is used to assess gene expression by normalization with reference genes. Most gene expression studies focused on honey bee select actin as an internal control (Chen et al., 2014), though RPS5 (Wheeler et al., 2006) and EIF S8 (Hunt et al., 2007) are also occasionally used. However, the identification of reference genes was not made in advance when these studies assessed viral dynamics and the changes in corresponding target genes.
As summarized by Shakeel et al. (2018), the identification and validation of reference genes for qRT-PCR have been performed in various insect species. For honey bee, the stability of reference genes has been investigated in different developmental stages and tissues (Lourenco et al., 2008;Reim et al., 2013;Moon et al., 2018a), and even after bacterial challenge (Scharlaken et al., 2008). However, suitable candidate reference genes have not been validated during honey bee viral infection. Moreover, most of these genes displayed variable expression levels under different experimental conditions (Thellin et al., 2009).

Collection of Virus-Infected Honey Bees
Israeli acute paralysis virus-infected A. mellifera adults with obvious paralysis symptom (Maori et al., 2007) were collected from three different apiaries in Beijing, Miyun, China, while CBPV-infected A. mellifera adults with paralysis symptom (Olivier et al., 2008) were collected from three different apiaries in Guangdong Province in China. They were tested by RT-PCR with special primers (Diao et al., 2018). Healthy A. cerana larvae were collected from three different apiaries in Guangdong Province in China. CSBV-infected larvae were taken from naturally infected colonies with obvious cystic phenotypes and symptoms. The samples were immediately transported to the laboratory on dry ice to avoid the degradation of the active substances therein.
We collected seemingly healthy bee samples from the experimental apiary at the Institute of Apicultural Research (IAR), Chinese Academy of Agricultural Sciences, Beijing, China. Newly emerged bees (A. mellifera) were obtained from brood frames taken from the experimental honey bee hives and maintained in an incubator at 30 • C and 60% relative humidity (RH) for approximately 12 h with access to a 50% sucrose solution. Honey bee samples were divided into five groups with three repetitions per group with 30 bees. The first, second, and third groups consisted of bees treated with IAPV, CBPV, and CSBV, respectively; the fourth group consisted of bees treated with dsRNA-pgrpsa; and the fifth group was used to test the reliability of the reference genes. Each group contained a blank control, except for the dsRNA experiment, which had dsRNA-GFP as a control.
To gain reliable samples, PCR detection was performed. Samples in IAPV group were collected every day (24, 48, 72, and 96 h) and followed by PCR using IAPV primers ( Supplementary  Table S1). Similarly, samples in CBPV group were collected every day (0, 24, 48, 72, and 96 h) and followed by PCR using CBPV primers (Supplementary Table S1). The naturally CSBV-infected larvae and the healthy larvae of A. cerana were collected from 2th to 6th-instar and followed by PCR using CSBV primers (Supplementary Table S1). Samples in dsRNA treatment group were collected every day (24, 48, 72, and 96 h) and followed by PCR using PGRP-SA primers (Supplementary Table S1).

Purification and Inoculation of IAPV and CBPV
To purify the viruses, the IAPV (IAPV-BJ, KX421583.1) and CBPV (CBPV-HB, MF175174.1, and MF175173.1) field strains originating from paralyzed bees showing clinical symptoms of viral infection were collected in May 2017 in China. IAPV and CBPV were screened for predominance in the samples, as well as the absence of other common viruses. IAPV was purified as described by Maori et al. (2007). CBPV was purified from the heads of bees by ultracentrifugation in a 10 to 40% (w/v) sucrose gradient as previously described by Olivier et al. (2008). To avoid interference by other viruses, PCR was run to identify whether there was infection by four other common honey bee viruses before virus purification (Supplementary Table S1), such as black queen cell virus (BQCV), acute bee paralysis virus (ABPV), deformed wing virus (DWV) and aphid lethal paralysis virus (ALPV) following the method of Diao et al. (2018) and Yang et al. (2019). The purified virus particles were dissolved in sucrose at −80 • C and aliquoted until use.
Purified virus particles (5 µl, 1 µg/µl) with about 4.0 × 10 4 genome equivalent copies was injected into the abdominal intersegment space of 30 individual newly emerged bees taken from a healthy colony. For the IAPV and CBPV infection groups, 5 µg of IAPV and CBPV, respectively, was injected into newly emerged honey bees (A. mellifera). The naturally infected CSBV and normal control A. cerana larvae were collected from the 2nd to 6th instars. For the dsRNA group, dsRNA-PGRP-SA and dsRNA-GFP (5 µg) were injected into newly emerged honey bees (A. mellifera) at 0 and 48 h, respectively. The honey bee samples were detected for six other viruses before they were used for the infection experiments. The infected bees were maintained in an incubator at 30 • C, and their mortality was calculated each day. Bees in control group were injected with PBS buffer.

RNA Extraction and RT-qPCR
Total RNA was extracted with the TRIol Kit (Ambion, Life Technologies, United States) following the manufacturer's instructions. A NanoDrop 2000 (Thermo Scientific, United States) was used to check the quality of each RNA sample. RNA samples with A260/280 ratios ranging from 1.8 to 2.2 were used for cDNA synthesis with the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) according to the manufacturer's instructions. The cDNAs of all samples were stored at −20 • C until use.
The specific primers for the reference genes were designed using Primer-BLAST ( Table 1). The RT-qPCR reaction system consisted of 12.5 µL 2 × SYBR Premix Ex Taq TM II (Takara, Dalian, China), 0.5 µL each of the forward and reverse primers (10 mM), 2 µL cDNA template, and 9.5 µL double-distilled H 2 O for a total volume of 25 µL. The RT-qPCR conditions followed the manufacturer's instructions. All reactions were performed using SYBR Green Premix and amplified under the following cycling conditions: an initial cycle at 95 • C; 40 cycles of 95 • C for denaturation, 25 s at 55 • C for annealing and 20 s at 72 • C for extension; and the generation of a melting curve consisting of a single peak to rule out non-specific products and primer dimers afterward. The RT-qPCR analysis was performed with three biological replicates for each sample and three technical replicates for each biological replicate and was measured as the mean Ct value. The results were analyzed using the 9600 plus Software.

PCR Detection for Viral and Immune Genes
The primer sequences, orientation and references are provided in Supplementary Table S1. The cDNAs of all samples were gained in the above, followed by the PCR cycles. The PCR reaction consisted of a total 20 µL volume containing 10 µL 2 × GoTaq reaction buffer (Promega, United States), 0.5 µL 10 µM of the sense and antisense primers, 1 µL of cDNA, and 8 µL nucleasefree water. The cycling conditions were as follows: 1 min at 95 • C; 33 cycles of 30 s at 94 • C, 30 s at 55 and 72 • C for 1 min; a final extension of 10 min at 72 • C; and cooling to 4 • C. The PCR amplification products were separated in a 2% agarose gel stained with GV II (BIOMEC, China) and photographed with a FR-200A luminescent and fluorescent biological image analysis system (Furi, China). The product size was determined using a 100-bp molecular size ladder.

dsRNA Synthesis
The specific primers for PGRP-SA and GFP tagged with T7 promoter sequence were designed using Primer-BLAST to amplify the templates for synthesis of dsRNA (Supplementary Table S2). Amplification was performed under normal PCR conditions as follows: 30 s at 95 • C; followed by 35 cycles of 8 s at 95 • C, 30 s at 56 • C, and 30 s at 72 • C; and a final extension incubation at 72 • C for 5 min. The PCR templates were purified using a Qiagen purification kit (Qiagen, Germantown, MD, United States). dsRNA-PGRP-SA and dsRNA-GFP were produced using a T7 RNAi Transcription Kit (Vazyme, China) according to the manufacturer's instructions. The reaction product was subjected to DNase and RNAse digestion, and then incubated at 37 • C for 8 h. Subsequently, the samples were dried at 37 • C for 10 min and resuspended in 20-40 µL of nuclease-free H 2 O. The quality of the dsRNA was checked by electrophoresis and quantified with a spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). Then, 5 µg each of PGRPSA dsRNA and GFP dsRNA was injected into newly emerged A. mellifera at 0 and 48 h, respectively.

Stability of Expression of Selected Reference Genes
The stability of expression of the candidate reference genes was evaluated by five statistical algorithms, namely, the delta CT method (Silver et al., 2006), BestKeeper (Pfaffl et al., 2004), geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and the online platform RefFinder (Xie et al., 2012). geNorm is based on the principle that the expression ratio of two reference genes should be constant in all samples, regardless of the experimental conditions or sampling time (Hellemans et al., 2007). The gene expression stability (M) is defined as the average pairwise variation compared with all other tested candidate genes by geNorm, and the reference genes with the lowest M-value remained (Vandesompele et al., 2002). Pairwise variation (Vn/Vn + 1) analysis between the normalization factors (NFn and NFn + 1) was performed by geNorm to determine the optimal number of reference genes (Vandesompele et al., 2002). NFn was also calculated by the stepwise inclusion of a less stable gene until the (n + 1)th gene had no significant contribution to the newly calculated NFn + 1 (Vandesompele et al., 2002). In particular, if the pairwise variation Vn/n + 1 between the two sequential normalization factors NFn and NFn + 1 was lower than the cut-off value of 0.15, it was considered that NFn + 1 is not required (Peng et al., 2018). NormFinder was used to estimate changes and intra-and intergroup expression variation among the candidate reference genes (Andersen et al., 2004). The most stable gene is the one with the lowest stability value according to the intra-and intergroup variability of each gene (Andersen et al., 2004). BestKeeper determines the most stably expressed genes according to the standard deviation (SD) and coefficient of variation (CV) of all the Ct values for each gene (Pfaffl et al., 2004). When the SD values of genes are <1, they are considered stable (Marcial-Quino et al., 2016). Similarly, the delta-Ct method was used to identify suitable housekeeping genes based on the average SD value of each gene (Silver et al., 2006). To identify the optimal reference genes, RefFinder, a user-friendly web-based comprehensive tool that provides a comprehensive ranking evaluation based on the geometric mean and integrates with other computational algorithms (geNorm, NormFinder, BestKeeper, and the comparative Delta-Ct method), was used to reduce bias or avoid contradictory results caused by using individuals (Peng et al., 2018). Individual genes were assigned an appropriate weight and the geometric mean of their weights was calculated for the overall final ranking based on each program.

Validation of the Selected Reference Genes
The three strategies for expression profile normalization using the reference genes in each treatment are as follows: (1) the optimal reference gene from all the samples, (2) the optimal reference gene from each treatment, and (3) the least stable reference gene from each treatment. To validate the candidate reference genes, the relative expression levels of two immune genes (defensin1 and prophenoloxidase) were used to evaluate the reference genes by RT-qPCR. The 2 − CT method was employed to validate and quantify the expression levels of def1 and ppo after viral infection. One-way ANOVA was used to analyze the significance and the means were compared using Tukey's multiple comparison test (P < 0.05) using the GraphPad Prism software.

Identification of Primer Specificity and the Efficiency of the Candidate Reference Genes
The reference gene candidates in this study were selected based on the following two criteria: (a) genes routinely used in honey bees or insects for transcript normalization and (b) genes described as stably expressed under environmental stress. Eleven candidate reference genes were selected for A. mellifera and A. cerana under different treatments. The amplicon sizes of the target genes were between 100 and 320 bp with primer lengths of 17-25 bp and a 40-60% GC content. Primer specificity for each gene was validated by melting curve analysis during qRT-PCR, and a single peak was observed in each reaction (Figures 1A,B). Primer specificity for all genes was further tested using 2.0% agarose gels. A single PCR amplification band of the expected size was observed, suggesting that primer dimers and nonspecific amplified products were not generated (Figures 1C,D). Additionally, the amplification efficiencies ranged from 92.2 to 107% based on the standard curves for the 18 genes. Moreover, correlation coefficients (R 2 ) were greater than 0.99, which fit the RT-qPCR requirements (Supplementary Table S3).

Characterization of Candidate Reference Genes
To analyze the expression profiles of the candidate reference genes, sample validation was performed by PCR (Supplementary Table S1). As shown in Supplementary Figure S1, the samples from the IAPV-treated group were positive but the control groups were not (Supplementary Figure S1A). Similarly, positive results were also observed for the CBPV and CSBV groups (Supplementary Figures S1B,C). The results for the dsRNA treatment groups showed that the expression of PGRP-SA had been markedly downregulated compared to the control group (Supplementary Figure S1D).
A reliable reference gene should show constant expressions level under different experimental conditions (Galli et al., 2015). We next evaluated the expression levels of the candidate reference genes and found that the average Ct values of the reference genes in different samples ranged from 16 to 32 (Figure 2). In the A. mellifera groups, RPS18 was the most abundant with the lowest mean Ct value at 19.43 ± 0.67 in the IAPV-treated group (Figure 2A), 17.31 ± 0.37 in the CBPV-treated group ( Figure 2B) and 18.52 ± 0.96 in the dsRNA-treated group ( Figure 2C). In contrast, Amα-tubulin transcription showed the lowest level with the highest mean Ct value at 30.99 ± 1.33 in the IAPV-treated group, 29.58 ± 0.62 in the CBPV-treated and 31.33 ± 0.80 in the dsRNA-treated group (Figure 2). In the A. cerana group, ache2 exhibited the smallest SD value (28.6 ± 1.47) but had the lowest expression level (Figure 2D). In contrast, rpl14 exhibited the smallest Ct value (19.59 ± 1.56) but had a larger SD value. Tif was a moderately stable gene with Ct average values at 25.25 ± 0.73 in the IAPV-treated group, 24.38 ± 0.26 in the CBPV-treated group, 26.10 ± 0.9 in the dsRNA-treated group and 24.43 ± 1.57 in the CSBV-infected A. cerana groups (Figure 2). However, reference genes from the CSBV-infected group exhibited larger SD values, suggesting that studies should carefully consider when only using one gene as an internal reference for the quantification of genes related to CSBV in A. cerana ( Figure 2D).

Stability of the Candidate Reference Genes Under Different Viral Infections and dsRNA Treatment
The rankings for the expression stability of these 11 genes in honey bee under treatment with IAPV, CBPV, CSBV or dsRNA was analyzed using different algorithms (Table 1 and Figure 3). In the IAPV treatment group, Amrps18, Amtif, and Amache2 were considered the most stable genes with the lowest M value according to the geNorm and NormFinder analysis ( Table 2). However, the ranking obtained with BestKeeper and Delta CT showed that Amache2, Amtbp, and Amβ-actin were the top three stable candidate reference genes ( Table 2). To identify the most suitable reference genes, the results of these four programs were integrated with RefFinder (He et al., 2019), which showed that the stability of these genes in A. mellifera treated with IAPV was as follows: Amache2, Amrps18, Amactin, Amtbp, Amtif, Amrpl32, Amgadph, Amubc, and Amα-tubulin ( Figure 3A). However, Amubc and Amα-tubulin were ranked as the least stable genes in the IAPV-infected group. As shown in Table 2, Amtif and Amtbp were the top two stable candidate reference genes in A. mellifera infected with CBPV based on the ranking with geNorm, Delta-CT and BestKeeper. However, the ranking obtained with NormFinder was Amtif, Amactin and Amrps18, which was different from that obtained with geNorm. Based on the comprehensive analysis by RefFinder, the stability of these genes in A. mellifera treated with CBPV was as follows: Amβ-actin, Amtif, Amtbp, Amrps18, Amgadph, Amache2, Amrpl32, Amubc, and Amα-tubulin ( Figure 3B). Because Amubc and Amα-tubulin were ranked as the least stable genes in IAPV-infected and CBPV-infected A. mellifera, they were not recommended as reference genes for qRT-PCR analyses in honey bee infected with viruses (Figures 3A,B).
For A. cerana larvae infected with CSBV, Actif, Acache2, and Acrpl14 were the top three stable candidate reference genes based on the geNorm algorithms (Table 2). Acrpl14, Actif, and Acrpsa were the top three stable candidate reference genes according to NormFinder. However, the ranking obtained with BestKeeper showed that Acache2, Acα-tubulin and Actbp were the top three stable candidate reference genes. Furthermore, the ranking obtained with the Delta CT method showed that Acache2, Actbp, and Acrpsa were the top three suitable genes ( Table 2). The outcome from integration with RefFinder showed that the stability of the genes in A. cerana infected with CSBV was as follows: Acrpl14, Actif, Acrpsa, Acubc, Acache2, Acα-tubulin, Acacetin, Actbp, and Acgadph ( Figure 3C). However, according to their larger M values (1.3 to 2.0), none of the reference genes used in CSBV-infected A. cerana was stable ( Table 2).
For A. mellifera treated with dsRNA-PGRP-SA and dsRNA-GFP, Amtif, Amrpl32 and Amubc were the top three stable candidate reference genes according to the geNorm algorithms ( Table 2). Amtif, Amrpl32 and Amactin were the top three stable FIGURE 4 | Determination of the optimal number of reference genes for normalization based on the pairwise variation values calculated by geNorm. The value of Vn/Vn + 1 < 0.15 means n should be the optimal number of reference genes selected for RT-qPCR analysis among IAPV-infected (black), CBPV-infected (green), and CSBV-infected (blue) honey bees, as well as honey bees exposed to dsRNA-PGRPSA (purple).
candidate reference genes based on NormFinder. However, the ranking obtained with BestKeeper showed that Amα-tubulin, Amache2, and Amrpl32 were the top three stable candidate reference genes. Moreover, the ranking by the Delta-CT method showed that Amα-tubulin, Amache2 and Amtif were the top three reference genes ( Table 2). According to the RefFinder results, the stability of the genes in A. mellifera treated with dsRNA was as follows: Amrpl32, Amtif, Amrps18, Amubc, Amα-tubulin, Amactin, Amache2, Amtbp, and Amgadph ( Figure 3D).

Determination of the Optimal Number of Reference Genes
To determine the optimal number of reference genes, we used geNorm to calculate the pairwise variation. In the IAPV treatment group, V 3 /V 4 < 0.15 indicated that three reference genes was sufficient to accurately normalize RT-qPCR data (Figure 4). Therefore, Amache2, Amrps18 and Amactin were considered to be a suitable combination of reference genes for RT-qPCR analyses in honey bee infected with IAPV. In the CBPV treatment group, V 2 /V 3 < 0.15 indicated that two reference genes was sufficient to normalize RT-qPCR data. Thus, Amactin and Amtif were considered to be a suitable combination of reference genes for RT-qPCR. In contrast, in the CSBV treatment group, Vn/Vn + 1 > 0.15 indicated that more reference genes are required and should be combined for RT-qPCR data normalization. In the dsRNA treatment group, V 4 /V 5 < 0.15 suggested that rpl32, tif, rps18, and ubc were a suitable combination of reference genes for RT-qPCR analyses in honey bee treated with dsRNA (Figure 4).

Validation of the Selected Reference Genes
As shown in Figure 5A, the relative expression levels of the defensin 1 and ppo genes varied based on the different reference genes. In the IAPV-infected group, the relative expression levels of defensin 1 were downregulated from 5.1-to 2.9-fold and from 3.2-to 2.6-fold when Amache2 and Amtif were used as reference genes, respectively, while the expression levels of defensin 1 were upregulated from 2.5-to 7-fold if Amα-tubulin was the reference gene. In the CBPV-infected group, the normalization of ppo expression levels to Amactin or Amtif as reference genes exhibited similar downregulated expression trends for both CBPV-infected and uninfected honey bees ( Figure 5B). However, the relative expression levels of ppo were upregulated to 2.25fold at 96 h after CBPV infection when Amα-tubulin was the reference gene ( Figure 5B). In the CSBV-infected group, the relative expression levels of defensin 1 were upregulated from 1.4to 5.6-fold and from 1.2-to 3.4-fold in 5th and 6th CSBV-infected honey bee larvae when we used Acrpl14 and Actif as reference genes, respectively, while the relative expression levels of defensin 1 were upregulated from 7.3-to 15.3-fold when Acgadph used as the reference gene ( Figure 5C). In contrast, when Acgadph was used for normalization, the expression levels of defensin 1 exhibited a significant difference from that of Acrpl14 (P < 0.01). In the dsRNA treatment group, the relative expression levels of pgrpsa were downregulated 0. 16-, 0. 06-, and 0.28-fold after 72 h of dsRNA-pgrpsa treatment when Amrpl32 (best), Amgadph (worst) and Amtif (modest) were used as the reference genes ( Figure 5D), respectively, while the relative expression levels of pgrpsa were upregulated 6. 5-, 2. 4-, and 7.2-fold at 72 h after dsRNA-gfp treatment if we used Amrpl32, Amgadph and Amtif as the reference genes, respectively ( Figure 5D). This result showed that normalization of pgrpsa expression using the best (Amrpl32) and modest (Amtif ) reference genes exhibited nearly similar outcomes after dsRNA treatment, though its expression exhibited significant differences (P < 0.05 and P < 0.01) compared to those of Amrpl32 and Amgadph as reference genes ( Figure 5D).

DISCUSSION
RT-qPCR is an important, simple and practical technique for assessing gene expression when compared to other quantitative methods such as Northern blotting, in situ hybridization and RNA-seq technology (Peng et al., 2018;Adeyinka et al., 2019;He et al., 2019). The application of accurate reference genes is crucial for the quantification of gene expression in A. mellifera and A. cerana. In our study, five algorithms were used to identify the stability of reference genes. Ache2, rps18, actin, tbp, and tif are suitable to normalize the gene expression levels in IAPVinfected A. mellifera. The combination of actin and tif is suitable for CBPV-infected A. mellifera, while the combination of rpl14, tif, rpsa, ubc, and ache2 as well as more reference genes is suitable for CSBV-infected A. cerana larvae. This result suggests that further studies should be performed to select candidate reference genes for CSBV infection. Moreover, rpl14, tif, rps18, ubc, and αtubulin showed the most stable expression in A. mellifera under dsRNA treatment.
Our results showed that reference genes were different under different viral infection conditions. An increasing number of studies have reported that using a single reference gene to normalize RT-qPCR might be insufficient and inaccurate for the quantification of gene expression (Yang et al., 2018). As shown in our four experiments, the number of suitable reference genes was greater than two (Figure 4). Additionally,  the stability of these reference genes was different under different experimental conditions. As determined by RefFinder (Figure 3), tif was a relatively stable reference gene under all of our experimental conditions. Actually, tif have been widely used as housekeeping genes for gene expression analysis during Bombyx mori Cytoplasmic Polyhedrosis Virus (BmCPV) and Bombyx mori Bidensovirus (BmBDV) infection (Guo et al., 2016). Ache2 was the most stable reference gene under IAPV infection, while β-actin was the most stable reference gene under CBPV infection. In contrast, α-tubulin exhibited the worst stability under IAPV and CBPV infection. AcRpl14 and Amrpl32 were the most stable reference genes under CSBV infection and dsRNA treatment, respectively, while Acgadph and Amgadph exhibited the worst stability under CSBV infection and dsRNA treatment, respectively. However, gapdh and rps18 were suggested to be the optimal reference genes for RT-qPCRbased under the labor-specific gene expression, bacterial infected and developmental gene expression in Western honey bee as shown in Table 2 (Scharlaken et al., 2008;Reim et al., 2013;Moon et al., 2018a). Previous studies have shown that some ribosome-associated genes have been used as stable internal reference genes for quantitative analysis (Wang et al., 2019). For instance, Freitas et al. (2019) found that rpl32, rps5, and rps18 were identified as suitable reference genes for the different tissues of stingless bees. Similarly, rpl23 was identified the reliable reference gene in bumblebees challenged by IAPV (Niu et al., 2014; Table 3). Although our study showed that rpl32 displayed stable expression in A. mellifera under dsRNA treatment, it exhibited the worst stability during IAPV and CBPV infection in A. mellifera (Figure 3), while rps18 was the better reference gene in A. mellifera under dsRNA treatment or viral infection (Table 3). Likewise, the stability of the same reference gene varied under different experimental conditions. Although actin is used as an internal control in most gene expression studies during honey bee viral infection, such as IAPV infection (Chen et al., 2014), CSBV infection (Shan et al., 2017), and CBPV infection (Coulon et al., 2017), actin exhibited poorer stability during CSBV infection (Figure 3). Gadph, the other most commonly used reference gene, also showed expression variability from the 2nd to 6th honey bee infected with CSBV infection as well as dsRNA treatment (Figure 3). Similarly, α-tubulin, another commonly used reference gene, exhibited worse stability under IAPV and CBPV infection. This might result due to interactions between α-tubulin and acetylation, as Zhang et al. (2016) found that acetylated α-tubulin will enhance viral titers by promoting fusion with viral inclusion bodies.
This study systematically analyzed 11 reference genes for RT-qPCR under viral infection in A. mellifera and A. cerana. These results will facilitate future studies in honey bee under different viral infection conditions.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
CH and YD conceived the manuscript. YD, HZ, SY, LiZ, and LinZ performed the experiments. YD analyzed the data and wrote the manuscript. CH revised the manuscript. All authors contributed to the article and approved the submitted version.