Validation of reference genes for expression analysis by quantitative real-time PCR in Leptinotarsa decemlineata (Say)

Background L. decemlineata is an exotic invasive insect pest, and invaded in Xinjiang Uygur autonomous region in China in the 1990s from Kazakhstan. It is a notorious defoliator of potato throughout most of the northern Xinjiang in current, and often causes extremely large yield losses of potato. Results The expression stability of nine L. decemlineata house-keeping genes (Actin, ACT1 and ACT2; ADP-ribosylation factor, ARF1 and ARF4; TATA box binding protein, TBP1 and TBP2; ribosomal protein RP4 and RP18; translation elongation factor 1α EF1α) was evaluated by quantitative real-time polymerase chain reaction (qRT-PCR) in seven developmental stages, three larval tissues and two insecticide treatments. The results were analyzed using three software programs: geNorm, NormFinder and BestKeeper. Although there was no consistent ranking observed among the house-keeping genes across the samples, the overall analysis revealed that RP18, RP4, ARF1, and ARF4 were the four most stable house-keeping genes. In contrast, ACT1 and ACT2, two of the most widely used reference genes, had the least stability. Our results suggest that the combined use of the four most stably expressed genes may produce optimal normalization for qRT-PCR. Conclusions The expression stability of the house-keeping genes varies among different developing stages, in different tissues and under different experimental conditions. Our results will enable a more accurate and reliable normalization of qRT-PCR data in L. decemlineata.


Background
The Colorado potato beetle, Leptinotarsa decemlineata (Say), is an exotic invasive insect pest. It invaded in Xinjiang Uygur autonomous region in China in the 1990s from Kazakhstan. It is a notorious defoliator of potato throughout most of the northern Xinjiang in current, and often causes extremely large yield losses of potato [1][2][3][4][5]. L. decemlineata has a complicated and diverse life cycle. Moreover, the beetle has developed resistance to many classes of insecticides, among them are organophosphates, carbamates, pyrethroids and neonicotinoids [6,7].
Understanding the molecular regulation mechanisms that underlie these ecological and physiological adaptations in L. decemlineata may provide insights into complex regulatory networks, and may help to develop intriguing targets for the control of this pest. Using the next generation sequencing method, we have obtained a transcriptomic database. Deciphering gene expression profiles and validation of mRNA levels for target genes via quantitative real-time polymerase chain reaction (qRT-PCR) have been crucial to on-going studies. For normalizing qRT-PCR results, the mRNA levels of the target genes need to normalize to internal control house-keeping genes (HKGs) [8][9][10].
In L. decemlineata, the expression stability of three HKGs, β-actin, RP4 and RP18, was evaluated. RP4 was the most stable among samples and was used as reference gene [29]. According to Vandesompele et al. (2002), however, at least two endogenous control genes are recommended. In the present paper, therefore, we selected 9 HKGs (Actin, ACT1 and ACT2; ADP-ribosylation factor, ARF1 and ARF4; TATA box binding protein, TBP1 and TBP2; ribosomal protein RP4 and RP18; translation elongation factor 1α EF1α) from the transcriptome and validated the stability of their expression. We have identified several suitable reference genes for gene expression studies.

Insect collection and rearing
Post-diapause L. decemlineata adults were collected from potato field in spring at Urumqi city (43.82 N, 87.61E), Xinjiang Uygur autonomous region in China. Insects were routinely reared in an insectary at 28 ± 1°C under a 14 h:10 h light-dark photoperiod and 50-60% relative humidity using fresh potato foliage as food. The adults deposited their eggs in batches of about 30 on the underside of potato leaves. After approximately 7 days, the eggs hatched into reddish-brown larvae. Larvae progressed through four distinct instars, with the average periods of the 1 st -, 2 nd -, 3 rd -, and 4 th -instar stages of 2.5, 2.5, 3.0 and 5.0 days, respectively. Upon reaching full size, the 4 th instars spent an additional 4-7 days as a non-feeding prepupae. The prepupae then dropped to the soil and burrowed to a depth of 3-5 cm to pupate. The pupae emerged in roughly 10 days. Both the male and female adults spent an average of 7 days to become sexually mature.
Total RNA extraction and cDNA synthesis L. decemlineata samples of developing eggs, 1 st -, 2 nd -, 3 rd -and 4 th -instar larvae, pupae (5 days after burrowing to soil), sexually mature adults (10 days after emergence) were collected from laboratory-rearing individuals. Moreover, the 4 th -instar larvae were dissected to obtain midgut, fat body and cuticle. Furthermore, the 4th-instar larvae feeding on foliage immersed 0.09 mg/L chlorantraniliprole [1] and 0.12 mg/L spinosad [2] represented insecticide treatments. The samples were frozen with liquid nitrogen and homogenized in a tissue grinder. Then, total RNAs were extracted with TRIzol reagent (Invitrogen Co., USA) according to the manufacturer's specifications, and were treated for 30 min at 37°C with RNase free DNase I (Ambion, Austin, TX) to eliminate traces of chromosomal DNA. The integrity of RNA was checked on a 1.5% agarose gel and visualized by ethidium bromide staining. The purity and amount of the total RNA samples were determined by NanoDrop ND-1000 spectrophotometer (Nanodrop Technologies, Rockland, DE, USA). First strand cDNA was synthesized from 1 μg total RNA using M-MLV reverse transcriptase (Takara Bio, Dalian, China) and Oligo (dT18) as the anchor primer. The reaction mixtures were incubated at 70°C for 10 min followed by 42°C for 1 h and 70°C for 15 min. The cDNA samples were used as templates for polymerase chain reaction (PCR).

Selection and authentication of candidate HKGs
To search for HKG sequences from L. decemlineata transcriptome data, a reciprocal BLAST hits approach was used. The HKGs from other insect species in GenBank were downloaded from NCBI (http://www. ncbi.nlm.nih.gov/), queried individually to L. decemlineata transcriptome using the TBLASTN program with a permissive E-value cutoff of 10 -3 to get the hits. And then, each of the queried hits was compared back against nonredundant database of NCBI by the BLASTX program (E-value <10 -3 ) to determine whether the original sequence was one of the hits. The selected HKG sequences were listed in Table 1.
The unigenes of selected nine HKGs were assembled and clustered from short reads, there were inevitably issues with clone contamination and mix-up. Reverse transcriptase PCR (RT-PCR) was performed to authenticate the HKGs using the primers listed in Table 1. The components of PCR reaction buffer were 2.5 mM of dNTP, 10 mM of each primer, 25 mM of MgCl 2 , 5 U/μL of Ex-Taq DNA polymerase (Takara Bio, Dalian, China), in a total volume of 25 μL. Thermal cycling conditions of RT-PCR were available from the authors upon request. The amplified products were separated by electrophoresis on 1.5% agarose gel and purified using the Wizard W PCR Preps DNA Purification System (Promega). Purified DNA was ligated into the pGEM W -T easy vector (Promega) and several independent subclones were sequenced from both directions. The nucleotide sequences obtained after the sequence analysis were submitted to GenBank database (Accession No. KC190026-KC190034).

Quantitative real-time PCR (qRT-PCR)
The qRT-PCR primers were designed using Beacon Designer 7 (Premier Biosoft International, Palo Alto, Calif., USA), and were given in Table 2. The qRT-PCR reactions were performed using SYBR Premix Ex Taq (Perfect Real Time) (Takara Co., Otsu, Japan) and ABI Real-Time 7300 PCR system (Applied Biosystems) according to the manufacturer's protocol. The reaction mixture consisted of 2 μL of cDNA template (corresponding to 0.9 ng of the starting amount of RNA), 10 μL of SYBR Premix Ex Taq (Takara), 1 μL of forward primer (10 μM), 1 μL of reverse primer (10 μM), 0.4 μL of Rox Reference Dye (50×) in a final reaction volume of 20 μL. A reverse transcription negative control (without reverse transcriptase) and a non-template negative control were included for each primer set to confirm the absence of genomic DNA and to check for primer-dimer or contamination in the reactions, respectively. The qRT-PCR protocol included an initial step of 95°C for 30 sec, followed by 40 cycles of 95°C for 5 sec and then annealed at 60°C for 31 sec, followed by one cycle of 95°C for 15 sec, 60°C for 60 sec, and 95°C for 15 sec. PCR amplicons were subjected to melting curve analysis. The specificity of the qRT-PCR reactions was monitored with melting curve, analyzing by SDS software (version 1.4) and gel electrophoresis. Amplification efficiencies were determined by a 10-fold dilution series of template. All experiments were repeated in triplicate.

Data processing
The raw Ct values were obtained using the SDS software of ABI 7300 (version 1.4). The algorithms including geNorm [8], BestKeeper [30] and NormFinder [9] were used to analyze the stability of selected HKGs, strictly following the manuals of the algorithms.

Selection of candidate housekeeping genes
According to the published results, nine HKGs were selected. They were given names as ACT1, ACT2, ARF1, ARF4, EF-1α, RP4, RP18, TBP1 and TBP2, corresponding to homologous genes in the T. castaneum genome ( Table 1). The sequence correctness of the nine HKGs was proven by RT-PCR and the obtained sequences were submitted to GenBank database.
The products from qRT-PCR were also confirmed by bi-direction sequencing. Primer specificities for qRT-PCR were verified by melting curve analysis. All primer pairs amplified a single PCR product with the expected sizes, showed a slope less than −3.0, and exhibited efficiency values ranging from 2.0-2.1 by a 10-fold dilution series of template (Table 2). These data indicate that amplification efficiencies of primers reach the standard requirements of conventional qRT-PCR.

Expression levels of the candidate reference genes
The temporal expression patterns of mRNAs encoding the nine candidate HKGs were analyzed by qRT-PCR. All of the nine genes were expressed in seven developmental stages, and among three larval tissues, indicated by the presence of a single amplicon of the expected size on an agarose gel (data not shown).
Putting all the data with the same HKG at all tested samples together, the raw expression levels varied dramatically. According to the variations of Ct values, RP18 showed the smallest gene expression variation (below ten cycles), followed by ARF1, ARF4 and RP4, whereas ACT2 and ACT1 had the highest expression variation (Figure 1). Thus, determining HKGs as reference genes requires careful confirmation of expression stability.

geNorm analysis
The geNorm defines two parameters to quantify the gene stability: M (the average expression stability) and V (the pairwise variation). The gene with the lowest M value is considered to have the most stable expression, while the one with the highest M value has the least stable expression. The V value should be below the default cut-off value of 0.15.
In summary, geNorm analysis indicates that the combined use of the four most stably expressed genes (ARF1, RP18, RP4 and ARF4) produces optimal normalization for qRT-PCR among different developing stages or tissues. Moreover, combining two most stable genes ARF1 and RP18 is sufficient to normalize the target genes among different treatments within the 4 th -instar stage.

NormFinder analysis
NormFinder is designed to calculate stability using the combined estimation of intra-and inter-group expression variations of the analyzed genes. According to the NormFinder, the genes that are more stably expressed are indicated by lower average expression stability values. For developmental stages, the most stable genes were ARF1 and RP18, the most unstable genes were ACT1 and ACT2 ( Figure 3A). For different tissues, the stabilities of selected HKGs were RP4 > RP18 > ARF1 > EF1α > ARF4 > TBP1 > TBP2 > ACT2 > ACT1 ( Figure 3B). Among the selected HKGs for different treatment, the four most stable genes were RP4, RP18, EF1α and ARF1. And again, the least stable genes were ACT2 and ACT1 ( Figure 3C). Thus, the NormFinder results validate the findings of the geNorm algorithm, and ARF1, RP18, RP4, EF1α and ARF4 are among the best combination of HKGs.

BestKeeper analysis
According to the BestKeeper analysis, the stabilities of selected HKGs were RP18 > ARF1 > RP4 > ARF4 > EF1α > TBP1 > TBP2 > ACT1 > ACT2 in developmental stages ( Figure 3D). The weighted index BestKeeper calculated for the nine candidates showed an SD of CP = ±0.98 cycles. The SD (±CP) value was higher for EF1α, TBP1, TBP2, ACT1 and ACT2. This constituted a reason to exclude these genes from the BestKeeper index calculation, because they were not reliable reference genes in this setting. After the exclusion of EF1α, TBP1, TBP2, ACT1 and ACT2 from the index, its variation decreased (SD = ±0.77 cycles). After further exclusion of ARF4 and RP4 subsequently, the variation SD was ±0.71 and ±0.65 cycles, respectively (Table 3).
For different tissues, RP18, ARF1 and RP4 were the most stable because they showed SD values lower than 1. EF1α, ARF4, TBP1 and TBP2 had SD values around 1.2. And ACT1 and ACT2 were the most unstable genes ( Figure 3E). Among the selected HKGs for different treatment, the four most stable genes were RP4, RP18, EF1α and ARF1 due to their SD values lower than 1. And again, the least stable genes were ACT2 and ACT1, with their SD values more than 1.2 ( Figure 3F).
In summary, BestKeeper analysis validates the findings of the geNorm and NormFinder algorithms.

Discussion
Detection and quantification of transcript abundance in different conditions are important tasks in molecular biology. qRT-PCR is an exceptional and trustworthy technique even for low abundant mRNA transcripts [31]. However, the variations of qRT-PCR will be unavoidably introduced during RNA preparation, cDNA synthesis, and PCR process. The use of HKGs as reference genes to normalize gene expression is a strategy to minimize the variations of qRT-PCR. A good reference gene should meet three criteria. Firstly, its amplification efficiency is similar to the target genes. Secondly, it is expressed at moderate level. Lastly, its expression is stable in all test samples [11]. Unfortunately, almost all HKGs are regulated by other 'regulators' , no gene is constitutively expressed in all cell types and under all experimental conditions [32]. Therefore, HKGs should be validated before using them as reference genes.
Evaluation of the expression stability of HKGs requires mathematical methods. Many algorithms such as geNorm [8], NormFinder [9], BestKeeper [30], ΔCt approach [33], and stability index [34] have been developed. Among them geNorm, BestKeeper and NormFinder are the most common statistical algorithms. Moreover, the three algorithms are easy to use and freely available for download [28,30,[35][36][37][38]. geNorm and NormFinder use different  Ranking, stability and determination of optimal number of reference genes using geNorm. A, B and C represent ranking based on average expression stability value ( M) for developmental stages, tissues and treatments. In the plots, genes were ordered from least (left) to most (right) stable. D, E and F display optimal number of genes required for the accurate estimation of the target gene mRNA in developmental stages, tissues and treatments calculated by the pairwise variation (V) of the normalization factors (NFn and NFn + 1). A large pairwise variation indicates that the added gene has a significant effect and should preferably be included in the normalization. The cut-off value for such significance should be 0.15. mathematical methods to estimate the expression stability. The results of two algorithms can be used for cross validation [8,9]. Thus, we used geNorm, NormFinder and BestKeeper to identify suitable reference genes for gene expression studies in the present paper.
Ribosomal proteins are involved in translation and protein synthesis. In this study, we found that ribosomal protein RP4 and RP18 were among the most stable reference genes in L. decemlineata. Consistent with our results, ribosomal proteins are reported to be the best reference genes in many insects [22,26]. In L. decemlineata, RP4 was the most stable among three HKGs, and was used as control gene [29]. In other coleopterans, RPS3, RPS18, and RPL13a were suitable reference genes in T. castaneum [20], and RPL7 was the most stable in A. planipennis [12]. In lepidopterans, RP49 had the greatest stability in different tissues of Chilo suppressalis, Bombyx mori and Spodoptera exigua [11]. In dipteran Drosophila melanogaster, L32 was the most stable across three different treatments (i.e. injury, heat-shock stress and diets) [15]. In hemipterans, L26, RPS18, RPL18 and RPS9 were suitable reference genes in Rhodnius prolixus [13], Delphacodes kuscheli [16], Cimex lectularius [17] and Aphis glycines [14], respectively.
ADP-ribosylation factors (ARFs) are ubiquitous in eukaryotic cells. ARFs are soluble, or associate with membranes due to their N-terminus myristoylation. They function as regulators of vesicular traffic and actin remodeling [39]. Previous results revealed that ARF was one of the most stable HKGs in plants such as wheat, barley, rye and citrus [40][41][42]. In the present paper, we provided the first line of evidence that ARF1 and ARF4 are also stably expressed in animal species such as L. decemlineata.
EF1α encodes a protein that is involved in protein synthesis. It is widely used as reference gene in gene expression studies in plants and animals, including insects [12,21,23]. Our data showed that EF1α exhibited moderate stability across the samples assayed in L. decemlineata. In contrast, EF1α was the most appropriate reference gene in another coleopteran A. planipennis [12], in dipteran D. melanogaster [15], in hemipteran A. glycines [14], in orthopteran Schistocerca gregaria [23] and Chortoicetes terminifera [18], and in hymenopteran Bombus terrestris and Bombus lucorum [22]. Meanwhile, E2F was the most stable gene in different tissues of lepidopteran B. mori, S. exigua and C. suppressalis [11].
The TATA-binding protein (TBP) is a general transcription factor that binds specifically to a DNA sequence called the TATA box. Our data showed that TBP1 and TBP2 were among the least stable across the samples in L. decemlineata. In contrast, TBP showed high stability across the samples assayed in A. glycines [14].
Actins play key roles in cell motility and cytoskeleton maintenance. Therefore, they are assumed to be constitutively expressed and widely used as reference genes. Actin was among the most stable in collembolan Folsomia candida and Orchesella cincta [24], orthopteran S. gregaria [23], hemipteran R. prolixus [13] and D. kuscheli [16], siphonapteran Lepeophtheirus salmonis [43], dipteran D. melanogaster [15] and hymenopteran Apis mellifera [26]. However, actin genes have been challenged for their suitability as the internal controls [44]. In lepidopteran insects, actin A1 was the most stable in Plutella xylostella and C. suppressalis, but it was the least stable in B. mori and S. exigua across different developmental stages [11]. In the present paper, we found ACT1 and ACT2 were the most unstable genes, irrespective of the software programs used. Similarly, ACT showed least stability among the candidate reference genes analyzed in A. planipennis [12].

Conclusion
In summary, our work presented a number of stable HKGs that are suitable to be used as the reference genes in L. decemlineata. This will enable a more accurate and reliable normalization of qRT-PCR data. Moreover, these validated reference genes could also serve as the basis for selection of candidate internal controls in any given insect species.

Competing interests
The authors declare that they have no competing interests.
Authors' contributions XQS and WCG performed all the experimental procedures, data analysis, and was the primary author of the manuscript. PJW, LTZ and XLR assisted in manuscript revising and provided helpful discussions. TA and KYF maintained the insects. GQL wrote the manuscript, conceived and supervised the research. All authors read and approved the final manuscript.  In the three last columns the BestKeeper index is computed together with the same descriptive parameters, for four genes (ARF4, ARF1, RP4 and RP18) or for three genes after removal of ARF4 (ARF1, RP4 and RP18) or for two genes after removal of RP4 (ARF1 and RP18).