Identification and Validation of Reference Genes and Their Impact on Normalized Gene Expression Studies across Cultivated and Wild Cicer Species

Quantitative Real-Time PCR (qPCR) is a preferred and reliable method for accurate quantification of gene expression to understand precise gene functions. A total of 25 candidate reference genes including traditional and new generation reference genes were selected and evaluated in a diverse set of chickpea samples. The samples used in this study included nine chickpea genotypes (Cicer spp.) comprising of cultivated and wild species, six abiotic stress treatments (drought, salinity, high vapor pressure deficit, abscisic acid, cold and heat shock), and five diverse tissues (leaf, root, flower, seedlings and seed). The geNorm, NormFinder and RefFinder algorithms used to identify stably expressed genes in four sample sets revealed stable expression of UCP and G6PD genes across genotypes, while TIP41 and CAC were highly stable under abiotic stress conditions. While PP2A and ABCT genes were ranked as best for different tissues, ABCT, UCP and CAC were most stable across all samples. This study demonstrated the usefulness of new generation reference genes for more accurate qPCR based gene expression quantification in cultivated as well as wild chickpea species. Validation of the best reference genes was carried out by studying their impact on normalization of aquaporin genes PIP1;4 and TIP3;1, in three contrasting chickpea genotypes under high vapor pressure deficit (VPD) treatment. The chickpea TIP3;1 gene got significantly up regulated under high VPD conditions with higher relative expression in the drought susceptible genotype, confirming the suitability of the selected reference genes for expression analysis. This is the first comprehensive study on the stability of the new generation reference genes for qPCR studies in chickpea across species, different tissues and abiotic stresses.


Introduction
Gene expression studies have become extremely important to obtain insights into gene function and to understand the molecular mechanisms. Among various techniques used for gene expression studies, quantitative Real-Time PCR (qPCR) has become the most important and reliable method owing to its accuracy and high-throughput analysis [1]. While appropriate applications of qPCR requires robust reference genes for accurate normalization, non-availability can result in experimental deviations or errors that inevitably occur during sample preparation leading to unreliable quantification of gene transcripts. Ideally, the endogenous genes selected as reference genes for qPCR for normalization of gene expression data should be expressed stably in all plant tissues under various experimental conditions [2]. Various housekeeping genes such as, glyceraldehydes-3-phosphate dehydrogenase (GAPDH), actin (ACT) tubulin (TUB) cyclophilin (CYP), elongation factor (EF1) and ribosomal RNA (18S and 28S rRNA) have been widely used as reference genes in normalization of qPCR data. However, several studies have revealed that the expression levels of many commonly used housekeeping genes vary across tissues, treatments and species [1,[3][4][5]. Hence, emphasis has been on the identification of new generation reference genes that are stable under different experimental conditions [3,6]. Since several studies have shown that no single universal gene has consistent expression under all experimental conditions, the evaluation of reference gene(s) under specific experimental conditions is essential for reliability of qPCR analysis [7,8]. Moreover, several algorithms such as geNorm [9], NormFinder [10], BestKeeper [11] and comparative ΔCt method [12] have been developed to evaluate the most stable reference gene(s) from set of candidate genes.
Chickpea is an important food legume of the semi-arid tropical (SAT) regions of the world, known to be a nutraceutical (or health benefiting food) because of its high nutritional value [36]. Despite growing demand and high yield potential, chickpea yield is unstable and productivity is stagnant at unacceptably low levels due to constraints such as abiotic stresses (drought, heat, cold and high-salinity) and biotic stresses (Ascochyta blight, Fusarium wilt and pod borer). With chickpea genome sequencing allowing high-powered functional genomics studies to proceed [37], these can significantly accelerate molecular breeding efforts for the discovery and introgression of stress tolerance genes into cultivated germplasm. Nevertheless, for crop improvement in the post genomic era, there is a need to understand the function of genes in response to various stresses and during stages of growth and developmental, thereby necessitating gene expression profiling to identify the candidate genes.
Keeping this in view, we have evaluated a set of reference genes including traditional (commonly used) and new generation reference genes in a diverse set of biological samples of chickpea, including nine genotypes representing cultivated and wild species across primary, secondary and tertiary gene pools, plant tissues from various developmental stages, and six abiotic stress treatments (drought, salt, high vapor pressure deficit, abscisic acid, cold and heat shock). To identify the most stable reference gene(s) for normalization of qPCR data, this study evaluated 25 reference genes, including ATP-binding cassette transporter (ABCT), alcohol dehydrogenase class-3 (ADH3), calcium-dependent protein kinase 4 (CDPK4), clathrin adaptor complexes medium (CAC), cyclophilin (CYP), eukaryotic elongation factor 1-alpha (ELF1a), eukaryotic elongation factor 1-beta (ELF1b), galactose oxidase/kelch repeat superfamily protein (FBOX), glucose-6-phosphate dehydrogenase (G6PD), glyceraldehyde3-phosphate dehydrogenase (GAPDH), heat shock cognate protein 80 (HSP80), translation initiation factor IF-3 (IF3), eukaryotic initiation factor 4A-15 (IF4a), peroxin4 (PEX4), protein phosphatase2A subunit A3 (PP2A), pentatrico peptide repeat superfamily protein (PPR), s-adenosyl methionine decarboxylase (SAMDM), SAND-family protein (SAND), F-box protein SKIP16-like (SKIP16), TIP41-like protein (TIP41), un-characterized conserved protein UCP022280 (UCP), unknown protein (UNK), ubiquitin-protein ligase 7 (UPL7), vacuolar protein sorting-associated protein 53 homolog (VPS) and yellow leaf specific protein 8 (YLS8). The expression stability of these genes was evaluated across gene pools, different tissues and stress treatments using geNorm, NormFinder and RefFinder algorithms. Further, we have selected the best and least stable reference genes from the all samples set and validated by normalizing the expression levels of two aquaporin genes PIP1;4 (Plasma membrane intrinsic protein) and TIP3;1 (Tonoplast intrinsic protein) of chickpea. These two aquaporin genes express differentially in susceptible and tolerant genotypes of chickpea under different abiotic stress conditions (Unpublished data). The expression levels of the selected aquaporin genes under vapor pressure deficit (VPD) stress were tested in three chickpea genotypes with different tolerance levels.

Plant materials and stress treatments
Seeds of the chickpea cultivars and wild species were obtained from the mini core collection of the International Crop Research Institute for the Semi-Arid Tropics (ICRISAT), Hyderabad, India. Chickpea plants were grown in 8-inch pots containing 4.5 kg of black clay soil (Vertisol, 20% water holding capacity) under glasshouse conditions with 28/20°C day/night temperature. Nine chickpea genotype representations from cultivated genotypes and wild species across primary, secondary and tertiary gene pools (Table 1) were used for leaf tissue sample collection under normal growth conditions after 28 d from sowing. Besides, different tissue samples from seedlings of var. JG11 including leaf, flower, seed and roots were collected at different growth

Selection of candidate reference genes and primer design
The candidate genes were selected based on a review of studies in model plant Arabidopsis and other grain legumes such as soybean, chickpea and peanut followed by an in silico analysis using BLAST tools of the NCBI database. A total of 28 candidate reference genes that showed highly stable expression in previous studies including CAC, FBOX, GAPC2, PEX4, PP2A, PPR, PTB1, SAMDM, SAND, TIP41, UCP, UNK, UPL7 and YLS8 in Arabidopsis [3]; ABCT, CDPK4, IF3 and SKIP16 in soybean [40]; ELF1a, GAPDH, HSP80 and IF4a in chickpea [35]; ACT, ADH3, CYP, ELF1b and G6PD in peanut [1] and the VPS gene from chickpea (unpublished) were selected for evaluation. To retrieve the orthologous mRNA and their corresponding DNA sequences in chickpea, mRNA sequences of the Arabidopsis genes, EST sequences of the chickpea, soybean and peanut were used to query databases of the National Center for Biotechnology Information (NCBI) using BLASTX. The retrieved mRNA sequences of the chickpea sequences were used to design PCR primers using Primer3 software [41] with GC content between 45 and 50%, primer length of 20-22 nucleotides, and an expected product size of 100-150 base pairs. Most of the primer pairs were designed from two adjacent exons, which were separated by an intron (Table 2).

Genomic DNA isolation and PCR
The genomic DNA from leaves of chickpea variety JG11 was isolated using NucleoSpin Plant II DNA isolation kit (Macherey-Nagel, Duren, Germany). PCR reactions for all the candidate genes were performed using genomic DNA as template and gene-specific primers designed for qPCR. The PCR was performed in a total volume of 25 μl using 100 ng of genomic DNA, 200 nM of each primer, 1.5 mM MgCl 2 , 200 μM dNTP and 1 U Taq polymerase (Invitrogen, life technologies, NY, USA). PCR samples were amplified in an Eppendorf Thermal Cycler with an initial denaturation at 95°C for 5 min followed by 35 cycles of 95°C for 1 min, 62-64°C for 1 min and 72°C for 3 min, followed by a final extension at 72°C for 10 min, and PCR products tested by using 1% agarose gel electrophoresis.

RNA extraction and cDNA synthesis
Total RNA was extracted from 100 mg tissue by using NucleoSpin RNA plant kit (Macherey-Nagel, Duren, Germany) including in-column DNAse1 treatment. The concentration and purity of all RNA samples was tested using NanoVue plus spectrophotometer (GE health care, USA) at 260/280 nm absorbance, selecting the ones that ranged from 1.8 to 2.0. Integrity of the RNA was tested by using 1.4% agarose gel electrophoresis with SYBR safe DNA gel stain (Invitrogen-life technologies, USA). The total RNA sample (2.5 μg) was reverse transcribed to cDNA using SuperScript 1 III first strand cDNA synthesis kit (Invitrogen, Life Technologies, NY, USA) and oligodT primer in 25 μL reaction, according to manufacturer's protocol. The cDNA preparations were diluted 12-times with nuclease-free water (Qiagen, Valencia, CA, USA) to use as template in qPCR analysis. To confirm the total absence of genomic DNA, cDNA was used as a template for PCR amplification using ADH3 and G6PD primer pairs spanning an intron and PCR was performed as mentioned above.

Quantitative real-time PCR analysis
All qPCRs were carried out in Realplex (Eppendorf, Germany) Real-Time PCR system using SYBR Green in 96 well optical reaction plates (Axygen, Union City, CA, USA) sealed with ultra-clear sealing film (Platemax). The PCR reaction was performed in a total volume of 10 μL containing 5 μl of 2X SensiFAST TM SYBR No-ROX (Bioline, UK) mix, 400 nM of each primer, 1.0 μL of diluted cDNA and nuclease-free water to make up the final volume. The reaction conditions were set as 2 min at 95°C (polymerase activation); 40   using a pooled cDNA sample of the chickpea variety JG11. The 12-times diluted pooled cDNA sample was taken and 2-fold serial dilutions were carried out. Five serially diluted cDNA samples were used as templates to construct the standard curves of each primer pair, where PCR composition and conditions were same as above. Standard curves were constructed using the RealPlex (Eppendorf, Germany) real time PCR instrument software by using linear regression based on the quantitative cycle (Cq) values for all dilution points in a series. The correlation coefficients (R 2 ) and slope values were obtained from the standard curve, and corresponding PCR amplification efficiencies (E) were calculated using the slope of the calibration curve according to the following equation: E = (10 −1/slope -1).

Data analysis
All experimental samples were divided into different sets based on their nature. The tissue set comprised of leaf, flower, root, seedling and seed of chickpea variety JG11 grown under normal growth conditions (five samples). The abiotic stress treatment set comprised of eight leaf samples of chickpea variety JG11 grown under drought, salt, high VPD, ABA, cold, heat shock, drought control and control. The genotypes sample set comprised of leaf sample of four cultivars and five wild species (total nine genotypes listed in Table 1). All the 20 samples (JG11 was a common control in all sample sets) were analyzed together as 'all samples set'. The geNorm and NormFinder algorithms of genEX Professional software (MultiD Analyses AB, Sweden) were used to identify stably expressed gene(s) in a set of sample analyzed. The raw Cq values of each gene were corrected according to their respective PCR efficiencies before converting them into relative quantities, and the mean values of the biological replicates were taken as the input data for the geNorm and NormFinder analysis. The pair-wise variation analysis was carried out to identify the optimal number of genes required for normalization in each sample set by using geNorm of qBase plus software (ver: 2.4; Biogazelle, Belgium) [42]. RefFinder, a webbased (http://www.leonxie.com/referencegene.php) tool that integrates the currently available four major computational programs {geNorm [9], NormFinder [10], BestKeeper [11] and comparative ΔCt method [12]} and calculates the geometric mean for the comprehensive ranking was used to calculate the comprehensive ranks.

Reference gene validation by aquaporin gene expression studies under drought
Two aquaporin genes PIP1;4 (Plasma membrane intrinsic protein) and TIP3;1 (Tonoplast intrinsic protein) of chickpea were selected as target genes for quantification of gene expression levels under high vapor pressure deficit (VPD) treatment. Three chickpea varieties with contrasting drought tolerance viz. a susceptible variety ICC1882 and two tolerant varieties ICC4958 and JG11 were selected for reference gene validations. Experimental conditions and sample collection followed were same as mentioned above. Gene expression levels of PIP1;4 and TIP3;1 were normalized using the two most stable reference genes (ABCT and UCP) and two least stable genes (CYP and SKIP16) individually and in combination. Relative expressions of these two aquaporin genes in drought stressed leaf samples were estimated by comparing with expression levels of leaf sample collected from GH grown control plants of same variety using the REST [43] software.

Primer specificity and PCR efficiencies
The candidate reference genes selected in this study represent different functional classes and gene families, including traditional and new generation reference genes based on reports in Arabidopsis, soybean, peanut and chickpea. Primers designed for 28 candidate reference genes were tested with cDNA of chickpea variety JG11. Since the primer pairs of ACT, GAPC2 and PTB1 genes yielded more than one band with cDNA as template, these genes were eliminated from further analysis. All the remaining 25 primer pairs except for PP2A (The PP2A forward primer was designed from the exon-exon junction region) yielded single PCR product of expected size with the genomic DNA (S1A Fig). The selected 25 candidate reference genes produced single band on agarose gel (S1B Fig) and a single peak in melting curve using the cDNA template (S2 Fig). The amplification products of cDNA of 25 primer pairs when sequenced and searched in GenBank (NCBI) using the BLASTN algorithm, showed matches with the retrieved mRNA sequences of the chickpea, thereby demonstrating the gene specificity of the primer pairs used for qPCR. The amplification efficiencies (E) of the candidate reference genes ranged from 0.90 to 1.09 (Table 2).

Expression profiling of the reference genes
A qPCR assay based on SYBR Green was used for transcript profiling of the 25 candidate reference genes across all 20 samples of chickpea including those from six abiotic stress conditions (drought, salt, high VPD, ABA, cold and heat shock), five different tissues (leaf, root, flower, seedlings and seed) and nine genotypes representing cultivated and wild species across the chickpea gene pools.
To minimize the variability associated with qPCR analysis, all RNA samples were taken in equal concentration and quality for conversion to cDNA. The expression levels of the all candidate reference genes were individually determined in each of the samples as the quantification cycle value (Cq). A relatively wide range of Cq values across all 20 samples of chickpea was observed for the 25 reference genes, suggesting various levels of transcript abundance of the genes analyzed (Fig 1).

Expression stability of the candidate reference genes
Since the stability value of the best reference gene or the best combination of genes may vary from one experimental set to another, the stability ranks of candidate genes in four sample sets was determined separately to identify the most suitable reference genes using geNorm and NormFinder algorithms. The geNorm program defines a stability measure (M) as the average pairwise variation between a gene and reference genes in a set of samples where genes with the lowest M values have the most stable expression. The lowest M value was observed for the pair CAC/ABCT (M = 0.46) corresponding to the most stable expression in all sample set, whereas the M values for SKIP16 and CYP was considerably higher than the rest of the candidate genes.
In tissue sample set, VPS and ABCT were most stable (M values 0.28) and PPR, CDPK were least stable (M values 1.20 and 1.12, respectively). Under abiotic stress conditions, while TIP41 and G6PD were most stably expressed than all other candidate reference genes tested, ELF1a and PPR were unstable. The genes YLS8 and PEX4 were most stable in the tested genotypes across various gene pools while CYP and SKIP16 were found to be least stable (Fig 2A-2D). Similarly, the stability ranks of candidate reference genes determined with NormFinder (Table 3) showed UCP as the most stable gene followed by ABCT and G6PD in the all sample set, whereas CYP, SKIP16 and ELF1a genes were least stable among the all candidate reference genes tested. The TIP41 expressed most stably under abiotic stress conditions followed by CAC and PEX4, whereas ELF1a and PPR remained least stable. Across genotypes and species, UCP, G6PD and VPS genes were identified as the most stable, whereas the CYP gene was identified as least stable followed by SKIP16 and ELF1a genes. In the tissue sample set, PP2A, CAC, ABCT were identified as highly stable genes, whereas PPR, CDPK, CYP were among the least stable ones, according to NormFinder analysis.

RefFinder analysis
The rankings of candidate reference genes in all four-sample sets as re-determined using the RefFinder were highly consistent with the ones obtained by using geNorm and NormFinder. The RefFinder analysis revealed that while UCP, G6PD, CAC and YLS8 reliably expressed in chickpea genotypes across species, ABCT, UCP, CAC and G6PD were most stable in all samples set. Similarly, while PP2A, ABCT, VPS and CAC genes were ranked as top four in differential tissue sample set, TIP41, CAC, G6PD and PEX4 were observed to be highly stable under abiotic stress conditions. The comprehensive ranking revealed that the CYP gene was most unstable gene in all the sample sets, except under abiotic stresses. Various other genes such as ELF1a and GAPDH in abiotic stress, CDPK in differential tissues, SKIP16 in genotypes across species and all samples were also found unreliable in this analysis (Table 4).

Optimal number of reference genes for normalization
The optimal number of reference genes required for accurate normalization was estimated using pairwise variation analysis of geNorm algorithm. Since pairwise variation analysis showed that V2/3 value (0.159) was higher than 0.15, an ideal normalization will require to include CAC, ABCT and TIP41 as reference genes to normalize gene expression data in all samples set. Nevertheless, in the abiotic stressed sample set, genotype samples set and different tissue sample set, only two genes would be sufficient since the V2/3 values in these three sample sets was inferior to the 0.15 cut-off level. However, different experimental sets clearly required a different set of reference genes for normalization of gene expression, based on pair wise variation analysis viz. G6PD and TIP41 genes for abiotic stresses sample set, and YLS8 and PEX4 genes for genotype sample set. The genes, ABCT and VPSP were required for normalization of gene expression data in different tissue sample set, where the V2/3 value was 0.125 (Fig 2E-2H). Reference Gene Identification and Validation across Cicer Gene Pools

Validation of the selected reference genes
To validate the selected reference genes in qPCR, the relative expression level of two aquaporin genes PIP1;4 and TIP3;1 were investigated in three chickpea varieties, ICC1882, ICC4958 and JG11 under VPD stress. For normalization of the relative expression level of two aquaporin genes, two most stable and two least stable reference genes as evaluated by RefFinder in all sample set were used (Table 4). This analysis revealed that the expression level of PIP1;4 did not shown significant difference between drought susceptible and tolerant varieties, whereas TIP3;1 shown difference, when normalized with ABCT and UCP reference genes individually or in combination. The relative expression levels of the two aquaporin genes were did not followed the above pattern, when normalized using CYP and SKIP16 as the two most unstable reference genes, compared to the expression levels obtained after normalization with TIP41 and CAC genes. Moreover, both genes PIP1;4 TIP3;1 did not shown significant up regulation in susceptible variety ICC1882, when relative expression levels were normalized with CYP and SKIP16 genes individually or in combination (Fig 3).

Discussion
Microarray and next generation RNA sequencing technology-based studies has indicated that transcription levels vary between species, between ploidy levels and between tissue types [44]. The qPCR has become the first choice for accurate quantification of gene expression profiles of selected genes in diverse biological samples due to its sensitivity, accuracy and high throughput [45]. Nevertheless, stably expressing internal reference gene(s) are required for accurate gene expression. However, based on the previous reports, no single gene expresses stably across the experimental conditions [7,8]. Hence making it, critical to identify a set of stable reference genes that show stable expression in distinct biological samples: across different genotypes, tissue, growth stages and under varied experimental conditions. Genome sequence for both desi and Kabuli type chickpea is now available to researchers and functional genomics will be the major focus area in chickpea crop improvement programs [46,37]. Since, qPCR will be a very useful tool in the identification of candidate genes that controls major traits by studying their expression profiles, it is critical to develop a reference gene tool-kit with stably expressing internal genes for transcript normalization under specific experimental conditions. While, a previous report has evaluated 12 commonly used reference genes in a desi type chickpea var. ICC 4958 [35], the present study evaluated the expression stability of 25 candidate reference genes including new generation candidate reference genes. Moreover, assessing these genes in diverse tissue samples including, genotypes of chickpea representing both cultivated and wild species across gene pools, different tissues/developmental stages and abiotic stress conditions, ensures the robustness of our results owing to a wide range of variability covered in our samples. The expression profiling showed that all the 25 candidate reference genes varied significantly across the 20 samples tested. These results indicated that the most stable gene(s) among all the tested genes must be determined statistically. The stability value of the most suitable reference gene or the best combination of genes may vary from one experimental set to another, and for that reason the choice for reference genes for optimal normalization was made from a set of candidate genes for each experiment and for all of them together. To identify most stable reference gene under different experimental conditions, we grouped the diverse chickpea tissue samples into four sample sets and validated the 25 candidate reference genes using geNorm and NormFinder, two distinct statistical algorithms. Although, several studies have reported differences between the outputs of geNorm and NormFinder [8,47,48], we observed a high correlation between these rankings. In the present study, geNorm or NormFinder analysis in all four sample sets revealed a high degree of similarity. From the top 10 stable genes identified, 6 to 8 genes were almost same, with slight changes in their ranking orders. Interestingly, three most unstable genes remained same in all sample sets with both the algorithms (Table 3). To determine the comprehensive ranking of candidate reference genes in a sample set, we also used RefFinder that considered together the results of the four algorithms (geNorm, NormFinder, bestkeeper and deltCt method). Based on the stability rankings, different candidate reference genes have been proposed for their best suitability as reference genes depending on the experimental conditions. This being an important consideration would imply that the reference genes used for biotic or abiotic stress studies in chickpea might not be suitable for gene expression studies across genotypes and species.
Interestingly, new generation candidate reference genes selected for the current study performed better than the traditionally employed reference genes. In particular, CAC (Arabidopsis homologue: AT5G46630) was identified among the top four reference genes in all four-sample sets, indicating that CAC transcript levels were stable in chickpea under different experimental conditions tested. In this study, the CAC gene was selected as a candidate reference gene based Reference Gene Identification and Validation across Cicer Gene Pools on its highly stable expression in the Arabidopsis [3]. Besides, CAC has been reported to be the best reference gene for normalization in tomato [6], drought stressed roots and genotypes of coffee spp. [8,17], during fruit development in cotton [18], across total developmental stages and cultivars in mustard [14], and at various developmental stages in bamboo [49]. Based on the present and previous studies, it can be concluded that the CAC gene is stably expressed in different experimental conditions.
The uncharacterized conserved protein UCP022280 (UCP, Arabidopsis homologue: AT4G26410) was most stably expressed in different genotypes of chickpea and ranked second in all sample set. Previously, UCP gene was identified as stably expressed gene in Arabidopsis microarray analysis [3], and different experimental conditions [50]. Interestingly, UCP gene expression has reportedly been moderately stable in Arabidopsis under metal stress [4] and buckwheat [51]. In contrast, the expression of UCP gene was most unstable in vegetative tissues and maturing embryos of rapeseed [52]. An ATP-binding cassette transporter (ABCT) gene known to play a variety of cellular roles such as auxin transport was also selected for validation in the present study based on the microarray and qPCR validation in soybean, where ABCT has been reported as most suitable for gene expression normalization [40]. In the present study, ABCT gene was most stable when all samples were analyzed together and second most stable gene in differential tissue sample set, which is in agreement with the study in soybean where the ABCT gene was expressed stably under various abiotic stress conditions [53]. Besides soybean, this study is the first to report ABCT as the stable reference gene in crop species. Similarly, the TIP41-like family protein (TIP41; Tonoplast Intrinsic Proteins like protein; At4g34270) coding gene that was selected based on the Arabidopsis microarray data [3] emerged as most stable gene under abiotic stress conditions, whereas ranked fifth in different tissues and sixth in all sample set. These observations are in accordance with previous studies where TIP41 gene stably expressed under different abiotic conditions in desert shrub [54] and Chinese celery [55]. TIP41 gene also showed stable and similar expression regulation across Brassica species regardless of experimental conditions [14,16,52] and has been reported to be one of the most stable reference gene in many other plant species including soybean [56], tomato [6], buckwheat [51], peanut [57], common hop [58], cucumber [59], bamboo [49], Chinese pear [60] and olive tree [61].
Since our previous studies indicated stable expression of glucose-6-phosphate dehydrogenase (G6PD) gene under different experimental conditions in peanut [1], this was selected as one of the candidate in the present study. Similar to peanut, G6PD gene ranked in top four most stably expressed genes across different chickpea genotypes under abiotic stress conditions and all samples set. However, these results are in contrast to the findings in soybean where G6PD gene expression was least stable under developmental stages and different photoperiodic treatments [56], under cadmium stress [62], and also varied most in all tested experimental conditions [63]. The vacuolar protein sorting-associated protein 53 homolog (VPS) gene selected as candidate reference gene in the present study has not been previously explored as reference gene for qPCR analysis, it was identified as low copy gene in the legumes (unpublished data) and subsequently cloned. Interestingly, VPS gene expressed stably in different tissues and showed moderate stability across all samples and genotypes, ranking 5 th and 6 th respectively, indicating its potential for inclusion as reference gene in expression studies in other crop species.
The expression of CYP, ELF1a and IF4a genes were highly variable among the genes evaluated under different experimental conditions. The CYP gene was least stable in all samples set, across genotypes and different tissues, and showed high variability under abiotic stress. These observations were intriguing since this was in contrast to our previous studies in peanut, where CYP gene expressed most stably in vegetative stages and under abiotic stress conditions [1], and various tissues of soybean [56] and common bean [32]. Similarly, in this study ELF1a gene was found least stable under abiotic stress conditions and showed highly variable expression in other three sample sets, whereas IF4a genes was highly variable in all the four sample sets. These observations are in contrast with the previous study in chickpea [35] where ELF1a and IF4a genes were reported among the top three stably expressed genes in all three experimental conditions. Besides, the other two genes GAPDH and HSP80 that were previously reported to have stable expression [35], showed high variability under abiotic stress, but were moderately stable in other three sample sets. These differences justify the evaluation of this new set of reference genes indicating that, these new generation reference genes performed better than the conventional reference genes.
According to Vandesompele et al. [9] and as per the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines [45], two or more reference genes were required for accurate normalization of gene expression data. Several other reports also suggested that the application of more than one reference gene could lead to more reliable data in qPCR assays [18,64]. The geNorm proposed 0.15 as the threshold of pairwise variation values, below which the inclusion of an additional reference gene is not required. In this study, three out of four sample sets (genotypes, different tissue and abiotic stress sample sets) showed pairwise variation V2/3 values below 0.15, thereby indicating that the two reference genes were enough for normalization under these experimental conditions, whereas for all sample set three reference genes were required for normalization. These high pairwise variation values in all sample set may have resulted due to a highly diverse sample set, i.e., highly diverse samples might have significant differences in gene transcription including those of the candidate reference genes and were likely related to the developmental and metabolic properties of each distinct tissue.
To validate the selected candidate reference genes based on RefFinder ranking, expression of two chickpea aquaporin genes (PIP1;4 and TIP3;1) that belong to two different aquaporin sub-families (unpublished data) were assessed. Plant aquaporin proteins play major role in water transport through cell membranes and these genes are expressed differentially in different tissues that are also altered under different abiotic stresses in higher plants. The chickpea TIP3;1 gene was significantly up regulated under high VPD conditions in leaf tissues of three chickpea genotypes with relatively higher expression in susceptible genotype ICC1882, when normalized with ABCT and UCP reference genes. These results were in accordance with our previous study in sorghum [65] and in chickpea (unpublished data), where up-regulation of TIP genes in leaf tissues was reported under different abiotic stress conditions, but the normalization was obscured when the least stable reference gene(s) (CYP and SKIP16) were used. The validation results clearly indicated the demerits of using unstable reference gene(s) for normalization and confirmed the reference gene stability of the selected candidate genes.

Conclusions
The new generation candidate reference genes (CAC, UCP, ABCT, G6PD, VPS and TIP41) selected for the present study performed significantly better than the traditionally employed reference genes that were stably expressed across different experimental conditions. However, the CYP and ELF1a genes that were most unstable should be avoided for use as reference genes in gene expression studies in chickpea. Pairwise variation analysis indicated that two genes are enough for normalization of gene expression data in different sample sets except when all samples were analyzed together, where three genes are required. The validation of selected reference gene by normalizing aquaporin gene expression data further confirmed the stability of the ABCT and UCP as reference genes under high VPD conditions. These results provide information that can ensure more accurate qPCR based gene expression quantification towards chickpea functional genomics.
Supporting Information S1 Fig. Amplification of specific PCR products: with genomic DNA (A) and cDNA (B) as templates using gene-specific primers for each candidate reference gene tested in the study. 1 to 25 indicates the loading order of the candidate reference genes as mentioned in Table 2, M-DNA size marker. All primer pairs except CYP, FBOX, HSP80 and PPR amplified a larger size PCR product with DNA template as compared to cDNA template, indicating the position of primer pairs spanning at least one intron.