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

1The 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. 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.

2Understanding 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).

3However, increasing evidence has suggested that HKGs varies among different insect species, in different tissue samples and under different experimental conditions. In coleopteran species Agrilus planipennis, for example, translation elongation factor 1α (EF1α) was the most stable gene, whereas glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and actin (ACT) showed least stability among six candidate HKGs for all the samples. In another coleopteran species Tribolium castaneum, the most stable ones were ribosomal protein genes, RPS3, RPS18, and RPL13a, whereas β-actin, α-tubulin, and RPS6 were not stable. Moreover, choosing unsuitable endogenous control genes resulted in low precision or misleading results. Therefore, HKGs should be validated to assure expression stability before using them as endogenous control genes in qRT-PCR.

4In 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. 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.

5Post-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 1st-, 2nd-, 3rd-, and 4th-instar stages of 2.5, 2.5, 3.0 and 5.0 days, respectively. Upon reaching full size, the 4th 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.

6L. decemlineata samples of developing eggs, 1st-, 2nd-, 3rd- and 4th-instar larvae, pupae (5 days after burrowing to soil), sexually mature adults (10 days after emergence) were collected from laboratory-rearing individuals. Moreover, the 4th-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 and 0.12 mg/L spinosad 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).

7To 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 non-redundant 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.

8The 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 MgCl2, 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® PCR Preps DNA Purification System (Promega). Purified DNA was ligated into the pGEM®-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).

9The 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.

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

11According 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.

12The 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.

13The 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).

14Putting 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.

15The 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.

16For developmental stages, RP18 and ARF1 were the best reference genes, with M-values below 0.5. RP4, ARF4, EF1α, TBP1 and TBP2 followed, that showed M-values between 0.5 and 1.0. ACT1 and ACT2 were the least stably expressed genes, with the M-values more than 1.0 (Figure 2A). The V2/3 and V3/4 values (the pairwise variation when the number of normalization factors is increased from two to three and from three to four) were above 0.15, whereas the V4/5 value was below 0.15 (Figure 2D). For different tissues, the stabilities of selected HKGs were RP18 > ARF1 > RP4 > EF1α > ARF4 > TBP1 > TBP2 > ACT1 > ACT2 (Figure 2B). The V2/3 and V3/4 values were above 0.15, whereas the V4/5, V5/6, V6/7, V7/8 and V8/9 values were below 0.15 (Figure 2E). For different treatments, the stabilities of selected HKGs were ARF1 > RP18 > RP4 > EF1α > ARF4 > TBP2 > TBP1 > ACT2 > ACT1 (Figure 2C). The V2/3, V3/4, V5/6, V6/7, V7/8 and V8/9 values were below 0.15, whereas the V4/5 value was above 0.15 (Figure 2F).

17In 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 4th-instar stage.

18NormFinder 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.

19For 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.

20According 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).

21For 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).

22In summary, BestKeeper analysis validates the findings of the geNorm and NormFinder algorithms.

23Detection 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. 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. Unfortunately, almost all HKGs are regulated by other ‘regulators’, no gene is constitutively expressed in all cell types and under all experimental conditions. Therefore, HKGs should be validated before using them as reference genes.

24Evaluation of the expression stability of HKGs requires mathematical methods. Many algorithms such as geNorm, NormFinder, BestKeeper, ΔCt approach, and stability index 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. geNorm and NormFinder use different mathematical methods to estimate the expression stability. The results of two algorithms can be used for cross validation. Thus, we used geNorm, NormFinder and BestKeeper to identify suitable reference genes for gene expression studies in the present paper.

25Ribosomal 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. In L. decemlineata, RP4 was the most stable among three HKGs, and was used as control gene. In other coleopterans, RPS3, RPS18, and RPL13a were suitable reference genes in T. castaneum, and RPL7 was the most stable in A. planipennis. In lepidopterans, RP49 had the greatest stability in different tissues of Chilo suppressalis, Bombyx mori and Spodoptera exigua. In dipteran Drosophila melanogaster, L32 was the most stable across three different treatments (i.e. injury, heat-shock stress and diets). In hemipterans, L26, RPS18, RPL18 and RPS9 were suitable reference genes in Rhodnius prolixus, Delphacodes kuscheli, Cimex lectularius and Aphis glycines, respectively.

26ADP-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. Previous results revealed that ARF was one of the most stable HKGs in plants such as wheat, barley, rye and citrus. 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.

27EF1α 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. 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, in dipteran D. melanogaster, in hemipteran A. glycines, in orthopteran Schistocerca gregaria and Chortoicetes terminifera, and in hymenopteran Bombus terrestris and Bombus lucorum. Meanwhile, E2F was the most stable gene in different tissues of lepidopteran B. mori, S. exigua and C. suppressalis.

28The 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.

29Actins 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, orthopteran S. gregaria, hemipteran R. prolixus and D. kuscheli, siphonapteran Lepeophtheirus salmonis, dipteran D. melanogaster and hymenopteran Apis mellifera. However, actin genes have been challenged for their suitability as the internal controls. 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. 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.

30In 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.

31PCR: Polymerase chain reaction; RT-PCR: Reverse transcriptase PCR; qRT-PCR: Quantitative real-time PCR; cDNA: Complementary DNA; HKG: House keeping gene; ACT: Actin; ARF: ADP-ribosylation factor; TBP: TATA box binding protein; RP: Ribosomal protein; EF: Translation elongation factor; SD: Standard deviation.

32The authors declare that they have no competing interests.

33XQS 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.