Introduction
Gene expression analysis is an important component of biological research and quantitative real-time PCR (qRT-PCR) is widely used for gene expression analysis. It can measure slight changes in individual gene expression because of its large dynamic range, high sensitivity and good reproducibility (VanGuilder et al., Reference VanGuilder, Vrana and Freeman2008; Lin & Lai, Reference Lin and Lai2010; Shen et al., Reference Shen, Jiang, Wang and Wang2010a ; Toutges et al., Reference Toutges, Hartzer, Lord and Oppert2010). Although qPCR is one of the most effective methods for analysis of gene expression, sample quantity, variations in efficiency of RNA extraction, cDNA concentration, primer performance, PCR efficiency and experimental precision are all factors that can introduce error (Udvardi et al., Reference Udvardi, Czechowski and Scheible2008; Bustin et al., Reference Bustin, Benes, Garson, Hellemans, Huggett, Kubista, Mueller, Nolan, Pfaffl, Shipley, Vandesompele and Wittwer2009). The conventional use of a single gene for normalization can lead to relatively large errors in a significant proportion of samples (Vandesompele et al., Reference Vandesompele, De Preter, Pattyn, Poppe, Van Roy, De Paepe and Speleman2002). In recent studies concerning reference gene selection, a single classic housekeeping gene was found to be inadequate for normalizing expression data of the other target genes (Lu et al., Reference Lu, Yuan, Gao, Kang, Zhan, Wan and Li2013; Liang et al., Reference Liang, Guo, Zhou and Gao2014; Sun et al., Reference Sun, Lu, Tang and Du2015).
Much insect research has attempted to validate and assess reference genes under a variety of biotic and abiotic conditions. The studied include Drosophila melanogaster (Ponton et al., Reference Ponton, Chapuis, Pernice, Sword and Simpson2011), Plutella xylostella (Fu et al., Reference Fu, Xie, Zhang, Wang, Wu, Liu, Zhou, Zhou and Zhang2013), Bemisia tabaci (Liang et al., Reference Liang, Guo, Zhou and Gao2014), Tribolium castaneum (Lord et al., Reference Lord, Hartzer, Toutges and Oppert2010), Spodoptera litura (Lu et al., Reference Lu, Yuan, Gao, Kang, Zhan, Wan and Li2013), Sesamia inferens (Sun et al., Reference Sun, Lu, Tang and Du2015), Spodoptera exigua (Zhu et al., Reference Zhu, Yuan, Shakeel, Zhang, Wang, Wang, Zhan, Kang and Li2014), Bactrocera minax (Wang et al., Reference Wang, Zhao and Liu2014) and Helicoverpa armigera (Zhang et al., Reference Zhang, An, Li, Wu, Yang, Liu, Cao, Zhang, Zhang and Liu2015b ). No reference gene is universally applicable under all conditions. It is therefore necessary to evaluate the expression profiles of candidate reference genes for each specific experiment. The comparative ΔC t method (Silver et al., Reference Silver, Best, Jiang and Thein2006) and various computational programs (NormFinder (Andersen et al., Reference Andersen, Jensen and Ørntoft2004), geNorm (Vandesompele et al., Reference Vandesompele, De Preter, Pattyn, Poppe, Van Roy, De Paepe and Speleman2002), BestKeeper (Pfaffl et al., Reference Pfaffl, Tichopad, Prgomet and Neuvians2004)) have been used to identify the most stably expressed reference genes within a given set of biological samples. RefFinder, an online tool, can integrate these methods to compare reference genes and recommend the best-suited candidate reference genes (Xie et al., Reference Xie, Sun, Stiller and Zhang2011).
Ideal reference genes should not be regulated or influenced by the experimental procedure, or by different conditions. They should also have high expression rates and exhibit similar, stable mRNA expression levels under diverse treatments and in different tissues (Radonic et al., Reference Radonic, Thulke, Mackay, Landt, Siegert and Nitsche2004). Housekeeping genes, such as β-actin (ACTB), glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and 18S ribosomal RNA (18S rRNA), are commonly used as reference genes in different organisms. No reference genes are stably expressed and suitable for all cell and tissue types, and all experimental conditions (Derveaux et al., Reference Derveaux, Vandesompele and Hellemans2010).
Galeruca daurica (Joannis) (Coleoptera: Chrysomelidae) is a serious insect pest occurring in the Inner Mongolia grasslands. It causes significant damage to Allium spp. (Zhang et al., Reference Zhang, Zhou, Pang, Chang, Shan and Zhang2015a ). Moreover, the pest would lead to a devastating disaster in turf seeding growth on grassland, especially on barren pasture. The physiological and biological characteristics of this species are known (Gao et al., Reference Gao, Zhou, Pang, Bao, Luo and Erdeng2015), but the mechanisms underlying population outbreaks of G. daurica are unclear. Expression analysis of relevant genes may provide insight into G. daurica physiology and biology. Prior to expression analysis study, effective reference gene combinations in G. daurica are required.
We selected a set of commonly used housekeeping genes from other insect species as candidate reference genes for normalization of gene expression. Eight housekeeping genes (actin, GAPDH, succinate dehydrogenase (SDHA), ribosomal protein L32, tubulin-alpha (TUB-α), tubulin-beta (TUB- β), glutathione S-transferase (GST) and TATA-box), and two target genes (a P450 CYP6 gene and Hsp70 gene) were evaluated. The stability of these candidate genes was investigated under four biotic conditions (developmental stage, gender, tissue type and diapauses) and one abiotic condition (temperature). Our objective was to determine suitable reference genes in G. daurica and to evaluate the importance of variations in relative qualification under a range of conditions.
Materials and methods
Insects
G. daurica has one generation per year in its natural habitats. The G. daurica strain used in this study was originally collected in the Arxant village of Xilinhot City in early May in 2014, and it has been maintained on garlic chives for one generation in a growth chamber (23 ± 2°C, 40 ± 10% RH, L16:D8) in our laboratory.
Treatments
Temperature treatment
To examine temperature influence, second-instar larvae were exposed to −14, 0, 10, 20, 30 and 40°C for 1 h, and then returned back to 25°C for 30 min. They were then snap frozen and stored at −80°C until qPCR testing.
Development stages
Different developmental stages of G. daurica were collected and pooled as follows: eggs (150–200 per pool), first-instar larvae (80–100 per pool), second-instar larvae (8–10 per pool), third-instar larvae (3–5 per pool), pupae (3–5 per pool) and mixed sex adults (2–3 per pool). All samples were snap frozen in liquid nitrogen, and stored at −80°C until qPCR testing.
Gender effects
Three to five male and three to five female adults were collected separately and placed in separate centrifuge tubes, and then snap frozen and stored as described.
Tissue effects
Dissection of body parts (head, thorax and abdomen) from male and female G. daurica adults was done under a stereo microscope (OLYMPUS, LG-PS2-2, Tokyo, Japan). The dissected parts were snap frozen and stored as previously described.
Diapause and non-diapause adults
Mixed sex G. daurica adults (7 days old) were selected as non-diapause individuals. At 20–30 days, adult G. daurica had entered diapause. A group of 30 days old mixed sex G. daurica adults were collected and evaluated as diapause test insects. All the samples were snap frozen in liquid nitrogen, and stored at −80°C until qPCR test.
Note: For all the collections, the garlic chives must be removed for half an hour from all the test bugs before they are frozen to prevent plant retention in bugs’ digestive system.
Total RNA extraction and cDNA synthesis
Total RNA was extracted using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer instructions. RNA concentration and quality were measured according to the optical density at 260 nm and the A260/A280 absorption ratio using a Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA). RNA samples used had an A260/A280 ratio ranging from 1.8 to 2.0. After adjusting the samples to equivalent concentration, one microgram of RNA was reverse transcribed into first-strand cDNA using the PrimeScript 1st Strand cDNA Synthesis Kit (Takara, Japan). The cDNA was stored at −20°C until use.
Primer design and qRT-PCR
Eight candidate reference genes, including five commonly used reference genes (actin-β (ACT), TUB-α, TUB-βGAPDH and SDHA), as well as three infrequently used reference genes (ribosomal protein L32 (RPL32), GST- δ and TATA-box in G. daurica, were evaluated for their expression stability. The primers were designed using online Primer3Web software (http://primer3.ut.ee/). The sequences, length of products, and source of these candidate genes are listed in table 1.
Candidate reference genes and target gene primers used for qRT-PCR analysis.
qRT-PCR was performed using an ABI 7500 PCR system (Applied Biosystems, StepOnePlus, USA). The 25 µl reactions contained 1 µl of cDNA template, 12 µl of SYBR Green Real-time PCR Master Mix (Takara, No. DRR420S) and 0.5 µl of each primer, and nuclease-free water was added for a total of 25 µl. Reactions were carried out under the following conditions: 95°C for 5 min followed by 40 cycles of 95°C for 20 s, 52°C for 20 s and 72°C for 30 s. Each treatment included three replicates, and each reaction was run in triplicate. A 10-fold dilution series of cDNA was employed to construct a standard curve to determine the PCR efficiency. Corresponding qRT-PCR efficiencies (E) were calculated according to the equation: E = (10[−1/slope]−1) × 100 (Pfaffl, Reference Pfaffl2001).
Statistical analysis
All biological replicates were used to calculate the average C t value. Stability values of the eight candidate reference genes were assessed using geNorm, NormFinder, BestKeeper and the comparative ΔC t method. RefFinder, a user-friendly web-based software (http://www.leonxie.com/referencegene.php), was used to rank the expression stabilities of candidate reference genes (Pfaffl et al., Reference Pfaffl, Tichopad, Prgomet and Neuvians2004; Silver et al., Reference Silver, Best, Jiang and Thein2006).
The geNorm provides a measure of gene expression stability (M), and genes with the lowest M values have the most stable expression. An M value below 1.5 indicates that the candidate reference gene is stable and appropriate for use as a normalizer. The geNorm also performs pairwise comparisons of one selected gene to others, and calculates a serial value of Vn/Vn + 1. A value above 0.15 indicates that an additional reference gene could be added to improve normalization. NormFinder provides a stability value for each gene and ranks the stability of tested candidate reference genes. BestKeeper determines the SD with the user selecting the best genes based on these variables. RefFinder integrates the currently available major tools (geNorm, Normfinder, BestKeeper and the comparative ΔC t method), assigns an appropriate weight to each candidate gene, and calculates the geometric mean of their weights for the entire ranking.
Validation of reference gene selection
A putative insecticide resistance-associated gene Cytochrome P450 (CYP6) and a heat shock-related gene (Hsp70) were used to assess the validity of selected reference genes under different biotic and abiotic conditions. Both target genes were tested. The most stable gene (Normalization factor 1, NF1), the least stable gene (NF8) (as determined by RefFinder) and the combined set of reference genes (NF1–2 or NF1–3) (as calculated and recommended by geNorm) were used for comparative purposes. Relative expression was conducted based on the 2−ΔΔCt value method. One-way ANOVA test was used for statistical analysis, student's t-tests were performed to compare target gene expression calculated with three sets of reference genes with significance reported for P < 0.05.
Results
Amplification efficiency of the primers
All of the eight candidate reference genes and two target genes were expressed in all G. daurica sample sets, as visualized by the presence of a single band of the expected size on a 1% agarose gel. A standard curve was generated for each gene, using 10-fold serial dilutions of the pooled cDNA generated from each experiment. All amplification efficiencies in the qRT-PCR analysis for the eight candidate genes and three target genes ranged from the lowest for TATA-box (90.661%) to the highest from SDHA (114.108%). Linear regression coefficients (r 2) for all ten genes and the slopes of the standard curve are shown in table 1.
Expression profiles of candidate reference genes
Expression analyses of the eight reference genes displayed a range of mean C t values, covering all the experimental conditions (fig. 1). The raw C t values ranged from 18.421 (TUB-α) to 27.679 (GST). The eight candidate reference genes exhibited variability in their C t values. The C t values of TUB-α, RPL32, ACT and GAPDH ranged from 18 to 22. The Ct values of TUB-β, TATA-box and SDHA ranged from 23 and 25. The least abundant transcript was GST with a mean C t value of 27.679.
Expression levels of candidate reference genes in different samples of G. daurica. The expression level of candidate G. dalerica reference genes in the test samples is shown in terms of the cycle threshold number (Ct-value). The black dot indicates the mean of duplicate samples, and the bars indicate the SD of the mean.
Expression stability analysis of the candidate reference genes
Among four programs, the gene stability rankings by BestKeeper analysis differed from the results generated by the other three methods. SDHA, SDHA and ACT and TUB-α were calculated as the most stable genes using the ΔC t method, geNorm and NormFinder, respectively (table 2). RefFinder, produced stability rankings (most stable to the least stable) genes in the temperature-stressed samples as: SDHA > TUB-α > ACT > GAPDH > RPL32 > TATA-box > TUB-β > GST (fig. 2a).
Expression stability of the candidate reference genes as calculated by RefFinder. A lower ranking value indicates more stable expression. Expression stability of reference genes in the following samples: (a) Different temperature of G dalerica; (b) different developmental stages of G. dalerica; (c) different sexes of G. dalerica; (d) different tissues of G. dalerica; (e) diapause condition.
Stability of reference gene expression under biotic and abiotic conditions.
For the development stage group, BestKeeper and geNorm analysis identified TUB-α, RPL32 and GST as the most stable genes. The stability rankings generated by the ΔC t method and NormFinder indicated that SDHA was the most stable gene(table 2). The order of gene expression stability (most stable to least stable) was calculated by RefFinder as: RPL32 > SDHA > GST > ACT > TUB-α > TATA-box > GAPDH > TUB-β (fig. 2b).
RefFinder showed that the order of gene expression stability for both of male and female adult group (most stable to least stable) was: ACT > TUB-α > SDHA > RPL32 > GST > TUB-β > TATA-box > GAPDH (fig. 2c). However, the most stable genes were completely different under different methods.
In different tissues of G. dalerica, the most stable genes were identified as SDHA, TUB-α, RPL32 and SDHA by geNorm, NormFinder, BestKeeper and ΔC t methods, respectively (table 2), and the eight candidate reference genes were ranked (highest to lowest stability) by the RefFinder as: SDHA > TUB-α > RPL32 > GAPDH > GST > TATA-box > ACT > TUB-β (fig. 2d).
The expression stability of the eight genes for diapause and non-diapause G. dalerica was ranked as SDHA > TUB-α > GAPDH > RPL32 > GST > ACT > TUB-β > TATA-box by RefFinder from the highest to the lowest (fig. 2e).
Optimum number of genes for normalization
The geNorm algorithm calculates an expression stability value (M) for each gene and then performs a pairwise comparison (Vn/Vn + 1) of this gene with others. A threshold of V value (Vn/Vn + 1) less than 1.5 was suggested for valid normalization. For the development stage group, the V2/3 and V3/4 variation values exceeded the proposed 0.15 cutoff, indicating that normalization with four stable reference genes was required, and the optimum gene combination number in other four groups was two (fig. 3). The most stable reference gene set according to the V values calculated with the geNorm is listed in table 3.
Determination of the optimal number of reference genes calculated by geNorm for normalization in G. daurica. The pairwise variation (Vn/Vn + 1) was analyzed by the geNorm software to determine the optimal number of reference genes included in the qRT-PCR analysis. A value less than 0.15 indicates that an additional reference gene will not significantly improve normalization.
Selected reference genes under different conditions.
Validation of selected reference genes in G. daurica
To demonstrate the effect of reference genes on target gene expression data, the relative expression levels of two target genes, CYP6 and Hsp70, were analyzed under all experimental conditions. For G. daurica exposed to different temperatures, CYP6 and Hsp70 in −14 and 10°C groups significantly overexpressed, compared to those in the 30°C group. Significant differences were found among the expression of CYP6 using the best reference gene (NF1: SDHA), the recommended normalization factors (NF1–2: SDHA, TUB-α), and the least stable gene (NF8: GST) (fig. 4a), but no significant differences were found in Hsp70 expression using the three sets of reference genes (fig. 4b). CYP6 expression levels exhibited significant differences among three larval stages, pupae and adults (P < 0.01) after normalization with the most stable gene (NF1: RPL32), the reference gene combination (NF1–3: RPL32, SDHA, GST), and the least stable gene (NF8: TUB-β) (fig. 4c). Similarly, expression levels of Hsp70 also exhibited differences except in the pupal stage (fig. 4d). When normalized using the best reference gene (NF1: ACT), the recommended normalization factors (NF1–2: ACT, TUB-α) and the least stable gene (NF8: GAPDH), CYP6 and Hsp70 expression in male adults were higher than in female adults, but no significant difference appeared among the three normalization factors (fig. 4e, f). CYP6 and Hsp70 expression levels were better normalized using the recommended normalization factors (NF1–2: SDHA, TUB-α) than normalized using the calculated best reference gene (NF1: SDHA) but significant differences were not evident in head, thorax and abdomen, whereas the expression levels were significantly less normalized using the unstable normalization factor (NF8: TUB-β), (fig. 4g, h). Normalized by the best reference gene (NF1: SDHA) and reference gene combination (NF1–2: ACT, TUB-α), the CYP6 expression in non-diapause adults increased by 1.61- and 1.56-fold compared with diapause adults, respectively. Reverse results were obtained when normalized against the gene with unstable normalization factor (NF8: TATA-box), and the CYP6 expression decreased to 0.59-fold (fig. 4i). Hsp70 in non-diapause adults was highly expressed compared with diapause adults after normalization using NF1, NF1–2 and NF8, but the expression levels using NF1 and NF1–2 were markedly higher than unstable reference gene TATA-box.
Validation of the gene stability measures. Expression profiles of Hsp70 and CYP6 in second-stage larvae at three different temperatures (a, b), in five development stages (c, d), in female and male adults (e, f), in three tissues (g, h), in diapause and non-diapause adults (i, j), were investigated under different normalization factors. Results are represented as mean ± SE.
Discussion
The qRT-PCR technique is superior to other conventional methods (northern hybridization and semi-quantitative PCR), and it is an essential procedure for analysis of gene expression (Hoogewijs et al., Reference Hoogewijs, Houthoofd, Matthijssens, Vandesompele and Vanfleteren2008; Huis et al., Reference Huis, Hawkins and Neutelings2010). Proper selection of reference genes is important for successful qRT-PCR analysis (Bustin et al., Reference Bustin, Benes, Garson, Hellemans, Huggett, Kubista, Mueller, Nolan, Pfaffl, Shipley, Vandesompele and Wittwer2009). However, using only one single endogenous control for normalization leads to a conclusion highly dependent on this single control gene, and can lead to inaccurate data interpretation (Ferguson et al., Reference Ferguson, Nam, Hopkins and Morrison2010). No single universal reference is available for different species under diverse conditions (Teng et al., Reference Teng, Zhang, He, Yang and Li2012). Therefore, it is important to select and validate reliable reference gene(s) stably expressed in different experimental conditions to minimize qRT-PCR analysis errors (Liang et al., Reference Liang, Guo, Zhou and Gao2014).
Our results demonstrate that it is difficult to find a universally applicable reference gene covering all conditions because gene expression is typically highly variable under different conditions. Therefore, we suggest that the stability of reference gene expression must be validated for each experimental condition under investigation. Ribosomal protein L32 (RPL32), a commonly used reference gene in gene expression studies (Scharlaken et al., Reference Scharlaken, De Graaf, Goossens, Brunain, Peelman and Jacobs2008; Shen et al., Reference Shen, Li, Ye, Wang, Lu and Xing2010b ), did not show good expression stability under all test conditions and it was not recommended as a proper normalized factor for G. daurica. Another widely used housekeeping gene, Actin (ACT), is a major component of the protein scaffolding that supports the cell and determines its shape. ACT has been used as a reference gene in many insect species (De Boer et al., Reference De Boer, De Boer, Marien, Timmermans, Nota, Straalen, Ellers and Roelofs2009; Hiel et al., Reference Hiel, Wielendaele, Temmerman, Soest, Vuerinckx, Huybrechts, Broeck and Simonet2009). In our research, ACT in G. daurica was qualified as the most stable reference gene under development stage conditions, but not for other conditions. Our results are consistent with the those of Li et al. (Reference Li, Xie, Wang, Wu, Yang, Yang, Pan, Zhou, Bai, Xu, Zhou and Zhang2013) and Zhu et al. (Reference Zhu, Yuan, Shakeel, Zhang, Wang, Wang, Zhan, Kang and Li2014), who found that actin was unsuitable for normalizing qRT-PCR data on sweetpotato whitefly, beet armyworm due to large errors and expression instability under different conditions.
Another important conclusion from our study is that multiple internal reference genes are necessary for studying gene expression under different experimental conditions. This is especially valid in the case of many samples, because more complex sample sets will exhibit higher reference gene variability. We found that two reference genes were sufficient for normalizing expression values of target genes in most of the samples, but four reference genes are recommended by pairwise variation (Vn/Vn + 1) and calculated by the geNorm in different developmental stage samples (table 3). Similarly, five reference genes were needed for all of the developmental stages samples; Zhu et al. (Reference Zhu, Yuan, Shakeel, Zhang, Wang, Wang, Zhan, Kang and Li2014) found that larger sample sizes require a higher number of reference genes for accurate normalization. However, we used reference genes in the top third to normalize the expression levels of target genes, because the threshold value of V < 0.15 is not absolute. Zhu et al. (Reference Zhu, Yuan, Shakeel, Zhang, Wang, Wang, Zhan, Kang and Li2014) considered additional reference genes to be required when adding more samples to a study, because it is more difficult to reach a minimum value of Vn/n + 1 when more unstable factors are introduced. Fu et al. (Reference Fu, Xie, Zhang, Wang, Wu, Liu, Zhou, Zhou and Zhang2013) demonstrated that the stability of multi-gene normalizer maybe may decline after adding a fourth, relatively unstable, reference gene. He recommended a combination of the three best reference genes for tissue samples as adequate.
Some researchers have used reference gene sets to analyze expression validation of target genes. We identified and compared the different expression levels of two target genes (CYP6 and Hsp70) with the most stable reference gene (NF1), the most unstable reference gene (NF8) and the recommended reference gene combination (NF1-2, or NF1-3). The target gene expression levels, in many samples, were better normalized using a reference gene combination than by using a single reference gene. In CYP6 expression analysis in different tissues, diapause and non-diapause adults, and Hsp70 in different tissues, it is clear that expression levels of target genes were completely reversed when normalized against NF1, NF1–2 and NF8. Following calculation by geNorm, the gene expression stability value (M) of all eight reference genes was <1.5, which illustrated that all of the reference genes can be used for normalization. However, using the most unstable normalized factor (NF8) would produce the opposite results. From these data validation tests, it is clear that extreme care must be taken for the selection of internal reference genes before their application of qRT-PCR. The stability of reference genes must be determined on a case-by-case basis.
This is the first report on establishing a standardized qRT-PCR procedure guideline for an important grassland insect pest. This study provides a foundation for advanced transcriptome validation tests and RNAi-based functional studies of G. daurica.
Acknowledgements
This research was supported by the National Natural Science Foundation of China (Grant No. 31360441). We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.
