Assessment of suitable reference genes for qRT-PCR analysis in Adelphocoris suturalis

LUO Jing , MA Chao LI Zhe ZHU Bang-qin, ZHANG Jiang LEI Chao-liang, JIN Shuang-xiaJ. Joe Hull, CHEN Li-zhen

1 National Key Laboratory of Crop Genetic Improvement and National Centre of Plant Gene Research, Huazhong Agricultural University, Wuhan 430070, P.R.China

2 Hubei Collaborative Innovation Center for Green Transformation of Bio-Resources, College of Life Science, Hubei University,Wuhan 430062, P.R.China

3 Hubei Insect Resources Utilization and Sustainable Pest Management Key Laboratory, College of Plant Science and Technology,Huazhong Agricultural University, Wuhan 430070, P.R.China

4 U.S. Arid Land Agricultural Research Center, Agricultural Research Service, U.S. Department of Agriculture, Maricopa 85138,USA

Abstract Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is the most commonly-used tool for measurement of gene expression, but its accuracy and reliability depend on appropriate data normalization with the use of one or more stable reference genes. Adelphocoris suturalis is one of the most destructive pests of cotton, but until recently knowledge of its underlying molecular physiology had been hindered by a lack of molecular resources. To facilitate research on this pest, we evaluated 12 common housekeeping genes studied in insects (GAPDH, ACT, βACT, TBP, SDH, βTUB, EF1γ,EF1α, EF1δ, RPL32, RPS15, and RPL27) for their expression stability in A. suturalis when subjected to various experimental treatments, including three biotic (developmental stage and sex, tissue type, and metathoracic scent gland for varying developmental stages and sexes) and one abiotic (RNA interference injection) conditions. Four dedicated algorithms (ΔCt method, geNorm, BestKeeper and NormFinder) were used to analyze gene expression stability. In addition, RefFinder provided an overall ranking of the stability/suitability of these candidates. This study is the flrst to provide a comprehensive list of suitable reference genes for gene expression analyses in A. suturalis, which can serve to facilitate transcript expression study of related biological processes in this and related species.

Keywords: Adelphocoris suturalis, reference gene, qRT-PCR, normalization, expression stability

1. Introduction

Adelphocoris suturalis(Hemiptera: Miridae) is a highly polyphagous pest that can attack a broad range of cultivated crops including cotton, pastures, vegetables, and fruit trees(Jianget al.2015). This plant bug was originally a secondary pest of cotton, but has since become a signiflcant problem for cotton growing regions in China due to the reduction in broad-spectrum insecticides that followed the widespread adoption of transgenicBacillus thuringiensis(Bt) cotton(Luet al.2008a, 2010; Liet al.2010). Currently, control of these plant bugs is largely dependent on the use of broadspectrum chemical insecticides (Lu and Wu 2008). Although such treatments can be effective, they often threaten human health, adversely affect the environment, and can give rise to resistant populations. Hence, the development of less hazardous, environmentally sound, and sustainable pest management strategies is needed. However, a deeper understanding of the biology of these species is necessary for developing novel alternatives to broad-spectrum chemical control approaches.

In recent years, the ecology and physiology ofA. suturalishas been expansively studied (Luet al.2008b; Lu and Wu 2008; Zhanget al.2014; Jianget al.2015; Zhanget al.2016), but the molecular mechanisms of the biology have largely remained unknown. With the advent of next generation sequencing approaches, there is now an unparalleled opportunity to investigate the genetic basis of its biology and physiology. For this investigation, both the silencing of gene expression to assess gene function as well as quantiflcation of gene expression are required.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) is an indispensable tool for measurement of gene expression (Bustinet al.2010). Older, traditional gene expression measurements such as competitive reverse transcription polymerase chain reaction (RT-PCR),Northern blot analysis, orin situhybridization lacks its high sensitivity, speciflcity, reproducibility, and high-throughput convenience (Ginzinger 2002; Wong and Medrano 2005;Espyet al.2006; Bustin 2010). However, the accuracy of qRT-PCR is limited by batch to batch variation in RNA extraction and variable efflciency of cDNA synthesis and of PCR reaction (Bustinet al.2005; Huggettet al.2005).To limit variability, a common technique in qRT-PCR is normalizing target gene expression data to reference genes(Vandesompeleet al.2002). Consequently, one of the most important considerations in a valid qRT-PCR analysis is the appropriateness of the reference gene(s).

Housekeeping genes, which are expressed constitutively and are essential for survival in all cells, are commonly used as reference genes, under the assumption that their transcript levels will remain constant regardless of experimental treatment and/or physiological condition (Thellinet al.1999;Butteet al.2001). However, these assumptions may not be valid in practice. Numerous reports have demonstrated that some of the most common reference genes (housekeeping genes) undergo signiflcant regulation in response to diverse experimental conditions (Leeet al.2002; Yuanet al.2014; Bansalet al.2016), indicating that these genes are inappropriate for normalization purposes. Indeed, there is no stably expressed reference gene suitable for all cells and experimental conditions. If the reference gene is not selected properly, it will cause erroneous interpretation of the qRT-PCR data (Leeet al.2002; Vandesompeleet al.2002; Radoni?et al.2004; Huggettet al.2005; Nolanet al.2006; Yuanet al.2014). Therefore, it is essential to select and validate a suite of best-suited qRT-PCR reference gene(s) prior to quantifying transcript levels. However,the expression stability of reference genes inA. suturalis,or in any plant bug, has not been systematically assessed to date. Hence, the identiflcation and validation of suitable endogenous reference genes inA. suturalisis urgently needed.

The overall goal of this study was to evaluate and select appropriate reference gene(s) for accurate quantiflcation of mRNA transcripts inA. suturalis. To achieve this goal, 12 housekeeping genes frequently used in gene expression study of other sap-sucking insects were selected as candidate reference genes:glyceraldehyde-3-phosphate dehydrogenase(GAPDH),actin(ACT),β-actin(βACT),TATA-box binding protein(TBP),succinate dehydrogenase(SDH),β-tubulin(βTUB),elongation factor-1γ(EF1γ),elongation factor-1α(EF1α),elongation factor-1δ(EF1δ),ribosomal protein L32(RPL32),ribosomal protein S15(RPS15), andribosomal protein L27(RPL27) (Liet al.2013; Bansalet al.2016; Ibanez and Tamborindeguy 2016; Koramutlaet al.2016; Zhanget al. 2018). These housekeeping genes were chosen from different functional classes and gene families to avoid the effect of co-regulation.We evaluated the expression stability of these candidate genes under varying abiotic (RNA interferenceviadsRNA injection) and biotic (developmental stage and sex, multiple tissue types, and a narrow focus on metathoracic scent glands (MTGs) from different developmental stages and sexes) conditions using the four statistical algorithms geNorm (Vandesompeleet al.2002), NormFinder (Andersenet al.2004), BestKeeper (Pfafflet al.2004), and the ΔCt method (Silveret al.2006). In addition, RefFinder, a comprehensive platform integrating the above-mentioned algorithms, provided an overall ranking of the stability/suitability of these candidates (Yanget al.2015a). Our data provide the flrst comprehensive assessment of reference genes inA. suturalis, and will beneflt gene expression studies in this plant bug species.

2. Materials and methods

2.1. Insect rearing

A. suturaliswas collected from a Bt cotton fleld located in Wuhan (Hubei Province, China) in July 2015. Nymphs and adults were reared in plastic cages (22.5 cm×15 cm×11 cm) and fed a diet of green beans and 5% sugar solution.The green beans were also used to collect eggs (Luet al.2008b). The newly emerged adults were separated by sex every day and reared in a separate plastic cage at a density of 30 adults per cage (Lu and Wu 2008). All insects were maintained for their entire life cycle at (26±2)°C, (75±5)%relative humidity, and a photoperiod cycle of 16 h L:8 h D.

2.2. Experimental treatments

Developmental stage and sexThe different developmental stages and sexes ofA. suturalisincluded eggs, second instar nymphs, flfth instar nymphs, sexually immature (1-d-old)male and female adults, and sexually mature (8-d-old) male and female adults. A total of 100 eggs, 30 second instar nymphs, 6 flfth instar nymphs, 3 male adults (1-d-old), 3 female adults (1-d-old), 3 male adults (8-d-old) and 3 female adults (8-d-old) were collected, then immediately frozen in liquid nitrogen and stored at -80°C until use. The samples were collected in three biological replicates.

TissueFive body regions, including head, MTG, gut, ovary and fat body, were dissected from 8-d-old female adults ofA. suturalis. Three independent biological replicates were collected and each tissue group was derived from a minimum of 25 insects. All the samples were handled and stored as described above.

MTG at different developmental stages and sexesThe MTGs of 1-d-old male, 1-d-old female, 8-d-old male, and 8-d-old female adults ofA. suturaliswere dissected. Three independent biological replicates were collected with each group corresponding to a minimum of 30 insects. All the samples were handled and stored as described above.

dsRNA injectionFor RNAi treatments, 1-d-old female adults ofA. suturaliswere microinjected with 100 nL dsRNA (10 μg μL–1) againstfatty acyl-CoA reductase(FAR), encoded an NAD(P)H-dependent oxidoreductase that catalyzes the reduction of fatty acyl-CoA precursors into fatty alcohols, using a micro-injector (World Precision Instruments, Sarasota, FL, USA). Control injections consisted of H2O anddsGFP(dsRNA ofgreen fluorescent proteingene). The dsRNA was synthetized using the corresponding primers (Appendix A) as described (Liuet al.2016). At 3 days post-injection, 3 individuals from each treatment were collected with 3 independent biological replicates performed. All the samples were handled and stored as described above.

2.3. Total RNA isolation and cDNA synthesis

Total RNA was extracted using RNAiso Plus reagent(TaKaRa, Kyoto, Japan) according to the manufacturer’s protocol. The RNA integrity was conflrmed by 1.5% agarose gels electrophoresis and quantifled on a Nano-Drop 2000(Thermo Scientiflc, Wilmington, DE, USA). Total RNA(1 μg) was depleted of residual genomic DNA and then reverse transcribed using the PrimeScript?RT Reagent Kit with gDNA Eraser (TaKaRa, Kyoto, Japan) according to the manufacturer’s protocol. The synthesized cDNA was stored at –20°C until use. For qRT-PCR analysis, each cDNA sample was diluted 20 times with nuclease-free water.

2.4. Candidate reference genes and primer design

UsingA. suturalistranscriptomic data (Luoet al.2014),sequences corresponding to 12 candidate reference genes commonly used in qRT-PCR analyses in other insect species were selected (Appendix B). All of the candidate reference genes were PCR amplifled fromA. suturaliscDNA using the corresponding primers (Appendix C), sub-cloned using the pEASY-T1 Simple Cloning Kit (TransGen, Beijing, China)and sequenced. The corresponding gene sequences were deposited in GenBank (the accession numbers are listed in Table 1). After conflrmation of candidate reference genes, primers for the subsequent qRT-PCR analyses were designed using an online tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The PCR ampliflcation efflciency and primer speciflcity were assessed using standard curves,melt curve analyses, and 2% agarose-gel electrophoresis.

2.5. qRT-PCR

All of the qRT-PCR reactions were performed using a Bio-Rad iQ2 Real-time PCR Detection System (Bio-Rad, Hercules, CA, USA) following the MIQE (Minimum Information for publication of Quantitative real time PCR Experiments) guidelines (Bustinet al.2010). Reactions were carried out in a flnal volume of 10 μL containing 2 μL diluted cDNA template, 5 μL 2×SYBR?PremixExTaq? II(TaKaRa, Kyoto, Japan) and 400 nmol L–1of each genespeciflc primer (Table 1). The reaction cocktails were set up in 96-well format Microseal PCR plates (Bio-Rad Hercules,CA, USA) in triplicate. The PCR program consisted of 95°C for 30 s, followed by 40 cycles of dissociation at 95°C for 5 s, and then annealing and extension at 62°C for 30 s. For melt curve analysis, continuous fluorescent measurements were made as the temperature was ramped up from 55 to 95°C in increments of 0.5°C every 6 s. A serial 5-fold dilution of cDNA template was used to generate standard curves and the gene speciflc PCR efflciency (E) of each gene were calculated according to previously described methods (Yanget al.2015a; Koramutlaet al.2016). Three biological replicates were performed for individual treatment.

2.6. Stability analysis of candidate reference genes

The Bio-Rad iQ5 Optical System software (ver. 2.1.94.617)(Bio-Rad, Hercules, CA, USA) was used to analyze the qRTPCR data. Three biological replicates were used to calculate the average cycle threshold (Ct) values. The stabilities of the 12 candidate reference genes were evaluated by geNorm (Vandesompeleet al.2002), NormFinder (Andersenet al.2004), BestKeeper (Pfafflet al.2004), and the ΔCt method (Silveret al.2006). GeNorm calculates an average expression stability value (M) in which lower values indicate more stable expression or lower variation (Vandesompeleet al.2002). M is calculated by a geometric averaging of the mean pairwise variation of a candidate reference gene to all the other candidate reference genes. An M value less than 1.5 is recommended to identify stably expressed gene. NormFinder determines the expression stability by considering intra- and inter-group variations for candidate reference genes (Andersenet al.2004). NormFinder provides the stability value (SV) for each candidate reference gene. Genes with a lower SV are considered to be more stably expressed and are ideal to select as reference gene for that particular experimental conditions (Andersenet al.2004). The BestKeeper program determines the stability of a candidate reference gene based on the standard deviations (SD) of the Ct values. SD values below 1 are recommended for stably expressed genes, and the lower the SD, the better the gene is as a reference (Pfafflet al.2004). In the ΔCt method, rank order is determined based on pair-wise comparisons of gene-sets using mean ΔCt values within a particular treatment. Therefore, the average standard deviation of each gene-set is inversely proportional to the gene-expression stability (Silveret al.2006). Finally,RefFinder, a comprehensive software platform integrating all four algorithms, provides an overall ranking of the stability/suitability of the candidates (Yanget al.2015a).Furthermore, geNorm performs a stepwise calculation of the pairwise variation (Vn/Vn+1) between sequential normalization factors (NFnand NFn+1) to determine the optimal number of reference genes required for accurate normalization. A threshold value below 0.15 suggests no additional reference genes are necessary for normalization (Vandesompeleet al.2002). BestKeeper and RefFinder use raw Ct values,whereas geNorm and Normflnder use expression values calculated as 2–ΔΔCT.

2.7. Validation of reference genes

Validation of the selected reference genes was performed using various tissues. TheA. suturalisvacuolar-type H+-ATPase(V-ATPase) (accession no. MF102282) gene and thefatty acid synthase(FAS) (accession no. MG520370)gene were selected as target genes for stability validation(V-ATPase, F: 5′-ACCCTCCATCAGCGTCCCAT-3′,R: 5′-AGGCGCCAAAGGAGTATCGAC-3′;FAS,F: 5′-ACTGGGG CGAATGTGGATGGTTAC-3′, R:5′-GGTCTCCTACCTTG GTTCCTGTTC-3′). The relative expression level ofV-ATPaseandFASwere normalized using the best reference gene pair (RPL32/RPS15identifled by geNorm), the single best reference gene (RPS15determined by RefFinder), and the least stable reference gene (βACTdetermined by all flve algorithms), respectively.The qRT-PCR reactions were carried out as described above and qRT-PCR data were analyzedviathe 2–ΔΔCTmethod(Schmittgen and Livak 2008). One-way ANOVA followed by Tukey’s HSD Multiple Comparison was used to determine statistical signiflcance. Statistical differences are shown as different letters.

3. Results

3.1. Performance of qRT-PCR primers

We conflrmed speciflc PCR ampliflcation by testing the primer speciflcity for each candidate reference gene with RT-PCR. PCR ampliflcations for each primer pair showed single bands of the anticipated sizes on the 2.0% agarose gel (Fig. 1-A). All amplicons were sequenced and conflrmed to exhibit 99–100% identity with the corresponding transcriptomic sequences. Melt curve analyses revealed single peaks for each primer pair, suggesting the absence of non-speciflc ampliflcation (Fig. 1-B). The PCR efflciency (E)and correlation coefflcient (R2) for each standard curve are shown in Table 1. The PCR efflciencies for all tested primer pairs varied between 91.0 and 103.4%, with associatedR2values of 0.996–1.000.

Fig. 1 Ampliflcation speciflcity of primers in RT-PCR and qRT-PCR. A, single amplicon of the expected size for each gene was visualized on a 2% agarose gel. M, marker. B, single peaks in melt curve analysis.

3.2. Expression profiles of candidate reference genes

To provide an overall representation of primer variability under varying experimental conditions, the expression proflles of the candidate reference genes were examined(Fig. 2). Ct values of the 12 candidate reference genes under the four experimental conditions spanned a range of 14.53–29.27 cycles.ACTandTBPhad average Ct values＞24 cycles, while the average Ct value of the other candidate reference genes (GAPDH,βACT,SDH,βTUB,EF1γ,EF1α,EF1δ,RPL32,RPS15andRPL27) ranged between 16–24 cycles (Fig. 2-A–D). Furthermore, we found that experimental treatment influenced the degree of variability in candidate reference gene expression. For example,βTUBvaried less (～1 cycle) between samples before and after dsRNA injection than across tissue types (＞6 cycles). When considering all the experimental conditions, we found that a limited number of the candidate genes (e.g.,EF1γ,EF1α,EF1δ,RPS15andRPL27) were relatively stable (＜4 cycle difference), whereas the others exhibited greater variation in expression (e.g.,βACTat ～11 cycles). Although variation was observed in all treatments, it was less pronounced following dsRNA injection (Fig. 2-D).

3.3. Stability of candidate reference genes

Performance of the 12 candidate genes was assessed in four experimental sets, including different developmental stages and sexes, across multiple tissues, in MTGs from different developmental stages and sexes, and after dsRNA injection. To identify the most stable reference gene(s) for these different experimental conditions, the expression stabilities were evaluated using the ΔCt method,BestKeeper, NormFinder, and geNorm. The overall stability ranking was obtained by RefFinder.

Fig. 2 Box-and-whisker plots of expression proflles of candidate reference genes under four experimental conditions. A, different developmental stages and sexes. B, different tissues. C, metathoracic scent glands (MTGs) from different developmental stages and sexes. D, dsRNA injection. The expression levels of candidate reference genes are shown as Ct values. Each data point represents the Ct values of each biological replicate in each treatment. The median is represented by the line in the box. The interquartile range is bordered by the upper and lower edges, which indicate the 75th and 25th percentiles, respectively. The whisker caps indicate the minimum and maxium values.

For different developmental stages and sexes, the ΔCt method and NormFinder indicated thatEF1δ,EF1γandRPS15were the most stable, whereasEF1αandβACTexhibited the greatest variation (Table 2). Based on BestKeeper,EF1γandGAPDHwere the most stable reference genes (Table 2). Similarly, GeNorm calculated the lowest M value for theEF1γ/GAPDHpair (0.051), suggesting that they are the most stable transcripts. The M value of all candidate reference genes exhibited little variation and remained ＞0.3 (Table 2).RefFinder ranked the genes from most to least stable as:EF1δ＞EF1γ＞RPS15＞GAPDH＞βTUB＞SDH＞RPL27＞RPL32＞TBP＞ACT＞βACT＞EF1α(Fig. 3-A).

In our analysis of multiple tissue types, both the ΔCt method and NormFinder suggestedEF1δ,EF1γ,andEF1αwere the most suitable reference genes (Table 2). In contrast,RPL32/RPS15andRPS15were considered the most stable genes by geNorm and BestKeeper, respectively(Table 2). Furthermore, all four algorithms identifledβACTas the least suitable reference gene (Table 2). For tissues,the overall RefFinder stability ranking (from most to least stable) was:RPS15＞EF1δ＞RPL32＞EF1γ＞EF1α＞RPL27＞SDH＞GAPDH＞βTUB＞ACT＞TBP＞βACT(Fig. 3-B).

For MTGs from different developmental stages and sexes, both the ΔCt method and NormFinder suggestedRPS15andRPL32were the most stable genes (Table 2),whereas BestKeeper rankedEF1αandRPL32as the most stable and geNorm rankedEF1γandEF1δas the best pair(Table 2).βACTwas again identifled by all four algorithms as the least stable reference gene (Table 2). For the MTG analyses, the RefFinder ranking was:RPL32＞RPS15＞EF1α＞EF1γ＞EF1δ＞GAPDH＞βTUB＞ACT＞TBP＞SDH＞RPL27＞βACT(Fig. 3-C).

In the last set of experiments, which assessed the effect of dsRNA-mediated RNAi on reference gene expression,ACTwas identifled by all four analyses as one of the most stable genes (Table 2). In addition,RPS15(ΔCt method and NormFinder),RPL32(BestKeeper), andβTUB(geNorm)were also identifled as having a calculated stability value equivalent to that ofACT(Table 2). The stability values(calculated by BestKeeper and geNorm) for all candidate genes were lower than the recommended threshold for reference gene suitability with SD＜0.460 for BestKeeper and M＜0.160 for geNorm (Table 2). This is consistent with the smaller variation in RNAi-injected treatments (Fig. 2-D).For the dsRNA injection study, the RefFinder ranking was:ACT＞RPS15＞RPL27＞βTUB＞RPL32＞EF1δ＞TBP＞SDH＞GAPDH＞EF1γ＞EF1α＞βACT(Fig. 3-D).

3.4. Determination of the optimal number of reference genes for normalization

Established methods of qRT-PCR often use a single reference gene with sufflcient expression for analysis,though use of more than one reference gene strengthens analysis (Vandesompeleet al.2002). Thus, we used geNorm to estimate the pairwise variation (Vn/Vn+1) to determine the optimal number of reference genes required for normalization of samples under a given experimental condition. All pairwise variations were determined to be below 0.15 (the recommended threshold of cut-off) for all treatments, indicating that no additional genes are required for the normalization (Fig. 4). Thus, the use of only the top two reference genes for each experimental set is sufflcient for normalization.

3.5. Validation of selected reference genes in A. suturalis

To validate the reference genes selected, we assessed the relative expression ofA. suturalisV-ATPaseandFASin various tissues. The best reference gene pair determined by geNorm,RPL32/RPS15, the single best reference gene determined by RefFinder,RPS15, and the most variable reference gene determined by all algorithms,βACT, were used to normalize the expression levels of those two target genes.

V-ATPaseis a highly conserved evolutionarily ancient enzyme that functions in cellular homeostasis, found at the plasma membrane of cells lining the gut and the Malpighian tubules of many insects, where it regulates pH, energizes ion transport, and modulates fluid secretion (Nelsonet al.2000;Wieczoreket al.2009). Normalization of transcripts usingRPL32/RPS15andRPS15alonerevealed peakV-ATPaseexpression in the gut (Fig. 5-A). In contrast, normalization withβACTsuggested the highest expression in MTG.FASis a multi-enzyme protein that catalyzes thede novosynthesis of long-chain fatty acids, and plays an important role in energy production and storage, cellular structure,and the biosynthesis of pheromone in insects (Volpe and Vagelos 1973; Tillmanet al.1999). We saw peak expression in the MTG whereA. suturalissynthesizes and releases pheromones (Zhanget al.2014) regardless of the reference gene used. Although the trend was similar, normalization with an unsuitable reference gene such asβACTnot only increased gene expression levels, but also resulted in larger standard error values (Fig. 5-B). Therefore, our results support the importance of the selection and validation of accurate reference genes RT-qPCR to avoid misinterpretation of the expression data.

4. Discussion

qRT-PCR is the most widely used molecular technique for gene expression analysis, but the accuracy and reliability of the results are critically dependent on appropriate data normalization with the use of stable reference gene(s)(Vandesompeleet al.2002; Bustinet al.2010). Using inappropriate reference gene(s) can signiflcantly impact quantiflcation results, leading to false inferences or misinterpretations (Radoni?et al.2004; Huggettet al.2005; Nolanet al.2006; Yuanet al.2014). Furthermore,numerous qRT-PCR studies have revealed that most reference/housekeeping genes exhibit variable expression depending on the organism and experimental conditions,which suggests that there is no suitable “universal” reference gene (Leeet al.2002). Therefore, it is essential to validate reference gene(s) prior to examining the effects of different experimental conditions on target gene expression.

Table 2 Rank order of the candidate Adelphocoris suturalis reference genes under different experimental conditions

Fig. 3 Expression stability and relative ranking of candidate reference genes calculated by the Geomean method of RefFinder.A, different developmental stages and sexes. B, different tissue. C, metathoracic scent glands (MTGs) from different developmental stages and sexes. D, dsRNA injection. A lower Geomean ranking indicates more stable expression.

AlthoughA. suturalisis one of the most destructive pests in major cotton growing regions (Jianget al.2015), its molecular physiology had not been actively explored due to a lack of molecular resources. However, recent transcriptomic advances (Luoet al.2014; Tianet al.2015) have opened the door for functional genomics research and associated gene-expression studies. However, the lack of dependable information regarding reference genes forA. suturaliscan limit data normalization and lead to false inferences or misinterpretations (Huggettet al.2005; Fergusonet al.2010). In this study, we used abiotic and biotic conditions to evaluate the expression stability of 12 housekeeping genes frequently used in expression studies of other sapsucking insects (Liet al.2013; Bansalet al.2016; Ibanez and Tamborindeguy 2016; Koramutlaet al.2016).

To better evaluate the candidate reference genes and to avoid analysis errors caused by selecting co-regulated transcripts, four statistical models (ΔCt method, geNorm,BestKeeper and NormFinder) were used to evaluate gene expression stability. We found that stability rankings frequently varied with the analysis package utilized. For example, when compared across multiple tissues,RPS15had the best BestKeeper and geNorm scores, but was not prioritized by NormFinder. These variations can be attributed to differences in the scaling system utilized by the algorithms (Zhaiet al.2014; Sagriet al.2017). Although the ranking orders fluctuated according to the analysis software used, the general trends were similar. Genes among the diverse tissues were largely split into two groups by all four algorithms (Table 2), those that were stably expressed(EF1δ,EF1γ,EF1α,RPL32,RPS15,RPL27andSDH) and those which exhibited variable expression patterns (GAPDH,ACT,TBP,βTUBandβACT). Thus, we used RefFinder, a comprehensive platform that integrates all four algorithms to rank the overall stability of candidate genes.

Fig. 4 Determination of the optimal number of reference genes for normalization by geNorm analysis. MTG, metathoracic scent gland. Average pairwise variations (Vn/Vn+1) were calculated by geNorm between the normalization factors NFn and NFn+1 to indicate whether inclusion of an extra reference gene would add to the stability of the normalization factor. A value below 0.15 indicates that an additional reference gene will not signiflcantly improve normalization.

Reference gene suitability can vary in response to diverse biotic and abiotic factors (Table 2; Fig. 3). This result is consistent with previous studies. For example,EF1αwas stably expressed inA. suturalisMTGs at different developmental stages and sexes, whereas its expression was highly variable at the whole-body level under similar conditions. The variation ofEF1αin different developmental stages and sexes may only occur in certain tissues, e.g.,eggs (Fig. 2). This has also been observed in many other species (Panet al.2015; Tanet al.2015; Bansalet al.2016; Koramutlaet al.2016). For example,the expression ofαTUB1inColaphellus bowringivaried across different developmental stages, though it is stably expressed under different photoperiods (Tanet al.2015).

Furthermore, our veriflcation results showed that normalization of transcripts usingRPL32/RPS15andRPS15alone revealed peakV-ATPaseexpression in the gut (Fig. 5-A), consistent with results fromCulex quinquefasciatus(Filippovaet al.1998),Holotrichia parallela(Weiet al.2016), andNasutitermes takasagoensis(Kumaraet al.2015). In contrast, normalization withβACTsuggested the highest expression in MTG. TheFASshowed a peak expression in the MTG regardless of the reference gene used. The same expression pattern was also observed inSpodoptera litura(Linet al.2018),Agrotis ipsilon(Guet al.2013) and Bumblebee (Zaceket al.2013). Although the trend was similar, normalization with an unsuitable reference gene such asβACTnot only increased gene expression levels, but also resulted in larger standard error values(Fig. 5-B). When the target gene was normalized to different reference genes that were either the most or least suitable for different tissue types, there was a dynamic shift in gene expression levels (Fig. 5). This shift in gene expression levels is driven by varied expression of the selected reference gene under varying experimental conditions,and contributes to biased conclusions that are based on calculations utilizing a faulty reference gene dataset(Vandesompeleet al.2002; Radoni?et al.2004; Huggettet al.2005; Nolanet al.2006; Ling and Salvaterra 2011).Previous studies have demonstrated that this variability under diverse conditions can affect transcript quantiflcation results, leading to false inferences or misinterpretations (Leeet al.2002; Yuanet al.2014; Wanget al.2015; Koramutlaet al.2016). All the above underscore the need for validation of reference genes prior to their experimental use.

Fig. 5 Validation of reference genes. The relative expression level of vacuolar-type H+-ATPase (V-ATPase, A) and fatty acid synthase (FAS, B) were normalized using the best reference pair (RPL32/RPS15 identifled by geNorm), the single best reference gene (RPS15 determined by RefFinder), and the least stable reference gene (βACT determined by all flve algorithms). The results are depicted as the mean±SE based on three independent biological replicates. Different letters indicate signiflcant differences(P＜0.05, one-way ANOVA followed by Tukey’s HSD Multiple Comparison).

Here, we found that the quantiflcation results of candidate reference genes varied across tissue types and developmental stages, and that exogenous dsRNA injection has minor impact on gene expression (Fig. 2). This type of variation has also been described inHalyomorpha halys(Bansalet al.2016) andColeomegilla maculate(Yanget al.2015b). We speculated that this variation of gene expression in the different tissue types and developmental stages is partly caused by variable RNA extraction efflciencies from different tissues of insects (Huggettet al.2005). To minimize the experimental errors caused by the variability of RNA extraction, a similar sample size with similar tissue volume and weight is needed. Then, RNA extraction steps should be optimized to accommodate different type of samples.Finally, it is essential to accurately quantify and assess the quality of RNA prior to reverse transcription (Bustin and Nolan 2004; Huggettet al.2005). RNAi is an effective tool to assess gene function that utilizes introduced dsRNAs to trigger a sequence-speciflc knockdown of gene expression at the post-transcriptional level (Hannon 2002). Measuring gene expression is thus critical for conflrming the reduction in target transcript levels following RNAi-mediated silencing.Many dsRNA delivery systems have been used successfully for this purpose. Among these, injection remains the most frequently used method due to its high efflciency and accuracy (Hughes and Kaufman 2000). The mechanical stress exerted by puncturing the cuticle, or the stress associated with the introduction of exogenous dsRNA, could alter housekeeping gene expression, and affect study of target gene silencing. Therefore, it is crucial to identify and validate reference gene expression stability under variable conditions in order to develop RNAi. Our study found that all used candidate reference genes remained stable under the experimental treatments (Table 2, Fig. 3). However, it is possible that the stability of these reference genes could change when injected with other dsRNAs, and further study on this topic is required.

5. Conclusion

Given the economic importance ofA. suturalisas an agricultural pest, the application of gene expression analyses and functional genomics can further facilitate basic and applied research on this species. Hence, it is critical to establish a standardized qRT-PCR analysis. The current study identifled stable reference genes inA. suturalisunder an array of biotic and abiotic conditions. Based on a comprehensive analysis integrating flve commonly used methods to compare and rank the expression stability of candidate reference genes, we recommend the following genes as the best suited for use inA. suturalis: 1)EF1δandEF1αfor gene expression studies across different developmental stages and sexes; 2)RPS15andEF1δfor tissue proflling; 3)RPL32andRPS15for MTGs at different developmental stages and sexes; and 4)ACTandRPS15for dsRNA-mediated RNAi studies. This is the flrst study to evaluate suitable reference genes for gene expression analyses inA. suturalisand the results will facilitate future research into the transcriptional changes associated with various biological processes.

Acknowledgements

This work was funded by the National Special Key Project for Transgenic Breeding, China (2016ZX08011002) and the Open Fundation of State Key Laboratory for Biology of Plant Diseases and Insect Pests, Institute of Plant Protection,Chinese Academy of Agricultural Sciences (SKLOF201415).Mention of trade names or commercial products in this article is solely for the purpose of providing speciflc information and does not imply recommendation or endorsement by the U.S.Department of Agriculture (USDA) is an equal opportunity provider and employer.

Appendicesassociated with this paper can be available on http://www.ChinaAgriSci.com/V2/En/appendix.htm