Identification and Evaluation of Suitable Reference Genes for Normalization of MicroRNA Expression in Helicoverpa armigera (Lepidoptera: Noctuidae) Using Quantitative Real-Time PCR.

More and more studies have focused on microRNAs (miRNAs) expression in the pest Helicoverpa armigera (Lepidoptera: Noctuidae) recently. Quantitative real-time PCR (qRT-PCR) is being widely used in miRNA expression studies. Suitable reference genes are necessary for the correct analysis of results. In this study, 10 candidate genes of H. armigera were selected and analyzed for their expression stability under different biotic and abiotic conditions with 3 statistical methods, including geNorm, NormFinder, and Bestkeeper. Combination the best number of reference genes was calculated by geNorm. One target gene, let-7, was used to validate the selection of reference genes. The suitable candidate reference genes were shown as follows: miR-9 and U6 snRNA for developmental stages, miR-100 and U6 snRNA for larval tissues, miR-100 and miR-305 for adult tissues, miR-9 and miR-279 for parasitic treatment, miR-998 and U6 snRNA for nuclear polyhedrosis virus infection, miR-9 and U6 snRNA for insecticide treatment, miR-92a, miR-100, and miR-279 for temperature treatment, miR-92a, miR-305, and miR-998 for starvation treatment, miR-9 and miR-279 for light treatment, miR-305 and miR-998 for hormone treatment, and there was not one reference gene suitable for all samples. This study could promote future research on miRNAs expression in H. armigera with optimal reference genes under different experimental conditions.

miRNAs are single-stranded, small (about 18-24 nucleotides) genome-encoded noncoding RNAs derived from much longer preprimary transcripts (Bartel 2004). They control the expression of target genes through binding to complementary target ‘‘seed match’’ sites within the 3’ or 5’ untranslated region (UTR) of mRNA targets (Lee et al. 1993, Wightman et al. 1993, Asgari 2013). These miRNAs play some crucial roles in regulating at the posttranscriptional gene expression level lead to mRNA degradation, translational repression, and even mRNA upregulation during almost all known physiological and pathophysiological processes (Bartel 2004, Kim et al. 2009, Hussain et al. 2011, Djuranovic et al. 2012).

To date, several different methods have been developed to measure miRNA expression (Kang et al. 2012). Among these methods, compared with northern blot and microarray technologies are shortage of stability, the fluorescence-based quantitative real time reverse transcriptase PCR (qRT-PCR) is one of the most powerful approaches (Shi and Chiang 2005) because of its high sensitivity, accuracy, specificity, good reproducibility, and linear dynamic range of quantification (Hellemans and Vandesompele 2011, D’haene et al. 2012, Zhang et al. 2012, Serafin et al. 2014). However, the quality of results could be influenced by several nonspecific variables, i.e., the extraction, purification and stability of RNA, and the efficiency of reverse transcription and PCR amplification (Mahoney et al. 2004, Bustin et al. 2009, Bustin et al. 2010). To avoid bias and produce accurate relative results, suitable reference genes are one of the crucial factors (Huggett et al. 2005, Lin and Lai 2013).

Recently, many studies have focused on selecting suitable reference genes for qRT-PCR analysis of miRNAs in plants and human diseases (Kulcheski et al. 2010, Feng et al. 2012, Torres et al. 2013). But, there are very few studies on reference genes of miRNA selection in insects. As we know, only one study on Plutella xylostella (Lepidoptera: Plutellidae) under developmental stages and insecticide-induced stress was reported (Feng et al. 2014). Many studies on miRNAs have adopted U6 snRNA or 5S rRNA (widely expressed in various tissues and cells, and relatively stable in eukarya) as reference genes to normalize miRNA expression without systematic selection (Puthiyakunnon et al. 2013, Wang et al. 2013, Singh et al. 2014). It has been proven that using unreliable reference genes for normalization may lead to large errors and therefore leading to fault results (Andersen et al. 2004, Shi and Chiang 2005). Some previous studies indicated that at least two or three reference genes were necessary for normalizing gene expression data (Vandesompele et al. 2002, Lin and Lai 2013).

The cotton bollworm, Helicoverpa armigera (Hübner) (Lepidoptera: Noctuidae), is one of the most serious polyphagous insect pests of numerous crop plants and has a widespread distribution in Asia, Australia and Africa (Fitt 1989, Matthews and Tunstall 1994). It has caused enormous economic loss in the cotton, corn, vegetable and other crop industries throughout Asia, including China (Wu et al. 2008, Lu et al. 2012). Because of its rapid reproduction and good breeding, H. armigera could be seen as one of the model insects in recent years. Nowadays, there are many publications on the transcriptomics, proteomics, functional gene validation, and miRNAs profiling on H. armigera, and almost all of the studies empirically selected U6 snRNA as reference gene (Ge et al. 2013, Agrawal et al. 2013). Therefore, to insure the validity of further studies on miRNAs in H. armigera, selecting uniform and suitable reference genes under different experimental conditions is very necessary.

In this paper, 10 relatively stable candidate reference genes, including 8 miRNAs and 2 ribosomal RNAs of H. armigera retrieved through high-throughput sequencing, were selected to normalize qRT-PCR data, and 1 development and innate immunity-related target gene (let-7) (Ling et al. 2014) was used to validate optimal reference genes. The expression stability of the reference genes was evaluated under a set of biotic conditions (developmental stages, larval tissues, adult tissues, parasitized by wasps, and nuclear polyhedrosis virus (NPV) infection) and five abiotic stresses (insecticide, temperature, starvation, light, and hormone). Three software packages, including geNorm, NormFinder, and BestKeeper, were then used to evaluate the stability of these genes to enable selection of the most stable normalizer for miRNA qRT-PCR analysis of H. armigera. The present work aims to provide some normalization genes for future expressional and functional research about miRNAs in H. armigera, and to provide references for other insects in selecting of suitable miRNA reference genes.

Materials and Methods

Insects

Helicoverpa armigera had been cultured in our laboratory (IPM laboratory of Entomology department in China Agricultural University, Beijing, China) for more than 10 years. Both the larvae and adults were reared at 26 ± 1°C, with a 14-h light:10-h dark (L14:D10) photoperiod and 70–80% relative humidity (RH). The larvae were reared individually on artificial diet (Bot 1966), and the adults were supplemented with 10% honey solution.

The samples treated with biotic factors and abiotic stresses were all in three biological replicates and immediately stored at −80°C after snap frozen in liquid nitrogen for subsequent total RNA extraction. Six treated insects were collected as one replicate for RNA extraction under different treatments, while more insects were needed in the developmental stages and tissues treatments.

Biotic Factors

Developmental Stages

The samples of our present study contained 400 first-day eggs, 60 first-instar larvae, 40 second-instar larvae, 20 third-instar larvae, 10 fourth-instar larvae, 6 fifth-instar larvae, 6 first-day male pupae, 6 first-day female pupae, 6 first-day male adults, and 6 first-day female adults for each replicate.

Tissues

Twelve tissues were obtained from larvae and adults using a dissection needle and a tweezer in ice-cold PBS solution (140 mM NaCl, 2.70 mM KCl, 10 mM Na2HPO4, 1.80 mM KH2PO4, pH 7.40) (Bear et al. 2010). Seven larval tissues, including brain, epidermis, fat body, midgut, Malpighian tubules, hemocytes and central nervous system were obtained from the second day of fifth instar larvae. Five adult tissues including head, thorax, abdomen, legs, and antenna were dissected from the second day after emerging of male and female, respectively. At least 20 insects were collected for each tissue.

Parasitic Treatment

The parasitic wasp Microplitis mediator (Haliday) (Hymenoptera: Braconidae) was obtained from the Institute of Plant Protection, Hebei Academy of Agriculture and Forestry. The adult wasps were reared at 21 ± 1 °C, 70–80% RH, under a L14:D10 photoperiod.

Paired adults were put into glass tubes and supplied with cotton balls containing 10% honey solution. Then, five females after mating of 3–4 days were put into a glass tube (2 cm × 8.5 cm) with one second instar larva of H. armigera. The larva was removed immediately after the parasitoid made one successful oviposition into it, then replaced by another larva. The parasitized host larvae were fed with artificial diet. Insects were then collected at 0, 24, 48, 72, and 96 h.

NPV Infection

The inoculation of H. armigera with NPV was conducted exactly as in our previous study (Zhang et al. 2014, 2015). Briefly, H. armigera NPV powder was diluted to 106 PIB/ml with sterile water. Then, 10 μl of the solution was pipetted on pieces of the artificial diet. The control group was added with an equal volume of sterile water. The piece of treated diet was provided to one newly molted fourth-instar larva pretreated with overnight starvation in a glass tube (2 cm × 8.5 cm), and normal diet was replenished when the diet with NPV was eaten up. The treated insects were collected at 24, 48, 72, and 96 h, respectively.

Abiotic Stresses

Insecticide-Induced Stress

Two stomach-toxicity insecticides used to control lepidoptera larvae, β-cypermethrin and Spinetoram, which are commonly applied in the management of H. armigera, were chosen in this study. Twenty fourth-instar larvae were treated with the LC50 concentration (Zhang et al. 2014) of each insecticide for 48 h, using similar methods as the NPV infection. Control groups were treated with an equal amount of sterile water.

Temperature-Induced Stress

Fourth-instar larvae were exposed to the temperatures of 4 °C (cold), 26 °C (suitable temperatures), and 40 °C (hot) in culture chambers (Wang and Kang 2003, Li et al. 2013) for 2, 6, and 12 h.

Starvation Treatments

Fourth-instar larvae were placed in glass test tubes without food for 6, 12, 24, 36, and 48 h, respectively.

Light Treatments

Fourth-instar larvae were treated under three light treatments, including all light (L24:D0), all dark (L0:D24) and normal photoperiod (L14:D10). Insects were assembled after 24, 48, 72, and 96 h treatment.

Hormone Treatments

A 0.3 μg of Juvenile hormone or ecdysone (Sigma, St. Louis, MO, USA) dissolved in dimethyl sulfoxide (DMSO, Sigma-Aldrich) per larva were injected into fourth-instar larvae by NANOJECT II (Drummond Scientific Co., Broomall, PA, USA). At 6, 12, and 24 h postinjection, six insects were sampled. The control group was injected with an equal volume of DMSO solution (Zhang et al. 2016).

Total RNA Extraction and Complementary DNA Synthesis

The TRIzol Reagent kit (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from the collected samples, and the quality was evaluated by agarose gel electrophoresis. The RNA samples were considered intact if an 18S band could be clearly observed. Then, an ultraviolet spectrophotometer (Nano-Drop-2000, Thermo Scientific) was used to measure the purity and concentration of total RNA, and only the RNA samples whose A260 as A280 ratios were between 1.9 and 2.1 and A260 as A230 ratios were higher than 2.0 were used for further analysis (Yang et al. 2014).

The resultant total RNA was reverse transcribed using anmiScriptII RT Kit (QIAGEN, Dusseldorf, Germany). First-strand complementary DNA (cDNA) was synthesized from 2 μg of total RNA. All cDNA samples were stored at −80 °C.

Candidate Reference Gene Selection and Primer Design

We got the sequences of miRNAs from the H. armigera parasitized by M. mediator and the control ones through high-throughput sequencing (LC-BIO, Hangzhou, China) in our laboratory. Ten putative stable reference miRNAs were selected, including six miRNAs tested from our Solexa sequencing (bmo-miR-2b-3p_R-2 (miR-2b), bmo-miR-92a (miR-92a), bmo-miR-305-5p_R + 1 (miR-305), bmo-miR-998_R + 2 (miR-998), pxy-miR-6497-p5_2ss7GC19CT (miR-6497), dme-miR-100-5p_R + 1_1ss9AG (miR-100)) which had approximate raw and normalized data, high absolute value of free energy, and medium expression levels (Supp Table 1 [online only]), all above meant the miRNAs were stable before parasitized treatment compared to controls, two miRNAs from previous report of other insects (ppc-miR-9_R + 2 (miR-9) and dme-miR-279-3p_R + 1_2ss20TC21AC (miR-279)) (Liang et al. 2013, Feng et al. 2014), and two widely used ribosomal RNAs (U6 snRNA and 5S rRNA).

The forward primers of all miRNAs were designed according to the fully tested miRNA sequence, and the reverse primer was a universal reverse primer for all miRNAs as liquid from the kit. All of the primer sequences were synthesized by Sangon Biotechnology Co., Ltd (Shanghai, China) (Supp Table 2 [online only]).

Quantitative Real-Time PCR Analysis

A Bio-Rad CFX Connect Real-Time PCR System (Bio-Rad, USA) was used to perform qRT-PCR experiments in 96-wells reaction plates using miScript SYBR Green PCR Kit (QIAGEN, Dusseldorf, Germany). The qRT-PCR reaction mixture was run in a 20-μl volume reaction, consisted of 10 μl Quanti Test SYBR Green PCR Master Mix, 1.0 μl of each specific primer, 1.0 μl of 10× diluted cDNA template, and 7.0 μl RNase-free water. The reaction program was as follows: 95 °C for 15 min, followed by 40 cycles of 94 °C for 15 s, 55 °C for 30 s, and 70 °C for 30 s. Each treatment was prepared in three technical and biological replicates.

A range of series 10-fold dilution of pooled cDNA was used to create the five-point standard curves using the linear regression model (Pfaffl et al. 2004). The regression equation was carried out to calculate the efficiency (E) and correlation coefficient (R2) of each primer pair. The efficiencies (E) of corresponding primers were estimated based on the equation: E = (10[−1/slope] − 1) × 100.

Validation of Reference Gene Selection

To evaluate the validity of the selected reference genes, the transcription levels of microRNA let-7 (Ling et al. 2014) was estimated. Let-7 expression levels were determined in different developmental stages, larval tissues, second instars of H. armigera parasitized by wasps and insecticide treatment of fourth instars with specific primers, respectively. The transcript levels of let-7 were compared among the results when the best (NF1), the worst (NF12) and the optimal recommended combination of genes (NF1-2)/(NF1-3) was used as normalization (Fu et al. 2013). The 2−ΔΔCT algorithm was used to calculate the relative expression levels of the target gene in different samples (Livak and Schmittgen 2001). All the treatments were performed in three biological replicates for RNA extraction and three technical replicates for the qRT-PCR test.

Statistical Analysis

The expression stability of 10 selected reference genes was evaluated with 3 commonly used software tools: geNorm v. 3.5 (Pfaffl et al. 2004), NormFinder version 0.953 (Andersen et al. 2004), and Bestkeeper (Vandesompele et al. 2002). Raw cycle threshold (Ct) values were loaded directly for analysis using Bestkeeper, and linear scaled expression quantities transformed from Ct values were used for the analysis with NormFinder and geNorm software programs (Yang et al. 2014). The stability value (M) calculated by geNorm software (Pfaffl et al. 2004), stability value (SV) calculated by Normfinder software (Andersen et al. 2004) and standard deviation (SD) calculated by BestKeeper (Vandesompele et al. 2002) were used to evaluate and rank the stability of suitable reference genes. The high stability gene has lower numerical values. t-Test was used for statistical analysis between the best combination of reference genes and the worst, using SPSS 17.0 (IBM, GeoHack, NY, USA), and statistical differences were denoted by * (P < 0.05).

Results

Expression Profiles of Selected Reference Genes

The PCR efficiency values of all the primer pairs ranging from 93.59% for miR-279 to 111.26% for miR-2b, and the correlation coefficient (R2) of all primer pairs showed an R2> 0.99 (Supp Table 2 [online only]).

The mean Ct values of the 10 reference genes varied from 17.88 for miR-6497 to 26.37 for 5S rRNA (Supp Table 3 [online only]). The remaining reference genes were expressed at moderate levels, with mean Ct values of 20.90, 26.09, 23.53, 26.18, 25.21, 22.01, 21.69, and 19.83 for miR-2b, miR-92a, miR-305, miR-998, miR-100, miR-9, miR-279, and U6 snRNA, respectively (Fig. 1; Supp Table 3 [online only]). The mean Ct value of the target gene let-7 was 28.93. miR-92a exhibited the lowest dispersion (4.29 cycles) over all samples, and miR-100 (9.91 cycles) displayed the highest dispersion (Supp Table 3 [online only]; Fig. 1).

Fig. 1.

Average qPCR cycle threshold (Ct) values of candidate reference genes tested in H. armigera under different treatments. The black dot represents the outliers of replicated samples, while whiskers indicate the standard deviation of the mean.

Expression Stability of Selected Reference Genes under Biotic Conditions

The most stable genes recommended by geNorm for gene profiling among developmental stages were miR-279, miR-998 and U6 snRNA, and Normfinder determined miR-9, miR-2b, miR-6497, while Bestkeeper allocated U6 snRNA, miR-9, and miR-92a as the three best-suited genes. The lowest three ranked genes were miR-100, 5S rRNA, and miR-305, and the three miRNAs were predicated as the least stable genes in all of the biological samples by geNorm (Table 1).

Table 1.

Expression stability of the candidate reference genes under different biotic conditions

Biotic condition	Reference gene	geNorm		Normfinder		Bsetkeeper
		M	Rank	SV	Rank	SD	R	Rank
Developmental stages	miR-2b	1.295	6	0.039	2	1.042	0.872*	7
	miR-92a	1.114	4	0.052	4	0.732	0.015	3
	miR-305	1.672	8	0.075	8	1.757	0.775*	8
	miR-998	0.654	1	0.064	7	0.839	0.070	4
	miR-6497	1.375	7	0.044	3	0.894	0.693*	6
	miR-100	2.140	10	0.095	9	2.525	0.737*	10
	miR-9	1.186	5	0.024	1	0.667	0.699*	2
	miR-279	0.654	1	0.056	6	0.854	0.218	5
	U6 snRNA	0.936	3	0.054	5	0.506	−0.102	1
	5S rRNA	1.917	9	0.097	10	2.220	0.800*	9
Larval tissues	miR-2b	0.553	4	0.017	1	0.464	0.704*	3
	miR-92a	1.017	6	0.043	5	1.068	0.500	7
	miR-305	0.458	3	0.017	2	0.546	0.703*	4
	miR-998	1.402	9	0.061	8	1.539	0.875*	9
	miR-6497	0.810	5	0.062	9	0.960	0.509	6
	miR-100	0.236	1	0.025	3	0.289	0.305	2
	miR-9	1.626	10	0.095	10	1.835	0.411	10
	miR-279	1.127	7	0.045	6	1.143	0.799*	8
	U6 snRNA	0.236	1	0.030	4	0.203	0.038	1
	5S rRNA	1.267	8	0.055	7	0.930	0.381	5
Adult tissues	miR-2b	0.417	1	0.005	1	0.858	0.947*	5
	miR-92a	0.964	6	0.044	5	0.806	0.343	4
	miR-305	0.484	3	0.006	2	0.784	0.836*	3
	miR-998	1.807	10	0.102	9	2.326	0.465	10
	miR-6497	1.550	9	0.117	10	1.338	0.082	8
	miR-100	0.417	1	0.009	3	0.707	0.920*	1
	miR-9	0.843	5	0.047	6	1.242	0.708*	7
	miR-279	1.361	8	0.100	8	1.787	0.827*	9
	U6 snRNA	1.161	7	0.068	7	1.075	0.387	6
	5S rRNA	0.637	4	0.028	4	0.755	0.614*	2
Parasitic	miR-2b	0.156	1	0.011	5	0.148	0.078	2
	miR-92a	0.388	8	0.018	8	0.415	0.623	9
	miR-305	0.217	4	0.008	3	0.178	0.467	4
	miR-998	0.269	5	0.000	1	0.260	0.894*	6
	miR-6497	0.571	10	0.056	10	0.869	0.948*	10
	miR-100	0.448	9	0.026	9	0.335	−0.683	8
	miR-9	0.179	3	0.005	2	0.123	0.671	1
	miR-279	0.156	1	0.014	6	0.158	−0.200	3
	U6 snRNA	0.347	7	0.017	7	0.255	0.236	5
	5S rRNA	0.320	6	0.009	4	0.288	0.629	7
NPV infection	miR-2b	0.487	3	0.018	2	0.466	0.680*	5
	miR-92a	0.600	6	0.020	6	0.429	0.575	3
	miR-305	0.562	5	0.020	5	0.550	0.851*	7
	miR-998	0.368	1	0.010	1	0.450	0.810*	4
	miR-6497	0.709	9	0.040	10	0.527	0.713*	6
	miR-100	0.642	7	0.027	7	0.568	0.359	8
	miR-9	0.548	4	0.019	4	0.593	0.882*	9
	miR-279	0.674	8	0.033	8	0.327	−0.241	1
	U6 snRNA	0.368	1	0.018	3	0.423	0.712*	2
	5S rRNA	0.788	10	0.037	9	0.855	0.763*	10

M, stability value; SD, standard deviation; SV, stability value; R, Pearson correlation coefficient.

P ≤ 0.001.

Larval Tissues

For the larval tissues study, the top two ranked genes, U6 snRNA and miR-100, were identified as the most stable genes by geNorm and Bestkeeper, while miR-2b was ranked in the top position according to Normfinder, and miR-9 was identified as the least stable gene by all three methods (Table 1).

Adult Tissues

The rankings of the best-suited reference genes determined by the geNorm and Normfinder programs were similar for adult tissues. The two methods identified the top three ranked genes as miR-2b, miR-100, and miR-305, while Bestkeeper recommended miR-100, 5S rRNA, and miR-305 as the three best-suited genes. Both geNorm and Bestkeeper identified miR-998, but Normfinder ranked miR-6497 as the worst stable gene (Table 1).

For parasitic treatment, miR-2b, miR-9, and miR-279 were recognized as the best-suited genes by geNorm and Bestkeeper, while Normfinder ranked miR-998, miR-9, and miR-305 as the best ones. MiR-9 and miR-279 were exactly confirmed as the most stable genes in all of the biological samples. Meanwhile, the least steady gene was ranked as miR-6397 by all the three programs (Table 1).

Both geNorm and Normfinder determined miR-998, U6 snRNA, and miR-2b as the best-suited genes, while Bestkeeper software obtained miR-279, U6 snRNA, and miR-92a as the best ones. 5S rRNA was ranked as the lowest stable gene by geNorm and Bestkeeper, while Normfinder ranked miR-6497 (Table 1).

Expression Stability of Selected Reference Genes under Abiotic Conditions

Insecticide Treatment

For insecticide treatment, the top ranked three genes were miR-305, U6 snRNA, and miR-9 by geNorm, miR-998, miR-92b, and miR-9 by NormFinder, while Bestkeeper ranked miR-9, U6 snRNA, and miR-998. The worst-suited gene calculated by the three software packages was miR-6497 (Table 2).

Table 2.

Expression stability of the candidate reference genes under different abiotic conditions

Abiotic condition	Reference gene	geNorm		Normfinder		Bsetkeeper
		M	Rank	SV	Rank	SD	R	Rank
Insecticide	miR-2b	0.435	7	0.029	8	0.517	0.404	8
	miR-92a	0.366	6	0.003	2	0.276	0.893*	5
	miR-305	0.187	1	0.016	7	0.252	−0.692	4
	miR-998	0.329	5	0.003	1	0.181	0.752	3
	miR-6497	0.759	10	0.069	10	1.027	0.034	10
	miR-100	0.281	4	0.012	5	0.277	0.111	6
	miR-9	0.222	3	0.007	3	0.134	−0.272	1
	miR-279	0.483	8	0.029	9	0.540	0.662	9
	U6 snRNA	0.187	1	0.009	4	0.136	−0.638	2
	5S rRNA	0.532	9	0.013	6	0.418	0.797	7
Temperature	miR-2b	0.276	1	0.023	7	0.439	0.567	4
	miR-92a	0.303	4	0.013	2	0.305	0.654*	1
	miR-305	0.409	7	0.010	1	0.490	0.860*	7
	miR-998	0.359	5	0.021	6	0.485	0.632*	6
	miR-6497	0.784	10	0.060	10	1.039	0.575	10
	miR-100	0.294	3	0.018	3	0.389	0.539	3
	miR-9	0.383	6	0.019	4	0.454	0.718*	5
	miR-279	0.276	1	0.020	5	0.328	0.492	2
	U6 snRNA	0.521	8	0.028	8	0.566	0.633*	8
	5S rRNA	0.664	9	0.032	9	0.919	0.756*	9
Starvation	miR-2b	0.416	5	0.031	7	0.777	0.738	7
	miR-92a	0.166	1	0.005	1	0.682	0.981*	4
	miR-305	0.166	1	0.009	3	0.731	0.959*	6
	miR-998	0.273	3	0.005	2	0.579	0.982*	2
	miR-6497	0.901	10	0.079	10	1.098	0.518	10
	miR-100	0.330	4	0.014	4	0.712	0.881*	5
	miR-9	0.486	6	0.028	6	0.457	0.585	1
	miR-279	0.517	7	0.035	8	0.591	0.486	3
	U6 snRNA	0.643	8	0.037	9	1.034	0.966*	9
	5S rRNA	0.730	9	0.026	5	0.926	0.792*	8
Light	miR-2b	0.335	3	0.023	5	0.378	0.326	4
	miR-92a	0.524	7	0.014	3	0.502	0.756*	8
	miR-305	0.468	6	0.013	2	0.461	0.860*	7
	miR-998	0.360	4	0.010	1	0.231	0.537*	1
	miR-6497	0.696	9	0.052	9	0.694	0.616*	9
	miR-100	0.381	5	0.028	8	0.442	−0.025	6
	miR-9	0.303	1	0.017	4	0.314	0.363	3
	miR-279	0.303	1	0.026	6	0.298	−0.169	2
	U6 snRNA	0.582	8	0.026	7	0.397	0.605*	5
	5S rRNA	0.870	10	0.056	10	1.297	0.656*	10
Hormone	miR-2b	0.146	1	0.025	7	0.354	0.015	8
	miR-92a	0.337	7	0.008	1	0.277	0.565	5
	miR-305	0.298	6	0.008	2	0.204	0.521	1
	miR-998	0.166	3	0.017	6	0.274	−0.064	3
	miR-6497	0.785	10	0.082	10	1.444	0.793*	10
	miR-100	0.181	4	0.014	4	0.221	−0.040	2
	miR-9	0.225	5	0.016	5	0.296	0.277	6
	miR-279	0.146	1	0.025	8	0.349	−0.197	7
	U6 snRNA	0.393	8	0.013	3	0.274	0.488	4
	5S rRNA	0.582	9	0.035	9	0.923	0.866*	9

M, stability value; SD, standard deviation; SV, stability value; R, Pearson correlation coefficient.

P ≤ 0.001.

Temperature Treatment

MiR-279, miR-2b, and miR-100 were identified as the most stable genes by geNorm, while NormFinder identified miR-305, miR-92a, and miR-100. BestKeeper identified miR-92a, miR-279, and miR-100 as the most stable genes. All three computational programs identified miR-6497 and 5S rRNA as the least stable genes (Table 2).

Starvation Treatment

For the starvation treatment, the three top ranked genes were similar between geNorm and NormFinder, although the rank order was slightly altered. The top three genes identified by geNorm were miR-92a, miR-305, and miR-998. MiR-92a was identified as the most stable gene by NormFinder, while miR-998 and miR-305 were ranked in the second and third position, respectively. However, Bestkeeper analysis found miR-9 was selected as the most suitable normalization factor for qRT-PCR normalization, followed by miR-998 and miR-279. MiR-6497 was identified as the lowest stable gene by all of the three programs (Table 2).

Light Treatment

The worst-suited genes exhibited by the three software packages were 5S rRNA and miR-6497 under light treatment, but the best-suited genes identified by the three programs were not the same. geNorm ranked the most stable genes as miR-9, miR-279, and miR-2b, NormFinder ranked miR-998, miR-305, and miR-92a, while Bestkeeper ranked miR-998, miR-279, and miR-9 (Table 2).

Hormone Treatment

For the hormone treatment, the most stable genes determined by geNorm were miR-2b, miR-279, and miR-998, Normfinder recommended miR-92a, miR-305, and U6 snRNA, while Bestkeeper recognized miR-305, miR-100, and miR-998 as the three best-suited genes. All the three software packages identified miR-6497 and 5S rRNA as the least stable genes (Table 2).

Combination the Best Number of Reference Genes

According to Figure 2, the V2/3 values of parasitic treatment and abiotic stress were all below 0.15. For the different larval tissues, three best selected reference genes should be used due to the fact that V3/4 was below the proposed 0.15 value. But for developmental stages and adult tissues, there was no V value that was below 0.15.

Fig. 2.

Pairwise variation analysis for an accurate normalization in H. armigera. The pairwise variation (V/V + 1) was analyzed by geNorm software between the normalization factors NFn and NFn + 1 to determine the optimal number of reference genes. Each pairwise variation value is compared with 0.15, below of which the inclusion of an additional reference gene is not required.

Summarized Best Genes under Different Conditions

According to the different results from three packages, we assumed that a combination of different mathematical models enabled a better evaluation of the most stable reference genes. Thus, we summarized most suitable genes under different conditions through the same coincident results provided by three different programs (Table 3).

Table 3.

Preferable reference genes in H. armigera under different experimental conditions

Experimental conditions		Preferable reference genes
Biotic factors	Developmental stages	miR-9	U6 snRNA
	Larval tissues	miR-100	U6 snRNA
	Adult tissues	miR-100	miR-305
	Parasitic	miR-9	miR-279
	NPV injection	miR-998	U6 snRNA
Abiotic factors	Insecticide	miR-9	U6 snRNA
	Temperature	miR-92a	miR-100	miR-279
	Starvation	miR-92a	miR-305	miR-998
	Light	miR-9	miR-279
	Hormone	miR-305	miR-998

For various development stages, when using the best reference gene or the recommended two and three most stable references, let-7 transcript levels were higher in pupa and adult periods compared with larval stages (Fig. 3A). Although using NF (1-2) (miR-279, miR-998) is significantly higher than the ones using NF1 (miR-279) and NF (1-3) (miR-279, miR-998, U6 snRNA), their expression tendency were consistent in different development. However, when the most unstable gene NF12 (miR-100) was used, the let-7 transcript levels were very low in all sample developmental stages, and no evident difference was detected. Moreover, the expression level of let-7 normalized against the combination of two best reference genes was significantly different from the least stable reference gene (P < 0.05; Fig. 3A).

Fig. 3.

Validation of reference gene selection. Relative Expression levels of target gene let-7 (A) in 10 developmental stages (B) different larval tissues (C) under parasitic (D) and two insecticides treatment. B: brain, FB: fatbody, MG: midgut, MT: Malpighian tubules, EP: epidermis, BC: β-cypermethrin, SP: Spinetoram. NF1, using the single one best reference gene for normalization, NF (1-2), using two best reference genes, NF (1-3), using three best reference genes, NF (12), using the least stable reference gene. Asterisks indicate significant differences in the expression of the target gene based on three biological replicates (P < 0.05, t-test).

The expression profiles of let-7 in different tissues were similar when normalized using NF1 (miR-100), NF (1-2) (miR-100, U6 snRNA), NF (1-3) (miR-100, U6 snRNA, miR-305). The transcript levels of let-7 normalized using NF (1-3) and NF12 (miR-9) were significantly different in larval tissues (P < 0.05; Fig. 3B).

For parasitic treatment, the expression profiles of let-7 were similar when normalized using NF1 (miR-9) and NF (1-2) (miR-9, miR-279). The expression levels normalized using NF12 (miR-6497) were 1.36-fold lower than using NF (1-2) in the first day after parasitizing, but they were increased in the other days after parasitizing (P < 0.05; Fig. 3C).

For insecticide treatment, when normalized against NF12 (miR-6497), compared with NF (1-2) (miR-9, U6 snRNA), the expression of let-7 was increased by 8.77-fold for the treatment of β-cypermethrin and 2.98-fold for Spinetoram (P < 0.05). Using the NF1 (miR-9) for normalization, similar relatively expression levels were observed just like using NF (1-2) (Fig. 3D).

Discussion

Rearing H. armigera larvae using artificial diet were convenient and economic under experiment condition, reference genes must have stable expression under different conditions, so we think our experiment using artificial diet would not affect the expression of reference genes.

In our study, the reaction efficiencies of reference genes indicated that the selected genes can be used for further analysis. And the distributions of the mean Ct values of the reference genes indicated that the quantity of expression of those genes was relatively intermediate across different samples.

It is difficult to apply a universal appropriate reference gene for all samples when analyzed by the three different programs. This is probably because that each program has a different algorithm (Andersen et al. 2004, Mallona et al. 2010, Mafra et al. 2012). The principle of geNorm software is that the transcript ratio was the same of the two best reference genes. It ranks the reference genes via calculating the expression stability value (M) and recommends the optimal number of reference genes via analyzing the pairwise variation (V). NormFinder calculated an expression stability value by using a different mathematical model. BestKeeper is an excel-based tool that is based on the geometric mean of Ct values and PCR amplification efficiency. It determines the best reference gene according to SD values, the probability (P) value and the Pearson correlation coefficient (r) (Pfaffl et al. 2004). Expression stability of the candidate reference genes under different biotic and abiotic conditions showed that the rank of most suitable genes obtained from BestKeeper was slightly different from other programs, this is commonly seen in these kind of reference gene selection (Yang et al. 2014, Zhang et al. 2014). The reason might be that, compared with another two software packages, BestKeeper only analyzes the stability of reference genes individually instead of considering the pairwise variation between two reference genes (Silver et al. 2006, Teng et al. 2012), and it usually represents the average of the best suitable reference genes.

Furthermore, previous studies have reported that compared to the use of a single reference gene, more reference genes could increase the accuracy of quantitation (Vandesompele et al. 2002, Haller et al. 2004). Vandesompele et al. (2002) recommended to determine the best combination of candidate housekeeping genes by calculating the normalization factor (NF) with geNorm. If the pairwise variation (V/) was below 0.15, adding n + 1 gene would have no great influence for the results (Vandesompele et al. 2002). However, other factors, for example cost and convenient operation, especially the limit of template, have to be considered to quantify the number of reference genes. In general, two or three stable genes should be used for further data analysis.

The test of relative expression of the target gene using different combination of reference genes proved that, inappropriate reference genes could significantly alter the outcome of miRNA quantitation and will result in false interpretation. Therefore, the appropriate reference gene for normalization was crucial for accurate estimation of target gene expression. Our study will provide useful information on miRNA target gene expression and function in H. armigera and other insects in future studies.

Conclusions

In this study, 10 candidate reference genes of H. armigera were systematically evaluated for their expression stability across different biotic and abiotic experimental conditions. Results were as follows: miR-9 and U6 snRNA for developmental stages; miR-100 and U6 snRNA for larvae tissues; miR-100 and miR-305 for adult tissues; miR-9 and miR-279 for parasitic treatment; miR-998 and U6 snRNA for nuclear polyhedrosis virus infection; miR-9 and U6 snRNA for insecticide treatment; miR-92a, miR-100, and miR-279 for temperature treatment; miR-92a, miR-305, and miR-998 for starvation treatment; miR-9 and miR-279 for light treatment; miR-305 and miR-998 for hormone treatment. Our work indicated that the expression stability of selected reference genes should be evaluated in different experimental conditions.

Supplementary Data

Supplementary data are available at Journal of Insect Science online.

Click here for additional data file.