- Research article
Dissection of the style’s response to pollination using transcriptome profiling in self-compatible (Solanum pimpinellifolium) and self-incompatible (Solanum chilense) tomato species
BMC Plant Biologyvolume 15, Article number: 119 (2015)
Tomato (Solanum lycopersicum) self-compatibility (SC) is defined as self-pollen tubes that can penetrate their own stigma, elongate in the style and fertilize their own ovules. Self-incompatibility (SI) is defined as self-pollen tubes that are prevented from developing in the style. To determine the influence of gene expression on style self-pollination, a transcriptome-wide comparative analysis of SC and SI tomato unpollinated/pollinated styles was performed using RNA-sequencing (RNA-seq) data.
Transcriptome profiles of 24-h unpollination (UP) and self-pollination (P) styles from SC and SI tomato species were generated using high-throughput next generation sequencing. From the comparison of SC self-pollinated and unpollinated styles, 1341 differentially expressed genes (DEGs) were identified, of which 753 were downregulated and 588 were upregulated. From the comparison of SI self-pollinated and unpollinated styles, 804 DEGs were identified, of which 215 were downregulated and 589 were upregulated. Nine gene ontology (GO) terms were enriched significantly in SC and 78 GO terms were enriched significantly in SI. A total of 105 enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified in SC and 80 enriched KEGG pathways were identified in SI, among which “Cysteine and methionine metabolism pathway” and “Plant hormone signal transduction pathway” were significantly enriched in SI.
This study is the first global transcriptome-wide comparative analysis of SC and SI tomato unpollinated/pollinated styles. Advanced bioinformatic analysis of DEGs uncovered the pathways of “Cysteine and methionine metabolism” and “Plant hormone signal transduction”, which are likely to play important roles in the control of pollen tubes growth in SI species.
In flowering plants, the male organ of the flower is the stamen and the female organ of the flower is pistil. The stamen comprises an anther generating pollen grains and a filament supporting the anther. The pistil comprises the stigma, the style and the ovary. Pollination is a process of pollen-pistil interaction during which pollen adheres, hydrates, and germinates on the stigma, the pollen tube elongates on an active extracellular matrix in the style and finally transports male gametes (sperm cells) to the ovary, releasing it into ovules to complete fertilization . Mate selection is crucial to successful reproduction and species survival . Self-compatibility (SC) and self-incompatibility (SI) are the two predominant forms of mate selection. SC is defined as self-pollen that can penetrate its own pistil and fertilize its own ovules ; SI is where self-pollen is prevented from developing on the pistil .
Tomatoes (Solanum lycopersicum) are one of the most important vegetable crops in the world, and possess genetic diversities in fruit color, size, and mating system. In particular, the mating systems play key roles to control the reproductive habits between intra-/interspecies in tomatoes. Generally, color-fruited species such as Solanum lycopersicum, S. pimpinellifolium and S. neorickii are SC species, while some green-fruit species, such as S. habrochaites and S. chilense, are SI species . However, the growth of pollen tubes within styles differs between SI and SC species. Pollen growth is arrested at the middle style in SI species, but not in SC. Some models were proposed for growth behavior of pollen tubes within styles that are related to pollen factors such as F-box protein and pistil factor of RNase [5,6]; however, the mechanism controlling the growth of pollen tubes remains unclear in tomatoes.
The transcriptome is the sum of all the RNA transcription for specific cells in a certain functional condition, including mRNAs, non-coding RNAs (ncRNA) and small RNAs [7,8]. RNA-Seq is a deep-sequencing technology [7,9] that has many advantages compared with Serial Analysis of Gene Expression (SAGE) , Expressed Sequence Tag (EST) , cDNA-amplified fragment length polymorphism (AFLP) , DNA microarrays  and massively parallel signature sequencing (MPSS) . RNA-seq has already been widely used for transcriptome research in Miscanthus sinensis , tomato , Wolfiporia cocos , Hevea brasiliensis , Populus tomentosa , Lolium rigidum  and wheat . It has also been applied to study pollination in maize [22,23], and to study SC/SI in Citrus clementina , lemon  and Leymus chinensis . To understand what occurs after pollination in the styles of tomatoes of different mating types at the transcriptome level, we compared the transcription profiles differences between tomato SI and SC species. The results provide valuable information for understanding the growth behavior of pollen tubes within styles.
At present, research into tomato SC and SI has mainly concentrated on the S-RNase aspect, with no comprehensive transcriptome-level studies. Thus, to the best of our knowledge, this is the first study to perform comparative transcriptome analyses of SC and SI tomato unpollinated/pollinated styles using RNA-seq. The results of RNA-seq were analyzed by mapping, differential gene expression analysis, GO and pathway analysis. The results revealed comprehensive information concerning SI and SC, and provided clues to the molecular mechanisms of SI and SC.
Summary of RNA-seq datasets
SC unpollination/self-pollination (SCUP/SCP) and SI unpollination/self-pollination (SIUP/SIP) styles (total of 12 samples) were performed RNA-seq. The raw sequence data yielded approximately 3.0 gigabases (GB) per sample and more than 96% of the raw read pairs obtained had a quality score of ≥ Q20. Total raw read pairs among the 12 samples ranged from 15 to 18 million. By later removing reads containing adapters, reads containing poly-N and low-quality reads from the raw data, high-quality read pairs were obtained. The number of high-quality read pairs among the 12 samples ranged from 14 to 17 million (about 98% of the raw read pairs). Approximately 90% of the high-quality read pairs from the SC samples and 70% of the SI samples could be mapped to the tomato reference genome sequence. In addition, unmapped read pairs ranged from 1 to 5 million and multiple mapped read pairs ranged from about 0.30% to 0.50% of mapped read pairs among the 12 samples (Table 1).
Differential gene expression profiles of unpollinated (UP) and self-pollinated (P) styles in SC and SI, and hierarchical cluster analysis
To quantify the expression levels of the transcripts, HT-seq was used to count the read numbers mapped to each gene, based on the 34,726 genes of the tomato reference genome. These data were then normalized to reads per kilobase of exon region in a given gene per million mapped reads (RPKM) values, which were calculated based on the length of the gene and read count mapped to this gene. The RPKM values for each gene are listed in Additional file 1. To determine differential expression genes (DEGs) of UP and P styles in SC and SI, we screened for DEGs between UP and P styles in SC, and between UP and P styles in SI using the following criteria: Log2 fold-change (FC) > 1 or Log2FC < −1 and P-value < 0.05. We identified 1341 DEGs between UP and P styles in SC, and 804 DEGs between UP and P styles in SI (Additional file 2). Of these DEGs, 753 genes were downregulated and 588 genes were upregulated after self-pollination in SC; 215 genes were downregulated and 589 genes were upregulated after self-pollination in SI (Figure 1). We used hierarchical cluster analysis to compare the DEGs between UP and P styles in SC, between UP and P styles in SI, and the similarity of the expression patterns of the three biological replicates (Figure 1).
GO annotation of all DEGs in SCP vs. SCUP and SIP vs. SIUP
To identify the functions of thee DEGs, we performed gene ontology (GO) analysis. A total of 798 DEGs of SC comparing UP and P styles were assigned GO annotations and 525 DEGs of SI comparing UP and P styles were assigned GO annotations. GO has three ontologies: molecular function, cellular component and biological process. In many cases, one gene was annotated with multiple GO terms. The GO terms of 798 DEGs of SCP vs. SCUP styles were categorized into 42 main functional groups belonging to the three categories and the GO terms of 525 DEGs of SIP vs. SIUP styles were categorized into 41 main functional groups belonging to the three categories (Figure 2).
Comparative analysis of GO terms assigned to SCP vs. SCUP DEGs and those assigned to SIP vs. SIUP DEGs
To better understand the distribution of gene functions at the macro level, the GO function classification of the DEGs in SCP vs. SCUP styles and SIP vs. SIUP styles were analyzed using the WEGO online tool. The comparative analysis showed that DEGs in SCP vs. SCUP styles and SIP vs. SIUP styles shared broad similarities in the proportion of genes in the three main categories, but differences were detected in many subcategories (Figure 2). Most GO terms of DEGs in SCP vs. SCUP styles and SIP vs. SIUP styles were categorized into the same biological processes, cellular components and molecular functions. Most GO subcategories terms were detected in both of SCP vs. SCUP styles and SIP vs. SIUP styles; however, GO subcategory terms, including membrane-enclosed lumen, organelle part, molecular transducer, transcription regulator, biological regulation, developmental process, multicellular organismal process, pigmentation, reproduction, reproductive process and response to stimulus showed significant (P-value < 0.05) differences in counts between SCP vs. SCUP styles and SIP vs. SIUP styles. These results suggested that despite certain mechanisms of SC and SI appear to be conserved, the regulation mechanisms appear to be different between these two reproductive systems.
GO enrichment analysis of all DEGs in SCP vs. SCUP and SIP vs. SIUP
Significantly enriched GO terms were identified using singular enrichment analysis (SEA). The results showed that nine GO terms were significant in DEGs of SCP vs. SCUP based on a P-value < 0.05 and the false discovery rate (FDR) < 0.05 cutoffs (Figure 3A), which comprised two, three and four terms for the cellular components, molecular functions, biological processes categories, respectively. Seventy-eight GO terms were significant in DEGs of SIP vs. SIUP based on a P-value < 0.05 and the FDR < 0.05 cutoffs (Figure 3B, only 9), which comprised eight and 70 terms for the molecular functions and biological processes categories, respectively. The detailed results of the SCP vs. SCUP and SIP vs. SIUP Go enrichment analysis are presented in Additional file 3.
KEGG pathway mapping of all DEGs in SCP vs. SCUP and SIP vs. SIUP
To further investigate the influence of the DEGs on pathways, statistical pathway enrichment analysis of DEGs in SCP vs. SCUP and SIP vs. SIUP were performed based on KEGG database, using Fisher’s exact test. The DEGs of SCP vs. SCUP were enriched in 105 KEGG metabolic pathways and the DEGs of SIP vs. SIUP were enriched in 80 KEGG metabolic pathways (Additional file 4). The top ten KEGG metabolic pathways and three P-value < 0.05 metabolic pathways of the DEGs in SCP vs. SCUP are shown in Figure 4A. Among these 105 pathways of SCP vs. SCUP, those containing the greatest numbers of DEGs transcripts were “Metabolic pathways” (containing 111 DEGs) and “Biosynthesis of secondary metabolites” (containing 75 DEGs). Other GO terms associated with higher numbers of DEGs were “Starch and sucrose metabolism” (16 DEGs), “Plant hormone signal transduction” (16 DEGs), “Biosynthesis of amino acids” (15 DEGs), “Carbon metabolism” (15 DEGs), “Plant-pathogen interaction” (12 DEGs), “Phenylpropanoid biosynthesis” (11 DEGs), “Glycolysis/Gluconeogenesis” (nine DEGs), and “Amino sugar and nucleotide sugar metabolism” (eight DEGs); The pathways of “Biosynthesis of secondary metabolites”, “Biotin metabolism”, “Brassinosteroid biosynthesis” and “Degradation of aromatic compounds” had P-values < 0.05 (Figure 4A). For SIP vs. SIUP, of 13 KEGG metabolic pathways were identified. The top 11 KEGG metabolic pathways and two P-value < 0.05 metabolic pathways of DEGs in SIP vs. SIUP are shown in Figure 4B. Among the 80 pathways of SIP vs. SIUP, those containing the greatest numbers of DEGs were “Metabolic pathways” (69 DEGs), “Biosynthesis of secondary metabolites” (40 DEGs), “Plant hormone signal transduction” (22 DEGs), “Plant-pathogen interaction” (10 DEGs), “Starch and sucrose metabolism” (9 DEGs), “Biosynthesis of amino acids” (nine DEGs), “Phenylpropanoid biosynthesis” (nine DEGs), “Carbon metabolism” (eight DEGs), “Pentose and glucuronate interconversions” (eight DEGs), “Phenylalanine metabolism” (seven DEGs). The pathways of “Cysteine and methionine metabolism”, “Plant hormone signal transduction”, “Pentose and glucuronate interconversions”, “Flavonoid biosynthesis” and “Stilbenoid, diarylheptanoid and gingerol biosynthesis” all had P-values < 0.05 (Figure 4B). In addition, the pathways of “Cysteine and methionine metabolism” and “Plant hormone signal transduction” were significant pathways in DEGs of SIP vs. SIUP, based on a P-value < 0.05 and the FDR < 0.05 cutoffs (Figure 4B). The detailed results of the SIP vs. SIUP significant pathways enrichment analysis are presented in Figures 5 and 6.
“Cysteine and methionine metabolism” is the ethylene biosynthesis pathway, which was significantly enriched in the SIP vs. SIUP analysis. DEGs were enriched in the step of O-Acetyl-L-serine conversion to L-Cysteine, L-Homocysteine conversion to L-Methionine, L-Methionine conversion to S-adenosyl-L-methionine (AdoMet), AdoMet conversion to 1-aminocyclopropane-1-carboxylate (ACC) and ACC production ethylene (Figure 5). L-Methionine conversion to AdoMet was the first step of ethylene biosynthesis, AdoMet conversion to ACC was the rate-limiting step in ethylene biosynthesis and ACC production ethylene was the last steps for ethylene biosynthesis. Plant hormone signal transduction is very important to hormone-instigated biochemical changes during plant growth, development, and environmental information processing pathways, which were also significantly enriched in the SIP vs. SIUP comparison. DEGs were also enriched in Auxin signal transduction, Abscisic acid (ABA) signal transduction, Ethylene signal transduction, Jasmonic acid (JA) signal transduction and Salicylic acid (SA) signal transduction (Figure 6).
Significant pathways enrichment analysis showed that cysteine and methionine metabolism and plant hormone signal transduction were the most important pathways in SIP vs. SIUP comparison, and plant hormone signal transduction was the key biological event. All the plant hormone signaling pathways pointed to it and the significant pathway of “Cysteine and methionine metabolism” also (Figure 7). This evidence indicated that plant hormone signal transduction plays important roles in tomato SI.
RNA-seq is a powerful tool that can provide a global overview of gene expression at the transcriptome level. With reductions in sequencing costs and the advance of technologies, RNA-seq will become more accessible to researchers to identify and track the expression changes of all genes . The present study identified 1341 significant (P-value < 0.05) DEGs after comparing UP and P styles in SC and 804 significant (P-value < 0.05) DEGs in the comparison of UP and P styles in SI, using RNA-seq analysis. The total number of gene changes demonstrated that SC self-pollination and SI self-pollination are complex processes. This finding is consistent with other plant pollination studies. For example, 1025 differentially expressed genes were potentially involved in the pollination response and SI mechanisms in sheepgrass . In a comparison of pollinated and unpollinated stigmas with styles, 4785 DEGs were identified in SI lemon . These data demonstrate the complex nature of the transcriptome changes in SC self-pollination and SI self-pollination.
Pollination shares striking similarities with fungal infection in terms of biological responses and processes that result in cell death [27,28]. Our transcriptome GO enrichment analysis identified several significant GO terms involved in pathogen invasion responses, such as defense response to fungus, response to fungus, immune response, and immune system process in the SCP vs. SCUP comparison. This result is consistent with other plant pollination studies, such as in Arabidopsis [29,30] and rice . However, GO terms involved in stimuli and hormones were the most important of the 78 significant GO terms in the SIP vs. SIUP comparison.
Pollination leads to senescence of petunia corollas by inducing many hormonal, physiological, and molecular changes . Ethylene is a gaseous plant hormone with a wide range of effects on plant growth and development . Ethylene is synthesized from L-Methionine via the intermediates AdoMet and ACC (Figure 5) [34-36]. AdoMet is made from L-Methionine by the enzyme S-adenosylmethionine synthase (SAM), representing the first step of ethylene biosynthesis (Figure 5). 1-aminocyclopropane-1-carboxylate synthase (ACS) gene family members and 1-aminocyclopropane-1-carboxylate oxidase (ACO) gene family members are two important enzymes for ethylene biosynthesis. ACS catalyzes the conversion of AdoMet to ACC, which is the rate-limiting step in ethylene biosynthesis. ACO then catalyzes the conversion of ACC to ethylene (Figure 5) . After SI self-pollination, one SAM gene (S-adenosylmethionine synthase 2-like) (Solyc10g083970), five ACS gene family members (Solyc00g095760, Solyc08g081550, Solyc08g008100, Solyc08g081540, Solyc00g095860) and four ACO gene family members (Solyc02g036350, Solyc07g026650, Solyc07g049530, Solyc07g049550) were significantly upregulated, which indicated that SI self-pollination is associated with results in significant upregulation of ethylene biosynthesis related genes and ethylene production. It has been reported that ethylene biosynthesis is induced by pollination in petunias . After SC self-pollination, although the pathway of “Cysteine and methionine metabolism” was not a significant enrichment pathway in the SCP vs. SCP comparison, two ACS gene family members (Solyc08g081540, Solyc00g095860) and one ACO gene family member (Solyc07g049530) were significantly upregulated, which indicated that SC self-pollination results in some upregulation of ethylene biosynthesis of partly related genes. The above results suggest that SI self-pollination induces more ethylene production than SC self-pollination.
Plant hormone signal transduction is very important to hormone triggered biochemical changes . Plant hormone signal transduction plays an important role in pollination of petunias pollination; for example, RNA-seq revealed that plant hormone signal transduction-related KEGG pathways were enriched in petunia corollas when comparing pollinated and unpollinated samples . After SI self-pollination, plant hormone signal transduction-related KEGG pathways were significantly enriched in the SIP vs. SIUP comparison, but not after SC self-pollination (Figure 6). This result indicated that plant hormone signal transduction might play an important role in tomato SI. Plants recognize and transduce the ethylene signal via ethylene receptors (ETR)  in the ethylene signal transduction pathway (Figure 6) . We identified two ethylene receptors, LeETR6 (Solyc06g053710) and tETR (Solyc09g089610), which were significantly upregulated in the SIP vs. SIUP comparison, both of which mapped to the plant hormone signal transduction KEGG pathway. LeETR2 (Solyc07g056580) was the only ethylene receptor identified from the SCP vs. SCUP comparison, and significantly downregulated in P styles compared with UP styles. This protein also mapped to the plant hormone signal transduction KEGG pathway, which was not a significantly enriched pathway in the SCP vs. SCUP comparison. The above results indicated that SI self-pollination not only involves the induction of ethylene production, but also enhanced the perception ethylene. Although SC self-pollination may involve some enhancement of ethylene production, the ability to perceive ethylene was weakened by the significant downregulation of LeETR2. Plant responses to ethylene initiates with ethylene binding to ETRs and terminates in a transcription cascade of plant-specific transcription factors families, especially the ethylene-insensitive protein 3 (EIN3/EIL) and ethylene-responsive transcription factor (ERF). EIN3 protein is a key transcription factor for mediating the expression of ethylene-regulated genes and morphological responses. EIN3 interacts physically with the Ein3-binding f-box protein1/2 (EBF1/EBF2) and is ultimately and quickly degraded through a ubiquitin/proteasome pathway mediated by the SCF complex, which comprises a RING-box protein 1 (RBX1), Cullin 1 (Cul1), S-phase kinase-associated protein 1 (Skp1), F-box protein (F-box) [42,43]. We identified one EBF1/2 (Solyc07g008250) from the SC and two EBF1/2 (Solyc07g008250, Solyc12g009560) from the SI, both of which were significantly upregulated in P compared with UP styles. In addition, we also identified one Skp1 (Solyc01g111640) and one Cul1 (Solyc01g067120) from SI, which were significantly upregulated in P compared with UP styles. This result indicated that key transcription factor EIN3 was negatively regulated by targeting EIN3 it for degradation through the ubiquitin/proteasome pathway after SI self-pollination, but not in SC pollination.
A previous study demonstrated that auxin was significantly increased after compatible pollination and ethylene was strongly increased after incompatible pollination [44,46]. The last step of indole-3-acetic acid (IAA) biosynthesis is performed by aldehyde dehydrogenase. We identified one aldehyde dehydrogenase (aldehyde dehydrogenase family 2 member B4, Solyc08g068190) from SC that was significantly upregulated in P compared with UP styles and one aldehyde dehydrogenase (aldehyde dehydrogenase family 3 member H1-like, Solyc06g060250) from SI that was significantly downregulated in P compared with UP styles. This result is consistent with the results of the previous study. Auxin is likely to be directly or indirectly involved in pollen-pistil recognition and pollen tube elongation in Nicotiana  and might have an important role in the SI response in plants such as Theobroma cacao , Petunia hybrida  and Olea europaea . Auxins regulate plant growth and development by a complex signal transduction network , which was included in the significantly enriched KEGG pathways of plant hormone signal transduction KEGG in the SIP vs. SIUP comparison. Auxin influx carrier (AUX1 LAX family) is a polar auxin transporter in cells that is involved in attaining a hormone maximum (Figure 6) . After SC self-pollination, LAX2 protein (auxin influx carrier, AUX1 LAX family) (Solyc01g111310) was significantly downregulated. Auxins alter three major gene families: auxin/indole-3-acetic acid (Aux/IAA), GH3 and small auxin-up RNA (SAUR) to direct plant growth and development (Figure 6) [49,51]. Aux/IAA gene families: IAA1 (Solyc09g083280), IAA2 (Solyc06g084070), IAA3 (Solyc09g065850), IAA19 (Solyc03g120380), IAA22 (Solyc06g008580), IAA26 (Solyc03g121060), IAA35 (Solyc07Vg008020) and IAA36 (Solyc06g066020) were significantly upregulated in the SIP vs. SIUP comparison, and only IAA2 (Solyc06g084070), IAA29 (Solyc08g021820) and IAA 35 (Solyc07g008020) were significantly upregulated in the SCP vs. SCUP comparison. For the GH3 gene families, only one probable indole-3-acetic acid-amido synthetase GH3.1-like gene (Solyc02g092820) was significantly upregulated in the SCP vs. SCUP comparison. For the SAUR gene families, small auxin-up protein 58 (Solyc06g053260), auxin-induced protein 10A5-like (Solyc03g033590), uncharacterized LOC101249064 (Solyc03g124020) and uncharacterized LOC101254455 (Solyc12g009280) were significantly upregulated, and auxin-induced protein 15A-like (Solyc01g110570) and auxin-induced protein 10A5-like (Solyc01g110560) were significantly downregulated in the SIP vs. SIUP comparison. Only auxin-induced protein 15A-like (Solyc09g009980) and indole-3-acetic acid-induced protein ARG7-like (Solyc04g081250) were significantly upregulated in the SCP vs. SCUP comparison. These results indicated that although auxin was strongly increased after compatible pollination, because the auxin influx carrier (AUX1 LAX family) (Solyc01g111310) was significantly downregulated, fewer auxin-responsive genes showed altered expressions. During SC pollination, the auxin influx carrier (AUX1 LAX family) was not affected, resulting in many auxin-responsive genes showing altered expression after incompatible pollination. A previous study indicated that auxin influx carriers (AUX1 LAX family) were involved in auxin-ethylene interactions in Arabidopsis thaliana ; however, whether auxin influx carriers (AUX1 LAX family) are also involved in auxin-ethylene interactions in tomato SI is unknown.
Ethylene and JA, as well as ABA and auxin, have direct or indirect interactions , but the roles of JA and ABA in tomato pollination, especially in SI self-pollination, were unknown. ABA is a phytohormone that acts in seed dormancy, plant development and environmental stress. The carotenoid biosynthesis pathway is an ABA biosynthesis pathway (Figure 6) that was enriched in SC and SI. Endogenous ABA levels are regulated by both ABA biosynthesis and ABA catabolism: xanthoxin dehydrogenase is a key enzyme for ABA biosynthesis and ABA 8′-hydroxylase is a key enzyme for ABA catabolism [53,54]. Xanthoxin dehydrogenase (Solyc12g056600) was significantly upregulated in SC not in SI and ABA 8′-hydroxylase 1-like (CYP707A2, Solyc08g005610) was significantly upregulated in SI but not in SC, which indicated that endogenous ABA levels increased in SC and decreased in SI styles during pollination. Pyrabactin resistance/pyrabactin resistance-like (PYR/PYL) family is an ABA receptor that is very important to ABA recognition and signaling [55,56]. We identified two genes of the PYR/PYL family: ABA receptor PYL8-like (Solyc03g007310) from SI and ABA receptor PYL6-like (Solyc06g050500) from SC. PYL8-like was significantly downregulated in SI and PYL6-like was significantly upregulated in SC styles during pollination, which indicated that the ability to perceive ABA was weakened in SI and enhance in SC. A previous study showed that PYR/PYLs are negative regulatory receptors, whereby ABA binds to PYR/PYLs, which in turn binds to type 2C protein phosphatases (PP2Cs) to inhibit PP2Cs. SNF1-related protein kinase subfamily 2 (SnRK2) is located downstream of PP2Cs and is negatively regulated by PP2Cs (Figure 6). SnRK2 (serine/threonine-protein kinase SAPK3-like, Solyc08g077780) was upregulated in SI (in which PP2Cs are not inhibited) and an SnRK2 (serine/threonine-protein kinase SAPK7-like, Solyc05g056550) was downregulated in SC, wherePP2Cs are inhibited. In addition, SnRK2s can phosphorylate b-ZIP transcription factors, which bind to the ABA-responsive element to activate ABA-responsive genes. Phosphorylated b-ZIP transcription factors are important to active ABA-responsive genes . One b-ZIP transcription factor (Solyc10g076920) was significantly downregulated in the SCP vs. SCUP comparison, but not in the SIP vs. SIUP comparison. This indicated that ABA might have important regulatory roles in SI. Jasmonates are phytohormones that are essential for plant development and survival, and can induce jasmonate ZIM-domain proteins (JAZs) to be degraded through the ubiquitin/proteasome pathway, mediated by the SCFCOI1 complex. In addition, JAZs negatively regulate MYC2, which is a key jasmonate responses transcriptional activator . We identified a JAZ (jasmonate ZIM-domain protein 1, Solyc12g009220) and a transcription factor MYC2 (Solyc08g076930), both of which were both significantly upregulated in the SIP vs. SIUP comparison. The TGA family comprises key transcription factors of the salicylic acid (SA)-mediated signal transduction pathway . After SI self-pollination, TGA family transcription factor (Solyc10g080410) was significantly upregulated.
This is the first global transcriptome-wide comparative analysis of styles from SC and SI tomatoes using a high-throughput RNA-seq. The enriched GO term analysis of the identified DEGs showed that nine GO terms were significantly enriched in the SCP vs. SCUP comparison and 78 GO terms were significantly enriched in the SIP vs. SIUP comparison. The ethylene biosynthesis pathway of the cysteine and methionine metabolism pathway and the plant hormone signal transduction pathway play an important role in tomato SI. Further GO and KEGG analyses showed that SI self-pollination induced more ethylene production and catabolism of ABA, and SC self-pollination induced more auxin production and ABA biosynthesis. Moreover, the phytohormones ethylene, auxin and ABA play important roles by plant hormone signal transduction in tomato SI.
Tomato seeds of S. chilense (LA0130, SI) and S. pimpinellifolium (LA1585, SC) were obtained from the Charles Rick Tomato Genetics Resource Center (UC, Davis http://tgrc.ucdavis.edu/index.aspx). The seeds were germinated in peat pellets and seedlings with three to four leaves were grown on medium containing the perlite: peat (1:1) under a thermoperiod of 26/20°C (day/night) in a greenhouse. Plants were supplied with a commercial fertilizer every week. During flowering, 24 h UP and P styles (containing stigmas) (Additional file 5) were collected from S. chilense (LA0130) (SIUP/SIP) and S. pimpinellifolium (LA1585) (SCUP/SCP), respectively, and immediately frozen in liquid nitrogen and stored at −80°C for RNA extraction. Three biological replicates of each sample were collected and used for RNA extraction.
RNA extraction and deep sequencing
Total RNA was extracted from each sample using an RNAprep pure Plant Kit (Tiangen, Beijing, China), according to the manufacturer’s protocol. The RNA concentration of each sample was measured using a NanoDrop 2000 (Thermo Scientific, Waltham, MA, USA). The RNA quality was assessed using an Agilent2200 (Agilent Technologies, Santa Clara, CA, USA).
The sequencing library for each RNA sample was prepared using a TruseqTM RNA sample prep Kit (Illumina, San Diego, CA, USA), following the manufacturer’s protocol. Briefly, mRNA was purified using poly-T oligo-attached magnetic beads (Invitrogen,Carlsbad, CA, USA) from 5 μg total RNA. The mRNA was fragmented, and the RNA fragments were reverse transcribed and amplified to double-stranded cDNA. Index adapters were then ligated to the cDNA according to the protocol of the TruseqTM RNA sample prep Kit (Illumina). The library was quantified using a TBS-380 mini-fluorometer (Picogreen, Cohasset, MA, USA). The clustering of the index-coded samples was performed on a cBot Cluster Generation System, using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and a sequence length of 2*101 bp paired-end reads were generated.
Filtering raw reads and mapping
The raw reads were pass-filtered using the Trimmomatic tool  and then used for mapping. The reference tomato genome and gene model annotation files were downloaded from the genome website (http://solgenomics.net/) directly. The paired-end clean reads were aligned to the reference tomato genome using Tophat  and the mapped reads were counted with using HT-seq .
Identification of DEGs
Gene expression levels were estimated as RPKM . Differential expression analysis of SCUP/SCP groups and SIUP/SIP groups was performed using the DESeq R package (1.10.1), which provides statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. After statistical analysis, the DEGs were identified using significance analysis by t-tests, with a P-value < 0.05 and at least two-fold changes (either up- or downregulation) being considered significant.
The blast2go  program was used to obtain GO annotations for all identified genes. GO functional classification was performed using the WEGO online tool  to gain an understanding of the distribution of gene functions at the macro level. GO is the key functional classification of NCBI, which was applied to analyze the functions of the DEGs [66,67]. GO enrichment analysis of DEGs was implemented using SEA , in which Fisher’s exact test and a χ2 test were used to classify the GO categories; the FDR was calculated to correct the P-value [69,70]. P-values for the GOs of all the DEGs were computed. The significant GO terms were defined as having a P-value < 0.05 and an FDR < 0.05.
KEGG is a database resource for understanding high-level functions and utilities of biological systems, such as cells, organisms and ecosystems, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-through put experimental technologies (http://www.genome.jp/kegg/). KEGG pathway analysis was used to identify the significant pathways involving the DCEGs [71-73]. Fisher’s exact test and a χ2 test were used to identify significant pathways (P-value < 0.05 and FDR < 0.05) [74-76]. We used the KEGG Orthology Based Annotation System (KOBAS) software to test the statistical enrichment of DEGs in KEGG pathways.
Availability of supporting data
The data sets supporting the results of this article are available in the Gene Expression Omnibus repository under accession no GSE67654 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67654) .
Serial Analysis of Gene Expression
Expressed sequence tag
Massively parallel signature sequencing
Reads per kilobase of exon region in a given gene per million mapped reads
Differentially expressed genes
Singular enrichment analysis
False discovery rate
Kyoto encyclopedia of genes and genomes
S-adenosyl methionine synthase
Ethylene-insensitive protein 3
Ethylene-responsive transcription factor
Ein3-binding f-box protein1/2
- AUX1 LAX:
Auxin influx carrier
Small auxin-up RNA
RING-box protein 1
S-phase kinase-associated protein 1
Pyrabactin resistance/pyrabactin resistance-like
Type 2C protein phosphatases
SNF1-related protein kinase subfamily 2
Jasmonate ZIM-domain protein
KEGG Orthology Based Annotation System
Lord EM, Russell SD. The mechanisms of pollination and fertilization in plants. Annu Rev Cell Dev Biol. 2002;18:81–105.
Samuel MA, Tang W, Jamshed M, Northey J, Patel D, Smith D, et al. Proteomic analysis of Brassica stigmatic proteins following the self-incompatibility reaction reveals a role for microtubule dynamics during pollen responses. Mol Cell Proteomics. 2011. Doi: 10.1074/mcp.M111.011338.
Kear PJ, McClure B. How did flowering plants learn to avoid blind date mistakes? Self-incompatibility in plants and comparisons with non self rejection in the immune response. Adv Exp Med Biol. 2012;738:108–23.
Li W, Chetelat RT. A pollen factor linking inter- and intraspecific pollen rejection in tomato. Science. 2010;330(6012):1827–30.
Qiao H, Wang F, Zhao L, Zhou J, Lai Z, Zhang Y, et al. The F-box protein AhSLF-S2 controls the pollen function of S-RNase-based self-incompatibility. Plant Cell. 2004a;16:2307–22.
Franklin-Tong VE, Franklin FC. Gametophytic self-incompatibility inhibits pollen tube growth using different mechanisms. Trends Plant Sci. 2003;8:598–605.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.
Costa V, Angelini C, De Feis I, Ciccodicola A. Uncovering the complexity of transcriptomes with RNA-Seq. J Biomed Biotechnol. 2010;2010:853916.
Haas BJ, Zody MC. Advancing RNA-Seq analysis. Nat Biotechnol. 2010;28:421–3.
Veleulescu VE, Zhang L, Zhou W, Vogelstein J, Basrai MA, Bassett Jr DE, et al. Characterization of the yeast transcriptome. Cell. 1997;88:243–51.
Adams MD, Kelley JM, Gocayne JD, Dubnick M, Polymeropoulos MH, Xiao H, et al. Complementary DNA sequencing: expressed sequence tags and human genome project. Science. 1991;252:1651–6.
Brugmans B, del Fernandez CA, Bachem CW, van Os H, van Eck HJ, Visser RG. A novel method for the construction of genome wide transcriptome maps. Plant J. 2002;31:211–22.
Schena M, Shalon D, Davis RW, Brown PO. Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science. 1995;270:467–70.
Brenner S, Johnson M, Bridgham J, Golda G, Lloyd DH, Johnson D, et al. Gene expression analysis by massively parallel signature sequencing (MPSS) on microbead arrays. Nat Biotechnol. 2000;18(6):630–4.
Swaminathan K, Chae WB, Mitros T, Varala K, Xie L, Barling A, et al. A framework genetic map for Miscanthus sinensis from RNAseq-based markers shows recent tetraploidy. BMC Genomics. 2012. doi: 10.1186/1471-2164-13-142.
Koeniga D, Jiménez-Gómez JM, Kimura S, Fulop D, Chitwood DH, Headland LR, et al. Comparative transcriptomics reveals patterns of selection in domesticated and wild tomato. Proc Natl Acad Sci U S A. 2013;110:E2655–62.
Shu S, Chen B, Zhou M, Zhao X, Xia H, Wang M. De novo sequencing and transcriptome analysis of Wolfiporia cocos to reveal genes related to biosynthesis of triterpenoids. PLoS One. 2013;8(8), e71350.
Salgado LR, Koop DM, Pinheiro DG, Rivallan R, Le Guen V, Nicolás MF, et al. De novo transcriptome analysis of Hevea brasiliensis tissues by RNA-seq and screening for molecular markers. BMC Genomics. 2014;15:236. doi:10.1186/1471-2164-15-236.
Liao W, Ji L, Wang J, Chen Z, Ye M, Ma H, et al. Identification of glutathione S-transferase genes responding to pathogen infestation in populus tomentosa. Funct Integr Genomics. 2014;14:517–29.
Gaines TA, Lorentz L, Figge A, Herrmann J, Maiwald F, Ott MC, et al. RNA-Seq transcriptome analysis to identify genes involved in metabolism-based diclofop resistance in Lolium rigidum. Plant J. 2014;78:865–76.
Li A, Liu D, Wu J, Zhao X, Hao M, Geng S, et al. mRNA and small RNA transcriptomes reveal insights into dynamic homoeolog regulation of allopolyploid heterosis in nascent hexaploid wheat. Plant Cell. 2014;26:1878–900.
Xu XH, Chen H, Sang YL, Wang F, Ma JP, Gao XQ, et al. Identification of genes specifically or preferentially expressed in maize silk reveals similarity and diversity in transcript abundance of different dry stigmas. BMC Genomics. 2012;13:294. doi: 10.1186/1471-2164-13-294.
Li G, Wang D, Yang R, Logan K, Chen H, Zhang S, et al. Temporal patterns of gene expression in developing maize endosperm identified through transcriptome sequencing. Proc Natl Acad Sci U S A. 2014;111:7582–7.
Caruso M, Merelo P, Distefano G, Malfa SL, Lo Piero AR, Tadeo FR, et al. Comparative transcriptome analysis of stylar canal cells identifies novel candidate genes implicated in the self-incompatibility response of citrus clementina. BMC Plant Biol. 2012;12:20. doi:10.1186/1471-2229-12-20.
Zhang S, Ding F, He X, Luo C, Huang G, Hu Y. Characterization of the ‘Xiangshui’ lemon transcriptome by de novo assembly to discover genes associated with self-incompatibility. Mol Genet Genomics. 2015;290(1):365–75.
Zhou Q, Jia J, Huang X, Yan X, Cheng L, Chen S, et al. The large-scale investigation of gene expression in Leymus chinensis stigmas provides a valuable resource for understanding the mechanisms of poaceae self-incompatibility. BMC Genomics. 2014. doi: 10.1186/1471-2164-15-399.
van Doorn WG, Woltering EJ. Physiology and molecular biology of petal senescence. J Exp Bot. 2008;59(3):453–80.
Tintor N, Ross A, Kanehara K, Yamada K, Fan L, Kemmerling B, et al. Layered pattern receptor signaling via ethylene and endogenous elicitor peptides during Arabidopsis immunity to bacterial infection. Proc Natl Acad Sci U S A. 2013;110(15):6211–6.
Tung CW, Dwyer KG, Nasrallah ME, Nasrallah JB. Genome-wide identification of genes expressed in arabidopsis pistils specifically along the path of pollen tube growth. Plant Physiol. 2005;138:977–89.
Swanson R, Clark T, Preuss D. Expression profiling of arabidopsis stigma tissue identifies stigma-specific genes. Sex Plant Reprod. 2005;18:163–71.
Li M, Xu W, Yang W, Kong Z, Xue Y. Genome-wide gene expression profiling reveals conserved and novel molecular functions of the stigma in rice (Oryza sativa L.). Plant Physiol. 2007;144:1797–812.
Broderick SR, Wijeratne S, Wijeratn AJ, Chapin LJ, Meulia T, Jones ML. RNA-sequencing reveals early, dynamic transcriptome changes in the corollas of pollinated petunias. BMC Plant Biol. 2014. doi:10.1186/s12870-014-0307-2.
Wilkinson JQ, Lanahan MB, Clark DG, Bleecker AB, Chang C, Meyerowitz EM, et al. A dominant mutant receptor from Arabidopsis confers ethylene insensitivity in heterologous plants. Nat Biotechnol. 1997;15(5):444–7.
Yang SF, Hoffman NE. Ethylene biosynthesis and its regulation in higher plants. Annu Rev Plant Physiol. 1984;35:155–89.
Kende H. Enzymes of ethylene biosynthesis. Plant Physiol. 1989;91(1):1–4.
Bleecker AB, Kende H. Ethylene: a gaseous signal molecule in plants. Annu Rev Cell Dev Biol. 2000;16:1–18.
Chae HS, Kieber JJ. Eto Brute? Role of ACS turnover in regulating ethylene biosynthesis. Trends Plant Sci. 2005;10(6):291–6.
Singh A, Evensen KB, Kao TH. Ethylene synthesis and floral senescence following compatible and incompatible pollinations in Petunia inflata. Plant Physiol. 1992;99(1):38–45.
Bowler C, Chua NH. Emerging themes of plant signal transduction. Plant Cell. 1994;6:1529–41.
Schaller GE, Bleecker AB. Ethylene-binding sites generated in yeast expressing the Arabidopsis EtRI K7gene. Science. 1995;270:1809–11.
Sakai H, Hua J, Chen QG, Chang C, Medrano LJ, Bleecker AB, et al. ETR2 is an ETR1-like gene involved in ethylene signaling in Arabidopsis. Proc Natl Acad Sci U S A. 1998;95:5812–7.
Potuschak T, Lechner E, Parmentier Y, Yanagisawa S, Grava S, Koncz C, et al. EIN3-dependent regulation of plant ethylene hormone signaling by two arabidopsis F box proteins: EBF1 and EBF2. Cell. 2003;115(6):679–89.
Guo H, Ecker JR. Plant responses to ethylene gas are mediated by SCF (EBF1/EBF2)-dependent proteolysis of EIN3 transcription factor. Cell. 2003;115(6):667–77.
Holden MJ, Marty JA, Singh-Cundy A. Pollination-induced ethylene promotes the early phase of pollen tube growth in Petunia inflata. J Plant Physiol. 2003;160(3):261–9.
Chen D, Zhao J. Free IAA in stigmas and styles during pollen germination and pollen tube growth of Nicotiana tabacum. Physiol Plant. 2008;134:202–15.
Hasenstein KH, Zavada MS. Auxin modification of the incompatibility response in Theobroma cacao. Physiol Plant. 2001;112:113–8.
Kovaleva L, Zakharova E. Hormonal status of the pollen-pistil system at the progamic phase of fertilization after compatible and incompatible pollination in Petunia hybrida L. Sex Plant Reprod. 2003;16:191–6.
Solfanelli C, Bartolini S, Vitagliano C, Lorenzi R. Immunolocalization and quantification of IAA after self- and free pollination in Olea europaea L. Sci Hortic. 2006;110:345–51.
Khan S, Stone JM. Arabidopsis thaliana GH3.9 influences primary root growth. Planta. 2007;226(1):21–34.
Grebe M, Friml J, Swarup R, Ljung K, Sandberg G, Terlou M, et al. Cell polarity signaling in Arabidopsis involves a BFA-sensitive auxin influx pathway. Curr Biol. 2002;12(4):329–34.
Jain M, Kaur N, Tyagi AK, Khurana JP. The auxin-responsive GH3 gene family in rice (Oryza sativa). Funct Integr Genomics. 2006;6(1):36–46.
Vandenbussche F, Petrásek J, Zádníková P, Hoyerová K, Pesek B, Raz V, et al. The auxin influx carriers AUX1 and LAX3 are involved in auxin-ethylene interactions during apical hook development in Arabidopsis thaliana seedlings. Development. 2010;137(4):597–606.
Gonzalez-Guzman M, Apostolova N, Belles JM, Barrero JM, Piqueras P, Ponce MR, et al. The short-chain alcohol dehydrogenase ABA2 catalyzes the conversion of xanthoxin to abscisic aldehyde. Plant Cell. 2002;14:1833–46.
Okamoto M, Kuwahara A, Seo M, Kushiro T, Asami T, Hirai N, et al. CYP707A1 and CYP707A2, which encode abscisic acid 8′-hydroxylases, are indispensable for proper control of seed dormancy and germination in Arabidopsis. Plant Physiol. 2006;141:97–107.
Nishimura N, Hitomi K, Arvai AS, Rambo RP, Hitomi C, Cutler SR, et al. Structural mechanism of abscisic acid binding and signaling by dimeric PYR1. Science. 2009;326:1373–9.
Park SY, Fung P, Nishimura N, Jensen DR, Fujii H, Zhao Y, et al. Abscisic acid inhibits type 2C protein phosphatases via the PYR/PYL family of START proteins. Science. 2009. doi: 10.1126/science.1173041.
Fujii H, Zhu JK. Arabidopsis mutant deficient in 3 abscisic acid-activated protein kinases reveals critical roles in growth, reproduction, and stress. Proc Natl Acad Sci U S A. 2009;106:8380–5.
Chini A, Fonseca S, Fernandez G, Adie B, Chico JM, Lorenzo O, et al. The JAZ family of repressors is the missing link in jasmonate signalling. Nature. 2007;448:666–71.
Zhou JM, Trifa Y, Silva H, Pontier D, Lam E, Shah J, et al. NPR1 differentially interacts with members of the TGA/OBF family of transcription factors that bind an element of the PR-1 gene required for induction by salicylic acid. Mol Plant Microbe Interact. 2000;13:191–202.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.
Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31(2):166–9.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNASeq. Nat Methods. 2008;5:621–8.
Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.
Ye J, Fang L, Zheng HK, Zhang Y, Chen J, Zhang ZJ, et al. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34:W293–7.
Gene Ontology Consortium. The Gene Ontology (GO) project in 2006. Nucleic Acids Res. 2006;34(Database issue):D322–6.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.
Du Z, Zhou X, Ling Y, Zhang Z, Su Z. AgriGO: a GO analysis tool kit for the agricultural community. Nucleic Acids Res. 2010;38:W64–70.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Sta Soc Ser B. 1995;57(1):289–300.
Pawitan Y, Michiels S, Koscielny S, Gusnanto A, Ploner A. False discovery rate, sensitivity and sample size for microarray studies. Bioinformatics. 2005;21:3017–24.
Wang K, Singh D, Zeng Z, Coleman SJ, Huang Y, Savich GL, et al. MapSplice: accurate mapping of RNA-seq reads for splice junction discovery. Nucleic Acids Res. 2010;38(18), e178.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.
Joshi-Tope G, Gillespie M, Vastrik I, D’Eustachio P, Schmidt E, de Bono B, et al. Reactome: a knowledgebase of biological pathways. Nucleic Acids Res. 2005;33:D428–32.
Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The KEGG resource for deciphering the genome. Nucleic Acids Res. 2004;32 (Database issue):D277–80.
Yi M, Horton JD, Cohen JC, Hobbs HH, Stephens RM. WholePathwayScope: a comprehensive pathway-based analysis tool for high-throughput data. BMC Bioinformatics. 2006;7:30. doi: 10.1186/1471-2105-7-30.
Draghici S, Khatri P, Tarca AL, Amin K, Done A, Voichita C, et al. A systems biology approach for pathway level analysis. Genome Res. 2007;17:1537–45.
Zhao P, Zhang L, Zhao L. Dissection of the style’s response to pollination using transcriptome profiling in self-compatible (Solanum pimpinellifolium) and self-incompatible (Solanum chilense) tomato species. Gene Expression Omnibus. http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE67654.
We thank the Charles Rick Tomato Genetics Resource Center at the University of California Davis for supplying the tomato seeds for this study. The research was supported by the Key Technology Research and Development Program of Shanghai Science and Technology Committee (No. 13391901202 and No. 14JC1403400), the National Natural Science Foundation of China (No. 31071810) and the China National ‘863’ High-Tech Program (No. 2011AA100607).
The authors declare that they have no competing interests.
ZP participated in the experimental design, collected the material, performed the RNA extraction, participated in the bioinformatics analyses and wrote the manuscript. ZL participated in the bioinformatics analyses. ZL designed and wrote this manuscript. All authors read and approved the final manuscript.
RNA-seq of data of all counts for SI and SC compared with self-pollinated and unpollinated styles.
List of differentially expressed genes for SI and SC compared with self-pollinated and unpollinated styles.
GO analysis of differentially expressed genes for SI and SC compared with self-pollinated and unpollinated styles.
Pathway analysis differentially expressed genes for SI and SC compared with self-pollinated and unpollinated styles.
The structure of the tomato pistil. Red lines show the cutting position of a style containing a stigma.