Stable Reference Gene Selection for RT-qPCR Analysis in Nonviruliferous and Viruliferous Frankliniella occidentalis

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is a reliable technique for measuring and evaluating gene expression during variable biological processes. To facilitate gene expression studies, normalization of genes of interest relative to stable reference genes is crucial. The western flower thrips Frankliniella occidentalis (Pergande) (Thysanoptera: Thripidae), the main vector of tomato spotted wilt virus (TSWV), is a destructive invasive species. In this study, the expression profiles of 11 candidate reference genes from nonviruliferous and viruliferous F. occidentalis were investigated. Five distinct algorithms, geNorm, NormFinder, BestKeeper, the ΔC t method, and RefFinder, were used to determine the performance of these genes. geNorm, NormFinder, BestKeeper, and RefFinder identified heat shock protein 70 (HSP70), heat shock protein 60 (HSP60), elongation factor 1 α, and ribosomal protein l32 (RPL32) as the most stable reference genes, and the ΔC t method identified HSP60, HSP70, RPL32, and heat shock protein 90 as the most stable reference genes. Additionally, two reference genes were sufficient for reliable normalization in nonviruliferous and viruliferous F. occidentalis. This work provides a foundation for investigating the molecular mechanisms of TSWV and F. occidentalis interactions.


Introduction
The western flower thrips, Frankliniella occidentalis (Pergande) (Thysanoptera: Thripidae), is a destructive species that is found globally. This insect damages hundreds of plant species through direct and indirect mechanisms via feeding and transmitting tospoviruses, respectively [1][2][3]. The tomato spotted wilt virus (Family: Bunyaviridae; Genus: Tospovirus; TSWV) is transmitted in a circulative, propagative manner by thrips, with F. occidentalis as the most efficient vector of TSWV [4]. Currently, TSWV ranks among the top ten most economically important plant viruses worldwide [5]. TSWV infects in as many as 1,000 plant species, causing disease in many vegetables and ornamentals [6]. Global epidemics of TSWV are likely related to the worldwide distribution of F. occidentalis [7].
To better understand the molecular mechanisms of TSWV and F. occidentalis interactions, the transcriptomes of F. occidentalis were generated [8][9][10][11]. Additionally, RNA interference (RNAi) tool that enables functional genomics assays was successfully developed for F. occidentalis [12]. Given the nature of RNAi mechanisms, the impact of RNAi will likely be a method to control gene expression profiles of pest organisms. To make better use of these genomic resources, the establishment of a standardized reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) procedure in F. occidentalis according to the MIQE (Minimum Information for publication of Quantitative real-time PCR Experiments) guidelines is instrumental for functional genomics studies of this pest [13].
RT-qPCR is a high-throughput technique for measuring and evaluating gene expression [14]; however, there remain limitations that can significantly influence the normalization of gene expression, including sample amount, RNA quality and quantity, efficiency of reverse transcription, cDNA quality, PCR efficiency, and experimental operation between different samples [15,16]. Previously, Actin was used as stable reference gene to investigate the TSWV titers in F. occidentalis [12,[17][18][19], whereas another reference gene, RP49, was used as an internal control gene to estimate the abundance of TSWV N RNA in F. occidentalis [10]. The 18S gene was used as an endogenous control to investigate the expression profiles of 36 selected genes in viruliferous and non-viruliferous F. occidentalis [8], however, only one reference gene for the RT-qPCR experiment was used in these studies [8][9][10][11][17][18][19]. In standardized RT-qPCR experiments, two or several reference genes are commonly used to normalize gene expression data, which are stably expressed across various experimental conditions and serve as the internal controls [14][15][16]. Recently, seven candidate genes from F. occidentalis including Actin, 18S, H3, Tubulin, GAPDH, EF-1A, and RPL32 were evaluated for their suitability as reference genes across different developmental stages and temperatures [20]. A previous study strongly suggests the necessity of conducting custom reference gene selection designed for all experimental conditions, even when examining the same abiotic or biotic factor [21].

Ethics Statement
The western flower thrips F. occidentalis (Pergande) (Thysanoptera: Thripidae) was collected from clover plants, Trifolium repens L., at the Experimental Station of Qingdao Agricultural University. No specific permit was required for the collection described.
Insect rearing, plant cultures, and TSWV inoculation F. occidentalis was reared on the common bean Phaseolus vulgaris and maintained in MGC-250BP-2 incubators (Shanghai Yiheng Instruments, China) at 55-60% relative humidity and under a light: dark cycle of 16: 8 h. F. occidentalis adults were allowed to lay eggs on P. vulgaris for 1 day, and the adults were subsequently removed. First instar larvae were obtained in a few days.
Pepper (Capsicum annuum L., cv. Zhongjiao 6) plants were grown in soil mixed with vermiculite and organic fertilizer in 1.5-L pots (one plant per pot) under natural light and controlled temperatures (30 ± 2°C) in a greenhouse.
TSWV was maintained on Datura stramonium L. (Solanaceae). TSWV-infected pepper plants were prepared according to the method described in our previous study [22]. Control plants were treated similarly, but were inoculated by applying ground healthy plant material.

Nonviruliferous and viruliferous F. occidentalis
Healthy and TSWV-infected pepper leaf discs (diameter, 26 mm) were obtained using a cork borer. Each leaf was kept in a ventilated vial containing 3 ml of 1.0% agar to keep them fresh. Thirty newly hatched first instar larvae (nymphs) as one replicate were maintained on healthy or TSWV-infected discs in respective vials for 24 h to obtain nonviruliferous and viruliferous F. occidentalis, respectively [23]. Each treatment consisted of four replicates. Samples were frozen in liquid nitrogen and stored at -80°C.

Total RNA extraction and cDNA synthesis
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) according to previously described methods [24]. First-strand cDNA was synthesized from 1.0 μg of total RNA using a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan) according to the manufacturer's recommendations. cDNA was diluted 10-fold for subsequent RT-qPCR studies.

Reference gene selection and primer design
Eleven candidate reference genes were selected that are commonly used in RT-qPCR and that have been verified as stable genes in other species (Table 1). Primers were designed based on the sequences obtained from GenBank ( Table 1). The primers used for RT-qPCR analysis were designed online (https://www.idtdna.com/Primerquest/Home/Index) using the following parameters: amplicon length 75-150 bp with the optimum amplicon at 100 bp, Tm 59-65°C with the optimum amplicon at 62°C, primer lengths 17-30 bp, optimized to 22 bp, and GC content 35-65% with the optimal content at 50%. PCR amplifications were performed in 25-μl reactions containing 2.5 μl of 10× PCR Buffer (Mg 2+ Plus), 0.5 μl of dNTP mix (10 mM of each nucleotide), 0.5 μl of each primer (10 μM each), and 0.25 μl of TaKaRa Taq (5 U/μl) (TaKaRa). Target genes were amplified using the following parameters: initial denaturation at 94°C for 3 min; 35 cycles of 94°C for 30 s, 59°C for 45 s, and 72°C for 1 min; and a final elongation step of 72°C for 10 min. Primer specificity was confirmed by melting curve analysis and agarose gel electrophoresis of the amplification product.

Validation of selected reference genes
One TSWV-receptor gene of F. occidentalis was used to evaluate the candidate reference genes (GenBank No. AF247969). TSWV-receptor gene expression levels were investigated in nonviruliferous and viruliferous F. occidentalis. Two normalization factors (NFs) were calculated based on (1) the geometric mean of genes with the lowest Geomean values, and (2) a single reference with the lowest or highest Geomean value (as determined by RefFinder). Relative quantification of the TSWV-receptor gene in different samples was performed using the 2 -ΔΔCt method [25].

Data analysis
All biological replicates were used to calculate the average C t value. The stability of candidate reference genes was evaluated by the algorithms geNorm [14], NormFinder [15], BestKeeper [26], and the ΔC t method [27]. Finally, the tested candidates were compared and ranked using the web-based comprehensive analysis tool RefFinder (http://www.leonxie.com/referencegene. php). Nonparametric tests (K independent samples) were used to compare the expression

Transcriptional profiling of candidate reference genes
All tested genes were visualized as a single amplicon of the expected size on a 2.0% agarose gel (S1 Fig). Furthermore, gene-specific amplification was confirmed by a single peak in real-time melting curve analysis (S2 Fig). The linear regression equation, correlation coefficient, and PCR efficiency for each standard curve are shown in Table 1.
The mean and standard deviation (SD) of C t values were calculated for all samples (S1 Table). EF1A (SD = 0.68) had the least variable expression and this was reflected in its low SD values. By contrast, 28S (SD = 1.35) had the most variable expression, as shown by its high SD values. Additionally, 18S had the lowest C t values (C tavg = 10.09), suggesting that it had the highest expression level, whereas Actin was the lowest expressed gene among the candidates (C tavg = 31.62) (Fig 1; S1 Table).

Stability of candidate reference genes
GeNorm determines its ranking based on the geometric mean of the SD of each transformed gene pair combination (M-value). The lower the M-value, the higher the ranking. EF1A and HSP70 were co-ranked as the most stable genes (M = 0.270). The overall order based on geNorm from the most stable to the least stable reference gene was as follows: EF1A = HSP70, RPL32, HSP60, ATPase, HSP90, NADH, 28S, 18S, Actin, and Tubulin ( Table 2).
The ΔC t method relies on relative pair-wise comparisons. Using raw C t values, the average SD of each gene set is inversely proportional to its stability. HSP60 (0.69) was the top-ranked gene ( Table 2). The overall order from the most stable to the least stable reference gene based on the ΔC t method was as follows: HSP60, HSP70, RPL32, HSP90, EF1A, ATPase, NADH, 28S, 18S, Actin, and Tubulin (Table 2; S2 Table).
A low stability value suggests a more stable gene by NormFinder. HSP60 (0.152) was the most reliable and stable reference gene. The overall order from the most stable to the least stable reference gene based on NormFinder was as follows: HSP60, HSP70, RPL32, EF1A, ATPase, HSP90, NADH, 28S, 18S, Actin, and Tubulin (Table 2).
A low SD value suggests a more stable gene by BestKeeper. The stability of a gene is inversely proportional to the SD value. EF1A (SD = 0.499) had the least variable expression levels across all samples ( Table 2). The overall order from the most stable to the least stable reference gene based on BestKeeper was as follows: EF1A, HSP70, RPL32, HSP60, Actin, 18S, Tubulin, ATPase, NADH, HSP90, and 28S (Table 2).

Quantitative analysis of candidate reference genes based on geNorm
Under plant virus stress, the first V-value less than 0.15 was after V2/3 (Fig 2). This means that two reference genes were sufficient for reliable normalization regardless of the virus infection status of the insect.
Comprehensive ranking of the best reference genes using RefFinder According to RefFinder, which integrates the above-mentioned four software tools to compare and rank candidates, the comprehensive ranking of candidate reference genes from the most to the least stable was as follows: HSP70, HSP60, EF1A, RPL32, ATPase, HSP90, NADH, 18S, Actin, 28S, and Tubulin (Table 2). Among these, Tubulin had a geometric mean of almost 10.0, and was determined to be the least suitable candidate to serve as a reliable reference gene for normalizing gene expression.

Validation of selected reference genes
Using one, two, or three of the best reference gene combinations for normalization, expression of the TSWV-receptor gene did not differ between viruliferous and nonviruliferous F. occidentalis (P > 0.05) (Fig 3).

Discussion
Recent studies suggest that there is no single reference gene that is suitable for all types of normalization experiments [14,[28][29][30]. Therefore, candidate reference genes should be evaluated before they are used for normalization [21,31]. Recently, a study found that the ranking of candidate reference genes in the sweet potato whitefly, Bemisia tabaci varied substantially among intra-and inter-classes of insecticides [21]. This strongly suggested the need for conducting reference gene selection specifically designed for all experimental conditions, even when examining the same abiotic or biotic factor [21]. Zheng et al. (2014) focused on the selection of reference genes across different developmental stages and temperatures in F. occidentalis [20], while this study mainly focused on the impact of virus infection on the stability of these internal controls. As expected, these selected reference genes varied across different  [20]. Interestingly, Actin has been used previously as reference gene to investigate the TSWV titer in F. occidentalis [12,[17][18][19]. Therefore, we suggest that custom reference gene selection should be conducted for each experimental condition.
In previous studies, only one reference gene was selected for the RT-qPCR experiment in F. occidentalis [8][9][10][11][17][18][19]. To avoid biased normalization, many researchers have started to advocate the use of multiple reference genes to analyze gene expression [14,24,[28][29][30]. The geNorm program first calculates an expression stability value (M) for each gene and compares the pair-wise variation (V) of this gene with the other genes. A threshold of V<0.15 was suggested for valid normalization. geNorm starts with a single gene pair, and tests whether the inclusion of a 3 rd gene adds significant variation. The pair-wise variation (V n /V n+1 ) was analyzed between the normalization factors NF n and NF n+1 by geNorm to determine the optimal number of references genes required for qRT-PCR data normalization. In our study, the first pair-wise variation value less than 0.15 was after V2/3 (Fig 2), suggesting that two reference genes are sufficient for studying gene expression in nonviruliferous and viruliferous F. occidentalis (Fig 2).
In short, 11 candidate reference genes were selected for RT-qPCR analysis in nonviruliferous and viruliferous F. occidentalis assessed by five algorithms (geNorm, NormFinder, Best-Keeper, ΔC t method, and RefFinder). Among them, HSP60, HSP70, and RPL32 are the three most stable reference genes under the impact of plant virus infection as found by all five algorithms. Recently, the introduction of RNAi as a tool for functional genomics assays for F. occidentalis was developed [28]. Therefore, our study not only provides a standardized protocol for the quantification of gene expression in F. occidentalis, but also provides a solid foundation for genomic and functional genomics research assessing the interactions between TSWV and F. occidentalis.  Table. The mean and standard deviation (SD) of the C t value for each candidate reference gene. (DOCX) S2 Table. Summary of mean and SD values of gene pairwise comparison using the ΔC t method.