Identification and Validation of Reference Genes for RT-qPCR Normalization in Mythimna separata (Lepidoptera: Noctuidae)

Mythimna separata is a major agricultural pest with seasonal migrating trait in China. Formation and regulation mechanism of migration behavior has resulted in a large number of fundamental researches involving quantitative studies of gene expression in this species. Using appropriate reference gene is critical in RT-qPCR data normalization. A comprehensive study on the reference genes in M. separata is lacking. In this paper, expression stabilities of ten candidate reference genes were evaluated in M. separata under various biotic and abiotic conditions by employing four different software geNorm, NormFinder, BestKeeper, and the comparative ΔCT method. The comprehensive stabilities ranking of these genes were suggested by RefFinder. PKG as a target gene was employed to justify the number of reference genes in four larval tissues and two photoperiod treatments. Results demonstrate that the first three most stable genes were as follows: EF, CypA, and β-TUB for developmental stages; EF, CypA, and RPL12 for larval tissues; EF, TBP, and β-TUB for adult tissues. RPL12, β-TUB, and EF for densities; EF, RPL12, and GAPDH for photoperiod treatments; β-TUB, EF, and ATPase for temperature treatments. Stable reference gene combinations may reduce bias in normalization. This work provides for the first time a comprehensive list of appropriate reference genes and facilitates future studies on gene function of M. separata.


Introduction
The oriental armyworm, Mythimna separata (Lepidoptera: Noctuidae), is an important agricultural pest in Asia and Australia [1]. In China, it is widely distributed except in the far northwest (Xinjiang). The larvae of the armyworm can inflict devastating damage to more than 104 species of plants in 16 families including food crops (wheat, rice, foxtail millet, and maize) [2]. Each year, the nocturnal moths engage in a seasonal long-distance migration between southern and northern China [3]. M. separata has been extensively studied as a model migratory species in China [3][4][5]. The research contents mainly involved environmental, physiological, hormonal, genetic, and molecular factors [3]. The fifth and sixth instars of the larvae (the gluttony period) are the very important stages in which larvae consume vast food to accumulate sufficient nutrition and energy to meet remaining development and migration requirements. Environmental factors such as temperature [6,7], photoperiod [8], and larval density [9][10][11] can trigger or inhibit the onset of migration. Moths response to larval environmental conditions [12] and are most likely to take off on a migratory flight the first or second night after adult eclosion [13]. Therefore, molecular studies towards environmental factors and specific developmental stages are the focus of migration research in M. separata. Exploring target gene expression profiles in tissues from various environmental conditions or developmental periods will promote a greater understanding of the migration regulation mechanism. Reports on molecular regulation mechanisms of migration for M. separata were relatively scarce.
The real-time quantitative polymerase chain reaction (RT-qPCR), as a reliable and rapid method, has been widely applied in gene profiling or expression analysis in organisms. Compared with classic quantitative PCR, this approach offers distinct advantages including high sensitivity, specificity, accuracy, reproducibility, and requirement for less post-PCR processing [14,15]. However, some variables may influence the accuracy of RT-qPCR, such as total RNA quality, the efficiency of reverse transcription, and primer transcription efficiency [16][17][18][19].
In the current paper, 10 commonly used reference genes were examined in M. separata to normalize RT-qPCR data. The expression profile of protein kinase G gene (PKG) was applied to further validate the selected reference genes. The aims are to identify and evaluate the expression stability of the reference genes under different biotic (development stage, tissue, and density) and abiotic (photoperiod and temperature) conditions. This will facilitate future researches on gene regulation involved in M. separata migration.

Insects.
The insects were initially obtained from a culture at the Biorational Pesticides Research and Development Center, Northwest A&F University, China, and maintained under a L14: D10 regime at 23 ± 1 ∘ C and 70 ± 10% relative humidity (RH). In each treatment, the colonies were derived from one batch of eggs laid by a single male/female pair. The larvae were immediately reared in one 850-ml jar from the day of hatching with fresh maize seedling. The samples were dissected under light CO2 aesthesia and immediately immersed in liquid nitrogen and stored at −80 ∘ C until used.

Developmental Stages.
Samples of brain were collected from 8 developmental stages, including 3 larval stages (from fourth instar to sixth instar; L4-L6), 1 prepupal stage (PP), 3 pupal stages (P1; P5; P9), and 1 adult (A0). In the larval stage, two-day-old larvae from each instar stage were selected in this paper. The newly formed pupa is designated as day P0. Subsequent days of pupal development are counted as P1-P11. The freshly formed adult is termed A0. The samples were in triplicate, 20 larvae in each repetition, and stored as just described.

Larval Tissues.
Four larvae tissues (brains, epidermises, fat bodies, and alimentary canals) were obtained from second day larvae of 5th instar stage. Each sample was repeated three times, 15 larvae in each repetition. All samples were stored as mentioned earlier.

Adult Tissues.
Six adult tissues, including male brains, female brains, testes, ovaries, male fat bodies, and female fat bodies, were dissected from one-day-old adults. For each tissue, 15 insects were collected. Each sample was repeated three times. All samples were stored as mentioned earlier.

Densities.
The larvae were immediately separated after hatching and formed four larval densities regimes including 1, 10, 20, and 30 larvae/850 ml jar and then were reared with fresh maize seedlings under a L14:D10 cycle at 23 ± 1 ∘ C and 70 ± 10% relative humidity (RH). Each density regime had three repetitions, 20 larvae in each repetition. Brains were obtained from the fifth instar larvae and stored as described above.

Photoperiod Treatments.
Colonies newly hatched were reared at a density of 10 larvae per 850-ml jar and subjected to two different photoperiod regimes (8L: 16D and 16L: 8D) at 23 ± 1 ∘ C and 70 ± 10% relative humidity (RH). Each photoperiod treatment had three repetitions, 20 larvae in each repetition. Brains were obtained from 5th instar larvae and stored as mentioned above.

Temperature Treatments.
As photoperiod treatment, larvae after hatching were reared at a density of 10 larvae/ 850 ml jar and were exposed at two temperature regimes (18 ∘ C and 25 ∘ C) under 70 ± 10% RH and L14: D10 photoperiod. Each temperature treatment had three repetitions, 20 larvae in each repetition. Brains were obtained from 5th instar larvae and stored as mentioned previously.

RNA Extraction and cDNA Synthesis.
Total RNA was extracted with RNAiso Plus (Takara Bio Inc., Dalian, China) according to the manufacturer's instructions. First-strand cDNA was synthesized from 1 g of total RNA using the PrimeScript5RT reagent Kit (Perfect Real-Time) (Takara Bio Inc.) with gDNA Eraser following the manufacturer's protocols. The synthesized cDNA was stored at −80 ∘ C.

Real-Time RT-qPCR and Expression Stability Analysis.
RT-qPCR was conducted in an iQ65 Multicolor Real-time PCR Detection system (Bio-Rad, Hercules, CA, USA). The amplification reactions were performed in a volume of 20 L containing 8 diluted cDNA, 10 L SYBR Premix EX Taq polymerase (Takara Bio Inc.), and 1 L of each 10 M primer. The cycling conditions were 95 ∘ C for 30s, 40 cycles of 95 ∘ C for 30 s, and 60 ∘ C for 30 s, followed by a dissociation-curve program from 60 ∘ C to 95 ∘ C for verification of PCR amplicons specificity. Standard curves were created with 10-fold dilution series of pooled cDNA for each treatment using the linear regression model. For both biotic and abiotic conditions, all reactions were performed in three technical replicates. All biological replicates were used to calculate the average CT value.
The expression stabilities of the candidate reference genes were evaluated by geNorm [34], NormFinder [35], BestKeeper [36], and the comparative ΔCT method [37]. RefFinder, a web-based analysis tool (http://150.216.56.64/ referencegene.php), was finally applied to evaluate and screen the optimal reference genes by integrating the results obtained from the above four programs. In addition to the stability evaluation, the optimal number of reference genes is determined with the geNorm software by calculating a normalization factor (NF). geNorm calculates the pairwise variation V n /V n+1 between two sequential normalization factors NF n and NF n+1 containing and increasing number of reference genes. If V n /V n+1 < 0.15 the inclusion of an additional reference gene is not required and the recommended number of reference genes is given by n [34]. Given the case that all pairwise variation values are above the proposed 0.15 cut-off value, the optimal number of reference genes for normalization, the lowest pairwise variation is recommended [34,38].
cGMP-dependent protein kinase (PKG) gene was used to evaluate the validity of selected reference genes under biotic (larval tissues) and abiotic (photoperiod) conditions. All the experiments were performed in triplicate. Its expression levels were determined according to the CT value based on the 2 −ΔΔCT method [39]. The relative expression levels in the four larval tissues (brains, epidermises, alimentary canals, and fat bodies) were analyzed using a one-way analysis of variance (ANOVA) followed by Tukey's post hoc test. Significance values were set to P < 0.05. (Table 1) were selected to normalize the gene expression levels in M. separata. The qPCR primer pairs for these reference genes were specific according to the analysis of a single sharp peak in melting curve and the presence of specific bands of the expected size in the agarose gel electrophoresis. Each primer pairs amplification efficiency ranged from 96.8% to 108.9% and correlation coefficients (R2) were not less than 0.99 (Table 1), which were within an acceptable range for qPCR.

Primer Specificity and PCR Efficiency. Ten candidate reference genes
A broad range of CT values from 10.71 (28S) to 25.91 (TBP) in RT-qPCR exhibited all reference genes expression levels across all treatments ( Figure 1). 18S (mean CT value) and 28S (mean CT value) were expressed at the highest levels and EF (mean CT value) and TBP (mean CT value) at the lowest levels. The six remaining reference genes were expressed at moderate levels (mean CT values of 18.87, 19.84, 19.83, 21.57, 15.82, and 18.7 for GAPDH, RPL12, -TUB, ATPase, -ACT, and CypA, respectively).

Stability of Candidate Reference Genes under
Biotic Conditions The stability ranking order of the first 3 most stable genes obtained from four programs were inconsistent ( Table 2). The least stable genes determined by four programs were 18S, 28S, and -ACT. According to RefFinder, the stability ranking of the reference genes from the most stable to the least stable across different developmental stages is as follows: (Figure 2(a)). The geNorm analysis showed that all pairwise variation values were above the proposed 0.15 cut-off value (Figure 3(a)). The lowest pairwise variation was shown at V 6/7  (the variation between the normalization factors of six genes in relation to seven genes); thus, six genes are recommended as reference genes for normalization. It is relatively impractical to use excessive numbers of endogenous control genes for normalization. The use of the three most stable control genes was considered to be adequate for normalization of RT-qPCR [34,38]. To verify this recommendation, the correlation of NF values between the three most stable genes and the optimal number of genes was calculated. As shown in Figure 3(b)-(A), there is a good correlation between two NF measures (proposed number, three; the theoretical optimal number, six) for developmental stages (r = 0.972). This result supports that the three most stable internal control genes are sufficient for an accurate normalization of developmental stages data. Thus, the normalization with the combination EF, CypA, and RPL12 should be suggested by geNorm for a suitable normalization in the different developmental stages. And geNorm is basically in line with RefFinder identifying three out of ten most stable genes (Figure 2(a)).

Larval Tissues. 18S
and -ACT were identified as the least stable genes by the four programs in different tissues (Table 2(a)). The four programs revealed that EF was the most stable gene (Table 2(a)). According to RefFinder, the stability ranking of reference genes from the most stable to the least stable gene under various tissues was as follows: (Figure 2(b)). In contrast with developmental stages, all pairwise variation values obtained by the geNorm were below the proposed 0.15 cut-off value (Figure 3(a)). The normalization with the combination EF and RPL12 was proposed by geNorm in the various tissues samples. It is roughly identical to the analysis result of RefFinder software calculating three most stable genes.

Adult
Tissues. -ACT was identified as the least stable genes by the three programs in different tissues (Table 2(a)) except BestKeeper. The four programs, except BestKeeper, revealed that EF and TBP were the most stable genes (Table 2(a)). According to RefFinder, the stability ranking of reference genes from the most stable to the least stable gene under various tissues was as follows: (Figure 2(c)). V 2/3 pairwise variation value obtained by the geNorm was below the proposed 0.15 cutoff value (Figure 3(a)). Therefore, the normalization with the combination EF and TBP was proposed by geNorm in the various tissues samples. It is consistent with the calculational result of RefFinder software proposing three out of ten most stable reference genes.

Densities.
It is similar to the results determined by NormFinder in which the stability ranking identified by geNorm revealed that RPL12, -TUB, and 28S were the most stable genes (Table 2(b)). CypA and TBP were the least stable genes identified by four programs (Table 2(b)). According to RefFinder, the stability ranking of the reference genes from the most stable to the least stable gene across density was as follows: (Figure 2(d)). The geNorm analysis showed that all pairwise variation values were below the proposed 0.15 cut-off value (Figure 3(a)). The combination of two genes, RPL12 and -TUB, is recommended by geNorm for a suitable normalization in the density samples. It is similar to the analysis result by RefFinder software recommending the two most stable control genes.

All Biotic Samples.
The stability ranking of the reference genes proposed by NormFinder was closely similar to the result revealed by the ΔCT method (Table 2(b)), which identified -TUB and EF as the most stable genes. Three genes ( -ACT, 18S, and 28S) were the least stable genes revealed by four programs (Table 2(b)). The geNorm showed RPL12 and CypA were the most stable genes, whereas, EF and TBP determined by BestKeeper were the most stable genes (Table 2(b)). According to RefFinder, the stability ranking of the reference genes from the most stable to the least stable gene under all biotic samples was as follows: (Figure 2(e)). All pairwise variation values calculated by the geNorm were above the proposed 0.15 cut-off value (Figure 3(a)). According to the lowest pairwise variation appearing at V 5/6 , five genes are proposed as reference genes for normalization. As for the developmental stages, the correlation of NF values between NF 3 and NF 5 was analyzed. As shown in Figure 3(b)-(B), there is a good correlation between two NF measures for all biotic samples (r = 0.936). This result supports that the three most stable reference genes may meet an accurate normalization for all biotic samples. Thus, the normalization with RPL12, CypA, and EF was proposed by geNorm across all biotic samples. This result is in general accordance with that of RefFinder.

Photoperiod Treatments.
Gene stability ranking order determined by BestKeeper was closely similar to the results revealed from the ΔCT method, which revealed EF, RPL12, and CypA to be the most stable genes ( Table 3). The most stable genes identified by geNorm were GAPDH, -TUB, and CypA, whereas those identified by NormFinder were EF, -ACT, and ATPase ( Table 3). The four programs revealed 18S and 28S to be the least stable genes (Table 3). According to RefFinder, the stability ranking of the reference genes from the most stable to the least stable gene across photoperiod was as follows: (Figure 2(f)). For this treatment, the first V-value < 0.15 showed at V 2/3 (Figure 3(a)), suggesting that two reference genes were sufficient for reliable normalization. Thus, the normalization with the combination of GAPDH and -ACT was proposed by geNorm in all photoperiod treatments.

Temperature.
Gene stability ranking results of most stable genes determined by geNorm was very similar to the results obtained from NormFinder, which showed -TUB, ATPase, and EF to be the most stable genes (Table 3). BestKeeper and the ΔCT method were closely similar in the ranking results, which exhibited that RPL12, EF, and -TUB were the most stable genes (Table 3). 28S and 18S determined by the four programs were the least stable genes. According to RefFinder, the gene stability ranking from the most stable to the least stable across temperature was as follows: (Figure 2(g)). The geNorm analysis showed that the pairwise variation values of V 5/6 , V 6/7 , and V 7/8 were the cut-off value of 0.15 (Figure 3(a)). According to the lowest pairwise variation shown at V 5/6 , five genes are proposed as reference genes for normalization. As previously described, the correlation of NF values between NF 3 and NF 5 was calculated. As shown in Figure 3(b)-(C), there is a good correlation between two NF measures for all biotic samples (r = 0.998). It supports that the three most stable reference genes may suit an accurate normalization for temperature treatments. Thus, the normalization with -TUB, EF, and ATPase was recommended by geNorm across temperature treatments. This result is consistent with that of RefFinder.

All Abiotic Samples.
The results of gene stability ranking obtained from BestKeeper and the ΔCT method were similar, which revealed that EF, RPL12, and -TUB were the most stable genes (Table 3). Four programs identified 28S and 18S as the least stable genes. -ACT and CypA determined by geNorm were the most stable genes; however, -TUB and ATPase recommended by NormFinder were the most stable genes (Table 3). According to RefFinder, the gene stability ranking from the most stable to the least stable across temperature was as follows: (Figure 2(h)). Pairwise variation values of V 2/3 , V 3/4 , and V 4/5 calculated by the geNorm were below the proposed 0.15 cut-off value (Figure 3(a)). Thus, the combination of two reference genes, CypA and -ACT, proposed by geNorm was suitable for normalizing RT-qPCR data in all abiotic samples.

All Samples.
The stability ranking results determined by NormFinder and the ΔCT method were similar. These programs revealed that -TUB, EF, and CypA were the most  stable genes (Table 4). 18S and 28S were identified by the four programs as the least stable genes ( Table 4). The stability ranking of the first four genes determined by BestKeeper was inconsistent with that of the ΔCT method (Table 4). According to RefFinder, the gene stability ranked from the most stable to the least stable across all samples was as follows: (Figure 2(i)). All pairwise variation values exhibited by geNorm were above the proposed 0.15 cut-off value (Figure 3(a)). The value at V 5/6 was lowest in all pairwise variations; five are proposed as reference genes for normalization. As described above, the correlation of NF values between NF 3 and NF 5 was calculated, and there is a good correlation between two NF measures for all biotic samples (r = 0.998) (Figure 3(b)-(D)). So, the three most stable reference genes (CypA, RPL12, and EF) may meet an accurate normalization for all sample. This result is basically consistent with that of RefFinder.

Validation of Reference Gene Selection
To assess the validity of selected reference genes, the expression profiles of PKG under four larval tissues and two photoperiod treatments as a sample were normalized using single reference gene or gene combinations recommended by geNorm (Figure 2) as follows: two most stable reference genes, respectively, the combination of two most stable reference genes, and the least stable gene. Among four various larval tissues, the results of normalizing PKG expression were similar when, respectively, using EF or RPL12 (the most stable reference gene) and the combination EF and RPL12 (the most stable two) (Figure 4(a)). However, the difference of PKG expression among four tissues was significantly different using -ACT normalizing PKG expression and the expression levels of PKG normalized using -ACT (least stable) were 7.7fold to 41.1-fold higher than those of PKG normalized using the other two reference genes or gene combinations in fat bodies and 3-fold to 9.2-fold in epidermises. PKG expression normalized by GAPDH or -ACT (the most stable reference gene), the combination of GAPDH and -ACT (the most stable most), and 18S (the least stable) (Figure 4(b)) was higher in 8L:16D than that in 16L:8D, but the difference was not statistically significant (P>0.05).

Discussion
RT-qPCR is a routine technique widely used to quantify mRNA levels of target genes. Its accuracy and reliability rest with the reference gene used as endogenous control. However, to date, there are no universal reference genes that are stably expressed in all types of samples under various treatment conditions. Therefore, the evaluation of reference genes must be conducted prior to gene expression profiling experiments. In this study, we firstly examined the validation of ten reference genes commonly applied in all kinds of organisms in M. separata with four algorithms (geNorm, NormFinder, BestKeeper, and the ΔCT method) under various biotic and abiotic conditions. The ranking orders of investigated reference genes across the four programs were varied based on various conditions (Tables 2(a), 2(b), 3, and 4). RefFinder, a web-based algorithm, integrates all data obtained from the four programs mentioned above and determines the stability of selected reference genes by calculation of Geometric Mean values, and those with the lower GM values are identified as more stable genes. RefFinder proposed that EF, -TUB, CypA, and RPL12 are the most stable reference genes for M. separata under all samples (Figure 2(i)). Elongation factor-1 alpha (EF), a GTP-binding protein that catalyzes the binding of aminoacyl-transfer RNAs to the ribosome [40], was the most stable reference gene for all biotic factors (developmental stages, tissues, and densities) proposed by RefFinder. The stability of EF in two biotic factors (developmental stages and tissues) is in accordance with reference gene analyses in Drosophila melanogaster [41], Orthoptera [42], Hymenoptera [43], and Plutella xylostella [44] and is counter to that in S. exigua [28]. Under temperature treatments, the ranking of EF in M. separata is inconsistent with that in P. xylostella [44] and H. armigera [29], in which EF was one of the least stable genes. EF was relatively stable across different treatments in M. separata; however, its expression level was low, with a mean CT of 24.49 ( Figure 1). Therefore, the expression level of target gene will decide whether EF may be used as a control gene.
CypA is a peptidyl-prolyl cis/trans isomerase originally identified as the target of the immunosuppressive drug cyclosporine A [45]. Its applications as reference gene in insects were relatively spares. In M. separata, the expression stabilities across developmental stages and tissues were higher than those in other conditions. This result was similar to the stabilities of CypA in Hippodamia convergens [46]. However, it was the least under the developmental stages and tissues in Danaus plexippus (L.) [47]. More remarkably, the expression of CypA was unstable among four densities (Figure 2(d)). Previous study has examined that M. separata larvae reared at high density show higher resistant to nucleopolyhedrovirus (NPV) than those reared at low density [48]. CypA was also reported to regulate several severe viruses in human [49]. So, further studies are needed to determine whether CypA is involved in the regulation of resistance to NPV in M. separata while controlling for larval density.
Ribosomal proteins (RPs) are a large group of proteins which, in conjunction with rRNA, make up the ribosomal subunits involved in the cellular process [50]. In mice (C57B16, 10-12 weeks old), ribosomal protein genes exhibited important tissue-dependent variation in mRNA expression and cannot be considered as true reference genes [51]. In the present study, RPL12 gene is identified as being the most stable gene under densities and photoperiod treatments and being of relative higher stable ranking in developmental stages, larval tissues, and temperature treatments ( Figure 2). Previous researched also exhibited that RP genes were stable under varied treatments in insects, as RP49 in the different developmental stages of B. mori [23], PRL10, and RPL7 in the all tissue sample sets of S. exigua [28], RPL15 under different larvae tissues and RPL32 under all tested samples in H. armigera [29], and RPL18 in developmental stages in A. obliqua [30]. Although the number of ribosomal proteins is large, no single gene has been identified to show expression stability across all biotic and abiotic treatments up to now.
18S and 28S, ribosomal RNAs, were suggested by related studies to be applicable reference genes, due to the fact that their transcript is less affected by treatments that significantly alter mRNA expression [52][53][54]. In this study, 18S and 28S exhibited the least stable expression at the different developmental stages, photoperiod treatments, and temperature treatments. This result conflicts with that in H. armigera [29], but identifies with the findings of P. xylostella [44]. These results further prove that it is necessary to validate the expression stability of reference genes.
Further, the expression of PKG was conducted in different tissues and in response to photoperiod stress for validating the reference genes. Although the results demonstrated that the expression trends in different tissues were accordant using various reference gene or gene combinations. PKG transcript was not induced by photoperiod. However, the results showed that using unstable reference genes may generate the wrong interpretation, and stable reference gene combinations may reduce bias in normalization (Figure 4(a)). Therefore, the validation of candidate reference genes should be conducted to accurately estimate target gene expression, and two or three gene combinations should be used for data normalization.

Conclusion
In the present study, ten candidate reference genes for normalizing RT-qPCR data in M. separata were selected and validated for their expression stability under various biotic and abiotic conditions. Results obtained from RefFinder exhibited that the first three most stable gene were as follows

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.