Reference Genes for Expression Analyses by qRT-PCR in Phthorimaea operculella (Lepidoptera: Gelechiidae)

Simple Summary Quantitative real-time fluorescent polymerase chain reaction (qRT-PCR) is a momentous tool for calculating the expression levels of targeted genes across various experimental conditions. The selection and evaluation of stable reference genes for qRT-PCR analysis is an essential precondition for reliable expression assessment. Phthorimaea operculella is one of the most serious Lepidopteran pests that attack potatoes around the world. In the present paper, a total of 10 commonly used reference genes, namely ACT, α-TUB, 18S, 28S, GAPDH, EF1α, RPL4, RPL13, RPL27 and SOD, were selected and validated for suitability under three treatments (developmental stages, tissues/organs and temperatures) using five methods (Ct value, geNorm, NormFinder, BestKeeper and RefFinder). These results indicated that EF1α and RPL13 were the best suitable reference genes for diverse backgrounds. The relative transcript levels of the target gene chitin synthase A gene (PoChSA) were abundantly expressed in epidermal cells, and lowly transcribed in the midgut. Our findings will be beneficial for improving the accuracy of qRT-PCR analysis for future functional analysis of the target gene expression in P. operculella. Abstract Due to a lack of effective internal references, studies on functional genes in Phthorimaea operculella, a serious Lepidopteran pest attacking potatoes worldwide, have been greatly limited. To select suitable endogenous controls, ten housekeeping genes of actin (ACT), α-tubulin (α-TUB), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1α (EF1α), 18S and 28S ribosomal RNA (18S, 28S), ribosomal protein genes RPL4, RPL13 and RPL27 and superoxide dismutase (SOD) were tested. Their expression levels were determined under three different experimental conditions (developmental stages, tissues/organs and temperatures) using qRT-PCR technology. The stability was evaluated with five methods (Ct value, geNorm, NormFinder, BestKeeper and RefFinder). The results clarified that RPL13, EF1α and RPL27 are ranked as the best reference gene combination for measuring gene expression levels among different developing stages and under various temperatures; EF1α and RPL13 are recommended to normalize the gene expression levels among diverse tissues. EF1α and RPL13 are the best reference genes in all the experimental conditions. To validate the utility of the selected reference pair, EF1α and RPL13, we estimated the tissue-biased expression level of chitin synthase A gene (PoChSA). As expected, PoChSA was abundantly expressed in ectodermally derived epidermal cells, and lowly transcribed in the midgut. These findings will lay the foundation for future research on the molecular physiology and biochemistry of P. operculella.


Introduction
Quantitative real-time fluorescent polymerase chain reaction (qRT-PCR) is a powerful tool for the quantification of nucleic acids owing to its advantages of high specificity,

Specimens among Various Tissues
Ten fully grown larvae were selected as a replicate. They were dissected and the head capsule, foregut, midgut, hindgut, fat body, hemocytes and epidermis were collected. The tissue collection was repeated three times. The tissue specimens were placed in RNAlater R (Thermo Fisher Scientific Inc., Waltham, MA, USA) and stored for several weeks at −80 • C before total RNA isolation.

Collections during Varied Temperature Incubation
The final instar larvae were transferred into three temperatures (4 • C, 26 • C and 35 • C). Ten larvae as a replicate were collected after 2, 6 and 12 h. A total of nine treatments were set. The collection was repeated three times and stored for several weeks at −80 • C before total RNA isolation.

Samples for the Expression Analysis of PoChSA
The ultimate instar larvae were dissected and the head capsule, foregut, midgut, hind gut and epidermis were collected. A total of 10 individuals were dissected for each replicate. The tissue collection was repeated three times. The tissue specimens were placed in RNAlater R (Thermo Fisher Scientific Inc., Waltham, MA, USA) and stored for several weeks at −80 • C before total RNA isolation.

Selection and Authentication of Candidate HKGs
Ten HKG sequences (actin, ACT; α-tubulin, α-TUB; glyceraldehyde-3-phosphate dehydrogenase, GAPDH; elongation factor 1α, EF1α; 18S and 28S ribosomal RNA, 18S and 28S; ribosomal proteins RPL4, RPL13 and RPL27; superoxide dismutase, SOD) were selected. The accession numbers of these genes are listed in Table 1. Reverse transcription PCR (RT-PCR) was performed to authenticate the HKGs using the primers listed in Table 1. The amplified products were separated by electrophoresis on 1.5% agarose gel and purified using the Wizard ® PCR Preps DNA Purification System (Promega, Madison, WI, USA). Purified DNA was ligated into the pGEM ® -T easy vector (Promega, Madison, WI, USA) and several independent subclones were sequenced from

Quantitative Real-Time PCR (qRT-PCR)
The qRT-PCR primers were designed using Beacon Designer 7 (Premier Biosoft International, Palo Alto, Santa Clara, CA, USA), and are given in Table 2. The qRT-PCR reactions were performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) and QuantStudio™ 7 Pro Real-Time PCR System (Applied Biosystems, Life Technologies, Carlsbad, CA, USA) according to the manufacturer's protocol. The reaction mixture consisted of 1 µL of cDNA template, 10 µL of 2× ChamQ Universal SYBR qPCR Master Mix, 0.4 µL of forward primer (10 µM), 0.4 µL of reverse primer (10 µM) 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 dimers or contamination in the reactions, respectively. The qRT-PCR protocol included an initial step of 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s and then annealed at 60 • C for 34 s, followed by one cycle of 95 • C for 15 s, 60 • C for 60 s and 95 • C for 1 s. PCR amplicons were subjected to melting curve analysis. The specificity of the qRT-PCR reactions was monitored by melting curve analysis using QuantStudio™ Design & Analysis Software (version 1.5.0) and gel electrophoresis. Amplification efficiencies were determined by a 10-fold dilution series of template. All experiments were repeated in triplicate.

Evaluation of Reference Gene Selection
ChSA of P. operculella was used to evaluate the stability of candidate reference genes. The primer sequence of the target gene was as follows: forward (5 -GCCTGGAGTTCACAGTCAGA-3 ) and reverse (5 -GCCGGTCTTTCTTAAGTTGC-3 ). The average relative levels of PoChSA in different tissues were computed based on 2 −∆∆CT method and from three replicates. We used SPSS for Windows (Chicago, IL, USA) for statistical analyses. The averages (±SE) were submitted to analysis of variance with the Tukey-Kramer test.

Data Processing
The raw Ct values were obtained using the QuantStudio™ Design & Analysis Software (version 1.5.0). The algorithms, including geNorm [26], BestKeeper [27] and Normfinder [28], were used to analyze the stability of selected HKGs, strictly following the manuals of the algorithms. Finally, the comprehensive ranking of each condition was obtained according to RefFinder [29,30].

Selection of Candidate HKGs
We selected ten HKG genes and designated them as ACT, α-TUB, 18S, 28S, GAPDH, EF1α, RPL4, RPL13, RPL27 and SOD. The resultant sequences were submitted to GenBank; the accession numbers were listed in Table 1. The correctness of the ten HKGs was proven by RT-PCR.
The products from qRT-PCR were confirmed by sequencing. The primer specificity for qRT-PCR was verified by melting curve analysis. All the primer pairs amplified a single PCR product with the expected sizes and sequences, showed a slope less than −3.0 and exhibited regression coefficient (R2) and efficacy values ranging from 0.991-0.999 and 93.80-103.89% (Table 2). These data indicate that the amplification efficiencies of the primers reached the standard requirements of conventional qRT-PCR [5].

Expression Variations of the Ten HKGs
The specimens were collected from four developmental stages (young and old larvae, pupae and adults), seven larval tissues (head capsule, foregut, midgut, hindgut, fat body, hemocytes and epidermis) and three temperature treatments (4 • C, 26 • C and 35 • C). Using the products obtained by qRT-PCR for agarose gel electrophoresis, we found that all ten genes had single amplicons of expected size (data not shown). Therefore, these ten genes were expressed during different developmental stages, among different larval tissues and under different temperatures.
The overall threshold cycle (Ct) values under different experimental conditions are shown in Figure 1 and Table S1. Across developing stages, EF1α and RPL13 had the smaller gene expression variation, whereas ACT and 18S had the higher expression difference ( Figure 1A). Among various tissues, except for GAPDH and SOD, the expression fluctuations were small in selected HKGs ( Figure 1B). Under different temperatures, the expression fluctuations were small in selected HKGs except for SOD and GAPDH ( Figure 1C). A combination of these results revealed that the variations in RPL13, EF1α, RPL27 and α-TUB were smaller, whereas the ranges in ACT, GAPDH, 18S and SOD were larger ( Figure 1D).

Expression Stability of the Ten HKGs during Developmental Stages
The geNorm algorithm evaluates the candidate reference genes based on their expression stability values (M-values) and pairwise variations (Vn/Vn+1). The expression stability values revealed that EF1α, RPL13 and 28S were the better reference genes during developing, with M-values below 0.5. The values of other genes were below 1, except for ACT, and their stability values were similar ( Figure 2A, Table 3). The pairwise variation analysis showed that the V3/4 value was near 0.15; indicating three different reference genes are needed for gene expression analysis during development ( Figure 2B).      Expression stability of ten house-keeping genes among various tissues in Phthorimaea operculella. Head capsule, foregut, midgut, hindgut, fat body, hemocytes and epidermis were dissected from the fourth instar larvae. The expression stability rankings were determined by geNorm, NormFinder and BestKeeper.

Expression Stability of the Ten HKGs during Developmental Stages
The geNorm algorithm evaluates the candidate reference genes based on their expression stability values (M-values) and pairwise variations (Vn/Vn+1). The expression stability values revealed that EF1α, RPL13 and 28S were the better reference genes during

Expression Stability of the Ten HKGs during Developmental Stages
The geNorm algorithm evaluates the candidate reference genes based on their expression stability values (M-values) and pairwise variations (Vn/Vn+1). The expression stability values revealed that EF1α, RPL13 and 28S were the better reference genes during developing, with M-values below 0.5. The values of other genes were below 1, except for ACT, and their stability values were similar (Figure 2A, Table 3). The pairwise variation analysis showed that the V3/4 value was near 0.15; indicating three different reference genes are needed for gene expression analysis during development ( Figure 2B).
According to the NormFinder, those genes with low stability values, based on intra- Expression stability of ten house-keeping genes in different samples of Phthorimaea operculella. The stability of the reference genes calculated by the Geomean method of RefFinder. A lower Geomean of ranking value denotes more stable expression. According to the NormFinder, those genes with low stability values, based on intraand inter-group expression variations, are considered to be the most stable reference genes. Across different development stages, the stable genes were RPL13, EF1α and RPL27, with the p value less than 1.0. The most unstable gene was ACT, with the p value of 3.8 ( Figure 2C, Table 3).
Based on the BestKeeper analysis, the stable orders of selected HKGs were RPL13, 28S, EF1α, α-TUB, RPL27, SOD, RPL4, GAPDH, 18S and ACT, from the most stable to the least. The last two genes, 18S and ACT, had Cp values of more than 1 ( Figure 2D, Table 3), indicating that they should be excluded as reference genes for qRT-PCR to test the expression level of the target gene.
The online tool RefFinder combines the three methods above to compare and rank the tested reference genes [29]. It ranks the selected HKGs in the following order from the most to least stable: Figure 5A). Therefore, RPL13, EF1α and RPL27 are ranked as the best reference gene combination for measuring target genes among different developing stages.

Expression Stability of the Ten HKGs among Different Tissues
Among the three tissues, the stability of the selected HKGs were EF1α = RPL13 > RPL4 > 28S > RPL27 > α-TUB > 18S > ACT > GAPDH > SOD; based on the geNorm algorithm, the M-values of EF1α, RPL4 and RPL13 were below 0.4 ( Figure 3A, Table 3). The pairwise variation analysis displayed that the V2/3 to V8/9 values were below 0.15, suggesting two reference genes are enough for gene expression determination within various tissues ( Figure 3B).
The BestKeeper data uncovered that EF1α, RPL13, 28S, RPL4 were the most stable because they showed Cp values of 0.200, 0.279, 0.328 and 0.373, respectively. The Cp values of α-TUB, RPL27, 18S, ACT and GAPDH were less than 1.0., and the Cp value of SOD was more than 2.0 ( Figure 3D, Table 3).
The RefFinder showed a comprehensive ranking order from the most to the least stable: Figure 5B). Thus, the two HKGs (EF1α and RPL13) are recommended to be used to test the target gene expression levels among various tissues.

Stability of the Ten HKGs under Different Temperatures
The geNorm algorithm results showed that the comprehensive reference gene rankings from the best to the least stable were EF1α, RPL13, RPL27, ACT, α-TUB, RPL4, 28S, 18S, GAPDH and SOD ( Figure 4A, Table 3). Except for SOD, the other genes in the selected HKGs showed values below 1, indicating their stabilities were similar. Moreover, the pairwise variation analysis showed that the V3/4 value was below 0.15, indicating three different reference genes are needed for gene expression analysis under different temperatures ( Figure 4B).
By the NormFinder analysis, the stable orders of the selected HKGs from the most stable to the least were RPL13, α-TUB, EF1α, RPL27, RPL4, ACT, 28S, 18S, GAPDH and SOD. Again, the p values revealed by the NormFinder analysis indicated that RPL13, α-TUB, EF1α, RPL27 were smaller, demonstrating that the genes have similar stability ( Figure 4C, Table 3).
The BestKeeper data unveiled that the steady orders were 28S, α-TUB, 18S, RPL13, EF1α, RPL27 and RPL4 ( Figure 4D, Table 3). Since the Cp values of ACT, GAPDH and SOD were more than 1, they cannot be used as reference genes for qRT-PCR to test the expression level of the target gene. The other genes showed values below 1, indicating their stability values were similar ( Figure 4D, Table 3).
According to the RefFinder results, the stability rankings were as follows: Figure 5C).
When the three different conditions were combined together, the RefFinder results indicated that the stability rankings from the most to the least were RPL13, EF1α, RPL27, α-TUB, 28S, RPL4, GAPDH, 18S, ACT and SOD ( Figure 5D). Thus, the two HKGs (EF1α and RPL13) can be selected as reference genes for measuring the target gene expression levels among diverse backgrounds.

Validation of the Selected Reference Genes after Gene Expression
To demonstrate the utility of EF1α and RPL13 in accurate gene expression analysis, the expressions of chitin synthase A gene (PoChSA) in the head capsules, epidermis, foregut, midgut and hindgut were calculated after normalization with a combination of EF1α and RPL13. The highest accumulated mRNA level of PoChSA was found in the head capsule and epidermis, followed by those in the foregut and hindgut; the lowest level was detected in the midgut (Figure 6). and RPL13. The highest accumulated mRNA level of PoChSA was found in the head capsule and epidermis, followed by those in the foregut and hindgut; the lowest level was detected in the midgut ( Figure 6).

Figure 6.
Tissue expression of the chitin synthase A gene (PoChSA) in Phthorimaea operculella. The head capsule (Head), foregut (FG), midgut (MG), hindgut (HG) and epidermis (EP) were dissected from the fourth-instar larvae. For each sample, 3 independent pools of 20-30 individuals were measured in technical triplicate using qRT-PCR. The values were calculated using the 2 −ΔΔCT method, using the selected reference genes EF1α and RPL13. The relative transcripts are the ratios of copy numbers in different developing stages relative to the head capsule, which is set as 1. The columns represent averages, with vertical lines indicating SE. Different letters indicate significant difference at p value < 0.05 using analysis of variance with the Tukey-Kramer test.

Discussion
In the present paper, we investigated the expression stability of ten HKGs in P. operculella. Out of the ten HKGs, ACT, RPL, TUB, GAPDH, 18S and EF1α are the top 10 most frequently used reference genes [25].
It is widely accepted that moderately expressed HKGs should be chosen as potential reference genes because genes with extremely high or low expression levels are less reliable [31]. According to the Ct value obtained in the present paper, ACT and SOD are the less expressed and GAPDH and 18S are the most expressed. Even though the results obtained using the BestKeeper, geNorm and NormFinder algorithms were not completely consistent, the data still revealed that the mRNA levels of ACT and SOD are changed dramatically throughout the developing stages among tissues and under different temperatures, respectively. It appears that the four genes should be excluded as reference genes for qRT-PCR.
Actin (ACT) plays an important role in cell secretion, motility cytoplasm flow and experimental cytoskeleton maintenance and is abundantly expressed in most cell types. Even though ACT is used for qRT-PCR studies when measuring the expression of target genes in P. operculella [3,4], it was verified to be one of the most unstable genes in the present paper. The transcript level of ACT is also less stable in several Coleopteran insect species, such as Phaedon brassicae, Henosepilachna vigintioctomaculata, Leptinotarsa decemlineata, Coleomegilla 11 aculate, Coccinella septempunctata and Hippodamia convergens [32][33][34][35][36][37], although ACT is one of the most stable reference genes across several developmental stages in Orthopteran (Schistocerca gregaria and Chortoicetes terminifera), Hemipteran (Diu- Figure 6. Tissue expression of the chitin synthase A gene (PoChSA) in Phthorimaea operculella. The head capsule (Head), foregut (FG), midgut (MG), hindgut (HG) and epidermis (EP) were dissected from the fourth-instar larvae. For each sample, 3 independent pools of 20-30 individuals were measured in technical triplicate using qRT-PCR. The values were calculated using the 2 −∆∆CT method, using the selected reference genes EF1α and RPL13. The relative transcripts are the ratios of copy numbers in different developing stages relative to the head capsule, which is set as 1. The columns represent averages, with vertical lines indicating SE. Different letters indicate significant difference at p value < 0.05 using analysis of variance with the Tukey-Kramer test.

Discussion
In the present paper, we investigated the expression stability of ten HKGs in P. operculella. Out of the ten HKGs, ACT, RPL, TUB, GAPDH, 18S and EF1α are the top 10 most frequently used reference genes [25].
It is widely accepted that moderately expressed HKGs should be chosen as potential reference genes because genes with extremely high or low expression levels are less reliable [31]. According to the Ct value obtained in the present paper, ACT and SOD are the less expressed and GAPDH and 18S are the most expressed. Even though the results obtained using the BestKeeper, geNorm and NormFinder algorithms were not completely consistent, the data still revealed that the mRNA levels of ACT and SOD are changed dramatically throughout the developing stages among tissues and under different temperatures, respectively. It appears that the four genes should be excluded as reference genes for qRT-PCR.
18S ribosomal RNA is a part of the ribosomal RNA (rRNA), which accounts for more than 80% of the total RNA pool [38], whereas mRNA accounts for only 3 to 5%. This is consistent with our data that 18S is the most expressed in P. operculella (this study). Therefore, the use of rRNAs as reference genes may mask subtle changes in target mRNAs [39]. Moreover, 18S shows a large variation in different development stages in Myzus persicae [40].
Superoxide dismutase (SOD) is known as an antioxidative stress protein by scavenging the superoxide radicals, used to defend against reactive oxygen species (ROS) damage caused by a variety of unfavorable environmental stressors in some insect species [41][42][43]. In this study, the SOD gene was verified to be one of the most unstable genes under three different conditions.
Similarly, the instability of GAPDH expression has been documented in Colaphellus bowringi [44], D. caesalis [22], Mythimna separata [45], Ophraella communa [46] and P. brassicae [37]. GAPDH functions as a glycolytic enzyme involved in glycolysis and is associated with cell proliferation under adverse conditions where its catalytic activity is impaired [47]. It is presumed that any perturbation toward energy metabolism or cell proliferation would have a potential impact on GAPDH expression. Considering these issues, it is inappropriate to adopt GAPDH as a reference gene. Therefore, we focus on five genes, i.e., α-TUB, EF1α, RPL4, RPL13 and RPL27, in P. operculella for the selection of a suitable reference gene combination.
It has been suggested that multiple reference genes should be used in order to avoid biased normalization [48]. Additionally, from the present study, we recommended two reference genes, EF1α and RPL13, to normalize the gene expression levels among diverse conditions in P. operculella. Consistent with our results, the conserved nuclear gene elongation factor 1 alpha (EF1α) plays an important role in translation by catalyzing the GTP-dependent binding of aminoacyl-tRNA to the acceptor site of the ribosome. EF1α is evaluated as the most stable gene under diverse conditions in D. melanogaster [15], C. terminifera [13], Bombus terrestris and Bombus lucorum [49], Frankliniella occidentalis [50] and Helicoverpa armigera [39].
To sum up, in this study, the genes RPL13, EF1α and RPL27 are indicated to be ranked as the best reference gene combination for measuring gene expression levels among different developing stages and under various temperatures, while EF1α and RPL13 are recommended to normalize gene expression levels among diverse tissues. EF1α and RPL13 are the best reference genes in all the experimental conditions in P. operculella (Lepidoptera: Gelechiidae). Interestingly, in another Lepidoptera insect Spodoptera frugiperda (Noctuidae), based on the online program RefFinder, SOD, RPL10 and RPS24 were reported to be the most stable reference genes for different developmental stages, while α-TUB, RPL10 and ATP were for various tissues, AK, RPL10 and 18S for mating status, 18S and AK for hormone treatment, 18S, RPL10 and SOD for diets treatment, and RPL10, 18S and RPS24 for temperature treatment [67]. The results verified that the expression stability of the reference genes varied under different treatments. Similarly, the ribosomal protein genes are also the most stable reference genes selected under almost all the experimental conditions. In addition, the difference in housekeeping genes under similar treatments may be related to the phylogenetic relationship and feeding habits of the two lepidoptera insects.
In order to demonstrate the utility of EF1α and RPL13 in accurate gene expression analysis in P. operculella, we evaluated the relative gene expression level of PoChSA in the head capsules, epidermis, foregut, midgut and hindgut. Our results showed that PoChSA was abundantly expressed in the head capsule and epidermis, moderately transcribed in the foregut and hindgut and lowly expressed in the midgut. Our expression data are consistent with the fact that ChSA encodes an enzyme that catalyzes the biosynthesis of chitin in the ectodermally derived epidermal cells forming epidermis, trachea, foregut and hindgut [3,[68][69][70][71]. Thus, the tissue-biased expression pattern of PoChSA demonstrates that EF1α and RPL13 can be used as endogenous controls to assess gene expression in P. operculella.

Conclusions
Our findings recommend EF1α and RPL13 as the optimal reference gene set under three different experimental conditions. EF1α and RPL13 combinations can be proposed as reference genes for measuring the target gene expression levels among diverse backgrounds in P. operculella. To date, this is the first study to screen out candidate reference genes for gene expression analysis in P. operculella. The results lay a foundation for molecular research. Nevertheless, the application of these loci as reference genes under other physiological or experimental conditions remains to be determined.