Reference genes validation in Phenacoccus solenopsis under various biotic and abiotic stress conditions

Real-time PCR (RT-qPCR) expression analysis is a powerful analytical technique, but for normalization of data requires the use of stable reference genes. However, suitable reference genes are still not known in the case of Phenacoccus solenopsis under variable experimental treatments. The present study focused on the identification of stable housekeeping genes as a reference for analysis under different abiotic and biotic factors in P. solenopsis. We analyzed the relative expression of six commonly used candidate reference genes in different developmental stages, host-feeding assay, temperature treatments and field distribution conditions. Expression stabilities were analyzed by geNorm, NormFinder, and RefFinder. Under developmental and field distribution conditions, β-Tubulin was found to be most stable reference genes followed by rpl32 and α-Tubulin. In the case host feeding treatment conditions, β-Tubulin and α-tubulin identified to be the most stable reference genes, while in temperature stress, a combination of α-Tubulin and rpl32 found to be suitable for normalizing the RT-qPCR data. Further, the above-identified genes were validated using RT-qPCR based gene expression analysis of four objective genes namely, Myoinhibitory peptides (MIPs), Zinc_metalloprotease (Zn_Mp), fatty acid synthase (fas) and alpha-glucosidase. Identified reference genes will facilitate gene expression studies in future under different stress treatments in P. solenopsis.

meaningful expression data 11,12 . The RT-qPCR analysis depends upon the appropriate selection of the reference gene(s) in the experimental set up which should be stably expressed across various tissue developmental stages and experimental conditions [12][13][14][15] . A good reference gene should have a ubiquitous expression, very low variance, and reasonable stability under different experimental conditions. Previous studies have shown that a single gene cannot stand stable under all biotic and abiotic stress conditions 16,17 . Therefore, it is imperative to use two or more reference genes for accurate normalization. For gene expression analysis through RT-qPCR in insects a few housekeeping genes such as Actin (ACT), GAPDH (Glyceraldehyde-3-phosphate dehydrogenase), 18S rRNA and EF-1a (Elongation factor-1a) have been commonly used as reference genes 16,[18][19][20][21][22][23][24] . However, several studies reported inconsistency in their expression that questioned their suitability as reference genes in RT-qPCR studies 14,25 . In fact, variable expression of housekeeping genes has also been reported depending on the organism, developmental stages, and various experimental conditions 26 . Therefore, none of the genes can be universally used as reference gene and hence it is necessary to validate the stability of reference gene under desired experimental conditions to ensure proper normalization during RT-qPCR analysis in the target species 14 . In most of the hemipteran species including P. solenopsis, many reference genes have been used without their prior validation. Studies on validation of the genes for RT-qPCR is now becoming prominent in other species but still lack for the mealybug species under different environmental conditions. As a result, considering the certainty of biotic and abiotic stress related internal control genes for accurate normalization of transcripts in RT-qPCR, the present work demonstrate the expression stability of six individual reference genes namely Actin (Act), Ribosomal protein (rpl32), β-tubulin, α-tubulin, Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and Succinate dehydrogenase (SDH) in P. solenopsis is subjected to individual biotic and abiotic stress treatments. Further, the identified genes were validated using RT-qPCR based gene expression analysis of three objective genes namely, MIPs, Zn_MPs and fatty acid Synthase and alpha-glucosidase.
Samples were taken into two biotic conditions (Sub group: developmental stages and field distribution) and two abiotic condition (Sub group: temperature variations and host feeding treatment). Further in each sub group, two variable conditions were taken as described in material and method section, and the detailed flow-chart of the experiment performed is shown in supplementary Fig. S1 (Total sample size, n = 48). RNA was isolated from all samples separately and treated with DNaseI to remove DNA contamination ( Supplementary Fig. S2), Further PCR reaction carrying DNaseI treated RNA as a template does not showed any amplification while desired band size was obtained in PCR (Used cDNA as template), further ruled out any possibility of DNA contamination ( Supplementary Fig. S2).
RT-qPCR primers specificity was analyzed using agarose gel electrophoresis and melt curve for the individual gene. Single desired size band and amplicon peak suggested the absence of non-specific PCR products and primer-dimers (Fig. 1). The PCR Amplification efficiencies (E) for all the primer pairs were calculated using a fivefold dilution series of pooled cDNA and measured in triplicate. PCR efficiencies (E) and correlation coefficients (R2) of the six candidate reference genes were evaluated by plotting the standard curves with their slopes obtained by serial dilution that fall under allowed range of 1.9 to 2.0 for E, and 0.97-1.00 for the R2. The description of six candidate reference genes, accession numbers, primer sequences, amplicon lengths, slope, R2 and amplification efficiency are enlisted in Supplementary Table S1.

Expression stability of the candidate reference genes
The expression profile of all six candidate reference genes was assessed by RT-qPCR, and the Ct values indicating differential transcript level were examined under four subgroups additionally all-samples subgroup (total five). The mean Ct values of these genes ranged from 18.63 to 27.76 (Table 1). The distribution of Ct values is listed in Supplementary Table S2.
In all-samples subgroup, the mean Ct values ranged from a minimum of 18.63 ± 0.58 to the maximum of 27.76 ± 5.68 with highest and lowest expression levels for ACT and SDH, respectively. In the field distribution subgroups, the mean Ct values showed a minimum of 16.27 ± 0.33, and the maximum of 25.12 ± 0.42 with highest and lowest expression levels for α-Tubulin and SDH, respectively. In developmental stages distribution mean Ct values for minimum and maximum expression were 17.60 ± 0.12 and 25.59 ± 0.36 for α-Tubulin, and SDH, respectively. Rpl32, α-Tubulin and SDH also showed a minimum and maximum average Ct values in abiotic stress including temperature treatments and host feeding assays sample set. The two most abundant genes were Actin (18.63 ± 0.58 cycles) and GAPDH (Ct = 19.11 ± 0.37 cycles) which showed moderate variation in Ct range of 0.91 and 0.96 cycles, respectively. Additionally, CV of the Ct values was also calculated to evaluate the expression levels of candidate reference genes under all five experimental subgroups, where higher values represent higher variability or minimum stability and vice-versa. The CV of six reference genes among all samples ranged between 1 to 2%. α-Tubulin was the least variable reference gene with a CV of 1.26% and GAPDH was the most variable with a CV of 1.81%. The stability ranking of all candidate reference genes on the basis of CV values was as follows: GAPDH < β-Tubulin < Actin < rpl32 < SDH < α-Tubulin (Table 2). During gene expression profiling across the sample sets in all the samples subgroups (field distribution; Fig. 2a, developmental stages; Fig. 2b, temperature variations; Fig. 2c, and different host feeding assay; Fig. 2d) separately, genes Actin, β-Tubulin, rpl32, SDH and GAPDH showed significant variation in Ct value (p-value < 0.05; Mann-Whitney U Non-parametric test) in biotic stress in comparison to the abiotic stress treatments (Supplementary Fig. S3). The results indicated that none of the candidate reference genes had constant expression levels across the different sample sets of P. solenopsis. Therefore, it was necessary to validate the expression stability of various candidate reference genes.

Expression stability ranking of candidate reference gene
Expression stability of the candidate reference genes was evaluated by geNorm (expression stability values M), NormFinder and RefFinder which calculated the stability values of each reference gene using the Ct values across all the experimental sets (Tables 3 and 4). The stability ranking analyses carried out by each of the three algorithms is described in the following sections.    Table 2. Details of C t and CV values of each of the selected candidate reference genes tested in Phenacoccus solenopsis.  Table 3 and Supplementary Table S3).   Table 4. Expression stability ranks of six candidate reference genes in abiotic stress treatments and all sample groups calculated using geNorm (GN), NormFinder (NF) and RefFinder methods.

Optimal number of internal candidate genes for normalization
To produce accurate and consistent results, two or more genes are required for the RT-qPCR experiment. To obtain the optimal number of reference genes for data normalization, we calculated the pairwise variation (Vn/n + 1) of serial log-transformed NF ratios using N relative to N + 1 reference genes (i.e. log2 (NFn/NFn + 1) as previously described 27 . Since the individual reference gene had considerably differential expression across the samples, NF could be sensitive to the stepwise inclusion of these reference genes resulting in an increase or decrease in Vn/n + 1 value. A cut-off value of V = 0.15 was recommended in deciding the addition of the next reference gene. A large pairwise variation with a significant effect means that the added gene is preferably included for calculation of a reliable normalization factor. In contrast, an extra reference gene is not necessary for normalization when the cut-off value is below 0.15. As shown in Fig. 4, the Pairwise variations of all biotic treatment sample set, V2/3 value were above the cut-off value 0.15, which indicates that two reference genes were not enough for accurate normalization. The addition of the third or fourth reference gene was required for significant normalization of the results. In our study, the V2/3 values of the developmental set of group were above the proposed 0.15 cut-off. But V3/4 values were below the cut-off and these results indicate that normalization with three (α-Tubulin, rpl32 and β-Tubulin) stable reference genes was required to produce valid expression data set. In all samples sets, each Pairwise variation values were below the proposed 0.15 cut-off. In the case of field distribution, the V3/4 was close to the proposed 0.15 cut-off. Moreover, the inclusion of additional reference genes did not lower the V value below the recommended 0.15 cut-off until the fourth or five genes was added. In temperature treatments, the V2/3 values were below the proposed 0.15 cut-off value (Fig. 4). According to geNorm, two reference genes (β-Tubulin and GAPDH) should be required for a suitable normalization in the high or low temperature-stressed samples. In the experimental set of host feeding assay, all the pair-wise variation V2/3 values were below the recommended 0.15 cut-off. These results indicate that normalization with additional reference genes did not lower the V value below the proposed 0.15 cut-off.
NormFinder analysis. The result of NormFinder analysis presented that the stability ranking of the six candidate reference genes was relatively consistent with the data array of geNorm findings. The results of gene expression stability analysis are shown in Supplementary Table S4. For the all samples set rpl32, β-Tubulin and SDH were identified as the best three optimal internal control genes with the stability value of 0.08, 0.16, and 0.21, while ACT with stability value (SV) 0.49 showed higher variation. In the case of field distribution, ACT and GAPDH were unstable with the stability value of (0.47), (0.55). But in the sample set of developmental stages, ACT and SDH showed stability value of (0.77), (0.78) and said to be very unstable. The rpl32, β-Tubulin and α-Tubulin found to be highly stable with the stability of 0.16, 0.17 and 0.70 in developmental conditions, respectively. The β-Tubulin and rpl32 were also very stable and showed the stability value of (0.38 and 0.42) at field distribution sample sets (Table 3). This means that α-Tubulin, β-Tubulin, and rpl32 were the best stably expressed gene under developmental subgroups, and β-Tublin and rpl32 performed better in the field conditions. In the sample set temperature variations, ACT and β-Tubulin were found to be very unstable with the stability values of (0.23 and 0.22, respectively) while α-Tubulin and rpl32 behave stable with the stability of (0.08 and 0.02). Under different host feeding assays, β-Tubulin and ACT found to be very unstable with stability values of (0.44, 0.63), while rpl32 and GAPDH were stable with stability value of (0.24, 0.56). The results shown with NormFinder was similar to the one obtained from geNorm to some extent (Tables 3 and 4). Overall, NormFinder predicted the reference genes α-Tubulin, β-Tubulin and rpl32 were the highly stable genes in development and host feeding assay sample, which was set similar to the results obtained by geNorm analysis. In the case of field distribution and temperature treatments, α-Tubulin and rpl32 were the highly stable genes and ACT, GAPDH were the unstable gene in the all sample set subgroups.
RefFinder analysis. The RefFinder analysis revealed that α-Tubulin and β-Tubulin were most stable genes in field distribution samples subgroup. However, β-Tubulin and rpl32 were ranked as top in different stages sample subset. The rpl32, and α-Tubulin genes were the most stable genes under different temperature conditions, whereas GAPDH, rpl32 and α-Tubulin were observed to be highly stable in different host feeding assays. The comprehensive ranking also revealed that the β-Tubulin and rpl32 were the most stable genes in the all biotic stress treatment samples, while rpl32 and α-Tubulin in abiotic stress treatments. Other candidate genes such as ACT in the field, different stages and host feeding assay sample subsets; and GAPDH in temperature treatments were found as the unstable genes in this analysis (Supplementary Table S5).

Expression analysis of MIP, Zn_Mp, alpha-glucosidase and fatty acid synthase gene from P. solenopsis under different biotic and abiotic stress conditions for reference genes validation.
In order to validate the selected reference gene, the differential expression of four objective genes namely Myoinhibitory peptide (MIP; for the developmental stage and host feeding), Zinc_metalloprotease (Zn_Mp; for developmental stage), fatty acid synthase (fas; for temperature stress) and alpha-glucosidase (for field distribution) were analyzed by RT-qPCR. The expression profiles of the genes MIP and Zn_Mp were normalized using the combination of most stable reference gene (β-tubulin, α-tubulin, and rpl32) and the least stable reference gene GAPDH and ACT recommended by geNorm, RefFinder, and NormFinder. Relative transcript accumulation of MIP was found to be unbiased when α-Tubulin, β-Tubulin, rpl32 or combination of α-Tubulin_β-tubulin_rpl32 were used confirming the high level of expression stability of these three genes. Using α-Tubulin, β-Tubulin, rpl32 or combination of it, the expression of MIP and Zn_Mp were increased under different developmental and host feeding treatment conditions (Fig. 5a,b). For field distribution and temperature treatment conditions, expression profiles of alpha glucosidase and Fas were normalized using the most stable reference genes α-Tubulin, β-Tubulin and rpl32 and combination of it (Fig. 5c,d). Fas, which plays a central role in lipid synthesis and in stress tolerance showed increase in their relative expression during any abiotic stress conditions 28 . The alpha-glucosidase has role in the metabolism of oligosaccharides and play major role in carbon nutrition and osmoregulation 29 . The relative expression pattern of MIP, Zn_Mp, alpha-glucosidase and Fas showed strong deviation and failed to achieve consistency when GAPDH and ACT for transcript normalization under development, host feeding treatments, field distribution and different temperature treatments.

Discussion
Gene expression analyses during different developmental stages and environmental conditions are extremely important to reveal the molecular mechanisms underlying various biological processes. RT-qPCR is the best technique even for low abundant mRNA transcripts and is becoming one of the most commonly utilized technology for the analysis of gene expression in recent years 7,13,26,30-33 . However, accurate normalization of gene expression remains a major concern for precise RT-qPCR analysis. The normalization of gene expression value reduces the variations occurred during sample preparations, various steps of RT-qPCR, cDNA preparations, PCR Figure 4. Determination of the optimal number of reference genes for geNorm analysis. The pairwise variation (Vn/Vn + 1) was analyzed for the normalization factors NFn and NFn + 1 by geNorm program to determine (V < 0.15) the optimal number of reference genes. efficiency, and various others. The majority of studies used one gene for their data normalization which produces an ambiguous result, and the gene was often selected on the basis of earlier reports rather than on their suitability for a particular experimental conditions. Utilization of right reference genes is a key strategy to minimize the variations in expression analysis. Earlier, it has been demonstrated the use of two or more reference genes during RT-qPCR provides more precise results 27 . In this context, there is an utmost need of suitable reference genes for expression analysis in economically devastating pest of cotton, P. solenopsis and other hemipteran pest of agricultural crops as well.
The PCR amplification efficiency is directly related to master mix performance, sample quality and assay performed. Our result suggested, the amplification efficiencies and correlation coefficients of the six candidate reference genes were ranges of 1.9 to 2.0, and 0.97-1.00, respectively, which is similar to the previous reports 34 . Numerous statistical algorithms such as geNorm 27 , NormFinder 11 , RefFinder (http://fulxie.0fees.us/?type=reference), BestKeeper 35 and delta Ct 36 have been developed to estimate the expression stability of individual reference genes, but the rankings of the best reference genes were different with the different programs. So, for accurate normalization in gene expression studies, a comparison of different algorithms allows better evaluation of gene expression patterns 37 . We found that the MV value of all analyzed six reference genes were less than 1.5 under different abiotic and biotic stress conditions, this indicating possible expression stability of selected reference genes 27 . All three of the algorithms (geNorm, NormFinder, and RefFinder) used in present study recommended α-tubulin, β-tubulin and rpl32 as the top three most stable reference genes in all the samples set across all the treatment conditions with some variability in their ranking order. α-Tubulin, β-Tubulin constituted the most stable reference gene pair with rpl32 according to geNorm in the biotic stress treatment conditions. We consistently found α-Tubulin, β-Tubulin as best reference genes this may be due to their structural requirement of the cell so their expression could be consistent in organism invariable to any treatment. Under developmental, β-Tubulin was found to be most stable reference genes followed by rpl32 and α-Tubulin. An earlier study on Frankliniella occidentalis suggested EF-1 and rpl32 are the most stable reference genes in different developmental stages 38 . Although rpl32 is common in both studies but we cannot compare these result simultaneously due to the difference in model insect and experimental result analysis. Further, α_tubulin and rpl32 are the best reference genes under different temperature conditions; these results are consistent with the earlier finding where α_tubulin and rpl32 selected as best reference genes under heat stress in P. solenopsis 39 . In field distribution, ACT and GAPDH were used as the stable gene to normalize the data in Aphis gossypii, but we found both genes were unstable in expression. This might be due to change in different climatic conditions of collection field 34 .
The applicability of the reported three most stable genes were assessed by analyzing their expression pattern of MIP, Zn_MP, alpha-glucosidase and fas under different biotic and abiotic stress conditions. The objective genes showed amplification with three of the internal control genes α-Tubulin, β-Tubulin, and rpl32 under different biotic and abiotic stress. However, they showed variations when gene expression was assessed by using the unstable genes, ACT, and GAPDH. The results showed that the expression stability of the six most commonly used reference genes in P. solenopsis was varied. Of these six reference genes; α-Tublin, β-Tublin, and rpl32 showed least changes in expression level in developmental stages, field distribution and host feeding assay treatment conditions. The ACT and GAPDH expression were found to be most variable during biotic stress and host feeding assay. None of the single reference genes showed constant expression in all the four different experimental conditions. Reportedly, expression of each gene was profoundly affected by the different experimental conditions on the chosen reference genes, though all of them were not changed by each treatment conditions. Collectively, the results re-confirmed the previous report that no single reference gene is showing the same expression across all the treatment conditions. The best reference gene(s) under different experimental conditions were indicated by rank order. However, our results showed that ACT and GAPDH are also not an ideal reference gene for all the studies of P. solenopsis, it is only suitable under certain experimental conditions and needed to cooperate with other reference genes to obtain the reliable results. For instance, in the temperature stress experiments, the combination of rpl32 and α-Tubulin could be suitable for normalizing the RT-qPCR data. In the developmental stage and field distribution condition, the use of combining the expression of three reference genes namely, α-tublin, β-Tubulin and rpl32 gave better results. Similarly, when the combination of α-Tubulin and β-Tubulin were used for normalizing the analytical RT-qPCR results of host feeding assays, an accurate result could be obtained. The results demonstrated that using unsuitable reference gene(s) for normalization might lead to deviated results. Therefore, it is a perquisite to select appropriate reference genes for accurate estimation of target gene expression.
In conclusion, we have investigated the expression of six candidate reference genes at different biotic and abiotic treatment conditions, an attempt to identify most suitable reference genes for normalizing gene expression in P. solenopsis. Our results suggested that we could use combination of reference genes for the different experimental set-up to avoid any variations in gene expression to achieve perfect normalization. To the best of our knowledge, this is one of the first reports on the evaluation of candidate reference genes across different experimental conditions in P. solenopsis and will provide guidance to other researchers to select reference genes against target genes by RT-qPCR in the related species. solenopsis were collected from reared population. Due to the lack of sufficient biological material, egg and first instar nymphs were pooled and considered as a single sample while adult insect was taken as separate sample (Total two samples). Further samples were collected from two different plant considers as biological duplicate, while the experiment was performed in triplicate with each sample.

Field distribution.
A field variation has a great impact on the growth parameters of the insect due to change in temperature and relative humidity at the place of maintenance of the insect populations and could be used as one of the factors for gene expression studies. Two geographical populations were collected from two physically separated fields (~800 Km) one from the cotton field of CSIR-NBRI, Lucknow (26°55′N, 80°59′E) and the other from Panjab University, Chandigarh (30.7601°N, 76.7663°E). Collected populations were multiplied in the net-house conditions and glass house conditions, separately. Hundreds of adult stage female populations were collected separately and stored in RNAse free micro-centrifuge tubes containing RNA later solution and stored in −80 °C until use.
Abiotic factors. Temperature treatments. Temperature has a great impact on the insect developmental cycles. To look out for variability in the growth parameters with a varied range of temperatures, adult mealybugs were exposed to two temperature regimes 15 °C and 42 °C for 1 h. In total, a hundred adults at each temperature points were collected for RNA extraction.
SCIeNtIfIC RepoRtS | 7: 13520 | DOI:10.1038/s41598-017-13925-9 Host plant-feeding assays. Feeding assay was designed to evaluate the growth of adult mealybugs on cotton and tobacco leaves. Leaves of both plants were collected from the glass house and washed thoroughly and dry on bloating sheet. Leaves were cut and placed in petridish (35 mm) pre-poured with 1% agarose. Adult female mealybugs were starved for 5 hours and released on a plate (5/dish). After 48 h, until it started secreting the honeydews, insects were gently brushed from leaves disc and stored in RNA later solution and stored in −80 °C until use.
Extraction of total RNA and cDNA synthesis. Total RNA was extracted using 500 μL of TRIzol ® Reagent (Sigma-Aldrich, St. Louis, USA) from pooled sample of Eggs and first instars nymph, and third instar and adult female mealybugs according to the manufacturer instructions. The extracted RNA was resuspended in 50 μL of Nuclease-free water (Life Technologies). RNA concentration of each sample was determined using NanoDrop-1000 spectrophotometer (NanoDrop Technologies). The OD260/OD280 nm absorption ratio (1.98-2.01) and OD260/OD230 (≥2.0), was used to determine the quality and purity of RNA preparations. RNA integrity was analyzed by means of electrophoresis in 1.5% agarose gels visualize after ethidium bromide staining using Bio-Rad Gel Doc systems (Bio-Rad, USA) ( Supplementary Fig. 1). Genomic DNA was eliminated using the Turbo-DNase kit (Ambion, USA) ( Supplementary Fig. 1). Further, the presence of DNA contamination was checked through PCR using DNase-treated RNA samples with −RT and +RT reaction. Minus-reverse transcriptase ("−RT") control is a mock reverse transcription containing all the RT-PCR reagents (no cDNA template), except the reverse transcriptase and +RT samples contained all RT-PCR reagents (cDNA template). RNA samples (1 µg) were immediately used to synthesize cDNA using the SuperScript III Reverse Transcriptase PCR and RT-qPCR. Amplification specificity of 6 primer pairs was confirmed by PCR. The PCR cocktail was having the volume of 10 μl contained 10× PCR buffer (ThermoFischer Scientific, USA), 10 mM dNTPs mix (ThermoFischer scientific, USA), 5 μM each of forward and reverse primers, 1 μl of cDNA, and 0.2 μl of Taq polymerase. The PCR program for amplification was as follows: an initial denaturation step of 3 min at 94 °C, and 30 cycles of 30 s at 94 °C, 45 sec at 60 °C and 45 sec at 72 °C followed by an extension step of 10 min at 72 °C. PCR products were separated on 1.5% agarose gels. The RT-qPCR was executed using SYBR Green detection chemistry on the 7500 fast Real-time PCR system v2.0.6 (Applied Biosystems). The reaction mixture contained 1 μl of 5 times diluted cDNA (equal to 50 ng of initial amount of RNA), 5 pico moles of each gene-specific forward and reverse primer, and 5 μl of 2 × FAST SYBR Green PCR Master Mix (Applied Biosystems, USA) in 10 μl of total reaction volume. The RT-qPCR cycling conditions were as follows: an initial denaturation step of 20 s at 50 °C, 10 min at 95 °C, and 40 cycles of 15 s at 95 °C, and 1 min at 60 °C followed by melt curve analysis using default parameters in order to check the PCR specificity by steady increase in temperature from 60 °C to 90 °C. The PCR efficiency of each primer set was determined through slope of the amplification curve in the exponential phase, obtained by five-fold dilution series of cDNA (1:5, 1:25, 1:125, 1:625, and 1:3125) using E = 10 (−1/slope)−1 (3.3 ≥ slope ≥ 3.1) by the software itself (Applied Biosystems). Triplicates of two biological replicates of each sample were used for real-time PCR analysis and three technical replicates were analyzed for each biological sample.
Data analysis for expression stability of candidate reference genes. The candidate reference genes expression levels were determined by threshold cycles (Ct). The stability and suitability of these genes were evaluated using three independent Excel based calculations and algorithms, i.e., geNorm 27 , NormFinder 11 , and RefFinder (http://fulxie.0fees.us/?type=reference). The Ct values were imported into the Excel sheet, and average expression stability value-M (AESV-M) for each candidate gene was calculated. Then, pair-wise variation (V) of each gene was calculated and compared with all other candidate reference genes. The geNorm stability value provides brief ideas about the expression stability of the genes on the basis of stability value, M. Through the NormFinder analysis, first ct values were converted into relative quantities and ranked the stability of the reference genes depending on the variation in gene expression. Genes with the lower stability values are considered to be most stable gene and are best to select as reference gene for that particular experimental conditions. RefFinder is a web-based (http://fulxie.0fees.us/?type=reference) online tool, which integrates the presently available four advanced computational programs geNorm 27 , NormFinder 11 , comparative ΔCt method 36 and BestKeeper 35,38 and calculates the geometric mean for the comprehensive ranking.
Validation of reference genes. On the basis of available literature 28,29,[40][41][42][43][44] and in-house transcriptome data (unpublished) four genes namely Myoinhibitory peptide (MIP), Zinc_metalloprotease (Zn_Mp), fatty acid synthase (Fas) and alpha-glucosidase were selected for validating three most stable reference gene α-Tubulin, β-Tubulin, and rpl32; and the two most varying reference genes namely GAPDH and ACT identified from this study. Sequence of the four objective genes was obtained from a P. solenopsis transcriptome in-house data (data not published yet) and real time primers were synthesized (supplementary Table 6). RT-qPCR was performed as described above. The relative mRNA expression data of objective genes normalized separately with internal control genes used was analyzed using 2 −ΔΔCt method 45 .