Comprehensive Selection of Reference Genes for Gene Expression Normalization in Sugarcane by Real Time Quantitative RT-PCR

The increasingly used real time quantitative reverse transcription-PCR (qRT-PCR) method for gene expression analysis requires one or several reference gene(s) acting as normalization factor(s). In order to facilitate gene expression studies in sugarcane (Saccharum officinarum), a non-model plant with limited genome information, the stability of 13 candidate reference genes was evaluated. The geNorm, NormFinder and deltaCt methods were used for selecting stably expressed internal controls across different tissues and under various experimental treatments. These results revealed that, among these 13 candidate reference genes, GAPDH, eEF-1a and eIF-4α were the most stable and suitable for use as normalization factors across all various experimental samples. In addition, APRT could be a candidate for examining the relationship between gene copy number and transcript levels in sugarcane tissue samples. According to the results evaluated by geNorm, combining CUL and eEF-1α in hormone treatment experiments; CAC and CUL in abiotic stress tests; GAPDH, eEF-1a and CUL in all treatment samples plus CAC, CUL, APRT and TIPS-41 in cultivar tissues as groups for normalization would lead to more accurate and reliable expression quantification in sugarcane. This is the first systematic validation of reference genes for quantification of transcript expression profiles in sugarcane. This study should provide useful information for selecting reference genes for more accurate quantification of gene expression in sugarcane and other plant species.


Introduction
Real time quantitative reverse transcription-PCR (qRT-PCR) is increasingly used in gene expression analysis owing to its simple, reproducible and high-throughput features. qRT-PCR provides a useful and rapid means of understanding gene expression in living organisms by measuring the expression of target genes across different samples. In addition, qRT-PCR is a low cost and widely accepted method in the tracking of gene expression levels in genetically modified organisms (GMO) as well as in molecular breeding and gene mining. When performing qRT-PCR analysis, several factors such as sample amount, RNA integrity, cDNA quality, as well as the tissues or cell activities, can affect the quantitative measurement of gene expression [1]. Thus, in order to obtain a reliable analysis of gene expression by qRT-PCR, one or several reference genes should serve as the internal control to normalize and monitor the expression variation between samples and reactions [2][3][4]. The expression of these reference genes should remain stable under various experimental treatments and/ or at different stages of development and growth periods [2][3][4]. Specifically, a suitable reference gene for performing qRT-PCR analysis should: (i) have stable expression across all or most of the samples analyzed; (ii) have no association with any pseudogene, to avoid the amplification of non-functional gene family members; (iii) reflect variations in RNA quality and quantity, as well as cDNA synthesis and amplification; (iv) possess the stability of transcription that is suitable for the target gene; and (v) exhibit moderate expression levels (i.e. a threshold cycle Ct of 15 to 30) [5].
Housekeeping genes (HKGs) related to basal cell activities and cellular structure components have been historically used as internal controls in medical science [6] and later in plant science [7]. The most commonly used HKGs are 25S rRNA, GAPDH (glyceraldehyde-3-phosphate dehydrogenase), ACT (b or c actin) and TUB (a or b tubulin). However, using a genome wide approach, hundreds of genes in Arabidopsis thaliana were shown to outperform traditional reference genes in terms of expression stability throughout development and under a range of environmental conditions [8]. Among them, genes encoding a protein phosphatase 2A subunit, a coatomer subunit and an ubiquitinconjugating enzyme were identified as novel reference genes [8].
Sugarcane (S. officinarum6S. spontaneum) is a widely grown sugar crop in the tropics and subtropics with increasing demand due to biofuel production and challenges with biomass production [31]. Identifying genes for sucrose accumulation and stress resistance, which requires a set of reference genes for gene expression normalization, could serve to increase sugarcane yield and sucrose content using both genetically modified strategies and molecular marker-assisted breeding. Sugarcane has limited available genomic information so only a few of genes, mainly GAPDH and 25S rRNA, have been verified as practicable reference genes [32,33]. However, there is expression variation of GAPDH and 25S rRNA under various experimental conditions [2,12,34].
In the present study, the 13 genes 25S rRNA, GAPDH, ACT (bactin), TUB (b-tubulin), APRT, PRR, 18S rRNA, eEF-1a, eIF-4a, CAC, TIPS-41, CUL and UBQ were selected as candidate reference genes for evaluation in sugarcane. This study utilized 57 sugarcane samples and aimed to reveal which reference genes should be used in experimental samples with different treatments or different tissues. A combination of reference genes was also introduced to evaluate their potential for more accurate and reliable qRT-PCR analysis of gene expression in sugarcane.

Ethics Statement
No specific permissions were required for these locations/ activities. The field studies did not involve endangered or protected species and the specific location of this study is longitude: 119.23E, latitude: 26.08N.

Plant Materials, Growth Conditions and Treatments
Five sugarcane cultivars, ''ROC''20, ''ROC''22, FN40, Liu-cheng03-182 and YC05-179, were used for tissue sample collection. Five tissue samples were harvested (leaf, leaf sheath, stem epidermis, stem pith and bud) from 7-to 8-month-old sugarcane grown in the field. For each cultivar, 6 plants from the same experimental plot were collected to provide 6 replicates. The leaf and leaf sheath samples were taken from the last fully expanded leaf (+1 leaf), while the stem samples (epidermis and pith) and the buds were harvested from the 6 th or 7 th stem internodes. All materials were cut into small pieces, wrapped in tinfoil, immediately snap-frozen in liquid nitrogen and then kept at 280uC until RNA isolation.
Four sugarcane cultivars, ''ROC''20, FN40, Liucheng03-182 and YC05-179, were used for the stress experiments. The single bud node shoots that used for in vitro disease-free plantlets regeneration, were incubated in 50uC water with the fungicide carbendazim (100 mg?L 21 ; Tianjin, China) for 40 min. The shoots were then planted in autoclaved soil before harvest for meristem excision, callus induction, shoot differentiation and rooting [35]. The plantlets were then placed in distilled water and kept for ten days in a tissue culture room under a constant temperature of 2561uC. Different sets of plants were then transferred into test tubes containing water solutions along with different treatments, including abscisic acid (100 mM, ABA), methyl jasmonate (25 mM, MeJA) and salicylic acid (5 mM, SA) for 6 h; hydrogen peroxide (500 mM, H 2 O 2 ), sodium chloride (250 mM, NaCl) and polyethylene glycol (25% w/v, PEG) for 12 h or copper chloride (100 mM, CuCl 2 ) and cadmium chloride (500 mM, CdCl 2 ) for 24 h. Each set of samples comprised three seedlings as biological replicates for expression analyses. The plants without treatment (kept in distilled water) were harvested as control. All materials were wrapped up in tinfoil, immediately snap-frozen in liquid nitrogen, and then kept at 280uC until RNA isolation.
A total of 57 samples (56 plus the control), including 25 tissue samples and 31 treated samples (the H 2 O 2 treated sample of YC05-179 was absent) exposed to various stress treatments involving four different cultivars, were employed in the experiments for evaluation of candidate reference genes in sugarcane. The first-strand cDNA was synthesized with a 10 mL reaction system according to the instructions of the TAKARA PrimeScrit RT Reagent Kit (Perfect for Real Time) (TaKaRa Biotechnology, Dalian, China). The quality and integrity of cDNA were determined in the same way as above. All the cDNA samples were diluted to 5 ng?mL 21 for the following qRT-PCR reaction and stored at 220uC until use.

qRT-PCR and Data Analyses
An ABI 7500 Real Time PCR machine (Applied Biosystems, Foster City, CA, USA) and its default program (2 min at 50uC and 10 min at 95uC followed by 40 cycles at 94uC for 15 s, and at 60uC for 60 s.) were employed for qRT-PCR with a reaction mixture volume of 20 mL in an optical 96-well plate. 10 mL of SYBR Green Master Mix (Roche), 10 pM of each primer, 10 ng of final cDNA and 6.4 mL of RNase-free water were added to the reaction mixture. A control was also included in each plate with 2.0 mL of RNase-free water as a template. Three technical replicates were contained in each plate. Specificity verification of the PCR amplification dissociation and the PCR efficiency curves were determined for each candidate reference gene prior to the qRT-PCR evaluation of these genes in sugarcane.

Statistical Analysis of Gene Expression and Comparison of Normalization Methods
The PCR amplification efficiencies of all gene assays were estimated based on the slope of the standard curve in conjunction with the following formula: Standard, PCR efficiency and the correlation coefficient (R 2 ) curves of each gene were generated in Microsoft Excel 2003 using a range of seven dilutions made in ten, five or three-fold decrements on the YC05-179 control sample. GAPDH and TUB were ten-fold; 25S rRNA, ACT, 18S rRNA, eEF-1a, CUL and eIF-4a were five-fold and CAC, TIPS-41, UBQ, APRT and PRR were three-fold.
The overall expression levels of the candidate genes were transformed into threshold cycle (Ct) values by the ABI 7500 Fast Real Time PCR System. Ct values over 40 indicated undetectable product and were considered as missing values (NA) for subsequent calculations. After collecting and converting the threshold PCR cycle data, Ct average values for the evaluated genes were inputted into software according to the corresponding manuals of geNorm (trial version; Biogazelle, Zwijnaarde, Belgium) [2] and NormFinder (ver. 0.953) [1]. After that, the values of stability of candidate genes achieved from geNorm, NormFinder and deltaCt (only the deltaCt results from RefFinder were employed [16]) were used to calculate Pearson correlation values (r value) by SAS S21.0. The values reflect the level of correlation of between the results from geNorm, NormFinder and deltaCt.

Screening of Candidate Reference Genes and Primer Design
A total of 13 candidate reference genes were selected in sugarcane or in other plant species for evaluation on the basis of their stable expression across developmental stages and/or abiotic stresses. These included four previously assessed sugarcane candidate genes, 25S rRNA, GAPDH, ACT and TUB [31,32], as well as nine new candidate reference genes, APRT, PRR, 18S rRNA, eEF-1a, eIF-4a, UBQ, CAC, TIPS-41 and CUL. Since there is limited available sugarcane genome sequence information, the publicly available gene sequences from O. sativa (18S rRNA, AK059783; eEF-1a, AK061464; eIF-4a, AK073620; UBQ5, AK061988), Z. mays (CUL, GRMZM2G166694_T04) and A. thaliana (CAC, TIPS-41) were used as the probes to search within a sugarcane expressed sequence tags (ESTs) database (www.ncbi. nlm.nih.gov/nucest/?term = sugarcane). Two candidate reference genes in sugarcane which performed well in B. juncea [12], named CAC and TIPS-41, were identified by querying homologous sugarcane sequences with A. thaliana genes complete CDS (CAC, At5G46630; TIPS-41, At4G34270). The remaining five new candidates, 18S rRNA, eEF-1a, eIF-4a, UBQ, and CUL, performed well in O. sativa [9] and Z. mays [11]. All seven EST sequences, including those of 18S rRNA, eEF-1a, eIF-4a, UBQ5, CAC, TIPS-41 and CUL, were acquired from the publicly available database in NCBI using candidate ESTs with the highest homology to the target sequences. Both reference sequences from O. sativa, Z. mays, B. juncea, T. aestivum or A. thaliana and the corresponding target sequences from sugarcane were aligned together in DNAMAN to identify the complete sequence identity. This information was used to design primers using the Primer-BLAST tool from the NCBI (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primer sequences for all 13 candidate reference genes are shown in Table 1.

Verification of Primer Specificity, Efficiency and Gene Expression Profile
Using the control sample of YC05-179 as the cDNA template, the specificity of primers used in qRT-PCR experiments was confirmed based on a melting curve analysis and agarose gel electrophoresis. The standard curve, PCR efficiency and the correlation coefficient (R 2 ) of each gene were generated in Microsoft Excel 2003 using a range of seven dilutions made in ten (GAPDH and TUB), five (25S rRNA, 18S rRNA, eEF-1a, eIF-4a and CUL) or three-fold (CAC, TIPS-41, UBQ, APRT and PRR) decrements on the YC05-179 control sample. The qRT-PCR efficiency formula (Eq. 1) was used in the calculation. The regression coefficient R 2 for all the primers varied between 0.9876-0.9999 over a serial dilution of cDNA. qRT-PCR efficiencies of primers ranged from 93.24% to 113.83% (Table 1).
Over all samples, the Ct values of the 13 candidate reference genes varied over a wide range, and the mean Ct values of these genes across all the samples ranged from 14.27 to 28.21 (Table 1). Among these candidate reference genes, 25S rRNA was the most abundantly expressed gene in all of the samples (mean Ct6SD = 14.2760.71) followed by 18S rRNA (mean Ct6SD = 15.3860.74), whereas PRR was the least abundantly expressed gene (mean Ct = 28.2162.04). ACT, GAPDH and eEF-1a were close in Ct values (25.03, 24.64 and 24.24, respectively) but differed in SD values (1.06, 1.23 and 1.33, respectively). The mean Ct6SD of the remaining genes (TUB, TIPS-41, CAC, CUL, eIF-4a, UBQ, APRT and PRR) varied from 26.4761.29 to 28.2162.04. The data also showed that the expression of ACT and eIF-4a were the least variable (Covariance (CV) of 4.23% and 4.42%, respectively). The CV values of 18S rRNA, UBQ, 25S rRNA, GAPDH, CAC and APRT ranged from 4.81% to 5.05%. The CV values of TUB, eEF-1a and CUL ranged from 5.31% to 5.52%. In contrast, the CV values of TIPS-41 and PRR were 7.67% and 7.25%, respectively (Table 1).
In order to make a wider expression analysis of candidate reference genes, the 57 diverse samples were divided into four experimental sets. The 1 st set consisted of leaf, leaf sheath, stem epidermis, stem pith and bud samples from ''ROC''20, ''ROC''22, FN40, Liucheng03-182 and YC05-179 ( Table 2). The 2 nd set was comprised of 12 samples from FN40, Liucheng03-182, ''ROC''20 and YC05-179 treated with ABA, MeJA and SA. The 3 rd set contained 19 samples treated with H 2 O 2 , NaCl, PEG, CuCl 2 and CdCl 2 (abiotic stresses). In the 4 th set, both the 2 nd and 3 rd sample sets were included. The variation of transcript levels indicated that the expression of the candidate genes was affected by tissue types and experimental conditions ( Table 2). In the 1 st set, 18S rRNA performed as the least variable gene, whereas 25S rRNA, CUL and PRR had a larger range of expression level variation. 25S rRNA and eIF-4a had the lowest expression variation in the 2 nd set and, 25S rRNA and ACT had the lowest variation in the 4 th set. ACT and TUB performed less variability in the 3 rd set. Additionally, TIPS-41&TUB and TIPS-41&PRR displayed the most variable expression profiles when treated with hormones and abiotic stresses, respectively (Table 2). A significant variation in the expression of PRR and APRT was observed in all treatment samples ( reference genes displayed constant expression levels throughout the various cultivars and among the different kinds of treatments, six genes among these 13 reference genes (25S rRNA, 18S rRNA, GAPDH, ACT, UBQ and eIF-4a) varied at a relatively lower level according to CV values (Table 2).

Expression Stability Analysis and Ranking of All 13 Candidate Reference Genes in Sugarcane
Based on the expression stability analysis of the 13 candidate reference genes by geNorm, the six top ranked genes in different cultivar tissues were CUL = CAC.APRT.TIPS-41.GAPDH. eEF-1a. The top ranked genes in the cultured plantlets under hormone treatments were CUL = eEF-1a.GAPDH.CAC, and the rank of the top six genes in the 3d set was CAC = CUL. APRT. GAPDH.eEF-1a (Fig. 1). The top six genes that were expressed relatively constantly in the 4 th set were eEF-1a = GAPDH.CUL. CAC.APRT.eIF-4a (Fig. 1). It is interesting to note that the most stable traditional housekeeping gene was GAPDH compared to 25S rRNA, ACT and TUB, which is similar to previous reports [32,33]. Three reference genes, the eEF-1a, CAC and CUL, were more stable than the other six new candidates (18S rRNA, eIF-4a, UBQ, TIPS-41, APRT and PRR). Also, in response to various abiotic stresses and hormone treatments, UBQ or TIPS-41 were identified as the least stable genes (Fig. 1).
NormFinder analysis revealed that eEF-1a and eIF-4a were always two of the top three most stable reference genes in the treated samples, while the third place was occupied by ACT, GAPDH or TUB in the 2 nd , 3 rd and 4 th sets, respectively (Table 3). Three reference genes, GAPDH, eEF-1a and 18S rRNA, were most stable in the different tissues from all the five tested cultivars (Table 3). UBQ had the worst stability of the 13 candidate reference genes in untreated samples and in samples under hormone treatments (Table 3). Likewise, TIPS-41 ranked at the bottom of these candidate genes in both the 3 rd and 4 th sets (Table 3). Though ACT, TUB and GAPDH expressed at less variable levels in the 1 st , 2 nd and 4 th sets, the expression of all the four traditional housekeeping genes fluctuated across all four experimental sets and was less stable than eIF-4a or eEF-1a in most of the samples ( Table 3).
The deltaCt method was used in RefFinder, which integrates geNorm, NormFinder, BestKeeper and deltaCt into a web-based program. The ranking of the 13 genes based on the deltaCt was consistent with that of NormFinder in different tissue samples except for the ranking of 25S rRNA (Table 3 and Table 4). Comparing the average expression stability of the top six genes valued by geNorm, NormFinder and deltaCt, the three candidate genes GAPDH, eEF-1a and eIF-4a performed better in all the 31 treated samples, and eEF-1a, APRT and GAPDH performed better in all the 25 tissue samples (Table 4). Conversely, using the three algorithms, TIPS-41 and UBQ were the worst performing genes of all treated samples and untreated samples, respectively (Table 4). Pearson correlations were calculated among the geNorm, NormFinder and deltaCt methods using the stability values in the 1 st and 4 th sets. Interestingly, the Pearson correlations among all three stability tests were positive. A significant correlation was observed between outcomes from NormFinder and deltaCt (r = 0.946), indicating that the ranking results of all 13 reference genes from the above two methods were nearly identical in the 1 st set ( Table 5). The lower correlations (Table 5) between geNorm and NormFinder in the 4 th set (r = 0.438) and the 1 st set (r = 0.476) were reflected by the results of the different rankings (Table 4).

Optimal Number of Reference Genes for Expression Normalization across Different Experimental Sets
The optimal number of reference genes in the normalizing experiment was determined by geNorm by calculating the pairwise variation (V n /V n+1 ) between the normalization factors (NF) across all the samples of the different experimental sets. The pairwise variation (V = V n /V n+1 ), which was counted between NF n and NF n+1 , was used to find the best combination of genes for reliable normalization [2]. Vandesompele et al. suggested a threshold of V = 0.15, which indicates that adding one more gene has little influence on the calculation of normalization factor [2]. As shown in Figure 2, the rank order of gene stability established by means of stepwise exclusion of the least stable gene suggested that the combination of CAC, CUL, APRT and TIPS-41 could provide a dependable result while normalizing the qRT-PCR data of the target gene in the 1 st sample set (Fig. 2). However, in the 2 nd set the two reference genes (CUL and eEF-1a) were enough to achieve a V 2/3 = 0.134, which is close to the V = 0.15 threshold proposed by Vandeompele et al [2] and is thus an efficient and economical strategy to quantify sugarcane samples from hormone treatments. Similarly, the use of CAC and CUL was enough to achieve a V 2/3 = 0.144 in samples treated with compounds eliciting abiotic stress (Fig. 2). The V 3/4 value in the 4 th set samples from four sugarcane cultivars was 0.101 which suggests that the combination of GAPDH, eEF-1a and CUL is the best choice for quantification of gene expression in qRT-PCR (Fig. 2).

Discussion
Sugarcane, sorghum, maize and rice belong to the Andropogoneae tribe and share many similarities in their genetic composition [36][37][38][39][40]. Rice, O. sativa, which diverged from a common Andropogoneae ancestor around 50-70 million years ago, performs as a model plant for grass species in several fields of modern molecular biology research [41]. Modern sugarcane cultivars are hybridized crosses between S. officinarum and S. spontaneum, resulting in a high degree of polyploidy and frequent aneuploidy [42]. Comparisons involving sugarcane as well as O. sativa, S. bicolor and Z. mays indicated significant common conservation of gene content and a few rearrangements [40]. The results of the present study also showed that the stable expression of candidate genes, such as eIF-4a, may follow the same pattern in S. officinarum, O. sativa, S. bicolor and Z. mays.
qRT-PCR is a quick, reliable and accurate tool to analyze gene expression. Accurate mRNA normalization requires one or two internal control genes which are stably expressed should be used.
Since no single control is appropriate for all experimental treatments, it is generally suggested to select suitable internal controls prior to use for normalization. Although it is generally considered that housekeeping genes are expressed constantly, Lilly et al. and Nicot et al. revealed that the expression of such genes can undergo significant stability changes during biotic or abiotic stresses [18,43]. Thus, several statistical algorithms, such as geNorm [2], NormFinder [1], BestKeeper [17], the deltaCt method [10] and the RefFinder WEB-based software [11], were developed to assess expression stability. Among these algorithms, geNorm, NormFinder and the deltaCt method were used in the Table 3. Expression stability of 13 reference genes in Saccharum officinarum as calculated by Normfinder.  present work. The deltaCt results were obtained from RefFinder, a web-based tool that integrates the most commonly used statistical algorithms. The normalization results were mostly affected by variation of the quantity and quality of introducing RNA and cDNA or even reaction-to-reaction variations [1]. When employing a variable reference gene, Gutierrez et al. revealed that nearly 100-fold variations could be found when quantifying target gene expression, thus leading to a huge potential scope for misinterpretation of the expression pattern of target genes [7]. Guo et al. observed much different expression values during a study of sugarcane dirigent protein gene expressed in sugarcane stems when using the different internal control genes of 25S rRNA and GAPDH [44]. Therefore, to achieve a reliable result with qRT-PCR, a systematically validated reference gene should be used and taken as an essential component of qRT-PCR analysis. In the present study, the three statistical algorithms, geNorm, Normfinder and deltaCt, were used and the Ct values were inputted into RefFinder. When correlation analysis based on the ranking order of the 13 candidate genes was applied for comparison of the three statistical algorithms, the results showed that the correlation coefficient between NormFinder and deltaCt was more positive (significantly in the cultivar set, r = 0.946) than the correlations between geNorm and Normfinder or between geNorm and deltaCt. These results were in accordance with the results of Jacob et al. [45]. The ranking order of the evaluated candidate genes by NormFinder and deltaCt in our study was generally consistent.
According to the CV values, 25S rRNA, 18S rRNA, GAPDH, ACT, CAC, eIF-4a, UBQ and APRT showed low variation in gene expression in all four experimental sets (Table 2). However, only GAPDH, eEF-1a and eIF-4a performed well in all three statistical algorithms. These three genes had relatively stable expression in different sugarcane genotypes and tissues. Our results also showed that the three genes CUL, eEF-1a and eIF-4a were suited for gene expression normalization in hormone-treated experiments, while GAPDH and eEF-1a were ideal when analyzing samples comparable to the 3 rd sample set. However, when these two sets (the 2 nd and 3 rd sets) were integrated into the 4 th set and analyzed together, eEF-1a, eIF-4a and GAPDH were the most stable reference genes.
In previous studies, GAPDH was found to possess the most stable mRNA expression in sugarcane [32] and 25S rRNA was the most stable gene in sugarcane infected with Ustilago scitaminea [33]. Similarly, this study identified GAPDH as one of the most stable reference genes for all types of hormone and abiotic treatments. 25S rRNA and 18S rRNA were two of the more abundantly expressed target genes under investigation, which is contrary to the reference gene selection principle that emphasizes moderate expression levels. It is difficult to detect the variation of 25S rRNA and 18S rRNA due to their rich content, so these two genes may not be ideal for expression normalization. It should be pointed out that some genes with stable but low expression levels, such as genes  coding for certain transcription factors may be used as internal references in qRT-PCR experiments. A suitable reference should have a minimum difference in Ct value since this would lessen influence on the quantification and hence be more accurate [8,54].
In the present study, four historically used housekeeping genes, GAPDH, 25S rRNA, ACT and TUB, were also included in the evaluation [32,33]. 25S rRNA, ACT and TUB had poor performance across all four experimental sets of sugarcane samples, and similar poor performance was seen when using some of the novel candidates, such as UBQ, TIPS-41, PRR and 18S rRNA.
In this study, eEF-1a, which has been proven to be a suitable reference gene for expression normalization in O. sativa and C. sativus [5,12], also ranked at the top when evaluated by geNorm, NormFinder and deltaCt across our four experimental sets. Being an important specific protein factor involved in the process of protein translation, eIF-4a had the same performance as eEF-1a in most experimental treatment conditions in our study. Our study showed that eIF-4a had a relatively high expression under most experimental conditions, which is in agreement with the study of Zhu et al. in C. papaya [46]. eIF-4a also performed well in Musa paradisiaca, Lycoris longituba, Hevea brasiliensis and Coffea spp. [47][48][49][50]. The low copy number of APRT in the sugarcane genome, which was reported by Casu et al. [30], provides an advantage in analyzing gene copy number. However, in the present study the expression of APRT was found to be easily affected by abiotic stress conditions (3 rd set) and thus its application is limited. With more copies in the sugarcane genome than APRT [30], PRR displayed variable performance under different stresses and in different sugarcane tissues.

Conclusion
In the present study, 13 candidate reference genes were evaluated by three independent analysis approaches according to stability of transcript profiles across sugarcane samples. The aim was to select the most suitable reference genes for further gene expression quantification in different tissues and stress-treated samples in sugarcane. The stability analysis of gene expression by geNorm, NormFinder and deltaCt revealed that GAPDH, eEF-1a and eIF-4a are the most suitable for normalization controls across different samples. This is the first systematic validation of reference genes for quantification of transcript profiles in sugarcane. In our study, combining different reference genes is advocated for reliable normalization in different experimental samples. In addition, we also provide the procedure for identifying suitable reference genes by qRT-PCR in detail, which will be helpful in other plant species. Therefore, this study should provide useful information for selecting reference genes for accurate quantification of gene expression in sugarcane and other plant species.