Evaluation of Housekeeping Genes for Quantitative Real-Time PCR Analysis of Bradysia odoriphaga (Diptera: Sciaridae)

The soil insect Bradysia odoriphaga (Diptera: Sciaridae) causes substantial damage to Chinese chive. Suitable reference genes in B. odoriphaga (Bradysia odoriphaga) have yet to be identified for normalizing target gene expression among samples by quantitative real-time PCR (qRT-PCR). This study was focused on identifying the expression stability of 12 candidate housekeeping genes in B. odoriphaga under various experiment conditions. The final stability ranking of 12 housekeeping genes was obtained with RefFinder, and the most suitable number of reference genes was analyzed by GeNorm. The results revealed that the most appropriate sets of internal controls were RPS15, RPL18, and RPS18 across developmental phases; RPS15, RPL28, and GAPDH across temperatures; RPS15 and RPL18 across pesticide treatments; RSP5, RPS18, and SDHA across photoperiods; ACTb, RPS18, and RPS15 across diets; RPS13 and RPL28 across populations; and RPS15, ACTb, and RPS18 across all samples. The use of the most suitable reference genes versus an arbitrarily selected reference gene resulted in significant differences in the analysis of a target gene expression. HSP23 in B. odoriphaga was found to be up-regulated under low temperatures. These results will contribute to the standardization of qRT-PCR and will also be valuable for further research on gene function in B. odoriphaga.


Introduction
Quantitative real-time PCR (qRT-PCR) is considered as a reliable technique for the gene quantification [1][2][3]. However, gene expression can be affected by many confounding factors, such as RNA extraction, reverse transcription, and qRT-PCR efficiency [4,5]. Therefore, housekeeping genes are commonly used as "reference genes" to decrease the effects due to confounding factors and to increase the accuracy of the quantification analysis related to the particular biological environment [6,7]. The reference genes overcome the whole steps of the analyses along with interest genes and suppress the variations within the treatment group to the lowest level. Determining the number and identity of the reference genes to be employed for count data of normalization factors (NF) among comparable

Expression Images of Candidate Reference Genes
To analyze mRNA expression level of the 12 candidate housekeeping genes, C t values were calculated for all samples in this study. As shown in Figure 2, the mean C t values of the 12 candidates were <30. The average C t value was lowest for RPL28 (15.95) and highest for TUB (25.32).

Expression Images of Candidate Reference Genes
To analyze mRNA expression level of the 12 candidate housekeeping genes, Ct values were calculated for all samples in this study. As shown in Figure 2, the mean Ct values of the 12 candidates were <30. The average Ct value was lowest for RPL28 (15.95) and highest for TUB (25.32).

Stability of Reference Genes
The following results are based on analyses across the range of each factor. For developmental stage, for example, stability is based on an analysis across all stages.

Developmental Stages
According to the four algorithms, TUB and EF1a were the least steady across developmental stage ( Table 2). The most stable genes (in order) were RPS15, RPL18, and ACTb according to the ΔCt method; RPS18, RPS13, and RPL28 according to BestKeeper; SDHA, ACTb, and GAPDH according to NormFinder; and RPL18, RPS15, and RPS18 according to GeNorm (Table 2).
According to RefFinder, the order of the reference gene stability across developmental stages was: Figure 3A). GeNorm analysis results showed that the pair-wise values of V2/3 to V6/7 were all above the cut-off value of 0.15 but that the pair-wise value of V7/8 was <0.15 ( Figure 4); a value <0.15 indicates that the supplemental reference genes will not evidently change the normalization. Based on the RefFinder recommendations for selection of reference genes and on convenience of operation, RPS15, RPL18, and RPS18 were considered suitable reference genes across developmental stages of B. odoriphaga (Table 3).

Stability of Reference Genes
The following results are based on analyses across the range of each factor. For developmental stage, for example, stability is based on an analysis across all stages.

Developmental Stages
According to the four algorithms, TUB and EF1a were the least steady across developmental stage ( Table 2). The most stable genes (in order) were RPS15, RPL18, and ACTb according to the ∆C t method; RPS18, RPS13, and RPL28 according to BestKeeper; SDHA, ACTb, and GAPDH according to NormFinder; and RPL18, RPS15, and RPS18 according to GeNorm (Table 2).
According to RefFinder, the order of the reference gene stability across developmental stages was: Figure 3A). GeNorm analysis results showed that the pair-wise values of V2/3 to V6/7 were all above the cut-off value of 0.15 but that the pair-wise value of V7/8 was <0.15 ( Figure 4); a value <0.15 indicates that the supplemental reference genes will not evidently change the normalization. Based on the RefFinder recommendations for selection of reference genes and on convenience of operation, RPS15, RPL18, and RPS18 were considered suitable reference genes across developmental stages of B. odoriphaga (Table 3).    Pair-wise variations (V) Figure 4. Pair-wise variation (Vn/Vn + 1) analysis of the number of candidate reference genes in B. odoriphaga. Pair-wise variation was analyzed by GeNorm software. A value <0.15 indicates that the normalization could not be dramatically changed by additional reference genes.

Temperatures
According to the ΔCt method and NormFinder, the most steady candidate genes across temperature treatments were RPS15, RPL28, and GAPDH, and the least stable were RPS18, SDHA, and TUB (Table 2). According to BestKeeper, the most stable candidate genes were RPL28, RPS15, and UBCE, and the least steady were RPS18, EF1a, and ACTb (Table 2). According to GeNorm, the most stable candidates were RPL18, RSP5, and RPL28, and the least stable were TUB, SDHA, and ACTb (Table 2).
According to RefFinder, the order of reference gene stability across temperatures was: Figure 3B). The GeNorm data predicted that the pair-wise values from V2/3 to V3/4 were <0.15 ( Figure 4). Therefore, RPS15, RPL28, and GAPDH were considered stable candidate genes across the tested temperatures (Table 3).

Pesticides
TUB, GAPDH, and EF1a were regarded as the least steady genes across pesticide treatments by the ΔCt method and by GeNorm and NormFinder but not by BestKeeper (Table 2). According to the comparative ΔCt method and GeNorm, the most stable candidates were RPS15, RPL18, and RPL28 (Table 2), while they were RPS15, RPL18, and RPS18 by using NormFinder and were SDHA, EF1a, and ACTb according to BestKeeper (Table 2).
Based on RefFinder, the order of reference gene stability across pesticide treatments was: Figure 3C). The GeNorm analysis showed that the pair-wise value of V2/3 was <0.15 ( Figure 4). Therefore, RPS15 and RPL18 were considered suitable candidate genes across the tested pesticide treatments (Table 3).
According to RefFinder, the order of reference gene stability across photoperiod treatments was: Figure 3D). The GeNorm analysis data showed that only the pair-wise value of V7/8 was

Temperatures
According to the ∆C t method and NormFinder, the most steady candidate genes across temperature treatments were RPS15, RPL28, and GAPDH, and the least stable were RPS18, SDHA, and TUB (Table 2). According to BestKeeper, the most stable candidate genes were RPL28, RPS15, and UBCE, and the least steady were RPS18, EF1a, and ACTb (Table 2). According to GeNorm, the most stable candidates were RPL18, RSP5, and RPL28, and the least stable were TUB, SDHA, and ACTb (Table 2).
According to RefFinder, the order of reference gene stability across temperatures was: Figure 3B). The GeNorm data predicted that the pair-wise values from V2/3 to V3/4 were <0.15 ( Figure 4). Therefore, RPS15, RPL28, and GAPDH were considered stable candidate genes across the tested temperatures (Table 3).

Pesticides
TUB, GAPDH, and EF1a were regarded as the least steady genes across pesticide treatments by the ∆C t method and by GeNorm and NormFinder but not by BestKeeper (Table 2). According to the comparative ∆C t method and GeNorm, the most stable candidates were RPS15, RPL18, and RPL28 (Table 2), while they were RPS15, RPL18, and RPS18 by using NormFinder and were SDHA, EF1a, and ACTb according to BestKeeper (Table 2).
Based on RefFinder, the order of reference gene stability across pesticide treatments was: Figure 3C). The GeNorm analysis showed that the pair-wise value of V2/3 was <0.15 ( Figure 4). Therefore, RPS15 and RPL18 were considered suitable candidate genes across the tested pesticide treatments (Table 3).
According to RefFinder, the order of reference gene stability across photoperiod treatments was: Figure 3D). The GeNorm analysis data showed that only the pair-wise value of V7/8 was below the cut-off value of 0.15 ( Figure 4). RSP5, RPS18, and SDHA were considered to be the most stable candidate genes across photoperiod treatments (Table 3).

Diets
Both NormFinder and ∆C t method results shared the same stable genes (ACTb, RPS18, and RPS15) across diets and confirmed SDHA, EF1a, and RSP5 as the least steady genes across diets ( Table 2). According to BestKeeper, the most steady genes were RPS15, EF1a, and GAPDH, and the least stable were RPL18, UBCE, and SDHA (Table 2). According to GeNorm, the most stable genes were RPL18, RPS18, and ACTb, and the least stable were SDHA, RSP5, and EF1a (Table 2).
According to RefFinder, the ranking order of reference gene stability across diets was: Figure 3E). The GeNorm analysis showed that the pair-wise value of V4/5 was <0.15 ( Figure 4). Therefore, ACTb, RPS18, and RPS15 were considered fitted reference genes across diets (Table 3).

Populations
Across B. odoriphaga populations, TUB, RPL18, and UBCE were identified as the least stable genes by all the four algorithms ( Table 2). The most stable genes were RPS13, RPS15, and GAPDH according to the comparative ∆C t method; RPL28, SDHA, and GAPDH according to BestKeeper; RPS13, RPS15 and RPL28 according to NormFinder; and EF1a, RSP5, and GAPDH according to GeNorm (Table 2).
According to RefFinder, the order of reference gene stability across populations was: Figure 3F). The GeNorm analysis showed that V2/3 value was <0.15 ( Figure 4). Therefore, RPS13 and RPL28 were considered suitable reference genes for gene expression (Table 3).

Ranking of Reference Genes for All Specimens
Across all samples, the three computational programs, and the comparative ∆C t method ranked RSP5, RPS13, and TUB as the least stable genes ( Table 2). The most stable genes were RPS15, ACTb, and RPL18 according to the ∆C t method; RPS18, ACTb, and RPL28 according to BestKeeper; ACTb, RPS15, and RPS18 according to NormFinder; and RPL28, RPS15, and RPL18 according to GeNorm (Table 2). Based on RefFinder, the order of reference gene stability across all samples was: Figure 3G). The GeNorm analysis showed that only the pair-wise values of V6/7 to V7/8 were less than the cut-off value of 0.15 ( Figure 4). Therefore, RPS15, ACTb, and RPS18 were regarded as the most suitable reference genes for qRT-PCR (Table 3).

Target Gene Expression
The selection failure of internal controls led to remarkable differences in quantification target genes. The relative expression level of HSP23 significantly differed among temperature treatments (4,´5, or´10˝C) when normalized by the most stable reference genes (such as RPS15) ( Figure 5). Similar changes observed in analyzing relative expression level of HSP23 with the normalization of two reference genes (such as RPS15 and RPL28) ( Figure 5) or three reference genes (such as RPS15, RPL28, and GAPDH) ( Figure 5). HSP23 in B. odoriphaga was found to be up-regulated under low temperatures, especially when the temperature was below´10˝C. However, HSP23 expression did not significantly differ among these treatments when expression was calculated with an arbitrary reference gene (such as ACTb) ( Figure 5). *** ** *** * Figure 5. Relative expression of a target gene, HSP23, was affected by three temperature treatments and standardized with different numbers, and kinds of reference genes. The expression level was separately normalized by: A (RPS15); B (RPS15 and RPL28); C (RPS15, RPL28 and GAPDH); or D (ACTb) reference genes. The reference genes were selected depending on the expression stability of the 12 housekeeping genes among the three temperature treatments. Values are means ± SD of three biology replications; the "*" means remarkable differences, * p < 0.05; ** p < 0.01; *** p < 0.001.

Discussion
Results obtained with qRT-PCR depend on several critical factors including RNA quantity, primer efficiency, and an internal control, i.e., a reference gene. When mRNA expression level is determined by qRT-PCR, the RNA must be intact, and primer efficiency must be determined [26]. Here, the OD ratio (A260/A280) of all RNA samples were between 1.8 and 2.0, and the amplification efficiency of the 12 candidates ranged from 90% to 110% (all R 2 > 0.990) ( Table 1). Thus, the quality of the RNA and amplification was sufficient for qRT-PCR.
Previous researches have reported that expression level of reference genes is not always stable under all experimental conditions [27][28][29] and that mRNA expression levels varied among several housekeeping genes [2,30]. These earlier findings were confirmed in the current study with B. odoriphaga (Table 2). In the current study, none of the candidate genes exhibited the same level of expression under all experiment conditions [31]. This indicates that reference genes need to be optimized and chosen depending on experimental parameters. Our data showed that, among the tested genes, mRNA expression of RPS15 was the most stable across development stages, temperatures, pesticide treatment, and all samples of B. odoriphaga, which is consistent with previous studies concerning development stage and temperature treatments for Nilaparvata lugens [9] and insecticide treatments for Helicoverpa armigera [32]. In B. odoriphaga, RSP5 was the most stable gene across photoperiod treatments, while RPS13 was the most stable across populations.
Previous studies have reported high expression stability for genes in the ribosomal protein (RP) genes family [27,33]. For example, among different organs, geographic populations, pesticide treatments, and starvation treatments, expression stability in Nilaparvata lugens was highest for RPS11 [9]; among different organs and developmental stages of Tetranychus cinnabarinus, Figure 5. Relative expression of a target gene, HSP23, was affected by three temperature treatments and standardized with different numbers, and kinds of reference genes. The expression level was separately normalized by: A (RPS15); B (RPS15 and RPL28); C (RPS15, RPL28 and GAPDH); or D (ACTb) reference genes. The reference genes were selected depending on the expression stability of the 12 housekeeping genes among the three temperature treatments. Values are means˘SD of three biology replications; the "*" means remarkable differences, * p < 0.05; ** p < 0.01; *** p < 0.001.

Discussion
Results obtained with qRT-PCR depend on several critical factors including RNA quantity, primer efficiency, and an internal control, i.e., a reference gene. When mRNA expression level is determined by qRT-PCR, the RNA must be intact, and primer efficiency must be determined [26]. Here, the OD ratio (A 260 /A 280 ) of all RNA samples were between 1.8 and 2.0, and the amplification efficiency of the 12 candidates ranged from 90% to 110% (all R 2 > 0.990) ( Table 1). Thus, the quality of the RNA and amplification was sufficient for qRT-PCR.
Previous researches have reported that expression level of reference genes is not always stable under all experimental conditions [27][28][29] and that mRNA expression levels varied among several housekeeping genes [2,30]. These earlier findings were confirmed in the current study with B. odoriphaga (Table 2). In the current study, none of the candidate genes exhibited the same level of expression under all experiment conditions [31]. This indicates that reference genes need to be optimized and chosen depending on experimental parameters. Our data showed that, among the tested genes, mRNA expression of RPS15 was the most stable across development stages, temperatures, pesticide treatment, and all samples of B. odoriphaga, which is consistent with previous studies concerning development stage and temperature treatments for Nilaparvata lugens [9] and insecticide treatments for Helicoverpa armigera [32]. In B. odoriphaga, RSP5 was the most stable gene across photoperiod treatments, while RPS13 was the most stable across populations.
Previous studies have reported high expression stability for genes in the ribosomal protein (RP) genes family [27,33]. For example, among different organs, geographic populations, pesticide treatments, and starvation treatments, expression stability in Nilaparvata lugens was highest for RPS11 [9]; among different organs and developmental stages of Tetranychus cinnabarinus, expression stability was highest for RPS18 [34]; in Phenacoccus solenopsis, expression stability among temperature treatments was highest for RPL32 [35]; among different developmental stages of Schistocerca gregaria, expression stability was highest for RPL49 [36]; among different organs and developmental stages of Cimex lectularius, expression stability was highest for RPL18 [37]; in Spodoptera litura, expression stability among different larval tissues, populations, and food treatments was highest for RPL10 [33]; in Plutella xylostella, expression stability among different developmental stages and photoperiods was highest for RPS13 [38]; in response to virus infection in Tribolium castaneum, expression stability was highest for RPS3 [39]; and in Helicoverpa armigera, expression stability among temperature treatments was highest for RPL28 [40]. As a principal component of ribosomes, ribosomal protein (RP) is important in intracellular protein biosynthesis, DNA repair, cell differentiation, etc. [31]. These results indicate that ribosomal protein genes might be useful as reference genes in interest gene expression studies.
In the current study with B. odoriphaga, however, an exception was that RPS13 showed the least steady expression across photoperiod treatments. Another exception was reported for Rhodnius prolixus: RPL26 was the most variable gene in the salivary glands of starved and non-starved specimens [41].
Because actin is the main structural protein of the cellular skeleton and is important for cell function [42], expression of the actin gene is substantial in most cell types [43]. The actin gene is the most stable gene in Chilo suppressalis [44], Schistocerca gregaria [36], and Apis mellifera [45]. Our study showed that ACTb is an ideal reference gene in B. odoriphaga subjected to diet treatments. In Helicoverpa armigera, however, ACTb exhibited the least stable expression in response to photoperiod and temperature treatments [40]. These results further confirmed that validating the stability of reference gene is very significant. The suitability of reference genes relative to both species and experimental conditions.
In addition to be affected by species and conditions [40], the ranking of reference gene stability is also affected by the tools used to perform the ranking. In the current study with B. odoriphaga, for example, the most stable genes across temperature treatments were RPS15, RPL28, and GAPDH by using NormFinder and ∆C t method but were RPL28, RPS15, and UBCE due to BestKeeper. This difference in ranking probably results from differences in the statistical algorithms: while BestKeeper individually analyzes the stability among candidate reference genes, NormFinder and the ∆C t method mainly think of the pair-wise variation between two candidate genes, and then confirm the stability of one of them [44,46]. Therefore, we used RefFinder software to comprehensively estimate the stability ranking of the 12 candidates. In addition, the optimal number of reference genes was confirmed by GeNorm, which calculates the pair-wise variation (V n /V n + 1 ) between the continuous standardization factors or NF (NF n and NF n + 1 ) [14] (Figure 3). If the first V value (V2/3) is <0.15, this indicates that two reference genes are enough for reliable normalization [14]. Nevertheless, the most appropriate number of reference genes also appears arbitrary without proper statistical verification under appropriate experimental condition. Some analyses, for example, failed to achieve V n / n + 1 <0.15, but could get relatively stable expression genes across final ranking estimated by GeNorm [47]. The most suitable number of reference genes conforms to the steadiest NF feasible with a unique sample set and a unique panel of candidates [48].
Random selection of reference genes reduces the accuracy of detecting interest genes expression because such a standardization strategy will be either under-estimate or over-estimate the expression differences among specimens. Such as the expression level of HSP23 among different temperature samples did not significantly differ using ACTb as internal control, but did significantly differ using other reference gene (such as RPS15) ( Figure 5). Normalization with two or more stable reference genes may be demanded, and researchers have recommended that multiple normalization genes were used to get more credible results [49][50][51]. Vandesompele et al. [14] recommended that reliable normalization needs at least three reference genes, and the pair-wise variation analysis in GeNorm hinted the need to include more than two genes in the current study. According to the ranking of expression stability among the 12 candidates evaluated by RefFinder in this work, we selected RPS15, RPL28, and GAPDH to assess the target gene HSP23 in B. odoriphaga under different temperatures; the results showed that HSP23 expression was up-regulated by low temperature, which was consistent with an earlier study that used RPS20 as reference gene [52]. In the current study, however, an arbitrarily selected reference gene (such as ACTb) failed to detect a significant effect of temperature on the expression profile of HSP23. Therefore, optimization of reference genes is critical for exact normalization of mRNA, especially for the subtle difference. To improve the accuracy of results, it is necessary to use the panel of selected housekeeping genes for any sample set.

Insects
B. odoriphaga was collected from a Chinese chive field on the Yang Town farm, ShunYi area (40˝1 1 N, 116˝6 1 E), Beijing, China. Individuals were reared for three generations with rhizomes of Chinese chive in an incubator (MLR-352H-PC) at 25˘1˝C, 70%˘5% relative humidity, and 12:12 (L:D). The specimens were promptly put into liquid nitrogen for further RNA isolation, and then screened following 12 candidate genes and amplification efficiencies.

Factors that Could Affect the Expression of Housekeeping Genes
The effects of the following factors on candidate reference genes mRNA were measured: developmental stage, temperature, population, pesticide exposure, diet, and photoperiod. After "exposure" to each factor (as described in the following sections), the specimens were placed in liquid nitrogen and then saved in´80˝C fridge for further study. Each factor was assessed in three independent experiments.

Temperatures
Groups of 20 4th-instar larvae were placed in a plastic Petri dish and exposed to 25, 4,´5, or´10˝C. After 4 h, they were exposed to 25˝C for another 24 h. Four living insects per group were then put in the tube (1.5-mL), frozen, and stored.

Pesticides
Groups of 40 4th-instar larvae were sprayed in culture dishes (φ = 60 mm) with the LC 90 value of allyl isothiocyanate, lime nitrogen, or thiamethoxam. An additional group of 40 larvae was sprayed with distilled water. After 24 h at 25˝C, four living larvae per group were saved in a 1.5-mL plastic tube, frozen, and stored.

Diets
Groups of four 4th-instar larvae were maintained in an incubator at 25˘1˝C, 70%˘5% relative humidity, and 12:12 (L:D) and were provided with one of the following: ginger slices, garlic bulbs, Chinese chive rhizomes, onion bulbs, or artificial diet [53]. After three generations, four larvae were placed into a 1.5-mL Eppendorf tube, frozen, and stored.

Populations
Larvae collected from three locations in China (Dezhou, Shandong; Baoding, Hebei; and Shunyi, Beijing) were reared on rhizomes of Chinese chive in an incubator at 25˘1˝C, 70˘5% relative humidity, and 12:12 (L:D). In the third generation, 12 4th-instar larvae from each population were placed in 1.5-mL micro centrifuge tubes (four larvae per tube), frozen, and stored.

Candidate Reference Genes
We assessed 12 "housekeeping" genes are known as reference genes selected from other insects. They were EF1a, UBCE, RSP5, GAPDH, RPS18, RPL18, ACTb, SDHA, RPL28, RPS13, RPS15, and TUB [33,34,36,40]. The sequences were obtained from our B. odoriphaga transcriptome data. The secondary structure of DNA template was predicted by the mfold web server [54], with the sets as follows: melting temperature for 60˝C; Na + concentration for 50 mM; Mg 2+ concentration for 3 mM; and linear DNA sequence. Other parameters were used as default. The primers used here were designed and checked by NCBI (National center for Biotechnology Information) Primer-BLAST, under the following conditions: primer GC content between 40% and 60%; primer melting temperature for 60˝C; and PCR products size of between 80 and 200 base pairs (Table 1).

Total RNA Abstraction and cDNA Synthesis
Total RNA was abstracted by the Trizol method. Each sample was homogenized with 1 mL of Trizol in a glass homogenizer following the manufacturer's protocol (TIANGEN, Beijing, China). The quality and quantity of RNA were assessed with a Thermo Scientific NanoDrop 2000c UV-Vis spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). The quality of the nucleic acid sample was considered good if the OD ratio (A 260 /A 280 ) was between 1.81 and 2.05. The cDNA was synthesized using the TransScript ® (TAKARA, Japan) All-in-One First-Strand cDNA Synthesis SuperMix in a 20 µL volume, with 4 µL 5ˆTransScript ® Buffer, 1 µg total RNA, and 1 µL gDNA Remover. Following the manufacturer's instruction, the 20-µL mixture was reacted in a Bio-rad PCR machine for 15 min at 42˝C before both the TransScript ® RT and gDNA remover were inactivated for 5 s at 85˝C. The cDNA was stored at´20˝C.

qRT-PCR
Each reaction was operated in a 20-µL solution including 0.4 µL cDNA, 10 µL 5ˆTransStart ® SuperMix, 0.4 µL forward primer, 0.4 µL reverse primer, and 0.4 µL 50ˆPassive Reference Dye. The amplification conditions for the qRT-PCR were set as follows: 94˝C for 30 s, followed by 40 cycles of 94˝C for 5 s, 60˝C for 15 s, and 72˝C for 34 s. Then, the 10-fold dilution series of cDNA was used for a standard curve. The melting curve analysis from 80 to 90˝C was used for assuring specificity of the amplified product [55]. The corresponding qRT-PCR efficiencies (E) were counted by means of the equation: E = (10 [´1/slope]´1 )ˆ100 [30,55].

Constancy of Gene Expression
The constancy of candidate genes was estimated by the ∆C t method [46] and with the following software: BestKeeper [56], GeNorm [14], and NormFinder [4]. The lower the value estimated by these algorithms, the greater the stability of expression. RefFinder [57], a useful web-based tool, was applied to estimate and screen the most suitable reference genes by combining the results of the four algorithms. Based on rankings from each algorithm, RefFinder assigned a suitable weight to each gene and counted the geometric mean of the overall ultimate ranking.

Evaluation of a Target Gene Expression
To select the suitable reference genes from 12 candidates, we estimated latent up-or down-regulation of the HSP23 gene in B. odoriphaga under different temperature treatments. Gene expression ratios were calculated by using the formula (2´∆ ∆Ct ) [58].
∆C t " C t ptarget geneq´C t preference geneq ∆∆C t " ∆C t psampleq´∆C t pcontrolq

Statistical Analysis
Results are showed as means˘SD. The means were calculated with Tukey's test at p < 0.05 by the software SPSS 19.0 for Windows (SPSS Inc., Chicago, IL, USA).

Conclusions
In summary, we first systematically evaluated 12 candidate reference genes in B. odoriphaga under various conditions. Four algorithms (NormFinder, BestKeeper, GeNorm, and the comparative ∆C t method) were used for evaluating the suitable reference genes. RefFinder, which was applied to combine the results of the different algorithms, then indicated that the most suitable reference genes were RPS15, RPL18, and RPS18 across developmental phases; RPS15, RPL28, and GAPDH across temperatures; RPS15 and RPL18 across pesticide treatments; RSP5, RPS18, and SDHA across photoperiods; ACTb, RPS18, and RPS15 across diets; RPS13 and RPL28 across populations; and RPS15, ACTb, and RPS18 across all samples. The use of the best reference genes vs. an arbitrarily selected reference gene resulted in substantial differences in the estimation of expression of a target gene. The results of this study will be valuable for research concerning gene function in B. odoriphaga.