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.

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 non-biological variations of different samples. Hence, the systematic 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; 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.

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.

Conceived and designed the experiments: HW and YC; Performed the experiments: YC, WB, XP, JY, MR, QY, RW, ZL, GZ, and ZY; Analyzed the data: JY, WB, XP and GZ; Wrote the paper: YC, WB, XP, GA, and HW. All authors have read and approved the manuscript.