Evaluation of New Reference Genes in Papaya for Accurate Transcript Normalization under Different Experimental Conditions

Real-time reverse transcription PCR (RT-qPCR) is a preferred method for rapid and accurate quantification of gene expression studies. Appropriate application of RT-qPCR requires accurate normalization though the use of reference genes. As no single reference gene is universally suitable for all experiments, thus reference gene(s) validation under different experimental conditions is crucial for RT-qPCR analysis. To date, only a few studies on reference genes have been done in other plants but none in papaya. In the present work, we selected 21 candidate reference genes, and evaluated their expression stability in 246 papaya fruit samples using three algorithms, geNorm, NormFinder and RefFinder. The samples consisted of 13 sets collected under different experimental conditions, including various tissues, different storage temperatures, different cultivars, developmental stages, postharvest ripening, modified atmosphere packaging, 1-methylcyclopropene (1-MCP) treatment, hot water treatment, biotic stress and hormone treatment. Our results demonstrated that expression stability varied greatly between reference genes and that different suitable reference gene(s) or combination of reference genes for normalization should be validated according to the experimental conditions. In general, the internal reference genes EIF (Eukaryotic initiation factor 4A), TBP1 (TATA binding protein 1) and TBP2 (TATA binding protein 2) genes had a good performance under most experimental conditions, whereas the most widely present used reference genes, ACTIN (Actin 2), 18S rRNA (18S ribosomal RNA) and GAPDH (Glyceraldehyde-3-phosphate dehydrogenase) were not suitable in many experimental conditions. In addition, two commonly used programs, geNorm and Normfinder, were proved sufficient for the validation. This work provides the first systematic analysis for the selection of superior reference genes for accurate transcript normalization in papaya under different experimental conditions.


Introduction
Gene expression analysis is an important step to understand the roles of genes in developmental and cellular processes, such as the signaling and metabolic pathways [1]. Real-time reverse transcription PCR (RT-qPCR) has emerged as the most widely used method to quantify changes in gene expression profiles in response to different environmental conditions. It has shown important attributes such as accuracy, precision and relative ease of use due to its speed, sensitivity and specificity [2,3]. Nevertheless, to accurately quantify gene expression, several experimental variations, such as quality and amount of starting material, presence of inhibitors in different sample materials, primer design, and RNA extraction and retro-transcription efficiencies, should be taken into account [4]. Therefore, selection of an appropriate normalization strategy is of crucial importance for the acquisition of biologically meaningful data. Among several methods proposed so far [4,5], the use of one or more reference genes is currently the preferred method of normalization [6]. An ideal reference gene should be expressed at a constant level across various conditions and its expression is assumed to be unaffected by experimental parameters [7,8]. Moreover, the reference gene and the target gene should have similar ranges of expression in the samples to be analyzed [9]. Genes involved in basic metabolism and maintenance of the cell, e.g. b-actin, glyceraldehyde-3-phosphatede hydrogenase (GAPDH), ribosomal subunits and ubiquitin are commonly used as reference genes [10,11,12]. However, several reports have demonstrated that there are no universally applicable reference genes with an invariant expression, and that the using of unstable reference gene will lead to inappropriate biological data interpretation [8,13,14]. Thus, there is an urgent need to systematically evaluate the stability of potential reference genes for every particular experimental condition prior to their use in RT-qPCR normalization. Meanwhile, several algorithms, such as geNorm [15], NormFinder [16], BestKeeper [17], qBasePlus [18], and RefFinder [19] have been well developed to validate the most stable reference gene(s) from a panel of potential genes or candidate genes under a given set of experimental conditions [20,21].
Papaya (Carica papaya L.) is the only species within the genus Carica and the most commercially important species within the family Caricaceae [46] and it has been widely cultivated in tropical and subtropical lowland regions for its nutritional benefits and medicinal applications. Papaya is also the first perennial transgenic fruit variety for commercial application of the world [47,48]. However, papaya fruit is subject to some problems such as rapid ripening, susceptible to biotic or abiotic stresses, which could result in a high percentage of product loss [46,47,49]. Due to all these matters, papaya has been the focus of many studies at physicochemical, biochemical, and molecular levels [46]. The postharvest biology of papaya fruit has been an important aspect of those studies as well [46,50,51,52]. The understanding of expression patterns of some key genes, especially for the genes associated with ripening and stress responses, will help us to gain insights into the mechanisms involved in these processes, and in turn, to improve fruit quality and storage potential. To date, almost all studies on gene expression in papaya fruit with RT-qPCR have used Actin or 18S rRNA as reference gene [10,11,53,54]. However, the stability of these two genes has not been verified yet and it is not clear whether they are the suitable reference genes in papaya. Therefore, the application of RT-qPCR analyses of gene expression in papaya fruit has been limited by the use of potentially inappropriate reference genes.
For further development of RT-qPCR in papaya, the present study aimed at defining reference genes suited for quantitative analysis of papaya genes under different experimental conditions. Here, we reported a systematic analysis of 21 genes to identify the internal reference gene(s) most suitable for normalization gene expression data obtained with RT-qPCR analysis in papaya. These genes have different roles in the cells, including those involved in cell structure, membrane proteins, transcription, protein translation, protein degradation and metabolic pathways. The data for each gene were obtained from a large set of biological samples representing different experimental conditions, including various tissues, fruit developmental stages, different storage temperatures, different cultivars, postharvest ripening, pathogen stress, 1-MCP treatment, hot water treatment, modified atmosphere packaging (MAP) and hormone treatment. Furthermore, in order to illustrate the usefulness of the newly identified reference genes, expression analysis of one interesting gene related to fruit softening, CpaEXY1, was presented. The result provided a superior set of validated reference genes that were suitable for RT-qPCR analysis in papaya fruit under different experimental conditions and clearly indicated that different reference genes should be validated according to the particular experimental conditions.

Plant Materials and Experimental Conditions
Roots, stems, leaves, flowers, and pre-climacteric papaya fruits at some hint yellow stage (three-line yellow) were freshly harvested from a local commercial plantation nearby Guangzhou, southeastern China. The vegetative tissue samples, such as root, leaves and stem, were taken from young tissue; flowers were harvested at full bloom. At each sampling time, plant materials except for fruits were immediately frozen in dry ice after harvesting, transported to the laboratory and then stored at 280uC until total RNA was isolated. For samples of different development stages, papaya fruits were sampled at 1, 2, 3, 4 and 5 months after anthesis.
For samples of different cultivars, different developmental stages, postharvest ripening, stresses, 1-methylcyclopropene (1-MCP) treatment, hormone treatment, MAP treatment and hot water treatments, pre-climacteric papaya fruit at the first sign of color break (,10% skin yellow stage) were harvested. Fruits free from visual symptoms of any disease or blemishes were randomly selected for uniformity of weight, shape, and maturity. The selected fruits were firstly cleaned, dipped in a 1% hypochloride solution for 1 min for contraction wounds and then soaked in 0.2% (w/v) Sporgon solution (Bayer, Leverkusen, Germany) for 10 min to eliminate potential microbes. They were then allowed to air-dry at 25uC for 3 h and treated as described below.
For biotic stress samples, the selected papaya fruits were inoculated with 20 ul (4610 6 spores ml 21 ) of Colletotrichum gloeosporioides Penz. spores in suspension as described by De Capdeville, et al. [52], covered with wet cotton and sealed with bag but open after 24 hours. Fruits were placed into unsealed plastic bags and stored at 25uC. Samples were taken at 0, 2, 4, 6, 8, 10 and 12 days after treatment.
For 1-MCP treatment, papaya fruits were sealed in the closed airtight containers, and 300 ml/L of 1-MCP were injected into the containers through a rubber septum. Fruits were incubated with 1-MCP for 16 h at 2561uC. The containers were then opened to allow ripening in air in the same temperature conditions and fruits were placed into unsealed plastic bags, as well. Samples were taken at 0, 2, 4, 6, 8, 10 and 12 days after treatment.
For hot-water treatment samples, papaya fruits were dipped into hot water (54uC) for 4 min and were taken out and then allowed to air-dry at 25uC. After that fruits were placed into unsealed plastic bags and stored at 25uC. Samples were taken at 0, 2, 4, 6, 8 and 10 days after treatment.
For samples of modified atmosphere packaging (MAP) treatment, the selected papaya fruits were packed with thick of 0.02 mm thick PE bags, sealed by capper and stored at 25uC. Samples were taken at 0, 3, 6, 9, 12, 15 and 17 days after treatment.
For samples of ethephon treatment, papaya fruits were dipped into ethephon aqueous solution (500 ml/L) for 3 min, and then taken out and allowed to air-dry at 25uC. After that the fruits were placed into unsealed plastic bags and stored at 25uC. Samples were taken at 0, 1, 2, 3, 4 and 5 days after treatment.
All experiments were performed using biological triplicates. The information about all the experimental conditions mentioned above are summarized in Table 1, which composed thirteen sample sets for data analysis (Table 1). For all fruit samples, fruit core was excluded and the peel and flesh were chop up, frozen in liquid nitrogen and stored at 280uC until further use.
Total RNA Isolation, Quality Control, and cDNA Synthesis All the frozen tissues were ground in liquid nitrogen for RNA isolation. Total RNA was extracted using the hot borate method of Wan and Wilkins [55] and then treated with DNAseI digestion using the RNAse-free kit (TaKaRa, Japan) to eliminate the potential DNA contamination. The RNA concentration and purity were evaluated by measuring absorbance at 230, 260 and 280 nm, respectively, using a BioPhotometer plus (Eppendorf, Germany). The integrity of the RNA samples was assessed on 1.2% agarose/formaldehyde gel electrophoresis. Only RNA samples with 260/280 ratio between 1.8 and 2.1 and 260/230 ratio higher than 2.0, as well as both 28S and 18S ribosomal RNA bands with a density ratio about 2:1 were used for further analyses. Two microgram of total RNA was reverse-transcribed using the ReverTra Ace qPCR RT kit (TOYOBO, Japan) according to the manufacturer's instructions. The final cDNA products were diluted 150-fold prior to using in RT-qPCR.

Design and Validation of Reference Gene Primers
Primer pairs were designed based on selected sequences of the 21 candidate reference genes using Primer Premier 5.0 and Primer Premier 6.0 software under default parameters. All primers were designed in 39-untranslated region (39-UTR), to ensure the specificity of amplification. Ordinarily two or more primer pairs were designed for each gene. Then primers were checked by oligo 6.0 software. All primer pairs were custom-ordered from a commercial supplier (Sangon, Guangzhou, China). Prior to the regular gene expression analysis with RT-qPCR, all primer pairs were tested by RT-qPCR to check for the specificity of the amplicon by the melting-curve after amplification with RT-qPCR analysis. Only primer pairs tested performed well which showed single product and no product in no template control (NTC) were selected for further use. The primer specificities were further confirmed with 2.5% agarose gel (TaKaRa, Japan) electrophoresis and ethidium bromide staining for a single product and the expected size. In addition, the target amplicons were sequenced to confirm specificity of the PCR products. A standard curve using a dilution series of the mixed cDNAs from all tested samples as the template (spanning five orders of magnitude) was made to calculate the gene-specific PCR amplification efficiency (E) and correlation coefficient (R 2 ) for each gene. The primer sequences and amplicon characteristics including Tm, length, amplification efficiency and correlation coefficient of the 21 candidate reference genes are listed in Table 2.

RT-qPCR Conditions
RT-qPCR was carried out in 96-well plates with Bio-Rad CFX96 Real-Time PCR System and Bio-Rad CFX96 Manager Software (Bio-Bad, USA) using SYBR Green-based PCR assay. Each reaction mix containing 5 ml diluted cDNAs, 10 ml of THUNDERBIRD SYBR qPCR Mix (TOYOBO, Japan), 0.25 mM of each primer to a final volume of 20 ml was subjected to the following conditions: 95uC for 1 min, followed by 40 cycles of 95uC for 15 s, 55uC for 30 s, and 72uC for 35 s in 96-well optical reaction plates (Bio-Rad, USA). The melting curves were analyzed at 65-95uC after 40 cycles. In addition, reverse transcription negative control was included to check for potential genomic DNA contamination. Each RT-qPCR analysis was performed in triplicate and the mean was used for RT-qPCR analysis.

Data Analysis
Two publicly available software tools, geNorm (Version3.5) and NormFinder, were used to evaluate the stability of the 21 candidate reference genes under different experimental conditions. The comprehensive ranking of these genes was generated according to a method reported previously [56,57]. An additional tool, RefFinder (http://www.leonxie.com/referencegene.php), was used to confirm the reliability of calculation. Expression levels of the tested reference genes were determined by CT values [58], the number of amplification threshold cycles needed to reach a specific threshold level of detection. Results were imported into Microsoft Excel and transformed to relative quantities. For each candidate gene, where the highest relative quantity (the minimum CT value) was set to 1, the other average CT value of each duplicate reaction of sample was converted to relative quantity data. Log-transformed data were then exported into geNorm (version 3.5) [15] and NormFinder [31], respectively, to analyze gene expression stability. The geNorm algorithm firstly calculates an expression stability value (M) for each gene and then the pairwise variation (V) of this gene with the others. All the tested genes are ranked according to their stability in the tested sample sets, and the number of reference genes necessary for an optimal normalization is indicated as well. The NormFinder program identifies the gene(s) with optimal normalization among a set of candidate genes. The lowest stability value indicates the most stable expression within the gene set examined. Therefore, it ranks the set of candidate normalization genes according to the stability of their expression patterns in a given sample set under a given experimental design.

Normalization of CpaEXY1
Endoxylanase acting on matrix polysaccharide xylan degradation is based upon gene expression that occurs during the papaya fruit ripening [51] and plays a role during papaya fruit softening [59]. CpaEXY1 encoding endoxylanase was obtained from NCBI (Genbank: AY138968.1), used as a target gene to demonstrate the usefulness of the validated candidate reference genes in RT-qPCR. Gene expression levels of CpaEXY1 were quantified during postharvest ripening using the one or two most stable reference gene(s) and the most unstable gene determined by geNorm and Normfinder in the same RT-qPCR conditions mentioned above. Primer pairs (forward:59TAATATGGTCAGCGTGGTC39,reverse:59GAGATGAGGAAGAAGGTAACT-39) of CpaEXY1 were also verified by melting curve analysis and sequencing as described for reference genes.

Selection of Candidate Reference Genes and Amplification Specificity
Twenty-one candidate reference genes based on previous reports were selected and cloned for investigation in present study. Primer pairs were designed and selected based on those candidate reference genes. Gene names, accession numbers, the primer sequences and amplicon characteristics including Tm, length, amplification efficiency and correlation coefficient of the 21 candidate reference genes are listed in Table 2. We observed that not all the initially designed primer pairs performed well in the melting curve obtained after 40 cycles of amplification, only those showing single product and no product in no template control (NTC) were selected for further experiment. Specificity of the amplifications was also confirmed by agarose gel electrophoresis, which revealed that the majority of primer pairs selected by melting curve analysis amplified a specific PCR product with the expected size. However, RT-qPCR with the primer pairs for a few candidate reference genes (ACTIN, UBQ, EF2) generated multiproducts with different size although with the same melting temperature. For these genes, we redesigned and selected and tested additional primer pairs until specific amplifications were detected by both means (Figure 1 and Figure S1). Furthermore, sequence analysis of the cloned amplicons revealed that the amplified fragments were identical or nearly identical (with 1-2 bp mismatched but the sequences of amino acids were fully identical) to the sequences used for primer design. The gene-specific PCR amplification efficiency (E) was calculated by the regression coefficient (R 2 ) of the slope of the standard curve. The PCR amplification efficiency for the 21 reference genes varied from 91.1% for APT to 106.9% for TBP2, and correlation coefficients ranged between 0.992 and 0.999 for APT and EF1, UBCE, UBQ or ACTIN, respectively (Table 2).

Expression Profile of the Reference Genes
Some variations amongst 21 reference genes were identified from the analysis of the raw expression levels across all samples ( Figure 2). The CT values of these genes ranged from 13.02 to 34.78 in all tested samples, with the majority of these CT values were between 22.51 and 29.18 ( Figure 2). The gene encoding 18S rRNA was highly expressed compared to the protein coding genes, reaching cycle threshold after only 13.02 amplification cycles, whereas the average CT value of all reference genes within the datasets was approximately 25.23 cycles. As a result, the 18S rRNA transcript levels were about 4100-fold more abundant than the dataset's average. The CT values of CHY and APT were 34.78 and 33.19, respectively, indicating the least abundant transcripts. The individual reference gene had different expression ranges across all studied samples sets. SAND and EIF showed smaller gene expression variation (below 6 cycles, with 5.91 and 5.9 cycles, respectively) among studied reference genes, while GAPDH, EF2 and CHY had much higher expression variations (above 8 cycles, with 9.06, 8.27 and 14.56 cycles, respectively). The wide expression ranges of the tested reference genes confirmed that no single reference gene has a constant expression in different sets of papaya samples. Therefore, it is of great importance to select the reliable reference gene(s) to normalize gene expression under a certain condition in papaya.

Expression Stability of Reference Genes
As no one candidate reference gene showed a constant expression in different sets of papaya samples, it was necessary to use statistical methods to rank the stabilities of the 21 genes and to determine the number of reference genes necessary for accurate gene-expression profiling under the given experimental conditions. The programs geNorm and NormFinder, two most widely used algorithms, were used in the following analysis. The RefFinder was used as verification tool.

GeNorm Analysis
The average expression stability (M) value for each candidate reference gene was calculated based on the average pair-wise variation between all genes tested. The results were presented in Table S1 and Figure 3. Stepwise exclusion of the least stable gene allowed the genes to be ranked according to their M value (the lower the M value, the higher the gene's expression stability). Among the 21 candidate reference genes used for analysis, not all of the most stable reference genes were identical in the different sample sets (Figure 3 and Table S1). For example, the EIF and RPS genes were ranked highest in different storage temperatures with an M value of 0.2786 (Figure 3a), whereas the EF1 and TBP2 genes were most stably expressed in papaya fruit samples in hot water treatment with an M value of 0.0981 (Figure 3b). For modified atmosphere packaging samples, the most stable genes were EIF and TBP1 with an M value of 0.1181 (Figure 3c), which were the same as the hormone-treated samples but with a different M value of 0.1346 (Figure 3d). The TBP2 and TBP1 genes performed best in 1-MCP treatment fruit samples with an M value of 0.1286 (Figure 3e). For the papaya samples at different development stages, the UBCE and TBP1 genes were ranked highest with an M value of 0.1734 (Figure 3f). The SAND and EIF genes were most stably expressed in different tissue samples with an M value of 0.1030 (Figure 3g). For biotic stress samples, the CYP and SAMDC genes were proved to be the best with an M value of 0.1508 (Figure 3h). For the postharvest ripening of different papaya cultivar samples, the UBCE and SAND genes were ranked highest in the cultivars sample of 'Shuiyou 29 with an M value of 0.1596 (Figure 3i), while EF1 and EF2 or TBP1 and TBP2 in cultivar samples of 'Hongri1' or 'Hongri 39 were most stably expressed, with an M value of 0.0598 or 0.0777, respectively (Figure 3 j, k). As all the three cultivars samples (for different cultivars) were taken together, the TBP1 and UBQ genes performed best with an M value of 0.2024 (Figure 3l). When all sample sets were analyzed together, the TBP1 and TBP2 were the most stably expressed genes with an M value of 0.3056 and might be widely used as a single reference gene for multiple samples (Figure 3m). In contrast, CHY, 18S rRNA and GAPDH were the three least stable among the genes examined. These results highlight the importance of choosing the appropriate reference genes according to the experimental conditions.
The geNorm program was also applied to calculate the optimal number of reference genes required for accurate normalization in the different sample sets. The software determines the pair-wise variation V n/n+1 , which measures the effect of adding further reference genes on the normalization factor (that is calculated as the geometric mean of the expression values of the selected reference genes). It is advisable to add additional reference genes to the normalization factor until the pair-wise variation V n/n+1 is inferior to a cutoff value (0.15) used by Vandesompele et al. [15], below which the added gene has no significant effect and the inclusion of an additional reference gene is not required. As shown in the Figure 4, pair-wise variation analysis suggested that normalization required the use of only two reference genes in all sample sets since the V 2/3 value of all sample sets were under the 0.15 cut-off level. Therefore, according to geNorm and Norm-Finder, the best combinations for all of the sample sets were summarized in Table 3.

NormFinder Analysis
NormFinder was also used to evaluate the expression stability of candidate reference genes. More stable gene expression is indicated by lower average expression stability values. In this mathematical model, estimation of both intra-and inter-group variation and a separate analysis of the sample subgroups in expression levels are included into the calculation of a gene expression stability value [16]. Therefore, thirteen sample sets were established as geNorm analysis. At the same time, all samples with no subgroups and the other three sample-subgroups series were analyzed using this approach as well. The results of the Normfinder analysis applied to our data sets were summarized in Table S2. It is noteworthy that definition of sample-subgroups had a notable effect on NormFinder output. However, the NormFinder output with different sample subgroups and no subgroups exhibited almost the same top six stable genes, but with the slight changes in ranking orders. When the outcomes of geNorm and NormFinder were compared, only few differences were observed except for the biotic stress samples, 'Hongri 39 samples and total samples set, which had some obvious nonconformity in the rankings (Table S1, S2).
A method previously described by Chen et al. [56] and Zhang et al. [57] was used to give a comprehensive ranking of candidate reference genes. We firstly assigned a series of continuous integers starting from 1 to 21 as weight to each reference gene, according to the reference genes ranking by each algorithm from the most stable gene to the least stable gene; then we calculated the geometric mean (GM) of each gene weights across the two methods and then re-ranked these reference genes. The gene with the less GM is viewed as more stable reference gene. The comprehensive ranking results were presented in Table S3, and the consensus of the results obtained by both geNorm and NormFinder analysis were summarize in Table 3 according to the performance. In most sample sets, geNorm or NormFinder analysis revealed almost the same top seven stable genes, although with the slight changes in ranking orders and genes. In addition, no matter how the order is changed, the most unstable gene almost remains the same in all sample sets.

RefFinder Analysis
RefFinder was used to confirm the results obtained from geNorm and NorFinder. RefFinder (http://www.leonxie.com/ referencegene.php) is a user-friendly web-based comprehensive tool developed for evaluating and screening reference genes from extensive experimental datasets. It integrates the currently available major computational programs (geNorm, Normfinder, BestKeeper, and the comparative DCt method) to compare and rank the tested candidate reference genes. Based on the rankings from each program, it assigns an appropriate weight to an individual gene and calculated the geometric mean of their weights for the overall final ranking. The CT values were input into the program directly and the ranking of the four programs and the comprehensive ranking were then calculated. The data of 1-MCP treatment and modified atmosphere packaging treatment were analyzed by RefFinder and the results were presented in Table S4, which shared a high consistency with the results evaluated by  Figure 4. Determination of the optimal number of reference genes. Pair-wise variation (V) calculated by geNorm to determine the minimum number of reference genes for accurate normalization in different experiment conditions. Arrow indicates the optimal number of genes for normalization in each sample sets. doi:10.1371/journal.pone.0044405.g004 geNorm and Normfinder described above. The top four genes and their ranking were exactly the same, although there was slight difference in the ranking of the other genes. Therefore, these results have proved that the results obtained from the two software, geNorm and Normfinder, were sufficient for our validation.

Validation for the usefulness of the Selected Reference Gene
It has been documented that the use of inappropriate references can dramatically change the interpretation of the expression pattern of a given target gene [60]. To demonstrate the usefulness of the validated candidate reference genes in RT-qPCR, the relative expression level of one papaya fruit gene, CpaEXY1, was investigated in two papaya cultivar fruit during postharvest ripening, using one or two of most stable reference genes, and the most unstable gene for normalization, which had been validated by geNorm or NormFinder as described above (Table S1, 4 and 5; Figure 3). The analysis revealed that the expression level of CpaEXY1 in 'Shuiyou2' was not obviously changed in the first four days, but was increased sharply from the 4th to the 6th days, then decreased in the later time. Similar change patterns with slight difference were also seen when UBCE alone and the combination of UBCE+CYP was used as reference gene(s) for normalization, respectively (Figure 5a). CpaEXY1 expression level in 'Hongri1' increased from the first two days, then decreased slightly and increased for a second peak during postharvest ripening and also showed similar change patterns when using CYP alone and the combination of CYP +UBQ as reference gene (s) for normalization (Figure 5b). These results were similar to those reported by a former study on CpaEXY1 gene, with northern-blot analysis [59]. However, these change patterns were completely obscured when the least stable reference gene (GAPDH) was used for normalization in 'Shuiyou 29 (Figure 5a) or CHY in 'Hongri 19 (Figure 5b). Except for the different expression trends, the expression level of CpaEXY1 normalized by CHY was 25-fold higher than that normalized by UBQ or UBQ+CYP (Figure 5b). This analysis illustrated the adverse effect of using an unsuitable reference gene for normalization and further confirmed the importance of validating reference gene stability to ensure that low precision or misleading results do not occur.

Discussion
The analysis of gene expression under different experimental conditions is a major aspect of the functional analysis of genes. Currently, one of the most commonly used technologies for gene expression analysis, which can provide more accurate data, is RT-qPCR, a method that combines high specificity and sensitivity [61]. However, quantification of gene expression is affected by several factors, such as the quantity of the initial material, the quality of the RNA, the efficiency of cDNA synthesis, primer performance, and the methods to be used for statistical analysis [38]. Among several normalization strategies that have been proposed, the use of one or more reference genes is currently the preferred way of normalization [6], and it represents a strategy that is simple to use and can control for every stage of the real-time PCR [5]. However, as no single gene has a stable expression under every experimental condition [14], it is advisable to validate the expression stability of candidate reference genes under specific experimental conditions prior to their use in RT-qPCR normalization, rather than using reference genes published elsewhere without validation [62].
Several genes including GAPDH, ACTIN, 18S rRNA, UBQ, EF, CYP and TUA have been commonly used as the reference genes for gene expression studies in many plant species [12,29,36,38,63,64]. However, recent studies have indicated that these traditionally used reference genes are not always stably expressed when tested in other species or in a wider range of experimental treatments [14,20]. For example, 18S rRNA and ACTIN have been demonstrated to be performed poorly and were less stable during the different treatments [12,65]. Therefore, it has been stated that the reference genes need to be validated for each plant species and for each specific experimental setup [66]. The results of present study further support this statement. We have demonstrated that GAPDH was not the best reference gene but the worst one for normalization during the different treatments in papaya. Furthermore, the most commonly used internal reference genes ACTIN, 18S rRNA and APT performed poorly when all the samples were taken together. Unfortunately, almost all current studies on gene expression in papaya fruit with RT-qPCR have used either Actin or 18S rRNA as the reference gene [10,11,53,54]. The poor performance of Actin or 18S rRNA as reference genes indicated the urgent need to identify the other more appropriate reference genes, which serve as the strong reasonable for the present study. To date, a large number of detailed studies have focused on reference gene selection for expression profiles in other kinds of plants [12,14,20,29,31,32,36,39,43,45,67] but none on papaya. Nevertheless, these studies have provided mass of potential reference genes and allowed the identification of suitable reference gene under a wide range of experimental conditions for papaya possible. In present work, we selected 21 candidate reference genes based on previous reports for identification of the most stably expressed reference gene(s) for normalization under thirteen different sets of experimental conditions in RT-qPCR studies. According to The MIQE Guidelines [58] and the Eleven Golden Rules of Quantitative RT-PCR [68], we tried to control all the sources of variation along the entire workflow of RT-qPCR analysis and used the two most commonly used software, geNorm and NormFinder, to evaluate the expression stability of those candidate reference genes under different experimental conditions. Thus, we conducted a comprehensive evaluation of candidate reference genes in papaya under a wide variety of conditions and treatments. This study allowed identification of the appropriate reference genes suitable for gene expression analyses under different experimental conditions. It should be pointed out that a single-peak and no product in no template control (NTC) in melting curve analysis do not necessarily mean the single product. Agarose gel electrophoresis should be also performed to confirm the specificity. In this study, we observed that a few selected primer pairs did perform well in melting curves analysis. However, multi-products were detected in the agarose gel electrophoresis analysis. In this case, additional primer pairs needed to be designed and tested until the appropriate primer pairs were found. Especially for the Actin gene obtained from NCBI, we hadn't found an appropriate primer pairs until we cloned a novel gene Actin 2, which belong to the Actin gene family. As described in the study on reference gene in banana [20], most of the reference genes tested in our study are probably members of large gene families and thus, it is difficult to obtain specific primers.
In the analysis of our datasets, we found that no single reference gene had an optimal performance across all of the experimental conditions tested. The geNorm program ranked TBP1 and TBP2 the most stable and best candidates for the normalization of general gene expression for papaya when all papaya samples were tested, and the most commonly used GAPDH and 18S rRNA were proved to be bad reference genes. Different sets of samples had their own best reference genes ( Figure 3, Table S1). For example, in the analysis of data with geNorm, reference genes EIF, RPS and ACTIN ranked higher in different storage temperature whereas EF1, TBP2 and TBP1 did better than EIF, RPS and ACTIN in the hot water treatment. EIF and TBP1 were the best reference genes for MAP and ethephon treatment. For fruit in different development stages, UBCE and TBP1 performed better than others but were not better than SAND and EIF for different tissues and CYP and SAMDC for biotic stress samples. All the differences were summarized in Table S1. Our analysis indicated that each experimental condition tested requires a specific set of reference genes. This result emphasizes the importance of reference genes validation for each experimental condition, especially when samples belong to very different sets. This idea is consistent with a number of studies by others [20,36,67].
In slight contrast to geNorm, NormFinder ranked EIF and SAND as the most stably expressed genes in the all samples' data set, which were also ranked high in the geNorm. For the other sample sets, however, there were some differences observed (Table  S1, 4), the ranked high genes almost the same although the ranking orders were somehow different slightly for most of sample sets except of the sets of biotic stress samples, 'Hongri 39 samples and the total samples, which some obviously differences were observed. In addition, no matter how the order was changed, the most unstable gene would almost remain the same in all sample sets, which had been also observed in other studies [14,20,31,60]. Several studies have also reported the similar same results with some minor changes in gene stability ranking [20,67,69]. However, others have observed relatively substantial changes between the two methods [32,60]. GeNorm and NormFinde depend on different mathematical approaches to calculate stability. GeNorm selects two genes with a low intra-group variation and approximately the same non-vanishing inter-group variation. In comparison, NormFinder selects the two best genes with minimal combined inter-and intra-group expression variation [16], which can have a notable effect on the subsequent gene stability ranking [31]. Therefore, the fact that the ranking of candidate reference genes by NormFinder is not always identical to that calculated by geNorm is not surprising.
Taking the two algorithms into account, the consensus of the results obtained by both geNorm and NormFinder analyses were listed in Table S3 and Table 3. From the result obtained by the two programs, EIF, TBP1, TBP2 and EF1 appeared to be suitable as the reference genes for papaya, due to the stability in most of sample sets obtained under different experimental conditions. These results are consistent with those reported on the whole developmental series of tomato for which the TBP exhibit a remarkable stability of expression levels [31]. As described in flax [70], both EF1 and EIF genes were the most stable reference genes in flax tissue samples. In contrast, the most commonly used reference gene such as ACTIN, 18S rRNA, GAPDH, TUA and APT were not suitable for most of experimental conditions. Several studies shared the similar results. For example, in petunia, GAPDH was considered the least stably expressed gene during leaf and flower development [28]. In tomato, GAPDH was poorly ranked as a good reference gene based on the analysis of EST data [30]. ACT2 was also found to be the least stably expressed gene among the 27 tested in Arabidopsis [13]. 18S rRNA was proved to the least reliable reference gene in peach study under different conditions [12]. TUA was found to be not stable for reference gene in the whole developmental series in tomato [31] and in different flax tissues [70], although it showed quite consistent stability in expression in several studies [35,71]. Thus the use of reference gene should be validated according to the special experimental conditions. On the other hand, considering that the reference gene and the target gene should have a similar range of expression [9], 18S rRNA might be a good choice as the reference gene for those target genes that have a relative high expression levels under some experimental conditions. However, it should be pointed out that validations of 21 reference genes by the same procedures used in present work do not always give support to their frequent use in other plants, as many studies have suggested that the reference genes are regulated differently in different plant species and might exhibit differential expression patterns [14,20]. For example, EIF or EF1 show highly stable expression in papaya (present study), Perennial ryegrass [67] and banana [20], whereas its putative homologue had been shown unsuitable for normalization in Petunia hybrida [28] and tomato [31].
Studies that fail to use appropriate reference genes could lead to bias gene expression profiles and low precision or misleading results [35,65]. To demonstrate the usefulness of the validated candidate reference genes in RT-qPCR, the relative expression level of CpaEXY1 was investigated in two papaya cultivar fruits during post-harvest ripening ( Figure 5). The results showed that normalization using the most stable reference genes (UBCE, CYP, UBQ) were coincident and similar to those reported in a former study on CpaEXY1 gene expression by northern-blot analysis [59], but the normalization was obscured when the least stable reference gene(s) (GAPDH, CHY) were used. These results further confirmed the importance of selection of stable reference genes for the correct normalization of RT-qPCR data.
To the best of our knowledge, this is the first report on the selection of appropriate reference genes using the magnitude of samples tested with different experimental conditions in papaya. Our results provide a foundation for the more accurate and widespread use of RT-qPCR in the analysis of gene expression in papaya. More importantly, our results suggest that ACTIN, 18S rRNA and GAPDH are not suitable to be used as reference genes for normalization in papaya under many experimental conditions, whereas EIF, TBP1 and TBP2 could instead serve well as reference genes due to their good performance in most of experimental conditions. Obviously, the appropriate use of these identified reference genes should be based on the given species and experimental conditions after validation. In addition, using a combination of two genes as reference genes might improve the reliability of gene expression by RT-qPCR in papaya. Figure S1 Dissociation curve data for the 21 reference genes and one target gene tested. Dissociation curves for twenty-one candidate reference genes and one target gene CpaEXY1 showed single peaks and no amplicon was observed in no template control (NTC) indicated by the pink lines. (DOC)