Metabolic detoxification is a common mechanism of insecticide resistance, in which detoxifying enzyme genes are overexpressed. Aphis glycines Matsumura (Hemiptera: Aphididae) is one of the major soybean pests in the United States and has developed resistance to pyrethroid insecticides after almost two decades of use. To date, there are no validated reference genes to normalize expression of detoxification genes for pyrethroid resistance in A. glycines. From a literature review, a list was compiled of genes from 36 gene families (68 sequences) frequently used as reference genes in gene expression analysis in Hemiptera. Exon–exon junction primers were designed for the best alignment matches to a draft A. glycines genome and were assayed in a three-phase screening. The first screen eliminated nonamplifying primers. The second screen used nine A. glycines populations varying in resistance to pyrethroids and eliminated primers with inconsistent amplification or low amplification efficiency, and quantitatively assessed the stability of expression in the 14 remaining candidates using NormFinder and a generalization of BestKeeper. The third screen quantitatively validated these results on the best candidates. Six genes were identified with the greatest stability across technical and biological replication and the nine populations. The genes identified as the most suitable reference genes for the study of detoxifying enzymes related to pyrethroid resistance in soybean aphid were: actin, RPL9 (ribosomal protein L9), RPS9 (ribosomal protein S9), AK (arginine kinase), RNAPol2 (RNA polymerase II), and RPL17 (ribosomal protein L17). Our findings will support studies related to insecticide resistance in A. glycines.
Soybean, Glycine max (L.) Merrill, is an important commodity crop in the U.S. economy. The major insect pest affecting soybean in the North Central Region of the United States and southeast Canada is soybean aphid, Aphis glycines Matsumura (Hemiptera: Aphididae) (Ragsdale et al. 2011). Since the detection of A. glycines in 2000 in the United States (Hartman et al. 2001), chemical control has been the most widely used control tactic (Hodgson et al. 2012, Koch et al. 2018). Reliance on the few insecticide groups available, primarily pyrethroids and organophosphates, favored the selection of A. glycines resistance to pyrethroids (Hanson et al. 2017, Koch et al. 2018, Menger et al. 2020). Several mechanisms of resistance to insecticides have been documented in insects. Among these mechanisms is metabolic resistance (Feyereisen 1995), which results from the overexpression of detoxification enzymes, such as cytochrome P450 monooxygenases, glutathione-S-transferases, and esterases (E4 and CES) (Coppin et al. 2012, Panini et al. 2016). For A. glycines, there is evidence of metabolic resistance to pyrethroids in China (Xi et al. 2015) and in the United States (Paula et al. 2020). Metabolic resistance is the most common mechanism of insecticide resistance, and it can present a considerable challenge in insecticide resistance management (IRAC - Insecticide Resistance Action Committee 2019). Therefore, the monitoring of metabolic resistance is imperative to improve integrated pest management (IPM) programs (Li et al. 2016) for a pest such as A. glycines with a pyrethroid resistance history.
Real-time reverse transcription–quantitative polymerase chain reaction (RT-qPCR) is a sensitive, practical, and low-cost tool for expression analysis of detoxification genes (Bansal et al. 2012, Kozera and Rapacz 2013). An essential component of gene expression analysis by RT-qPCR is the selection and validation of appropriate reference genes (Huggett et al. 2005, Koramutla et al. 2016). Reference genes serve as endogenous or internal controls to normalize variability in the RT-qPCR signal across samples, usually introduced in the RNA extraction, complementary DNA (cDNA) synthesis, or PCR stages (Huggett et al. 2005, Kozera and Rapacz 2013), because they have constitutive expression that is not affected by the experimental conditions or treatments (Yang et al. 2014a) or the population to which the organism belongs (Lu et al. 2013, Sun et al. 2010). They are generally housekeeping genes (HKGs), so called because they are responsible for basic cell metabolism, among other functions (Butte et al. 2001, Thellin et al. 1999). Nevertheless, since it has been demonstrated that expression of HKGs may vary in certain circumstances in response to diverse biotic or abiotic factors (Thellin et al. 1999), it became advisable to use multiple reference genes (Lü et al. 2018, Radonić et al. 2004, Vandesompele et al. 2002).
According to Lü et al. (2018), the most common reference genes used in gene expression analysis in insects are: actin (main groups: alpha, beta, and gamma), ribosomal protein L (RPL), tubulin, elongation factor 1 alpha (EF1a), glyceraldehyde-3 phosphate dehydrogenase (GAPDH), ribosomal protein S (RPS), TATA-box binding protein (TBP), 18S ribosomal RNA (18S), heat shock protein (HSP), and succinate dehydrogenase subunit A (SDHA). However, their common use does not assure their validity as reference genes as their expression stability may vary in different insect species and experimental conditions (Gutierrez et al. 2008). Several reference genes were validated for many gene expression analyses in aphids, e.g., Acyrthosiphon pisum (Harris) (Yang et al. 2014a), Aphis gossypii Glover (Ma et al. 2016), and Myzus persicae (Sulzer) (Kang et al. 2017), including A. glycines (Bansal et al. 2012). However, reference genes have not yet been identified for the study of expression of detoxification genes related to pyrethroid resistance in A. glycines.
To accurately quantify and monitor the incidence of metabolic resistance to pyrethroids in field populations of A. glycines, this work aimed to identify reference genes across nine populations varying in resistance to pyrethroids. This assessment was performed as a three-phase screening for stability of expression of candidate reference genes from 36 gene families (68 sequences). Our results will provide a fundamental base for additional study of gene expression analysis in A. glycines under different levels of pyrethroid resistance.
Materials and Methods
Soybean aphids. Nine populations of A. glycines were studied. The populations were: (1) a laboratory susceptible control population (Biotype 1, previously shown to be insecticide susceptible [Hanson et al. 2017]) and field-collected populations from (2) Sutherland (IA), (3) Lamberton (MN), (4) Howard Lake (MN), (5) Rosemount (MN), (6) Fairfax (MN), (7) St. Paul (MN), (8) Rochester (MN) in 2019, and (9) Hancock (MN) in 2018. The field-collected populations from 2019 were collected from one infested soybean plant at each of five locations per field, with locations spaced at least 20 m apart. The Hancock population was collected in 2018 from a few plants within 0.5 m of each other at one location in the field. The pyrethroid (lambda-cyhalothrin) resistance of each field-collected population was characterized by a LC99 glass-vial bioassay (Menger et al. 2020). Clonal populations were started from survivors of the initial bioassays for each field-collected population, except for Howard Lake. All these populations were maintained in a greenhouse (University of Minnesota) in 60-cm2 cages containing healthy soybean plants (SD01-76R cultivar) at 25 ± 2°C, 18–22% relative humidity, and 16:8 h (L:D) photoperiod. For each clonal population, the level of pyrethroid (lambda-cyhalothrin) resistance was assessed with the abovementioned LC99 glass-vial bioassay. Tukey's test (P, 0.05) was performed to compare the mean proportion mortality among populations. From the clonal populations, apterous adult aphids were collected from the plants in the colonies and preserved for further work. For the Howard Lake population, apterous adult aphids were collected and preserved directly from the initial bioassay (i.e., a clonal population was not developed for this population). Three biological replicates from each of the nine populations (27 A. glycines) were studied for gene expression. Prior to preservation, aphids were inspected with a dissecting microscope, using RNase-free materials, to select only intact apterous adults. Individuals were transferred to 1.5-ml microtubes, flash-frozen in liquid nitrogen, submerged in RNAlater, and stored at –80°C.
Selection of the candidate reference genes and primer design. The literature was searched for candidate reference genes for RT-qPCR studies, and 36 gene families (68 sequences) were selected (Table 1) from nine aphids and one psyllid species: Aphis craccivora Koch (Yang et al. 2015), A. glycines (Bansal et al. 2012), A. gossypii (Ma et al. 2016), Diaphorina citri Kuwayama (Bassan et al. 2017), Diuraphis noxia (Mordvilko) (Sinha and Smith 2014), Lipaphis erysimi (Kaltenbach) (Koramutla et al. 2016), Megoura viciae Buckton (Cristiano et al. 2016), Myzus persicae (Kang et al. 2017), Rhopalosiphum padi (L.) (Wu et al. 2014), and Toxoptera citricida (Kirkaldy) (Aphis citricidus (Kirkaldy)) (Shang et al. 2015). Using their GenBank accession numbers, the nucleotide sequences were retrieved and used as queries to search for orthologous A. glycines sequences within the database consisted of contigs from the “Genome Assembly v1.0 of Aphis glycines, Biotype 4 (Ag_bt4)” (Wenger et al. 2017) available in the “A. glycines blast server” supported by the Bioinformatics Platform for Agroecosystem Arthropods (bipaa. genouest.org/sp/aphis_glycines/blast). For the search, BLASTx 2.6.0+ (Altschul et al. 1997) was used with default parameters. Considering only the best alignment matches, 32 gene candidate sequences were obtained (from 26 gene families) (Table 2). Exon and intron positions were identified using AUGUSTUS (http://bioinf.uni-greifswald.de/augustus/submission) to design exon–exon junction primer-pairs. Two exon–exon junction primer-pairs were designed for almost all of them, totaling 58 primer-pairs, using IDT's PrimerQuest® Tool 2012 (https://www.idtdna.com/pages/tools/primerquest) with the parameters: melting temperature 58°C (minimum), 60°C (ideal), and 63°C (maximum); GC content 35% (minimum), 48% (ideal), and 50% (maximum); primer length 18 bp (minimum), 22 bp (ideal), and 25 bp (maximum); amplicon length 80 bp (minimum), 130 bp (ideal), and 200 bp (maximum); and to target exon–exon junctions (3′–5′). The GenBank accession numbers of the 32 A. glycines nucleotide sequences and the 58 primer-pair sequences are presented in Table 2.
RNA extraction and cDNA synthesis. Single apterous adult aphids were individually transferred to 2.0-ml screw-cap tubes containing two 5-mm borosilicate beads and 100 ll of RNeasy Lysis Buffer (RNeasy Mini Kit Qiagen, Venlo, Netherlands). Aphids were homogenized at 4 m/s for 20 s in a FastPrep-24™ homogenizer (MP Biomedicals, Irvine, CA). Total RNA was extracted using the RNeasy Mini Kit according to the manufacturer instructions. The RNA yield was measured using Qubit™ RNA HS Assay Kit and Qubit™ 3 Fluorometer (Thermo Fisher Scientific, Waltham, MA). Twenty nanograms of RNA from each sample was used to synthesize first-strand cDNA using SuperScript® IV RT (Invitrogen–ThermoScientific, Waltham, MA), according to the manufacturer instructions. A rough estimate of the cDNA obtained per sample was performed in a NanoDrop 2000 spectrometer (Thermo Fisher Scientific) to normalize the amount of cDNA across samples. A single 20-fold dilution was performed for each cDNA sample with nuclease-free water for the RT-qPCR analysis.
RT-qPCR. The RT-qPCR was performed in a LightCycler® 480 Instrument II (Roche, Basel, Switzerland) using Maxima SYBR Green/ROX qPCR (Thermo Fisher Scientific) and 384-well plates. Each sample had three technical replicates and each primer-pair had no template controls (NTCs). The RT-qPCR program for all the primers consisted of: one cycle of initial denaturation at 95°C for 10 min, followed by 45 cycles of two-step amplification process (denaturation at 95°C for 10 s, annealing and extension at 60°C for 60 s), and subsequent melting curve with temperature increase of 1°C/s starting at 40°C for 1 min and going to 95°C, then final cooling at 40°C for 10 s. We performed three-phase screening with RT-qPCR analysis. In the first screening, we used cDNA from the three biological replicates from the Hancock population to test primer-specificity and efficiency, as well as to check for the possibility of nonexpressing genes (nonfunctional copies or pseudogenes), for the 58 primer-pairs of the 32 candidate reference genes (Table 2). In the second screening, genes/primer-pairs retained from the first screening (Table 2) had expression tested in the three biological replicates from the nine A. glycines populations to select genes with more stable expression across replicates and populations. In the third screening, the most promising candidate reference genes from the second screening results were analyzed again as to validate the results of the second screening (Table 2).
Stability of gene expression. The raw fluorescence from each well was used to generate melting curves using the MBmca package (Nucleic Acid Melting Curve Analysis on Microbead Surfaces) (Rödiger et al. 2014) and the function diffQ to calculate the melting temperatures (Tm) from the first derivative (Rödiger et al. 2014). After verification of presence of single peaks and no amplification in the NTCs, we used LinRegPCR version 2014.5 (Ruijter et al. 2009) to analyze the RT-qPCR raw fluorescence data and estimate Log10 (N0), which is equivalent to efficiency-corrected Cq (quantification cycle) for all samples. In the first screening, primers with no amplification, nonspecific primers with multiple peaks in the melting curve, and primers with low amplification efficiency were eliminated. In the second screening, primers that did not amplify any one of the biological replicates or averaged more than one nonamplifying technical replicate were disregarded, leaving 14 candidate reference genes for quantitative analysis. Missing values were imputed for subsequent analysis using multilevel multiple imputation (50 times) with the jomo package (Quartagno and Carpenter 2019) in R. For each imputation, we calculated five measures of gene stability and the Pearson correlation coefficient for all pair-wise gene expression levels across the populations and biological and technical replicates. The results from each imputation were averaged across the 50 imputations. To represent the clustering of candidate reference gene groups a principal components analysis (PCA) was performed (Pearson 1901) using the average correlation matrix. This reveals potential reference genes that provide complementary measures of stability and correlated genes that reinforce each other. Three of the stability measures were calculated from a generalization of the BestKeeper method (Pfaffl et al. 2004). The BestKeeper method uses the standard deviation of gene expression across samples to estimate stability, with the smallest standard deviation being the most stable. Our data were structured to be able to estimate variation among technical replicates, biological replicates (individual aphids), and populations. We estimated the standard deviation (SD) among technical replicates within aphids for each gene (pooled across aphids within a population), the SD among aphids (biological replicates) within a population for each gene (pooled across populations), and the SD among populations for each gene as three measures of stability. We also used NormFinder version 2015 (Andersen et al. 2004) to estimate the stability of expression for individual genes and for pairs of genes. We did not use geNorm (Vandesompele et al. 2002), because, as shown by Andersen et al. (2004) in their supplementary information, geNorm assumes independence of expression among candidate genes and selects the gene that is most similar to the other genes tested and, therefore, can give erroneous results. The five measures of stability were normalized by transforming each to standard normal deviates across the candidate genes. To determine the most stable genes, the five normalized measures were averaged and the genes with the lowest values were selected as the most stable genes. Finally, using the results from the PCA, we selected genes that were independently expressed to have multiple measures of the reference genes, as well as correlated genes (similar to BestKeeper) to reduce variation in expression associated with one gene.
Pyrethroid susceptibility among soybean aphid populations. The soybean aphids from Biotype 1 (susceptible control), Sutherland, and Lamberton showed the highest mean mortalities (high susceptibility) of 1.00, 0.82, and 0.97, respectively (Fig. 1). Soybean aphids from Howard Lake and Hancock presented an intermediate susceptibility (0.60 and 0.45, respectively) and those from St. Paul, Rochester, Rosemount, and Fairfax were the least susceptible (0.15, 0.02, 0.02, and 0.00, respectively) (Fig. 1).
Amplification of the candidate reference genes (Screen 1). We chose the best primer-pair for each of the candidate reference genes based on the amplification efficiency showing a single peak in the melt curve analysis, and the R2 of the regression to estimate Cq (Schmittgen and Livak 2008) in RT-qPCR. Across all 58 primer-pairs (from 32 candidate reference genes) assessed in the initial screening, amplification efficiency ranged from 0.000 to 2.071 and the R2 of the regression to estimate Cq ranged from 0.499 to 1.000. A total of 30 primer-pairs for 30 sequences of the candidate reference genes were selected from the first RT-qPCR screening to proceed to the next screening (Table 2). The primer-pairs selected to advance to the second screening had amplification efficiencies ranging from 1.720 to 1.982 and R2 ranging from 0.683 to 0.987. The candidate reference genes helicase and RPL27 were eliminated in the first screening.
Expression stability (Screen 2). Out of the 30 sequences of the candidate reference genes analyzed in the second screening (Table 2), 16 were removed from consideration because they had at least one biological replicate with nonamplification or averaged more than one nonamplifying technical replicate. The remaining 14 candidate sequences, with an average 12.08% nonamplification of technical replicates, were statistically analyzed for stability (Fig. 2). The NormFinder single gene stability analysis indicated that β-actin, AK (arginine kinase), RPL17 (ribosomal protein L17), RNAPol2 (RNA polymerase II), and RPL5 (ribosomal protein L5) had the highest stability (Fig. 3). Better stability (<0.25) was obtained with NormFinder paired gene stability analysis. Actin, RPS18 (ribosomal protein S18), RPL5, β-actin, AK, RPL9 (ribosomal mitochondrial protein L9), RPL17, and RPS9 (ribosomal protein S9) contributed to the highest stability in gene pairs (Fig. 4). With BestKeeper, stable candidate reference genes exhibit a standard deviation <1 (Sundaram et al. 2019). The most stable candidate reference genes from the generalized BestKeeper analysis differed by source of variation (i.e., technical, biological, and between-population variation) (Fig. 5). In terms of technical variation, actin, RPS18, β-actin, AK, and RPL5 were the most stable genes. In terms of biological variation (i.e., variation among aphids within populations), AK, actin, RPL9, RNAPol2, and RPL17 were the most stable genes. In terms of variation among populations, actin, β-tubulin, α-tubulin, RPL9, and AK were the most stable genes. The average of the normalized values of the five stability measures (Table 3) indicated that AK, actin, β-actin, and RPL9 were the most stable genes. In addition, RNAPol2, RPL5, and RPL17 were more stable than the average candidate reference genes.
Two groups of genes were identified based on the PCA of the 14 sequences of the candidate reference genes (Fig. 6). The first PC axis (PC1) explained 88.68% of the variance in expression and the first two axes explained 94.05% of the variation. PC1 separated the candidate reference genes into two groups. The first group of four genes was actin, β-tubulin, α-tubulin, and RPL9. Three of these, actin, β-tubulin, and α-tubulin, were highly correlated with each other (0.70–0.85), and RPL9 was less correlated with them (0.67–0.71). These four genes were not highly correlated with the genes in the second group (0.004–0.60). Of these four genes, actin and RPL9 had the highest expression stability (Table 3) and were selected for third level of screening (i.e., validation analysis). The other two genes, α-tubulin and β-tubulin, had poor stability and were disregarded from further consideration as a reference gene.
The second group comprised 10 highly correlated candidate reference genes. The six genes forming the core of this group, RPS18, RPS9, β-actin, RPL5, RPL7, and 60S, were highly correlated (0.50–0.60) with stabilities ranging from 0.79 to 0.97. Of these, β-actin was the most stable (Table 3) and was retained for validation. RPS9 and 60S were the most unrelated to other candidates in this core group (Table 4), and as RPS9 had higher stability than 60S, it was retained for validation. The remaining four genes in this group were less correlated with the core group (0.52–0.82), and therefore could provide partially complementary information to the core. RNAPol2 and RPL17 had the highest stability and were selected for validation. AK was weakly associated with the core group (0.52–0.76) and was retained because it might provide the most complementary information to the core. These four genes could provide complementary information to the core group and to each other.
In summary, the second screening of candidate reference genes resulted in seven genes for validation: actin, β-actin, RPL9, RPS9, AK, RNAPol2, and RPL17. These genes were distributed throughout the PCA space, as would be expected for complementary gene expression.
Validation of the selected candidate reference genes (Screen 3). In the third level of screening, β-actin had >50% of wells without amplification and, therefore, this gene was disregarded as a reference gene. The candidate reference genes actin, RPL9, and RPS9 had better stability characteristics than in the second screening; however, AK and RNAPol2 had poorer stability, and RPL17 had similar stability (Tables 3, 5). Correlation analysis at this level of screening showed that actin and RPL9 were highly correlated, but they were only moderately correlated in the second level of screening. The other genes were not highly correlated, and most were less correlated than in the previous screening, indicating that they may provide complementary expression.
Even though there were differences in the collection methods and in the level of resistance of the A. glycines populations (different phenotypes), we found good stability in the gene expression of the candidate reference genes (Table 6). Finally, after the third level of screening, six genes (actin, RPL9, RPS9, AK, RNAPol2, and RPL17) were validated as suitable candidate reference genes for gene expression analysis in different A. glycines populations with contrasting levels of pyrethroid resistance.
The accuracy of the RT-qPCR data analysis depends on an adequate selection of reference genes (Everaert et al. 2011). It has been found that for trustworthy results more than one reference gene should be used in an experiment because together they improve stability among samples (Radonić et al. 2004, Vandesompele et al. 2002). Reference genes that are highly but not perfectly correlated provide measures of the same or similar biological process. Having multiple highly correlated reference genes would reduce the influence of methodological and other sources of technical variation but would not reduce variation associated with biological replication or among populations. To reduce error associated with these sources of variation, the expression profiles of the candidate genes should not be highly correlated. Reference genes that are not highly correlated provide complementary information and will reduce the influence of variation in reference gene expression across biological replicates and populations. Our introduction of the correlation matrix and PCA for the selection of appropriate reference genes provides a statistical method for identifying complementary reference genes. This allowed us to select reference genes that were individually stable and complementary.
It is also important to evaluate gene expression stability across populations or treatments (Lu et al. 2013, Mamidala et al. 2011, Zhai et al. 2014). Both NormFinder and BestKeeper, two of the most commonly used programs for identifying best reference genes, do not clearly address this issue, because they only consider stability for a single source of variation. Our data were highly structured, which allowed us to partition the observed variation in expression and estimate stability for three sources of variation, technical variation, biological replication, and among populations. This allowed us to generalize BestKeeper and required the best reference genes to be stable for all three sources of variation.
Adequate reference genes have been identified in many insect species (Bansal et al. 2012, Chang et al. 2017, Cristiano et al. 2016, Lu et al. 2013, Ma et al. 2016, Mamidala et al. 2011, Paim et al. 2012, Rodrigues et al. 2014, Van Hiel et al. 2009, Yang et al. 2014a, Zhai et al. 2014). For our work, nine aphids and one psyllid (A. craccivora, A. glycines, A. gossypii, Diaphorina citri, Diuraphis noxia, L. erysimi, Megoura viciae, Myzus persicae, R. padi, and T. citricida [Aphis citricidus]) served as the source of potential reference genes for A. glycines for analysis of detoxification gene expression. From our three-phase screening, the HKGs identified as suitable reference genes for pyrethroid resistance expression were actin, RPL9, RPS9, RPL17, AK, and RNAPol2. Of the 10 most frequently used reference genes in insects (i.e., Actin, RPL, Tubulin, GAPDH, RPS, 18S, EF1α, TATA, HSP, and SDHA) (Lü et al. 2018), our work includes three of these genes (Actin, RPL, and RPS). Other studies on aphid species also validated these three genes as robust reference genes as detailed below (Bansal et al. 2012, Kang et al. 2017, Koramutla et al. 2016, Sinha and Smith 2014).
Actin plays many roles in cell function (Perrin and Ervasti 2010) and was identified to have good expression stability in our study. Actin also showed good stability in other aphid species such as Diuraphis noxia under exposure to host plant resistance (Sinha and Smith 2014), R. padi with viral infection (Wu et al. 2014), and Myzus persicae across different tissues, wing dimorphism, photoperiod, and temperature (Kang et al. 2017). Additionally, actin was validated in other insects under different experimental conditions such as in Apis mellifera L. subjected to bacterial infection (Scharlaken et al. 2008), Drosophila melanogaster Meigen under virus challenge, heat stress, and various diets (Ponton et al. 2011), Anopheles sinensis Wiedemann in different life stages, tissues, and levels of insecticide (neonicotinoid) resistance, and Liriomyza trifolii (Burgess) in different development stages and temperatures (Chang et al. 2017). Jiang et al. (2010) and Yang et al. (2014b) also found actin to be suitable as a reference gene for Liposcelis bostsrychophila Badonnel under pyrethroid stress and immatures of Locusta migratoria (L.). However, actin was found to have low stability under different biotic and abiotic conditions in Drosophila suzukii (Matsumura) (Zhai et al. 2014), suggesting that even commonly used HKGs must be validated on a case-by-case basis.
Ribosomal proteins are involved in multiple processes in the genome (Plocik and Guthrie 2012), and translation is one of the most important (Hoffman et al. 1996). We found three ribosomal genes with good stability. The gene RPL9 had the best stability across populations of A. glycines with varying levels of insecticide resistance. It has been reported as a good reference gene for different nymphal stages and tissues of other insects (An et al. 2016, He et al. 2014). The gene RPS9 has been established as a suitable reference gene in A. glycines and Diabrotica virgifera virgifera LeConte subjected to host plant resistance and exposure to dsRNA (Bansal et al. 2012, Rodrigues et al. 2014). Our results demonstrate the stability of this gene across populations of A. glycines with a range of pyrethroid resistance. The gene RPL17 was established as a good reference gene in our work. This gene also showed stable expression in T. citricida undergoing temperature variation (Ma et al. 2016) and Myzus persicae under different photoperiods, temperatures, and levels of insecticide susceptibility (Kang et al. 2017).
The enzyme AK catalyzes phosphorylation in cells (Zhou et al. 1998) and was stable in our experiment. Expression of this gene was found to be stable across different body tissues in Bombus terrestris (L.) (Horňáková et al. 2010); insecticide-induced stress in Spodoptera litura (Fabricius) (Lu et al. 2013); different populations, developmental stages, and tissue in Drosophila suzukii (Zhai et al. 2014).
The multiprotein RNAPol2 transcribes DNA into mRNA (Hirose et al. 1999). RNAPol2 showed the highest levels of stability across populations of A. glycines in our study. The expression of RNAPol2 was also found to be stable under starvation and UV irradiation stress treatments in other insects (Shang et al. 2015).
Interestingly, some genes (e.g., TBP, EF1α, GAPDH, β-Actin, and RPS18) that are widely used in studies related to pyrethroid resistance were among the least stable across our populations of A. glycines with different level of pyrethroid resistance. The genes TBP and EF1α have been used as reference genes in various insect studies, including Schistocerca gregaria (Forskal) (Van Hiel et al. 2009), Drosophila melanogaster (Ponton et al. 2011), A. glycines (Bansal et al. 2012), Plutella xylostella (L.) (Fu et al. 2013), T. citricida (Shang et al. 2015), and A. gossypii (Ma et al. 2016). In the aforementioned studies, TBP and EF1α were found to be stable across different experimental conditions. Additionally, EF1α was found to be stable between populations of Bemisia tabaci (Gennadius) varying in neonicotinoid resistance (Li et al. 2013) and between groups of Ctenocephalides felis (Bouché) exposed and not exposed to avermectin (McIntosh et al. 2016). However, we did not find that these genes were expressed consistently across A. glycines populations in this study. GAPDH, which is involved in energy production (Nicholls et al. 2012), is commonly used in insect species as a reference gene under disease stress, different tissues, temperature fluctuations, and insecticide exposure (avermectin) (Chang et al. 2017, Ma et al. 2016, McIntosh et al. 2016, Paim et al. 2012, Scharlaken et al. 2008, Van Hiel et al. 2009). However, GAPDH did not show adequate stability for A. glycines populations with different levels of pyrethroid resistance. The genes β-Actin and RPS18 have been used several times as reference gene in insects, being mostly successful in Rhodnius prolixus Stål in experiments with different tissues (Paim et al. 2012), T. citricida under heat stress and across various life stages (Shang et al. 2015), A. gossypii undergoing different developmental studies, geographical distribution and feeding conditions for β-Actin, and temperature oscillation for RPS18 (Ma et al. 2016), and Tetranychus cinnabarinus (Boisduval) strains with different levels of resistance to insecticides in different developmental stages for RPS18 (Sun et al. 2010). Although at the initial screening stages of our study β-Actin and RPS18 showed good stability characteristics when paired with other genes, their stability was not good by themselves. Then the stability of these two genes decreased in the second level of screening; therefore, they were not chosen to be final candidate reference genes. These results demonstrated that even commonly used reference genes need to be validated in specific species and experimental condition combinations.
This is the first study that sought to assess stability of candidate reference genes across A. glycines populations for examination of expression of detoxification genes related to pyrethroid resistance. To understand if differences in pyrethroid susceptibility are due to overexpression of detoxification genes in important agricultural pests, such as A. glycines, it is imperative to have adequate reference genes to normalize the expression across populations. Furthermore, the work presented here brings three innovations to the process of mining the best reference genes for gene expression analysis. First, a three-phase screening process was used to identify potential genes based on expression and stability, and to validate the identified genes. Second, the commonly used software for selecting reference genes, BestKeeper, generalized to handle nested data structures and independently evaluate sources of variation in expression. Third, PCA was used to identify nonredundant genes. PCA provided an objective method to select a stable set of genes. The methods and results of this robust study should contribute to gene expression studies not only in A. glycines, but also to other aphid species.
We are grateful to Amelia Lindsey for providing a review of an earlier version of this paper, Lisa Behnken and Bruce Potter for helping to locate fields with soybean aphid infestations, and James Menger, Arthur Vieira, Mads Bartz, and Pheylan Anderson for assistance in the field and laboratory. This work was supported by the Minnesota Rapid Agricultural Response Program (RARF) and the USDA National Institute of Food & Agriculture (NIFA) (award no. 2019-68008-29892).