The ys mutant of tomato was generated in the background of the inbred line TTD302A through EMS mutagenesis using Saito’s method (Saito et al., 2011), and stabilized via six generations of selfing prior to this study. The seeds of ys mutant and WT were germinated on wet filter paper in a Petri dish at 28°C in dark for 3 days. Then the resulting seedlings were grown in a growth chamber under a 16 h/8 h (light/dark) photoperiod with 25/16°C temperatures, respectively. Upon four true-leaf stage, plants were transferred to a greenhouse in the experimental field of the Northwest A&F University. Pest control and water management were carried out according to standard practices.

The fresh stigmas at the anthesis stage were collected from WT and ys mutant at the same time on the same day, and then immediately used to pigment measurement. Total chlorophylls and carotenoids were extracted and quantified by spectrophotometric method as described previously (Lichtenthaler, 1987; Bou-Torrent et al., 2015).

Flavonoids and other polyphenols extraction and analyses were carried out according to Bino et al. (2005) with some modifications. The 100 mg frozen tomato stigmas powder was extracted with 80% methanol for 12 h at 4°C followed by 30 min sonication. The mixtures were then centrifuged at 12,000 rpm for 10 min at 4°C and the supernatants were filtered (0.2 μm). The extraction was repeated once by the above method. The samples were analyzed and identified using a Waters Alliance 2695 HT HPLC system coupled with a Q-TOF Mass Spectrometer (Waters-Micromass, Milford, MA, United States). The 20 μL of sample extract was injected and separated using a Luna C18 reverse-phase analytical column (3 μm, 150 mm × 2.1 mm; Phenomenex) at 40°C and a gradient of 5–50% acetonitrile in 0.1% formic acid at a flow rate of 0.3 mL min^-1. The eluted compounds were detected at 200–600 nm using a 996 PDA detector (Waters, Milford, MA, United States), followed by the mass spectrometer equipped with an electrospray ionization (ESI) source. The conditions for LC-MS runs were as follows: desolvation temperature was 250°C with a nitrogen gas flow rate of 10 L min^-1 and capillary voltage was 3 kV; source temperature was 120°C; cone voltage was 35 eV with 1 L min^-1 gas flow and collision energy was 5 eV in positive ion mode or 10 eV in ion mode. Ions in the m/z range 100–1500 were detected using a scan time of 0.9 s and an interscan delay of 0.1 s. The flavonoid compounds and other polyphenols were identified using retention times obtained by authentic standards and mass calculated for the ion (M + H)⁺ (Supplementary Table S1), and quantified by comparing the area of each individual peak with the standard curves obtained from the pure compounds. All measurements were repeated with three independent biological samples, and each sample was assayed in triplicate.

The stigmas at the anthesis stage from WT and ys mutant were harvested and freeze-dried in liquid nitrogen. The lignin content was determined using the acetyl bromide method as described previously with some modifications (Yan et al., 2012). The samples were ground into powder, then 10 mg of powder was rinsed four times with 95% ethanol and twice with the mixture of 100% ethanol and n-hexane (1:2 in volume). The precipitate was collected, dried at 60°C and then suspended in 2 mL of 25% acetyl bromide in glacial acetic acid. After incubation at 70°C for 30 min, 0.9 mL 2 M NaOH was added, followed by 2 mL glacial acetic acid and 0.1 mL 7.5 M hydroxylamine hydrochloride. The mixture was centrifuged at 4,500 rpm for 5 min and the supernatant was collected, then the absorbance was measured at 280 nm. The lignin content was expressed as Absorption 280 on a fresh weight basis. All measurements were performed with three biological samples, and each sample was assayed in triplicate.

Young stigmas of about 1 mm in length were collected from WT and ys mutant at the same time on the same day. Samples were immediately frozen in liquid nitrogen and stored at -80°C for RNA-Seq analyses. Total RNA was isolated using the RNA extraction kit (Promega, United States). RNA quality was checked by RNase-free agarose gel electrophoresis to avoid possible degradation and contamination, and then verified using Agilent 2100 Bio-analyzer (Agilent Technologies, Santa Clara, CA, United States). Next, Poly (A) mRNA was isolated using oligo-dT beads (Qiagen, Germany), and then broken into short fragments by adding fragmentation buffer. First-strand cDNA was synthesized using random hexamer-primed reverse transcription, followed by the synthesis of the second-strand cDNA using RNase H and DNA polymerase I. The cDNA fragments were purified using a QIA quick PCR extraction kit, and then washed with EB buffer for end reparation poly (A) addition and ligated to sequencing adapters. Following agarose gel electrophoresis and extraction of cDNA from gels, the cDNA fragments were purified and enriched by PCR to construct the final cDNA library, which was then sequenced on the Illumina sequencing platform (Illumina HiSeq^TM 2500) using the paired-end technology. Three biological replicates were performed for each line from WT and ys mutant, thus six DGE libraries were generated and sequenced.

Raw reads were filtered to remove low quality sequences (there were more than 50% bases with quality lower than 20 in one sequence), reads with more than 5% N bases (bases unknown) and reads containing adaptor sequences through the Perl program. Then the clean reads were mapped to the tomato reference genome using TopHat2 (Consortium, 2012; Kim D. et al., 2013), allowing up to one mismatch. Unigenes mapped by at least one read, in at least one sample, were identified for further analysis. The DEGs were identified using the R package edgeR (Robinson et al., 2010). The expression level of each unigene was calculated and normalized to FPKM. The FDR was used to determine the threshold of the P-value in multiple tests. In our study, the FDR < 0.05 and fold change > 2 were used as significance cut-offs of the gene expression differences.

Sequencing data were deposited to the Short Read Archive (SRA) database at the National Center for Biotechnology Information (NCBI) under the accession number SRP080654.

Further, the DEGs were used for GO and KEGG enrichment analyses according to Zhang et al. (2013). GO terms with corrected P-value < 0.05 and KEGG pathways with P-value < 0.05 were considered significantly enriched by differential expressed genes.

Quantitative real-time PCR assays were performed using the independent tomato stigmas in the same developmental stage as those used for DGE analysis. Total RNA was isolated using the RNA extraction kit (Promega, United States), and cDNA was synthesized using MultiScribe reverse transcriptase (Applied Biosystems, United States). qRT-PCR was carried out using SYBR^® Premix Ex Taq^TM from TaKaRa (China) on an Applied Biosystems 7500 real-time PCR system (Applied Biosystems, United States). The tomato EF-1α gene was used as reference gene to normalize the expression data. Each qRT-PCR experiment was repeated with three biological samples, and each sample was assayed in triplicate. The relative expression of genes was calculated using the 2^-ΔΔCt method and standard error was calculated between three biological replicates. The gene specific primers for qRT-PCR are listed in Supplementary Table S2.

While the WT tomato plant was characterized by typical green stigma, the ys mutant, generated through EMS mutagenesis, showed unusual “yellow stigma” phenotype (Figures 2A–C). Notably, this stigma color difference between WT and ys mutant was apparent even in the very early stage of stigma development (stigma length around 1 mm) (Figure 2D).

To explore important pigments responsible for the “yellow stigma” phenotype, we quantified the contents of chlorophylls, carotenoids and flavonoids in the stigmas of WT and ys mutant during anthesis stage. To our surprise, no significant changes were noticed in both total chlorophyll and carotenoid levels between WT and ys mutant (Figure 3A and Supplementary Table S3). However, remarkable differences were found for some polyphenol contents in the flavonoid biosynthetic pathway. For instance, the yellow stigma accumulated 4.3- and 8.9-fold increased p-coumaric acid and naringenin chalcone, respectively. However, the levels of kaempferol-3-rutinoside and quercetin-3-rutinoside (rutin) in stigmas were not changed between WT and ys mutant (Figure 3B and Supplementary Table S3). It is worth mentioning that the p-coumaric acid is a biosynthetic precursor of naringenin chalcone (Figure 1) (Petrussa et al., 2013; Saito et al., 2013; Falcone Ferreyra et al., 2015). Therefore, it is highly plausible that the higher accumulation of p-coumaric acid in the ys mutant led to an increased naringenin chalcone level. These results suggested that increased accumulation of yellow-colored flavonoid naringenin chalcone may be responsible for the yellow stigma phenotype in the ys mutant.

To identify genes and gene networks that are involved in the formation of yellow stigma in tomato, we performed genome-wide expression analyses to compare the transcriptome profiles of the stigmas between WT and ys mutant through the DGE approach (Eveland et al., 2010). Given that the phenotypic changes occurred in the early stage (Figure 2D), we chose young stigmas (length around 1 mm) for RNA-Seq analyses. In this project, 42.3–50.8 million raw reads from each DGE library were generated. After removal of adapter sequences and low-quality reads, we obtained 41.8–50.1 million high-quality clean reads with a total of 5.2–6.3 billion nucleotides. Among these clean reads, the percentage of Q20 (base quality more than 20) and GC was 94.4–94.8% and 42.8–43.0%, respectively (Supplementary Figure S1 and Table S4). Further, we clustered the clean reads into unique reads, which were mapped to the tomato genome using TopHat2 (Consortium, 2012; Kim D. et al., 2013). In the six DGE libraries, about 91.2–92.9% of clean reads from RNA-Seq data were mapped uniquely to the tomato genome (Supplementary Table S4).

Based on deep sequencing, 23,978 genes were detected in all libraries (Supplementary Table S5). We used the R package edgeR to identify the DEGs (Robinson et al., 2010). The expression level of each gene was normalized to FPKM. Using FDR < 0.05 and fold change > 2 as the significance cut-offs, we obtained 507 DEGs, in which 84 genes were significantly up-regulated and 423 genes were dramatically down-regulated in the stigmas of ys mutant as compared to those of WT plants (Supplementary Table S6). To validate the RNA-Seq data, we performed qRT-PCR assays using the independent tomato stigmas in the same developmental stage as those used for DGE analysis. Sixteen DEGs, including 6 up-regulated genes and 10 down-regulated genes, were randomly chosen for qRT-PCR analysis. As shown in Table 1, the qRT-PCR data were in full agreement with the RNA-Seq data in terms of relative fold change in the expression of these 16 genes between WT and ys mutant (Pearson correlation coefficient 0.969, P = 6.3E - 10), suggesting that the RNA-Seq results were highly reliable.

In order to understand the expression profiles and potential functions of DEGs identified by DGE, GO term enrichment analysis (corrected P-value < 0.05) was performed. The DEGs were classified into cellular component, biological process and molecular function categories (Supplementary Table S7). As shown in Figure 4, there was only one group “extracellular region” (P = 8.8E - 09) in the cellular component category. The GO terms of catabolic and metabolic processes were dominantly enriched in the biological process group. In the molecular function category, the top three significantly enriched GO terms were “oxidoreductase activity” (P = 1.7E - 11), “heme binding” (P = 1.7E - 08) and “iron ion binding” (P = 3.0E - 08). Interestingly, these three groups all shared 23 genes encoding cytochrome P450s (CYP450), a large family of enzymes that play a crucial role in the biosynthesis of flavonoids that are responsible for flower colors (Tanaka, 2006; Tanaka and Brugliera, 2013), and 10 genes encoding PODs (Supplementary Table S7), however, most of those genes were down-regulated in the stigmas of ys mutant compared to WT (Table 2). For example, the expression levels of two CYP450 genes, Solyc10g078220 and Solyc10g078230, and two POD genes, Solyc07g017880 and Solyc02g084790, decreased 139.8-, 130.8-, 54.9- and 28.2-fold, respectively, in the ys mutant compared with those in the WT, and qRT-PCR verification revealed the same expression pattern (Table 1). These data indicated that the cytochrome P450s and PODs may play important role in the formation of yellow stigma in tomato.

To further identify the biological pathways that are responsible for the formation of yellow stigma in tomato, we mapped the detected DEGs to reference canonical pathways in the KEGG (Kanehisa et al., 2004), and then compared these with the whole transcriptome background to search for genes involved in the metabolic or signal transduction pathways that were significantly enriched. Using the P-value < 0.05 as the significance cut-off, 29 DEGs, including 3 up-regulated and 26 down-regulated genes, were assigned to 8 KEGG pathways (Supplementary Table S8). Among them, the “Phenylpropanoid biosynthesis” (P = 1.4E - 03), “Biosynthesis of secondary metabolites” (P = 6.5E - 03) and “Biosynthesis of unsaturated fatty acids” (P = 1.7E - 02) were the top three significantly enriched KEGG pathways (Figure 5). However, the “Flavonoid biosynthesis” pathway (P = 4.4E - 02), which was putatively associated with the yellow stigma phenotype in the ys mutant (Figure 3B), was also enriched, despite the P-value showed higher. Notably, two CYP450 genes (Solyc10g078220 and Solyc10g078230) and six POD genes (Solyc07g017880, Solyc02g084790, Solyc10g076240, Solyc03g006700, Solyc02g079500, and Solyc04g071890), whose expressions were dramatically down-regulated in the stigma of ys mutant (Table 2), were mapped into the “Phenylpropanoid biosynthesis” and “Biosynthesis of secondary metabolites” (Supplementary Table S8). Moreover, the two CYP450 genes were also involved in the “Flavonoid biosynthesis” pathway (Supplementary Table S8). In these pathways, the two CYP450 genes encoded the C3H which catalyze the p-coumaric acid, p-coumaroyl-CoA and p-coumaraldehyde to caffeic acid, caffeoyl-CoA and caffeyl aldehyde, respectively (Figure 6 and Supplementary Figure S2) (Xue et al., 2014; Fornaléa et al., 2015), so, we named these two genes as SlC3H1 (Solyc10g078220) and SlC3H2 (Solyc10g078230). In addition, the six POD genes were involved in the last step of biosynthesis of lignin, an insoluble cell wall-associated polymer (Figure 6) (Quiroga et al., 2000; Higuchi, 2006).

It is worth mentioning that the p-coumaric acid acts as the common precursor for both flavonoid and lignin biosyntheses, and C3H and POD are two key enzymes in lignin biosynthetic pathway (Higuchi, 2006; van der Rest et al., 2006; Petrussa et al., 2013; Renault et al., 2017). Given the repression of two SlC3Hs and six PODs in yellow stigma of the ys mutant (Table 2), we speculated that the higher accumulation of p-coumaric acid in the ys mutant (Figure 3B) may result from the reduced lignin biosynthesis. To test this hypothesis, we compared the levels of caffeic acid, the downstream production of p-coumaric acid catalyzed by C3H and also a biosynthetic precursor of lignin (Figure 6) (Higuchi, 2006; Renault et al., 2017), and lignin in the stigmas at anthesis stage between WT and ys mutant. As expected, the caffeic acid and lignin contents were 7.6- and 8.5-fold decreased in the ys mutant compared to those of WT, respectively (Figure 7 and Supplementary Table S3), suggesting that the reduced expressions of two SlC3Hs and six PODs might inhibit the biosynthesis of lignin, leading to a change in the direction of the metabolic flow of p-coumaric acid from lignin to flavonoid biosynthetic pathway, which eventually resulted in a higher naringenin chalcone level (Figure 3B). These observations indicated that the two SlC3Hs and six PODs may function as negative regulators in naringenin chalcone biosynthesis, which could play a positive role in yellow stigma formation in the ys mutant.