DUAN Yao-ke ,HAN Rong ,SU Yan ,WANG Ai-ying ,LI Shuang ,SUN Hao ,GONG Hai-jun

1 Shaanxi Engineering Research Center for Vegetables/College of Horticulture,Northwest A&F University,Yangling 712100,P.R.China

2 Henan Key Laboratory of Ion-Beam Bioengineering/College of Agricultural Sciences,Zhengzhou University,Zhengzhou 450000,P.R.China

Abstract Tomato (Solanum lycopersicum) is a model plant for research on fruit development and stress response,in which gene expression analysis is frequently conducted. Quantitative PCR (qPCR) is a widely used technique for gene expression analysis,and the selection of reference genes may affect the accuracy of results and even conclusions. Although there have been some frequently used reference genes in tomato,it has been shown that the expressions of some of these genes are not constant in different tissues and environmental conditions. Moreover,little information on genomic identification of reference genes is available in tomato. Here,we mined the publicly available transcriptional sequencing data and screened out fifteen candidate reference genes,and the expression stability of these candidate genes and seven traditionally used ones were evaluated under stress and hormone treatment. The results showed that over half of the selected candidate references were housekeeping genes in tomato cells. Among the candidate reference genes and the traditionally used ones,the most stably expressed genes varied under different treatments,and most of these genes were recommended as preferred reference genes at least once except Solyc04g009030 and Solyc07g066610,two traditionally used reference genes. This study provides some novel reference genes in tomato,and the preferred reference genes under different environmental stimuli,which may be useful for future research. Our study suggests that excavating stably expressed genes from transcriptome sequencing data is a reliable approach to screening reference genes for qPCR analysis.

Keywords: tomato (Solanum lycopersicum),gene expression,qPCR,reference gene,expression stability

1.Introduction

Gene is the basic unit for genetic information storage in plants,and the transcriptional products initially determine complex signaling and metabolic pathways in the cell,which forms the basis for plant development and environmental responses (Shinozaki and Dennis 2003;Aslamet al.2022;Wanget al.2022). To determine the transcriptional abundance of genes,several approaches have been developed. These include Northern Blot(Alwineet al.1977),semi-quantitative PCR (semi-qPCR;Wanget al.1989),quantitative real-time PCR (qPCR;Gibsonet al.1996),digital quantitative real-time PCR(dPCR;Vogelstein and Kinzler 1999),microarray (Fodoret al.1991),and RNA-seq (Wanget al.2009). Among these approaches,qPCR has been most widely used due to its advantages of high sensitivity and specificity,cost-effect,and a broad quantification range (Gachonet al.2004). However,problems like potential inherent variability of RNA,extraction protocols,and reverse transcription efficiency mean that proper normalization is critical to objectively reflect the gene expression level (Bustin and Nolan 2004). Strategies including normalization to sample size,total RNA,or relative expression with internal reference genes have been developed,with the relative expression to one or more reference genes being the most popular practice (Huggettet al.2005;Wong and Medrano 2005).

The selection of proper reference genes determines the accuracy of the qPCR results. An ideal reference gene should be expressed constantly regardless of tissue or cell types,developmental stage or sample treatment,and the expression should not be too high or too low(Huggettet al.2006;Bustinet al.2013). Previously used reference genes are mostly introduced from wellstudied housekeeping genes as they are assumed to be expressed in all tissues or cells,and their expression is unaffected by environmental conditions (Gutierrezet al.2008). However,it has been shown that the expressions of these traditional reference genes are not constant in different tissues,organs,or environmental conditions (Masciaet al.2010;Hruzet al.2011). For instance,Expósito-Rodríguezet al.(2008) reported that the frequently used glyceraldehyde-3-phosphate dehydrogenase (GAPDH),elongation factor 1-alpha(EFα1),and alpha-tubulin (TUA) demonstrated unstable expressions during tomato development. InSetariaviridis,bothGAPDHandactin(ACT) were considered the least stable reference genes during development and abiotic stress response (Martinset al.2016). Using inappropriate reference genes may cause ignorance of slight distinction or even incorrect conclusions (Vandesompeleet al.2002;Huggettet al.2005).

Traditionally,the selection of reference genes always depended on evaluating previously used candidate genes,introducing reference genes from other species,or excavating stably expressed genes from expressed sequence tag (EST) (Coker and Davies 2003;Expósito-Rodríguezet al.2008;Wanet al.2010;Zhanget al.2020;Liet al.2021;Wuet al.2021;Yanget al.2021;Yinet al.2021;Zhuet al.2021),which overwhelmingly relied on previous research and thus restricted the excavation of reference genes. The rapid development of whole transcriptome sequencing techniques like RNAseq not only makes a high-throughput analysis of gene expression available,but also provides a mass of data in different treatments or experiments,where screening and evaluating reference genes for qPCR analysis is possible (Smithaet al.2019). However,up to date,this strategy has only been adopted in very few studies. For instance,using this approach,Gaoet al.(2018),Gonget al.(2016),Smithaet al.(2019),and Tonget al.(2021)have successfully identified some novel reference genes inPyropiayezoensis,goji,Gossypiumhirsutum,and Chinese chive,respectively;and they found that the expressions of these genes were more stable than those of the traditionally used ones,indicating that the selection of stably expressed genes from transcriptome sequencing data is a powerful and reliable approach to excavating reference genes.

Tomato (Solanumlycopersicum) is an important model plant in genetic and molecular research,and it is also a widely planted vegetable crop worldwide. Its growth is inhibited,and the yield is reduced by various biotic and abiotic stresses like drought,salinity,pathogen,and pest(Gerszberg and Hnatuszko-Konka 2017;Baiet al.2018;Biondiet al.2018;Nemeskéri and Helyes 2019;Jahanet al.2021). Various genes that consist of ingenious and complex signal pathways have been found to be involved in the responses of tomato plants to environmental stresses (Redmanet al.2001;Rivas and Thomas 2005;Shinozaki and Yamaguchi-Shinozaki 2007;Sunet al.2010;Habibet al.2021). In gene expression analysis in tomato using qPCR,previously often used reference genes included the genes encoding phosphatase (GAPDH and calcineurin B-like protein),kinase (phosphoglycerate kinase and UMP-CMP kinase 3),membrane protein(TIP41-like protein),and receptors (cyclophilin) (Expósito-Rodríguezet al.2008;Løvdal and Lillo 2009;Masciaet al.2010). These reference genes were mostly introduced from other species or excavated from tomato EST data.For instance,TIP41-like protein gene was introduced from grapevine (Reidet al.2006),UMP-CMP kinase 3 and Cyclophilin (Cyp) genes were previously employed inArabidopsis(Czechowskiet al.2005),while GAPDH and phosphoglycerate kinase genes were excavated from tomato EST data (Coker and Davies 2003). Recently,Pomboet al.(2017) screened reference genes from the RNA-seq data on tomato leaves infiltrated with different immunity inducers and bacterial strains. They found thatARD2andVIN3,but not other previous frequently used genes,were the most stably expressed genes.Mülleret al.(2015) evaluated reference genes fromXcvinfected tomato leaves based on microarray analyses and found that the most suitable reference genes werePHDandLSm7. After investigating the transcriptome data of fruit samples from different developmental stages,Chenget al.(2017) discovered thatSlFRG05(Solyc01g104170),SlFRG12(Solyc04g009770),SlFRG16(Solyc10g081190),SlFRG27(Solyc06g007510),andSlFRG37(Solyc11g005330) were the most stable reference genes for gene expression analysis during tomato fruit development.

However,stable reference genes in tomato under various environmental conditions remains to be excavated. In this study,we firstly analyzed tomato transcriptome with large-scale sequencing data and screened out fifteen candidate reference genes;then,we evaluated the candidates under stress and hormone treatment conditions. This study not only provides novel reference genes in tomato,but also provides the preferred reference genes under different environmental stimuli,which may be used in future research.

2.Materials and methods

2.1.ldentification of candidate reference genes

Both traditionally used reference genes and new ones from high-throughput data were evaluated in this research. All the previously evaluated or used reference genes were retrieved from Web of Science with the following searching parameters: TI=(tomato ORSolanumlycopersicumORLycopersiconesculentum)AND (reference gene OR internal control). The articles evaluating reference genes for qPCR in tomato were manually screened out. The traditional reference genes were selected from these articles based on their expression stability in the literature (Expósito-Rodríguezet al.2008;Løvdal and Cathrine 2009;Masciaet al.2010;Rezzonicoet al.2018). To screen candidate tomato reference genes from high-throughput transcriptome sequencing data,a software called Genevestigator (NEBIOB,Swiss,https://genevestigator.com) was employed (accessed in December 2020).The organism was limited to tomato and filtered to wildtype genetic background. The selection was refined by conditions including ‘biotic’,‘chemical’,‘elicitor’,‘hormone’,‘light intensity’,‘light quality’,‘nutrient’,‘photoperiod’,‘stress’,‘temperature’,and ‘others’,which resulted in 529 samples being selected in total. A total of 20 genes that evidenced a median expression level(signal intensity at the median of all the samples) and a low standard deviation were screened out,of which 15 ones containing a specific sequence region (suitable for primer design) were identified as the candidate reference genes.

2.2.GO and KEGG analyses of reference genes

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of all the reference genes were performed with an online tool(https://david.ncifcrf.gov/).

2.3.Plant materials preparation and treatment

After sufficient water absorption,the seeds of tomato(‘Alisa Craig’) were sterilized in water at 55°C for 15 min and then washed with sterile water. The seeds were placed in a Petri dish and cultured at 30°C in dark for 2 d,after which the seeds were sown in plugs filled with substrate. Two-leaf-stage seedlings were transferred to black-painted plastic boxes containing 1/4 Hoagland nutrient solution (pH 5.8,Hoagland and Arnon 1950),which contained 1.25 mmol L-1KNO3,1.25 mmol L-1Ca(NO3)2,0.5 mmol L-1MgSO4,0.25 mmol L-1KH2PO4,71 μmol L-1FeNa-EDTA,46.1 μmol L-1H3BO3,0.24 μmol L-1(NH4)6Mo7O24,9.1 μmol L-1MnCl2,0.76 μmol L-1ZnSO4,and 0.32 μmol L-1CuSO4. After 7 d of culture,the nutrient solution was changed to 1/2 Hoagland solution(CT) or solution containing 20% (w/w) polyethylene glycol 6000 (PEG),150 mmol L-1sodium chloride (NaCl),50 μmol L-1abscisic acid (ABA),1.5 mmol L-1hydrogen peroxide (H2O2),and 200 μmol L-1salicylic acid (SA).The employed concentrations of ABA,H2O2,and SA were based on Duet al.(2013),Linet al.(2018),and Mandalet al.(2009),respectively. After 3,12,and 24 h of treatments,the roots and leaves were collected and frozen in liquid nitrogen for RNA extraction and gene expression analysis.

2.4.mRNA isolation and synthesis of first-strand cDNA

DNA-free mRNA of tomato roots and leaves was extracted with E.Z.N.A.®Plant RNA Kit (Omega,GA,USA) under the guidance of the instructions,and the RNA quality and concentration were determined by Nanodrop (Thermo,WI,USA). The first strand of cDNA was synthesized with equivalent RNA by HiScript®II Q RT SuperMix for qPCR (+gDNA wiper,Vazyme,Nanjing,China) according to the instructions. Briefly,1 μg mRNA was added as the template in a 20-μL reaction system,and the cDNA was synthetized at 37°C for 15 min,followed by incubation at 85°C for 5 s.

2.5.Primer amplification efficiency of candidate reference genes

Primers of the previously used reference genes were adopted from the literature (Table 1). Primers of the new candidate reference genes from high-throughput transcriptional data were designed with the online Primer-Blast tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome),and the primers were designed on different exons with an intron of more than 500 bp if possible. All the primers are listed in Table 1.

To determine the primer amplification efficiency of candidate reference genes,we first performed PCR and purified the products. The PCR products were then respectively diluted at 1:10 000,1:80 000,1:640 000,and 1:5 240 000,and equivalent dilutions were employed as templates in qPCR. The qPCR was performed on a quantitative PCR system (K8930AA,AriaDx Real-Time PCR,CA,USA). The amplification efficiency of primers was calculated based on the linear fitting of quantification cycle (Cq) values and the logarithm (base two) of template concentration. To evaluate the amplification specificity of primers,melting curves of the qPCR products were run and checked.

2.6.Expression stability of candidate reference genes

The expression stability of candidate reference genes under different stress and hormone treatments was detected with qPCR. The reaction system was prepared under the instructions of ChamQTMSYBR®qPCR Master Mix (Vazyme,Nanjing,China). Briefly,the reaction mix was prepared with 10 μL of 2× ChamQTMSYBR®qPCR Master Mix,0.4 μL forward primer,0.4 μL reverse primer,0.5 μL cDNA,and ddH2O added to 20 μL. The PCR procedure was set as follows: pre-denaturation at 95°C for 30 s,40 cycles of denaturation at 95°C for 10 s,and annealing-extension at 60°C for 30 s. The melting curve was set as the default setting of the quantitative PCR system (K8930AA,AriaDx Real-Time PCR,CA,USA).

To evaluate the expression stability of candidate reference genes,the Cq values were assigned to the following four well-established algorithms for analysis: the comparative ΔCq method (Silveret al.2006),BestKeeper(Pfafflet al.2004),NormFinder (Andersenet al.2004),and GeNorm (Vandesompeleet al.2002). The integrative analyses of expression stability were performed with the publicly available web-based platform (https://www.heartcure.com.au/reffinder/?type=reference#),which builds up a combined stability index based on the ranking of the above four algorithms (Xieet al.2012).

3.Result

3.1.ldentification of candidate reference genes

The qPCR technology has been more and more frequently employed to investigate gene expression,as seen from the continuous increase of relevant published papers(Fig.1). The selection of suitable reference genes is critical for qPCR analysis. To evaluate the traditional reference genes for qPCR in tomato,30 reference genes that demonstrated relatively stable expression in the literature were collected (Table 1). A good number of them are housekeeping genes like those encoding actin (Solyc03g078400,Solyc04g011500,and Solyc11g005330),ubiquitin or ubiquitin-conjugating enzyme(Solyc01g068045,Solyc07g064130,Solyc06g007510,and Solyc11g071260),elongation factor 1-alpha(Solyc06g005060 and Solyc11g069700),and protein phosphatase subunit (Solyc05g006590,Solyc05g009600,and Solyc06g069180). After the expression stability was evaluated based on previous research and plant tissues in this study,seven traditional reference genes were selected (Table 2;Appendix A). These reference genes include Solyc04g009030 (encoding GAPDH) (Løvdal and Lillo 2009),Solyc07g066610 (encoding phosphoglycerate kinase) (Løvdal and Lillo 2009),Solyc10g049850 (encoding TIP41-like protein) (Expósito-Rodríguezet al.2008;Dekkerset al.2012;Lacerdaet al.2015),Solyc07g025390(encoding an uncharacterized protein LOC101263039)(Expósito-Rodríguezet al.2008;Lacerdaet al.2015;Gonzalez-Aguileraet al.2016),Solyc12g015870 (encoding calcineurin B-like protein) (Rezzonicoet al.2018),Solyc01g088480 (encoding UMP-CMP kinase 3) (Masciaet al.2010),and Solyc01g111170 (encoding cyclophilin)(Masciaet al.2010).

To excavate novel reference genes from tomato transcriptome sequencing data,the software Genevestigator was employed. With parameters set as described in Section 2,20 genes that demonstrated low standard deviations and medium expressions were screened out from 529 samples in total. Since Solyc09g091945,Solyc07g053525,Solyc10g038005,and Solyc09g098445 could not be accessed in the Sol Genomics Network (https://www.sgn.cornell.edu/),and the expression of Solyc06g009970 was very low,these five genes were excluded. Therefore,15 novel candidates with specific primers were introduced in this research(Table 2). Of these new candidates,Solyc06g009210,Solyc10g078960,and Solyc02g064670 encoded 60S ribosomal protein;Solyc01g100870,Solyc01g008000,and Solyc01g100860 encoded ADP-ribosylation factor;Solyc07g025400 and Solyc01g088040 werepredicted to be non-coding RNA;Solyc01g095050 and Solyc12g100060 respectively encoded a negatively light-regulated protein and a stress-associated protein;Solyc07g064150 encoded a protein translation factor SUI1 homolog;Solyc11g039840 encoded a ubiquinolcytochrome c reductase iron-sulfur subunit,and Solyc01g007860 encoded an ubiquitin-conjugating enzyme (Table 2). In addition,Solyc03g115490 and Solyc04g078070 encoded a LYR motif-containing protein and an uncharacterized protein,respectively (Table 2).

Table 1 Previously evaluated reference genes in tomato

3.2.Performance of amplification primers

The performance of primers was evaluated using amplification efficiency and amplification specificity.To determine the amplification efficiency,qPCR was conducted. The results showed that the amplification efficiency of all the primers was between 94.4 and 105.0%,with the means being 100.2% (Table 2). To evaluate the amplification specificity of primers,the melting curves of amplification products were run and checked. The results showed that there was one single peak in every qPCR analysis,indicating that there was only one amplification product (data not shown). These results suggest that all the primers had high amplification efficiency and specificity.

Table 2 Detailed information of traditionally used reference genes and selected novel candidate reference genes based on transcriptional sequencing data

3.3.GO function and KEGG pathway analysis of reference genes

GO function and KEGG pathway analysis of the candidate reference genes revealed that Solyc04g009030,Solyc11g039840,Solyc01g088480,and Solyc07g066610 were involved in various metabolic pathways,including carbon metabolism,pyrimidine metabolism,biosynthesis,and oxidative phosphorylation;Solyc02g064670,Solyc10g078960,and Solyc06g009210 encoded large and small ribosomal subunits;Solyc01g100870,Solyc01g100860,and Solyc01g008000 might function in cell endocytosis;Solyc07g064150 might help transport RNA;Solyc01g111170 was involved in protein folding;Solyc01g007860 was involved in postreplication repair and was associated with protein K63-linked ubiquitination(Fig.2). It can be seen that these 13 genes (over half of 22) were required for the maintenance of basic cellular functions and thus could be considered housekeeping genes.

However,there were nine genes that could not be enriched into a KEGG pathway or whose GO functions could not be identified. Among these,Solyc04g078070 and Solyc07g025390 encoded uncharacterized proteins,Solyc01g088040 and Solyc07g025400 encoded ncRNA,Solyc01g095050 and Solyc12g100060 were light-regulated and stress-associated protein-encoding genes,respectively. The functions of Solyc03g115490,Solyc10g049850,and Solyc12g015870 are still unclear.

3.4.Expression stability of candidate reference genes

In order to evaluate the stability of the candidate reference genes,we analyzed their expressions in both leaves and roots using qPCR under water deficit (simulated by PEG),salt stress,oxidative stress (treated by H2O2),ABA,and salicylic acid treatment. The candidate reference genes displayed a relatively wide range of expression levels.For instance,Solyc01g111170,Solyc07g064150,and Solyc10g078960 were the three genes with the lowest Cq values,suggesting their high expressions;while Solyc07g025400,Solyc01g095050,and Solyc01g088040 were the three genes with top Cq values,implying their lower expressions (Fig.3).

The candidate reference genes not only had differences in expression levels,but also demonstrated distinct expression stability under various stress and hormone treatments (Fig.3). To evaluate the gene expression stability quantitatively,four well-established algorithms -the comparative ΔCq method (Silveret al.2006),BestKeeper (Pfafflet al.2004),NormFinder(Andersenet al.2004),and GeNorm (Vandesompeleet al.2002) were introduced. The ΔCq method was used to compare the relative expression of ‘pairs of genes’ to identify the stability of reference genes (Silveret al.2006).BestKeeper determines the expression stability for each reference gene using the coefficient of variation (CV) of the relative expression levels after normalization(Pfafflet al.2004). Both NormFinder and GeNorm evaluate the variation of each gene pair along the sample series,and the former also takes intraand inter-group variation into consideration(Vandesompeleet al.2002;Andersenet al.2004). The expression stability coefficients of the candidates under various stress and hormone treatments are shown in Appendix B. Lower coefficients imply higher expression stability. The results demonstrated that the expression stability calculated by different alg orith ms differed.However,two traditionally used reference genes,Solyc04g009030 and Solyc07g066610,mostly demonstrated relatively poorer expression stability among the 22 genes under various treatments (Appendix B).

RefFinder,based on the rankings from each program above,calculated the overall final ranking of the references(Xieet al.2012). The number of reference genes recommended was determined with the help of GeNorm in qBase+Software.Based on the levels of variation in average reference gene stability with the sequential addition of each reference gene to the equation (for calculation of the normalization factor),the optimal number of reference targets in each or all treatments in this experiment was two(Fig.4). In different treatment conditions and in different tissues,the best reference combinations were different(Appendix B;Fig.4). For instance,in the control conditions,Solyc12g015870 and Solyc01g007860 were the best reference combination in leaves,and Solyc01g100870 and Solyc12g015870 the best in roots;when both leaves and roots were sampled and analyzed,Solyc11g039840 and Solyc01g100870 could be adopted as reference genes. When tomato seedlings were treated with PEG,NaCl,ABA,H2O2,or SA,the best reference combinations for gene expression analysis in leaves were Solyc10g078960+Solyc11g039840,Solyc04g078070+Solyc01g008000,Solyc01g008000+Solyc12g015870,Solyc01g008000+Solyc10g078960,and Solyc11g039840+Solyc06g009210,respectively;and the best reference combinations in roots were Solyc10g049850+Solyc01g100870,Solyc01g008000+Solyc10g078960,Solyc01g008000+Solyc06g009210,Solyc03g115490+Solyc12g015870,and Solyc01g111170+Solyc10g078960,respectively. When both leaves and roots were sampled for expression analysis,the best reference combinations were Solyc01g008000+Solyc03g115490 under water stress,Solyc01g008000+Solyc03g115490 under salt stress,Solyc01g008000+Solyc01g111170 under ABA treatment,Solyc12g100060+Solyc01g008000 under H2O2treatment,and Solyc01g008000+Solyc10g078960 under SA treatment. It is noticed that Solyc04g009030 and Solyc07g066610 were not recommended as preferred reference genes in any of the treatment conditions or tissues (Fig.4).

4.Discussion

Gene expression analysis has been an important part of biological research since it initially reflects how plants react to various internal and external signals (Lyndon and Francis 1992;Mele and Hake 2003). qPCR has been a commonly used technology for gene expression analysis due to its high sensibility and specificity (Gachonet al.2004;Zhanget al.2021). Normalizing qPCR raw data to one or more reference genes is the most commonly utilized approach (Huggettet al.2005),and thus introducing optimal reference genes is critical for the data normalization,while the use of unsuitable reference genes may dramatically change the interpretation of expression pattern (Gutierrezet al.2008). Mining transcriptome sequencing data is an approach to screening reference genes. However,relevant work is very limited. This study excavated 15 new candidate reference genes from publicly available transcriptional sequencing data in tomato;then,it evaluated the expression stability of these candidate reference genes together with seven traditionally used ones under different stress or hormone treatment conditions.

Quite a few previous researchers have evaluated the expression stability of reference genes in tomato in distinctive development processes or under various environmental stimuli (Table 1),and the candidates are mostly introduced from well-studied housekeeping genes.Based on the evaluation in the literature,seven reference genes were introduced in this study (Table 2;Appendix A).At least four genes of these,including Solyc01g088480,Solyc01g111170,Solyc04g009030,and Solyc07g066610 belonged to housekeeping genes,which are required for cell survival.

Identification of stably expressed genes from RNA-seq data provides a new approach to excavating reference genes for qPCR analysis,although most researchers are interested in differentially expressed genes. As RNAseq data can be seriously influenced by many factors,such as sample preparation and sequencing platform(Zhaoet al.2017),it is necessary to make integrative or comparative analyses among experiments to see if the same genes were stably expressed in different datasets.In this study,15 new candidate reference genes were screened out from a combined dataset of 529 samples,which included most public high-throughput transcriptional data of tomato as far as known. Among these new candidates,three genes encoded ADPribosylation factor,three encoded 60S ribosomal protein,two encoded ncRNA,two encoded enzymes,and others encoded translation factor,light-regulated or stressassociated protein (Table 2). Among these,the genes encoding ADP-ribosylation factor and 60S ribosomal protein have also been identified as ideal internal controls for expression analysis inCoffea arabica(Barsalobres-Cavallariet al.2009),castor bean (Cassolet al.2016),wheat (Longet al.2011),and potato (Tanget al.2017). In this study,nine of the 15 novel candidate reference genes were associated with or involved in basic biological processes,cellular components,or metabolism (Fig.2),indicating that they were housekeeping genes. When evaluating reference genes from transcriptomic sequencing data,previous scholars have also reported that some housekeeping genes were stably expressed. For instance,Zhouet al.(2017)found that MDP0000095375 (encoding casein kinase II subunit beta-4),MDP0000147424 (encoding a glutamate receptor interacting protein),MDP0000173025 (encoding membrane-related protein),and MDP0000233640(encoding ribosome biogenesis protein) were expressed stably in apple roots. Pomboet al.(2017) observed thatARD2(encoding acireductone dioxygenase) andVIN3(encoding vernalization insensitive protein 3) were the most stably expressed reference genes for studying tomato-Pseudomonaspathosystem. Tonget al.(2021)found thatGLY1andSKP1were the most stable genes across diverseAlliumtuberosumtissues. Therefore,our findings are consistent with previous studies that most reference genes were introduced from housekeeping genes (Expósito-Rodríguezet al.2008).

Given that high-throughput sequencing accuracy may be affected by the sequencing depth and coverage(Simset al.2014;Limaet al.2020),in this study,qPCR technology was used to verify the expression stability of the candidate reference genes under various stress conditions which are closely associated with plant growth and production. The conditions included water deficit,salt stress,oxidative stress,and chemical signals treatment like abscisic acid and salicylic acid. The results showed that the expression stability of the candidate reference genes differed when calculated by the four methods -comparative ΔCq method (Silveret al.2006),BestKeeper(Pfafflet al.2004),NormFinder (Andersenet al.2004),and GeNorm (Vandesompeleet al.2002). The differences may be associated with the calculation principles of these methods. Since the stability coefficients calculated by the above four methods were different,the comprehensive tool RefFinder was employed. RefFinder was based on the rankings from the above four programs,and it assigned an appropriate weight to an individual gene and calculated the geometric mean of their weights for the overall ranking (Xieet al.2012). The results demonstrated that the best reference gene combination differed between tissues and under different conditions (Fig.4). In addition,Solyc04g009030 and Solyc07g066610,two traditionally used genes,were not recommended as preferred reference genes in any of the treatment conditions or tissues (Fig.4),suggesting that the novel reference genes performed better than these two traditional ones. Similar results have also been reported previously. For instance,inGossypiumhirsutum,Smithaet al.(2019) found that Transmembrane 9 superfamily member 5 (TMN5) and Protein trichome birefringence-like 6 (TBL6) from publicly available RNA-seq dataset demonstrated more stable expression than the traditionally usedPP2A1(a protein phosphatase) andUBQ14(a polyubiquitin). Our study not only provides some new reference genes,but also provides preferred reference genes under various stress or hormone treatments,which may be useful for future research.

Some researchers have paid attention to the relationship between expression stability and gene length.Vinogradov (2004) and Eisenberg and Levanon (2013)found that the more stably expressed housekeeping genes were shorter than the tissue-specific expressed genes,which was driven by evolutionary selection for the economy;while Huanget al.(2019) did not find this relationship inCionasavignyi. In this study,no significant relationship between gene length and expression stability was observed,suggesting that evaluating reference genes based on gene length may not be a solid criterion.

5.Conclusion

Fifteen novel candidate reference genes for qPCR analysis were screened from tomato transcriptome sequencing data. All the novel candidate reference genes were expressed more stably than Solyc04g009030 and Solyc07g066610,two traditionally used genes,under different stress and hormone treatments. Our study suggests that excavating stably expressed genes from transcriptome sequencing data is a reliable approach to screening reference genes. This study provides some novel reference genes in tomato and the preferred reference genes under different environmental stimuli,which may be useful for future research.

Acknowledgements

This work was supported by the National Key Research and Development Program of China (2018YFD1000800),the National Natural Science Foundation of China(32072561) and the Natural Science Foundation of Henan,China (222300420282).

Declaration of competing interests

The authors declare that they have no conflict of interest.

Appendicesassociated with this paper are available on http://www.ChinaAgriSci.com/V2/En/appendix.htm