Reference gene selection to determine differences in mitochondrial gene expressions in phosphine-susceptible and phosphine-resistant strains of Cryptolestes ferrugineus, using qRT-PCR

Cryptolestes ferrugineus is a serious pest of stored grain and has developed high levels of resistance to phosphine fumigants in many countries. Measuring differences in expression levels of certain ‘resistant’ genes by quantitative real-time PCR (qRT-PCR) may provide insights into molecular mechanisms underlying resistance to phosphine in C. ferrugineus, but reliable qRT-PCR results depend on suitable reference genes (RGs). We evaluated the stability of nine candidate RGs across different developmental stages and phosphine strains of C. ferrugineus, using four softwares. The results showed that RPS13 and EF1α were the most stable RGs, whereas α-TUB was the least under developmental stages. Across the different strains, RPS13 and γ-TUB were the most stable RGs, whereas CycA and GAPDH were the least. We confirmed the reliability of the selected RGs by qRT-PCR analyses of the mitochondrial cox1 gene. Expression of cox1 was not significantly different in the phosphine-resistant strain compared with the phosphine-susceptible strain, but three mitochondrial genes (nad3, atp6 and cob) were significantly down-regulated. These results suggest that alterations in the expressions of these three genes may be associated with phosphine resistance in C. ferrugineus. The findings will facilitate future functional genomics studies on the development and phosphine resistance in C. ferrugineus.

The earliest reports of resistance to phosphine in the field were made in 1972-1973 7 , but since that time, phosphine resistance has been reported in a whole host of the most common stored-product insect pests, including Tribolium castaneum [8][9][10][11][12] , Rhyzopertha dominica 8,9,11,12 , Sitophilus oryzae [13][14][15] , Lasioderma serricorne 16 , Sitophilus zeamais 17 , Oryzaephilus surinamensis 12,18 , C. ferrugineus 4,5 and Liposcelis bostrychophila 19 . Despite several decades of study, the biological mechanisms underlying resistance to phosphine and its toxicity are not well understood 20 . Mitochondria are broadly recognised as the intracelluluar site of biological action of phosphine, where it inhibits aerobic respiration 21 . Phosphine resistance has consistently been associated with carbon dioxide production, metabolic rate, respiration rate and even with walking activity 12,18,[22][23][24] . Mitochondrial DNA is known to encode a limited number (<20) of the polypeptide components of respiratory complexes I, III, IV, and V. Genes controlling expressions of the components of complex II are conspicuously lacking in the mitochondrial genomes that have so far been characterized 25 , whilst the following genes of the mitochondrial genome are known to encode the subunits of respiratory complexes I, III, IV and V, respectively: NADH dehydrogenase subunit 3 (nad3), cytochrome b (cob), cytochrome c oxidase subunit I (cox1) and ATP synthase F0 subunit 6 (atp6). To date, we have used next generation gene sequencing techniques to establish a C. ferrugineus transcriptome database (GenBank accession number SRA245468) to identify potential molecular targets for the control of C. ferrugineus. The sequences in this database can be used in gene expression analyses to reveal the molecular mechanisms underlying biological processes such as insecticide resistance. However, despite a number of transcriptome analyses related to insecticide resistance having been performed with the aim of detecting and managing resistant insect strains, very limited information is currently available for the expression profiles of target genes associated with phosphine resistance in C. ferrugineus.
Quantitative real-time PCR (qRT-PCR) is a valuable tool for gene expression studies due to its high sensitivity, specificity and convenience in high throughput analyses 26 . In general, when the technique is used to analyze gene expressions, reference genes (RGs) are needed to normalize nonspecific variations or errors that can be caused by sample quantity, variations in efficiency of RNA extraction, cDNA concentration, primer performance, PCR efficiency and experimental precision 27,28 . The optimal number of RGs is another normalization factor because one single reference gene is not enough to normalize gene expression data 28,29 . The ideal RGs used in quantitative gene expression studies should be stable under the experimental conditions employed 27 . In the last two decades, several housekeeping genes, such as alpha-tubulin (α-TUB), beta actin (β-ACT), elongation factor 1 alpha (EF1α), glyceraldehydes 3 phosphate dehydrogenase (GAPDH), 18 S ribosomal RNA (18 S), ribosomal protein S13 (RPS13) and ribosomal protein L13a (RPL13a), have been widely used as RGs because they are considered to have a stable and uniform expression pattern as a result of their basic, ubiquitous, cellular functions and consistencyof expressions under different experimental conditions 26,30 . However, it has now been shown that many of the conventional RGs (e.g. 18 S, β-ACT and GAPDH) exhibit variable expression levels depending on the organism, its developmental stage, its response to external stimuli and the experimental conditions 27,[30][31][32] . Thus, if the RGs are not selected appropriately, it could result in inaccurate calculation of the normalization factor, which could obscure the detection of actual biological differences between samples. Therefore, it is necessary to validate the expression stability of RGs under specific experimental conditions before using them for normalization. To date, many studies have been conducted in various insect groups to select suitable RGs under different experimental conditions, such as Lepidoptera 29 , Hymenoptera 33 , Diptera 34 , Corrodentia 35 , Coleoptera 28 , Homoptera 30, 36, 37 , Thysanoptera 38 and Hemiptera 31 . However, none of these previous studies has focussed on C. ferrugineus, which limits the current opportunities for functional genomics studies on this important pest.
In this study, the expression stability of nine candidate RGs, namely succinate dehydrogenase complex subunit A (SDA), cyclin A (CycA), γ-tubulin (γ-TUB), α-TUB, EF1α, GAPDH, RPL13a, RPS13 and β-ACT, was evaluated in different developmental stages and phosphine resistant strains of C. ferrugineus using four statistical algorithm methods (geNorm 39 , NormFinder 40 , BestKeeper 41 and delta Ct 42 ). The optimal number of RGs (normalization factors, NFs) was also determined. The reliability of the selected RGs was confirmed by qRT-PCR analyses of the mitochondrial cox1 gene. The expression profiles of three other mitochondrial genes (nad3, atp6 and cob) were aso analyzed by qRT-PCR for different phosphine resistant strains of C. ferrugineus. Previous studies have shown the importance of these mitochondrial genes in phospine resistance 43 . The present study sought to identify reliable RGs that could be used to accurately measure gene expressions under different developmental stages and phosphine resistant strains, which could therefore benefit future functional genomics studies of the development and phosphine resistance in C. ferrugineus.

Results
Amplification efficiency and expression profiles of the candidate RGs and target genes. The primer sets of the nine candidate RGs and four target genes for qRT-PCR used in the study were shown in Table 1. The primer specificities were validated by the presence of a single band of the expected molecular weight on a 3% agarose gel ( Supplementary Fig. S1) and a sharp peak in the melting curve analysis ( Supplementary Fig. S2). A standard curve was generated for each gene using five-fold serial dilutions of cDNA ( Supplementary Fig. S3). The amplification efficiencies (E%) of the primer pairs ranged from 94.71% for SDA to 110.06% for cox1, with the correlation coefficient (R 2 ) values varying from 0.993 to 0.999. On this basis, all primers were considered suitable for use in the accurate amplification of the corresponding genes ( Table 1).
The expression profiles of the nine candidate RGs and four target genes were analyzed by calculating the cycle threshold (Ct) value, which represents the number of cycles required for the amplification to exceed a fixed threshold. As shown in Fig. 1 Expression stability of the candidate RGs in different developmental stages. We calculated and ranked the expression stability of nine candidate RGs across the different developmental stages of the insect using four statistical algorithms. The stability varied considerably between the nine candidates according to all four algorithms ( Table 2). The most stable RGs recommended by geNorm were largely similar with those identified by the delta Ct method, showing that RPS13, RPL13a and EF1α were the most stable RGs. The three most stable RGs determined by NormFinder were SDA, GAPDH and EF1α, whereas γ-TUB, β-ACT and RPS13 were judged to be the most stable by BestKeeper (Table 2). Thus, the different algorithms did not all predict the same RGs as the most stable. However, we found that RPS13 and EF1α were ranked among the top three most stable RGs by at least three algorithms. Interestingly, all four algorithm programs consistently identified α-TUB as the least stable RG. The geNorm analysis showed that the pairwise variation value of V 2/3 (0.106) was below the cut-off  value of 0.15 ( Fig. 2A), indicating that two RGs are sufficient for accurate normalization across developmental stages. Therefore, it was considered that the two most stable RGs (RPS13 and EF1α) could be normalized across the developmental stages.

Expression stability of candidate RGs in different phosphine strains.
Among the different insect strains tested, the expression stability rankings generated by three of the four algorithms (geNorm, NormFinder, and delta Ct) were highly similar: RPS13 was identified as the highest expression stability among all the RGs,  Table 2. Ranking of RGs based on their expression stability across different developmental stages of C. ferrugineus according to geNorm, NormFinder, BestKeeper, and delta Ct method.  Table 3. Ranking of RGs based on their expression stability in different strains according to geNorm, NormFinder, BestKeeper, and delta Ct algorithm methods.
while CycA was the least stable gene ( Table 3). The BestKeeper method identified RPS13 as the most stable of the RGs and CycA as one of the least stable genes. γ-TUB was identified as the most stable RG by BestKeeper. This gene was also one of the genes with the highest stability according to the NormFinder and delta Ct methods. Taking account of the predictions of all four algorithms, RPS13 and γ-TUB were the most stable RGs. The geNorm analysis revealed that the pairwise variation values of V 2/3 to V 8/9 were all below the cut-off value of 0.15 (Fig. 2B). Thus, RPS13 and γ-TUB were the most suitable combination of RGs to use for measuring gene expression in different C. ferrugineus strains.
Validation of the selected RGs in C. ferrugineus. To demonstrate the effect of the selected RGs on target gene expressions, the relative expression levels of cox1 were analyzed using different sets of RGs and compared between the four insect strains, i.e. phosphine-susceptible (SS0), -resistant strains (RS0) and each of these two exposed to phosphine for 4 h (SS4 and RS4) ( Table 4). Similar expression profiles were obtained when the expression levels were normalized using the most stable RG (RPS13) and the combination of the two most stable RGs (RPS13 + γ-TUB), showing that the transcriptional levels of cox1 were not significantly different among the different strains (Fig. 3). However, when the expression levels were normalized with the least suitable RG (CycA), the expression levels of cox1 were higher in three strains (SS4, RS0 and RS4) than that in the susceptible strain   (SS0). Similarly, the expression levels normalized using another less suitable RG (GAPDH) were also higher in the SS4 and RS4 strains than that in the susceptible strain (SS0). When the cox1 transcript expression levels were normalized by the combination of the two most suitable RGs, they were significantly different from those calculated with the least suitable RG (CycA or GAPDH), in both the SS4 and RS4 strains (P < 0.01) (Fig. 3).
Expression profiles of four mitochondrial genes in the different C. ferrugineus strains. The relative expression profiles of four mitochondrial genes (nad3, atp6, cob and cox1) were assessed in the different C. ferrugineus strains using a combination of RPS13 and γ-TUB as the RGs for the normalization of the data (Fig. 4).
The results showed that the expression levels of three mitochondrial genes (nad3, atp6 and cob) were markedly lower in the phosphine-resistant strains compared to the phosphine-susceptible strains, by approximately 1479, 753, and 5756 fold, respectively ( Fig. 4A-C). When exposed to phosphine, the transcriptional expression levels of these three genes were significantly lower than those in the unexposed strains, regardless of the duration of exposure (P < 0.01), and there was an overall decreasing trend in expression levels with increasing duration of exposure in all strains. However, the extent of this decline for any given exposure duration was different between the susceptible and resistant strains, especially for the shorter exposure times (0.5 h, 1 h and 2 h). After 0.5 h exposure to phosphine, the expression levels of cob had decreased from 1 (RS0) to 0.32 (RS0.5) in the resistant strain, but had decreased by only 0.55 times (SS0.5/SS0) in the susceptible strain (Fig. 4C). Clearly, the resistant strains responded faster to conditions of phosphine stressors than the susceptible strains. For cox1 under the same exposure duration, the expression levels changed from 1.07-(SS8) to 1.27-times (RS8) in the resistant strains, but there was no significant difference among any of all the strains in this respect (P > 0.05) (Fig. 4D).

Discussion
In this study, we performed a comprehensive analysis of nine RGs in C. ferrugineus to aid in the normalization of gene expression analyses by qRT-PCR since qRT-PCR has proven effective for exploring the underlying mechanisms of insecticide resistance in many insect pests [44][45][46] . We found that the most stable RGs across different developmental stages differed among different algorithms, whereas α-TUB was consistently identified as the least stable RG by all four algorithms. SDA and GAPDH were found to be the most stable RGs by NormFinder; however, they were not determined as stable by the other three algorithms. For different insect strains, slight differences in the stability ranking were observed among the four algorithms. RPS13, γ-TUB, β-ACT and EF1α were identified as the top four most stable RGs across strains, whereas CycA and GAPDH were ranked as the two least stable RGs, by all four algorithms. These variations have also been obtained in many previous studies of insects, such as Lipaphis erysimi 30 , Galeruca daurica 28 , Aphis gossypii 37 and Helicoverpa armigera 29 , as different algorithms The asterisks indicate a statistically significant difference compared with the control, by the Student's t-test (SS0 and RS0 were the controls, respectively; *P < 0.05, **P < 0.01).
and analytical procedures are used in the four algorithm programmes and these do not seem to affect the overall validation quality 47 . RPS13 has been identified as the most stable gene in different tissues of Sesamia inferens 48 , and RPS13 mRNA expression was shown to be unaffected by alcohol administration, SIV infection, or antiviral therapy in rhesus macaques 49 . Our results demonstrated that RPS13 was the most stable gene both among different developmental stages and among different strains of C. ferrugineus, indicating that RPS13 exhibits a wide degree of stability in a range of conditions. EF1α is a commonly used housekeeping gene in qRT-PCR experiments 50 . Previously, EF1α has been reported to exhibit a stable level of expression in A. gossypii under biotic and abiotic stressors 37 . Being an important specific protein factor involved in the process of protein translation, it was not suprising that EF1α was identified in the present study as one of the most stable RGs in terms of its expression across the developmental stages of C. ferrugineus. However, it was not the most stable RG among different phosphine strains of C. ferrugineus. γ-TUB, a gene that has not previously been identified as a stable RG, was determined to be one of the most stable RGs across different strains in the present study, but it was not the most appropriate RG in terms of its stability across different developmental stages in C. ferrugineus.
GAPDH has previously been confirmed to be one of the most suitable RGs in a number of other studies [35][36][37] , although low GAPDH expression stablitity has been observed in a few studies under certain conditions, such as in different tissues, developmental stages and treatments of Agrilus planipennis 50 . In another study, GAPDH was found to be the least stable gene in phosphine-induced strains of T. castaneum 51 . In the present study, we identified GAPDH as an unsuitable RG across different strains in C. ferrugineus 51 . RPL13a was calculated as one of the most stable RGs to phosphine induction in T. castaneum 51 , but it did not exhibit such a high degree of stability across different phosphine strains of C. ferrugineus in the present study, indicating that the stability of the same RG may be species-specific. α-TUB has also been widely used as a RG in gene expression studies in G. daurica 28 and Sogatella furcifera 36 . However, other researcher have revealed that α-TUB does not satisfy certain basic requirements for use as an internal control 37 . In the present study, the results showed that the expression levels of α-TUB were not stable across different developmental stages of C. ferrugineus, therefore this gene would not be the most reliable internal control for comparative gene expression analyses. These differences may partly be explained by the fact that α/β-TUB not only acts as one of the major components of the cytoskeletal structure of cells, but also participates in other cellular functions 52 . The structural protein β-ACT is expressed at moderately abundant levels in most cell types and has been considered an ideal RG for measuring gene expressions in many insects, including Drosophila melanogaster 34 and Apis mellifera 53 . However, β-ACT was not the most stable gene in the present study, which is consistent with some previous studies on A. planipennis 50 and Frankiniella occidentails 32 . These differences highlighted the need to validate RGs for different species, even under similar experimental conditions. Respiratory complexes I, III and IV constitute the main electron transport chain in aerobic organisms and generate the proton gradient in the intermembrane spaces of mitochondria. Complex V, or ATP synthase, uses the energy stored in a proton gradient across a membrane to drive the synthesis of ATP from ADP and phosphate (Pi) 54 . In the present study, the transcriptional profiles of four mitochondrial genes (nad3, cob, cox1 and atp6), encoding the subunits of complexes I, III, IV and V, respectively), were different in various phosphine strains of C. ferrugineus, indicating that these subunits may have different roles in protecting against phosphine stressors. No significant expression differences were observed in cox1 among the phosphine strains in the present study even though Complex IV has previously been proposed as the principal mitochondrial site of action of phosphine 20,22,55 . However, a study on cox1 and cox2 expressions in S. zeamais found no significant association with phosphine susceptibility, in agreement with our findings 56 . Our results suggested that modifications in the expressions of the cox1 gene encoding complex IV may be not associated with phosphine resistance in C. ferrugineus. In contrast, three other mitochondrial genes (nad3, cob and atp6) were significantly suppressed in C. ferrugineus strains exposed to phosphine, and showed extremely low expressions in the phosphine-resistant strain compared to the phosphine-susceptible strain. This could indicate that phosphine may inhibit the insect respiratory system and that a lower respiratory rate may be an important mechanism for phosphine resistance in C. ferrugineus. Our results in this respect are largely in agreement with those of previous studies conducted in T. castaneum 12 , R. dominica 12 , O. surinamensis 12 , S. zeamais 55 , Musca domestica 22 and Sitophilus granaries 22 . Previous studies also showed that the suppression of mitochondrial respiratory chain genes caused an increase in phosphine resistance in Caenorhabditis elegans, and the strongest resistance towards phosphine was identified when the complex III genes were suppressed 57 . These results could explain our observations regarding the magnitude of the decrease in the expression of cob (5756 fold), compared to that of nad3 (1479 fold) and atp6 (753 fold), suggesting that complex III may play a more important role in phosphine resistance in C. ferrugineus. Additionally, we found that when C. ferrugineus were exposed to phosphine, the transcriptional expression levels of nad3, atp6 and cob were altered more significantly in the phosphine-resistant strain than in the phosphine-susceptible strain, indicating a difference in the response to phosphine between the two strains. A compensation mechanism of ATP synthesis should exist in phosphine resistant insects, considering that no fitness cost has been reported in phosphine-resistant strains of the red flour beetle (T. castaneum) 10,58 .
Overall, our results have shown that expression levels in the mitochondrial genes encoded by mitochondria were negatively related to phosphine resistance in C. ferrugineus, supporting the existing hypothesis that phosphine directly targets the mitochondria. There were marked differences in the expression levels of mitochondrial genes between the susceptible and resistant strains of C. ferrugineus.

Methods
Insect sampling and rearing. Insect rearing. Phosphine-susceptible strain of C. ferrugineus was originally collected from farmer's barns in Zhangjiagang County, Suzhou City, Jiangsu Province, China, in July 2013 (LC 50 = 0.035 mg/L, SS). Phosphine-resistant strain was originally collected from the State Grain Reserve Depot in Taicang County, Suzhou City, Jiangsu Province, China, in July 2013 (LC 50 = 18.76 mg/L, RS). Both strains have been maintained in our laboratory since collection, without any special methods to maintain the resistance levels. When performing the present study in 2016, we also measured the LC 50 values for these two strains and found that the LC 50 values almost unchanged compared to that in 2013. All the insects used in this study were reared on whole wheat flour and yeast powder (10:1) in a growth chamber under conditions of 32 ± 1 °C, and 70-80% relative humidity, in total darkness. Experimental treatments. Developmental stages. Five different developmental stages, first instar larvae, fourth instar larvae, pupae, newly emerged adults and aged adults were collected from the susceptible strain ( Table 4). Hundreds of adults (of both sexes) were released on to a culture dish containing a small amount of feed, allowed to lay eggs and removed after 48 h. First instar larvae, fourth instar larvae and pupae were collected after approximately 6, 20 and 26 days, respectively. A proportion of the pupae were transferred to a new culture dish, from which newly emerged adults were collected after approximately 8 days, and aged adults were collected after approximately 40 days. Enough insect samples were collected for three biological replicates and these were immediately frozen in liquid nitrogen and then stored at −80 °C for RNA extraction.
Phosphine-induced stress. Two-week old adults of both susceptible and resistant strains of C. ferrugineus were exposed to phosphine (0.026 mg/L, LC 30 of the susceptible strain) in sealed gas-tight glass jars with a 1 L volume at 30 °C. The control groups were treated with an equal volume of air. The surviving insects were collected after 0.5, 1, 2, 4, 8, 12 and 48 h, respectively (Table 4). Enough samples were collected to provide three biological replicates and these were immediately frozen in liquid nitrogen and stored at −80 °C until used.
Total RNA isolation and cDNA synthesis. Total RNA was isolated by the RNeasy Plus Micro Kit (Qiagen, Germany) following the manufacturer's instructions. Genomic DNA was removed using a genomic DNA elimination column supplied with the kit. The purity and concentration of total RNA were assessed with a BioMate 3 S (Thermo Scientific, USA). RNA samples with an OD 260 /OD 280 ratio ranging from 1.8 to 2.2 and an OD 260 / OD 230 ratio >2.0 were stored at −80 °C for further processing. First-strand cDNA synthesis was carried out with 1 μg total RNA in a 20 μl reaction volume using a PrimeScript Tm II 1st Strand cDNA Synthesis Kit (Takara, Japan). cDNA was stored at −20 °C until use. Each cDNA sample was diluted 10 times with nuclease free water before being used as a template for qRT-PCR.
Selection of reference genes and primer design. We selected nine candidate RGs to investigate their stability as internal controls for qRT-PCR (Table 1). The sequences of the selected RGs were obtained from the C. ferrugineus transcriptome sequence (accession number: SRA245468). Four target genes (nad3, atp6, cob and cox1), obtained from the mitochondrial genome of C. ferrugineus (unpublished data), were selected to determine their phosphine resistance status. The primers used in qRT-PCR were designed using the online primer design software Primer 3.0 (http://frodo.wi.mit.edu/primer3/). Quantitative real-time PCR and amplification efficiency test. qRT-PCR was performed using an ABI 7500 PCR system (Applied Biosystems, USA). The 25 μl reactions contained 1 μl of cDNA template, 12.5 μl of QuantiFast ® SYBR ® Green PCR Master Mix (Qiagen, Germany) and 0.5 μl of each primer. Nuclease-free water was added to make the volume up to 25 μl. Thermal cycling conditions were 95 °C for 5 min followed by 40 cycles of 95 °C for 10 sec followed by 60 °C for 30 sec. At the end of each PCR run, a melting curve analysis from 55 to 95 °C was applied to all reactions to ensure the specificity of the amplified product. Each qRT-PCR reaction was performed in triplicate (technical replicates) on three individual samples (biological replicates). Standard curves were constructed to determine PCR efficiency (E%) and correlation co-efficient (R 2 ) based on a 5-fold dilution series of cDNA. The PCR efficiency (E%) was calculated according to the equation: E = (10 [−1/slope] − 1) × 100 59 .
Expression stability analysis. All biological replicates were used to calculate the average Ct value. Stability values of the nine candidate RGs were assessed using four commonly used software tools: geNorm version 3.5 39 , NormFinder version 0.953 40 , BestKeeper 41 and delta Ct 42 , which utilize different algorithms to calculate the stability in expression levels across samples.
The geNorm program estimates the stability of each gene by calculating an expression stability value "M" and ranks them in an order for a given set of samples. Values of M that are less than 1.5 are considered to identify a stable level of gene expression across treatments. The lower the M value, the more stable the expression. Furthermore, geNorm provides an algorithm to calculate the optimal number of RGs required for better normalization, because use of a single RG may result in inaccurate normalization, leading to significant errors in gene expression analyses 39 . The pairwise variation (V) of these genes is compared with that of others, and the pairwise variation values (V n /V n+1 ) between two sequential normalization factors were used to determine the optimal number of RGs. A pairwise variation value of less than 0.15 implies that the addition of another RG will not make any significant contribution to the accuracy of normalization 39 .
NormFinder calculates the gene expression stability values by taking into account variations between and within defined sample groups or treatments, and provides a gene rank order depending on the values obtained. The candidate gene with the lowest value is considered to be the most stable RG 40 . BestKeeper determines the best-suited RGs based on the geometric mean of the Ct values and PCR amplification efficiency, and combines them into an index 41 . The comparative delta Ct method compares basic Ct values and the relative expression of 'gene pairs' within each sample, to confidently identify reliable RGs 42 .
Validation of reference gene selection. To evaluate the validity of the optimized selection of RGs, the expression levels of one target gene (cox1) were analyzed in 4 different C. ferrugineus strains, i.e. phosphine-susceptible (SS0), phosphine-resistant (RS0) and each of these two exposed to phosphine for 4 h (SS4 and RS4). The least stable genes were CycA (as determined by geNorm, NormFinder and delta Ct method) and GAPDH (as determined by BestKeeper). The most stable gene was RPS13 (as determined by geNorm, NormFinder and delta Ct method) and the combined set of RGs RPS13 + γ-TUB (as calculated and recommended from the findings of this study) were used for comparative analyses. The relative expression levels were calculated according to the 2 −ΔΔCT method 60 . The target gene expression was normalized by only one RG (CycA, GAPDH or RPS13), and the recommended combination of RGs were calculated using a Student's t-test, in SPSS 17.0 software, with significance levels denoted by *(0.01 < P < 0.05) and **(P < 0.01). All experiments were performed in triplicate, and the results were expressed as means ± standard deviation (SD).
Determination of the expression profiles of selected mitochondrial genes. The expression profiles of nad3, atp6, cob and cox1 were calculated in different strains, comprising both susceptible and resistant strains exposed to phosphine for 0.5, 1, 2, 4, 8, 12, 48 h in a concentration of 0.026 mg/L (LC 30 of susceptible strain), with the combined set of RGs (RPS13 + γ-TUB). The relative expression of the four target genes was calculated using the 2 −ΔΔCt method 60 . All the experiments were performed in triplicate, and the results were expressed as means ± SD. The target gene expression of phosphine exposed and unexposed strains were calculated using a Student's t-test, in SPSS 17.0 software, with significance levels denoted by * (0.01 < P < 0.05) and ** (P < 0.01).

Conclusions
This study evaluated the stability of expression of nine RGs for normalization in qRT-PCR of C. ferrugineus using four algorithms, and then explored whether the differences in expression of four mitochondrial genes were associated with phosphine resistance in this pest. Our results showed that RPS13 and EF1α were the most stable RGs across different C. ferrugineus stages, whereas RPS13 and γ-TUB were the most stable across the different insect strains. Expression changes in mitochondrial genes may be an important mechanism of phosphine resistance in C. ferrugineus. Our results will benefit future functional genomics studies of the development and phosphine resistance in C. ferrugineus.