Identification and Validation of Reference Genes for Quantitative Real-Time PCR in Drosophila suzukii (Diptera: Drosophilidae)

To accurately evaluate gene expression levels and obtain more accurate quantitative real-time RT-PCR (qRT-PCR) data, normalization relative to reliable reference gene(s) is required. Drosophila suzukii, is an invasive fruit pest native to East Asia, and recently invaded Europe and North America, the stability of its reference genes have not been previously investigated. In this study, ten candidate reference genes (RPL18, RPS3, AK, EF-1β, TBP, NADH, HSP22, GAPDH, Actin, α-Tubulin), were evaluated for their suitability as normalization genes under different biotic (developmental stage, tissue and population), and abiotic (photoperiod, temperature) conditions. The three statistical approaches (geNorm, NormFinder and BestKeeper) and one web-based comprehensive tool (RefFinder) were used to normalize analysis of the ten candidate reference genes identified α-Tubulin, TBP and AK as the most stable candidates, while HSP22 and Actin showed the lowest expression stability. We used three most stable genes (α-Tubulin, TBP and AK) and one unstably expressed gene to analyze the expression of P-glycoprotein in abamectin-resistant and sensitive strains, and the results were similar to reference genes α-Tubulin, TBP and AK, which show good stability, while the result of HSP22 has a certain bias. The three validated reference genes can be widely used for quantification of target gene expression with qRT-PCR technology in D.suzukii.


Introduction
Quantitative real-time PCR (qRT-PCR) is a fast, sensitive, repeatable and accurate method for quantifying gene transcript levels [1,2]. However, several studies have revealed that interpretation would produce appreciable errors owing to different samples, treatments, RNA extraction techniques, amplification efficiencies, etc. Therefore, it requires an appropriate candidate reference gene to normalize the target gene expression, and the reference genes must have stable expression existing features in different cell types and experimental conditions [3][4][5].
Spotted wing drosophila (SWD), Drosophila suzukii (Diptera: Drosophilidae), is one of a handful of species that live on undamaged ripening fruits, and using its peculiar serrated ovipositor to break the skin of fresh ripening fruits and lay eggs in it [6]. D. suzukii was first recorded in Eastern Asia, and its western invasion started in 2008 with synchronous outbreaks in North America and Europe, which seriously threatened the local fruit industry, and has inflicted huge financial losses [7]. At present, there are few research reports on D. suzukii, and so far no qRT-PCR studies on it have been reported, so there is especially important to validated reference gene(s) for it. So far, several housekeeping genes, such as ribosomal protein L18 (RPL18) [8], ribosomal protein S3 (RPS3) [9], arginine kinase (AK) [10], elongation factor 1 beta (EF-1b) [11], TATA binding protein (TBP) [12], NADH dehydrogenase (NADH), heat shock protein 22 (HSP22) [13], glyceraldehyde-3-phosphate dehydrogenase (GAPDH) [14], Actin [15] and a-Tubulin [16,17] have been used widely as reference genes for qRT-PCR analysis in Insecta and other species. However, based on previous studies, it will not have a single universal reference gene that is suited for all experimental conditions [18][19][20]. Therefore, it is especially important to select the reliable reference genes in qRT-PCR based transcriptome studies.
In order to select the most suitable reference gene(s) for gene expression quantification by qRT-PCR, we examined the expressions of ten candidate reference genes (RPL18, RPS3, AK, EF-1b, TBP, NADH, HSP22, GAPDH, Actin, a-Tubulin) in four developmental stages (egg, larvae, pupa, and adult), three body regions (head, thorax, and abdomen), three different populations, and two abiotic conditions (photoperiod, temperature) treated samples. The gene expression quantification data were assessed using three statistical approaches (geNorm, Norm-Finder and BestKeeper) to normalized analysis.

Insects
The laboratory strain of D. suzukii was collected in cherry fields at Tai'an in Shandong Province, China in 2012. The solution of sucrose was serving as a food and water source for adults, and larvae were reared on a banana medium, which also served as ovipositional substrate. The rearing process was maintained in a greenhouse at a temperature of 2561uC, 16: 8 h light: dark photoperiod and 70-80% relative humidity. The abamectinresistant strain was derived from the laboratory (abamectinsensitive) strain, which feed on medium containing progressively higher doses of abamectin for more than 18 months. The species are common agricultural pests and not included in the ''List of Protected Animals in China''. No specific permissions were required as these fields are experimental plots that belong to Shandong Academy of Agricultural Sciences, Jinan, Shandong in China.

Reference Gene Selection and Primer Design
Ten commonly used reference genes were selected, based on the described insect reference genes in literature, the Spotted Wing Fly Base (http://spottedwingflybase.oregonstate.edu) was searched for available D. suzukii sequences: RPL18, RPS3, AK, EF-1b, TBP, NADH, HSP22, GAPDH, Actin, a-Tubulin. Primer 5.0 (http:// www.premierbiosoft.com/) was used to design primers for qRT-PCR analysis. The gene characteristics and primer sequences were summarized in Table 1.

Biotic Factors
Developmental stages. Four developmental stages (egg, larvae, pupa, and adult) were collected from D. suzukii. Samples used comprised 300 eggs, 30 first-instar larvae, 30 second-instar larvae, 20 third-instar larvae, 20 pupae, 20 first-day male and female adults, 20 seven-day old male and female adults and 20 fifteen-day old male and female adults for each replication. All the samples were immediately frozen in liquid nitrogen and stored at 280uC.
Tissue. Three body regions (head, thorax, and abdomen) were obtained from seven-day old adults using dissection needle and a stereo microscope in PBS solution on ice.
Population. One laboratory D. suzukii strain and two field collected populations from Zhejiang and Yunnan provinces were used.

Abiotic Conditions
Photoperiod. Each group of 10 seven-day old adults was kept in glass chambers with photoperiods (L/D) of 8:16, 12:12, and 16:8 for 72 h, and then the live adults stored for RNA extraction.
Temperature. A total of 30 seven -day-old adults were collected and placed individually into glass tubes. The tubes were then placed in water bath at 15.0, 25.0, and 35.0uC for 1 h, and stored as described earlier.

Total RNA Isolation and cDNA Synthesis
Total RNA was isolated with Total RNA Kit II (Omega, USA), RNA quantity was evaluated using a microvolume spectrophotometer NanoDrop 2000 (Thermo, USA), and the integrity was checked with 1% agarose gel electrophoresis. First-strand cDNA was synthesized with the PrimeScript RT reagent Kit (TaKaRa, Japan), the synthesized cDNAs were stored at 220uC for further qRT-PCR.

Quantitative Real-time PCR
The synthesized cDNA was amplified by PCR in 10 mL reaction mixtures using a Light Cycler 480 system (Roche, USA) and SYBR Premix Ex Taq (Takara, Japan) with the following procedure: 94uC for 5 min, followed by 45 cycles of 94uC for 15 s, 60uC for 20 s, and 72uC for 15 s. After the amplifications, a melting curve analysis was performed in triplicate, and the results were averaged. A 10-fold dilution series of cDNA was used to create the standard curve, and the qRT-PCR efficiency was determined for each gene and each treatment used to convert the Ct-values into the relative quantities.

Reference Gene Validation
P-glycoprotein, an ATP-dependent drug-efflux pump transmembrane protein, which closely related to abamectin resistance in Drosophila [21]. P-glycoprotein (DS10_00005769) was selected as the target gene to validate the reference genes. (P-glycoprotein F: AGCGTCACGACAAGAGGA; R: ACCACCCACAACGAG-GAA). The relative expression level of the target gene was calculated with different normalization factors, the selected most/ least stable genes, and which was determined in abamectinresistant and sensitive strains.

Statistical Analysis
Three software programs, geNorm (version 3.4), NormFinder (version 0.953) and BestKeeper (version 1), were used to evaluate the stability of the ten candidate reference genes. The geNorm program calculates an expression stability value (M) for each gene, the lower the M values, and the more stable reference genes. And then compares the pair-wise variation (V) of this gene with the others, we used pairwise variation V n/n+1 to estimate the optimal number of reference genes. A value below 0.15 indicates that an additional reference gene will not significantly improve normalization [14]. NormFinder uses a model-based approach to evaluate the overall variation of the candidate reference genes, the lower values, the more stable reference genes [22]. BestKeeper uses the geometric mean of the Ct values of the candidate reference genes with a standard deviation (SD), the lower index scores, and the more stable reference genes [23].
RefFinder (http://www.leonxie.com/referencegene.php?type = reference) is a web-based comprehensive tool that was used to evaluate optimal reference genes by integrating the results of the previous three analyses. Finally, we selected the ideal reference genes determined by RefFinder to choose the most stable genes [24].

Total RNA Quality and PCR Amplification Efficiencies
The concentration and purity of RNA isolated from different samples were determined using the microvolume spectrophotometer NanoDrop 2000, the absorbance ratios at 260/280 above 1.80 for all RNA samples. The integrity of all total RNA samples was confirmed using 1.0% agarose gel electrophoresis.
For each reference gene, a melting curve with a single peak suggested that each of the primer pairs amplified a unique product ranging from 100 (RPS3) to 227 bp (NADH; Fig. S1). The primer efficiencies of each standard curve for PCR efficiency, which was generated from third-instar larvae of the laboratory strain, were ranging from the lowest for RPL18 (92.6%) to the highest for AK (109.4%). The linear regression coefficients were larger than 0.99 for candidate reference genes except Actin (0.9899; Table 1).
Altogether, these results confirmed that the selected primers accurately amplify candidate reference genes.

Expression Profiles of Reference Genes
The cycle threshold values analysis of candidate reference genes in all different samples ranged from 11.01 to 27

Stability of Reference Gene Expression
geNorm Analysis. The geNorm program was used to calculate the average expression stability values (M values). The M values of tested genes was lower, the stability expression was higher, and vice versa. Considering the data obtained from various biotic conditions. For developmental stages, a-Tubulin and TBP were found to be the two most stable genes with the lowest M values. Similarly, HSP22 was found to be the least stable gene with the highest M value (Fig. S2A). For tissue, the AK and RPL18 genes showed the greatest expression stability, while Actin showed the least expression stability (Fig. S2B). For population, a-Tubulin and TBP were expressed more stable than other candidate reference genes, while Actin was expressed less stably than all others (Fig. S2C). Under various abiotic conditions, for photoperiod, a-Tubulin and RPL18 were the most stably expressed genes, while HSP22 was the least stably expressed gene (Fig. S2D). For temperature, the ranking of reference gene stability among those with M values ,0.5 is RPS3, TBP and a-Tubulin (Fig. S2E). Based on data obtained with five biotic and abiotic factors, a-Tubulin and TBP may be suitable for gene expression analysis in this study (Fig. 2).
In addition, geNorm was used to calculate the optimal number of control genes for normalization. If the pairwise variation of Vn/ n+1,0.15, it is not necessary to use $ n+1 genes as internal controls. Accordingly, the pairwise variation of V2/3 was lower than 0.15 in most of the conditions, except different tissues (Fig. 3). This indicates that combined use of the two stable control genes would be suitable for this research.
NormFinder Analysis. As with the geNorm method, the gene with the lowest M value has the most stable expression, and the NormFinder analysis result was similar to the result of geNorm analysis. NormFinder indicated that a-Tubulin, TBP and RPL18 are the most stable reference genes, and Actin and HSP22 are the least stable expression for five biotic and abiotic conditions ( Table 2).
BestKeeper Analysis. The BestKeeper program is another software tool to calculate the stability of a candidate reference gene, based on the standard deviation (SD) of the Ct values. The most stable reference genes were identified based on having the lowest standard deviation (SD). In this study, the BestKeeper analyses indicated that TBP and a-Tubulin were the most stably expressed genes, and HSP22 was the least stably expressed genes for different developmental stages (Table 3). In different tissues samples, the most stable genes were Actin and a-Tubulin, while HSP22 had the highest SD of all of the selected genes. Among the different populations, the most stable reference genes were TBP and TBPa-Tubulin, with the lowest SD values. In the different photoperiod-treated models, the most stable genes were found to be TBP and AK, and the least stable gene was found to be HSP22. For temperature, TBP and NADH were expressed more stable than other candidate reference genes.
RefFinder Analysis. Based on the geometric mean (GM) of the rankings obtained from three complementary statistical approaches, the five most reliable reference genes were the same for biotic or abiotic conditions samples. However, their ranks were a little different between biotic and abiotic samples. a-Tubulin (GM = 1.32)was the preferred candidate in developmental stages. The most stable to the least stable in the different tissues samples was TBP, AK, a-Tubulin, NADH and EF-1b. Among the different populations, a-Tubulin (GM = 1.41) was the most stably expressed gene, the remaining ranked candidates were TBP, EF-1b, NADH and AK with GM values ranging from 2.54 to 4.66, respectively. For photoperiod, a-Tubulin (GM = 1.41) and TBP (GM = 2.52) were expressed more stable than other candidate reference genes. For temperature, the order of best reference genes was as follows: TBP (GM = 1.00), AK (GM = 1.54), a-Tubulin (GM = 2.42), EF-1b (GM = 2.61), and NADH (GM = 3.54) ( Table 4).

Reference Gene Validation
P-glycoprotein expression levels was higher in resistant than in susceptible strain had been previously found by Lou et al [21] using Western blotting and Vanadate-sensitive ATPase assay. Pglycoprotein gene was further assessed by using the reference genes a-Tubulin, TBP and AK in comparison to HSP22. The expression profiles of the P-glycoprotein gene obtained using these three stable reference genes showed a similar trend in resistant or susceptible strain, but when HSP22 was used as the reference for normalization, the expression profile of P-glycoprotein gene changed. At the same time, there was significant expression difference of P-glycoprotein, which was higher in resistant than in susceptible strain, when calculated with a-Tubulin, TBP and AK as internal control, further demonstrating that a-Tubulin, TBP and AK were adequate internal controls On the contrary, there was no significant difference of the target gene, when used the least stable gene HSP22 (Fig. 4). This indicated the questionable results would be to produce by using unstable reference gene.

Discussion
Due to its high sensitivity, accuracy, specificity and repeatability, qRT-PCR is the first choice tool for gene expression analysis. Gene expression can vary across different physiological stages or cell types, so an appropriate internal control gene is required. The ideal reference gene should have relatively stable expression in distinct samples and different experimental conditions. In this study, we used three algorithms (geNorm, NormFinder and BestKeeper) and one web-based comprehensive tool (RefFinder) to evaluate the stability of 10 candidate reference genes in D.suzukii in response to different biotic (developmental stage, tissue and population), and abiotic (photoperiod, temperature) conditions.
Considering the developmental stages, the most stable genes were a-Tubulin and TBP (geNorm, NormFinder and Best-Keeper). For different tissues, AK and RPL18 were identified as the two most stable genes by geNorm and NormFinder, while the  two most stable genes by BestKeeper analysis were a-Tubulin and Actin. Among different populations, TBP and a-Tubulin had constant expression in all data sets produced by the three algorithms. In the different photoperiod-treated models, a-Tubulin and RPL18 were identified more stable genes by geNorm and NormFinder, while TBP and AK were considered the optimal reference genes by BestKeeper. For temperature, the most stable genes were RPS3 and TBP by geNorm and NormFinder, while the two optimal genes by BestKeeper analysis were TBP and NADH. Based on the three major statistic algorithms, a-Tubulin and TBP genes had a good performance under different conditions, while HSP22 and Actin were the lowest expression stability genes according to the stability values. Our results indicate that the ranking of these reference genes were identical among treatments based on the analysis results of various software programs (geNorm, NormFinder and BestKeeper), these programs have different algorithms and different sensitivities toward coregulated reference genes [19]. Therefore, each experimental investigation should establish which at least three reference genes are the most appropriate for the specific conditions [14]. According to the results of RefFinder, a-Tubulin, TBP and AK were the optimal reference genes.
Tubulin is one of several members of a small family of globular proteins, and is the major building block of microtubules in almost all eukaryotic cells, exhibited the most stable expression in the different conditions. Moreover, Tubulin had also been proved to be the normalized reference protein for Western blotting in Drosophila [25][26][27]. In this study, it was found to be the most stably expressed gene in the different treatment. TBP composed of transcription factor IID with TBP-associated factors, and is frequently used as a reference gene. Bansal et al recommend the TBP as a suitable HKG for efficient normalization among treatments, tissues, and developmental stages of A. glycines [28]. AK is the major phosphagen kinase in invertebrate groups, which has rarely been used as a reference gene in previous studies. In Bombus terrestris, AK was identified the most stable gene in the labial gland and fat body [8], and we found which the most stable gene was in different tissue ( Figure S2). RPL18 is located in the cytoplasm, which belongs to the L18E family of ribosomal proteins. Mamidala et al found that it is the most stable gene Table 4. Ranking of candidate reference genes in order of their expression stability as calculated by RefFinder.   expression across all the tissue and developmental-stage in Cimex lectularius [8]. It is surprising that most of housekeeping genes performed poorly as reference genes in this study. Previously, Actin has been considered to be one of the best reference genes for assessing gene expression in many insects [29,30], and GAPDH was also frequently used as a reference gene [11,31]. However, they were not the best reference genes in the present analysis and several studies have demonstrated that the stability of GAPDH expression was low in certain conditions [32].These results further suggest that the expression stability of reference genes is affected by different experimental conditions. This study was conducted to identify the stable reference gene(s) for qRT-PCR analyses of D.suzukii for studies of different biotic (developmental stage, tissue and population), and abiotic (photoperiod, temperature) conditions. After comprehensive consideration, we suggested that a-Tubulin, TBP and AK were the optimal reference genes. To our knowledge, this is the first systematic study to validate a set of different candidate reference genes for qRT-PCR in D.suzukii. The validation of reference genes in our study provides new information that will be useful for the accurate elucidation of the expression profiles of target genes. Furthermore, this study may also be useful for RNAi-based functional study and RT-PCR techniques that require reference gene for normalization in D.suzukii. Figure S1 Primer positions and ampliconic sequences are used for qRT-PCR. The DNA sequences are shown from the 59 to 39 end, and the primer positions are underlined. The products were first amplified by underlined PCR and then sent to Invitrogen for sequencing. (TIF) Figure S2 geNorm analysis of the expression stability of the 10 reference genes. Average expression stability values (M) and ranking of the candidate reference genes as calculated by geNorm software. A lower average stability value indicates more stable expression. (TIF)