Evaluation and Validation of Reference Genes for qRT-PCR Normalization in Frankliniella occidentalis (Thysanoptera:Thripidae)

Quantitative real time PCR (qRT-PCR) has emerged as a reliable and reproducible technique for studying gene expression analysis. For accurate results, the normalization of data with reference genes is particularly essential. Once the transcriptome sequencing of Frankliniella occidentalis was completed, numerous unigenes were identified and annotated. Unfortunately, there are no studies on the stability of reference genes used in F. occidentalis. In this work, seven candidate reference genes, including actin, 18S rRNA, H3, tubulin, GAPDH, EF-1 and RPL32, were evaluated for their suitability as normalization genes under different experimental conditions using the statistical software programs BestKeeper, geNorm, Normfinder and the comparative ΔCt method. Because the rankings of the reference genes provided by each of the four programs were different, we chose a user-friendly web-based comprehensive tool RefFinder to get the final ranking. The result demonstrated that EF-1 and RPL32 displayed the most stable expression in different developmental stages; RPL32 and GAPDH showed the most stable expression at high temperatures, while 18S and EF-1 exhibited the most stable expression at low temperatures. In this study, we validated the suitable reference genes in F. occidentalis for gene expression profiling under different experimental conditions. The choice of internal standard is very important in the normalization of the target gene expression levels, thus validating and selecting the best genes will help improve the quality of gene expression data of F. occidentalis. What is more, these validated reference genes could serve as the basis for the selection of candidate reference genes in other insects.


Introduction
Quantitative measurements of gene expression are increasingly important for understanding biological processes. Knowledge of a gene's expression profile can provide evidence about its regulation and function [1]. Because of its accuracy, high sensitivity and reproducibility, quantitative real-time PCR (qRT-PCR) [2], is generally considered the best method for use in a variety of experimental and clinical conditions to analyze gene expression [3][4]. For accurate gene quantification analysis, normalization of qRT-PCR data is essential in order to eliminate variations in initial mRNA samples, sampling methods, reverse-transcription, and amplification efficiencies among different samples [5][6]. The most common normalization method is to relate mRNA levels of target genes to reference genes whose expression levels are highly stable in living organisms during various phases of development or under different environmental/experimental conditions [7].
Many genes involved in basic, ubiquitous cellular functions are housekeeping genes [8], such as actin, 18S rRNA, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) are commonly used as reference genes due to presumed uniform levels of expression. It is believed that there are few fluctuations in the transcription level of these genes compared to other genes [7,9]. However, many studies have shown that the expression levels of housekeeping genes do not always remain constant, thereby compromising their use as stable internal standards [10]. Moreover, errors in expression data of up to 20-fold have been demonstrated in instances where a single reference gene was used [11]. Therefore, it is critical to validate the expression stability of reference genes under different experimental conditions before using them for normalization.
The western flower thrips, Frankliniella occidentalis, is a notorious invasive species in many countries, causing serious economic damage to vegetables and ornamental plants [12]. It is both a direct pest of crops and an efficient vector of plant viruses such as tomato-spotted wilt virus [13]. In China, the western flower thrips was first recorded in Beijing in 2003 [14], and a recent investigation reported that it has spread to more than 10 provinces [15]. Sequencing of the western flower thrips genome was recently included in the 5,000 insect genome initiative, reflecting the rising importance of this invasive pest (http:// arthropodgenomes.org/wiki/i5K). Meanwhile, rapid progress has been made in the sequencing of western flower thrips by means of transcriptome analysis [16,17]. These data provide comprehensive gene expression information at the transcriptional level that may facilitate our understanding of the molecular mechanisms underlying various physiological aspects, including the development and response of western flower thrips to environmental stress [18].
To examine the changes of gene expression in F. occidentalis, Li et al. used 18S rRNA, the most frequently used reference gene in qRT-PCR analyses [19]. Because reference genes for F. occidentalis were selected without a companion validation study to evaluate their suitability under specific experimental conditions, we have evaluated the stability and performance of the following seven candidate reference genes: 18S ribosomal RNA (18S rRNA), b-actin (actin), tubulin, elongation factor-1 alpha (EF-1), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), histone 3 (H3) and ribosomal protein L32 (RPL32) under different experimental conditions (developmental stage, high temperature, and low temperature). The reference genes identified here will help future research to more precisely assess gene expression, to better understand the functionality of target genes and to elucidate specific molecular mechanisms underpinning the physiological processes in F. occidentalis.

Materials and Methods
Insects F. occidentalis was originally collected in Hangzhou, China in 2008 from a field maintained by our group and reared in the laboratory as described by Li et al. [20]. Individuals collected at the first day of each developmental stage were used in the experiment. We declare that no specific permissions were required for these activities and that the field studies did not involve endangered or protected species.
Samples from different developmental stages F. occidentalis from different developmental stages were collected separately and pooled by stage for gene analysis with more individuals pooled from the earlier life stages: 1st instar nymphs (300 per pool), 2nd instar nymphs (200 per pool), pupae (200 per pool) and adults (150 per pool). Four replicate pools for each stage were then collected for RNA extraction.

High temperature-induced stress
Each pool of 200 2nd instar nymphs were exposed to temperatures of 31, 33, 35, 37, 39 and 41uC for 2 h in a glass tube placed in a metal bath. A pool of nymphs were held at 26uC as a control. Four samples from each temperature were then collected for RNA extraction.

Low temperature-induced stress
Each group of 200 2nd instar nymphs were exposed to temperatures of 2, 0, 22, 24, 26, 28 and 210uC for 2 h in a glass tube placed in a metal bath. A pool of nymphs were held at 26uC as a control. Four samples from each temperature were then collected for RNA extraction.

Selection of gene sequences and primer design
Primer sequences and the associated amplicons characteristics for seven commonly used reference genes are summarized in Table 1. Based on reference genes described in the literature, the NCBI database (http://www.ncbi.nlm.nih.gov) was searched for homologous F. occidentalis sequences: 18S, actin, EF-1, RPL32, tubulin. And GAPDH and H3 sequences were amplified based on homologies with the hemipteran insects. All gene sequences have been submitted to GenBank. The Primer Premier 5 software (http://www.premierbiosoft.com/primerdesign/index.html) was used to design the primers. The parameters in Primer 5 were set as follows: amplicon length 70-150 bp, melting temperature 58-62uC, primer lengths 20-26 bp, and GC content 40-60%.

Quantitative Real-time PCR analysis
Total RNA was extracted using the SV Total RNA isolation system (Promega Z3100 (Promega, USA)), followed by DNase treatment to eliminate DNA contamination. The integrity of the RNA in all samples was verified by comparing the ribosomal RNA bands in ethidium bromide-stained gels. RNA sample purity was estimated using spectrophotometric measurements at 260 and 280 nm (Eppendorf Biophotometer plus (Eppendorf, Germany)). The OD260/280 of all samples was 1.8-2.2. Real-time PCR reactions were performed in a 20 ml total reaction volume comprised of 10 ml of iTaq universal SYBR Green supermix(2x) (Bio-Rad Laboratories Inc., USA), 1 ml of each gene specific primer (Table 1), 2 ml of cDNA templates, and 6 ml of PCR-grade water. Reactions were carried out on a CFX-96 real-time PCR system (Bio-Rad Laboratories Inc., USA). The efficiencies of the target and reference genes were similar [5,21]. The PCR parameters for all genes were as follows: 95uC for 3 min, 40 cycles of 95uC for 5 sec, 30 sec at the Tm value of primer pairs ( Table 1) with melting curve analysis performed to determine the specificity of PCR products. Every treatment included four replicates, and each reaction was run in triplicate.

Evaluation of target gene expression
Hop, also known as stress-induced protein-1, is a co-chaperone that usually acts as an adaptor, mediating the association of molecular chaperones Hsp90 and Hsp70 [19]. Hop was used as a target gene to evaluate the candidate reference genes. Relative quantification of Hop in different samples was conducted according to threshold cycle (Ct) value based on 2 2DDCt method.

Statistical analysis
Expression levels were determined as the number of cycles needed for the amplification to reach a fixed threshold in the exponential phase of the PCR reaction [22]. This threshold was set at 400 for all genes to determine the Ct values. The stability of candidate genes was evaluated by three commonly used software tools: BestKeeper (http://www.wzw.tum.de/genequantification/ bestkeeper.html) [23,24], geNorm version 3.5 (http://medgen. ugent.be/,jvdesomp/genorm/index.php) [8], and NormFinder-0953 (http://www.mdl.dk/publicationsnormfinder.htm) [25]. Finally, we compared and ranked the tested candidates based on a web-based analysis tool, RefFinder (http://www.leonxie.com/ referencegene.php?type=reference), including the comparative DCt method, to compare and rank the tested candidate reference genes [26]. The Bestkeeper index is based on the raw data and PCR amplification efficiency to determine the best-suited standards and combines them to create an index. Lower index scores denote greater transcriptional stability and thus better suitability as a reference gene. Ct values were converted into relative quantities and imported into the geNorm and NormFinder software programs. The geNorm algorithm first calculates an expression stability value (M) for each gene and then compares the pair-wise variation (V) of this gene with the others. The value of Vn/Vn+1 indicates the pairwise variation between two sequential normalization factors and determines the optimal number of reference genes required for accurate normalization [27]. Using microarray data as a training set for the algorithm, a threshold of V ,0.15 was suggested for valid normalization [8]. NormFinder uses a model-based approach to estimate expression variation in the selection of suitable reference genes [25]. Genes with the lowest values are the most stable. Based on the rankings from each program, RefFinder assigns an appropriate weight to an individual gene and calculates the geometric mean of their weights for the overall final ranking. Multiple comparisons of Ct values were performed by analysis of variance (One-Way ANOVA) followed by Tukey s-b(k) at the 95% confidence level using PASW Statistics version 18 (IBM Corp., Somers, NY).

Total RNA Quality and PCR Amplification Efficiencies
The concentration and purity levels of total RNA isolated from different samples were determined using the Eppendorf Biophotometer plus. The A260/A280 ratios ranged from 1.80 to 2.20 for RNA samples, indicating a high purity of total RNA for all samples. The integrity of all total RNA samples was confirmed using 1.0% agarose gel electrophoresis. For each of the primer pairs, single peak qRT-PCR melting curves suggested that each of the primer pairs amplified a unique product. A standard curve was generated for each gene, using a ten-fold serial dilution of the pooled cDNAs. The correlation coefficient and PCR efficiency for each standard curve are shown in Table 1. The PCR efficiency of the seven candidate reference genes and one target gene were similar, ranging from the lowest for EF-1 (95.4%) to the highest for 18S (108.9%). Linear regression coefficients (R 2 ) for all eight genes were $0.990. The melting temperatures (Tm) of all PCR products ranged from 55.0uC for RPL32 to 60.6uC for tubulin.

Expression profile of candidate reference genes
Expression levels were determined as the number of cycles needed for amplification to reach a fixed threshold in the exponential phase of the qPCR reaction [21]. Except for 18S, the gene expression analysis of all seven candidate reference genes displayed a narrow range of mean Ct values across all experimental samples (Fig. 1)

Analysis of gene expression stability
Developmental stages. The stability rankings generated by the comparative DCt method were largely similar with the results obtained from geNorm. However, the most stable genes as ranked by BestKeeper and NormFinder analysis were different from the results generated by geNorm and the comparative DCt method. The geNorm and comparative DCt method identified actin as the least stable gene and EF-1 as the most stable gene. In contrast, BestKeeper and NormFinder identified EF-1 as the second most stable gene. According to the results of RefFinder, the stability ranking from the most to the least stable in the developmental stages was EF-1.RPL32. tubulin .18S.GAPDH.H3.actin (Fig. 2). However, GeNorm analysis revealed that all the pair-wise variation values were above the proposed 0.15 cut-off (Fig. 3). These results indicate that normalization with three stable reference genes (EF-1, RPL32, tubulin ) was required (as suggested by the geNorm manual).
High temperature. Most of the programs with the exception of BestKeeper identified RPL32 and GAPDH as the most stable genes and actin as the least stable gene. However, RPL32 was also idenitified as the most stable gene, but GAPDH ranked fifth in BestKeeper ( Table 2). According to the results of RefFinder, the stability ranking from the most stable to the least stable gene in high temperature-stressed samples was RPL32.GAPDH. tubulin .EF-1.18S.actin.H3 (Fig. 2). GeNorm analysis revealed that all the pairwise variation values were below the proposed 0.15 cut-off value (Fig. 3). According to geNorm, two reference genes (RPL32 and GAPDH) should be required for a suitable normalization in the high temperature-stressed samples (Fig. 2).
Low temperature. The overall stability ranking generated by the comparative DCt method was largely similar with the results obtained by NormFinder and BestKeeper. All four programs identified 18S as the most stable gene and actin as the least stable gene in low temperature-stressed samples ( Table 2). RefFinder analysis ranked gene stability as 18S.EF-1.GAPDH.H3. tubulin .RPL32.actin. GeNorm analysis revealed that the pairwise variation values were below the proposed 0.15 cut-off value (Fig. 3). According to geNorm, two reference genes (18S and EF-1) should be required for a suitable normalization in the low temperature-stressed samples (Fig. 2).

Validation of reference gene selection
To assess the validity of selected reference genes, the relative expression of the target gene Hop was estimated for different experimental conditions. Target expression analyses further showed that differences in quantification were detected when normalizing with arbitrary reference genes relative to the best reference genes. We compared the mRNA transcript level of Hop when using two best reference genes RPL32 and GAPDH for high temperature-stressed samples, 18S and EF-1 for low temperaturestressed samples and the three best reference genes (EF-1, RPL32 and tubulin) for the different developmental stages along with the most unstable gene actin, which is usually recommended for normalization by RefFinder.
When using arbitrary genes to determine the relative expression levels of target genes in different samples the results can differ. For example, when using the most unstable gene for normalization, Hop transcript levels were higher in the nymphal stage compared to adult and pupae in the different development stages. However, after the best reference gene or the recommended two most stable references were used to normalize, no evident difference was detected (Fig. 4A). Similar results also occurred when calculating the relative expression levels of Hop after normalization with the unstable reference genes among the samples exposed to low temperatures (Fig. 4C). Therefore, it is important to determine the optimal reference genes for accurate normalization of qRT-PCR data, especially when differences in expression levels are subtle, as arbitrary selection of reference genes may decrease the accuracy of determining target gene expression. For example, the relative expression level of Hop showed no significant differences among the samples exposed to high temperatures when calculated using actin as the reference gene; however, the expression was significantly different when normalized by RPL32 and GAPDH (Fig. 4B). Therefore, in order to obtain accurate expression data, the expression stability of putative reference genes needs to be verified before each qRT-PCR experiment.

Discussion
Western flower thrips is considered one of the most economically important pests of agricultural crops worldwide [28,29]. Due to these thrips small size and secluded behavior, there is no effective way to control them. Previous studies concerning the transcriptome sequencing for F. occidentalis have been conducted, and numerous unigenes have been annotated according to their putative functional categories. While these sequence and putative functional data serve as a valuable resource for the control of this pest, it will be essential to validate the gene functions, including discovering their expression profile [16,17]. qRT-PCR is a Figure 2. Expression stability of the candidate reference genes as calculated by the Geomean method of RefFinder. A lower Geomean ranking indicates more stable expression. Expression stability of reference genes in the following samples: A) different developmental stages of Frankliniella occidentalis; B) F. occidentalis exposed to high temperatures; C) F. occidentalis exposed to low temperatures. doi:10.1371/journal.pone.0111369.g002 Validation of Reference Genes in Frankliniella occidentalis PLOS ONE | www.plosone.org sensitive, accurate and operable method to comparatively quantify transcript expression levels and profile gene expression levels [5,30,31]. Several studies on reference gene validation have emphasized that multiple internal genes must be evaluated in order to improve the accuracy of the qRT-PCR analysis and interpretation of gene expression [8,32,33]. Because few reference genes have been systematically evaluated in F. occidentalis, we identified suitable reference genes from different experimental conditions.
This study was conducted to identify the optimal reference genes for gene expression analyses of F. occidentalis from different developmental stages and temperature stress conditions. Seven candidate reference genes were analyzed by qRT-PCR. To our knowledge, this is the first study to evaluate the expression stability of different candidate reference genes for qRT-PCR in F. occidentalis. As expected, the most stable reference genes varied across the different experimental conditions (Fig. 2). Similar to Schistocerca gregaria [1], EF-1 and RPL32 displayed the most stable expression in development. RPL32 and GAPDH were the most stable under high temperature stress conditions, while 18S and EF-1 exhibited the most stable under low temperature stress conditions. These results indicate that the stability of reference genes expression in F. occidentalis needs to be determined for each experimental condition.
Previously, the structural protein actin, which is expressed at moderately abundant levels in most cell types, has been considered an ideal reference gene. actin was found to be one of the most suitable reference genes in studies of gene expression in Apis mellifera [34], Schistocerca gregaria [1], Drosophila melanogaster [35], Plutella xylostella and Chilo suppressalis [36], Chortoicetes terminifera [37] and Diuraphis noxia [38]. However, in our study, actin ranked as one of the least stable genes in almost all of the experimental conditions, actin is not a suitable reference gene for F. occidentalis under those experimental conditions. This result is consistent with previous studies that actin is not stable in Nilaparvata lugens [18] or Solenopsis invicta [39].
Although only a single reference gene with high expression stability may be appropriate for normalization of gene expression under some experimental conditions, in most experimental conditions two or more reference genes are required for accurate and reliable results [8]. In this study, we determined the optimal number of reference genes under different experimental conditions as calculated by geNorm, and calculated the two sequential normalization factors (NF n and NF n +1) needed to determine the minimum number of reference genes. Vandesompele et al. [8] proposed 0.15 as a cut-off value for V, below this value the inclusion of additional control genes is not required. In our studies, V values in high temperatures and low temperatures were both lower than the threshold value (0.15); therefore, we suggest that the two best reference genes be used to obtain accurate and reliable results under these conditions. Since, V values for the developmental stages were higher than the threshold value (0.15), we suggest that the three best reference genes be used for these conditions as the 0.15 recommended by geNorm (2007) is a suggested starting point.
Previous reference gene validation studies combined high and low temperature treatments, collectively referring to them as temperature treatment [18,27,40,41]. Interestingly, we found that the gene rankings for high and low temperature were different (Fig. 2). The most stable reference genes under high temperaturestress conditions were RPL32 and GAPDH, whereas RPL32 was one of the least stable reference genes under low temperaturestress conditions (Fig. 2). In contrast, the most stable reference genes under low temperature stress were 18S and EF-1. In general, expression of insect genes may vary under different temperature stesses. For example, hsp19.7 and hsp20.7 of Spodoptera litura responded only to heat treatment, not cold stress, while hsp21.5 of Chilo suppressalis only responded to cold stress [42,43]. Therefore, we believe that it is more reasonable to validate reference genes separately for high and low temperature.
Although the rankings of the different reference genes provided by three of the programs (BestKeeper, NormFinder and geNorm) were similar, our results showed that the determination of the most appropriate reference gene can vary depending on the software used (Table 2). For the developmental stages, geNorm recommended the use of EF-1 and RPL32, NormFinder recommended tubulin and BestKeeper recommended 18S. To determine which was the best reference gene, we chose a user-friendly web-based comprehensive tool called RefFinder. Based on the rankings from the computational programs (geNorm, Normfinder, BestKeeper, and the comparative DCt method), it assigned an appropriate weight to an individual gene and calculated the geometric mean of their weights to generate an overall final ranking. To validate the RefFinder-determined results, the expression of Hop was investi-  The expression stability was also measured using the DCt method, BestKeeper, NormFinder, and geNorm and ranked from the most stable to the least stable.
doi:10.1371/journal.pone.0111369.t002 gated (Fig. 2). Previous studies have demonstrated that Hop is significantly upregulated under high temperature stress conditions [21]. Compared to the more suitable reference genes, the use of actin as an internal control resulted in erroneous expression profiles (Fig. 4). The quantification of Hop using different control genes showed that the choice of the internal standard is very important in the normalization of the target gene expression levels.

Conclusions
In conclusion, this work is the first study aimed at validating a set of candidate reference genes for gene expression normalization using qRT-PCR in F. occidentalis. In the present research, seven candidate reference genes were identified. We concluded that EF-1 and RPL32 are the most suitable reference genes for the analysis of developmental stages, RPL32 and GAPDH are the most suitable reference genes under high temperature stress, and 18S and EF-1 are the most suitable reference genes under low temperature stress. These results may constitute a starting point for selecting reference genes in the future for more accurate normalization under other experimental conditions in F. occidentalis.