Comparative transcriptome analysis provides insights into dwarfism in cherry tomato (Solanum lycopersicum var. cerasiforme)

Tomato, which can be eaten as a vegetable or fruit, is one of the most popular and nutritionally important crops around the world. Although most plants of the cherry tomato cultivar ‘Minichal’ have a normal phenotype, some plants have a stunted phenotype with reduced plant height, leaf size, and fruit size, as well as altered leaf and fruit shape. To investigate the molecular mechanisms underlying these differences, we generated RNA-seq libraries from pooled leaf samples of 10 normal (N) and 10 stunted (S) plants. Using the Illumina sequencing platform, we obtained a total of 115.45 million high-quality clean reads assembled into 35,216 genes and 35,216 transcripts. A total of 661 genes were differentially expressed between N and S plants. Of these, 420 differentially expressed genes (DEGs) were up-regulated, and 221 DEGs were down-regulated. The RNA-seq data were validated using quantitative reverse-transcription PCR. Enrichment analysis of DEGs using the Kyoto Encyclopedia of Genes and Genomes (KEGG) showed that the enriched pathways were involved in steroid biosynthesis, homologous recombination, and mismatch repair. Among these, three genes related to steroid biosynthesis, including 3BETAHSD/D2, DIM and DWF5 were down-regulated in S compared to N. Of these, DIM and DWF5 are known to be involved in brassinosteroid biosynthesis. Our results thus provide a useful insight into dwarfism in cherry tomato, and offer a platform for evaluating related species.


Introduction
Cultivated tomato (Solanum lycopersicum L.) is nutritionally rich, economically important, and widely grown around the world. It is ranked as the second most-consumed vegetable after the potato [1]. Tomato can be consumed fresh or in processed food items such as ketchup, paste, juice, pizza sauce, and soup. The ripe tomato fruit is abundant in lycopene, a red carotenoid pigment that has antioxidant properties, which help to protect against heart diseases, and lung and prostate cancer [2][3][4][5]. It also contains other carotenoids, including beta-carotene, neurosporene, lutein, and zeaxanthin, which support the human immune system [6]. PLOS  development, were used (Fig 1). These lines were grown in a glasshouse at the Department of Horticulture, Sunchon National University, Suncheon, Republic of Korea. Young leaves were sampled from 10 individual plants for each phenotypic category (N and S). Leaves were pooled and frozen in liquid nitrogen before storing them at -80˚C until required.

RNA extraction and library construction for transcriptome analysis
Total RNA was isolated from 100 mg finely powdered leaf tissues using the RNeasy Mini Kit (Qiagen, USA). The quantity and integrity were checked with a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, Delaware, USA) and an Agilent 2100 BioAnalyzer (Agilent Technologies, Palo Alto, CA, USA). Two RNA-seq libraries were constructed by Theragen Bio Institute (Suwon, South Korea) using the TruSeq RNA Library Prep Kit (Illumina Inc.) and RNA samples with a RIN (RNA integrity number) greater than 7. RNA sequencing was performed using an Illumina HiSeq 2000 platform (Illumina Inc.). RNA sequencing data were analyzed according to the method described by Trapnell et al. [29].

Quantification of expression patterns and differentially expressed genes
Reference genome and gene model annotation files for tomato (Solanum lycopersicum) were retrieved from the Ensembl database (https://plants.ensembl.org/). Clean reads were mapped to the reference genome using TopHat (v.2.1.1; http://ccb.jhu.edu/). Assembled genes were searched against the Swiss-Prot database and Gene Ontology (GO) categories. Gene expression patterns and differential expression were determined using Cufflinks (v.2.0.1; http:// cufflinks.cbcb.umd.edu/), as previously reported by Trapnell et al. [29]. The expression level was normalized by the number of fragments per kilobase of exon per million mapped reads (FPKM). Differentially expressed genes (DEGs) were detected using DEGseq [30] with an adjusted p < 0.005 and q < 0.05. All DEGs were subjected to GO analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using WebGestalt [31] and DAVID (https://david.ncifcrf.gov/).

Validation of RNA-seq data by qRT-PCR
The expression patterns of eight genes were selected for further validation by quantitative reverse transcription PCR (qRT-PCR). cDNA was synthesized from 2 μg of high-quality total RNA using SuperScript III (Invitrogen, Gaithersburg, MD). The qRT-PCR reaction was carried out using 50 ng cDNA and a gene-specific primer (S1 Table)

Overview of RNA sequencing
Leaves at similar stages of growth were collected and pooled from normal (N) and stunted (S) tomato plants for RNA isolation. Their transcriptomes were profiled using the Illumina sequencing platform. We obtained 117.995 million paired-end raw reads (Table 1). Subsequently, adapters, low-quality reads, and ambiguous reads were removed (Fig 2A). A final total of 115.450 million high-quality clean reads were obtained (54.755 and 60.695 million for N and S, respectively) ( Table 1). A total of 97.6% reads from N plants, and 98.0% reads from S plants, were mapped to the S. lycopersicum reference genome (Ensembl). On average, 91.6% and 91.5% reads, respectively, uniquely mapped to the reference database. High-quality clean reads were assembled into 35,216 transcripts and 35,216 genes (Table 1). Among the annotated genes, 72% and 71% had 90-100% coverage in N and S tomato libraries, respectively (Fig 2B and 2C), indicating that the distributions of reads were similar between tomato libraries. Both the samples had Q20 scores (indicating Phred-like quality) greater than 96%, indicating the high quality of the RNA sequencing. These high quality transcriptomic data from N and S plants therefore provide a basis for identifying the key genes involved in tomato dwarfism.

Differentially expressed genes between normal and stunted tomato pools
A total of 661 differentially expressed genes (DEGs) between N and S tomato plants were identified using the R package DEGseq (S2 Table) [30]. Of these DEGs, 420 genes were up-regulated, and 241 genes were down-regulated in S versus N (Fig 3). However, 32 DEGs were found to be expressed in S only, and 108 in N only (S2 Table). The distribution of up-regulated and down-regulated genes is shown using a volcano plot (Fig 4).
To obtain insight into the biological significance of identified DEGs, GO enrichment analysis was performed using the Gene Ontology database (http://www.geneontology.org/). Enriched GO terms for genes that were up-regulated and down-regulated between N and S tomato plants are shown in Fig 6. GO enrichment analysis revealed that 'catalytic activity' and 'metabolic process' were the most often enriched GO terms for both up-regulated and downregulated genes (S3 Table).

Validation of RNA-seq data
To test the reliability of RNA-seq results, qRT-PCR was used to measure the expression of eight genes with the same RNA samples used for RNA-seq. Relative expression of the tested genes was consistent with the RNA-seq data (Fig 10), confirming the efficiency and accuracy of the RNA-seq experiments.

Discussion
Deep sequencing-based RNA-seq technology has made it possible to rapidly analyze large genomic datasets and quantify transcriptomes [33]. This high-throughput, next-generation sequencing technology has become a powerful tool for analyzing transcriptomes, and has been successfully used for both human and plant transcriptomes [34]. Global gene expression patterns can be determined using RNA-seq in samples of tissues at different developmental stages, with contrasting characteristics, or in response to different environmental stimuli [33,35,36].
In this study, we observed a stunted phenotype of the cherry tomato cv. 'Minichal', which is characterized by alterations in plant height, leaf size/shape, and fruit size/shape compared to the normal phenotype. We used RNA-seq to profile the transcriptomes of normal (N) and stunted (S) cherry tomatoes of this cultivar. We obtained almost 115.450 million high-quality clean reads, which were assembled into 35,216 transcripts (Table 1). Our results identified 661 DEGs between the pooled RNA of N and S tomato plants (S2 Table). Subsequently, GO enrichment revealed that 'metabolic process' and 'catalytic activity' were the most enriched GO terms for both up-regulated and down-regulated genes between N and S (S3 Table).
To obtain further insight into the biological functions of these DEGs, GO functional annotation and KEGG pathway enrichment analysis were performed using the DAVID tool. Functional annotation clustering of DEGs revealed that the most enriched GO terms were associated with the sterol biosynthesis process (GO:0016126; enrichment score 2.37) ( Table 2). KEGG pathway enrichment also indicated that 'steroid biosynthesis' and 'homologous recombination' were the most enriched pathways (Table 3). These results clearly suggest that genes related to steroid biosynthesis might be involved in dwarfism in S tomatoes.
Several studies have been conducted on plant dwarfism. Dwarfism is sometimes advantageous; for example in cereal crops-specifically rice and wheat, where lodging decreases crop productivity [37]. However, in tomato, dwarfism is deleterious because it reduces both quality and productivity. Plant dwarfism results from many genetic defects, mostly associated with hormone biosynthesis and perception [12]. Functional analysis of several genes related to dwarfism has previously been reported, including genes related to BR and GA biosynthesis and perception in different plant species [22,28,[38][39][40][41].
BRs play significant roles in plant growth and development, and are biosynthesized via multiple parallel pathways starting with the precursor campesterol [15,24,42,43]. Defects in the BR biosynthesis/signaling cause dwarfism in plants [13]. For example, in Arabidopsis, dwarf5 (dwf5) mutants have a mutation in the gene for the enzyme 7-dehydrocholesterol reductase, which disrupts the sterol Δ 7 reduction step and leads to dwarfism [25]. Likewise in Arabidopsis,

Fig 6. GO enrichment of up-regulated (A) and down-regulated (B) genes in leaves of stunted (S) and normal (N) cherry tomato plants (cv. 'Minichal').
The P-value was corrected as P < 0.5.
https://doi.org/10.1371/journal.pone.0208770.g006 Transcriptome analysis provides insights into dwarfism in cherry tomato dwarfism occurs in dwf4 mutants, which have a mutation in the gene encoding steroid 22α hydroxylase (CYP90B1), which is involved in 22α-hydroxylation of the BR pathway [24]. In this study, three DEGs, 3beta-hydroxysteroid-dehydrogenase (3BETAHSD/D2, Solyc02g081730.2), 7-dehydrocholesterol reductase (DWF5, Solyc06g074090.2), and delta(24)sterol reductase (DIM, Solyc02g069490.2)-all related to the steroid hormone biosynthesis pathway-were down-regulated in S compared to N plants (Figs 7 and 10). We also checked the expression patterns of these three genes in leaf, inflorescence, and fruit tissues of the cherry tomato cv. 'Minichal' (Fig 11). The result indicated that the expression of these steroid biosynthesis genes was higher in N than S plants. In N, expression was highest in leaf and lowest in fruit, while in S, expression was similar in leaves and inflorescences, but was drastically reduced in fruits.
The protein interaction network of three steroid biosynthesis-related genes, 3beta-hydroxysteroid-dehydrogenase (3BETAHSD/D2), 7-dehydrocholesterol reductase (DWF5), and delta (24)-sterol reductase (DIM), highlighted their possible contribution to dwarfism of tomato plants (Fig 12). In S cherry tomatoes, these genes might be involved in dwarfism by directly or indirectly affecting steroid biosynthesis. In apple plants (Malus × domestica), colchicineinduced autotetraploid plants showed dwarfism, with decreased levels of indole-3-acetic acid (IAA) and BR compared to diploid plants. Furthermore, digital gene expression analysis of these apple plants revealed that DEGs between them were mostly related to IAA and BR biosynthesis pathways [44]. In Arabidopsis, a biosynthetic defect in dwf1, which encodes delta (24)-sterol reductase, resulted in dwarfism with reduced levels of BR synthesis compared to the wild type [45]. A similar dwf1 dwarf mutant has been reported in pea [46]. The dwf5 mutant, which is defective in BR biosynthesis, also showed a dwarf phenotype in Arabidopsis [25]. In rice, the dwarf mutant ebisu dwarf (d2) is deficient in BR biosynthesis and caused dwarfism, but exogenous application of BL restored the normal phenotype [39]. In tomato, the BR biosynthesis-defective mutant Dwarf (D), which harbors a mutation in cytochrome P450 (P450), exhibits dwarfism, while complementation 35S::D lines restore the normal phenotype [12,40]. Similar dwarf mutant dumpy (dpy) resulted from the mutation of mutation in the C-23 steroid hydroxylase (cpd) gene has also been reported in tomato by Kaka et al. [47]. However, our reported genes (3BETAHSD/D2, DWF5 and DIM) for dwarfism of cherry tomato are different from those previously reported mutants like D and dpy. Up-regulation of cytokinin dehydrogenase 3 (CKX3) in S tomatoes led to a higher rate of cytokinin degradation in these plants. Reid et al. [48] reported that the cytokinin content was negatively regulated by the activity of CKX3 in the root of Lotus japonicus ckx3 mutants. The auxin signaling genes AIR12 (Solyc09g056390.1), and AXX15 (Solyc04g053000.1), but not IAA14, were down-regulated in N (Fig 8). This suggests that auxin signaling genes might be affected, leading to defective plant development. A similar result has been reported in tetraploid apple [44]. Previous studies have revealed that GA has an effect on plant growth and development. For example, exogenous treatment of dwarf pea and dwarf maize seedlings with GA3 enhanced longitudinal growth rates [49]. However, we found no DEGs related to GA biosynthesis in this study.  The up-regulation of two 1-aminocyclopropane-1-carboxylic acid oxidase genes, ACO1 and ACO3, which are involved in the final step of ethylene biosynthesis, suggests higher levels of Transcriptome analysis provides insights into dwarfism in cherry tomato Transcriptome analysis provides insights into dwarfism in cherry tomato ethylene production in S tomatoes, which might affect plant growth and development. Ethylene overproduction has been shown to inhibit plant growth in Arabidopsis [50,51].
The 'short vegetative phase' (SVP) group of MADS-box genes, such as OsMADS22, OsMADS47, and OsMADS55, have been shown to act as negative regulators of BR responses in rice [52]. The double and triple RNAi plants (OsMADS22-OsMADS55 and OsMADS22-OsMADS47-OsMADS55, respectively) showed reduced stem elongation. Unexpectedly, in this study, we also found that the expression of the MADS-box genes SVP (Solyc04g076280.2) and AGL19 (Solyc08g080100.2) was down-regulated, and AGL36 (Solyc01g103550.1) and SEPAL-LATA2 (Solyc02g089200.2) were up-regulated in S compared to N plants (Fig 9). Overexpression of OsMADS1 causes dwarfism in rice via irregular activation of BR and GA synthesis pathways [53]. Transcriptome analysis provides insights into dwarfism in cherry tomato Kim et al. [54] demonstrated that BR controls stomatal development by activating mitogenactivated protein kinase kinase kinase (MAPKKK) in Arabidopsis. Likewise, we found that MAPKK was up-regulated in S compared to N tomatoes (Fig 9).

Conclusions
We conducted comparative transcriptome analysis using normal and stunted plants of the cherry tomato cv. 'Minichal'. DEGs related to steroid biosynthesis may be involved in dwarfism in this tomato cultivar. To best of our knowledge, this is the first comparative transcriptome analysis for plant dwarfism in tomato. Our results provide insight into the molecular mechanism of dwarfism and lay the foundation for future studies in related species.
Supporting information S1