Comprehensive evaluation of candidate reference genes for real-time quantitative PCR (RT-qPCR) data normalization in nutri-cereal finger millet [Eleusine Coracana (L.)]

Finger millet (Eleusine coracana L.) is an annual herbaceous self-pollinating C4 cereal crop of the arid and semi-arid regions of the world. Finger millet is a food security crop proven to have resilience to changing climate and scores very high in nutrition. In the current study, we have assessed sixteen candidate reference genes for their appropriateness for the normalization studies in finger millet subjected to experimental regimes and treatments. Ten candidate reference genes (GAPDH, β-TUB, CYP, EIF4α, TIP41, UBC, G6PD, S24, MACP and MDH) were cloned and six (ACT, ELF1α, PP2A, PT, S21 and TFIID) were mined from the NCBI database as well as from the literature. Expression stability ranking of the finger millet reference genes was validated using four different statistical tools i.e., geNorm, NormFinder, BestKeeper, ΔCt and RefFinder. From the study, we endorse MACP, CYP, EIF4α to be most stable candidate reference genes in all ‘tissues’, whereas PT, TFIID, MACP ranked high across genotypes, β-TUB, CYP, ELF1α were found to be best under abiotic stresses and ‘all samples set’. The study recommends using minimum of two reference genes for RT-qPCR data normalizations in finger millet. All in all, CYP, β-TUB, and EF1α, in combination were found to be best for robust normalizations under most experimental conditions. The best and the least stable genes were validated for confirmation by assessing their appropriateness for normalization studies using EcNAC1 gene. The report provides the first comprehensive list of suitable stable candidate reference genes for nutritional rich cereal finger millet that will be advantageous to gene expression studies in this crop.


Introduction
Finger millet, Eleusine coracana (L.) Gaertn, is a nutri-cereal grown for food and feed in Africa and South Asia regions of the world. This tiny grain displays high variability in the mineral composition and has superior nutritional qualities compared to other cereal crops including rice and wheat used as a health food, and in bakery [1][2][3]. Besides  nutrition to millions of people in semi-arids, finger millet is quite a resilient crop, and hence has attracted lot of attention of researchers for studying its genetics, genomics for its improvement [4][5][6][7][8].
Finger millet is vulnerable to both abiotic and biotic stresses, with blast disease being a primary constraint [9,10] and drought and salinity stresses affecting the crop production systems and economics [11,12]. To overcome these stresses, there is a need to deploy beneficiary regulatory and structural genes through functional genomics approaches. Towards this RT-qPCR technology offers promise for studying the function of desired genes with high sensitivity, precision, simplicity and robustness [13][14][15]. Nonetheless, there are certain limitations to this technology essentially due to lack of the appropriate reference gene (s), which further effects the threshold (Cq) values and eventually affect the precision of the expression [16,17]. Experiment-to-experiment difference depends on the reference genes expression, which is quashed through the process of 'normalization'. Across most species, the most commonly used reference genes (RG) have been housekeeping genes (HKG) with the fundamental supposition that their expression levels remain unchanged regardless of the condition or nature of the sample during the course of the experiments [18,19]. The trustworthiness of the RT-qPCR data trusts on stable expression of the candidate reference genes across the conditions irrespective of the samples [20][21][22].
Considering that until now there is no work done in this direction for finger millet, the present study was undertaken for assessment of sixteen reference genes, including, Actin (ACT), Beta Tubulin (β-TUB), Cyclophilin (CYP), Eukaryotic Initiation factor 4A (EIF4α), Elongation factor 1-alpha (EF1α), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), Glucose-6-phosphate 1-dehydrogenase (G6PD), MalonylCoA-Acyl Carrier protein (MACP), Malate dehydrogenase (MDH), 40S ribosomal protein (S24), Serine/threonine-Protein Phosphatase (PP2A), Phosphate transporter protein (PT), Ribosomal protein (S21), Transcription initiation factor (TFIID), Tonoplast intrinsic proteins-like protein (TIP41) and Ubiquitin Protein Isoform C (UBC) as reference genes for RT-qPCR in finger millet. The expression stability of the sixteen reference genes across the regime of diverse experiments was evaluated using geNorm [19], NormFinder [47], BestKeeper [48] and ΔCt [49] statistical tools. To our knowledge it is the first attempt on a systematic evaluation of the reference genes in Eleusine coracana (L.). The conclusion of the study definitely will advantage those experiments which involve gene expression studies in finger millet species and also in other closely related millets.

Plant material and abiotic stress treatments
Finger millet (variety GPU 28) has been used for the different abiotic stress treatments and tissue/organ collection. GPU 28 seeds were sown in pots comprising red soil mixture (3:2:1 clay: sand:manure) and grown in a greenhouse with day/night average temperatures of 27/22˚C and relative humidity of 70-80%. Five major abiotic stresses (salt, cold, heat, drought and ABA) and different tissues (seedling, leaf, root, panicle, and mature seed) were harvested [20,21]. Finger millet cultivars contrasting for drought stress response (Tolerant-IE 4073, IE  4797 and GPU 28; sensitive-IE 5106 and IE 2572) [50] were grown under greenhouse conditions and imposed progressive drought stress after 28 days and leaf tissues were collected when normalized transpiration ratio (NTR) reached at 0.1. All the samples were collected in triplicates and straightaway snap-frozen in liquid nitrogen and stored at -80˚C till RNA isolation.

Sequence mining, cloning and RT-qPCR primers designing
Sixteen candidate reference genes, including ACT, ELF1α, PP2A, PT, S21 and TFIID were retrieved from the available sequence information of finger millet deposited in the NCBI database. Remaining ten-candidate reference genes were cloned from the sequence information of different plant species, including pearl millet (β-TUB, S24, CYP and GAPDH), chickpea (EIF4α, TIP41 and UBC), groundnut (G6PD) and sorghum (MACP and MDH). Two micrograms of finger millet total RNA was used for cDNA synthesis (Invitrogen, USA) and PCR amplification was carried out with respective gene specific primers according to the manufacturer's instructions (Invitrogen). The amplified PCR products were cloned into the pCR4.0-TOPO vector (Invitrogen) and sequenced. RT-qPCR primers were designed using primer 3.0 software (http://bioinfo.ut.ee/primer3-0.4.0/) [51] with default settings with the following considerations: (a) product size: 90-170 bp; (b) primer length: 18-24 bp and (c) GC of 45-55%. The primer details are listed in Table 1. Primer specificity was evaluated by 2.0% agarose gel electrophoresis and as well as with the melt curve analysis.

Calculation of PCR efficiency
Ten-fold serial dilution of cDNA was used as template for calculating the amplification efficiency (E) of the primer pairs in RT-qPCR with minimum five dilution points. The amplification efficiency (E) and correlation coefficients (R 2 ) for each primer set were estimated according to the equation: E = 10 −1/slope .

RNA isolation and RT-qPCR
Total RNA of the finger millet samples was isolated from 100 mg of tissue by RNeasy Plant Mini kit (Qiagen, Germany) according to manufacturer's guidelines. Quantity and quality of RNA were determined by NanoVue plus spectrophotometer (GE health care, USA) and BioAnalyzer (Agilent). Total RNA samples with an absorbance ratio OD 260/280 ranged from 1.9-2.2 were used directly for RT-qPCR analysis. Integrity of RNA was confirmed by running the samples on 1.4% denatured agarose gel electrophoresis. Further, total RNA of all the finger millet experimental samples was diluted to 100 ηg/μl and it used for RT-qPCR assays. All the RT-qPCR assays were accomplished using SYBR green based quantification assay in a Realplex real-time PCR machine (Eppendrof). A reaction mixture was constituted of 1 μl-RNA (100 ηg), 5-μl one step SYBR RT-PCR buffer 4 (Takara, Japan), 0.4 μl of the prime script one step Enzyme Mix 2 (Takara, Japan) and 400nm of each primer and total volume made to 10 μl with RNase-free H 2 O. The one step RT-qPCR including reverse transcription cycling were as follows: 42˚C for 5 min and 95˚C for 10 s, followed by 40 cycles of denaturing at 95˚C for 15 s, annealing at 62˚C for 15 s with fluorescent signal recording. The dissociation (melt) curve analysis was included after 40 cycles of amplification cycles are completed by heating from 58 C to 95˚C with fluorescence measured within 20 min. All the RT-qPCR assays were repeated at least three times.

Samples size and grouping
Experimental samples used in the current study were classified into four sample sets based on their sample nature. The 'tissue set' included of five different tissues of the plant development, i.e., seedling, root, leaf, panicle and mature seed of finger millet variety GPU 28 is grown under greenhouse conditions. The 'abiotic stress set' comprised of five samples of finger millet grown under different abiotic stresses (Heat, cold, salt, drought and ABA). The 'genotypes set' comprised a leaf sample of drought stress and controls of five finger millet genotypes, including three drought tolerant (IE 4073, IE 4797 and GPU 28) and two drought sensitive (IE 5106 and IE 2572). Total 60 samples considering three biological replicates of each sample, from above three sample sets were considered jointly as 'all sample set'. Table 1. Comprehensive details of the finger millet candidate reference genes, primer sequences, amplicon size, melting temperature (Mt), amplification efficiency (E); regression coefficient (R 2 ), coefficient of variation (CV). Ah-Arachis hypogaea; Ca-Cicer arietinum; Ec-Eleusine coracana; Pg-Pennisetum glaucum, Sb-Sorghum bicolor.

Statistical programs for normalization
Statistical tools named geNorm, NormFinder, BestKeeper and ΔCt were adopted for identifying the best stable candidate reference genes in finger millet. The Cq values of each gene were converted into relative quantities after adjusting them according to their respective PCR efficiencies. The mean values of the relative quantities of the replicates were acquired as the input data for the geNorm and NormFinder tools and data was analysed using genEX Professional software (MultiD Analyses AB). The geNorm tool calculates the expression stability (M) and ranks the genes stability in an order as the lower the M value indicates the higher the stability [19]. The geNorm tool of qBase plus software (ver: 2.4; Biogazelle, Belgium) also relates the pairwise variation (V) of the most stable genes with the rest of the candidate reference genes for efficient normalization in each sample set. A threshold value of 0.15 and less than that indicates no additional reference gene required for normalization in a particular sample set. The NormFinder statistical tool calculates intra-and inter-group variations in gene expression stability and provides ranks accordingly [47]. Genes with the lowermost rank values were considered to be most stably expressed reference gene(s). BestKeeper is an excel-based tool that calculates a Pearson's correlation coefficient for each reference gene, values of p closer to 1.0 indicating greater stability [48]. In the ΔCt tool, the rank order is determined based on pairwise comparisons of gene-sets and lowest standard deviation indicates highest expression stability of the reference gene [49]. RefFinder, a web-based tool (http://150.216.56.64/ referencegene.php) combines all four major statistical tools (geNorm, NormFinder, Best-Keeper and comparative ΔCt method) for calculation of the comprehensive ranks.

Reference gene validation
Abiotic stress inducible EcNAC1 gene [24] was selected for RT-qPCR data normalization in different genotypes under progressive drought stress. Four finger millet genotypes contrasting with drought stress tolerance viz. susceptible (IE 5106 and IE 2572) and tolerant (IE 4073 and IE 4797) were selected for quantification of the EcNAC1 gene. Treatments and sample collection were done as mentioned in the plant material and abiotic stress treatments section. Expression of EcNAC1 gene was normalized with two best (PT and TFIID) and two least stable finger millet reference genes (UBC and MDH) selected from the "genotypes set". The relative expression of EcNAC1 gene in progressive drought stressed leaf samples was assessed by comparing with respective control samples of same genotype and as well as with the selected combinations of the reference genes using the REST software [52].

Selection and cloning of candidate genes
Six candidate reference genes, including ACT, ELF1α, PP2A, PT, S21 and TFIID extracted from the finger millet genome sequence available in the NCBI database were used for the primer designing and further in RT-qPCR study. Remaining ten genes (β-TUB, S24, CYP, GAPDH, EIF4α, TIP41, UBC, G6PD, MACP and MDH) were cloned from finger millet cDNA by respective gene specific primers from the various plant sources (Table 1). Amplicons were further verified using agarose gel electrophoresis (Fig 1a) and confirmed by sequencing before being used for RT-qPCR primer designing.

Primer specificity and PCR efficiency analysis
The amplification specificity of the sixteen finger millet candidate reference genes was studied using regular PCR. The PCR amplification results revealed that all sixteen genes showed distinct and individual amplification of predictable product sizes when resolved on agarose gel (Fig 1a) and melt curve analysis by RT-qPCR (Fig 1b). Predictable product size, distinct PCR amplified fragment, single and sharp melt curve peak and sequencing data representing that the primers had high specificity and were appropriate for RT-qPCR assays. Linear regression coefficient (R 2 ) values ranged between 0.922 (CYP)-0.999 (G6PD and β-TUB) and the PCR amplification efficiency in different samples varied from 0.87 (β-TUB)-1.08 (CYP) ( Table 1). Linear regression coefficient (R 2 ) and PCR amplification efficiency (E) values were within the acceptable range demonstrating that the primers of the sixteen finger millet reference genes are very specific and applicable for the further analysis.

Expression analysis of the finger millet reference genes
The expression analysis of sixteen finger millet reference genes was studied in 21 different experimental samples collected from tissues, abiotic stresses and different genotypes of finger millet. The variable Cq values of sixteen candidate reference genes throughout the experimental samples suggested that their expression levels are highly diverse by the treatments and conditions (Fig 2). The expression of S21 (CV = 13.96) followed by GAPDH (CV = 9.95) were highly affected from sample to sample. In a relative evaluation, a slight range of variable Cq values was detected in case of TIP41 (CV = 2.47) and G6PD (CV = 2.88) empirically suggesting their stable expression under different experimental samples. Regardless of the condition, S21 was found to be abundant with the lowermost mean Cq value 15.57, and while TIP41 had highest mean of Cq 30.33 among the tested reference genes throughout the experiments (Fig 2). These implied that the expression of the finger millet reference genes are inconsistent and without any specific pattern across all the experimental conditions. Therefore, there is a need to select a best stable reference gene (s) to normalize the gene expression in a set of samples grouped based on their experimental nature in finger millet.

Expression stability analysis of the finger millet reference genes
Expression stability of the finger millet reference genes was studied using four different statistical algorithms having distinct principles for stability rankings in order to select the best stable reference genes. In the current study, we detected almost similar tendencies for each condition with subtle variations, which may be attributed by variances in the tools. Because each tool has their respective advantages we adopted four statistical packages in choice of the most and the Selecting reference genes for finger millet least stable candidate reference genes. Further precision in the interpretation of the stability ranking made by each statistical tool has been presented individually. geNorm. In geNorm, stability of the sixteen candidate reference genes notably changed in ranking order from sample to sample (Table 2 and S1 Table). In 'all sample set' geNorm ranking revealed MACP and PT (0.85) to be most stable followed by CYP (0.96), where TIP41 (1.47), S21 (1.42) and G6PD (1.37) were the least stable reference genes (S1 Table). Under 'abiotic stress set', CYP, S21 (0.25) showed higher stability and PT (1.14) and ACT (1.10) were the least stable ones. In the 'tissues set', CYP and MACP (0.33) were the most stable and S21 (1.46) and G6PD (1.39) were the least stable genes, whereas in 'genotypes set' MACP and S21 (0.45) were the highest in terms of their stability and MDH (1.25) and UBC (1.17) as the least stable amongst all. The pairwise variation V2/3 value was greater than 0.15 in 'all-samples set' and less than 0.15 in the case of 'abiotic stress set' (V2/3 = 0.13) (Fig 3). These results implied that use of more than two best stable finger millet reference genes together would be required for normalization studies in genotypes (V4/5 = 0.15), tissues (V4/5 = 0.13) and all sample sets (V6/7 = 0.15) (Fig 3).
BestKeeper. BestKeeper determines rankings based on the standard deviation (SD) values, which is inversely proportional to the expression stability of the genes. TIP41 (SD = 0.6) followed by β-TUB and G6PD (SD = 0.65) were most stably expressed and ranked high across all experimental samples, in contrast to, MDH expression that revealed significant variation (SD = 1.72) ( Table 2 and S3 Table). Intriguingly, gene-expression stability of β-TUB (SD = 0.07) was least affected in 'abiotic stress set', whereas S24 and PT (SD = 0.99) were the least stable genes. In 'tissue set', G6PD and β-TUB (SD = 0.47) displayed most stable expression, while TFIID (SD = 1.83) and MDH (SD = 1.77) were found to be the least stable ones. In the 'genotype' set, TIP41 (SD = 0.31) and G6PD (SD = 0.42) were placed as the best reference genes followed by β-TUB (SD = 0.53), while MDH (SD = 1.74) and ACT (SD = 1.58) as least useful (S3 Table).

Validation of the reference genes
To analyze the transcript level of EcNAC1 gene under drought stress in different genotypes, we selected two best stable (PT and TFIID) and least stable (UBC and MDH) finger millet reference genes identified in the present study from the 'genotype set' and were used for validation of the normalized results. PT and TFIID individually or combined showed higher stability in RT-qPCR data compared to least stable UBC and MDH both independently and in a combination (Fig 5). Normalization with PT and TFIID genes produced more consistent and comparable results in drought susceptible (IE 5106 and IE 2572) and tolerant (IE 4073, IE 4797) genotypes which clearly showed differential expression of transcripts of EcNAC1 gene (Fig 5). EcNAC1 gene expression was higher in drought tolerant genotypes (IE 4073 and IE 4797) when compared with susceptible ones (IE 5106 and IE 2572) and followed a similar trend when normalized with best stable reference genes. In contrast, normalizations using the lesser stable reference genes showed significant variation and inconsistency in the results with discrepancies in the drought tolerant and susceptible genotypes (Fig 5).

Discussion
The recent draft genome sequence of Eleusine coracana (L.) provides an excellent genomic resource for the research community that should facilitate a holistic understanding of the genetic basis of its innate nutritional potentials and drought stress tolerance [8]. A total of 85,243 genes from this genome sequencing data, predicted as stress related transcription factors/genes, calcium transporters accumulation genes and C 4 photosynthetic pathway genes, make it worthwhile to validate the function of the genes for understanding their role in various biological processes. So far most of the RT-qPCR analyses in finger millet have used reference genes either from heterologous plant species or from the native finger millet without any systematic normalization validation studies [23][24][25][26][27][28][29][30][31][32][33][34]. This becomes a serious hindrance in precise analysis of the RT-qPCR data and requires a comprehensive evaluation of the reference genes in finger millet under various experimental conditions.
In the this study, sixteen candidate reference genes (Table 1) were used to analyze their expression-stability under various experimental conditions referred to as 'genotypes', 'tissues', 'abiotic stress' and 'all samples' sets. The gene-expression data clearly indicate the effect of treatments on the stability of reference genes. The expression stability of all finger millet genes was not affected by each treatment, few genes got affected and others have not irrespective of the conditions (Fig 2). Together, the results reported earlier imply that no single reference gene is stable across all the conditions and treatments [20,21]. The best stable reference gene (s) for each condition were indicated by its order of rank using different statistical algorithms ( Table 2 and S1-S5 Tables). Although expression stability ranking order of the finger millet reference genes was not same in each method, a realistic consent undoubtedly implied the top order finger millet reference genes for each condition (Fig 4). For example, CYP and EF1α were positioned as best stable reference genes by NormFinder and ΔCt method for studying gene expression in all samples set, while geNorm ranked MACP as the best and BestKeeper placed TIP41 as most stable gene. The subtle variation in ranking order of the top listed reference genes could be endorsed to difference in tools of the software and sensitivities towards the co-expressed reference genes [43]. It was also evident that the top-ranked stable reference gene for one condition may not appropriate to the other condition. Hence, it was vital to find the stable reference genes that show adequate for RT-qPCR assay, if not the uppermost stability in expression within and across the conditions can be considered.
RefFinder study suggested CYP, β-TUB, and EF1α as the top stable finger millet reference genes for their stable expression in most of the conditions. These interpretations are in correlation with earlier findings where CYP and EF1α were described as most steadily expressed reference genes under diverse experimental conditions of Petunia and Vicia faba [31,53,54]. In alternative study, CYP gene was revealed to express stably in different developmental stages and EF1α was revealed to be least variant in expression in the samples collected at different time intervals in soybean [55]. Additionally, CYP gene was recommended as best stable gene in combination with other three reference genes for normalization studies in soybean [56]. The EF1α gene was also reported as best stable reference gene under various experimental conditions in different plant species including chinese cabbage [57], pearl millet [20], and potato [58]. The third most stable reference gene β-TUB was also recommended as best stable in many plant species. The genes β-TUB and EF1α were suggested in combinations with ACT gene for gene expression analysis under biotic and abiotic stress treated samples of Vigna mungo [59], and Musa [60]. The β-TUB gene was also reported as the most stable under leaf senescence conditions in sunflower [61]. There has been an extensive discussion about the optimal number of most stable reference genes required for RT-qPCR data analysis. For preventing errors and increasing the accuracy in normalization process, investigators have proved use of more than one reference genes instead of single gene [20,21,41,42]. In our present study, geNorm pair-wise analysis implies use of more than two stably reference genes in achieving Selecting reference genes for finger millet accuracy during RT-qPCR data normalization for all the experimental samples except abiotic stress, where two reference genes would be beneficial (Fig 3). Therefore, we propose the use of two and more reference genes in combinations for normalizing of gene expression assays under different experimental samples in finger millet.
In the present study EcNAC1 a drought and salinity stress responsive gene [24], was selected as experimental gene for the validation of normalization results (Fig 5). EcNAC1 gene was expressed high in tolerant genotypes (IE 4073 and IE 4797) in comparison with susceptible genotypes (IE 5106 and IE 2572) and followed distinctive pattern in expression when normalized with the stable reference gene (s) (PT and TFIID). In contrary, conflict and inconsistency was observed in expression levels when uses least stable reference gene (s) (UBC and MDH) (Fig 5). In summary, we conclude that the reference genes validated in the present study were appropriate for data normalization in RT-qPCR studies under various experimental conditions in finger millet.

Conclusion
Owing to the agronomic significance of Eleusine coracana as future food security crop, gene expression studies would endure to represent a significant part of basic and functional genomics research in finger millet. Therefore, establishing standardized reference genes for RT-qPCR studies in E. coracana would assist the peer researchers working in finger millet functional genomics. The present study reveals that the traditional and heterologous plant source candidate reference genes may not be appropriate for their direct use in RT-qPCR data normalization studies without systematic investigational validation across the experimental conditions. The stability rank order of the finger millet reference genes from all the experimental conditions implied that no single reference gene could be used perfectly for all the experimental conditions. In summary, we recommend the use of CYP, β-TUB and EF1α, preferably in combination for robust normalization of RT-qPCR data under most of the experimental conditions. The present study is helpful for undertaking the future RT-qPCR based expression studies in the finger millet.
Supporting information S1 Table.