Evaluation of reference genes for quantitative real-time PCR normalization in the scarab beetle Holotrichia oblita

Quantitative real-time polymerase chain reaction (qPT-PCR) is commonly used to analyze gene expression, however, the accuracy of the normalized results is affected by the expression stability of reference genes. Holotrichia oblita (Coleoptera: Scarabaeidae) causes serious damage to crops. Reliable reference genes in H. oblita are needed for qRT-PCR analysis. Therefore, we evaluated 13 reference genes under biotic and abiotic conditions. RefFinder provided a comprehensive stability ranking, and geNorm suggested the optimal number of reference genes for normalization. RPL13a and RPL18 were the most suitable reference genes for developmental stages, tissues, and temperature treatments; RPL13a and RPS3 were the most suitable for pesticide and photoperiod treatments; RPS18 and RPL18 were the most suitable for the two sexes. We validated the normalized results using odorant-binding protein genes as target genes in different tissues. Compared with the selected suitable reference genes, the expression of OBP1 in antennae, abdomen, and wings, and OBP2 in antennae and wings were overestimated due to the instability of ACTb. These results identified several reliable reference genes in H. oblita for normalization, and are valuable for future molecular studies.


Introduction
Quantitative real-time polymerase chain reaction (qRT-PCR), based on fluorescent signal monitoring, is commonly used for quantitative analysis of genes [1][2][3]. In most molecular studies, such as RNA sequencing (RNA-Seq) or RNA interference (RNAi), qRT-PCR is required to confirm accurate transcript changes of the target genes [4]. The reliability of qRT-PCR results is influenced by the availability of the reference genes [5]. The minimum requirements for qRT-PCR indicate that the effectiveness of reference genes as internal controls must be verified by corresponding experimental design [6]. Housekeeping genes used as reference genes without experimental validation can lead to poor normalization. Shi et al. [7] found that significant differences in the expression of HSP23 in Bradysia odoriphaga could not be detected under different temperatures when using the reference gene ACTb. The  polymerase Ⅱ gene, which is not a housekeeping gene, had more stable expression than classical housekeeping genes in different human tissues [8]. An idealized reference gene with absolutely stable expression in all conditions could not be found. Expression levels can fluctuate because of various factors [7][8][9][10]. The use of inappropriate reference genes can bias quantitative results. Specific genes normalized by a single reference gene are not credible without proper validation [11,12]. The stability of reference genes under specific experimental conditions must be evaluated to prevent nonbiological variations or errors [13]. Holotrichia oblita Faldermann (Coleoptera: Scarabaeidae) is widely distributed in China and causes serious damage to crops, forests, and lawns [14,15]. H. oblita adults stay underground during the day, while flying, feeding, and mating occur at night. Larvae feed on plant roots and remain underground during their development. H. oblita has a broad host range, extended feeding period, and cryptic habits [16]. Chemical pesticides are often used to kill larvae, but these can lead to soil pollution. The use of pheromone-baited traps can reduce larval populations by trapping adults. qRT-PCR has been applied to quantitative studies of olfactory genes to increase understanding of the odor recognition mechanism of H. oblita [17][18][19][20][21]. The reference genes for H. oblita, however, have not been evaluated under corresponding experimental conditions. Previously, a single housekeeping gene was randomly used for normalization.
The stability and effectiveness of reference genes for H. oblita need to be systematically evaluated. We assessed 13 candidate reference genes that involved several factors. These included developmental stage, tissue, sex, temperature, pesticide treatment, and photoperiod. We used the geNorm, NormFinder, BestKeeper, ΔCt method, and RefFinder to identify the stability of candidates. We studied the expression profile of two tissue-specific genes across tissues, which were standardized by the selected reference genes and a commonly used reference gene.

Insects rearing
We collected H. oblita in Feixi County (117˚60 0 E, 31˚39 0 N), Hefei, Anhui, China, in May 2018. No permits were required for the described study, which complied with all relevant regulations. Adults of both sexes were caught from fields at night and then placed in plastic boxes (60 × 50 × 50 cm) [16]. The bottom of the box was covered with a 20-cm deep soil layer with a moisture content of 15−18%. About 200 adults were housed in each box and the sex ratio was approximately 1:1 [22]. The adults were fed with fresh elm leaves (Ulmus pumila). We collected eggs every week and placed them in the box with the same soil. After the eggs hatched, larvae were fed with slices of fresh potato. When the larvae reached the 2nd instar, each larva was put in a separate cup. All of the insects were reared in a walk-in chamber under constant conditions of 25 ± 1˚C with 60 ± 5% relative humidity and a 8:16 h (L:D) photoperiod.

Sample collection and treatment
We evaluated the candidate reference genes under the following settings: developmental stage, tissue, sex, temperature, pesticide treatment, and photoperiod. After being processed under each experimental condition, all samples were immediately put into liquid nitrogen and stored at -70˚C. Each treatment had three biological replicates.
Tissues. Male and female adults were dissected into six body parts (antenna, head without antenna, thorax, abdomen, leg, and wing) respectively. One hundred pairs of antenna were separated and pooled as one sample. Samples of the head without antenna, abdomen, and thorax were obtained from two individuals. Leg samples were obtained from 12 individuals and wing samples were obtained from 6 individuals.
Temperatures. The 3rd instar larvae of H.oblita were exposed to 4˚C, 10˚C, 20˚C, or 30˚C for 2 h. The surviving individuals at each temperature treatment were collected and frozen. There was one individual larva per sample.
Pesticide treatments. We used clothianidin and bifenthrin insecticides in this study. They were dissolved in acetone at 2000 mg/L and 100 mg/L to produce stock solutions. The stock solutions were diluted with deionized water and used to treat soil containing 3rd instar larvae and the tests were scored at 48 h for larval mortality. At 48 h, the LC 50 concentrations of clothianidin and bifenthrin were 44.668 and 0.875 mg per kg of soil, respectively. The 3rd instar larvae were then placed into soil with the LC 50 concentrations for 48 h at 25˚C. We collected the surviving larvae and used one individual for each sample.
Sexes. Three pairs of male and female adults were caught on the night of May 2, which is the local early emergence period. We used one individual for each sample.
Photoperiods. We used the laboratory-reared adults for the photoperiod experiments. Newly emerged H.oblita adults were immediately placed at five photoperiods, including 0:24 h (L:D), 6:18, 12:12, 16:8, and 24:0. One pair of mated adults was put into a transparent plastic box and exposed to one of the five photoperiods randomly. After 7 d, one pair of adults was taken from each photoperiod and constituted one sample.

Total RNA extraction and cDNA synthesis
We used the MiniBEST Universal RNA Extraction Kit (TaKaRa, Dalian, China) to extract total RNA for all of the noted samples after being ground in liquid nitrogen. The purity and concentration of each RNA sample were checked by NanoVue Plus (GE Company, Fairfield, CT, USA). We used the OD value at a 260/280 nm ratio between 1.85 and 2.10 was used for further cDNA synthesis. Total RNA (1 μg per sample) was reverse transcribed following the manufacturer's instructions for the PrimeScript 1 RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). We checked the cDNA concentrations on the CFX96 System using RPS6 as a reference and adjusted the Ct value of all samples to approximately 18. Then, the cDNA was stored at -20˚C until used.
The melting curve and standard curve were drawn to check the specificity and amplification efficiency. We generated the standard curve by a serial 10-fold dilution of cDNA and calculated the efficiency value (E) of all primers by the formula: E = (10 [−1/slope] −1)×100 [27,28].

qRT-PCR
Each amplification reaction (25 μL) contained 12.5 μL SYBR Premix (Takara Bio, Dalian, China), 2 μL cDNA, 1 μL of each primer (10 μM), and 8.5 μL ddH 2 O. According to the MIQE guidelines, we performed qRT-PCR on CFX96 System (Bio-Rad, Hercules, CA, USA) [6]. The qRT-PCR amplification conditions were set as follow: 95˚C for 30 s, followed by 45 cycles at 95˚C for 30 s, 58˚C for 30 s, and 72˚C for 30 s. We performed each treatment with three biological samples and each sample had two technical replicates.

Stability of candidate reference genes
We used five algorithms to evaluate the stability of the candidate reference genes: RefFinder [29], ΔCt method [30], BestKeeper [28], geNorm [31], and NormFinder [32]. The RefFinder program provided a proper weight for each gene and generated a comprehensive ranking. The geNorm determined the optimal reference gene number by calculating the pairwise variation (V n /V n+1 ). The stability measure (M value) of gene expression calculated by geNorm proposed 1.5 as a cut-off line. A value less than 1.5 meant that this reference gene was stably expressed. Generally, lower values calculated by these algorithms indicated higher stability.

Validation by two target genes
We used two odorant-binding protein genes (OBP1 and OBP2) to verify the stability of the reference genes [33]. We used the optimum single reference gene RPS13 (RefFinder), the best reference gene pair RPL13a/RPL18 (geNorm), and a normally used reference gene ACTb to  Table 1 provides descriptions of the gene name, designed primer pair, product length, and primer amplification efficiency of all genes. We calculated the amplification efficiency by the slope of the standard curves. The efficiency values ranged between 91.33% (TUBb) and 103.03% (SYN6), with all regression coefficient (R 2 ) values > 0.99. We evaluated the primer specificity of all genes by the melting curve. The melting temperatures ranged from 79.00˚C (RPL28) to 85.00˚C (TUBb) with a single sharp peak, which confirmed gene-specific amplification.

Cycle threshold (Ct) values and variations in candidate reference genes
For each candidate reference gene, we analyzed Ct values to reveal the level of transcription (Fig 1). Evaluation of candidate reference genes Biotic experimental conditions. As shown in Table 2, geNorm ranked RPL13a and RPL18 as the most stable reference genes for normalization across different development stages. The M values, provided by geNorm, of ACTb, TUBa, and UBC were greater than the critical value of 1.5, and they were considered unsuitable as reference genes (Fig 2A). RPL13a and RPS6 were recommended as the most stable reference genes by NormFinder, while in the BestKeeper ranking, they were RPS3 and RPL18 (Table 2). RPS6 and RPL13a had the greatest stability by the ΔCt method ( Table 2). The comprehensive ranking provided by RefFinder from the highest to lowest was as follows: RPL18, RPL13a, GAPDH, TUBb, SDHA, RPS6, RPS3, RPS18, RPL28, SYN6, UBC, ACTb, and TUBa. The pairwise values of V 2 /V 3 calculated by geNorm below 0.15 indicated that two reference genes were enough for accurate normalization (Fig 2B). Therefore, RPL18 and RPL13a were demonstrated to be the best reference genes across the developmental stages of H. oblita (Figs 2B and 3A). For tissue samples, RPL13a and RPL18 were considered suitable reference genes using geNorm, and the M values for all reference genes were below 1.5 except for TUBa (Fig 2A, Table 2). The most stable reference genes were RPL13a according to NormFinder and the ΔCt method (Table 2). RPS3 ranked first in BestKeeper was and ranked second in NormFinder ( Table 2). The general ranking, according to RefFinder, was as follows: RPL13a, RPL18, RPS3, RPS6, RPS18, RPL28, TUBb, GAPDH, UBC, SYN6, SDHA, ACTb, and TUBa (Fig 3B). On the basis of the pairwise results calculated by geNorm, RPL13a and RPL18 were selected for accurate normalization (Figs 2B and 3B).
For different sexes of adults, the top four of the stable ranking were the same using Norm-Finder and ΔCt (Table 2). RPS18, RPL13a, and RPL18 were recommended as the best reference genes by geNorm, NormFinder, and the ΔCt method; however, in BestKeeper, they were RPS3, TUBb, and RPS6 (Table 2). On the basis of the four statistical formulas, TUBa was confirmed as the most unstable gene under the three biotic factors ( Table 2). GeNorm indicated that TUBa and ACTb were the most unstable genes because the M values were greater than 1.5 (Fig 2A). As shown in Fig 3C, RefFinder provided the ranking as RPS18, RPL18, RPL13a, TUBb, RPS3, RPL28, SYN6, RPS6, SDHA, ACTb, UBC, GAPDH, and TUBa from most stable to most unstable. According to the pairwise value of geNorm, RPS18 and RPL18 were the most suitable reference genes across sexes (Figs 2B and 3C).
Abiotic experimental conditions. For different temperatures, RPL18 and RPL13a were ranked first and third, respectively, by BestKeeper and ΔCt, and they were the top two by geNorm (Table 3). However, NormFinder showed that SDHA was the most stable reference gene ( Table 3). The least stable five reference genes were the same as those calculated by the four algorithms with the sequence being RPL28, SYN6, UBC, ACTb, and TUBa (Table 3). According to the result of geNorm, the pairwise values from V 2 /V 3 to V 11 /V 12 were below the 0.15 cutoff line and all gene M values were below 1.5 except for TUBa (Fig 2A and 2B). RefFinder ranked RPL18 and RPL13a as the top two, so they were evaluated as the genes most suitable for temperature treatments (Figs 2B and 3D).
For pesticide treatments, RPL13a was ranked first by NormFinder and ΔCt, whereas in geNorm the best gene was RPL18 (Table 3). GeNorm results indicated that RPL18 and RPL28 were the best reference genes ( Table 3). According to geNorm, the results of pairwise values and M values were similar to those of the temperature treatments. ACTb and TUBa were demonstrated to be the most unstable genes using four analysis tools (Fig 2A). According to RefFinder's online tool, the integrated ranking was RPL13a, RPS3, RPL18, RPS18, RPL28, SYN6, RPS6, GAPDH, SDHA, UBC, TUBb, ACTb, and TUBa (Fig 3E). Combining the results of pairwise values by geNorm, we selected RPL13a and RPS3 to be the most reliable reference genes (Figs 2B and 3E).
For different photoperiods, the top two ranked genes were RPL13a and RPL18 by geNorm and ΔCt, whereas RPL13a and RPS3 were ranked as the top two by NormFinder, RPS18 and RPS3 were ranked as the top two by BestKeeper (Table 3). According to the four algorithms, ACTb, SYN6, and TUBa were the most unstable candidate reference genes, and their M values were all higher than 1.5 (Table 3, Fig 2A). In addition, geNorm provided the pairwise value that V n /V n+1 were all below the 0.15 cutoff, except for V 12 /V 13 , which indicated that two reference genes were sufficient across the photoperiod treatments (Fig 2B). RefFinder indicated the ranking as RPL13a, RPS3, RPL18, RPS18, RPS6, SDHA, GAPDH, RPL28, UBC, TUBb, ACTb, SYN6, and TUBa (Fig 3F). RPL13a and RPS3 were considered to be the most suitable reference genes across different photoperiods (Figs 2B and 3F).
Total samples. For all of the experimental samples, RPL13a and RPL18 were the best reference genes in the stable ranking based on geNorm and ΔCt methods (Table 4). RPL13a and RPS6 were ranked the top two by NormFinder, whereas they were RPS3 and RPS18 by Best-Keeper (Table 4). ACTb and TUBa had M values above 1.5 by geNorm analysis, but were ranked as least stable according to the four methods (Fig 2A). The general ranking generated by RefFinder was RPL13a, RPL18, RPS3, RPS18, RPS6, RPL28, TUBb, GAPDH, SDHA, UBC, SYN6, ACTb, and TUBa (Fig 3G). According to the results of the pairwise values, the two genes RPL13a and RPL18 were the best reference genes for all samples (Figs 2B and 3G).

Effect of reference gene selection
To evaluate the reliability of the chosen reference genes in H. oblita, we selected two odorantbinding protein genes (OBP1 and OBP2) as the target genes across different tissues (Fig 4). Despite its frequent use as a reference gene, ACTb ranked next to last among the 13 candidate genes in the different tissues ( Fig 3B). The expression profiles of the two target genes normalized with RPL13a were consistent with the results of RPL13a/RPL18 (Fig 4). When calculated with ACTb, however, the expression levels of OBP1 in the antennae, abdomen, and wings and OBP2 in the antennae and wings were overestimated (Fig 4). The unstable reference genes led to more than 70-fold (head without antenna) and 60-fold (wing) errors in the quantification of OBP1 and OBP2 (Fig 4). Compared with the normalization results of RPL13a and RPL13a/ RPL18, OBP1 expression in the abdomen normalized with ACTb increased by 2.83-fold and 3.14-fold, respectively.

Discussion
The reliability of qRT-PCR results is determined by components, including sample quality, primer specificity and suitable reference genes [11]. In this study, we obtained all of the samples from fresh insects that were well preserved after processing. We checked all RNA samples by A260/A280, and eliminated the DNA residues by the gDNA eraser in the reagent kit. The primers of the selected genes had good efficiencies, ranging from 91.33% to 103.03% (R 2 >0.99), and the corresponding melting curve contained a single sharp peak (Table 1). Our data suggested that RNA quality and primer design were adequate for qRT-PCR. The expression level of the reference gene should be within a reliable range [34]. The median Ct values were in the range of 17.23 and 25.22, which indicated that they conformed to the basic requirements for qualified reference genes (Fig 1). The evaluation of reference genes is an important procedure for data normalization [35][36][37]. Tribolium castaneum was the first coleopteran insect in which reference genes were validated [2,26]. Studies on the selection of reference genes have now been reported for several Coleoptera, including Agrilus planipennis [38], Mylabris cichorii [39], Harmonia axyridis [40], and Anomala corpulenta [18]. However, the best reference genes may vary significantly among different insect species. Also, the reference gene can be affected by biotic or abiotic factors, and no reference gene is likely to be suitable for all experiment treatments [41,42]. ACTb is an important component of the cytoskeleton and widely distributed in cells. It has been the most common reference gene used for expression analysis [9,43,44]. ACTb was selected as a standardized reference gene in H. oblita [17,[19][20][21]. Our results showed that RPL18 was the best reference gene across different tissues of H. oblita adults, but among the tested genes, ACTb was next to last in stability (Fig 3B). We used OBP1 and OBP2 to validate the effectiveness of the candidate reference genes. Because of the instability of ACTb, the normalized results of the target genes in antennae, abdomen, or wings were overestimated (Fig 4). ACTb expression was significantly lower in antennae and wings than in other tissues. Therefore, up to 70-fold (head without antenna) and 60-fold (wing) changes were calculated in the relative expression of OBP1 and OBP2. The results show that the inappropriate selection of reference genes could lead to substantial errors in the quantitation of the target genes. Several other studies have challenged the applicability of ACTb as a reference gene for qRT-PCR [45][46][47]. These earlier results, and those of the present study, indicate that this "classic" gene is variable and needs to be assessed before further use as a reference gene. Previous studies demonstrated that ribosomal protein genes exhibit high expression stability and should be considered as a source of candidate reference genes. In Cimex lectularius, RPL18 was the best reference gene among different developmental stages and tissues [48]. The expression of RPL32 was the best reference gene among developmental stages of Bactrocera (Tetradacus) minax [49]. RPL28 and RPS15 were the best reference genes under temperature stress [50]. After fungal infection, RPS3, RPS18, and RPL13a expressed stably in T. castaneum [2]. RPL10 expression was stable within populations in Spodoptera litura [23]. Ribosomal proteins are the main components of the ribosome, which is involved in cell metabolism and regulation [51]. Our results showed that RPL18 and RPL13a were always in the top three of the stability ranking, and they were the best pair at different temperatures, developmental stages, and tissues (Fig 3). In some studies, however, ribosomal protein genes had unstable expression [7,52]. Our study showed that the gene with the most variable expression in all of the experimental conditions was TUBa (Fig 3). The M value of TUBa, determined by geNorm, was above the cutoff value of 1.5 under all experimental conditions. This finding indicates that it is unsuitable as a reference gene for H. oblita (Fig 2B). TUBa was identified as an available reference gene, however, in Tetranychus cinnabarinus [53], Drosophila melanogaster [9], and Nilaparvata lugens [54].
We used five statistical programs to analyze the stability rankings. The ΔCt method analyzes stability by calculating the pairwise variation within the candidate genes [30]. Ct raw data was used by BestKeeper to determine the stability index by ranking the standard deviation (SD) and coefficient of variance (CV) of each candidate gene [27,28]. In contrast to BestKeeper, the Ct value should be transformed to linear relative expression for NormFinder and geNorm [31,32,55]. For each candidate gene, both programs measure the expression variation in pairs to provide a stability value. Because of the unique analytical methods, the statistical results of the best-ranked genes showed variation among the different tools [56]. RPS18 and RPL18 ranked in the top three most stable, across sexes, according to the three methods except for BestKeeper ( Table 2). Among the tested temperatures, RPL13a and RPL18 exhibited the great stability using ΔCt, BestKeeper, and geNorm, but in NormFinder, the most stable genes were SDHA, TUBb, and GADPH (Table 2). Despite these slight differences, the overall trends of the four methods were similar in the stability rankings. GeNorm analyzed whether the pairwise variation (V n /V n+1 ) was above 0.15 and determined the optimal number for reliable normalization. In this study, the V 2/3 value within all factors was below 0.15, which meant that two reference genes were sufficient as internal standards (Fig 2B). Combined with the comprehensive ranking of RefFinder, RPL13a and RPL18 were considered to be the best normalizing genes across all developmental stages, tissues, and temperature treatments; RPL13a and RPS3 across pesticide and photoperiod treatments; and RPS18 and RPL18 across sexes (Fig 3).

Conclusions
We investigated 13 candidate reference genes of H. oblita under abiotic and biotic stresses and used five algorithms to evaluate their expression stability. We found reasonable variation in stability ranking as a result of the different formulas. RefFinder synthesized the results of these algorithms and provided the final results. Combined with the optimal reference gene number provided by geNorm, RPL13a and RPL18 were the most suitable reference genes in developmental stages, tissues, and temperature treatments; RPL13a and RPS3 in pesticide and photoperiod exposure; and RPS18 and RPL18 between the two sexes. OBP1 and OBP2 were used as target genes to validate the difference between the selected optimum reference genes and a classic reference gene. The results identified suitable reference genes in H. oblita under various experimental conditions. These data will benefit future molecular studies.
Supporting information S1 Fig. Melting curve of 13 candidate reference genes and two target genes. (TIF)