Identification and Evaluation of Suitable Reference Genes for RT-qPCR Analysis in Hippodamia variegata (Coleoptera: Coccinellidae) Under Different Biotic and Abiotic Conditions

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is an accurate and convenient technique for quantifying expression levels of the target genes. Selection of the appropriate reference gene is of the vital importance for RT-qPCR analysis. Hippodamia variegata is one of the most important predatory natural enemies of aphids. Recently, transcriptome and genome sequencings of H. variegata facilitate the gene functional studies. However, there has been rare investigation on the detection of stably expressed reference genes in H. variegata. In the current study, by using five analytical tools (Delta Ct, geNorm, NormFinder, BestKeeper, and RefFinder), eight candidate reference genes, namely, Actin, EF1α, RPL7, RPL18, RPS23, Tubulin-α, Tubulin-β, and TufA, were evaluated under four experimental conditions including developmental stages, tissues, temperatures, and diets. As a result, a specific set of reference genes were recommended for each experimental condition. These findings will help to improve the accuracy and reliability of RT-qPCR data, and lay a foundation for further exploration on the gene function of H. variegata.


INTRODUCTION
The reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is one of the most sensitive, accurate, and convenient methods for detecting quantitative analysis of gene expression in biological samples (Kubista et al., 2006;Provenzano and Mocellin, 2007;Derveaux et al., 2010). For RT-qPCR studies, it is necessary to select appropriate reference genes to correct and standardize the expression level of target genes. This will reduce the impact of RNA quality and cDNA synthesis efficiency on the results. Therefore, the validity of candidate internal reference gene needs to be evaluated before conducting RT-qPCR assay.
Reference genes, also called housekeeping genes, are supposed to be stably expressed in organisms under various biotic or abiotic conditions (Thellin et al., 1999;de Jonge et al., 2007;Zhu et al., 2014). The ideal internal reference gene should be stably expressed in different types of tissues and in different treatments of the same tissue, and its expression level is not affected by any endogenous or exogenous factors (Brym et al., 2013;Janská et al., 2013). Nevertheless more and more researches proved that the internal reference genes show significant various expression levels in different experimental treatments such as development stages, tissues, cells, and temperatures (Sun et al., 2008;Huis et al., 2010). Therefore, it is necessary to identify reference genes for their expression under different experimental conditions (Guo et al., 2016;Wan et al., 2017;Renard et al., 2018). Common reference genes include actin, glyceraldehyde-3phosphate dehydrogenase (GAPDH), ribosomal protein S3 (RPS3), 18S ribosomal RNA (18S), 25S ribosomal RNA (25S), tubulin, translation elongation factor 1-alpha (EF1α), and others have been used extensively for RT-qPCR analysis (Lu et al., 2013(Lu et al., , 2015Yuan et al., 2014;Zhang et al., 2015;Ma et al., 2016;Yan et al., 2016;Gao et al., 2017;Hu et al., 2018;Yin et al., 2020).

Insect Rearing
The H. variegata were obtained from Langfang Experimental Station of Chinese Academy of Agricultural Sciences, Hebei Province, China. Larvae and adults were maintained in plastic containers (8 × 8 × 11 cm) with Myzus persicae as a food source. The colonies were reared under conditions of 25 ± 1 • C, 65 ± 10% relative humidity, and a 14:10-h light/dark photoperiod.
To evaluate the effects of temperature, third-instar larvae of H. variegata were exposed to 15, 25, and 35 • C for 5 h. In diet treatment groups, 23-instar larvae of H. variegata were fed 3-dayold larvae of Spodoptera frugiperda (FAW), nymphs and adults of M. persicae, and canola pollen separately after 12 h of starvation (Schuldiner-Harpaz and Coll, 2017a,b).
Each experiment was replicated three times independently. Samples were preserved in 1.5-ml centrifuge tubes and snap frozen immediately in liquid nitrogen before storage at -80 • C.

Total RNA Extraction and Reverse Transcription
Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, United States) following the manufacturer's instructions. RNA purity was checked on a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). Each sample of 2.0 µg total RNA was used to synthesize single-stranded cDNA using the FastQuant RT kit (Tiangen Biotech, Beijing Co., Ltd., Beijing, China).

Selection of Candidate Reference Genes
The reciprocal BLAST hits approach was used to screen reference genes from H. variegata transcriptome and genome data (unpublished data). The reference genes from other insect species in GenBank were downloaded from NCBI 1 queried individually to H. variegata transcriptome and genome using the TBLASTN program with a permissive E-value cutoff of 10 −5 to get the hits. And then, each of the queried hits was compared back against non-redundant database of NCBI by the BLASTX program (E < 10 −5 ) to determine whether the original sequence was one of the hits. All the selected reference genes are listed in Table 1.
The primers of eight candidate reference genes were designed using Beacon Designer TM 7.9 software (Premier Biosoft International, CA, United States) according to sequences obtained from our recently sequenced transcriptomes and genome for H. variegata. Primer parameters were as follows: optimal temperature 60 ± 2 • C, GC content 40-60%, and 18-24 bp in length. The amplification size was between 100 and 200 bp ( Table 1).
Gene-specific primers (Supplementary Table 1) were designed using Primer Premier 6.0 software (Premier Biosoft International, CA, United States) to amplify ORF of Orco gene. PCR amplifications were performed in 50-µl reactions containing 25 µl 2 × Phanta Max Master Mix, 2 µl of each primer (10 µM each), 4 µl first-strand cDNA, and 17 µl ddH 2 O. The PCR parameters were as follows: one cycle of 94 • C for 3 min; 35 cycles of 94 • C for 30 s, 55 • C for 30 s, and 72 • C for 1 min; a final cycle of 72 • C for 10 min. PCR products were purified and

RT-qPCR Analysis
The RT-qPCR measurements were performed on an ABI Prism 7500 system (Applied Biosystems, Carlsbad, CA, United States).

Analysis of the Stability of Reference Gene Expression
Five tools (Delta Ct, geNorm, NormFinder, BestKeeper, RefFinder) were used to evaluate the stability of each reference gene. The Delta Ct method was used to select the optimal reference gene (Silver et al., 2006). geNorm was used to calculate the normalization factor to determine the optimal number of reference genes (Vandesompele et al., 2002); the proposed Vn/Vn + 1 ratio below V ≤ 0.15 does not significantly affect the normalization. The NormFinder algorithm, a model-based method, was used to rank the candidate reference genes by the stability value, and the gene with the lowest value was considered as the most stable reference gene (Andersen et al., 2004). The BestKeeper was used to rank the candidate reference genes based on the SDs and the coefficients of variation (CV) values (Pfaffl et al., 2004). Finally, RefFinder assigned an appropriate weight of the four methods to a single gene and calculated the geometric average of their weights for the overall final ranking (Xie et al., 2012).

Validation of the Selected Reference Genes
The Orco is highly conserved among insect species and is essential for localization of ORs in olfactory sensory neuron (ORN) dendrites and reception of odor signals (Leal, 2013). Thus, Orco was used as the target gene to assess the stability of candidate reference genes. Based on the selected candidate reference gene set, Orco expression levels in different tissues were calculated. The primers are shown in Supplementary Table 1. The gene expression was measured in various tissues and normalized by the optimal reference genes (TufA and RPS23) and the least stable reference genes (Tubulin-α and Tubulin-β). RT-qPCR analysis was performed as described previously. Data were calculated using the 2 − Ct method (Schmittgen and Livak, 2008) and the expression level of target gene (Orco) has been normalized with the CT-value average of two best and least stable reference genes.
All the experiments were performed in three replications, and the results are expressed as means ± SD.

RT-qPCR Analysis
Before evaluating the applicability of the reference genes, the specificity and efficiency of PCR amplification should be first ascertained. PCR amplifications for each primer pair showed a single peak in melt curves (Supplementary Figure 1). The primer efficiency (E) ranged between 87.98 and 103.42% with linear regression coefficient (R 2 ) values of 0.9851-1 ( Table 1).
Combining the three aforementioned factors, the results showed that the designed primers can precisely amplify candidate reference genes. Ct-values of the eight candidate reference genes ranged from 13.2 to 29.6, covering all the experimental conditions (Figure 1). EF1α was the most abundant reference gene with the lowest Ct-value, followed by RPL18, RPL7, Tubulin-β, RPS23, TufA, Actin, and Tubulin-α.
For the tissue-specific experiment, geNorm, NormFinder, and Delta Ct identified TufA and RPS23 as the most suitable reference genes, and BestKeeper identified EF1α and RPL18 as the most suitable reference genes ( Table 2). Through the calculation of RefFinder, the overall ranking of the eight candidate internal reference genes was as follows: TufA > RPS23 > RPL18 > EF1α > RPL7 > Actin > Tubulinβ > Tubulin-α ( Figure 2B and Table 2). The geNorm pairwise variation analysis also revealed that the first V < 0.15 emerged at V2/3 (Figure 3). Thus, the best combination of reference genes for tissue-specific samples of H. variegata was TufA and RPS23.
According to the comprehensive ranking of RefFinder, under biological conditions, the most stable to the least stable candidate reference genes were RPL18, RPL7, TufA, RPS23, EF1α, Actin, Tubulin-β, and Tubulin-α ( Figure 2C and Table 2). The geNorm results identified an initial V < 0.015 at V2/3, indicating that only two references genes were needed to normalize the target gene data (Figure 3). Thus, the best combination of reference genes for tissue-specific samples of H. variegata was RPL18 and RPL7.

Stability of Candidate Reference Genes Under Abiotic Conditions
Under different temperature conditions, combining the results of four software analysis, the RefFinder rankings were as follows: EF1α, Tubulin-β, TufA, RPL18, Tubulin-α, RPL7, RPS23, and Actin ( Figure 2D and Table 2). Through geNorm pairwise variation analysis, the results showed that the two reference genes were sufficient to normalize the target genes in different temperatures (Figure 3). Therefore, for H. variegata at different temperatures, the best combination of reference genes is EF1α and Tubulin-β.
For the diet treatment, geNorm and Delta Ct identified EF1α and RPS23 as the most suitable reference genes, NormFinder identified EF1α and RPL7, whereas BestKeeper identified EF1α and RPL18 as the most suitable reference genes ( Table 2). Based on RefFinder analysis, the comprehensive reference genes are ranked from the most stable to the most unstable as follows: EF1α, RPS23, RPL7, RPL18, Tubulin-α, TufA, Tubulin-β, and Actin ( Figure 2E and Table 2). In addition, geNorm pairwise variation analysis indicated that two reference genes were sufficient to normalize the target gene in dietary treatments (Figure 3). Therefore, we suggest that the best combination of reference genes in dietary treatments is EF1α and RPS23.
According to comprehensive ranking of RefFinder, under abiotic conditions, the order of the most stable to unstable candidate reference genes was as follows: EF1α, RPL18, RPL7, TufA, RPS23, Tubulin-β, Tubulin-α, and Actin ( Figure 2F and Table 2). The geNorm results identified an initial V < 0.015 at V2/3, indicating that only two references genes were required to standardize the target gene data (Figure 3). Therefore, under different abiotic status, the best combination of reference genes of H. variegata was EF1α and RPL18.

Validation of Selected Reference Genes
When standardized with the two most stable reference genes (TufA and RPL23), the expression pattern of Orco in various tissues were similar. The expression level of Orco was as follows: head > leg > wing > abdomen > thorax. However, the expression level of Orco in the head was almost 600-fold higher than in the leg when the two most stable reference genes were used for normalization ( Figure 4A). In contrast, when the two  FIGURE 3 | Optimal number of reference genes required for accurate normalization of gene expression by geNorm. Based on geNorm analysis, average pairwise variations were calculated between the normalization factors NFn and NFn + 1. Values <0.15 indicate that n + 1 genes were not required for the normalization of gene expression. most unstable reference genes (Tubulin-α and Tubulin-β) were used for normalization, the Orco expression in the head was tremendously high. The calculated relative expression of Orco in head with two unstable genes was almost 10 times higher than that of normalized with two stable genes ( Figure 4B).

DISCUSSION
H. variegata is one of the most important predatory natural enemies of aphids in the Xinjiang Uygur Autonomous Region, the main cotton-producing area of China, which shows a strong control effect on cotton aphids (Feng et al., 2000). The genome and transcriptomes of both sexes of H. variegata have been sequenced recently (unpublished data). To further investigate the biological roles of any specific gene in H. variegata, quantitative gene expression through RT-qPCR is essential. Following the "Minimum Information for Publication of Quantitative Real-Time PCR Experiments" (MIQE) guideline, reference genes with high stability are necessary (Bustin et al., 2009). However, to our knowledge, the best or combination of best reference genes under various treatment conditions have not been identified in H. variegata so far.
In this study, BestKeeper (Pfaffl et al., 2004), geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and Delta values (Silver et al., 2006) were used to evaluate the expression stability of the candidate internal reference genes. The overall stability of selected reference genes were also evaluated using RefFinder (Xie et al., 2012), which is a web-based analysis tool that integrates all four major calculation programs. Because of the different algorithms, the stability levels derived from the four analysis tools may be different. Optimal number of reference genes could be determined using geNorm analysis by calculating the paired mutation value (Vn/n + 1). When Vn/n + 1 is <0.15, there is no need to add additional reference genes to improve accuracy (Vandesompele et al., 2002).
The expression stability of eight candidate reference genes showed that the most suitable candidate combinations of reference genes were as follows: RPL7 and EF1α for developmental stages, TufA and RPS23 for tissues, EF1α and Tubulin-β for temperature treatment, EF1α and RPS23 for diet treatment, RPL18 and RPL7 for biotic conditions, and EF1α and RPL18 for abiotic conditions. The EF1α, encoding a protein associated with translation elongation, is the most abundant protein in the cell and highly conserved in different species (Yin et al., 2020). Several research works have been reported that EF1α is the most stable reference gene in different group of insect species such as Lymantria dispar (Lepidoptera: Lymantriidae) (Yin et al., 2020), Cydia pommnella (Lepidoptera: Tortricidae) (Wei et al., 2020), Aphis gossypii (Hemiptera: Aphididae) (Ma et al., 2016), Locusta migratoria manilensis (Orthoptera: Acrididae) (Yang et al., 2014), and H. convergens (Pan et al., 2015). The ICG website 2 counts the top 10 reported internal reference genes, among which, EF1α is in the first rank (Sang et al., 2017). The RPL18 encodes a ribosomal protein that is a component of the 60S subunit and belongs to the L18E family of ribosomal proteins (Fagerberg et al., 2014). It has been demonstrated that under different experimental conditions, RPL18 is also one of the most stable reference genes in different group of insects such as Solenopsis invicta (Hymenoptera, Formicidae) (Cheng et al., 2013), Anastrepha obliqua (Diptera, Tephritidae) (Nakamura et al., 2016), and Cimex lectularius (Hemiptera, Cimicidae) (Mamidala et al., 2011). Consistently, RPL18 is stably expressed under different biotic and abiotic conditions and could be considered as a proper reference gene for normalization of RT-qPCR data in H. variegata. Similar to EF1α, RPL family genes are also ranked among the top five stable reference genes for RT-qPCR analysis by ICG website (Sang et al., 2017).
Ribosomal protein S family including RPS23 has also been reported as the most appropriate reference gene in insects (Fu et al., 2013) and other organisms (Yadav et al., 2012) and has been ranked seventh by the ICG website. In contrast, TufA is not a commonly used reference gene in insects (Lü et al., 2018b). However, based on our preliminary evaluation and main experimental data, it has been shown that TufA was stably expressed in different tissues and could be considered as a reference gene in H. variegata. Previous studies demonstrated that in different experimental conditions, TufA is one of the most stable reference genes in different species of bacteria (Savard and Roy, 2009;Jiang et al., 2019).
The Orco is highly conserved in all insects and is widely expressed in the most olfactory receptor neurons (ORNs) of olfactory sensilla, distributed on the head (mainly expressed in antennae and mouthparts), legs, wings, and abdomen (Jacquin-Joly and Merlin, 2004). To validate our obtained results, the two most stable reference genes and the two most unstable reference genes were applied for normalization of Orco. Using the two most stable internal reference genes to examine the expression pattern of Orco, the results were stable and almost similar in different tissues. On the other hand, by using the two most unstable internal reference genes, the relative expression level of Orco normalized with Tubulinβ was approximately double lower in the head than when normalized with Tubulin-α. A similar phenomenon was also observed when OBP20 was used as the target gene to verify the stability of the internal reference gene of Apolygus lucorum (Hemiptera, Miridae). The relative expression level of OBP20, normalized with the stable reference genes (RPL32 and RPL27), was approximately threefold higher in the head than when normalized with two unstable reference genes (ACT and βACT) (Luo et al., 2019). Under two normalization conditions, the expression of target genes in same tissue (head) was different, indicating that the inappropriate use of reference genes in any experimental condition may lead to inaccurate results. Therefore, selection and validation of the best reference genes are crucial to determine the accuracy of the expression pattern of different genes. Taken together, our findings provide a more rigorous method to normalize RT-qPCR data in H. variegata, which lays a foundation for the functional genomics research in the future.