Selection and validation of experimental condition-specific reference genes for qRT-PCR in Metopolophium dirhodum (Walker) (Hemiptera: Aphididae)

Metopolophium dirhodum (Walker) (Hemiptera: Aphididae) is one of the most common aphid pests of winter cereals. To facilitate accurate gene expression analyses with qRT-PCR assays, the expression stability of candidate reference genes under specific experimental conditions must be verified before they can be used to normalize target gene expression levels. In this study, 10 candidate reference genes in M. dirhodum were analyzed by qRT-PCR under various experimental conditions. Their expression stability was evaluated with delta Ct, BestKeeper, geNorm, and NormFinder methods, and the final stability ranking was determined with RefFinder. The results indicate that the most appropriate sets of internal controls were SDHB and RPL8 across geographic population; RPL8, Actin, and GAPDH across developmental stage; SDHB and NADH across body part; RPL8 and Actin across wing dimorphism and temperature; RPL4 and EF1A across starvation stress; AK and RPL4 across insecticide treatments; RPL8 and NADH across antibiotic treatments; RPL8, RPL4, Actin, and NADH across all samples. The results of this study provide useful insights for establishing a standardized qRT-PCR procedure for M. dirhodum and may be relevant for identifying appropriate reference genes for molecular analyses of related insects.

. Candidate reference gene expression levels. Candidate reference gene expression levels in the whole M. dirhodum sample set are expressed in terms of the threshold cycle number (Ct value). Data are presented as whisker box plots. The box represents the 25th-75th percentiles, the median is indicated by a bar across the box, and the whiskers on each box represent the minimum and maximum values. Wing dimorphism. The delta Ct and BestKeeper analyses identified Actin and RPL4 as the most stable genes, whereas both NormFinder and GeNorm identified RPL8 and EF1A as the most stable genes. All four analyses indicated that 18S, GAPDH, AK, and SDHB were the least stable genes ( Table 1). The RefFinder data for the wing dimorphism revealed a rank order (most to least stable expression) of RPL8, Actin, RPL4, EF1A, HSP68, NADH, SDHB, AK, GAPDH, and 18S (Fig. 2). On the basis of the GeNorm analysis, all pairwise variation values were below the 0.15 cut-off value (Fig. 3). According to RefFinder, RPL8 and Actin are required for normalizing target gene expression levels in wing-dimorphic insects.
Temperature-induced stress. The delta Ct method identified Actin and RPL8 as the most stable genes. Both BestKeeper and NormFinder identified RPL8 and RPL4 as the most stable genes, whereas GeNorm identified Actin and NADH as the most stable genes. All four analyses indicated HSP68, 18S, GAPDH, and AK were the   2). On the basis of the GeNorm analysis, all pairwise variation values were below the 0.15 cut-off value, except for V9/10 ( Fig. 3). The RefFinder analysis indicated RPL8 and Actin are required for normalizing target gene expression levels in M. dirhodum exposed to different temperatures.
Starvation-induced stress. The delta Ct method and the NormFinder algorithm identified EF1A and RPL4 as the most stable genes and Actin, 18S, and HSP68 as the least stable genes (Table 1). However, BestKeeper identified 18S and Actin as the most stable genes and SDHB and HSP68 as the least stable genes ( Table 1). The GeNorm algorithm identified NADH and AK as the most stable genes and Actin, 18S, and HSP68 as the least stable genes ( Table 1). The RefFinder results for the starvation treatment revealed a rank order (most to least stable expression) of RPL4, EF1A, AK, NADH, RPL8, 18S, GAPDH, Actin, SDHB, and HSP68 (Fig. 2). The GeNorm analysis indicated that the pairwise variation values for V2/3 were less than the proposed 0.15 cut-off (Fig. 3). The RefFinder analysis indicated RPL4 and EF1A are required for normalizing target gene expression levels in starvation-stressed M. dirhodum.
Insecticide-induced stress. The delta Ct and NormFinder data revealed AK and RPL4 as the most stable genes, whereas the BestKeeper results identified HSP68 and SDHB as the most stable genes. In contrast, Actin and AK were the most stable genes according to GeNorm. All four analyses identified 18S and EF1A as the least stable genes ( Table 1). The RefFinder data for the insecticide treatment revealed a rank order (most to least stable expression) of AK, RPL4, Actin, RPL8, HSP68, NADH, SDHB, GAPDH, 18S, and EF1A ( Fig. 2). Based on the GeNorm analysis, all the pairwise variation values were below 0.15 cut-off value (Fig. 3). Thus, AK and RPL4 are required for normalizing target gene expression levels in insecticide-treated M. dirhodum.
Antibiotic-induced stress. The delta Ct method identified RPL8 and RPL4 as the most stable genes. The Best-Keeper algorithm identified SDHB and Actin as the most stable genes, whereas NormFinder indicated NADH and RPL8 were the most stable genes. The GeNorm algorithm identified GAPDH and 18S as the most stable genes. All four analyses identified EF1A, SDHB, and HSP68 as the least stable genes ( Table 1). The RefFinder data for the antibiotic treatment revealed a rank order (most to least stable expression) of RPL8, NADH, RPL4, 18S, GAPDH, AK, Actin, SDHB, EF1A, and HSP68 (Fig. 2). According to the GeNorm analysis, all pairwise variation values were less than the proposed 0.15 cut-off, except for V9/10 (Fig. 3). The RefFinder analysis suggested RPL8 and NADH are required for normalizing the target gene expression levels in antibiotic-treated M. dirhodum.
Overall ranking of M. dirhodum candidate reference genes. An examination of the candidate reference gene expression stability for all treatments and conditions with the four methods used in this study produced similar rank orders, with RPL4 and RPL8 identified as the most stable genes and AK, 18S, and HSP68 revealed as the least stable genes ( Table 1). The RefFinder results for all treatments and conditions revealed a rank order (most to least stable expression) of RPL8, RPL4, Actin, NADH, EF1A, SDHB, GAPDH, 18S, AK, and HSP68 (Fig. 2). The GeNorm analysis indicated that the pairwise variation values for V4/5 were less than the proposed 0.15 cut-off (Fig. 3). Thus, an analysis of all treatments and conditions suggested that RPL8, RPL4, Actin, and NADH are suitable internal reference genes for normalizing target gene expression levels in M. dirhodum. Figure 3. Determination of the optimal number of reference genes for accurate normalization calculated by geNorm. The V n/n+1 value indicates the pairwise variation (Y axis) between two sequential normalization factors and determines the optimal number of reference genes required for an accurate data normalization. A value below 0.15 indicates that an additional reference gene will not significantly improve the normalization.

Discussion
There are several reports describing the application of qRT-PCR assays to clarify the gene expression levels associated with diverse biological processes [36][37][38][39] . Reference genes used for molecular investigations can influence the accuracy of target gene expression levels 6,[40][41][42] . Therefore, a stable reference gene is an important prerequisite for gene expression investigations. Housekeeping genes, which are constitutively expressed to maintain basic cellular functions, have traditionally been used as internal reference controls 6,10,11 . However, there is no universal reference gene that is stably expressed in all cell and tissue types under different experimental conditions 10,11,[43][44][45][46][47] . Therefore, every stable reference gene used to normalize gene expression data should be evaluated under each experimental condition 43,48 .
In this study, qRT-PCR was used to evaluate the expression-level stability of 10 candidate reference genes in M. dirhodum across specific conditions. The best reference genes varied among conditions. Specifically, RPL8 (mean Ct value ± SD, 21.10 ± 0.35) and Actin (26.79 ± 0.42) had the least variable expression levels, whereas HSP68 (24.82 ± 1.86) produced the most variable expression levels among the examined candidate reference genes (Fig. 1). Similarly, RPL8, RPL4, and Actin were the most stable reference genes, whereas HSP68 and 18S were the least stable reference genes under most conditions (Fig. 2).
Ribosomal proteins (RPs), which are the principal components of ribosomes, are one of the most highly conserved proteins in all life forms. Earlier research proved that RP-encoding genes are among the most stably expressed reference genes, and have been widely used to normalize gene expression levels in insect molecular investigations during the past 10 years 49 . For example, in Bradysia odoriphaga 50 , RPS15 was the most stably expressed gene in response to various temperature treatments. However, another study indicated that the expression levels of RP-encoding genes may vary under some conditions 49 . Moreover, RPS20 was detected as the least stably expressed gene for analyzing Plutella xylostella geographic populations as well as the effects of the temperature, photoperiod, and insecticides 10 . Consistent with these earlier findings, we identified RPL8 as the most stable gene in M. dirhodum across various conditions (except for analyses of different body parts, starvation stress, and insecticide treatments) (Fig. 2). Additionally, RPL4 was detected as the most stable gene in response to starvation and insecticide treatments, but was also almost the least stable gene during analyses of various M. dirhodum body parts (Fig. 2).
Actin, which encodes a major structural protein, is important for cell secretion, motility, cytoplasm flow, and cytoskeleton maintenance. Moreover, Actin is expressed at various levels in many cell types, and is considered the ideal reference gene for qRT-PCR, which may explain its frequent use 15,26 . For example, it has been used to study the effects of diet on B. odoriphaga gene expression 50 and for investigating M. persicae gene expression in different tissues and in response to the temperature, photoperiod, and wing dimorphism 26 . However, in Helicoverpa armigera, Actin was revealed to be the least stable reference gene following temperature and photoperiod treatments 51 . In our study, Actin was identified as one of the most stable reference genes for analyzing developmental stages, temperature effects, and wing dimorphism (Fig. 2).
The GAPDH gene has been commonly used as a reference gene in the studies of gene expression 7,52,53 . However, unstable GAPDH expression has been detected in Tetranychus cinnabarinus developmental stages 54 , in the labial glands and fat bodies of Bombus terrestris and Bombus lucorum 55 , and in various Sogatella furcifera body parts 56 . In the current study, GAPDH was revealed as a stably expressed candidate reference gene for analyses of developmental stages (Fig. 2). These results imply that the mechanism underlying the expression stability of endogenous reference genes is complex. Furthermore, the stability of potential reference genes in different biological samples should be tested prior to their use.
The protein encoded by EF1A affects translation by catalyzing the GTP-dependent binding of aminoacyl-tRNA to the acceptor site of the ribosome. The EF1A gene was recently used as a reference gene in multiple insect gene expression studies 55,57,58 . Our results suggest that EF1A is an appropriate reference gene only for analyzing the effects of starvation stress on M. dirhodum gene expression ( Table 2).
The AK gene encodes the phosphagen kinase in invertebrates, and it has rarely been used as a reference gene 59 . An earlier study revealed that AK is the most stably expressed gene in the B. terrestris labial gland and fat body 60 . In this study, AK was identified as the most stable gene following insecticide treatments (Fig. 2). In A. pisum, SDHB and NADH are reportedly the most stable housekeeping genes in developmental stages and in response to various temperatures 11 . However, we determined that SDHB and NADH are the most stable housekeeping genes only during examinations of different M. dirhodum body parts (Fig. 2). These results further suggest that reference gene expression stability is influenced by the experimental conditions.
The 18S rRNA gene is considered to be an ideal reference control because of its relatively stable expression levels 61 . Accordingly, it has been applied in previous studies involving Lucilia cuprina 62 , Rhodnius prolixus 63,64 , and Delphacodes kuscheli 65 . However, in this study, 18S was revealed as one of the least stable genes in almost all www.nature.com/scientificreports/ sample sets, implying it is an inappropriate reference gene for M. dirhodum (Fig. 2). This observation is consistent with the results of previous studies that indicated 18S rRNA is not a stable reference gene in Bactrocera dorsalis and Nilaparvata lugens under specific experimental conditions 66 . It is transcribed by a separate RNA polymerase, which may explain why rRNA is not a suitable reference control 67 . Moreover, the utility of 18S for normalizing target gene expression levels in a qRT-PCR assay is limited by the potential imbalance between rRNA and mRNA fractions among samples 61 . The HSP68 gene, which belongs to the HSP70 family, encodes a highly conserved chaperone involved in protein assembly, folding, and transport as well as in antigen processing and presentation. The expression of genes encoding HSPs can be affected by high temperatures or other stresses (e.g., due to chemicals) 68 . In the current study, HSP68 was the least stable gene for all conditions (Fig. 2). In a previous study on Coleomegilla maculata, HSP70 was identified as the most stably expressed gene for sexes, but was the least stably expressed gene for analyses of different tissues, and dsRNA exposure 44 .
It is becoming common for researchers to use multiple reference genes to normalize target gene expression levels in diverse studies because a single gene is usually insufficient for analyzing gene expression 69 . An earlier investigation indicated that too many or too few reference genes may adversely affect the robustness of data normalizations 70 . However, the simultaneous application of multiple reference genes in a given experiment may decrease the probability of biased normalizations. The optimal number of reference genes under specific experimental conditions can be determined with the geNorm algorithm, which calculates the pairwise variation V n/n+1 based on the normalization factors NF n and NF n+1 , with n ≥ 2. If V n/n+1 is below 0.15, n is the optimal number of reference genes. The results of this study indicate that the most appropriate number of reference genes varies under diverse experimental conditions (Fig. 3). This implies that the stability of reference genes must be evaluated before every qRT-PCR experiment.

Conclusions
To the best of our knowledge, this study is the first to evaluate and validate experimental condition-specific candidate reference genes for M. dirhodum gene expression analyses. We identified reference genes applicable for elucidating functional gene expression profiles. In this study, we examined 10 candidate reference genes under diverse conditions. Notably, the stability of candidate gene expression levels in M. dirhodum varies depending on the experimental conditions. Moreover, we identified internal reference genes suitable for normalizing and quantifying gene expression in M. dirhodum (Table 2). Our findings may be useful for establishing a more accurate and reliable method for normalizing M. dirhodum qRT-PCR data. They may also provide the basis for future investigations on RNA interference and gene transcription in M. dirhodum and other insect pests.

Materials and methods
Insects. Our  Body part. We used a dissection needle and a tweezer to separate the head, thorax, and abdomen from wingless M. dirhodum adults. These body parts as well as whole adult bodies were stored as described earlier.
Wing dimorphism. Three samples of 20 winged and wingless M. dirhodum adults were collected, flash frozen in liquid nitrogen, and stored at − 80 °C until total RNA extraction.
Temperature-induced stress. Potted wheat seedlings infested with M. dirhodum were divided into five groups for a 24-h exposure to one of the following five temperatures: 4, 10, 15, 20, and 25 °C. For each temperature, three samples of 20 adults were collected, flash frozen in liquid nitrogen, and stored at − 80 °C until total RNA extraction. None of the temperature treatments were lethal to the aphids.
Starvation-induced stress. Adult aphids were placed on moistened filter paper in a Petri dish (9 cm diameter) with no food for a 32-h incubation in a thermostatic chamber (20 ± 2 °C and 60% relative humidity, with a 16-h light:8-h dark cycle). The control (satiated) group comprised aphids able to feed on wheat seedlings in the same conditions. For the control and treatment groups, three samples of 20 adults were collected, flash frozen in liquid nitrogen, and stored at − 80 °C until total RNA extraction. The mortality rate among the starved aphids was approximately 10%.

Scientific Reports
| (2020) 10:21951 | https://doi.org/10.1038/s41598-020-78974-z www.nature.com/scientificreports/ Insecticide-induced stress. The effects of insecticides on the stability of candidate reference genes were assessed in M. dirhodum subjected to one of the following three insecticide treatments: imidacloprid (9.87 mg/L), thiamethoxam (122.00 mg/L), and beta-cypermethrin (17.28 mg/L). These concentrations were selected because a bioassay indicated they are 30% to the mortality of the population (LC 30 ) (Table S1). Aphids were treated with the insecticides via the leaf dip method 71 . The 1% insecticide stock solutions prepared in acetone were serially diluted with water (containing 0.1% Tween-80) to produce five concentrations. Water (containing 0.1% Tween-80) was used as a control solution. Wheat leaves with M. dirhodum were immersed in the prepared solutions for 3-5 s and then placed on moistened filter paper in a Petri dish (9 cm diameter). The samples were incubated for 24 h at 20 ± 2 °C and 60% relative humidity, with a 16-h light:8-h dark cycle. For each concentration, the mortality rate based on three replicates of 30 aphids was calculated. Additionally, for the control and treatment groups, three samples of 20 adults were collected, flash frozen in liquid nitrogen, and stored at − 80 °C until total RNA extraction.
Antibiotic-induced stress. The M. dirhodum adults were fed a 30% sucrose solution containing 50 µg/mL rifampicin or an antibiotic-free sucrose solution (control) ( Table 3. Functions, primer sequences, and amplicon characteristics of the candidate reference genes analyzed in this study. a Amplicon length. b qRT-PCR efficiency (based on a standard curve). c Reproducibility of the qRT-PCR. www.nature.com/scientificreports/ was set to 500 to obtain all the threshold cycle (Ct) values of tested candidate reference genes. The qRT-PCR analyses were completed with three biological replicates and three technical replicates.
Data analysis. The stability of the 10 candidate reference housekeeping genes was evaluated with the geNorm 40 , NormFinder 73 , and BestKeeper 74 algorithms and the comparative delta Ct method 75 . Finally, we compared and ranked the tested candidate reference genes with the web-based RefFinder analytical tool (https :// www.heart cure.com.au/for-resea rcher s).