Selection of appropriate reference genes for RT-qPCR analysis in Propylea japonica (Coleoptera: Coccinellidae).

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is a reliable technique commonly used in molecular biology to analyze RNA expression. The selection of suitable reference genes for data normalization is a precondition for credible measurements of gene expression levels using RT-qPCR. Propylea japonica is one of the most common pests of many crop systems throughout East Asia, and has often been used in the testing of non-target impacts during environmental risk assessments of genetically engineered plants. The present study assessed the suitability of nine frequently used reference genes for comparisons of P. japonica gene expression. Expression stability was compared across developmental stages, sex, a range of tissues, and following exposure to different temperatures. Data were analyzed using RefFinder, which integrated the results obtained using NormFinder, geNorm, BestKeeper, and the ΔCt method. This led to the identification of unique sets of reference genes for each experimental condition: ribosomal protein S18 (RPS18) and elongation factor 1 α (EF1A) for developmental stage comparisons, RPS18 and EF1A for sex comparisons, EF1A and ribosomal protein L4 for tissue comparisons, and RPS18 and EF1A for analyses of temperature-mediated effects. These reference genes will help to enhance the accuracy of RT-qPCR analyses of P. japonica gene expression. This work represents an initial move towards building a standardized system for RT-qPCR analysis of P. japonica, providing a basis for the ecological risk assessment of RNAi-based insect control products.

Introduction

Reverse transcriptase-quantitative PCR (RT-qPCR) is frequently employed as a powerful method for the quantification of gene expression. However, various factors, including RNA quantity and quality, cDNA quantity and quality, and PCR efficiency can significantly influence the quantification cycle (C) values obtained using this method [1,2]. RT-qPCR data are generally analyzed by normalizing target gene expression to that of one or more suitable reference genes. Although reference gene expression profiles should ideally be stable under the relevant experimental conditions, previous studies have indicated that the expression of many frequently used reference genes can be markedly affected by different treatments [3-11]. Therefore, preliminary evaluations should be conducted to identify stable reference genes for RT-qPCR analysis in a given species under the proposed experimental conditions.

Reference gene selection for RT-qPCR normalization in insect gene expression studies over the past 10 years was recently reviewed for studies that employed the most widely used SYBR Green method [12]. A total of 39 experimental factors were investigated in these papers [12]. Many of these studies employed RNA interference (RNAi), which has been widely used to investigate insect gene function under a range of experimental conditions. RNAi is a biological process in which RNA molecules inhibit protein production by neutralizing targeted mRNA molecules [13]. The mode of action of RNAi-based insecticides and/or RNAi genetically engineered plants suggests that unintended off-target impacts may occur due to altered gene expression in non-target organisms [13,14], and RT-qPCR offers an convenient instrument to detect these gene expression changes. The lady beetle, Propylaea japonica (Thunberg) (Coleoptera: Coccinellidae) is one of the most widespread indigenous natural pests in many planting systems in East Asia. As a representative species of predator lady beetles, P. japonica has been widely used to assess the latent risks of Bacillus thuringiensis crops [15-20]. This species is also likely to be selected as a surrogate species to evaluate the effects of RNAi-based insect control products. Because these products may cause lethal or sub-lethal effects on P. japonica by altering gene expression, it is important to identify suitable reference genes for RT-qPCR analyses of this species.

In this study, we aimed to identify stable reference genes for RT-qPCR analysis in P. japonica. Nine frequently-used reference genes were investigated: β-actin (Actin), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), α-tubulin (TUBA), elongation factor 1 α (EF1A), ribosomal protein L4 (RPL4), ribosomal protein S18 (RPS18), arginine kinase (ArgK), heat shock protein 90 (HSP90), and vacuolar-type H-ATPase subunit A (V-ATPase A). All of these reference genes have been used frequently for RT-qPCR analyses in other insects (S1 Table). The stability of each candidate was assessed for four experiments evaluating the effects of developmental stage, sex, tissue, and temperature on gene expression.

Materials and methods

Insects

Propylea japonica adults were collected in the Mengshan Mountain region (Shandong Province, China) during June 2010 [15]. Since then, the colony has been maintained in the laboratory at a temperature of 26 ± 1°C and a relative humidity of 60–80%, with a 14: 10 h light: dark cycle. They were supplied with Aphis craccivora Koch, which were reared on fava bean (Vicia faba) in a greenhouse at 20–28°C.

Experimental conditions

Each P. japonica developmental stage was sampled on the first day of each stage; this included eggs, four larval instars, pupae, and female and male adults. The numbers of sampled individuals for each replicate in each stage was as follows: 15 eggs; five individuals for the 1st instar; five individuals for the 2nd instar; three individuals for the 3rd instar; one individual for the 4th instar; one pupa; and one female or male individual for adult female or male stages. Different body tissues, including head, midgut, Malpighian tubule, and carcass (body except for the above tissues) were dissected from the 4th instar larvae and female and male adults; about 15 individuals were dissected per replicate. The tissues were stored in RNAlater (Thermo Fisher Scientific Inc., Waltham, MA, USA) at 4°C until total RNA isolation.

To investigate temperature-mediated effects, three replicate samples of five 1st instars were maintained at 8, 25, or 35°C for 3 h. The samples were then placed in 1.5-ml centrifuge tubes, snap-frozen in liquid nitrogen, and stored at -80°C until total RNA isolation.

Total RNA extraction and cDNA synthesis

Total RNA samples were extracted from eggs and Malpighian tubules using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), as previously described [20]. Total RNAs were isolated from the other samples using the HiPure Total RNA Micro Kit (Magen, Shanghai, China), in accordance with the manufacturer’s instructions. Gel electrophoresis and the NanoDrop One spectrophotometer (Thermo Fisher Scientific) were used to determine the quantity of total RNA. This was then dissolved in 10–70 μl ddH2O to obtain the following RNA concentrations (mean ± standard error of the mean): 200.7 ± 29.4 ng/μl for eggs; 284.4 ± 16.4 ng/μl for first instars; 637.8 ± 52.5 ng/μl for second instars; 563.0 ± 96.8 ng/μl for third instars; 864.3 ± 177.3 ng/μl for fourth instars; 804.9 ± 34.4 ng/μl for pupae; 831.4 ± 88.1 ng/μl for male adults; 866.6 ± 129.6 ng/μl for female adults; 279.8 ± 18.5 ng/μl for heads; 821.2 ± 140.4 ng/μl for carcasses; 359.5 ± 104.1 ng/μl for midguts; 553.2 ± 149.1 ng/μl for Malpighian tubules; 290.67 ± 28.71 ng/μl for first instars at 8°C; 240.57 ± 53.29 ng/μl for first instars at 25°C; and 277.10 ± 53.19 ng/μl for first instars at 35°C. The 260/280 nm optical density ratios were between 1.9 and 2.1 for all samples. The PrimeScript RT kit (containing gDNA Eraser, Perfect Real Time; TaKaRa, Dalian, China) was used to prepare first-strand cDNA for gene expression analysis. The cDNA was diluted tenfold prior to the following RT-qPCR investigations.

Gene cloning and primer design

A total of nine reference genes were assessed (S1 Table, Table 1). Degenerate primers for ArgK had been designed previously [10]. Primers for the other eight genes were designed using the sequences obtained from previous transcriptome datasets [14] (GenBank accession: SRX554957).

Table 1

Primers used for RT-qPCR.

Gene	Primer sequences (5’-3’)	Length(bp)	Efficiency(%)	R2	Linear regression
Actin	F: TGTGCTATGTCGCTTTGG	129	95.4	0.9988	y = -3.4376x+18.887
	R: CTGGGCAACGGAATCTTT
GAPDH	F: GTATCGGTCGTCTTGTACTG	121	93.5	0.9983	y = -3.4879x+22.93
	R: CCATGGGTGGAGTCATATTT
EF1A	F: CTGGAAAGACCACAGAAGAAA	114	91.7	0.9999	y = -3.5381x+18.601
	R: GAGGAGGGAATTCTTGGAAAG
TUBA	F: TGGTTGATAATGAAGCCATCTA	117	99.8	0.9813	y = -3.3267x+23.937
	R: GAGAAGCAGTGATTGAAGAAAC
RPL4	F: CGTCGTCTTAACCCACTTAC	118	91.9	0.9965	y = -3.532x+20.291
	R: CTTCTTCTCTGGCCAACTG
RPS18	F: CGCTGGTGATTCCAGATAAA	111	102.1	0.9992	y = -3.2836x+23.82
	R: GACGACCTACACCTTTGATG
HSP90	F: GTTACCAATCCCTCACCAATC	132	90.2	0.9994	y = -3.5803x+20.565
	R: CTAAATCGGCCTTGGTCATAC
ArgK	F: GACGTTCTTTGGAGGGATAC	102	91.8	0.9999	y = -3.5349x+21.624
	R: CATCGTCGAGTCCAGATAAAG
V-ATPase A	F: CATCTGCCACTCTTGGTATC	120	100.1	0.9987	y = -3.3187x+24.889
	R: CCAAAGCTCTCGTGTACTTC

PCR reactions were performed using a total volume of 20 μl, as described previously [21]. Amplicons of the expected lengths were purified using the TIANgel Midi Purification Kit (TIANGEN, Beijing, China), and subcloned into the pClone007 Blunt vector before transformation into Escherichia coli DH5α competent cells (TSINGKE, Beijing, China) for sequencing by TSINGKE company. Reference gene sequences were confirmed by comparison with the NCBI database.

RT-qPCR analysis

The RT-qPCR reactions were conducted in accordance with our previous study [21]. The melting curve and standard curve for each candidate gene was also generated as described previously [21]. The RT-qPCR efficiencies (E) were calculated using the following equation: E = (10[-1/slope] - 1) × 100.

Determination of reference gene expression stability

The stabilities of the nine reference genes were assessed using the following four approaches: geNorm [22], NormFinder [23], BestKeeper [24], and the ΔCt method [25]. Finally, the findings of these four analytical tools were integrated by RefFinder (http://150.216.56.64/referencegene.php), providing a stability ranking of the candidates. The optimal number of reference genes for target gene normalization was determined by pairwise variation (Vn/Vn+1) using V-values calculated by geNorm [22]. A Vn/Vn+1 cutoff value of ≤ 0.15 signified that the additional n + 1 reference gene was unnecessary; this indicated the appropriate number of reference genes for RT-qPCR data normalization;

Results

Candidate gene cloning and performance

All reference genes were expressed in P. japonica and each was visualized as a single amplicon (S1 Fig). The specific amplification of all reference genes was confirmed by melting curve analyses (Fig 1).

Fig 1

Melting curves of the nine reference genes examined for Propylea japonica.

β-actin (Actin), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1 α (EF1A), α-tubulin (TUBA), ribosomal protein L4 (RPL4), ribosomal protein S18 (RPS18), heat shock protein 90 (HSP90), arginine kinase (ArgK), and vacuolar-type H-ATPase subunit A (V-ATPase A).

Table 1 shows the E of each PCR, the linear regression equation, and the correlation coefficient (R2) for each standard curve. The standard curve for each gene is also shown (S2 Fig). The C values for these reference genes under the four experimental situations ranged from 19 to 26. EF1A and Actin had the highest expression levels, whereas TUBA and ArgK showed the lowest levels of expression (Fig 2).

Fig 2

Expression profiles of the nine Propylea japonica reference genes.

The Cq values for each gene are shown for the four experimental conditions. β-actin (Actin), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1 α (EF1A), α-tubulin (TUBA), ribosomal protein L4 (RPL4), ribosomal protein S18 (RPS18), heat shock protein 90 (HSP90), arginine kinase (ArgK), and vacuolar-type H-ATPase subunit A (V-ATPase A).

Reference gene expression stability for each experimental condition

Table 2 shows the overall order of gene expression stability determined using geNorm, NormFinder, BestKeeper, and the ΔC method, from the most consistent to the least consistent reference genes, under each experimental condition. The major results obtained using each method are indicated below.

Table 2

Stability of the expression of nine candidate reference genes in Propylea japonica under different experimental conditions calculated by the four different analytical tools geNorm, NormFinder, BestKeeper, and the ΔCt method, respectively.

Conditions	CRGs*	geNorm		Normfider		BestKeeper		ΔCt
		Stability	Rank	Stability	Rank	Stability	Rank	Stability	Rank
Developmental	EF1A	0.494	2	0.861	4	0.636	3	1.861	1
Stage	RPS18	0.366	1	0.641	1	0.582	1	1.871	2
	Actin	2.656	8	4.371	9	2.977	8	4.708	9
	RPL4	0.366	1	0.907	5	0.730	4	1.936	3
	ArgK	2.069	7	3.953	8	3.322	9	4.388	8
	HSP90	0.890	5	1.267	6	0.912	6	2.121	6
	GAPDH	0.808	4	0.834	3	0.790	5	2.038	5
	V-ATPase A	0.713	3	0.719	2	0.599	2	2.005	4
	TUBA	1.280	6	2.346	7	1.571	7	2.975	7
Sex	EF1A	0.179	1	0.090	1	0.588	5	0.881	1
	RPS18	0.260	2	0.123	2	0.487	3	0.887	2
	Actin	0.972	7	2.188	7	1.288	8	2.304	8
	RPL4	0.350	4	0.185	4	0.705	6	0.971	5
	ArgK	0.590	6	0.714	5	0.161	1	1.209	6
	HSP90	0.179	1	0.090	1	0.505	4	0.903	3
	GAPDH	0.308	3	0.170	3	0.452	2	0.920	4
	V-ATPase A	0.468	5	0.870	6	0.983	7	1.232	7
	TUBA	1.312	8	2.397	8	1.841	9	2.504	9
Tissue	EF1A	0.299	1	0.149	1	0.357	2	1.082	1
	RPS18	0.339	2	0.161	2	0.409	3	1.088	2
	Actin	1.206	7	1.493	7	1.077	7	1.850	7
	RPL4	0.299	1	0.239	3	0.273	1	1.148	3
	ArgK	1.557	8	2.637	9	2.199	9	2.786	9
	HSP90	0.457	3	0.350	4	0.683	5	1.168	4
	GAPDH	0.828	5	1.125	6	1.015	6	1.532	6
	V-ATPase A	1.022	6	1.616	8	1.293	8	1.880	8
	TUBA	0.670	4	1.024	5	0.580	4	1.477	5
Temperature	EF1A	0.082	1	0.031	2	0.149	2	0.326	1
	RPS18	0.082	1	0.124	4	0.137	1	0.343	3
	Actin	0.152	2	0.051	3	0.201	4	0.342	2
	RPL4	0.179	3	0.028	1	0.247	5	0.350	4
	ArgK	0.228	4	0.305	5	0.200	3	0.431	5
	HSP90	0.275	5	0.308	6	0.298	6	0.453	6
	GAPDH	0.316	6	0.412	7	0.323	7	0.506	7
	V-ATPase A	0.381	7	0.550	8	0.460	8	0.617	8
	TUBA	0.452	8	0.649	9	0.464	9	0.700	9

* Candidate reference gene

geNorm

Across different developmental stages, RPS18 and RPL4 were both ranked as the most stable genes, while EF1A and HSP90 were ranked together as the most stable genes in the sex comparison. For the tissue comparisons, EF1A and RPL4 were both ranked as the most stable genes, while EF1A and RPS18 were ranked together as the most stable genes in the temperature experiment.

NormFinder

Across different developmental stages, RPS18 was the most stable gene. In females and males, EF1A and HSP90 were ranked together as the most stable genes. Among different tissues, EF1A ranked as the most stable gene, while RPL4 showed the most stable expression in the temperature experiment.

BestKeeper

Across different developmental stages, RPS18 was the most stable gene, while ArgK showed the most stable expression in females and males. For the tissue comparisons, RPL4 was the most stable gene, whereas RPS18 showed the most stable expression at different temperatures.

The ΔCt method

The average standard deviation of each gene set was negative correlated with its stability. The most stable gene under each of the four experimental conditions was EF1A.

The overall RefFinder ranking of reference gene expression stability

The comprehensive reference gene rankings for expression stability under each experimental condition are shown in Fig 3.

Fig 3

Stability of the nine Propylea japonica reference genes according to RefFinder.

A lower Geomean value indicates a more stable expression. β-actin (Actin), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1 α (EF1A), α-tubulin (TUBA), ribosomal protein L4 (RPL4), ribosomal protein S18 (RPS18), heat shock protein 90 (HSP90), arginine kinase (ArgK), and vacuolar-type H-ATPase subunit A (V-ATPase A).

Optimal number of reference genes for target gene normalization based on geNorm

Although V-values for the analyses of developmental stage were never < 0.15, V2/3 was lower than V3/4 (Fig 4A). In addition, the expression of Actin and ArgK was obviously lower at the egg stage than at other stages. When the egg stage data were removed, and data from the remaining stages were analyzed, we found that V2/3 was < 0.15. Therefore, the two most stable candidates, RPS18 and EF1A, were recommended for data normalization across different developmental stages (Fig 4B). For the sex, tissue, and temperature comparisons, the first V-values < 0.15 emerged at V2/3, indicating that two reference genes were adequate for normalization under these experimental conditions. Thus, the nominated reference genes were EF1A and HSP90 for sex, EF1A and RPL4 for tissue comparisons, and EF1A and RPS18 for analyses of temperature-mediated effects on gene expression.

Fig 4

Pairwise variation (V) values determined using geNorm.

(A) Across the four indicated comparisons, and (B) for the comparison across developmental stage only.

Discussion

Taken together with previous studies, the present findings clearly showed that reference gene expression stability is condition-specific and, accordingly, no one gene can be used to normalize all RT-qPCR data. This strongly indicates the need to conduct customized reference gene selection for RT-qPCR analyses under the relevant experimental conditions, even for the same species. For example, six papers relating to reference gene selection for the whitefly, Bemisia tabaci, have been published over the past five years [5, 26–30].

The key recommendations for reference gene selection in the MIQE guidelines state that at least two reference genes should be employed in order to avoid biased normalization [31]. The optimal number of reference genes is typically determined by geNorm; based on these results, two reference genes were adequate for the experimental conditions employed in the present study. These results were partly in accordance with those reported previously for other lady beetle species [9–11, 32].

Previous investigations have demonstrated that there is no single reference gene that is applicable under all experimental conditions [1, 2, 5–11]. For example, Actin is a major structural protein that is often used as an internal control without prior validation. However, the present study found that Actin was one of the least stable reference genes under the test conditions; this is consistent with previous studies of four other Coccinellidae species [9–11, 32]. Therefore, we suggest that researchers should initially select reference genes that have been verified within the same family or genus as their target species.

As far as we know, the present study is the first to identify stable RT-qPCR reference genes in P. japonica. Sets of two reference genes were nominated for each experimental condition: EF1A and RPS18 for comparisons between different developmental stages; EF1A and HSP90 for comparisons of female and male adults; EF1A and RPL4 for comparisons of different tissues; and EF1A and RPS18 for investigation of temperature-mediated effects. This study represents an initial move towards building a standardized system for RT-qPCR analyses in P. japonica. This will inform ecological risk assessments of RNAi-based insect control products on P. japonica, and facilitate in-depth functional genomic studies of P. japonica.

Agarose gel electrophoresis of the nine candidate reference genes.

PCR reaction templates: 1) β-actin (Actin); 2) glyceraldehyde-3-phosphate dehydrogenase (GAPDH); 3) elongation factor 1 α (EF1A); 4) α-tubulin (TUBA); 5) ribosomal protein L4 (RPL4); 6) ribosomal protein S18 (RPS18); 7) heat shock protein 90 (HSP90); 8) arginine kinase (ArgK); 9) vacuolar-type H-ATPase subunit A (V-ATPase A). M, DL100 DNA marker.

(TIFF)

Click here for additional data file.

Standard curves of the nine candidate reference genes.

(TIFF)

Click here for additional data file.

Nucleotide sequences of the nine reference genes.

(DOC)

Click here for additional data file.