Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 7047.
Published online 2017 August 1. doi:  10.1038/s41598-017-07430-2
PMCID: PMC5539111

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 rusty grain beetle Cryptolestes ferrugineus (Coleoptera: Laemophloeidae) is a serious pest of stored cereal grains and processed commodities throughout the world1, 2. In China, C. ferrugineus is considered to be the most important pest of the Cryptolestes genus and has caused serious economic losses in several southern provinces in recent years3. Currently, phosphine (PH3) is the most widely used fumigant for controlling this species in grain stores worldwide4. With the continuing phase-out of methyl bromide because of its contribution to ozone depletion, it is likely that phosphine will continue to play an important role in the protection of stored products for the foreseeable future5. However, phosphine resistance occurs worldwide and is a major challenge to its continued effective use46. Recent studies from a variety of different countries have indicated that C. ferrugineus is one of the most highly resistant insect pests to phosphine4, 5. However, the molecular mechanisms of this resistance are currently poorly understood.

The earliest reports of resistance to phosphine in the field were made in 1972-19737, but since that time, phosphine resistance has been reported in a whole host of the most common stored-product insect pests, including Tribolium castaneum 812, Rhyzopertha dominica 8, 9, 11, 12, Sitophilus oryzae 1315, 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 understood20. Mitochondria are broadly recognised as the intracelluluar site of biological action of phosphine, where it inhibits aerobic respiration21. Phosphine resistance has consistently been associated with carbon dioxide production, metabolic rate, respiration rate and even with walking activity12, 18, 2224. 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 characterized25, 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 analyses26. 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 precision27, 28. The optimal number of RGs is another normalization factor because one single reference gene is not enough to normalize gene expression data28, 29. The ideal RGs used in quantitative gene expression studies should be stable under the experimental conditions employed27. 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 conditions26, 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, 3032. 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 Lepidoptera29, Hymenoptera33, Diptera34, Corrodentia35, Coleoptera28, Homoptera30, 36, 37, Thysanoptera38 and Hemiptera31. 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 (geNorm39, NormFinder40, BestKeeper41 and delta Ct42). 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 resistance43. 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.


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).

Table 1
Gene primer sequences and amplicon characteristics.

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, the distributions of the mean Ct values for each gene across all samples varied significantly. The mean Ct values of the 9 candidate RGs ranged from 17.82 to 29.81, with the lowest and highest Ct values being obtained from β-ACT (15.46) and SDA (32.85), respectively. Among the target genes, cox1 was the most abundant transcript with the lowest mean Ct value (24.96) and cob was the least abundant mRNA with the highest mean Ct value (27.25).

Figure 1
Cycle threshold (Ct) values of the candidate RGs across the experimental samples. The line in the box indicates the median value. The box represents the values between the 25th and 75th percentiles. The maximum and minimum values are represented by up ...

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 V2/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.

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.
Figure 2
Pairwise variation of RGs in C. ferrugineus, according to the geNorm algorithm. Line graphs show the average expression stability value (M) for (A) different developmental stages and (B) different strains of the insect. Bar graphs show the pairwise variation ...

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, 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 V2/3 to V8/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.

Table 3
Ranking of RGs based on their expression stability in different strains according to geNorm, NormFinder, BestKeeper, and delta Ct algorithm methods.

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).

Table 4
Sample collection information.
Figure 3
Validation of the selected RGs in C. ferrugineus. Expression levels of the target gene (cox1) were tested using different normalization factors. The data represent mean values ± SD. Comparison of means is carried out using the ...

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).

Figure 4
Relative expression of 4 mitochondrial genes according to the selected RGs (RPS13 + γ-TUB). (A) nad3, (B) atp6, (C) cob and (D) cox1 in sample sets across the phosphine-susceptible strain (SS) and phosphine-resistant strain (RS) ...


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 pests4446. 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 and analytical procedures are used in the four algorithm programmes and these do not seem to affect the overall validation quality47.

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 macaques49. 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 experiments50. Previously, EF1α has been reported to exhibit a stable level of expression in A. gossypii under biotic and abiotic stressors37. 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 studies3537, 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 control37. 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 functions52. 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 phosphine20, 22, 55. However, a study on cox1 and cox2 expressions in S. zeamais found no significant association with phosphine susceptibility, in agreement with our findings56. 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 suppressed57. 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.


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 (LC50 = 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 (LC50 = 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 LC50 values for these two strains and found that the LC50 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, LC30 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 OD260/OD280 ratio ranging from 1.8 to 2.2 and an OD260/OD230 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 PrimeScriptTmII 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 (

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) × 10059.

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.539, NormFinder version 0.95340, BestKeeper41 and delta Ct42, 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 analyses39. The pairwise variation (V) of these genes is compared with that of others, and the pairwise variation values (Vn/Vn+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 normalization39.

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 RG40. BestKeeper determines the best-suited RGs based on the geometric mean of the Ct values and PCR amplification efficiency, and combines them into an index41. The comparative delta Ct method compares basic Ct values and the relative expression of ‘gene pairs’ within each sample, to confidently identify reliable RGs42.

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 method60. 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 (LC30 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 method60. 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).


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.

Electronic supplementary material


This work was supported in part by the Special Fund for Grain Scientific Research in the Public Interest of China (201413007-02; 201513002-05), the National key Research and Development Program of China (2016YFD0401004), the National Natural Science Foundation of China (31000828) and the Priority Academic Program Development of Jiangsu Higher Education Institutions.

Author Contributions

Author Contributions

P.A.T. and M.L.Y. designed the study. J.Y.D., H.J.W. performed the experiment. J.Y.D. and H.J.W. analysed the data. P.A.T. provided reagents. P.A.T., J.Y.D., H.J.W. drafted the manuscript. M.L.Y. and X.R.J. revised the manuscript.


Competing Interests

The authors declare that they have no competing interests.


Electronic supplementary material

Supplementary information accompanies this paper at doi:10.1038/s41598-017-07430-2

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Contributor Information

Pei-An Tang, moc.361@naiepgnat.

Ming-Long Yuan, nc.ude.uzl@lmnauy.


1. Cox PD, Atkinson JE, Bannon KL, Hatfield J. Studies of refuge-seeking behaviour and locomotory activity in Cryptolestes ferrugineus (Stephens) (Coleoptera: Cucujidae) using time-lapse video recordings. J. Stored Prod. Res. 1993;29:283–289. doi: 10.1016/0022-474X(93)90043-4. [Cross Ref]
2. Jagadeesan R, Nayak MK, Dawson K, Byrne V, Collins PJ. Dietary media for mass rearing of rusty grain beetle, Cryptolestes ferrugineus (Stephens) and flat grain beetle, Cryptolestes pusillus (Schonherr) (Coleoptera: Cucujidae) J. Stored Prod. Res. 2013;55:68–72. doi: 10.1016/j.jspr.2013.08.004. [Cross Ref]
3. Wu Y, et al. Microsatellite markers for Cryptolestes ferrugineus (Coleoptera: Laemophloeidae) and other Cryptolestes species. Bull. Entomol. Res. 2016;106:154–160. doi: 10.1017/S0007485315000899. [PubMed] [Cross Ref]
4. Tay WT, Beckett SJ, De Barro PJ. Phosphine resistance in Australian Cryptolestes species (Coleoptera: Laemophloeidae): perspectives from mitochondrial DNA cytochrome oxidase I analysis. Pest Manag. Sci. 2016;72:1250–1259. doi: 10.1002/ps.3805. [PubMed] [Cross Ref]
5. Nayak MK, et al. Strong resistance to phosphine in the rusty grain beetle, Cryptolestes ferrugineus (Stephens) (Coleoptera: Laemophloeidae): its characterisation, a rapid assay for diagnosis and its distribution in Australia. Pest Manag. Sci. 2013;69:48–53. doi: 10.1002/ps.3360. [PubMed] [Cross Ref]
6. Manivannan S, Koshy GE, Patil SA. Response of phosphine-resistant mixed-age cultures of lesser grain borer, Rhyzopertha dominica (F.) to different phosphine-carbon dioxide mixtures. J. Stored Prod. Res. 2016;69:175–178. doi: 10.1016/j.jspr.2016.08.005. [Cross Ref]
7. Champ, B. R. & Dyte, C. E. Report on the FAO global survey of pesticide susceptibility of stored grain pests. 297 (FAO, Rome, 1976).
8. Bajracharya NS, Opit GP, Talley J, Gautam SG, Payton ME. Assessment of Fitness Effects Associated with Phosphine Resistance in Rhyzopertha dominica (F.) (Coleoptera: Bostrichidae) and Tribolium castaneum (Herbst) (Coleoptera: Tenebrionidae) Afr. Entomol. 2016;24:39–49. doi: 10.4001/003.024.0039. [Cross Ref]
9. Kaur R, et al. Phosphine resistance in India is characterised by a dihydrolipoamide dehydrogenase variant that is otherwise unobserved in eukaryotes. Heredity. 2015;115:188–194. doi: 10.1038/hdy.2015.24. [PMC free article] [PubMed] [Cross Ref]
10. Daglish GJ, Nayak MK, Pavic H, Smith LW. Prevalence and potential fitness cost of weak phosphine resistance in Tribolium castaneum (Herbst) in eastern Australia. J. Stored Prod. Res. 2015;61:54–58. doi: 10.1016/j.jspr.2014.11.005. [Cross Ref]
11. Chen Z, Schlipalius D, Opit G, Subramanyam B, Phillips TW. Diagnostic Molecular Markers for Phosphine Resistance in U.S. Populations of Tribolium castaneum and Rhyzopertha dominica. PLoS ONE. 2015;10:e0121343. doi: 10.1371/journal.pone.0121343. [PMC free article] [PubMed] [Cross Ref]
12. Pimentel MAG, Faroni LRDA, Tótola MR, Guedes RNC. Phosphine resistance, respiration rate and fitness consequences in stored-product insects. Pest Manag. Sci. 2007;63:876–881. doi: 10.1002/ps.1416. [PubMed] [Cross Ref]
13. Nguyen TT, Collins PJ, Duong TM, Schlipalius DI, Ebert PR. Genetic conservation of phosphine resistance in the rice weevil Sitophilus oryzae (L.) J. Hered. 2016;107:228–237. doi: 10.1093/jhered/esw001. [PMC free article] [PubMed] [Cross Ref]
14. Daglish GJ, Nayak MK, Pavic H. Phosphine resistance in Sitophilus oryzae (L.) from eastern Australia: Inheritance, fitness and prevalence. J. Stored Prod. Res. 2014;59:237–244. doi: 10.1016/j.jspr.2014.03.007. [Cross Ref]
15. Benhalima H, Chaudhry MQ, Mills KA, Price NR. Phosphine resistance in stored-product insects collected from various grain storage facilities in Morocco. J. Stored Prod. Res. 2004;40:241–249. doi: 10.1016/S0022-474X(03)00012-2. [Cross Ref]
16. Zettler JL, Keever DW. Phosphine resistance in cigarette beetle (Coleoptera: Anobidae) associated with tobacco storage in the south eastern United States. J. Econ. Entomol. 1994;87:546–550. doi: 10.1093/jee/87.3.546. [Cross Ref]
17. Pimentel MAG, Faroni LRDA, Guedes RNC, Sousa AH, Tótola MR. Phosphine resistance in Brazilian populations of Sitophilus zeamais Motschulsky (Coleoptera: Curculionidae) J. Stored Prod. Res. 2009;45:71–74. doi: 10.1016/j.jspr.2008.09.001. [Cross Ref]
18. Pimentel MAG, Faroni LR, Silva FH, Batista MD, Guedes RN. Spread of phosphine resistance among brazilian populations of three species of stored product insects. Neotrop Entomol. 2010;39:101–107. doi: 10.1590/S1519-566X2010000100014. [PubMed] [Cross Ref]
19. Nayak MK, Collins PJ. Influence of concentration, temperature and humidity on the toxicity of phosphine to the strongly phosphine-resistant psocid Liposcelis bostrychophila Badonnel (Psocoptera: Liposcelididae) Pest Manag. Sci. 2008;64:971–976. doi: 10.1002/ps.1586. [PubMed] [Cross Ref]
20. Nath NS, Bhattacharya I, Tuck AG, Schlipalius DI, Ebert PR. Mechanisms of phosphine toxicity. J. Toxicol. 2011;2011:494168. doi: 10.1155/2011/494168. [PMC free article] [PubMed] [Cross Ref]
21. Chaudhry MQ. A review of the mechanisms involved in the action of phosphine as an insecticide and phosphine resistance in stored-product insects. Pestic. Sci. 1997;49:213–228. doi: 10.1002/(SICI)1096-9063(199703)49:3<213::AID-PS516>3.0.CO;2-#. [Cross Ref]
22. Chefurka W, Kashi KP, Bond EJ. The effect of phosphine on electron transport in mitochondria. Pestic. Biochem. Physiol. 1976;6:65–84. doi: 10.1016/0048-3575(76)90010-9. [Cross Ref]
23. Jian F, Jayas DS, White NDG. Toxic action of phosphine on the adults of the copra mite Tyrophagus putrescentiae (Astigmata: Acaridae) Phytoprotection. 2000;81:23–28. doi: 10.7202/706197ar. [Cross Ref]
24. Pimentel MAG, Faroni LRA, Corrêa AS, Guedes RNC. Phosphine-induced walking response of the lesser grain borer (Rhyzopertha dominica) Pest Manag. Sci. 2012;68:1368–1373. doi: 10.1002/ps.3314. [PubMed] [Cross Ref]
25. Burger G, Lang BF, Reith M, Gray MW. Genes encoding the same three subunits of respiratory complex II are present in the mitochondrial DNA of two phylogenetically distant eukaryotes. Proc. Natl. Acad. Sci. USA. 1996;93:2328–2332. doi: 10.1073/pnas.93.6.2328. [PubMed] [Cross Ref]
26. Bustin S. Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems. J. Mol. Endocrinol. 2002;29:23–39. doi: 10.1677/jme.0.0290023. [PubMed] [Cross Ref]
27. Kozera B, Rapacz M. Reference genes in real-time PCR. J. Appl. Genet. 2013;54:391–406. doi: 10.1007/s13353-013-0173-x. [PMC free article] [PubMed] [Cross Ref]
28. Tan Y, Zhou XR, Pang BP. Reference gene selection and evaluation for expression analysis using qRT-PCR in Galeruca daurica (Joannis) Bull. Entomol. Res. 2017;107:359–368. doi: 10.1017/S0007485316000948. [PubMed] [Cross Ref]
29. Zhang S, et al. Identification and validation of reference genes for normalization of gene expression analysis using qRT-PCR in Helicoverpa armigera (Lepidoptera: Noctuidae) Gene. 2015;555:393–402. doi: 10.1016/j.gene.2014.11.038. [PubMed] [Cross Ref]
30. Koramutla MK, Aminedi R, Bhattacharya R. Comprehensive evaluation of candidate reference genes for qRT-PCR studies of gene expression in mustard aphid, Lipaphis erysimi (Kalt) Sci. Rep. 2016;6:25883. doi: 10.1038/srep25883. [PMC free article] [PubMed] [Cross Ref]
31. Yuan M, et al. Selection and evaluation of potential reference genes for gene expression analysis in the brown planthopper, Nilaparvata lugens (Hemiptera: Delphacidae) using reverse-transcription quantitative PCR. PLoS ONE. 2014;9:e86503. doi: 10.1371/journal.pone.0086503. [PMC free article] [PubMed] [Cross Ref]
32. Shen GM, Jiang HB, Wang XN, Wang JJ. Evaluation of endogenous references for gene expression profiling in different tissues of the oriental fruit fly Bactrocera dorsalis (Diptera: Tephritidae) BMC Mol. Biol. 2010;11:76–76. doi: 10.1186/1471-2199-11-76. [PMC free article] [PubMed] [Cross Ref]
33. Lourenço AP, Mackert A, Cristino AdS, Simões ZLP. Validation of reference genes for gene expression studies in the honey bee, Apis mellifera, by quantitative real-time RT-PCR. Apidologie. 2008;39:372–385. doi: 10.1051/apido:2008015. [Cross Ref]
34. Ponton F, Chapuis M-P, Pernice M, Sword GA, Simpson SJ. Evaluation of potential reference genes for reverse transcription-qPCR studies of physiological responses in Drosophila melanogaster. J. Insect Physiol. 2011;57:840–850. doi: 10.1016/j.jinsphys.2011.03.014. [PubMed] [Cross Ref]
35. Jiang HB, Liu YH, Tang PA, Zhou AW, Wang JJ. Validation of endogenous reference genes for insecticide-induced and developmental expression profiling of Liposcelis bostsrychophila (Psocoptera: Liposcelididae) Mol. Biol. Rep. 2009;37:1019. doi: 10.1007/s11033-009-9803-0. [PubMed] [Cross Ref]
36. An XK, Hou ML, Liu YD. Reference gene selection and evaluation for gene expression studies using qRT-PCR in the white-backed planthopper, Sogatella furcifera (Hemiptera: Delphacidae) J. Econ. Entomol. 2016;109:879–886. doi: 10.1093/jee/tov333. [PubMed] [Cross Ref]
37. Ma KS, et al. Identification and validation of reference genes for the normalization of gene expression data in qRT-PCR analysis in Aphis gossypii (Hemiptera: Aphididae) J. Insect Sci. 2016;16:17. doi: 10.1093/jisesa/iew003. [PubMed] [Cross Ref]
38. Zheng YT, Li HB, Lu MX, Du YZ. Evaluation and validation of reference genes for qRT-PCR normalization in Frankliniella occidentalis (Thysanoptera:Thripidae) PLoS ONE. 2014;9:e111369. doi: 10.1371/journal.pone.0111369. [PMC free article] [PubMed] [Cross Ref]
39. Vandesompele, J. et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3, research0034.1–research0034.11 (2002). [PMC free article] [PubMed]
40. Andersen CL, Jensen JL, Ørntoft TF. Normalization of Real-Time Quantitative Reverse Transcription-PCR Data: A Model-Based Variance Estimation Approach to Identify Genes Suited for Normalization, Applied to Bladder and Colon Cancer Data Sets. Cancer Res. 2004;64:5245–5250. doi: 10.1158/0008-5472.CAN-04-0496. [PubMed] [Cross Ref]
41. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper – Excel-based tool using pair-wise correlations. Biotechnol. Lett. 2004;26:509–515. doi: 10.1023/B:BILE.0000019559.84305.47. [PubMed] [Cross Ref]
42. Silver N, Best S, Jiang J, Thein SL. Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR. BMC Mol. Biol. 2006;7:33–33. doi: 10.1186/1471-2199-7-33. [PMC free article] [PubMed] [Cross Ref]
43. Oppert B, et al. Genes related to mitochondrial functions are differentially expressed in phosphine-resistant and -susceptible Tribolium castaneum. BMC Genomics. 2015;16:968. doi: 10.1186/s12864-015-2121-0. [PMC free article] [PubMed] [Cross Ref]
44. Wang K, Peng X, Zuo Y, Li Y, Chen M. Molecular cloning, expression pattern and polymorphisms of NADPH-cytochrome P450 reductase in the bird cherry-oat aphid Rhopalosiphum padi (L.) PLoS ONE. 2016;11:e0154633. doi: 10.1371/journal.pone.0154633. [PMC free article] [PubMed] [Cross Ref]
45. Pang R, et al. Functional analysis of CYP6ER1, a P450 gene associated with imidacloprid resistance in Nilaparvata lugens. Sci. Rep. 2016;6:34992. doi: 10.1038/srep34992. [PMC free article] [PubMed] [Cross Ref]
46. Elzaki MEA, et al. Constitutive overexpression of cytochrome P450 associated with imidacloprid resistance in Laodelphax striatellus (Fallén) Pest Manag. Sci. 2016;72:1051–1058. doi: 10.1002/ps.4155. [PubMed] [Cross Ref]
47. Zhang QL, et al. Selection of reliable reference genes for normalization of quantitative RT-PCR from different developmental stages and tissues in amphioxus. Sci. Rep. 2016;6:37549. doi: 10.1038/srep37549. [PMC free article] [PubMed] [Cross Ref]
48. Sun M, Lu MX, Tang XT, Du YZ. Exploring valid reference genes for quantitative real-time PCR analysis in Sesamia inferens (Lepidoptera: Noctuidae) PLoS ONE. 2015;10:e0115979. doi: 10.1371/journal.pone.0115979. [PMC free article] [PubMed] [Cross Ref]
49. Robichaux S, Lacour N, Bagby GJ, Amedee AM. Validation of RPS13 as a reference gene for absolute quantification of SIV RNA in tissue of rhesus macaques. J. Virol. Methods. 2016;236:245–251. doi: 10.1016/j.jviromet.2016.08.004. [PMC free article] [PubMed] [Cross Ref]
50. Rajarapu SP, Mamidala P, Mittapalli O. Validation of reference genes for gene expression studies in the emerald ash borer (Agrilus planipennis) Insect Sci. 2012;19:41–46. doi: 10.1111/j.1744-7917.2011.01447.x. [Cross Ref]
51. Liu J-B, et al. Identification of appropriate reference genes for gene expression studies by quantitative real-time PCR in Tribolium castaneum after exposure to phosphine. J. Agric. Biotechol. 2014;22:257–264.
52. Hammond JW, Cai D, Verhey KJ. Tubulin modifications and their cellular functions. Curr. Opin. Cell Biol. 2008;20:71–76. doi: 10.1016/ [PMC free article] [PubMed] [Cross Ref]
53. Scharlaken B, et al. Reference gene selection for insect expression studies using quantitative real-time PCR: The head of the honeybee, Apis mellifera, after a bacterial challenge. J. Insect Sci. 2008;8:33. doi: 10.1673/031.008.3301. [Cross Ref]
54. Ballard JWO, Melvin RG. Linking the mitochondrial genotype to the organismal phenotype. Mol. Ecol. 2010;19:1523–1539. doi: 10.1111/j.1365-294X.2010.04594.x. [PubMed] [Cross Ref]
55. Nakakita H. The inhibitory site of phosphine. J. Pestic. Sci. 1976;1:235–238. doi: 10.1584/jpestics.1.235. [Cross Ref]
56. Corrêa AS, et al. Are mitochondrial lineages, mitochondrial lysis and respiration rate associated with phosphine susceptibility in the maize weevil Sitophilus zeamais? Ann. Appl. Biol. 2014;165:137–146. doi: 10.1111/aab.12127. [Cross Ref]
57. Zuryn S, Kuang J, Ebert P. Mitochondrial modulation of phosphine toxicity and resistance in Caenorhabditis elegans. Toxicol. Sci. 2008;102:179–186. doi: 10.1093/toxsci/kfm278. [PubMed] [Cross Ref]
58. Jagadeesan R, Collins PJ, Daglish GJ, Ebert PR, Schlipalius DI. Phosphine resistance in the rust red flour beetle, Tribolium castaneum (Coleoptera: Tenebrionidae): Inheritance, gene interactions and fitness costs. PLoS ONE. 2012;7:e31582. doi: 10.1371/journal.pone.0031582. [PMC free article] [PubMed] [Cross Ref]
59. Pfaffl MW. A new mathematical model for relative quantification in real-time RT–PCR. Nucleic Acids Res. 2001;29:e45. doi: 10.1093/nar/29.9.e45. [PMC free article] [PubMed] [Cross Ref]
60. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25:402–408. doi: 10.1006/meth.2001.1262. [PubMed] [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group