Selection and evaluation of reference genes for expression analysis using quantitative real-time PCR in the Asian Ladybird Harmonia axyridis (Coleoptera: Coccinellidae)

Harmonia axyridis (Coleoptera: Coccinellidae) is a polyphagous insect that is an important biological agent used to control agricultural and forestry pests. The role of functional genes in H. axyridis based on quantitative real-time PCR (qRT-PCR) is increasingly well understood to investigate biology, physiology, feeding behavior and the role of important genes in physiological processes. Quantitative real-time PCR (qRT-PCR) is a powerful and reliable technique to quantify gene expression. Using qRT-PCR, expression levels of target genes are determined based on the levels of internal reference genes; therefore, reference genes need to be stably expressed under specific experimental conditions. However, there have been no studies on the stability of reference genes used in H. axyridis. In this study, we systematically investigated expression profiles of nine candidate reference genes from H. axyridis, including β-actin (ACTIN); elongation factor 1 α (EF1A); ribosomal proteins L10, L18, L28, S13, and S15 (RPL10, RPL18, RPL28, RPS13 and RPS15); glyceraldehyde-3-phosphate dehydrogenase (GAPDH); and superoxide dismutase (SOD). Four analytical methods (geNorm, BestKeeper, NormFinder, and the ΔCt method) were used to evaluate the suitability of these genes as internal reference genes for three biotic factors (developmental stage, tissue, and sex) and two abiotic treatments (temperature and photoperiod). RefFinder, a comprehensive evaluation platform integrating the four analytical methods, was used to rank the overall stability of these reference genes. Among the nine candidate genes, different reference genes were identified as having the most stable expression across biotic and abiotic factors. Genes encoding ribosomal proteins typically had the most stable expression, though EF1A was the most stable across developmental stages and photoperiods. To validate the suitability of these reference genes, heat shock protein 90 (HSP90) was chosen as a target. Significant up-regulation in HSP90 expression level in response to both low and high temperature was observed when using the most suitable reference genes but not when using an arbitrarily selected reference gene. The reference genes identified in this study will provide the basis for future functional genomics research in H. axyridis and will also facilitate the establishment of a standardized qRT-PCR program for other related insects.

Introduction ACTIN, EF1A, ribosomal proteins L10, L18, L28, S13, and S15 (RPL10, RPL18, RPL28, RPS13, and RPS15), GAPDH, and SOD from H. axyridis were tested. The effectiveness of these genes for expression normalization was further validated by qRT-PCR analysis of the well-studied target heat shock protein 90 (HSP90) gene. Many studies have demonstrated that heat-shock proteins act as molecular chaperones and play an important role in cellular responses to environmental stressors, including sublethal heat and cold shocks, infections, environmental contaminants [6,[27][28][29]. H. axyridis is the focus of agriculture and forestry production pest control strategies worldwide. And studying the temperature adaptability of H. axyridis is very important. Thus, it is of interest to study HSPs like HSP90 that may confer thermal tolerance. Based on our analysis of HSP90 and commonly used reference genes, we recommend specific combinations of reference genes to use for qRT-PCR analysis of different biotic and abiotic experimental conditions.

Insects
Harmonia axyridis (Coleoptera: Coccinellidae) was purchased from a commercial company in Beijing (Beijing Kuoye Tianyuan Biological Technology Co., Ltd., http://www.kuoye.com/). H. axyridis larvae and adults were fed with the aphid Aphis craccivora Koch (Hemiptera: Aphididae). Ladybeetles were reared in a growth chamber located at Beijing Academy of Agriculture and Forestry Sciences at a temperature of 23±1˚C, with a 16L:8D photoperiod and 70% relative humidity.

Factors
Effects of the following factors on reference gene expression were measured: development stage, tissue, sex, temperature, and photoperiod. The different developmental stages included eggs (30), first-instar larvae (15), second-instar larvae (10), third-instar larvae (2), fourth-instar larvae (2), pupae (1), and sex (one male and one female adult). The head, midgut, and carcass (body with head and viscera removed) were dissected from fourth instar larvae (10 for each replication) under a microscope (Invitrogen, Carlsbad, CA) and stored in TRIzol reagent. To determine the effect of sex on reference gene expression, one adult female and male were collected separately and placed in 1.5 ml centrifuge tubes. To examine the influence of temperature, groups of two third-instar larvae were separately exposed to 5˚C, 20˚C, and 35˚C for 3 h in a constant-temperature incubator (16L:8D photoperiod, 70±10% RH). To test the effect of photoperiod, groups of two third-instar larvae were kept in a constant-temperature incubator (23±1˚C, 70±10% RH) with a 16L:8D, 12L:12D, or 8L:16D photoperiod for 2 d. Each experiment was repeated three times. All samples were quickly frozen in liquid nitrogen after collection and stored at -80˚C in 1.5 ml centrifuge tubes for subsequent total RNA extraction.

Selection of gene sequences and primer design
We selected nine housekeeping genes and one target gene from our H. axyridis transcriptome data: ACTIN, EF1A, GAPDH, RPL10, RPL18, RPL28, RPS13, RPS15, SOD, and HSP90. The mfold web server (http://unafold.rna.albany.edu/?q=mfold/DNA-Folding-Form) was used to predict the secondary structure of the DNA template using the following settings: melting temperature, 60˚C; DNA sequence, linear; Na+ concentration, 50 mM; and Mg++ concentration, 3 mM. The default settings were used for the remaining parameters [30]. Primers were designed using NCBI Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/ index.cgi?LINK_LOC=BlastHome) with the following settings: GC content between 40-60%; melting temperature of 60˚C; and PCR product size between 100-200 base pairs. Excluded regions were defined based on the results of mfold analysis, and the default settings were used for the remaining parameters. PCR primer sequences used for quantification of the ten genes are shown in Table 1.

RNA extraction and cDNA synthesis
Total RNA was extracted using the Trizol method. Each sample was homogenized with 1 ml TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's protocol. The quality and quantity of RNA were assessed with a Thermo Scientific NanoDrop 2000 UV-Vis spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). The quality of the nucleic acid sample was considered good if the OD ratio (A260/A280) was between 1.8 and 2.05 [31]. Complementary DNA (cDNA) was synthesized from 1 μg total RNA using the Prime-Script™RT reagent Kit with gDNA Eraser (perfect Real Time) (TAKARA, Japan) according to the manufacturer's protocol. The cDNA was diluted 10-fold for subsequent qRT-PCR studies.

qRT-PCR
Real-time qPCR was carried out in 20 μl reactions containing 2.0 μl cDNA, 10 μl SYBR Premix Ex TaqTM II (Takara, Japan), 1 μl forward primer (10 μM), 1 μl reverse primer (10 μM), 0.4 μl Rox Reference Dye II (Takara, Japan) and 5.6 μl nuclease free water using an ABI PRISM 7500 Real-time PCR System (Applied Biosystems, USA). The amplification conditions for qRT-PCR were as follows: 95˚C for 30 s, followed by 40 cycles of 95˚C for 5 s, and 60˚C for 34 s. After all reactions, melting curve analysis (from 60 to 95˚C) was done to ensure consistency and specificity of the amplified product. A 10-fold dilution series of cDNA from the whole body of adults was used for a standard curve. The corresponding qRT-PCR efficiencies (E) were calculated using the equation: E = (10 [−1/slope] − 1) × 100 [32].

Constancy analysis of candidate reference genes
The constancy of candidate genes was evaluated using the ΔCT method [33], GeNorm [8], NormFinder [34], and BestKeeper [35]. We also used the online software RefFinder (http:// fulxie.0fees.us/?type=reference&ckattempt=1&i=1) to further evaluate the suitability of reference genes by analyzing the results of the four algorithms. Both geNorm and Normfinder require the conversion of raw cycle threshold (Ct) values into relative quantities. Initially, geNorm calculates a gene expression stability value (M) and then compares the pairwise variation (V) with other genes. Using microarray data as a training set for the algorithm, the pairwise variation between two sequential normalization factors, Vn/ Vn+1, is calculated. The optimal number of reference genes required for accurate normalization is determined based on a cutoff Vn/Vn+1 value of 0.15. NormFinder calculates gene expression stability for all samples in any number of groups based on intra-and inter-group variation and combines these values to provide a gene rank order based on the variation in gene expression. BestKeeper uses raw data (Ct values) and PCR efficiency (E) to compute bestsuited standards and combines them into an index. The comparative ΔCT method compares the relative expression of pairs of genes within each sample to select the optimal reference gene. Finally, RefFinder evaluates and screens reference genes by integrating the results of the above four major software programs.

Evaluation of target gene expression
The target gene H. axyridis HSP90 was used to evaluate the performance of nine candidate reference genes. We estimated up-or down-regulation of the HSP90 gene in H. axyridis across different temperatures. Relative expression of HSP90 was calculated using the formula (2 −ΔΔCT ) [36].

Evaluation of primer specificity and amplification efficiency
The primer specificity of nine candidate reference genes and one target gene (HSP90) was evaluated by PCR. Visualization of PCR products by 2% agarose gel electrophoresis revealed a single amplicon of expected size for each primer pair (S1 Fig). Furthermore, gene-specific amplification was confirmed by a single peak in real-time melting curve analysis (S2 Fig). A 5-point standard with known RNA standard concentrations was used to estimate the amplification efficiencies, which ranged from 94.7% to 106%. The coefficients of all 10 genes based on the linear regression were >0.990 (Table 1).

Expression profiling of reference genes
Mean Ct values of ten candidate genes of developmental stage, tissue, sex, temperature, photoperiod and total ranged from 17.

Stability of reference genes across biotic conditions
Developmental stage. The most stable genes based on developmental stage were RPL18, RPL28, and RPL10 according to geNorm. EF1A, RPS13, and RPS15 were the most stable according to Normfinder, BestKeeper, and the 4Ct method (Table 2).
Based on RefFinder, the ranking of reference genes from the most to the least stable was: EF1A, RPS13, RPL28, RPS15, RPL18, RPL10, SOD, ACTIN, and GAPDH (Fig 2). From GeNorm analysis, the pair-wise values of V2/3 and V3/4 were both above the cut-off value of 0.15 while the pair-wise value of V4/5 was <0.15 (Fig 3). A value <0. 15 indicates that supplemental reference genes are unnecessary. Based on RefFinder and convenient operation, EF1A, RPS13, and RPL28 were determined to be the best reference genes across different developmental stages (Table 3).
Tissue. Across tissue samples, RPL18, RPS13, and RPL28 were the most stable reference genes based on geNorm, Normfinder, BestKeeper, and the 4Ct method (Table 2), and the ranking of reference gene stability across different tissues based on RefFinder was: RPS13, RPL18, RPL28, EF1A, RPL10, RPS15, ACTIN, SOD, and GAPDH (Fig 2). From geNorm analysis, the pair-wise value of V3/4 was <0.15 (Fig 3). Thus, RPS13, RPL18, and RPL28 were considered the most suitable reference genes for comparisons across different tissues (Table 4). Evaluation of reference genes for qRT-PCR expression analysis using qRT-PCR in Harmonia axyridis Sex. According to geNorm, RPL18, RPL28, and EF1A were the most stable reference genes based on sex. Normfinder, BestKeeper and the 4Ct method all identified RPL18, RPL28, and ACTIN as the most stable genes across sexes (Table 2).

Stability of reference genes under abiotic conditions
Temperature. Based on geNorm, Normfinder, and the 4Ct method, RPL18, RPL28, and RPS15 were the most stable reference genes across temperature treatments. According to Best-Keeper, RPS15, EF1A, and RPL28 were the most stable (Table 3).
According to RefFinder, the ranking of reference genes from the most to least stable was: RPS15, RPL28, RPL18, SOD, EF1A, ACTIN, RPL10, RPS13, and GAPDH (Fig 2). The value of V2/3 from geNorm analysis was <0.15 (Fig 3). Therefore, RPS15 and RPL28 were considered the most stable reference genes across temperature treatments (Table 4). Photoperiod. EF1A, SOD, and RPL28 were the most stable reference genes across photoperiod treatments according to Normfinder and the 4Ct method; RPL18, RPL28, and SOD were the most stable according to geNorm; and SOD, ACTIN, EF1A, and RPL28 were the most stable according to BestKeeper (Table 3).

Target gene expression
To test the effect of reference genes on the calculation of target gene expression, the relative expression level of the target gene HSP90 under different temperature treatments (5, 20, 35˚C) was normalized to the most and least stable reference genes. The relative expression level of   Evaluation of reference genes for qRT-PCR expression analysis using qRT-PCR in Harmonia axyridis HSP90 significantly differed between temperature treatments. When the most stable reference gene, RPS15, was used, HSP90 was significantly up-regulated under low temperature (5˚C) and high temperature (35˚C), and the relative expression at 5˚C was significantly higher than at 35˚C (Fig 4). Similar results were obtained using RPS13 and RPL28. However, the expression level of HSP90 was not significantly different between the 5˚C and 35˚C samples when using an unstable reference gene (such as RPS13). Similarly, the expression of HSP90 was not significantly different between 35˚C and 20˚C samples when the least stable reference gene, GADPH, was used (Fig 4).

Discussion
In the current study, the stabilities of nine reference genes were evaluated across different biotic and abiotic conditions. We found that the best reference genes varied among conditions.  The expression level was normalized by different candidate reference genes: RPS15; RPS15 and RPL28; RPS13; and GAPDH. The reference genes were selected based on expression stability across three temperature treatments. Data represent the means±SEM of three biological replicates. The comparisons were analyzed using one-way ANOVA, and different letters (a,b) denote significant differences between normalization strategies determined by Tukey test (P < 0.05). https://doi.org/10.1371/journal.pone.0192521.g004 Evaluation of reference genes for qRT-PCR expression analysis using qRT-PCR in Harmonia axyridis The ribosomal proteins typically exhibited the highest stability across different experimental conditions. Across tissues, temperatures, and sexes, RPS13, RPS15, and RPL18, respectively, were the most stable reference genes. Consistent with our findings, ribosomal proteins have been previously reported to have the highest expression stability in several insects. For instance, in Bradysia odoriphaga, expression stability of RPS15 across temperature treatments was highest [31], and in Lucilia sericata, RPS3 and RPLP0 showed high stability across specific larval tissues [37]. in Tetranychus cinnabarinus, the highest expression stability across tissues was observed for RPS18 [38]. But expression stability was highest for RPS11 in Nilaparvata lugens [39], and for RPS13 in Sesamia inferens [40]. Our results and the results of previous studies indicate that in general ribosomal proteins are good reference genes for gene expression studies. However, an exception in our study was RPS13, which showed the least stable expression across sex and photoperiod treatment. RPS13 was also found to be the least stable reference gene across photoperiod treatments in B. odoriphaga [31].
In our study, among the nine reference genes, the expression level of EF1A was the most stable across developmental stages and photoperiod treatments. This result is consistent with previous studies concerning developmental stage in Frankliniella occidentalis [41] and developmental stage and photoperiod treatment in Danaus plexippus (L.) [42].
In our study, SOD ranked among the top three reference genes for photoperiod treatments (Table 4). In Spodoptera exigua, the expression stability of SOD was high across developmental stages and between sexest but was low across certain tissues [15].
GAPDH and ACTIN have been commonly used as internal controls in many gene expression studies [11,43,44]. In our study, however, ACTIN ranked low for most factors except for sex, where it was the third most stable gene.   [31] also found that ACTIN was a good reference gene in B. odoriphaga for diet treatments. In Helicoverpa armigera, however, ACTIN was the least stable reference gene across temperature and photoperiod treatments [45], and in Tribolium castaneum exposed to Beauveria bassiana, ACTIN was not stably expressed [46]. Based on our results, GAPDH is also not an ideal reference gene in H. axyridis. Several studies in other insects have also demonstrated that the expression stability of GAPDH is low in some circumstances, such as across developmental stages in Tetranychus cinnabarinus [19]; between the labial gland and fat body in Bombus terrestris and Bombus lucorum [47]; and across different body parts in Sogatella furcifera [48]. Taken together, these results suggest that the stability of reference genes may be dependent on the insect species and other characteristics such as instar stage or tissue type.
The ranking of reference genes is not only affected by different experimental conditions or factors, but also by the tools used for ranking. In this study, for example, the most stable reference genes under different temperature treatments were RPS15, RPL28, and RPL18 using NormFinder and the 4Ct method. BestKeeper also ranked RPS15 and RPS28 as the most stable, but EF1A was ranked the second most stable. The differences in ranking may result from differences in statistical algorithms. NormFinder and the 4Ct method mainly analyze the pairwise variation between reference genes, and then confirm the stability of one gene in each pair. BestKeeper individually analyzes the stability among reference genes [33,49]. RefFinder, a comprehensive evaluation platform that integrates the four algorithms we tested, was used to estimate the stability rankings of the nine reference genes. We also used GeNorm, which calculates the pairwise variation (Vn/Vn + 1) between the continuous standardization factors or NF (NFn and NFn + 1), to determine the optimal number of reference genes [8]. Based on GeNorm analysis, two reference genes were found to be sufficient for normalizing target gene expression values for sex, photoperiod, and temperature, but three reference genes were needed to normalize across developmental stages and different tissues [Fig 3]. These results suggest that it is necessary to use different combinations of reference genes to study changes in gene expression in H. axyridis in response to different factors.
In recent years, several reference genes have been used as internal controls for studying gene expression in H. axyridis under diverse experimental conditions. In this study, we found that nine candidate reference genes have different strengths and weaknesses in various conditions. For example, we recommend RPS15 and RPL28 to study gene expression in response to temperature in H. axyridis to temperature. Previously, Wang et al. [50] selected RP49 as a reference gene to study the expression of six small heat shock proteins mediating cold-hardiness in H. axyridis. Moreover, RP49 was also used as internal control to study the expression of genes in the H. axyridis trehalase and glycogen metabolic pathways [51,52]. Vilcinskas et al. [53] also used RPS3, another member of ribosomal protein family, as a reference gene in a study of genes encoding antimicrobial peptides and proteins, while study of immunity-related genes was exclusive of our setting conditions. In addition, Tang et al. [54] selected another 18s RNA to verify the expression of cold-resistance response genes, such as E3 ubiquitin-protein ligase, transketolase, trehalase, serine/arginine repetitive matrix protein 2, glycerol kinase and sugar transporter SWEET1-like. One of the most frequently used reference genes, ACTIN, was used as the internal control to study development-related genes in the ovaries of adult H. axyridis [55], and our study showed that this gene may not be as good as other reference genes for some factors such as developmental stage, tissue, temperature, and photoperiods.
Taken together our results indicate that commonly used reference genes are often not wellsuited for normalization in all qRT-PCR experiments, and the simultaneous measurement of a panel of candidate reference genes is critical for the accuracy of qRT-PCR quantification. The results of our study represent an important step for establishing a standardized gene analysis framework for H. axyridis.