Identification and validation of reference genes for quantitative RT-PCR normalization in wheat

Background Usually the reference genes used in gene expression analysis have been chosen for their known or suspected housekeeping roles, however the variation observed in most of them hinders their effective use. The assessed lack of validated reference genes emphasizes the importance of a systematic study for their identification. For selecting candidate reference genes we have developed a simple in silico method based on the data publicly available in the wheat databases Unigene and TIGR. Results The expression stability of 32 genes was assessed by qRT-PCR using a set of cDNAs from 24 different plant samples, which included different tissues, developmental stages and temperature stresses. The selected sequences included 12 well-known HKGs representing different functional classes and 20 genes novel with reference to the normalization issue. The expression stability of the 32 candidate genes was tested by the computer programs geNorm and NormFinder using five different data-sets. Some discrepancies were detected in the ranking of the candidate reference genes, but there was substantial agreement between the groups of genes with the most and least stable expression. Three new identified reference genes appear more effective than the well-known and frequently used HKGs to normalize gene expression in wheat. Finally, the expression study of a gene encoding a PDI-like protein showed that its correct evaluation relies on the adoption of suitable normalization genes and can be negatively affected by the use of traditional HKGs with unstable expression, such as actin and α-tubulin. Conclusion The present research represents the first wide screening aimed to the identification of reference genes and of the corresponding primer pairs specifically designed for gene expression studies in wheat, in particular for qRT-PCR analyses. Several of the new identified reference genes outperformed the traditional HKGs in terms of expression stability under all the tested conditions. The new reference genes will enable more accurate normalization and quantification of gene expression in wheat and will be helpful for designing primer pairs targeting orthologous genes in other plant species.


Background
Transcriptome and gene expression analyses are contributing to substantially improve our understanding of the signalling and metabolic pathways underlying developmental and cellular processes. Quantitative RT-PCR (qRT-PCR) is currently one of the most powerful and sensitive techniques for analyzing gene expression and, among its several applications, it is often used for validating output data produced by micro-and macro-arrays of wholegenomes and as a primary source for detecting specific gene expression patterns [1][2][3]. Reliable quantification by qRT-PCR analysis of gene expression levels requires the standardization and fine-tuning of several parameters, such as: amount of initial sample, RNA recovery and integrity, enzymatic efficiencies of cDNA synthesis and PCR amplification, overall transcriptional activity of the tissues or cells analyzed [4,5]. Among several proposed methods [5,6], internal control genes (reference genes) are most commonly used to normalize qRT-PCR and to reduce possible errors generated in the quantification of gene expression, which is obtained by comparing the expression levels in the analysed samples of the gene of interest and of stable constitutive control genes. Obviously, the success of this procedure relies on the choice of appropriate control genes, which ideally would be those showing stable expression under various experimental conditions and in different tissues types. Housekeeping genes (HKGs), whose protein products are involved in basic cellular processes and are supposed to have stable and uniform expression across different tissues and developmental stages, are commonly exploited as internal controls for normalization in gene expression analyses. The reference genes most widely used for qRT-PCR in different species are those encoding 18S rRNA, actin, tubulin, polyubiquitin and GAPDH, which have commonly been employed in traditional methods of gene expression detection, such as northern blot, RNase protection test and conventional semi-quantitative RT-PCR [6]. However, several reports have shown that frequently the most widely used HKGs are not reliable controls, since their expression level varies in different tissues [7][8][9][10][11][12][13], and it would be essential a preliminary evaluation for identifying the most stable HKGs in each species. It is rapidly increasing the number of articles on the identification and validation of novel and more stable reference genes in mammalians, as well as the programmes of statistical software for evaluating the stability of selected HKGs [13][14][15]. On the contrary, there are few specific studies on expression analysis of HKGs in plants, most of them evaluating the expression in different species and experimental conditions of the same well-known reference genes [16][17][18][19][20]. Recently a novel set of reference genes has been identified in Arabidopsis (using the large public collection of data from Affymetrix GeneChip experiments) [12], and in barley (on the basis of ESTs analysis and relative expression calculation for all the cDNA libraries available through the TIGR Barley Gene index database) [21]. These studies indicate that many new reference genes outperforming the traditional ones in terms of expression stability can be found by different approaches. However, the increased awareness of the importance of systematic validation has not permeated fully throughout the community of plant molecular biologists, and although the potentially highly misleading effects of using inappropriate references for normalization are widely known, they are also widely disregarded. Gutierrez et al. [22] have shown that genes validated as references have been used in only 3.2% of 188 real-time RT-PCR analyses published during a 6-month period from July through December 2007 in the three leading primary research journals in plant biology, according to the ISI Web of Knowledge (The Plant Cell, Plant Physiology, and The Plant Journal). In another recent paper Gutierrez et al. [23] have also shown that genes commonly used as references may be expressed unstably during the development of Arabidopsis plants. Using validated and non-validated reference genes, the authors illustrated to what extent the use of arbitrarily chosen genes could affect the quantification of target gene expression levels. Up to 100-fold variations were found in the expression of a target gene, which could be attributed only to variations in the expression of the reference genes, with consequently huge potential scope for misinterpretation of the results. Therefore, there is an urgent need to regard the systematic validation of reference genes as an essential component of real-time RT-PCR analysis to improve the reliability of published results and retain the accuracy of this powerful technique.
Wheat is the most widely grown crop in the world and has paramount importance for human nutrition. Cultivated bread wheat (Triticum aestivum L.) is an allohexaploid species (AABBDD) with three very large homoeologous genomes, each comprising seven pairs of chromosomes. The whole hexaploid wheat genome is about 40 times larger than that of rice, this huge difference is largely due to the inclusion of repetitive sequences [24]. The close synteny between the rice and wheat genomes shown by genetic maps based on cDNA markers [25,26] indicates that their physiological and developmental differences are mostly due to regulation and expression patterns of similar structural genes. Since the wheat genome is too large to be entirely sequenced in the next future, structural and functional analysis of the wheat transcriptome is particularly important. The screening of the whole wheat transcriptome is made possible by the recent availability of an Affymetrix GeneChip Wheat Genome Array containing 61,127 probe sets (55,052 transcripts), likely covering half of the wheat expressed genes [27][28][29], and of several application-specific cDNA arrays [30][31][32]. Since qRT-PCR is necessary for array validation and accurate expression studies, it is important to identify reliable internal control genes suitable for many experimental conditions. Using the terms "wheat", "gene expression", and "real time RT-PCR" linked by the Boolean operator "AND", we performed a PubMed search of articles published from January 1996 to March 2008 and retrieved 26 articles that used 16 different reference genes (Additional file 1). Remarkably, genes encoding 18S rRNA (8 times, 30%), actin (7 times, 27%) and alpha-tubulin (5 times, 19%) have been used in about 3/4 of the studies, whereas each of the remaining reference genes has been cited only 1-3 times (4-11% of the papers). Only four of the studies were based on multiple reference genes, whereas 22 researches used single reference genes, and 15 of them presumed their expression stability without any preliminary validation. Although 18S rRNA has been frequently used as a reference gene, it is far from ideal for several reasons: i) need of using total RNA and random primers for the RT reaction; ii) very high expression levels compared to the mRNA transcripts of the gene of interest; iii) imbalances in rRNA and mRNA fractions between different samples [33,34]; iv) expression affected by tissue and experimental conditions [4,33]; v) reduced degradation compared to mRNA [16]. Also the constitutive genes actin and alpha-tubulin, which are commonly used as reference genes in plants, do not represent the best choice because they are present in multiple copies showing high variation in their expression levels [35,36]. Our literature survey evidenced the lack of validated reference genes for expression studies in wheat and clearly emphasized the value of the present systematic study for identifying more reliable genes.
Microarray datasets can be a rich source of information for selecting qRT-PCR reference genes, as shown in Arabdopsis [12]; the same approach can not be followed in wheat for the limited availability of public microarray datasets [27][28][29][30][31][32]. As transcript abundance can be inferred by EST frequency, the computational analysis of ESTs from various cDNA libraries represents an alternative strategy. As of January 2008, about 1,000,000 wheat ESTs accumulated in the dbEST (NCBI) are freely available and represent an opportunity for extensive investigation to identify reference genes with stable expression across a wide range of developmental and environmental conditions. Wheat EST database contains sequences of cDNA libraries obtained from several tissues, including cell cultures, crown, spikes at different stage of development, flower organs (lemma, palea, anthers and pistils), seedlings, leaves, roots, stems, sheath, seeds at different stages of development, endosperm and embryo. Moreover, the same database includes EST sequences from cDNA libraries of wheat plants exposed to abiotic (salt, drought, high and low temperatures, aluminium) and biotic (Fusarium graminearum, powdery mildew, stripe rust, Hessian fly) stresses. The value of this resource is further enhanced by the availability of two public wheat-specific databases, Triticum aestivum UniGene and TIGR wheat Gene index, in which a set of non-redundant EST-index is produced. UniGene (NCBI) [37] is a system for partitioning GenBank sequences, including ESTs, into a non-redundant set of orientated clusters. Each UniGene cluster contains sequences representing a unique gene, which is linked to the tissue types where it is expressed and to characterized proteins of model organisms encoded by similar sequences. The over 947,303 wheat ESTs contained in GeneBank were reduced about 22-fold to 41,256 sequence clusters in T. aestivum UniGene April 2008 release (build 52). However, NCBI does not generate contigs and/or consensus sequences for UniGene clusters, whereas TIGR produces Gene Indices which are constructed by first clustering and then assembling ESTs sequenced at TIGR, ESTs from dbEST and annotated gene sequences from GenBank [38]. Each TIGR cluster contains a fasta formatted consensus sequence with a unique accession number (TC, Tentative Consensus), as well as additional information including details of the assembly, cDNA libraries from which the single assembled ESTs are derived and putative gene identification. The last release of T. aestivum Gene Index (release 11, July 13, 2008) contains 216,452 unique sequences, including 91,464 assembled TC sequences, 124,732 singleton ESTs and 256 singleton ETs (annotated sequences).
In the present research several novel candidate reference genes suitable for gene expression normalization in wheat were identified by a cross search for stable expression in Unigene and TIGR databases. The expression of 32 genes representing different functional classes was assessed by qRT-PCR of RNAs from 18 wheat tissues and floral organs and from seedlings exposed to low and high temperatures. Additionally, we carried out a comprehensive evaluation of the expression patterns of the actin and alphatubulin gene families, whose gene members have commonly been used as controls for the normalization of gene expression in wheat. The expression stability of the analysed genes was estimated by three statistical approaches: 1) threshold cycle (C t ) variation range and coefficient of variation; 2) geNorm and 3) Norm Finder programs. Finally, by studying the expression of a gene encoding PDI (Protein Disulfide Isomerase) in different tissues and developmental stages of wheat it was shown that the correct evaluation of its expression relies on the adoption of suitable normalization genes and can be negatively affected by the use of genes with unstable expression.

In silico identification of reference gene candidates
The in silico approach was applied to identify genes with stable expression levels among different wheat tissues and included two subsequent steps. First the UniGene EST ProfileViewer resource was used to analyse the expression of 41,256 sequence clusters of UniGene T. aestivum database (build 52, released 11 April 2008). The ProfileViewer estimates the expression pattern of each Unigene cluster by relating the number of clustered ESTs to the total ESTs obtained from the cDNA libraries of ten different tissues: callus, cell culture, crown, flower, inflorescence, leaf, root, seed, sheath and stem; the gene transcription profile of each cluster in each tissue is reported as a value representing the number of transcripts per million (TPM). Three criteria were adopted for selecting wheat candidate reference genes from UniGene: 1) each UniGene cluster would contain at least 60 single ESTs; 2) each UniGene cluster would be expressed in all ten tissues; 3) the difference between the highest and the lowest expression value (TPMs) for each UniGene cluster would not be below 40% in the six most representative tissues: flower, inflorescence, leaf, root, seed and stem. Using such criteria, 177 out of the 42,256 sequence clusters were selected and mean value (Mv), variance (V), standard deviation (SD) and coefficient of variation (CV) were calculated on their expression values (TPM) for the six most representative tissues. An independent search for members of the αtubulin and actin gene families carried out in UniGene identified ten more clusters, five for each gene family.
The second step consisted in the identification of the TC sequence corresponding to each of the 187 (177+10) selected UniGene clusters; the longest sequence, including part or the whole coding region and the 3' UTR, was identified, downloaded and exploited to BLAST search the TIGR wheat gene index database (TaGI version 11). This database is in many ways complementary to other gene indexing databases and its clustering and high stringency approach offer a number of significant advantages over alternatives such as UniGene. First, assembly provides a high confidence consensus representing each transcript, whereas low quality, misclustered or chimeric sequences are identified and discarded; moreover, closely related but distinct transcript isoforms can be recognized when sufficient sequence extent is available. Second, the produced consensus sequence is generally longer than its individual ESTs, providing a resource that can be used more effectively for functional annotation and for designing specific primers for its expression analysis. The TIGR Gene Indices compare and cluster ESTs and annotated cDNA sequences when they meet the following criteria: I) at least 40 base pair match; II) identity higher than 94% in the overlapping region; III) maximum 30 base pairs of unmatched overhang. These clusters are then assembled into a con-sensus sequence using Paracel Transcript Assembler [38]. UniGene applies different parameters: EST sequences are linked in a cluster if there is a 50 bp overlapping (100% identity) within the 3' UTR. However the identified clusters are not further checked through the more stringent criteria of assembly process and consensus sequences, this explains why sometimes several closely related but not matching TCs were included into the same UniGene cluster. Therefore we selected only the TCs showing identity higher than 95% with the longest sequence representative of each UniGene cluster. Since the more stringent assembly process used by TIGR enhances the detection of splicing variants, in the allohexaploid genome of T. aestivum it most likely allows the identification of closely related but distinct transcripts of homoeologous genes. For each TC or group of TCs linked to the 187 selected UniGene clusters, we determined their frequency (number of hosting libraries/total number of cDNA libraries) to find out the TCs reported in a significant number of different cDNA libraries and, consequently, more stable and suitable to be included in a preliminary extended set of candidate reference genes. Out of the 250 cDNA libraries used to construct the TIGR Wheat Gene Index (release 11), TC frequency counting was restricted to the libraries (197) containing at least 100 ESTs, because the remaining 53 include a limited number of clones deriving from differential display and cDNA-AFLP analyses or from SSH libraries. Functional classification of each UniGene cluster and linked TC or TCs was performed on the basis of its rice orthologous genes and relative GO (Gene Onthology) annotations.
For temperature stress experiments seeds of T. aestivum cv CS were germinated on cotton pads 6 days at 18°C, 16 h light. Seedlings were transplanted into small pots with soil and grown for 20 days (until reaching the third leaf unfolding stage, Feekes scale 1.3) at 18°C with 16 h light at 65% relative humidity. White fluorescent and incandescent lighting was combined to provide a light intensity of 196 μmol m -2 s -1 . Temperature treatments were performed exposing pots to 33°C and 4°C for 48 h. Two biological replications (20 plants) were exposed to each treatment. The aerial parts of the plants were harvested just before the start of temperature treatments (control plants), after 24 h and 48 h of exposition to 33°C and 4°C, frozen in liquid nitrogen and kept at -80°C until use.

RNA isolation and cDNA synthesis
From most samples total RNA was extracted using the TRIzol reagent (Invitrogen) according to manufacturer's instructions, but adding an additional chloroform extraction. Total RNA from caryopses was isolated by a LiCl based method [39]. The resulting RNA was treated with Rnase-free DNase I (Promega) according to the manufacturer's protocol. Following digestion, nucleotides were removed from RNA using a G50 sepharose buffer exchange column (Amersham). Absence of genomic DNA contamination in DNase I-treated samples was checked by qRT-PCR of 0.125 μg of RNA template using a primer pair (5'-CTTATGCAATCCAATGATGG-3' and 5'-TCTAT-GGTTCTGAAGAGGACC-3') designed to amplify an intron sequence of a gene encoding PDI and located in chromosome 4A [40]. When a single DNase treatment did not completely remove interfering genomic DNA, a second DNase incubation was performed to eliminate any detectable DNA. RNA concentration and integrity were checked with a UV/VIS spectrophotometer Lambda 3B (Perkin Elmer) before and after DNase I digestion. Only RNA samples with 260/280 wavelength ratio between 1.9 and 2.1 and 260/230 wavelength ratio greater than 2.0 before and after DNase I digestion were used for cDNA synthesis. The quality of RNA samples was also assessed by electrophoresis on 1% formaldehyde agarose gels. First-strand cDNA was synthesized from 3 μg of total RNA by Expand™ Reverse Transcriptase (RT) (Roche) and diluted 1:5 before use in qRT-PCR assays.

RT-PCR analysis of cold-and heat-responsive genes
The effect of temperature stresses was checked by the expression analyses of the wheat genes wcor14 and TaHSP101B, which are specifically induced by cold and heat stresses, respectively [41][42][43][44]. Expression analyses were carried out by semi-quantitative RT-PCR because transcripts of both genes were not detected in plants grown at 18°C (control plants in our temperature experiments), but were immediately induced by heat (TaHSP101B) and cold (wcor14) stress treatments.
Proteins of the HSP100 family are involved in the process of thermotolerance acquisition, in particular in renaturing proteins denatured by heat. They are ATPase involved in assembly/disassembly of protein complexes such as the ATP-dependent dissolution of cytosolic or nuclear protein aggregates forming during heat stress [45]. In wheat three gene sequences of the HSP100 family have been isolated and characterised: TaHSP101A, TaHSP101B and TaHSP101C [43]. To study the effect of high temperature treatment we chose the gene TaHSP101B because it showed the highest induction by heat shock [43,44]. A search in the NCBI database identified two sequences of Triticum aestivum encoding TaHSP101B: HSP101 (Accession number AF083344) [46] and HSP101B (Accession number AF097363) [43]. The same search found also two additional sequences isolated from T. durum (tetraploid wheat) and encoding TaHSP101B: TdHSP101B-A (Accession number AJ970533) and TdHSP101B-B (AJ970534), corresponding to homoeologous genes located in the long arm of chromosomes 1A and 1B, respectively [44]. The primers used for RT-PCR analysis of TaHSP1001B (forward: 5'-GACGCAGCTGTCCAAGAT-3'; reverse: 5'-GCACCTGGATGAGGATGT-3') were designed in conserved regions of the four cDNA sequences and amplified a product of 161 bp.
First-strand cDNA was synthesized as described above and the PCR reactions were performed by the TripleMaster PCR system (Eppendorf) using 1 μl of the RT reaction. After initial denaturation at 94°C for 2 min, amplification conditions were 35 cycles each at 94°C for 30 s, 60°C for 1 min and 72°C for 3 min, followed by a final extension step at 72°C for 7 min. RT-PCR of the constitutive actin gene (UniGene cluster Ta54825), one of the most stable genes in temperature treatments (see below), was performed as control. Samples of the amplification products (3 μl) were collected after 28, 32 and 35 PCR cycles and analysed by electrophoresis on 2% agarose gel. Each RT-PCR experiment was independently repeated twice to test amplification reproducibility. The specificity of the amplicons was checked by sequencing of the PCR products in order to confirm that its sequence corresponded to the target gene.

Design and validation of qRT-PCR primers
Specific primer pairs were designed for 32 candidate reference genes and for five distinct members for each of the αtubulin and actin gene families on the basis of their TC sequences. When multiple TC sequences linked to a single UniGene cluster were identified, the primers were designed in their conserved regions. Primer design and optimization were done by the Beacon Designer 6 software (STRATAGENE) imposing the following stringent criteria: T M of 55°C ± 2°C, PCR amplicon length between 60 and 280 bp, primer length of 20 ± 2 nt, and 40 to 60% guanine-cytosine content. Primers were also designed at the 3' end region of each TC sequence to encompass all potential splice variants and to ensure equal RT efficiencies. The complete set of primer pairs and their amplicon length, along with the corresponding UniGene clusters and linked TC sequences, are listed in Additional file 2. Single sequence amplifications were performed using gene-specific primer pairs, whereas multiple genes of the same family were also amplified by primer pairs designed in conserved regions [TEF-1α (m), ADP-RF (m), Actin (m) and α-tubulin (m)]; m representing amplification of multiple genes). Each of the primer pairs TEF-1α (m) and ADP-RF factor (m) amplified the transcripts of two identified homologous genes (Ta659/Ta53964 and Ta45379/ Ta2291, respectively) within their gene families, whereas each of the primer pairs targeting conserved regions of the Actin (m) and α-tubulin (m) families amplified the transcripts of five genes. Preliminary qRT-PCR assays carried out on a pool of all available cDNAs using six different primer concentrations ranging from 50 to 300 nM showed that the optimal primer concentration was 150 nM, because it generated the lowest C t value, a sharp peak by melting curve analysis and absence of non-specific products or primer-dimer artefacts. The specificity of the amplicons was also checked by electrophoresis on 2% agarose gel and sequencing of the PCR products in order to confirm that the product sequence was the same as the target gene. Five-point standard curves of a 5-fold dilution series (1:1-1:625) from pooled cDNA were used for PCR efficiency calculation of each primer pair. The PCR efficiency (E) is given by the equation E = (10 [-1/m] -1) × 100 [34], where m is the slope of linear regression model fitted over log-transformed data of the input cDNA concentration versus C t values according to the linear equation y = m*log(x) + b.

qRT-PCR conditions
qRT-PCR analysis was performed using an Mx3000PTM real time PCR system with Brillant SYBR green QPCR mas-ter mix (STRATAGENE) according to manufacturer's protocols. qRT-PCR assays were carried out in 25 μl reaction volumes containg 1 μl of each diluted cDNA (1/5) and 150 nM forward and reverse primers. No template and RTminus controls were run to detect contamination, dimer formation and presence of genomic DNA. As described previously, for each primer pair PCR efficiency was calculated in each run from a pool of all available cDNAs. The thermal profile comprised three segments: 1) 95°C for 10 min; 2) 40 cycles of 30 s denaturation at 95°C, 1 min annealing at 55°C and 30 s extension at 72°C (amplification data collected at the end of each extension step); 3) dissociation curve consisting of 1 min incubation at 95°C, 30 s incubation at 55°C, a ramp up to 95°C. Two biological replicates, resulting from two different RNA extractions, RT and qRT-PCR reactions, were used for quantification analysis and three technical replicates were analysed for each biological replicate. Gene expression levels were recorded as C t values, which are inversely related to the initial DNA concentration. The built in amplification-based proprietary algorithm (STRATA-GENE) was used to set the fluorescence threshold value for each primer pair reaction. This algorithm determines the portion of the amplification plots where all data curves display an exponential increase in fluorescence and calculate the threshold value that minimizes the standard deviation in C t values for each replicate set at a point which falls within 5-60% of the fluorescence shift for all curves. The number of cycles (C t ) at which the amplification-corrected normalized fluorescence (dRn) for each reaction crossed the threshold value was exported to Excel (Microsoft) for further analyses.

Assessment of expression stability
The expression stability of the analysed genes was estimated by three different statistical approaches. The first one consisting of single factor ANOVA on raw C t values of each gene was performed using the Excel Analysis Tool-Pak; genes showing high variance between the analysed samples should be avoided as controls. In this approach the C t difference (C t max-C t min) and CV for each gene were also computed to allow for a straightforward assessment of the most suitable reference genes that had a narrow C t range over all tissues and stress conditions analysed.
In the second and third approaches the stabilities among the samples of the candidate reference genes were evaluated by the software programs geNorm, version 3.4 [13] (Visual Basic application tool for Microsoft Excel available on the Internet) and NormFinder [15] (a Microsoft Excel Add-in available on the Internet) according to the author's recommendations. For both programs raw C t values were transformed to relative quantities using the delta-C t formula Q = E ΔCt , where E is the efficiency of the primer pair used in the amplification of a particular gene and ΔC t is the difference between the sample with the lowest C t (highest expression) from the data set and the C t value of the sample in question.
The statistical algorithm GeNorm determines a measure of gene expression stability (M value) corresponding to the average pairwise variation of a single candidate reference gene to all other genes under investigation. The Mvalue measure relies on the principle that the expression ratio of two hypothetical ideal internal control genes is identical in all samples, regardless of the experimental conditions or treatments. In this way, variation of the expression ratios of two candidate reference genes reflects the fact that one (or both) of the genes is (are) not constantly expressed, with increasing variation in ratio corresponding to decreasing expression stability (M value). Genes with the lowest M values have the most stable expression and therefore would be selected as ideal reference genes. In iterative steps, genes with the lowest expression stability (i.e. the highest M value) are removed. A new M value for each of the remaining genes is calculated until only two genes remain. Because these calculations are based on ratios, the final 2 genes cannot be resolved.
NormFinder assesses the expression stability of a gene by evaluating its expression variation within tissues or treatments ("groups" in NormFinder terminology) compared to variation among tissues/treatments. On the basis of a given "group identifier" the program can in fact discriminate between different groups, e.g. stem, leaf, spike or seedliings exposed to low or high temperatures. The program algorithm implies the estimation of intra-and intergroup variation and combines both results in a stability value for each investigated gene. The candidate gene with the lowest stability value is the most stable gene within the groups studied. The best combination of two genes is also indicated. In our analyses each tissue and temperature treatment is considered an experimental group, so we had 24 groups (18 different tissues and 6 temperature treatments, including controls and 24 and 48 h at 4°C and 33°C), each composed of the mean values of the two biological replicates.

Reference gene validation
To determine how the adoption of different reference genes can affect the normalization of the expression data for a gene of interest, the same 24 plant cDNA samples used for the stability analyses of reference genes (see Plant material) were also analysed by qRT-PCR for the expression of TaPDIL1-1 (a gene encoding the Protein Disulfide Isomerase). Isolation and characterization in wheat of the three homoeologous gene sequences encoding PDI (TaPDIL1-1) and of their full-length transcripts have been reported previously [40]. The primer pair used for TaPDIL1-1 qRT-PCR analysis (forward: 5'-CGTGGTCT-TCAAATCCTG-3'; reverse 5'-GTAACCCTGGACAT-CAAAC-3') was designed in conserved regions of the three homoeologous cDNA sequences at their 3' ends (Accession numbers: AJ868105, AJ868106, AJ868107) with annealing temperature of 55°C. The PCR efficiency was 100.35 ± 0,365 with a coefficient of determination (R 2 ) of 0,998; the amplicon size was 187 bp. The TaPDIL1-1 expressions were normalized using five different strategies: 1) geometric average of the three references genes selected as the most stable genes by geNorm; 2) geometric average of the two references genes selected as the most stable genes by NormFinder; 3) the single most stable gene identified by NormFinder; 4-5) the two genes related to the Unigene clusters Ta54825 (Actin) and Ta25534 (α-tubulin) used alone. Normalized TaPDIL1-1 relative values are given as mean value ± SD. Standard deviations on normalized expression levels were computed according to the geNorm user manual [47].

Identification of reference gene candidates
As described in the section "Methods", the approach adopted for the identification of wheat reference candidate genes included two subsequent steps. First the expression of 41,256 sequence clusters of UniGene T. aestivum database was analysed using the EST ProfileViewer resource (Additional file 3 reports the 177 selected Uni-Gene clusters sorted in ascending order on the basis of their TPM CVs). Since UniGene clusters can not generate contigs and/or consensus sequences, the second step consisted in the identification of TC or group of TCs corresponding to each cluster using the TIGR wheat gene index database (TaGI version 11), whose outputs are more effective for gene functional annotation and for designing qRT-PCR primers. On the basis of the frequency (number of hosting libraries/total number of cDNA libraries) of each TC or group of TCs linked to the 177 selected UniGene clusters it was possible to pinpoint the TCs which, being represented in a significant number of cDNA libraries from different tissues, were the most suitable candidate reference genes to include in a preliminary screening. Functional classification of single TCs or groups of TCs linked to UniGene clusters was performed on the basis of their homology to rice orthologous genes and relative GO (Gene Onthology) annotations (Additional file 4). Of the 177 initially selected TCs or groups of TCs, 34 (19.21%) were homologous to rice genes encoding putative proteins whose biological functions are either unknown or not clear (Table 1); based on their presumed biological functions (GO biological process annotations), the remaining 143 selected sequences related to characterized rice genes were assigned to ten functional categories ( Table 1). The categories "Cellular organization" and "Metabolism" contained the highest number of genes (18.08% and 17.51%, respectively), followed by "Posttranscriptional modification" (15.82%), "Protein synthesis" (9.04%) and "Signal transduction" (7.34%). These five functional categories were also reported among the most represented in sets of candidate human housekeeping or maintenance genes selected through microarray analysis using different Affymetrix platforms [48][49][50]. The remaining five categories ("RNA processing", "Energy", "Transcription", "Biotic and abiotic stress related" and "Development") were under-represented in the 177 selected sequences, accounting for only 12.99% (Table 1). The GO terms relative to biological processes over-represented among the different functional categories were determined to gain insight on the functional roles of the 177 selected sequences (Table 1). They included genes mediating a variety of basic cellular functions, such as: protein transport, ubiquitin-dependent protein degradation, proteolysis, protein translation, primary metabolism (amino acid biosynthesis, glycolysis and tricarboxylic acid cycle). Moreover, two of the over-represented GO terms referred to genes involved in the regulation of transcription (transcription factors) and in the intracellular signalling cascade (protein kinases). On the basis of the number of EST sequences included in each UniGene cluster the selected sequences were also assigned to four classes according to their level of expression: 1) high (>800 ESTs, 8.47%); 2) moderate (<800/>300 ESTs, 18.08%); 3) low (<300/>100 ESTs, 47.46%); 4) very low (<100 ESTs, 25.99%) (Additional file 5).
Most sequences identified as reference gene candidates were new, but they included also most of the known HKGs used as internal controls in gene expression analyses, such as GAPDH (Ta30768), malate dehydrogenase (Ta53937), β-tubulin (Ta44405), histone (Ta38797), cyclophylin (Ta35497), ribosomal proteins (Ta55103, Ta27771, Ta55901, Ta54208, Ta55238, Ta54141), members of the translation initiation and elongation factors gene families (Ta15761, Ta746, Ta14227, Ta1092, Ta54575, Ta53964, Ta659), ubiquitins (Ta50503, Ta55221), ubiquitin-conjugating enzymes (Ta56067, Ta54238, Ta24654) and proteasome subunits (Ta55216, Ta21374, Ta20028, Ta22845, Ta55217, Ta54768, Ta14021) (Additional file 4). Surprisingly, none of the genes belonging to the families encoding actin and αtubulin met the three criteria adopted for selecting wheat reference candidate genes from UniGene. Since these genes have commonly been used as controls for normalization of gene expression in wheat, we decided to perform a detailed search for members of these two gene families in UniGene simply using the terms "actin" or "alphatubulin", and "wheat" combined by the Boolean operator "AND". Out of the 35 and 13 UniGene clusters identified for actin and α-tubulin, respectively, using the UniGene EST ProfileViewer resource we selected those containing at least 60 single ESTs and whose expression was detected at least in the six most representative wheat tissues: flower, inflorescence, leaf, root, seed and stem. Adopting such criteria five distinct UniGene clusters were identified for each of the two gene families and the Mv, V, SD and CV were NS = Number of sequences assigned to a GO term for biological process %FC = Percentage of sequences assigned to a GO term for biological process relative to a particular functional class %TS = Percentage of sequences assigned to a GO term for biological process relative to the total (177) selected sequences.
calculated on their expression values (TPMs) for the six most representative tissues (Additional file 6). The CV values of three members of actin and five of α-tubulin gene families were higher than those of all the 177 UniGene clusters selected on the basis of the three stringent criteria adopted for selecting wheat reference candidate genes (compare Additional files 3 and 6), whereas the two lowest CV values (37.25 and 38.38) obtained for the two actin UniGene clusters Ta54825 and Ta53908 performed better than those of only 18 and 14, respectively, of the 177 selected sequences. Additionally, in all the members of actin and α-tubulin gene families the difference between the highest and lowest expression value (TPM) detected in the six most representative wheat tissues was inferior to 40% (Additional file 6). However, we decided to include all the 10 sequences in the qRT-PCR analysis to select the most stable member for each gene family and to compare them with all the other candidates selected as reference genes. For each UniGene cluster we identified the corresponding TC or group of TCs and determined their frequency (number of hosting libraries/total number of cDNA libraries). The five UniGene clusters and related TCs identified for the α-tubulin corresponded to the five non-homoeologous gene family members recently cloned in wheat [36]. Table 2 lists the 32 UniGene clusters and corresponding TCs selected for the analysis of their expression by qRT-PCR of RNA from 18 different wheat tissues and developmental stages and from seedlings exposed to low and high temperatures. These gene sequences were chosen considering three main criteria: 1) the ranking of the 177 reference gene candidates on the basis of the CV of their expression values (TPMs) in six wheat representative tissues (UniGene Profile Viewer); 2) the ranking of the TCs or groups of TCs linked to each of 177 UniGene clusters on the basis of their frequency (number of hosting libraries/total number of cDNA libraries); 3) their level of expression estimated as the number of EST sequences comprised in each UniGene cluster. Because one aim of the present research was to compare the effectiveness of traditional and novel reference genes, we included some of the most common HKGs. Finally, special attention was paid in order to select genes involved in distinct biological processes, this would theoretically reduce the chance of including co-regulated genes. However, for the purpose of testing the expression stability of sequences belonging to the same gene family, we also included multiple genes encoding ADP-ribosylation factors and elongation factors with very different expression levels ( Table 2). Out of the 32 sequences selected as potential reference genes for wheat, 12 were related to well-known HKGs representing different functional classes, whereas the remaining 20 (in italics in Table 2) included genes that, according to literature searches, were not related to any housekeeping activity or genes with unknown functions, hence they can be considered novel with respect to the normalization issue. The first 11 Unigene clusters and corresponding TCs reported in Table 2 were selected because they ranked among the top 25 sequences with the lowest CV for their expression values (TPM), whereas 16 additional clusters (12-27 in Table 2) were chosen because they ranked among the 25 top sequences with the highest frequencies.

Selection of the best candidate reference genes for qRT-PCR analyses
Three more UniGene clusters (28-30 in Table 2) were selected for their known or suspected housekeeping roles in basic cellular processes. Eventually, the UniGene cluster Ta55512 (31 in Table 2) was included because it was related to a cathepsin B-like cysteine protease, the only chosen gene assigned to the GO term for the biological process "Proteolys", whereas Ta2291, encoding for an ADP-ribosylation factor (32 in Table 2), was selected because it was homologous to Ta45379, but had a different expression level (18 in Table 2). The selected sequences covered a wide range of expression levels: 11 high (>800 ESTs in each UniGene cluster), 12 moderate (<800/>300 ESTs), 8 low (<300/>100 ESTs) and 1 very low (<100 ESTs) ( Table 2).

Specificity and primer efficiency of the qRT-PCR reactions
The expression stabilities of the 32 selected candidate reference genes and of the five distinct homologous genes for each α-tubulin and actin gene families were assessed by qRT-PCR in a set of 24 tissue samples subdivided into two subsets: the first one included 18 samples representing different tissues and developmental stages of wheat, the second one included six samples and consisted of two temperature treatments (4°C and 33°C), each for 24 and 48 h, and controls. The effect of temperature stresses was checked by RT-PCR analyses of the cold-and heat-responsive genes wcor14 and TaHSP101B (Additional file 7). Accumulation of wcor14 transcripts was undetectable in seedlings grown at 18°C (control samples), whereas strong induction was detected after 24 h of exposition to 4°C, which increased further after 48 h. These results are in perfect agreement with those obtained by Tsetanov et al. [41], who showed by northern analyses of RNA from 23 day old seedlings of T. aestivum that transcripts of wcor14 accumulated within 3-6 h of cold acclimation at 4°C and reached a maximum after three days. Exposure to 33°C for 24 and 48 h did not induce detectable amount of wcor14 transcripts (Additional file 7), confirming its specific induction by low temperatures [41]. In agreement with the results obtained by Campbell et al. [43] and Gulli et al. [44], who analysed the expression of TaHSP101B in bread and durum wheats, respectively, our RT-PCR analyses showed that transcripts of these gene were not present in control seedlings grown at 18°C, strong expression was detected after exposition to 33°C for 24 h, then it  7). The absence of TaHSP101B transcripts in seedlings exposed to cold temperature (Additional file 7) confirmed that its expression was specifically induced by heat shock.
Total RNA was isolated from the 24 samples, DNase I digested and reverse transcribed.
Since we chose a two-step qRT-PCR protocol, reverse transcription and PCR-mediated cDNA amplification were carried out by subsequent steps in separate tubes. The two-step protocol was preferred because it would reduce unwanted dimer formation between primers when SYBR Green is used as a detection dye [51]. The RT reaction was primed by oligo-(dT), rather than using random-sequence primers which anneal preferentially to abundant mRNA species. This choice was critical because, as deduced on the basis of the number of EST sequences comprised in the corresponding UniGene cluster, some of the selected candidate reference genes would be expressed at low level. PCR primer pairs for the 32 selected reference candidate genes and for the five members for each of the α-tubulin and actin gene families were designed and preliminarily tested in qRT-PCR reactions from a pool of all available cDNAs (see Methods). In order to avoid unreliable amplification of splice variants and assure uniform RT efficiencies, typically primer pairs were designed to target an amplicon located within the 3' end region of the transcript of interest. For each biological replicate and tissue sample, the same cDNA pool was used for qRT-PCR amplification of the 42 (32+10) analysed genes using the previously tested gene-specific primers. Moreover, qRT-PCR analysis of every gene was performed in triplicate for each of the 24 cDNA pools, along with no template and RT-minus controls. The 42 primer pairs used to amplify the candidate reference genes generated single amplicons of the expected size from the various cDNA pools, as shown by the presence of single bands in agarose gel electrophoresis and by single-peak melting curves of the PCR products obtained after 40 amplification cycles ( Fig. 1 and Additional file 8). A detailed analysis of the dissociation curves confirmed the absence of primer dimers and of other products resulting from non-specific amplification (Additional file 8). A more stringent test of PCR specificity was performed by sequencing the amplification products of each of the 42 candidate reference genes. Always the sequence of the PCR product matched that of the target cDNA, thereby confirming the exquisite PCR specificity of the developed primer pairs.
In addition to the primer pairs targeting single transcripts, we set out the evaluation of a further subset of primer pairs, designed in conserved regions identified by sequence comparison, to amplify in a single PCR reaction the transcripts of multiple members of some gene families. Primer pairs designed as TEF-1α (m), ADP-RF (m), Actin (m) and α-tubulin (m), where (m) denotes multiple amplifications, targeted conserved sequences within the transcripts of at least two homologous genes for each family. The simultaneous amplification by conserved primer pairs of sequences transcribed from multiple homologous genes with balanced expression levels has been proposed as a suitable strategy for qRT-PCR normalization. Such an approach has been exploited for the Arabidopsis genes AtACT2 and AtACT8, which display complementary patterns of expression, making their combined expression quasi-constitutive [52]. Gel electrophoresis and meltingcurve analysis showed that each of the four conserved primer pairs designed to amplify transcripts from multiple genes of the four analysed families of wheat amplified a PCR product of the expected size from the various cDNA pools (Fig. 2). Single PCR products purified from gels were cloned into the pGEM plasmid vector (Promega) and ten different clones were sequenced for each of the four multiple amplifications. Two sequences of 276 bp and with identity of 90.58% were observed among the ten clones obtained by the PCR with the primer pair ADP-RF (m); also among the ten clones produced by the amplification using the primer pair TEF-1α (m) there were two different sequences of 233 bp with 95.28% identity. In both cases the sequences corresponded to those used for designing the primer pairs within the conserved regions, Ta45379/Ta2291 for ADP-RF and Ta659/Ta53964 for TEF-1α. In each group of ten clones obtained by the primer pairs Actin (m) and α-tubulin (m) three similar homologous sequences were found, which matched three of the five members of their respective gene families (Ta54825/Ta53908/Ta20863 for Actin and Ta53981/ Ta33558/Ta54519 for α-tubulin), showing that multiple amplifications were produced also by these conserved primers.
Although fluorescence intensity at a given stage during qRT-PCR is primarily determined by the initial concentration of the cDNA template, its increase is strongly affected by the amplification efficiency. When the efficiency is 100%, the amount of a cDNA amplicon targeted by a given primer pair is doubled at every PCR cycle during the exponential phase; however, several factors can significantly affect the efficiency of individual PCR reactions. In our study the PCR efficiency of each primer pair in each of the two biological replicates was determined by five-point standard curves of a 5-fold dilution series (1:1-1:625) from pooled cDNA. By plotting the C t s against the logarithm of the relative amounts of template (serial dilutions from pooled cDNA) a standard curve was generated by the Moreover, the Mx3000PTM software system computed the R 2 value (Coefficient of determination), which represents an indicator of the quality of the fit of the standard curve to the plotted data points. The R 2 values calculated for all the standard curves in each of the two biological replicates were higher than 0.994, indicating a strong linear relationship between the detected C t values and the correspondent relative amount of template in all the amplification reactions (Additional file 9). The efficiencies of all the 46 primer pairs were comprised between 94% and 110%, about 44% (29/46) of them with values close to 100% (97-103%) (Additional file 9), reflecting the high quality of the amplification reactions. Generally the standard deviation of primer efficiencies was low (<1.4), indicating comparable amplification efficiencies between the two biological replicates, only 14 primer pairs showed values higher than the 2% (SD > 1. 4). The C t s of the three replicates and the efficiency values (Additional file 9) were used to compute the relative expression level of the selected candidate reference genes for each of the two biological replicates.

Selection of the best candidate reference gene within the same family
In order to select the most stable gene encoding actin, αtubulin, ADP-ribosylation factor and translation elongation factor, the expression stability of genes belonging to each family was analyzed by single factor ANOVA on the C t values of five different data-sets: 1) all 24 samples including tissues, developmental stages and temperature treatments; 2) 18 samples relative to different tissues and developmental stages; 3) six samples consisting of two temperature treatments (4°C and 33°C), each for 24 and 48 h and their controls; 4) six samples relative to the single floral organs from fully emerged spikes; 5) six samples referring to vegetative tissues and developmental stages (shoots, stems and leaves). These data sets were used to evaluate the expression stability of the genes belonging to the four families by 18 primer pairs, 14 of them targeting specific gene transcripts of the gene families (five each for actin and α-tubulin and two each for ADP-ribosylation factor and translation elongation factor) and four primer pairs which amplified multiple transcripts for each of the four gene families [Actin (m), α-tubulin (m), ADP-RF (m) and TEF-1α (m)]. Table 3 reports the CV values and the C t differences (ΔC t = C t max-C t min) obtained by the 18 primer pairs for each of the five data-sets considered, whereas the complete statistics obtained by ANOVA analysis are reported in Additional file 10. The CV of the C t values gives a straightforward indication of the expression stability of a particular gene; most genes with the lowest CV values showed also the lowest values of C t range (ΔC t ) and MST (Mean Square of Treatments) (Additional file 10). On the other hand, the very low values of the MSE (Mean Square of Error) for all the analysed genes, showed that the C t values obtained in different replicates were highly reliable and the qRT-PCR experiments generate highly reproducible results (Additional file 10). Table 3 and Additional file 10 show that the most stable genes (lowest CV values) for each family were: the gene corresponding to the UniGene cluster Ta25534 for αtubulin, Ta54825 for actin, Ta2291 for ADP-ribosylation factor and Ta53964 for translation elongation factor. These four genes were chosen as representative of the corresponding families and used for further comparison with the other candidate reference genes.

Expression levels of the selected candidate reference genes
In order to avoid possible errors due to the interaction of co-regulated genes, further analyses were carried out on 32 genes involved in distinct biological processes and metabolic pathways; they comprised 28 genes selected on the basis of their uniform expression detected in UniGene and TIGR and a single gene for each of the four gene families. The expression levels could be grouped into two arbitrary categories: 19 highly expressed genes, whose mean C t values were below the general mean (C t = 22.31 cycles) and 13 lowly expressed genes, with mean C t values above this value (Fig. 3). The highest expression was shown by Ta38797 (Histone), Ta53919 (S-adenosylmethionine decarboxylase), Ta50503 (Ubiquitin) and Ta54825 (Actin), with mean C t values between 20 and 21 cycles, the lowest by Ta25534 (α-tubulin), Ta2776 (RNase L inhibitor-like protein), Ta54448 and Ta54733 (Proteins of unknown function), and Ta44405 (β-tubulin), with mean C t values between 26 and 24 cycles (Fig. 3). On the basis of the number of EST sequences included in each Unigene cluster, the 19 reference genes with low C t values had been previously classified either as highly expressed (9 genes with more than 800 EST) and moderately expressed (10 genes with more than 300 EST). Of the 13 reference genes with low C t values, only three including more than 300 ESTs (Ta54963, Ta54447 and 54948) were considered moderately expressed in UniGene, whereas the expression level of the remaining 10 genes was low (<300/>100 ESTs) or very low (<100 ESTs). In our opinion, these data indicate a substantial agreement between the gene expression levels evaluated by qRT-PCR analyses and those estimated by the number of ESTs included in each corresponding UniGene cluster. Ta54963 (GABAreceptor associated protein), Ta54238 (Ubiquitin-conjugating enzyme), Ta35497 (Cyclophilin), Ta53967 (Vacuolar ATP synthase 16 kDa proteolipid subunit) exhibited the lowest gene expression variation (less than 3 cycles), while Ta44405 (β-tubulin), Ta38797 (Histone) and Ta1698 (NADP-isocitrate dehydrogenase) showed the most variable expression (more than 5 cycles).

Expression stability of the selected candidate reference genes
As discussed in the previous paragraph, it was evident that none of the 32 selected candidate reference genes showed a constant expression under all the tested conditions (tissues, developmental stages and temperature stress) in wheat. Hence, as expected, it was not possible to find an all purpose reference gene, and it would be necessary to select the most suitable one for normalizing gene expression considering the peculiar and specific experimental conditions. The relationships between expression stability of the 32 candidate genes and different types of tissues and/or experimental conditions were tested by running the programs geNorm [13] and NormFinder [15] using the five data-sets previously described. genes. The ranking order of the following eight most stable genes, with M values comprised between 0.39 and 0.55, was: Ta54227 (3 rd , Cell division control protein, AAA-superfamily of ATPases), Ta2291 (4 th , ADP-ribosylation factor), Ta53899 (5 th , Hypothetical protein), Ta4045 (6 th , Ubiquinol-cytochrome C reductase iron-sulfur subunit), Ta2776 (7 th , RNase L inhibitor-like protein), Ta53967 (8 th , Vacuolar ATP synthase 16 kDa subunit), Ta54963 (9 th , GABA-receptor-associated protein) and Ta53919 (10 th , S-adenosylmethionine decarboxylase) (Fig. 4A). The ranking remained very similar when the set of 18 samples, including different tissues and developmental stages, was analysed, with the lowest M value (0.33) for Ta35284 and Ta22845 ( Figure A in additional  file 11). The rank order of the following seven most stable genes was the same obtained analysing the whole data set of 24 samples (compare figures 4A and A of additional file 11), and the least stable gene with the highest M value (0.87) was Ta44405 (β-tubulin). In the temperature treatment series Ta54825 (Actin) and Ta53889 (Hypothetical protein), both with M = 0.071, were the most stable genes, followed by Ta22845 (3 rd , ATP-dependent 26S proteasome regulatory subunit), Ta2291 (4 th , ADP-ribosylation factor), Ta54238 (5 th , Ubiquitin-conjugating enzyme), Ta54948 (6 th , Rab GTPase homolog, Rab7 subfamily), Ta2776 (7 th , RNase L inhibitor-like protein), Ta30797 (8 th , Phosphogluconate dehydrogenase), Ta54227 (9 th , Cell division control protein, AAA-superfamily of ATPases) and Ta4045 (10 th , Ubiquinol-cytochrome C reductase iron-sulfur subunit) (Fig. B in additional file  11). Seven of the top ten most stable candidate reference genes in the data set of 18 samples, (Ta53889, Ta22845, Ta2291, Ta54948, Ta2776, Ta54227 and Ta4045) were also classified among the ten most stable genes in the temperature treatment series, where the two least stable genes were Ta53891 (Zn-finger, A20-like domain containing protein) and Ta44405 (β-tubulin) (Fig. B in    value 0.09), followed by Ta53899 (3 rd , Hypothetical protein), Ta54227 (4 th , Cell division control protein, AAAsuperfamily of ATPases), Ta54238 (5 th , Ubiquitin-conjugating enzyme E2), Ta53891 (6 th , Zn-finger, A20-like domain containing protein), Ta54447 (7 th , Acetyl-CoA acyltransferase), Ta2291 (8 th , ADP-ribosylation factor), Ta53919 (9 th , S-adenosylmethionine decarboxylase) and Ta55512 (10 th , Cathepsin B-like cysteine protease), whereas Ta1698 (NADP-isocitrate dehydrogenase) and Ta44405 (β-tubulin), with M values of 0.63 and 0.59, were the least stable genes (Fig. D in additional file 11).

Unigene cluster
The high expression stability of a gene indicates that the use of a single internal control is appropriate, but some studies require the normalization based on two or more stable internal control genes. The geNorm program, in addition to the gene stability measure M, computes a normalization factor (NF) and assesses the optimal number of reference genes for generating that factor. The pairwise variations Vn/Vn+1 between two sequential normalization factors (NFn and NFn+1) are used to determine the necessity of adding the next more stable control gene for reliable normalization. Pairwise variations were calculated using geNorm for each data set to determine the optimal number of internal control genes for normalization ( Fig. 4B and Additional file 12). In all the five data sets analysed the two most stable candidate reference genes yielded V values lower than 1.5. That cut-off was suggested as the limit beneath which it would not be necessary to include additional reference genes for normalization [13]. Since the V values did not significantly decrease when more than four reference genes were included ( Fig. 4B and Additional file 12), it can be concluded that normalization using either the two or three most stable reference genes would be an adequate normalization criterion for gene profiling studies in the five data sets considered.
NormFinder (the second adopted program based on a model approach) ranks the best candidate references genes according to their minimal combined inter-and intra-group variation of expression. The results obtained by NormFinder analysis for the five data sets are shown in Table 4. When the complete data set of 24 samples was considered, Ta54227 (Cell division control protein, AAAsuperfamily of ATPases) was identified as the most stable gene, with a stability value of 0.035, followed by Ta2291 (2 nd , ADP-ribosylation factor), Ta2776 (3 rd , RNase L inhibitor-like protein), Ta35284 (4 th , Protein transport protein Sec23A), Ta53919 (5 th , S-adenosylmethionine decarboxylase), Ta53967 (6 th , Vacuolar ATP synthase 16 kDa proteolipid subunit), Ta54171 (7 th , Superoxide dismutase), Ta54963 (8 th , GABA-receptor-associated protein), Ta54733 (9 th , Protein of unknown function, DUF 866) and Ta4045 (10 th , Ubiquinol-cytochrome C reductase iron-sulfur subunit) (stability values between 0.049 and 0.069); the two least stable genes were Ta44405 (βtubulin) and Ta38797 (Histone), with stability values of 0.150 and 0.198 respectively, the same identified as the worst genes by geNorm. The best combination of two genes suggested by NormFinder was that of Ta54227 and Ta2291, which improved the stability value to 0.030, indicating a more reliable normalization than that based on single genes. Similar results were obtained for 18 samples Expression level of 32 candidate reference genes in 24 tissues, developmental stages and seedlings exposed to high and low temperatures Figure 3 Expression level of 32 candidate reference genes in 24 tissues, developmental stages and seedlings exposed to high and low temperatures. The expression levels (C t values) of each candidate reference gene (UniGene cluster) are shown as average (diagonal cross), median (horizontal lines), 25 th to 75 th percentile (boxes) and expression ranges (whiskers) for the 24 cDNA pools analysed. The line at C t 22.31, corresponding to the general mean, discriminates the reference genes with high and low expression.
including tissues/organs at different developmental stages of wheat, with the lowest stability value (0.032) for Ta54227 (Table 4). In the temperature treatment series NormFinder selected Ta30797 (Phosphogluconate dehydrogenase), with a stability value of 0.030, as the most stable single gene, whereas, as in geNorm analysis, Ta44405 (β-tubulin) and Ta53891 (Zn-finger, A20-like domain containing protein) were identified as the two least stable genes ( Table 4). The best combination of two genes was that of Ta30797 and Ta4045, with a stability value of 0.023 (Table 4). In the six samples relative to the single floral organs NormFinder identified as the best combination of two genes that of Ta30797 and Ta54227, with stability value of 0.045 which was significantly lower than that (0.061) of the most stable gene (Ta30797) taken alone. Finally, in the six samples including only vegetative tissues Ta54227 showed again the lowest stability value (0.022), the same value calculated for the best combination of two genes (Ta54227 and Ta53889).
The stability ranking of the reference genes across the five analysed data sets obtained by geNorm and NormFinder was not identical for the genes with the most stable expression (Table 5). In particular, the best normalization genes selected by the two programs were not the same, but there was substantial agreement when the groups of genes with the most and least stable expression were considered. For instance, in the whole data set of 24 samples (I Data set in Table 5), the most stable genes identified by NormFinder (Ta54227 and Ta2291) were classified as the third and fourth most stable genes by geNorm, and out of the top ten most stable reference genes determined by NormFinder, eight were also classified among the ten most stable genes by geNorm (Table 5). Moreover, the ranking of the eleven least stable genes determined by NormFinder was the same of geNorm (Table 5). Similar relationships among the two methods were also observed by comparing the ranking of the 32 candidate reference genes in the remaining four data sets ( Table 5).

The choice of the reference genes affects the normalization of a gene of interest
The expression levels of the target gene TaPDIL1-1 were used as an example to show the effect of using different normalization genes on estimating its relative expression. TaPDIL1-1 encodes the typical PDI (Protein Disulfide Isomerase) which accomplishes several metabolic functions, the most important consisting in disulfide bond formation and isomerization during the folding of secretory proteins [39,40]. Previous Northern and RT-PCR analyses had shown the constitutive expression of PDI in several wheat tissues, but a much higher expression during the early stages of wheat seed development [39,40]. The transcript quantity of TaPDIL1-1 was determined by qRT-PCR in the same 24 cDNA samples used for selection of refer-ence genes. Moreover the relative expression level of the gene of interest was calculated in the samples of four data sets: 1) 18 different tissues and developmental stages of wheat; 2) six vegetative tissues and developmental stages (shoots, stems and leaves); 3) six floral organs from fully emerged spikes; 4) six samples consisting of two temperature treatments (4°C and 33°C), each for 24 and 48 h, and their controls. The TaPDIL1-1 expression was normalized in each of the four data sets following five different methods: i) geometric average of the three genes selected by geNorm; ii) geometric average of the two genes selected by NormFinder; iii) top-ranking genes identified by NormFinder; iv) gene related to the Unigene cluster Ta54825 (Actin); v) gene related to the Unigene cluster Ta25534 (α-tubulin). Since members the actin and αtubulin gene families have been commonly used as controls for the normalization of gene expression in wheat, Ta54825 and Ta25534 were used as examples of the consequences of using unstable reference genes.
Normalization of the expression of TaPDIL-1-1 using the reference genes selected by the three methods based on geNorm and NormFinder programmes did not detect significant differences in the relative expression of the gene of interest in all the four analysed data sets (Figs. 5 and 6). As expected, the mRNA of TaPDIL1-1 was constitutively expressed in all analysed tissues, with very high transcription at the early stage of seed development (Seeds-1) and moderately high expression also in the middle and late stages of spike development (Spikes 2-3) and in floral organs such as lodicules and pistil. The lowest transcription was detected in vegetative tissues, such as mature leaves and stems, and in the non-reproductive floral organs (glumes, lemma and palea) (Fig. 5). On the contrary, normalization of the TaPDIL1-1 transcription levels by the commonly used reference genes Ta54825 (Actin) and Ta25534 (α-tubulin) led to significant over-or underestimation of transcription of the gene of interest (Fig. 5). Normalization using Ta25534 (α-tubulin) of the TaPDIL-1 expression in spikes collected at different stages of development, for instance, led to an estimation of its transcription level about 2-6 times higher than those obtained using the most stable genes indicated by geNorm and NormFinder (Fig. 5). The misleading effect of using unsuitable reference genes for the normalization of expression data was even more evident when the analysis was performed on the data sets relative to single groups of tissues (vegetative or floral organs, Figs. 6A and 6B) or temperature treatments (Fig. 6C). In the data set including only the vegetative tissues normalization with Ta54825 (Actin) estimated a significantly lower expression level of TaPDIL1-1 (approximately 1.3-3.6 times) in four of the six tissues analysed (Shoots1, Shoots3, Shoots4 and leaves) (Fig. 6A). Moreover, the normalization procedure with Ta25534 (α-tubulin) showed a weak but significant geNorm output charts of M and V values for the 32 selected reference genes in 24 samples  Fig. 6A). Among the floral organs, the use of Ta54825 (Actin) and Ta25534 (α-tubulin) as single reference genes led to significant over-estimation of TaPDIL1-1 expression in lodicules and pistil (Fig. 6B). Moreover, when TaPDIL1-1 expression was normalized with Ta54825 (Actin) in stamens its transcription level was about 3.7 times higher than the estimation obtained using the best genes indicated by geNorm and NormFinder (Fig. 6B). Finally, when the normalization of the TaPDIL1-1 gene was carried out using the geometric mean of the reference genes selected by NormFinder and by geNorm, the single best gene identified by NormFinder and Ta54825 (Actin), the comparison of seedlings exposed for 48 h at low (4°C) and high (33°C) temperatures with the control samples exposed to 18°C detected a significant up-regulation of TaPDIL1-1 mRNA (Fig. 6C). It is noteworthy that Ta54825 (Actin) was one of the two most stable genes selected by geNorm in the temperature treatments (Table  3). On the contrary, the normalization of TaPDIL1-1 expression data using Ta25534 (α-tubulin) as reference gene resulted in a significant decrement of its transcription level after 24 and 48 h exposure to low and high temperatures (Fig. 6C). These contrasting results can be explained by the significant up-regulation of Ta25534 (αtubulin) induced by the cold and heat shocks.
Among the 177 in silico identified sequences, specific and efficient primers were developed for 32 genes, whose expression level was assessed by qRT-PCR using a set of cDNAs from 24 plant samples, which included different tissues, developmental stages and temperature stresses. The 32 selected sequences consisted of 12 well-known HKGs representing distinct functional classes and 20 genes that either were not related to any housekeeping activity or had unknown functions, therefore they can be assumed as novel in relation to the normalization issue. The UniGene clusters Ta54825, Ta25534, Ta2291 and Ta53964 were included among the 32 selected candidate reference genes because preliminary qRT-PCR analyses showed that they were the most stable of the genes belonging to four analysed gene families: actin, α-tubulin, ADP-ribosylation factor and translation elongation factor, respectively. The simultaneous amplification of two or more reference genes of the same family by a single primer pair designed in the conserved regions has been suggested by some Authors as a strategy for reducing the cost of using multiple and independent control genes [16,20,52]. This approach assumes the balanced expression of different amplified genes of the same family, but this is not always true. For example, UBQ5 (UBIQUITIN 5) has been one of the most suitable reference genes in a set of tissue samples in rice, whereas the expression of UBQ10 was unstable [19]. A similar situation has been observed for the actin gene family in a set of samples collected at different development stages of soybean, where ACT2/7 was stably expressed, whereas ACT11 showed variable profiling [53]. In wheat and in our experimental conditions the approach based on the use of conserved primer pairs amplifying multiple gene sequences was not successful, in fact the primer pairs amplifying the single best gene for each of the four analysed families performed better than those targeting multiple genes (Table 3).
Since the search for suitable reference genes is both timeconsuming and cost intensive, several computer programs based on different statistical approaches have been developed to assist the selection of a reference gene among several HKGs. The software geNorm has been reported in 2002 [13] and has been followed by NormFinder in 2004 [15], they have been very successful and widely used. Also BestKeeper is a freely available program that determines the best-suited reference gene by pair-wise correlation analysis of raw C t values [14]. According to Lyng et al. [54], this approach may be useful to narrow down a search if no specific genes are known to be plausible candidates, whereas more advanced statistics, as those provided by geNorm and NormFinder, are needed to rank the genes if several of them are identified as good candidates; moreover BestKeeper can analyse a maximum of ten reference genes. Some more statistical models have been developed for identifying optimal reference genes, but either they are not freely available, their application is complex [55,56], or they do not take the PCR efficiency into account [16]. The choice of the program using the most suitable statistical method represents a crucial but difficult task, because each algorithm has the potential of delivering different results, as demonstrated by several studies [20,53,57,58]. Since there is no universally and univocally accepted method to rank candidate reference genes on the basis of their expression stability, the present study analysed a set of wheat samples using GeNorm [13] and NormFinder [15] and compared the results obtained by the two programs. This comparison showed some discrepancies in the ranking of the candidate reference genes and in the identification of the best ones obtained by the two programs, but there was substantial agreement if the grouping of the genes with the most and least stable expression was considered ( Table 5).
As discussed by different authors [15,[57][58][59], since geNorm determines the reference gene stability by pairwise comparison of the variation of expression ratios, it needs a careful choice of the candidate genes included in the analysis. In particular, the candidate genes should not be co-regulated to avoid the selection of genes with the highest similarity of expression pattern but not stably expressed, consequently not suitable as reference genes.
To reduce the chance of using co-regulated genes, in our study special attention was paid in order to select for qRT-PCR analysis genes involved in distinct biological processes and metabolic pathways, even though for some of the new reference genes selected (e.g. Ta54733, Ta54448, Ta53889 and Ta53891 in Table 2) the functional information associated with their annotation was inadequate. Moreover, it should be remarked that besides similarities of gene expression pathways other factors can contribute to co-regulation. Yu et al. [60], for instance, showed that genes targeted by similar transcription factors have complex relationships across the co-regulated genes. Since some studies have shown that the removal of genes from the analysis modified the ranking of the candidate genes [58,61], we checked for possible co-regulation between the most stably expressed genes by removing single high ranking genes (as determined by pair wise comparisons), and reassessing the rank ordering of the remaining reference genes. Except for the floral organs (IV Data set), changes in rank ordering of the ten most stable genes were observed in the remaining four data sets when single high ranking genes were removed from geNorm analysis (data not shown).
The mathematical model of gene expression used by NormFinder enables to estimate both the overall variation of the candidate normalization genes and the variation between subgroups of the sample set. This approach has been reported to perform in a more robust manner than geNorm and has been shown to be less sensitive to the presence of co-regulated genes [15], therefore the choice of candidate reference genes without functional relationships is less critical than in geNorm. A very useful information given by geNorm is the ideal number of reference genes which should be included in a NF (Normalization Factor), whereas NormFinder indicates only the best single gene and the best combination of two genes. As shown by several studies [15,62,63], the normalization based on multiple reference genes gives more accurate evaluation of gene expression, especially when no single optimal reference gene is available. The cost/benefit evaluation is an important issue in the decision on the number of reference genes to be used for normalization, in particular Relative quantification of TaPDIL1-1 expression using selected reference genes as internal controls in 18 tissues and develop-mental stages of wheat Relative quantity Genorm NormFind Ta54227  Ta54825  Ta25534 Relative quantification of TaPDIL1-1 transcripts in three data sets using selected reference genes as internal controls Figure 6 Relative quantification of TaPDIL1-1 transcripts in three data sets using selected reference genes as internal controls. A) Six vegetative tissues and shoot developmental stages (shoots, stems and leaves); B) six floral organs from fully emerged spikes; C) six samples consisting of two temperature treatments (4°C and 33°C) for 24 and 48 h and their controls (LT = low temperature; HT = high temperature). The TaPDIL1-1 expression levels were normalized using the five methods described in Fig. 5. Normalized values of TaPDIL1-1 relative expression are given as average ± SD.
when large experiments, based on the analysis of a number of tissues, are planned. In such cases a compromise needs to be achieved between the instances of economical adequacy and the required degree of accuracy. In all the four data sets analysed to study the expression of TaPDIL-1-1, we observed that the normalization obtained using the single most stable reference gene identified by NormFinder was comparable to that based on the NF computed using the geometric mean of the two or three most stable genes selected by NormFinder and geNorm, respectively (Figs. 5 and 6). It might be concluded that in our experimental conditions a single reference gene was sufficient to obtain an adequate normalization of the qRT-PCR data and that NormFinder was more effective than geNorm in selecting the most effective reference genes.
Expression stability analyses showed that the most stable reference gene identified by NormFinder and geNorm was not the same for the five data sets analysed (Table 5). Even though our data confirm, as shown by most papers on the subject, that it is not possible to find an ideal all purpose reference gene, they highlight the convenience of pilot studies to select the most suitable gene for specific experimental conditions, as shown also by Gutierrez et al. [23]. When the complete data set including 24 samples was analysed, NormFinder identified Ta54227 (Cell division control prot., AAA-superfamily of ATPases), Ta2291 (ADP-ribosylation factor) and Ta2776 (RNase L inhibitorlike protein) as the most stably expressed genes (Table 4), whereas in geNorm the same genes ranked among the seven most stable genes (Table 5). Ta54227 can be identified as the best reference gene to normalize expression data from different wheat tissues and developmental stages, in fact NormFinder included Ta54227 among the two most stable genes in three of the analysed data sets (II, IV and V data sets, Table 4). Also geNorm included Ta54227 among the four most stable genes in three of the analysed data sets (Figs A, C and D in Additional file 11). However, Ta54227 is not suitable as control gene for analysing the samples subjected to cold and heat shocks, because there was considerable variation of its expression in this data set, wherein Ta2291 and Ta2776 performed much better (Table 5). Moreover, since the expression level of Ta2776 is low (mean C t 24.79), whereas that of Ta2291 (mean C t 21.97) and of Ta54227 (mean C t 21.56) is high, they are suitable as reference genes in a large data set of wheat samples, but also for the normalization of target genes with contrasting expression levels.
It is noteworthy that most of the novel reference genes identified in this study were more stably expressed than well-known and frequently used HKGs, such as Actin (Ta54825), α-tubulin (Ta25534), β-tubulin (Ta44405), Ubiquitin (Ta50503), GAPDH (Ta30768), Cyclophilin (Ta35497), Translation elongation factor (Ta53964), Translation initiation factor (Ta54280), Ribosomal protein (Ta27771) and Histone (Ta38797), in fact most of them classified as the least stable genes both by geNorm and NormFinder in the five data sets analysed (Table 5). When the complete data set including 24 samples was analysed by NormFinder, Ta22845 (ATP-dependent 26S proteasome regulatory subunit), which resulted the gene with the highest expression stability among the selected traditional HKGs, ranked only 12 th among the 32 candidate reference genes (I Data set in Table 4). In the 18 samples including different tissues/organs at different developmental stages of wheat the best HKG identified by NormFinder was Ta35497 (Cyclophilin) which ranked only tenth (II Data set in Table 4), whereas in the temperature treatment series the most stable HKG was again Ta22845 (position 8 th in III Data set of Table 4). Finally, Ta44405 (β-tubulin) was identified as the most stable HKG among the single floral organs (position 10 in the IV Data set of Table 5), whereas only in the six samples including vegetative tissues Ta50503 (ubiquitin), a traditional HKG, was classified as one of the most stable genes both by NormFinder and geNorm (V data set in Table 5).
The example of expression level analysis of TaPDIL1-1 shows that the selection of stable reference genes represents a crucial issue for the correct normalization of qRT-PCR data. The comparison of the relative TaPDIL1-1 mRNA expression using different normalization approaches (Figs. 5 and 6) showed that Actin (Ta54825) and α-tubulin (Ta25534), most frequently used in wheat gene profiling studies (Additional file 1), and most traditional HKGs (data not shown) analysed in our research were less effective than the best reference genes identified in this study.

Conclusion
In the present research we compared 20 novel candidate reference genes and 12 commonly used HKGs suitable for gene expression normalization in wheat. Both new and traditional reference gene sequences were fetched by an in silico cross search for stable expression in Unigene and TIGR databases. The expression stability of the 32 selected genes was assessed by qRT-PCR in 24 different plant samples, which included different tissues, developmental stages and plants exposed to temperature stresses. The computer programs geNorm and NormFinder were used to compare the expression patterns of the 32 candidate genes and to identify the best reference genes, which were not the same. However, when single genes were not considered, compatible results and substantial agreement was observed between the two programs in discriminating between groups with the most and least stable expression.
Our results indicate that many of the new identified reference genes outperform the traditional HKGs in terms of expression stability under all the tested conditions. Ta54227 (Cell division control prot., AAA-superfamily of ATPases), Ta2291 (ADP-ribosylation factor) and Ta2776 (RNase L inhibitor-like protein) resulted the best and most stable reference genes to normalize gene expression in different tissue and development stages of wheat; they outperformed all traditional HKGs such as α-tubulin (Ta25534), β-tubulin (Ta44405), Ubiquitin (Ta50503), Actin (Ta54825), GAPDH (Ta30768), Ribosomal protein (Ta27771) and Histone (Ta38797). Although in our experimental conditions a single reference gene was adequate to normalize the qRT-PCR data, the use of normalization factors generated from the expression data of two or three most stable genes should be considered in future studies to improve the accuracy and reliability of expression analysis. The set of qRT-PCR specific primers for the novel and more stable reference genes developed in this study will enable more accurate normalization and quantification of gene expression in wheat. Moreover, the design of suitable primer pairs for their orthologous genes could contribute to the development of reference genes with stable expression in other plant species.