Temperature and Development Impacts on Housekeeping Gene Expression in Cowpea Aphid, Aphis craccivora (Hemiptera: Aphidiae)

Quantitative real-time PCR (qRT-PCR) is a powerful technique to quantify gene expression. To standardize gene expression studies and obtain more accurate qRT-PCR analysis, normalization relative to consistently expressed housekeeping genes (HKGs) is required. In this study, ten candidate HKGs including elongation factor 1 α (EF1A), ribosomal protein L11 (RPL11), ribosomal protein L14 (RPL14), ribosomal protein S8 (RPS8), ribosomal protein S23 (RPS23), NADH-ubiquinone oxidoreductase (NADH), vacuolar-type H+-ATPase (ATPase), heat shock protein 70 (HSP70), 18S ribosomal RNA (18S), and 12S ribosomal RNA (12S) from the cowpea aphid, Aphis craccivora Koch were selected. Four algorithms, geNorm, Normfinder, BestKeeper, and the ΔCt method were employed to evaluate the expression profiles of these HKGs as endogenous controls across different developmental stages and temperature regimes. Based on RefFinder, which integrates all four analytical algorithms to compare and rank the candidate HKGs, RPS8, RPL14, and RPL11 were the three most stable HKGs across different developmental stages and temperature conditions. This study is the first step to establish a standardized qRT-PCR analysis in A. craccivora following the MIQE guideline. Results from this study lay a foundation for the genomics and functional genomics research in this sap-sucking insect pest with substantial economic impact.


Introduction
Quantitative real-time PCR (qRT-PCR) is a powerful technique to quantify gene expressions during different biological processes [1]. Although qRT-PCR is one of the premier research tools, limitations still exist, several factors can influence the threshold cycle values including RNA quality, cDNA concentration, and PCR efficiency [2,3]. The most extensively adopted approach in qRT-PCR analysis is to normalize the expressions of target genes through measuring in parallel the expression of a housekeeping gene (HKG). Housekeeping genes, involved in basic cellular functions, are typically believed to possess inherent stable and constitutive expression in different samples under various biotic and abiotic conditions [1].
The cowpea aphid, Aphis craccivora Koch (Hemiptera, Aphidiae), is an important pest of cowpea, Vigna unguiculata (L.), one of the most important food crops in the semiarid tropical regions, including Asia, Africa, southern Europe, and Central and South America. Aphis craccivora typically feeds on several species of legumes (family Fabaceae) worldwide, including alfalfa, beans, chickpea, lentils, lupins, and peanuts. Aphids can infest cowpea through direct feeding on leaves, pods and other aerial tissues of the plant, or indirectly through the transmission of virus diseases [4][5][6]. A. craccivora can cause great damage even at low population densities because of its ability to transmit at least 14 viruses including the potyviruses, the cowpea aphid-borne mosaic virus and the blackeye cowpea mosaic virus [6,7]. In order to better understanding the molecular basis and facilitate the development of integrated pest management strategies of A. craccivora, Roche 454 pyrosequencing technology was used to generate the transcriptome of A. craccivora [7]. To take advantage of these genomics resources, establishing a standardized qRT-PCR procedure in A. craccivora following the MIQE (Minimum Information for publication of Quantitative real time PCR Experiments) guidelines [8] will be instrumental for the subsequent functional and epi-genomic research.
The objective of this research was to address an important aspect of gene expression studies in A. craccivora. Here, ten candidate HKGs including elongation factor 1 α (EF1A), ribosomal protein L11 (RPL11), ribosomal protein L14 (RPL14), ribosomal protein S8 (RPS8), ribosomal protein S23 (RPS23), NADH-ubiquinone oxidoreductase (NADH), vacuolar-type H+-ATPase (ATPase), heat shock protein 70 (HSP70), 18S ribosomal RNA (18S), and 12S ribosomal RNA (12S) were selected from the publically available A. craccivora transcriptome resources and the sequence obtained from GenBank [7]. The expression profile of these candidate HKGs was investigated across different developmental stages and under various temperature regimes. As a result, a suite of reference genes were recommended for the qRT-PCR analysis in A. craccivora.

Ethics statement
The cowpea aphid, Aphis craccivora Koch (Hemiptera, Aphidiae), was collected from a greenhouse on fava bean, Vicia faba (Fabales, Fabaceae), at the University of Kentucky. Aphis craccivora colony was maintained on seedlings of fava bean in a growth chamber at 23°C with a photoperiod of 12: 12 (L: D) and 50% relative humidity. No specific permit was required for the described collection. A. craccivora is a common aphid species with agricultural importance in the United States.

Samples preparation
For the developmental stage treatment, 10 A. craccivora adults (only unwinged individuals) and 20 nymphs (mixed nymphal stages) were, respectively, placed on fava bean leaves resting on a wet filter paper in a petri dish (9 cm diameter) for 2 d at 22°C. There are six replicates for the adult and nymph stages, respectively; therefore, there were 12 biological samples in total. For the temperature treatment, 10 A. craccivora adults and 20 nymphs (mixed nymphal stages) were, respectively, exposed to 10°C, 22°C, and 30°C, respectively, for 2 d. Each treatment was repeated three times independently, therefore, there were 18 biological samples for the temperature experiment. All the experiments were conducted in a growth chamber with a photoperiod of 14: 10 (L: D) and 50% relative humidity. After treatments, aphids were initially snap frozen in liquid nitrogen in a 1.5 ml microcentrifuge tube and then stored at -80°C for the subsequent total RNA extraction.

Total RNA extraction and cDNA synthesis
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) according to previously described methods [9,10]. First strand cDNA was synthesized from 1 μg of total RNA with M-MLV reverse transcription kit (Invitrogen, USA) using a random N primer according the manufacturer's recommendations.
Reference gene selection and primer design and quantitative real-time PCR A total of 10 housekeeping genes that are commonly used in qRT-PCR analysis were selected as the candidate ( Table 1). Majority of these genes have been previously used as the reference genes in other insect species [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Primers for EF1A was designed based on the sequences obtained from GenBank, and the others were obtained from the transcriptome of A. craccivora [7]. Primers for the qRT-PCR analysis were designed online, https://www.idtdna.com/ Primerquest/Home/Index. The information of qRT-PCR amplifications and programs were described in detail in our previous study [9,10]. The standard curve and PCR efficiency of each candidate were constructed and calculated according to previously described methods [9,10].

Stability of gene expression
All biological replicates were used to calculate the average C t value. The stability of the ten HKGs was evaluated by algorithms geNorm [1], NormFinder [26], BestKeeper [27], and the comparative ΔC t method [28]. Finally, we compared and ranked the tested candidate HKGs based on a web-based analysis tool RefFinder (http://www.leonxie.com/referencegene.php).

Transcriptional profiling of candidate reference genes
The entire candidate HKGs were visualized as a single amplicon with expected size on a 1.5% agarose gel (S1 Fig). Furthermore, gene-specific amplification was confirmed by a single peak in real-time melting-curve analysis (S2 Fig). Standard curves were created for all the candidate HKGs, and the PCR efficiency and correlation coefficient for each standard curve were shown in Table 1.
The mean and the standard deviation (SD) of the C t values were calculated for all the samples (S1 Table). RPL11 (SD = 0.61) had the least variable expression level and it was reflected in its low SD values. On the contrary, EF1A (SD = 1.09) had the most variable expression levels, and it was shown in its high SD values. Additionally, 18S had the lowest C t values (C tavg = 8.50), suggesting that it had the highest expression level, whereas, NADH was the least expressed gene among the candidates (C tavg = 27.82) (Fig 1, S1 Table).

Selection of the best candidate reference genes
Based on geNorm, under the impact of temperature, RPL14 and RPS8 were co-ranked as the most stable genes. The overall order from the most stable to the least stable reference genes was: RPL14 = RPS8, RPL11, RPS23, ATPase, 12S, HSP70, NADH, EF1A, 18S (Table 2). Under the impact of development, RPL14 and RPS8 were co-ranked as the most stable genes. The overall order from the most stable to the least stable reference genes was: RPL14 = RPS8, RPL11, ATPase, RPS23, 12S, HSP70, NADH, EF1A, 18S (Table 3).
According to the ΔC t method, under the impact of temperature, RPS8 was the top-ranked gene. The overall order from the most stable to the least stable reference genes was: RPS8, RPL14, RPL11, ATPase, HSP70, NADH, 12S, RPS23, EF1A, 18S (Table 2, S2 Table). Under the impact of development, RPS8 was also the top-ranked gene. The overall order from the most stable to the least stable reference genes was: RPS8, RPL14, RPL11, ATPase, HSP70, NADH, 12S, RPS23, EF1A, 18S (Table 3, S3 Table). Based on NormFinder, under the impact of temperature, RPS8 was the most reliable and stable reference gene. The overall order from the most stable to the least stable reference genes was: RPS8, RPL11, RPL14, ATPase, 12S, RPS23, HSP70, NADH, EF1A, 18S (Table 2). Under the impact of development, RPS8 was also the top-ranked gene. The overall order from the most stable to the least stable reference genes was: RPS8, RPL14, RPL11, ATPase, HSP70, NADH, 12S, RPS23, 18S, EF1A (Table 3).
According to BestKeeper, the stability of a gene is directly proportional to the [r] value, while it is inversely proportional to the SD value. Under the impact of temperature, RPL11 had the highest [r] value, and RPL14 had the lowest SD value across all the samples ( Table 2, S4  Table). Under the impact of development, RPL14 had the highest [r] value, and HSP70 had the least variable expression levels across all the samples (Table 3, S5 Table) Comprehensive ranking of best reference genes using RefFinder Under the impact of temperature, according to RefFinder, the comprehensive ranking of candidate reference genes from the most to the least stable was: RPS8, RPL14, RPL11, ATPase, RPS23, 12S, HSP70, NADH, EF1A, 18S (Table 2). Under the impact of development, the comprehensive ranking of candidate reference genes from the most to the least stable was: RPS8, RPL14, RPL11, HSP70, ATPase, RPS23, 12S, NADH, EF1A, 18S (Table 3). Interestingly, RPL11, RPS8, and RPL14 were the three most stable HKGs throughout different developmental stages and temperature conditions.

Quantitative analysis of candidate reference genes based on geNorm
To decide the minimal number of genes mandatory for normalization, the V-value was computed by geNorm. geNorm analysis revealed that the pair-wise variation value V6/7 is higher than V5/6 (Fig 2). Increasing variation in this ratio corresponds to decreasing expression stability, due to the inclusion of a relatively unstable sixth gene. Therefore, five genes (PRL14, RPS8, RPL11, ATPase, and RPS23) are necessary for accurate normalization. Including a sixth reference gene has no significant effect on the normalization factor (Fig 2).

Discussion
qRT-PCR quantification demands a comprehensive normalization by housekeeping genes to counteract confounding variations in experimental data. Housekeeping genes have been considered to be expressed in all cell types of the organism at a constant level to maintain basic cellular functions. However, there are no "universal" reference genes that are stably expressed and appropriate for the entire cell and tissue, and all kinds of test conditions [1]. Most gene expression studies in the literature use one single housekeeping gene; this will deeply influence the outcome of the statistical analysis and may bring about inaccurate data interpretation [29]. Therefore, customized reference gene selection under specific experimental conditions is highly recommended [11].
Recently, there is an influx of reference gene selection studies in insects, including convergent lady beetle, Hippodamia convergens; sweet potato whitefly, Bemisia tabaci; diamondback moth, Plutella xylostella; brown planthopper, Nilaparvata lugens; beet armyworm, Spodoptera exigua; oriental leafworm moth, Spodoptera litura; oriental fruit fly, Bactrocera dorsalis; Colorado potato beetle, Leptinotarsa decemlineata; soybean aphid, Aphis glycines; Russian wheat aphid, Diuraphis noxia; bird cherry-oat aphid, Rhopalosiphum padi; pea aphid, Acyrthosiphon pisum; bumblebees, Bombus terrestris and Bombus lucorum; western flower thrips, Frankliniella occidentalis; and honeybee, Apis mellifera [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Here, the expression profiles of ten HKGs from A. craccivora were evaluated across different developmental stages and temperature conditions. Our results are largely consistent with previous studies. For example, RPS8 (the component of the 40S ribosomal subunit) and RPL14 (60S ribosomal subunit) were the most stable HKGs across different developmental stage and temperature conditions, whereas the expression of 18S varied under the two conditions [14,16]. Not surprisingly, the comprehensive rankings (RefFinder) of these candidate reference genes under the two experimental conditions were, in principal, comparable to the rankings complied by the four algorithms, geNorm, Normfinder, BestKeeper, and the ΔC t method, respectively (Tables 2 and 3). Based on the comprehensive analyses, RPS8, RPL14, and RPL11 were the most stable A. craccivora HKGs under different developmental stages and temperature conditions. There has been ongoing discussion about the optimal number of reference genes warrant for qRT-PCR analysis [9,14]. To prevent biased normalization, multiple instead of a single reference gene have been gradually adopted to normalize the expression of target genes under test conditions [30]. Our results showed that the pair-wise variation value of V6/7 is higher than that of V5/6 (Fig 2), suggesting that five reference genes are warranted for the accurate normalization in A. craccivora under different developmental stages and temperature conditions. A phloem-feeding cowpea aphid, A. craccivora, is one of the key pests of cowpea, a major protein source for people in West Africa. Most recently, Roche 454-based pyrosequencing generated 176,262 raw reads from an A. craccivora transcriptome, and de novo assembly produced 7,647 transcripts [7]. Building on this newly developed genomic resource, we carried out the first reference gene selection study in one of the major pest species of cowpea. Although studies involving different developmental stages and /or temperature regimes have been limited [31][32][33][34], the advent of the Genomics Era will facilitate our understanding of A. craccivora, and eventually will lead to the development of integrated pest management strategies. Therefore, this study not only sheds light on establishing a standardized qRT-PCR procedure for quantification of gene expression in A. craccivora, but also lays a solid foundation for the genomics and functional genomics research in this sap-sucking insect pest.