Genome-Wide Identification and Evaluation of Reference Genes for Quantitative RT-PCR Analysis during Tomato Fruit Development

Gene expression analysis in tomato fruit has drawn increasing attention nowadays. Quantitative real-time PCR (qPCR) is a routine technique for gene expression analysis. In qPCR operation, reliability of results largely depends on the choice of appropriate reference genes (RGs). Although tomato is a model for fruit biology study, few RGs for qPCR analysis in tomato fruit had yet been developed. In this study, we initially identified 38 most stably expressed genes based on tomato transcriptome data set, and their expression stabilities were further determined in a set of tomato fruit samples of four different fruit developmental stages (Immature, mature green, breaker, mature red) using qPCR analysis. Two statistical algorithms, geNorm and Normfinder, concordantly determined the superiority of these identified putative RGs. Notably, SlFRG05 (Solyc01g104170), SlFRG12 (Solyc04g009770), SlFRG16 (Solyc10g081190), SlFRG27 (Solyc06g007510), and SlFRG37 (Solyc11g005330) were proved to be suitable RGs for tomato fruit development study. Further analysis using geNorm indicate that the combined use of SlFRG03 (Solyc02g063070) and SlFRG27 would provide more reliable normalization results in qPCR experiments. The identified RGs in this study will be beneficial for future qPCR analysis of tomato fruit developmental study, as well as for the potential identification of optimal normalization controls in other plant species.


INTRODUCTION
Tomato (Solanum lycopersicum) is an economically important horticultural crop in terms of production, flavor, and nutritional value of fruits. During the course of development and ripening, tomato fruits undergo a number of physiological and biochemical processes that bring forth the overall changes in fruit size, color, texture, and aroma (Klee and Giovannoni, 2011;Ruiz-May and Rose, 2012). Moreover, tomato is considered as an important model for genetic and molecular studies, partly due to its typical climacteric fruit property (Colombiet et al., 2016). A number of studies had been carried out to improve agronomic traits of tomato fruits, including size, pigment content, and flavor substances focusing on the metabolic and regulatory networks (Klee and Giovannoni, 2011;Ruiz-May and Rose, 2012;Tieman et al., 2017). Recent developments in genomic resources and bioinformatics tools (e.g., Genome-wide association study, GWAS) have enabled rapid elucidation of the complicated biological processes that occur during fruit development. Moreover, relative gene expression profiles during fruit development provide valuable clues for understanding the biological functions of the corresponding genes. So far, quantitative real-time PCR (qPCR) is considered as one of the most efficient tools for the measurement of transcript abundance of a gene due to its high accuracy, sensitivity, and reproducibility (Ginzinger, 2002;Bustin and Nolan, 2004;Gachon et al., 2004;Bustin et al., 2005).
In qPCR experiments, the reliability of the results predominantly depends on the appropriateness of RGs used for normalization, which should be stably expressed under the given experimental conditions (Pfaffl, 2004;Huggett et al., 2005). Highly stable expression of RGs could effectively remove non-biological variations including the difference in amounts, variability in enzymatic efficiency of reverse transcriptase, and sample differences in the overall transcriptional activity (Suzuki et al., 2000;Bustin et al., 2005;Huggett et al., 2005;Exposito-Rodriguez et al., 2008;Gutierrez et al., 2008). Generally speaking, an ideal RG should be a gene that is stably expressed under various experimental conditions or among different tissues (Czechowski et al., 2005;Huggett et al., 2005;Exposito-Rodriguez et al., 2008;Dekkers et al., 2012;Wang et al., 2012). Housekeeping genes (HKGs) encoding, e.g., GAPDH, Actin, UBI, and 18 sRNAs, are usually regarded as suitable normalization controls (Stürzenbaum and Kille, 2001). However, some previous studies reported that the transcription of several HKGs can be fluctuated considerably under certain conditions (Czechowski et al., 2005;Jain et al., 2006;Gutierrez et al., 2008;Jian et al., 2008;Jarosova and Kundu, 2010), which illustrates the importance of systematic identification or validation of optimal RGs in order to avoid inaccurate results (Gutierrez et al., 2008;Guenin et al., 2009). In practice, the expression levels of most RGs are proved to be dependent on the specific conditions, including experimental treatments, tissue types, or developmental stages (Czechowski et al., 2005;Jian et al., 2008;Jarosova and Kundu, 2010). Hence, no single RG is widely applicable for different experimental conditions. Systematic evaluation of RGs must be conducted on each qPCR experiment before their use (Bustin et al., 2009;Jacob et al., 2013). Furthermore, it has been well-recognized that in some cases, one single RG may not be adequate for reliable normalization in gene expression analysis (Yoo et al., 2009;Cassan-Wang et al., 2012). To date, some common statistical algorithms, including geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and Bestkeeper (Pfaffl, 2004), have been developed to determine the expression stabilities of RGs, which effectively simplify the selection of appropriate RGs for qPCR analysis.
Over the decades, a good number of stably expressed RGs have specifically been identified for normalization in the fruits of several fruit crops, such as papaya (Zhu et al., 2012), blueberry (Die and Rowland, 2013), and watermelon (Kong et al., 2015). For tomato, although suitable RGs have been identified under different experimental conditions including biotic/abiotic stresses (Løvdal and Lillo, 2009;Mascia et al., 2010) and various tissues of different developmental stages (Suzuki et al., 2000;Dekkers et al., 2012), very limited number of RGs in tomato fruit have been characterized so far (Coker and Davies, 2003;Baldassarre et al., 2015). Moreover, we noticed that most studies involving RG identification, including those relevant to tomato fruit developmental studies, were mainly based on the evaluation or validation of some already known candidate RGs (Mostly HKGs), which are convenient for implementation but also greatly limit the choice of best RGs. With the availability of tomato genome sequence and subsequent transcriptome data (SGN:Sol genomics network, https://solgenomics.net/; TFGD: Tomato functional genomics database,http://ted.bti.cornell.edu/), our study was aimed to identify some novel RGs for qPCR analysis of tomato fruit development within the entire genome level.
In this study, we initially evaluated the expression stabilities of all the tomato (S. lycopersicum L.) genes during various fruit developmental stages based on the RNA-seq data. A total of 38 novel genes stably expressed were identified as putative RGs and were further evaluated through qPCR analysis. Moreover, using two different statistical algorithms (geNorm and Normfinder), five optimal RGs were identified as optimal RGs for normalization during different stages of tomato fruit development. Furthermore, we also found that the combined use of two top-ranked RGs (SlFRG03 and SlFRG27) would potentially improve the accuracy of the qPCR results. Thus, based on the analysis of the entire tomato genome database, we comprehensively identified and evaluated the optimal RGs through large-scale biological information mining and qPCR techniques. These results not only provide useful RG resources for accurate gene expression studies in tomato fruit, but also shed light on the RG identification in fruit developmental study of other plant species.

Collection and Evaluation of the Previously Reported RGs
Firstly, the potential tomato RGs reported in previous studies (Coker and Davies, 2003;Exposito-Rodriguez et al., 2008;Løvdal and Lillo, 2009;Mascia et al., 2010;Müller et al., 2015) were selected to evaluate the expression stabilities during different stages of tomato (S. lycopersicum L.) fruit development based on the RNA-seq data downloaded from the TFGD (http://ted.bti. cornell.edu/). Furthermore, the orthologous genes of 11 potential RGs for watermelon fruit development (Kong et al., 2015) were identified. All the potential RGs selected from RNA-seq data were evaluated for expression stability. Details including accession number, gene locus, gene description, and RNA-seq values were listed in Table 1 and Supplementary Table 1. The corresponding gene sequences of these candidate RGs were collected from NCBI (National Center for Biotechnology Information: https:// www.ncbi.nlm.nih.gov/), SGN (https://solgenomics.net/), TFGD (http://ted.bti.cornell.edu/) (Fraser et al., 1994), and CuGenDB (http://www.icugi.org/). Through Blastn search, the orthologous genes (E-value was set at 1e −5 ) were collected in tomato genome   Figure 1A). CV (co-efficient variation) value of each gene during fruit development was calculated as the ratio of the SD to the AVE (Supplementary Table 3).

Identification of Stably Expressed Genes through Mining Fruit RNA-Seq Data Set
Through mining the entire set of RNA-seq data file (TFGD), genes with medium expression (200 < RPKM value < 2,000) at all stages of fruit development were selected (Supplementary Table 2). The CV value of each gene was calculated (Supplementary Table 2) and the genes with CV < 0.35 were chosen as putative RGs in the following qPCR analysis ( Table 2; Supplementary Table 3).

Plant Materials
The tomato (S. lycopersicum L.) inbred line "S14" was used in the current experiments. Generally, the fruits of this genotype become completely mature (red flash) 45 days after pollination. On month old tomato seedlings were transplanted to the greenhouse at Zhejiang Academy of Agricultural Sciences, Hangzhou, China (east longitude 120 • 2 ′′ , north latitude 30 • 27 ′′ ) on August 5th (late summer), 2016. Field management was implemented following the standard commercial practices. Tomato fruits were harvested at four developmental stages

Total RNA Isolation and cDNA Synthesis
Total RNA from collected samples was isolated using TRIZOL reagent according to the manufacturer's protocol (Tiangen, Beijing, China) as described previously (Cheng et al., 2016). The concentration and purity of extracted RNA were measured using a BioDropULite spectrophotometer (Biochrom, England), and RNA samples with A260/A280>1.8 and A260/A230>2.0 (indicating good RNA quality) were used for following experiments. All RNA samples were adjusted to the same concentration in order to ensure that the RNA input was homogenized for subsequent reverse transcription reactions using a mix of random primers. Then, according to the manufacturer's instruction (TIANGEN, Beijing, China), genomic DNAs (gDNA) were eliminated from the RNA samples and single-stranded cDNAs were synthesized.

Primer Design and qPCR Analysis
Gene-specific primers were designed using Real-time PCR (TaqMan) primer design (https://www.genscript.com/) as listed in  Figure 1). The melting curves were created and exhibited for all the investigated qPCR products in the qPCR experiments (Supplementary Figure 2).

Evaluation of RG Expression Stability Using geNorm and Normfinder
The expression levels of the detected genes were obtained through the qPCR analysis and the results were demonstrated as Ct values (Supplementary Table 4). The amplification efficiency (E) and correlation coefficient (R 2 ) for each gene were calculated using the standard curve method by amplifying the 10-fold serial dilution of cDNA samples. The amplification efficiency (E) was calculated with the formula: E = (10 −1/slope -1). The geNorm and NormFinder software packages were used to evaluate the gene expression stability as described before in this study (Vandesompele et al., 2002;Andersen et al., 2004). The geNorm applet not only provides a measure of gene expression stability value (M), but also creates pairwise variation values (V) to determine the minimum number of RGs required for reliable normalization, no additional genes were required when the pairwise variation (Vn/n+1) was lower than 0.15 (Vandesompele et al., 2002). The NormFinder measured the variations across groups and determine the expression stability of each tested gene (Andersen et al., 2004). Lower stability values (M) in both geNorm and Normfinder implied the higher expression stability of the genes. The stability values (M) obtained from geNorm and NormFinder were listed in Supplementary Table 5.

Evaluation of Previously Reported RGs during Tomato Fruit Development
In our study, tomato (S. lycopersicum L.) genome database in the SGN (https://solgenomics.net/) and transcriptome data derived from the TFGD (http://ted.bti.cornell.edu/) were used for  analysis. Based on the previously reported RGs in tomato (Coker and Davies, 2003;Exposito-Rodriguez et al., 2008;Løvdal and Lillo, 2009;Mascia et al., 2010;Müller et al., 2015), 73 reported potential RGs of tomato were identified (Table 1). Moreover, by collecting the cDNA sequences of 11 previously reported candidate RGs for watermelon fruit study (From CuGenDB: Cucurbit genomics database, http://www.icugi.org/) (Kong et al., 2015), we subsequently collected their corresponding orthologs in tomato through blastN in the SGN (Table 1). Thus, a total of 84 tomato potential RGs were collected in the present study. Among all the reported RGs and their corresponding orthologs, three genes (Solyc01g028930, Solyc10g049850, and Solyc03g115810) were found to be redundant. Thus, a total of 81 previously reported RGs were collected and listed in (1088.0), Solyc11g005330 (581.1), Solyc01g028810 (509.5), Solyc04g009770 (357.6), and Solyc05g023800 (359.9), were high enough (>200) to be considered as candidate RGs in tomato fruit. Hence, we came to the conclusion that most of the previously reported candidate RGs identified so far were not well-qualified for normalization during tomato fruit development.

Identification of Putative RGs Based on RNA-Seq Data Mining
To comprehensively identify qualified RGs for tomato fruit, we searched the entire set of RNA-seq data derived from TFGD, and a total of 56 genes with RPKM values ranged between 200 and 2,000 among all developmental stages of tomato fruit were identified, by searching the derived RNA-seq data (Supplementary Table 2). The CV values of these genes were calculated, and more than 70% (40/56) of them were shown to be lower than 0.35 (Supplementary Table 3). The 40 genes were listed in details in Table 2. Further analysis revealed that the genes identified from the entire genome level ( Figure 1B) were generally more stably expressed than those previously reported as candidate RGs ( Figure 1A) during tomato fruit development.

qPCR Analysis of the Putative RGs in Tomato Fruit Development Process
Next, we intended to validate the expression stabilities of the selected 40 putative RGs by qPCR analysis. When designing the primers using genescript online tool (https://www.genscript. com/), we found that proper primers of two candidate RGs (Solyc01g095050 [321 bp] and Solyc10g074860 [72 bp]) for qPCR analysis cannot be designed due to their short cDNA sequences or high homologies with other genes in tomato. Thus, a total of 38 candidate RGs, which are designated as SlFRG01 to SlFRG38, were eventually chosen for further expression validation in the qPCR experiments ( Table 2). qPCR amplification of the 38 candidate RGs were carried out using specific primers listed in Table 2, and the amplicon lengths ranged from 67 bp (SlFRG21) to 149 bp (SlFRG07). PCR-amplification specificity of each primer pair was verified by 1.5% agarose gel electrophoresis using cDNA templates (Supplementary Figure 1), and the melting curve analysis also showed single product peak (Supplementary Figure 2), which both confirm the specificities of the primer pairs. The amplification efficiencies (E) of these candidate RGs were found to vary from 0.76 (SlFRG01) to 1.39 (SlFRG22), and E values of more than 50% primer pairs (20/38) were ranged from 0.9 to 1.1, indicating their qualifications as primer pairs (Tichopad et al., 2002;Chung et al., 2004). Notably, the amplification efficiencies (E) of SlFRG07 and SlFRG13 could not be calculated due to their low transcript level in tomato fruit (Supplementary Table 4). The correlation coefficients (R 2 ) of the 38 candidate RGs ranged from 0.915 (SlFRG37/SlFRG38) to 1 (SlFRG06/SlFRG33) ( Table 2). The
Although the two different assessing systems (geNorm and Normfinder) came up with different results, there were still seven putative RGs (SlFRG27, SlFRG04, SlFRG35, SlFRG05, SlFRG37, SlFRG16, SlFRG12) that were found to be commonly top-ranked in both statistical algorithms ( Table 3). Generally speaking, the primer pairs with amplification efficiency (E) between 0.9 and 1.1 (Tiangen, China) possess the lowest variability in qPCR analysis (Tichopad et al., 2002;Chung et al., 2004). However, we found that among the seven top-ranked RGs, the amplification efficiencies (E) of two primer pairs (SlFRG04 [0.78] and SlFRG35 [0.87]) were lower than 0.9, suggesting that the primer pairs of these two genes were not recommended for subsequent RG application in this study. Ct values of the remaining five RGs were all between 20 and 25 (qualified as RG). Therefore, SlFRG27, SlFRG05, SlFRG37, SlFRG16, and SlFRG12 were finally identified as qualified and optimal RGs for normalization in the tomato fruit developmental process.
Previously, some researchers have reported that the use of more than one internal control genes in normalization could effectively improve the reliability of qPCR results (Reid et al., 2006;Exposito-Rodriguez et al., 2008;Gutierrez et al., 2008). Thus, we applied the geNorm software to calculate the pairwise variation values (V) of the 38 putative RGs (described in detailed in the Materials and Methods section). The pairwise variation revealed that the V2/3 value was 0.06 (significantly < 1.5) (Figure 4), which indicated that the combined use of two most stably expressed RGs as reflected in geNorm, SlFRG03 and SlFRG27, was potentially sufficient for better normalization in qPCR experiments of tomato fruit developmental studies (Figure 4).

DISCUSSION
The advent of qPCR technology has brought a new revolution in the gene expression analysis area. Accurate interpretation of qRCR results mainly depends on the use of stable RGs for normalization, which can potentially minimize nonbiological variations of different samples. Hence, the systematic FIGURE 2 | Tomato fruit samples of four representative developmental stages. IM, Immature, 15 days after fertilization, 5 cm diameter fruit; MG, Mature green, 30 days after fertilization, 9 cm diameter fruit; B, Breaker, 35 days after fertilization, 10 cm diameter fruit; MR, Mature red, 45 days after fertilization, 10 cm diameter fruit.

FIGURE 3 | Expression stability of the 38 newly identified RGs evaluated by geNorm (A) and NormFinder (B). (A)
Ranking of geNorm is based on the principle that logarithmically transformed gene expression ratio between two ideal internal control genes should be identical if both genes are stably expressed in the tested sample set. Expression stability values (M) of the 38 candidate RGs are shown, RGs with a higher M value are less stably expressed. (B) NormFinder is a model-based approach that evaluates expression variation by comparing the variation within and between a certain number of sample groups and RGs with lower combined levels of intra and intra-group variation were regarded to be more stably expressed.
identification of appropriate RGs is essential for obtaining reliable results in qPCR experiments (Udvardi et al., 2008;Bustin et al., 2009;Guenin et al., 2009). Nowadays, some HKGs (e.g., Actin, Ubiquitin, and 18s rRNA), are usually used as RGs under various experimental conditions, or across a broad range of tissue samples (Bustin, 2002;Kong et al., 2014). However, an increasing number of evidence showed that optimal RGs varied depending on the experimental conditions or organs/tissues assayed, and it seems to be impossible to acquire a list of RGs universally practicable across a wide range of experimental conditions (Guenin et al., 2009;Warzybok and Migocka, 2013;Kong et al., 2014). Therefore, the identification of suitable RGs for specific experimental conditions is essential for avoiding unnecessary error in the qRCR experimental results. So far, many studies involving the identification or evaluation of tomato RGs under various experimental conditions, including biotic/abiotic stresses (Cucumber mosaic virus, tobacco mosaic virus, bacterium Xanthomonas campestris, nitrogen stress, cold, light stress) (Coker and Davies, 2003;Løvdal and Lillo, 2009;Mascia et al., 2010;Wieczorek et al., 2013;Müller et al., 2015) and different organs/tissues (Leaf, fruit, flower, and seed) (Exposito-Rodriguez et al., 2008;Dekkers et al., 2012;Gao et al., 2012;Baldassarre et al., 2015), had been conducted. For example, SlACT, SlCAC, and SlEF1α were validated to be suitable RGs in studies of host-virus interactions in tomato (Wieczorek et al., 2013). Exposito-Rodriguez et al. (2008) found that the widely used RGs, such as SlCAC, SlTIP41, Expressed, and SlSAND, provide superior transcript normalization in various tissues of tomato (Exposito-Rodriguez et al., 2008). When studying the changes of gene expression in the wounded ripening-stage tomato fruit, Baldassarre et al. (2015) selected two most stably expressed RGs (EF1-α and GADPH) from seven routine used HKGs for normalization in the subsequent qPCR analysis. Thus, it occurred to us that although tomato is regarded as model plant for fruit development study, little attention has yet been paid to screen the best RGs specifically for normalization during the development of tomato fruit. So far, most studies involving RG identification were based on the evaluation or validation of the expression stabilities of traditional or novel RGs under corresponding conditions (Czechowski et al., 2005;Løvdal and Lillo, 2009;Schmidt and Delaney, 2010;Dekkers et al., 2012;Baldassarre et al., 2015). In the present study, we collected 70 putative RGs that had been previously reported in tomato (Coker and Davies, 2003;Exposito-Rodriguez et al., 2008;Løvdal and Lillo, 2009;Mascia et al., 2010;Müller et al., 2015) and 11 orthologs of reported RGs in watermelon fruit study (Kong et al., 2015), and subsequently validated their expression stabilities during different stages of fruit development according to the RPKM values derived from RNA-seq data sets. Out of expectation, the majority of these putative RGs identified previously were not well-qualified for normalization as internal control genes due to either their unstable expressions or low transcript levels (Table 1; Figure 1A; Supplementary Table 1). Therefore, we next intended to identify some novel RGs that are stably expressed during the whole developmental process of tomato fruit.
Nowadays, the open sources of the SGN and TFGD allowed us to search for the most stably expressed genes on a comprehensive evaluation system. In the present study, we conducted a data mining based on the RNA-seq data set of different tomato fruit developmental stages, and initially collected 38 most stably expressed genes as putative RGs for normalization in tomato fruit developmental study ( Table 2). Further evaluation of these 38 putative RGs was conducted using qPCR analysis in four developmental stages of tomato fruit (IM, MG, B, and MG) (Figure 2). Next, we used two popular statistical algorithms for RG ranking, geNorm, and Normfiner (Vandesompele et al., 2002;Andersen et al., 2004), for the RG evaluation based on the qPCR results (Supplementary Table 4). We found different ranking results from the evaluation results of geNorm and Normfinder (Table 3), which were explicable as these two algorithms are based on different models and assumptions (Schmidt and Delaney, 2010). The geNorm algorithm is based on the principle that logarithmically transformed expression ratio of two genes should be constant if both of them are stably expressed in the tested sample set. The relative stability of each gene (M) is defined as the mean pairwise variation (reflected by standard deviation of the expression ratios of two genes) of the gene in the sample set. Furthermore, as normalization with single RG can cause inevitable errors, geNorm is also used to determine the minimum number of RGs required for more reliable normalization (Vandesompele et al., 2002). Normfinder measures gene expression stability by comparing the variation within and between a certain number of sample groups. The genes with the lowest combined levels of intra and inter-group variation were regarded as most stably expressed (Andersen et al., 2004;Schmidt and Delaney, 2010). Taken together, Normfinder is based on analyzing the variation level of each tested gene rather than pairwise analysis of gene stability relative to a set of potential RGs (Schmidt and Delaney, 2010). So far, numerous ranking differences of RGs derived by these two algorithms had been found in many previous studies (Schmidt and Delaney, 2010;Cassan-Wang et al., 2012;Kong et al., 2015). Nevertheless, we identified 7 putative RGs (SlFRG04, SlFRG35,SlFRG27, SlFRG05, SlFRG37, SlFRG16, and SlFRG12) that were common in geNorm and Normfinder. Considering the unqualified primer amplification efficiencies (E) of SlFRG04 and SlFRG35, these two genes were excluded from our recommendation list, and the remaining five genes (SlFRG05, SlFRG12, SlFRG16, SlFRG27, and SlFRG37) were finally identified as suitable internal controls for normalization in tomato fruit development. Notably, we believe that some alternative primer pairs of SlFRG04 and SlFRG35 with improved amplification efficiencies might be redesigned for RG use in tomato fruit developmental study.
In practice, it is believed that the use of more than one RG in the normalization can efficiently improve the reliability of qPCR results (Alba and Giovannoni, 2005;Exposito-Rodriguez et al., 2008;Gutierrez et al., 2008). Thus, in order to explore the minimun number of RGs needed, the pairwise variation (V) values were calculated in geNorm (Figure 4). According to the evaluation, the combined application of two RGs, SlFRG03, and SlFRG27, would be a better choice than the use of only one RG for normalization when more reliable qPCR results are needed. It is also worth noting that due to the multiple sections of tomato fruit and the complex biological processes of fruit development, gene expression analysis has been extended to more precise tissue parts (e.g., pericarp, flesh, and even seeds) or longer developmental stages of fruits (Fraser et al., 1994;Carrari and Fernie, 2006; FIGURE 4 | Analysis of best RG association based on geNorm algorithm. The optimum number of RGs is the lowest number of genes with an acceptably low Vn/n+1. Vandesompele et al. (2002) suggested 0.15 (15% variation in normalization factors) to be an upper limit for Vn/n+1. According to variations (V-value) calculation, V2/3 is 0.06 (<0.15), which means the most stable expressed RGs identified ingeNorm, SlFRG03 and SlFRG27, are well-qualified as RG combination for normalization. Fei et al., 2011;Cheng et al., 2016). Therefore, we propose that the RGs identified in this study should be further validated in different tissue sections or earlier developmental stages (e.g., 1, 2 cm green fruits) of tomato fruit in the future according to specific experimental requirements.

CONCLUSION
To our knowledge, this study is the first systematic identification and evaluation of putative RGs as internal controls for normalization of qPCR analysis in tomato fruit developmental process. According to our extensive evaluation, five identified RGs-SlFRG05, SlFRG12, SlFRG16, SlFRG27, and SlFRG37 could be recommended for normalization of qPCR experiments in tomato fruits. Furthermore, according to geNorm analysis, a combination of two most stably expressed RGs, SlFRG03 and SlFRG27, were recommended when more reliable qPCR results were needed. Moreover, by comparative analysis of the previously published materials involving RG identification for fruit developmental study in other plants, we found that two RGs identified in this study were also chosen as optimal RGs for fruit developmental study in other plants (Zhu et al., 2012;Die and Rowland, 2013;Kong et al., 2015), which are ubiquitin conjugating enzyme (UBI) encoding genes (SlFRG27 in tomato, PEX4 and UBC28 in blueberry, UBCE in papaya) and actin encoding genes (SlFRG37 in tomato, ClACT in watermelon, ACTIN in papaya). Thus, SlFRG27/SlFRG37 and their corresponding orthologs seem to be universally applicable as RGs among plants of different families including Cucurbitaceae, Rosaceae, Vacciniaceae, and Solanacea. Taken together, the results presented here not only unveil optimal RGs for qPCR analysis in tomato fruit development, but also provide referable guidelines for identification of RGs in other plant species.