Transcriptomic and proteomic analysis reveals mechanisms of low pollen-pistil compatibility during water lily cross breeding
BMC Plant Biology volume 19, Article number: 542 (2019)
In water lily (Nymphaea) hybrid breeding, breeders often encounter non-viable seeds, which make it difficult to transfer desired or targeted genes of different Nymphaea germplasm. We found that pre-fertilization barriers were the main factor in the failure of the hybridization of Nymphaea. The mechanism of low compatibility between the pollen and stigma remains unclear; therefore, we studied the differences of stigma transcripts and proteomes at 0, 2, and 6 h after pollination (HAP). Moreover, some regulatory genes and functional proteins that may cause low pollen-pistil compatibility in Nymphaea were identified.
RNA-seq was performed for three comparisons (2 vs 0 HAP, 6 vs 2 HAP, 6 vs 0 HAP), and the number of differentially expressed genes (DEGs) was 8789 (4680 were up-regulated), 6401 (3020 were up-regulated), and 11,284 (6148 were up-regulated), respectively. Using label-free analysis, 75 (2 vs 0 HAP) proteins (43 increased and 32 decreased), nine (6 vs 2 HAP) proteins (three increased and six decreased), and 90 (6 vs 0 HAP) proteins (52 increased and 38 decreased) were defined as differentially expressed proteins (DEPs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that the DEGs and DEPs were mainly involved in cell wall organization or biogenesis, S-adenosylmethionine (SAM) metabolism, hydrogen peroxide decomposition and metabolism, reactive oxygen species (ROS) metabolism, secondary metabolism, secondary metabolite biosynthesis, and phenylpropanoid biosynthesis.
Our transcriptomic and proteomic analysis highlighted specific genes, incuding those in ROS metabolism, biosynthesis of flavonoids, SAM metabolism, cell wall organization or biogenesis and phenylpropanoid biosynthesis that warrant further study in investigations of the pollen-stigma interaction of water lily. This study strengthens our understanding of the mechanism of low pollen-pistil compatibility in Nymphaea at the molecular level, and provides a theoretical basis for overcoming the pre-fertilization barriers in Nymphaea in the future.
Water lilies (Nymphaea) are important flowering plants that are distributed worldwide from the tropics to temperate regions . With the rapid improvement of China’s economy and the overall quality of life, the demand is increasing for new water lily hybrids with different characteristics. Therefore, it is necessary to breed new water lily hybrids with excellent ornamental characteristics. However, in the breeding of water lily hybrids, breeders often encounter non-viable seeds, which makes it difficult to transfer the desired or targeted genes of various water lily germplasm [2, 3]. For example, many breeders have hoped to transfer the colorful flowers of tropical water lilies to hardy water lilies through crossbreeding; however, viable hardy water lily varieties with blue flowers have not yet been developed . In recent years, we have carried out interspecific hybridization between the female N. ‘Peter Slocum’ and male N. micrantha for three consecutive years, aiming at transferring the color gene of male parent to female parent. However, we also did not obtain seeds, so we carried out a thorough and systematic study from the aspect of plant reproductive biology, and found that the main reason for the failure of the hybrid combination was the low compatibility between pollen and stigma before fertilization . Therefore, in this study, an interspecific cross between the female ‘Peter Slocum’ and male N. micrantha was performed. Our aim was to further reveal the reasons of low compatibility between pollen and stigma at the molecular level on the basis of previous studies.
Low compatibility between the pollen and stigma is a common issue that negatively impacts the efficiency of plant breeding and the yield of seeds or fruit [5, 6]. Therefore, over the past several decades, many researchers have conducted studies to investigate factors that cause low compatibility between the pollen and stigma [7,8,9,10]. However, the mechanisms underlying low compatibility between the pollen and stigma in Nymphaea remain poorly understood.
With the development of molecular biology technology, the use of transcriptome and proteomics technology may provide a new way to find the genes and proteins related to low compatibility between pollen and stigma [11,12,13]. In particular, transcriptome sequencing is a useful method for identifying novel transcripts and analyzing gene expression [14, 15]. Transcriptomic and proteomic analyses have been extensively applied to many plant species, but limited transcriptome and proteome data exists regarding pre-fertilization barriers in water lily [16, 17]. To understand the mechanism of low pollen-pistil compatibility in water lily at the genomic level, Illumina paired-end sequencing and a label-free analysis of the stigma after pollination were conducted. This comprehensive analysis of the transcriptome and proteome may substantially improve the overall understanding of the potential molecular mechanisms involved in low pollen-pistil compatibility in water lily and pave the way for further analyses. This study aimed to provide important molecular data supporting a deep understanding of low compatibility between the pollen and stigma in water lily and also provides an important clue to overcome hybridization barriers.
Pollen germination on stigmas after pollination
Previous studies showed that pollen began to germinate at 2 HAP, and abnormal growth of pollen tubes was observed at 6 HAP . The number of pollen tubes did not change at 6 and 12 HAP, indicating that the related genes in stigma had been expressed and new proteins were synthesized, while no new proteins were synthesized after 6 h. For this reason, the stigmas of 0, 2, and 6 HAP were collected.
In the ‘Peter Slocum’ × N. micrantha cross, no pollen tubes penetrated stigmatic tissue between 2 and 6 HAP. In addition, the accumulation of wax between the stigma and the surface of the pollen grains was commonly observed (Fig. 1a). In the self-pollinated ‘Peter Slocum’, the stigma and the surface of the pollen grains showed no wax at 6 HAP (Fig. 1b).
Overview of the transcriptomic analysis and proteomics analysis
Using Fragments Per Kilobase Million (FPKM), we explored the gene expression levels in the stigmas 0, 2, and 6 HAP. In three comparisons (2 vs 0 HAP, 6 vs 2 HAP, 6 vs 0 HAP), the number of DEGs was 8789 (4680 were up-regulated), 6401 (3020 were up-regulated), and 11,284 (6148 were up-regulated), respectively. Further details of the DEGs are presented in Additional file 1. Using the label-free analysis, a total of 3176 proteins were identified within a false discover rate (FDR) of 1% (Additional file 2). Following the proteins (2 vs 0 HAP; 43 increased and 32 decreased), nine proteins (6 vs 2 HAP; three increased and six decreased), and 90 proteins (6 vs 0 HAP; 52 increased and 38 decreased) were defined as DEPs (Additional file 3).
Comparison analysis of transcriptome and proteome data
To identify robust pathways that were supported by both datasets, we integrated DEGs and DEPs to find the corresponding genes and proteins, and the results are listed in Additional file 4. Overlaps between DEPs and DEGs are shown by Venn diagrams in Fig. 2. Specifically, there were considerable non-overlaps between DEPs and DEGs, probably due to the relatively low sensitivity of proteome detection. For instance, among the 234 differentially regulated proteins in the 2 vs 0 HAP comparison, only 96 genes and their corresponding proteins were differentially expressed (Fig 2a). Similarly, for the 6 vs 2 HAP and 6 vs 0 HAP comparisons, 12 and 127 of the DEPs, respectively, were correlated to the corresponding DEGs (Fig. 2b, c).
All expression data associated with protein level and transcription level were analyzed and Spearman correlation coefficient was calculated. Globally, the correlation coefficients of all quantitative proteins and their corresponding genes at 2 vs 0 HAP, 6 vs 2 HAP, and 6 vs 0 HAP were 0.2236, 0.02 and 0.123, respectively (Fig. 3a). However, there is a high correlation between the DEGs and their corresponding DEPs (r = 0.8178, 0.4, and 0.6985, respectively; Fig. 3b). The correlation between proteins and their corresponding mRNAs with the same or opposite trend was analysed, and the comparative group 2 vs 0 HAP and 6 vs 0 HAP had higher positive or negative correlations (Fig. 3c and Fig. 3d). However, we found poor correlations between proteins at 6 vs 2 HAP and their corresponding mRNAs with the same or opposite trend (Fig. 3c and Fig. 3d). Among the Correlation-DEGs-DEPs (cor-DEGs-DEPs) genes, 39 (2 vs 0 HAP), three (6 vs 2 HAP), and 41 (6 vs 0 HAP) genes had the same trend, while five (2 vs 0 HAP), one (6 vs 2 HAP), and seven (6 vs 0 HAP) genes had the opposite trend (Additional file 5). Thus, we suggest that some of these cor-DEGs-DEPs genes might play important roles in causing low pollen-pistil compatibility during water lily breeding.
Cluster analysis of expression patterns in the cor-DEGs-DEPs genes
Cluster analysis of the DEPs and their corresponding DEGs can visually show their expression patterns, and the results are shown in Fig. 4. Cluster analysis showed that 44 (2 vs 0 HAP), four (6 vs 2 HAP), and 48 (6 vs 0 HAP) DEPs were correlated with the change of mRNA abundance, and 39 (2 vs 0 HAP), three (6 vs 2 HAP), and 41 (6 vs 0 HAP) DEPs were matched with corresponding DEGS. However, 5 (2 vs 0 HAP), one (6 vs 2 HAP) and seven (6 vs 0 HAP) DEPS were opposite to their mRNA expression pattern. In total, 25 of the cor-DEGs-DEPs genes at 2 vs 0 HAP and 6 vs 0 HAP showed the same expression pattern; thus, we infer that these genes may cause low pollen-pistil compatibility in interspecific hybridization of water lily.
GO and pathway enrichment analysis of the cor-DEGs-DEPs genes
After obtaining the expression information of the cor-DEGs-DEPs genes at 2 vs 0 HAP, 6 vs 2 HAP, and 6 vs 0 HAP, GO functional annotation analysis of these genes was carried out (Table 1; Table 2). The results showed that 27 (2 vs 0 HAP), 0 (6 vs 2 HAP), 19 (6 vs 0 HAP) GO terms were highly enriched at both mRNA and protein levels (Fig. 5; Additional file 6). The subcategory identified in the cellular component category was extracellular region in both 2 vs 0 HAP and 6 vs 0 HAP. For the molecular function category, peroxidase activity, heme binding, antioxidant activity, and oxidoreductase activity, acting on peroxide as acceptor were the most abundant categories in both 2 vs 0 HAP and 6 vs 0 HAP. The most abundant biological processes categories identified in both 2 vs 0 HAP and 6 vs 0 HAP were cell wall organization or biogenesis, phenylpropanoid metabolic process, sulfur compound biosynthetic process, hydrogen peroxide catabolic process, and ROS metabolic process. In addition, no GO terms were significantly enriched in the cor-DEGs-DEPs genes at 6 vs 2 HAP.
KEGG pathway analysis can help us further understand the biochemical metabolic pathways and signal transduction pathways involved in the cor-DEGs-DEPs genes (Fig. 6; Table 1; Table 2; Additional file 7). The results showed that two KEGG pathways were highly enriched at both mRNA and protein levels including phenylpropanoid biosynthesis (ko00940) and stilbenoid, diarylheptanoid, and gingerol biosynthesis (ko00945) in both 2 vs 0 HAP and 6 vs 0 HAP. For the cor-DEGs-DEPs genes at 2 vs 0 HAP, cutin, suberin, and wax biosynthesis (ko00073) and flavonoid biosynthesis (ko00941) were significant pathways in both the proteome and transcriptome. In addition, no KEGG pathways were significantly enriched in the cor-DEGs-DEPs genes at 6 vs 2 HAP.
Parallel reaction monitoring (PRM) analysis
Four differentially expressed proteins (mainly related to flavonoid biosynthesis, peroxidase activity and phenylpropanoid biosynthesis) were chosen for PRM analysis. According to the relative expression quantity of the corresponding peptide fragment of four target proteins in peptide mixtures of the stigmas of 0, 2, and 6 HAP, the relative expression quantity differences of target proteins were obtained (Table 3). Detailed protein quantitative information and significant difference analysis results are shown in Additional file 8. The results from this analysis indicated that expression quantities of the four target proteins in the 2 vs 0 HAP and 6 vs 0 HAP comparisons were markedly up-regulated, whereas the expression quantity of the four target proteins in the 6 vs 2 HAP comparison was not significantly changed. The results of the PRM analysis indicated that the four candidate proteins show similar trends as the label-free results, which supported the credibility of the proteomics data.
Low pollen-pistil compatibility during water lily breeding are associated with the metabolism of ROS
In this study, the combined transcriptome and proteome analysis showed that the expression of genes and proteins related to the metabolism of ROS on the stigma increased significantly in the 2 vs 0 HAP and 6 vs 0 HAP comparisons, suggesting that ROS may be involved in regulating the interaction between the pollen and stigma of water lily after pollination. ROS participate in many pollen-related processes, such as tapetum and pollen development [18,19,20], in vitro pollen germination , growth of the pollen tube apex [22,23,24], the rupture of the pollen tube to release sperm , and self-incompatibility . The role of ROS in pollen tube growth has been well established, but little is known about its involvement in the early stage of pollen germination. The biological function of ROS and hydrogen peroxide on the stigma may be involved in some signal crosslinks in the interaction between the pollen and stigma [27, 28]. Numerous experiments have shown that mature pollen grains produce a large amount of Nitric Oxide (NO), which inhibits ROS production in stigma papilla cells [27, 29, 30]. The mutual exclusion of ROS and NO during pollen tube growth suggests that there may be a coordination mechanism between these signaling molecules during pollen tube growth ; this which indicates that ROS from the stigma and NO from the pollen participate in the pollen-stigma interaction as signaling molecules . In addition, ROS are mainly composed of hydrogen peroxide on the stigma, which is considered the most important redox signaling molecule because of its unique physical and chemical properties as well as its stability in cells. Hydrogen peroxide can oxidize the thiol group of target protein cysteine, thus changing the structure and function of proteins .
However, the regulatory mechanism of ROS in the interaction between pollen and stigma of water lily after pollination is unclear. We infer that low pollen-pistil compatibility after pollination is due to the change of the level of ROS on the stigma. ROS may act as a signaling molecule to oxidize downstream target proteins. The oxidized target proteins cannot function properly, affecting the pollen germination on the stigma surface. Therefore, the metabolic process of ROS on the stigma may be related to the interaction between the pollen and the stigma of water lily.
Effects of stigma flavonoids on low pollen-pistil compatibility during water lily breeding
Flavonoids, which can affect plant physiology, growth, and development, are common secondary metabolites in vascular plants [34, 35]. Flavonoids are mainly involved in the physiological processes of plant symbiosis, defense against disease and insect pests, auxin transport, seed and pollen germination, and root development. Flavonoids are needed for pollen germination and stimulate pollen tube growth in some plant species [36, 37]. Using pollen produced by in vitro culture of immature pollen of tobacco (Nicotiana tabacum L.) and petunia (Petunia hybrida) as materials, adding flavonol (quercetin, kaempferol, myricetin) to germination medium could significantly promote pollen germination frequency and pollen tube growth in vitro . Similarly, the pollen of a flavonoid-deficient mutant of petunia could not germinate on the stigma, but the addition of the exogenous flavonoid kaempferol could induce pollen germination on the stigma. This indicates that kaempferol may play a role in pollen germination . In this study, our transcriptome and proteome analyses showed that there were significant differences in genes and proteins related to flavonoid biosynthesis between stigmas at 0 HAP and stigmas at 2 and 6 HAP. For example, the expression and content of flavonol synthase in stigmas after pollination are lower than those of unpollinated stigmas, which indicates that the flavonoid content on the stigma after pollination is greatly reduced. In addition, flavonoids can affect plant reproductive and developmental processes and participate in the interaction between the pollen and stigma [38, 39]. Other transcriptomic studies on rice stigmas revealed that numerous genes encoding flavonols were expressed, and these genes were expressed on rice stigmas, suggesting that flavonoids play an important role in the interaction between pollen and stigma . Therefore, we infer that the biosynthesis of flavonoids is closely related to the low pollen-pistil compatibility during water lily breeding.
Effects of SAM metabolism on low pollen-pistil compatibility during water lily breeding
The SAM participates in many important physiological processes, such as transamination of propyl, methyl, and sulfur in plants, and is the main hub of methionine metabolism [40, 41]. The results in the present study showed that the SAM synthase gene was up-regulated by pollination. SAM can also be used as a precursor of ethylene and polyamine . Moreover, polyamines with appropriate concentrations are important for pollen germination and pollen tube growth . Similarly, a previous study has shown that ethylene plays a role in pollen tube growth after pollination in tobacco . We infer that up-regulation of genes related to S-adenosylmethionine synthesis and metabolism induced by pollination may facilitate the synthesis of polyamines and ethylene by the stigma papilla cells, and then regulate pollen germination and pollen tube growth. Thus, it is speculated that the metabolic pathway of SAM participates in the interaction between the pollen and stigma of water lily and plays an important role in the regulation mechanism of the pollen-stigma interaction.
Cell wall organization or biogenesis is associated with low pollen-pistil compatibility during water lily breeding
During the interaction between the pollen and stigma, some enzymes in stigma papilla cells are activated and released by certain signals. These enzymes are mainly involved in modifying cell walls, such as enzymes that degrade pectin, cellulose, and hemicellulose [45,46,47]. In our study, the combined transcriptome and proteome analysis showed that many DEGs involved in cell wall tissue metabolism were found on stigma at 2 and 6 h after pollination, suggesting that genes involved in cell wall synthesis were induced by pollination. Transcriptome studies on the stigma of many species have confirmed that there are a large number of cell wall metabolism-related genes in the stigma [48,49,50], and the products of these genes may be secreted by stigma papillae to help pollen tubes penetrate the stigma. In addition, the GO annotation showed that the expression of endoglucanase involved in cellulose hydrolysis was significantly down-regulated, while beta-galactosidase and xyloglucan endotransglucosylase/hydrolase genes involved in cell wall modification were significantly up-regulated. Thus, we infer that beta-galactosidase may be involved in cellulose synthesis and cell wall elongation during cell wall metabolism of Nymphaea stigmas , whereas xyloglucan endotransglucosylase/hydrolase is mainly involved in cell wall reinforcement . The KEGG analysis showed that the gene expression abundances involved in cutin, suberin, and wax biosynthesis were significantly increased, which resulted in cutin, suberin, and wax deposition on the cell walls of stigma cells, thereby increasing cell wall strength . Significant up-regulation of key genes involved in phenylpropanoid biosynthesis further confirmed the accumulation of lignin in stigma the cell wall of water lily after pollination because suberin biosynthesis is closely related to phenylpropanoid biosynthesis . We infer that the accumulation of cutin, suberin, and wax on the cell wall enhances the strength and thickness of the cell wall, thus hindering the pollen germination on the stigma. Therefore, cell wall organization or biogenesis is closely related to low pollen-pistil compatibility of Nymphaea species.
Other enzymes involved in low pollen-pistil compatibility of Nymphaea spp
In this study, the combined analysis of transcriptome and proteome showed that the activity of peroxidase increased significantly after pollination. Peroxidase mainly removes peroxides and participates in stress response, auxin metabolism, and signal transduction . Previous research has shown that the peroxidase activity of a mature stigma increased significantly and reached the highest value when the stigma developed was the most receptive to pollen, which is a common method used to judge stigma maturity in the field . Mc Innis was the first to discover stigma-specific peroxidase, which has cell specificity and specific expression patterns, and is specifically expressed in the cytoplasm and cell surface of stigma epidermal cells . Based on our results and previously published papers, we infer that in the process of interaction between the pollen and stigma of Nymphaea, peroxidase may directly participate in the process of mutual recognition between the stigma and pollen, perhaps guiding the pollen tube.
Glutathione S-transferase (GST) is a ubiquitous enzyme that regulates plant secondary metabolism, detoxification, and defense, and plays an important role in cell response to stress . A previous study on maize showed that the expression of GST was up-regulated at the early stage of silk-pollen interaction and pollen tube germination . GST was also upregulated in this study. It is possible that GST participates in the pollen-stigma interaction of Nymphaea, but its mechanism of involvement needs further study.
Cytochrome P450 has a wide range of catalytic activities. It mainly participates in the synthesis and metabolism of terpenoids, alkaloids, sterols, fatty acids, plant hormones, signal molecules, phenylpropane, flavonoids, and isoflavones . Our results showed that the expression of cytochrome P450 increased significantly. We deduce that cytochrome P450 is involved in the processes of phenylpropanoid and flavonoid biosynthesis, and indirectly participates in the pollen-stigma interaction of water lily.
Mitogen-activated protein kinase (MAPK) is an evolutionarily conserved serine/threonine protein kinase in eukaryotic organisms. It is responsible for regulating signal transduction between cells and amplifying stimulus signals from outside cells to induce appropriate physiological and biochemical reactions in receptor cells . In this study, we found that the expression of MAPK on the stigma of water lily was significantly increased after pollination. We infer that MAPK plays a signaling role in the pollen-stigma interaction and plays an important role in the regulation of the complex network of low pollen-pistil compatibility in water lily.
In this paper, the differences of stigma transcripts and proteomes at 0, 2, and 6 h after pollination were compared, resulting in the identification of some regulatory genes and functional proteins that may cause pre-fertilization barriers in water lily. The functional analysis showed that differential transcripts were mainly involved in cell wall organization or biogenesis, SAM, hydrogen peroxide decomposition and metabolism, ROS metabolism, secondary metabolism, secondary metabolite biosynthesis, and phenylpropanoid biosynthesis. Our transcriptomic and proteomic analysis highlighted specific genes, incuding those in ROS metabolism, biosynthesis of flavonoids, SAM metabolism, cell wall organization or biogenesis and phenylpropanoid biosynthesis that warrant further study in investigations of the pollen-stigma interaction of water lily. This study strengthens our understanding of the mechanism of low pollen-pistil compatibility in Nymphaea at the molecular level, and provides a theoretical basis for overcoming the pre-fertilization barriers in Nymphaea in the future.
An interspecific cross between the female N. ‘Peter Slocum’ and male N. micrantha was performed as described in a previous report . Self-pollinated ‘Peter Slocum’ and N. micrantha were successful. In addition, ‘Peter Slocum’ has a chromosome number of 84 (2n = 6x), whereas N. micrantha has a chromosome number of 56 (2n = 4x) . These plants were grown in ponds in Xingxiang, Zhenjiang, Jiangsu Province, China. These plants were collected from shanghai chenshan plant science research center, Chinese Academy of Sciences (Shanghai, China).
The stigmas of 0, 2, and 6 HAP were collected. Each treatment had three biological repeats. The stigmas of 0 HAP (fresh male parent pollen was placed on stigma) were used as the control, and the stigmas from 2 and 6 HAP were used as the treatment to dynamically study the interaction between the pollen and stigma after pollination. After collection, the three samples, the stigmas of 0 HAP, the stigmas of 2 HAP, and the stigmas of 6 HAP, were immediately frozen in liquid nitrogen and stored at − 80 °C.
High-throughput RNA-seq and data processing
Total RNA was extracted using Trizol reagent according to the manufacturer’s protocol (Takara Bio Inc., Otsu, Japan). The total RNA was checked for quality and quantity using an Agilent 2100 bioanalyzer (Agilent Technologies, CA, USA). The mRNA samples were enriched by oligo (dT) magnetic beads and then cut into fragments with fragmentation buffer at 80 °C. The first strand of cDNA was synthesized with random hexamer as primer. RNase H, dNTPs and DNA polymerase I were used to synthesize the second strand of cDNA using the first strand as template. After purification and end repair, adaptor sequences and double-stranded DNA poly A were attached to the ends of cDNAs. After selecting the fragment size and using Agilent 2100 biological analyzer for quality detection, the cDNA library was constructed by PCR amplification. Finally, the Illumina HiSeq 2500 system was used to sequence the qualified cDNA libraries. Three biological replicates were used in the RNA-seq experiments involving each sample, and sequenced as separate libraries.
Raw data (raw reads) of fasta format were firstly filtered using fastp software. In this step, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated. All the downstream analyses were based on the clean data with high quality. The transcripts were assembled by Trinity (http://trinityrnaseq.github.io). Clean reads with a certain overlap length were initially merged into long fragments without N (called contigs). TGICL software was used to cluster the related contigs to obtain unigenes (without N) which could not be extended at both ends, and nonredundant unigenes were acquired after removing redundancies . The functional annotation of all-unigenes was performed using a BLAST search (http://blast.ncbi.nlm.nih.gov/Blast.cgi) against the GO, Pfam, KEGG, Nr, and Swiss-Prot databases.
The unigenes obtained using Trinity were used as reference sequences, and the clean reads of each sample were compared to the reference sequences by using the default parameters of bowtie 2 in RSEM software. The results of bowtie2 comparison were statistically analyzed by RSEM, and the number of read counts of each sample compared to each gene was further obtained, and the number of read counts was converted to FPKM, then the gene expression level was analyzed. The FPKM (fragments per kilobase per million mapped reads) method eliminated the effect of gene length and sequencing depth on gene expression calculation. The difference of gene expression between samples was directly compared with FPKM values. The “base mean” value for DEGs identification was obtained by using the DESeq package. The absolute value of log2ratio ≥ 1 and FDR ≤0.01 were used as the threshold of significant difference in gene expression between the two samples.
Functional annotation and classification of the DEGs were conducted using the Blast 2 GO program (https://www.blast2go.com/) . Additionally, a KEGG pathway analysis (https://www.genome.jp/kegg/) was performed. The heat map was produced using Cluster 3.0 and treeview.
Label free analysis of the stigma proteome of water lily
Protein extraction and peptide enzymolysis
Protein extraction was performed using the SDT (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0) method. The protein concentrations were quantified with the BCA Protein Assay Kit (Bio-Rad, USA), and the samples were stored at − 80 °C. Protein digestion (200 μg for each sample) was performed using the filter aided proteome preparation procedure described by Wisniewski . The peptides from each sample were desalted on C18 cartridges, concentrated by vacuum centrifugation, and reconstituted in 40 μL of 0.1% (v/v) formic acid.
MS/MS protein identification and quantification
NanoLC-MS/MS analysis was performed on each peptide fraction. We first loaded the peptide mixture onto reversed-phase trap column (Thermo Fisher Scientific Acclaim PepMap100, 100 μm × 2 cm, nanoViper C18) and connected it to C18-reversed phase analytical column (Thermo Fisher Scientific Easy Column, 75 μm inner diameter, 10 cm long, 3 μm, C18-A2) in buffer A (0.1% formic acid), and then separated it with the linear gradient of buffer B (0.1% formic acid, 84% acetonitrile) under the flow rate of 300 nL / min controlled by intelliFlow technology.
The LC-MS/MS analysis was carried out on the Q Exactive mass spectrometer (Thermo Fisher Scientific) coupled with a Easy nLC (Proxeon Biosystems, now Thermo Fisher Scientific) for 240 min. Mass spectrometer was operated under positive ion mode. The MS data was obtained through the data-dependent top10 method and dynamically selected the most abundant precursor ions from the measurement scan (300–1800 m/z) to obtain higher collision energy dissociation (HCD) fragments. The automatic gain control target, the maximum injection time and dynamic exclusion duration was set to 1e6, 50 ms and 60.0 s respectively. Measurement scans were acquired at 70,000 resolution at 200 m/z; HCD spectra resolution was set to 17,500 at 200 m/z, the normalized collision energy was 30 eV, the isolation width was 2 m/z, and the under-fill ratio was 0.1%.
The MS raw files were dealed with using Maxquant22.214.171.124 software for protein identification . In this study, we searched the acquired MS/MS spectra based on the predicted protein databases translated from the above transcriptome databases. The maximum FDR and the minimum peptide length was set to 1% and six amino acids for both peptides and proteins respectively. Other parameters were set as described below: enzyme = trypsin; max missed cleavage = 2; fixed modification: carbamidomethyl (C); peptide mass tolerance = ± 20 ppm; variable modification: oxidation (M), acetyl (protein N-term). Protein quantification was determined by unique peptides and ‘razor’ [62, 63], and label free quantitation calculation was performed . For each fraction, the peptide was matched in different LC-MS/MS operations based on the mass and retention time (setting matching options between runs in MaxQuant) using a 2 min time window.
DEPs were analyzed as significantly up-regulated or down-regulated. For quantitative changes, a critical value of 2.0-fold was set to determine the up-regulated and down-regulated proteins, at least two replicates with a p-value of <0.05.
Each sample was repeated three times in label free analysis. In order to study the effect of differentially expressed proteins on different biological processes, we carried out enrichment analysis of GO and KEGG. Go and KEGG were enriched by Fisher’s exact test, and the whole quantitative protein annotation was used as the background data set. Benjamini-Hochberg was modified by multiple testing, so that the derived p-values could be adjusted . Only pathways and functional classifications with p-values <0.05 were identified as significant.
Correlation analysis of transcriptomics and proteomics
In a comparison group, if a gene and its corresponding protein were expressed, it is considered that the gene and its corresponding protein were correlated. Next, we determined the significant expression of correlated genes and proteins in the comparison group. If a gene and its corresponding protein were significantly expressed in a comparison group, they were named cor-DEGs-DEPs. GO and KEGG enrichment analysis for cor-DEGs-DEPs were then conducted, and the Spearman correlation coefficient was calculated.
Selection of target peptides for PRM analysis
Peptide mixture of the nine samples analyzed by label free analysis was prepared by trypsin. The equivalent peptides in each sample was collected, and 2 μg of the collected sample was introduced into the HPLC system through a capture column (5 μm-C18, 100 μm × 50 mm) and an analysis column (3 μm-C18, 75 μm × 200 mm). Next, the separated peptides were analyzed by Q-Exactive mass spectrometer (Thermo Fisher Scientific). Maxquant 126.96.36.199 software was used to analyze the raw files (missed cleavage = 0, enzyme = trypsin/P). The peptide with a score of more than 40 was considered to be the target peptide.
Quantitative PRM analysis for target proteins
We selected eight target peptides from the four DEPs for PRM quantitative analysis. We added the peptide retention time calibration mixture to the peptide mixture and used the labeled peptide “TASEFDSAIAQDK” (bold “K” for heavy isotope labeling) as the internal standard. We first separated 2 μg of peptide mixture with 20 fmol labeled peptides by HPLC, and then analyzed them by Q-Exactive mass spectrometer. Three times of repeated quantitative analysis were carried out, and raw data was calculated with skyline 3.5.0 software.
Availability of data and materials
The raw data from the three samples have been submitted separately to the National Center for Biotechnology Information (NCBI) under the accession number PRJNA548276.
Differentially expressed genes
Differentially expressed proteins
Fragments Per Kilobase Million
Hours after pollination
Kyoto Encyclopedia of Genes and Genomes
Liquid Chromatography-Parallel Reaction Monitoring/Mass Spectrometry
Mitogen-activated protein kinase
No significant difference in genes
No significant difference in proteins expression
Parallel reaction monitoring
Reactive oxygen species
Huang GZ, Deng HQ, Li Z, Li G. Water lily. Beijing: China Forestry Press; 2009.
Sun CQ, Ma ZH, . Sun GS, Dai ZL, Teng NJ, Pan YP. Cellular mechanisms of reproductive barriers in some crosses of water lily (Nymphaea spp.) cultivars. Hortic Sci. 2015; 50:30–35.
Sun CQ, Ma ZH, Zhang ZC, Sun GS, Dai ZL. Factors influencing cross barriers in interspecific hybridizations of water lily. J Am Soc Hortic Sci. 2018;143(2):1–6.
Pairat S, Vipa H. Intersubgeneric cross in Nymphaea spp. L. To develop a blue hardy waterlily. Sci Hortic. 2010;124:475–81.
Sun CQ, Chen FD, Teng NJ, Liu ZL, Fang WM, Hou XL. Factors affecting seed set in the crosses between Dendranthema grandiflorum (Ramat) Kitamura and its wild species. Euphytica. 2010;171:181–92.
Teng NJ, Wang YL, Sun CQ, Fang WM, Chen FD. Factors influencing fecundity in experimental crosses of water lotus (Nelumbo nucifera Gaertn.) cultivars. BMC Plant Biol. 2012;12:82.
De Graaf BHJ, Derksen JWM, Mariani C. Pollen and pistil in the progamic phase. Sex Plant Reprod. 2001;14:41–55.
Ram SG, Ramakrishnan SH, Thiruvengadam V, Bapu JRK. Prefertilization barriers to interspecific hybridization involving Gossypium hirsutum and four diploid wild species. Plant Breed. 2008;127:295–300.
Ram SG, Sundaravelpandian K, Kumar M, Vinod KK, Bapu JRK, Raveendran TS. Pollen-pistil interaction in the inter-specific crosses of Sesamum sp. Euphytica. 2006;152:379–85.
Mazzucato A, Olimpieri I, Ciampolini F, Cresti M, Soressi GP. A defective pollen-pistil interaction contributes to hamper seed set in the parthenocarpic fruit tomato mutant. Sex Plant Reprod. 2003;16:157–64.
Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57–63.
Sun QQ, Zhang N, Wang JF, Cao YY, Li XS, Zhang HJ, Zhang L, Tan DX, Guo YD. A label-free differential proteomics analysis reveals the effect of melatonin on promoting fruit ripening and anthocyanin accumulation upon postharvest in tomato. J Pineal Res. 2016;61(2):138–53.
Ma DY, Huang X, Hou JF, Ma Y, Han QX, Hou GG, Wang CY, Guo TC. Quantitative analysis of the grain amyloplast proteome reveals differences in metabolism between two wheat cultivars at two stages of grain development. BMC Genomics. 2018;19:768.
Raherison E, Rigault P, Caron S, Poulin PL, Boyle B, Verta JP, Giguère I, Bomal C, Bohlmann J, MacKay J. Transcriptome profiling in conifers and the PiceaGenExpress database show patterns of diversification within gene families and interspecific conservation in vascular gene expression. BMC Genomics. 2012;13:434.
Qiu ZB, Wan LC, Chen T, Wan YL, He XQ, Lu SF, Wang YW, Liu JX. The regulation of cambial activity in Chinese fir (Cunninghamia lanceolata) nvolves extensive transcriptome remodeling. New Phytol. 2013;199:708–19.
Sheoran IS, Pedersen EJ, Ross AR, Sawhney VK. Dynamics of protein expression during pollen germination in canola (Brassica napus). Planta. 2009;230(4):779–93.
Samuel MA, Tang W, Jamshed M, Northey J, Patel D, Siu KW, Muench DG, Wang ZY, Goring DR. Proteomic analysis of Brassica stigmatic proteins following the self-incompatibility reaction reveals a role for microtubule dynamics during pollen responses. Mol Cell Proteomics. 2011;10(12):111–3.
Hu L, Liang W, Yin C, Cui X, Zong J, Wang X, Hu J, Zhang D. Rice MADS3 regulates ROS homeostasis during late anther development. Plant Cell. 2011;23:515–33.
Luo D, Xu H, Liu Z, Guo J, Li H, Chen L, Fang C, Zhang Q, Bai M, Yao N. A detrimental mitochondrial-nuclear interaction causes cytoplasmic male sterility in rice. Nat Genet. 2013;45:573–7.
Xie HT, Wan ZY, Li S, Zhang Y. Spatiotemporal production of reactive oxygen species by NADPH oxidase is critical for tapetal programmed cell death and pollen development in Arabidopsis. Plant Cell. 2014;26:2007–23.
Speranza A, Crinelli R, Scoccianti V, Geitmann A. Reactive oxygen species are involved in pollen tube initiation in kiwifruit. Plant Biol. 2012;14:64–76.
Potocky M, Jones MA, Bezvoda R, Smirnoff N, Zarsky V. Reactive oxygen species produced by NADPH oxidase are involved in pollen tube growth. New Phytol. 2007;174:742–51.
Kaya H, Nakajima R, Iwano M, Kanaoka MM, Kimura S, Takeda S, Kawarazaki T, Senzaki E, Hamamura Y, Higashiyama T. Ca2+−activated reactive oxygen species production by Arabidopsis RbohH and RbohJ is essential for proper pollen tube tip growth. Plant Cell. 2014;26:1069–80.
Lassig R, Gutermuth T, Bey TD, Konrad KR, Romeis T. Pollen tube NAD (P) H oxidases act as a speed control to dampen growth rate oscillations during polarized cell growth. Plant J. 2014;78:94–106.
Duan Q, Kita D, Johnson EA, Aggarwal M, Gates L, Wu HM, Cheung AY. Reactive oxygen species mediate pollen tube rupture to release sperm for fertilization in Arabidopsis. Nat Commun. 2014;5:3129.
Wilkins KA, Bancroft J, Bosch M, Ings J, Smirnoff N, Franklin-Tong VE. Reactive oxygen species and nitric oxide mediate actin reorganization and programmed cell death in the self-incompatibility response of papaver. Plant Physiol. 2011;156:404–16.
Mc Innis SM, Emery DC, Porter R, Desikan R, Hancock JT, Hiscock SJ. The role of stigma peroxidases in flowering plants: insights from further characterization of a stigma-specific peroxidase (SSP) from Senecio squalidus (Asteraceae). J Exp Bot. 2006;57(8):1835–46.
Hiscock SJ, Bright J, Mcinnis SM, Desikan R, Hancock JT. Signaling on the stigma: potential new roles for ROS and NO in plant cell signaling. Plant Signal Behav. 2007;2(1):23–4.
Bright J, Hiscock SJ, James PE, Hancock JT. Pollen generates nitric oxide and nitrite: a possible link to pollen-induced allergic responses. Plant Physiol Biochem. 2009;47(1):49–55.
Zafra A, Rodríguezgarcía MI, Alché JD. Cellular localization of ROS and NO in olive reproductive tissues during flower development. BMC Plant Biol. 2010;10(1):36.
Popescu SC, Popescu GV, Bachan S, Zhang Z, Seay M, Gerstein M, Snyder M, Dinesh-Kumar S. Differential binding of calmodulin-related proteins to their targets revealed through high-density Arabidopsis protein microarrays. Proc Natl Acad Sci. 2007;104:4730–5.
Sharma B, Bhatla SC. Accumulation and scavenging of reactive oxygen species and nitric oxide correlate with stigma maturation and pollen-stigma interaction in sunflower. Acta Physiol Plant. 2013;35(9):2777–87.
Tian Y, Fan M, Qin Z, Lv H, Wang M, Zhang Z, Zhou W, Zhao N, Li X, Han C, Ding Z, Wang W, Wang ZY. Hydrogen peroxide positively regulates brassinosteroid signaling through oxidation of the BRASSINAZOLE-RESISTANT1 transcription factor. Nat Commun. 2018;9:1063.
Winkel-shirley B. Flavonoid biosynthesis. A colorful model for genetics, biochemistry, cell biology, and biotechnology. Plant Physiol. 2001;126(2):485–93.
Tsai CJ, Harding SA, Tschaplinski TJ, Lindroth RL, Yuan Y. Genome-wide analysis of the structural genes regulating defense phenylpropanoid metabolism in Populus. New Phytol. 2006;172(1):47–62.
Ylstra B, Touraev A, Moreno RM, Stöger E, van Tunen AJ, Vicente O, Mol JN, Heberle-Bors E. Flavonols stimulate development, germination, and tube growth of tobacco pollen. Plant Physiol. 1992;100:902–7.
Mo Y, Nagel C, Taylor LP. Biochemical complementation of chalcone synthase mutants defines a role for flavonols in functional pollen. Proc Natl Acad Sci. 1992;89:7213–7.
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. Plant Physiol. 2007;144:1797–812.
Cheung AY, Wu HM, Di Stilio V, Glaven R, Chen C, Wong E, Ogdahl J, Estavillo A. Pollen-pistil interantiongs in Nicotiana tabacum. Ann Bot. 2000;85(S1):29–37.
Shen B, Li C, Tarczynski MC. High free-methionine and decreased lignin content result from a mutation in the Arabidopsis S-adenosyl-L-methionnine synthetase 3 gene. Plant J. 2002;29(3):371–80.
Martinez-López N, Varela-Rey M, Ariz U, Embade N, Vazquez-Chantada M, Fernandez-Ramos D, Gomez-Santos L, Lu SC, Mato JM, Martinez-Chantar ML. S-adenosylmethionine and proliferation: new pathways, new targets. Biochem Soc Trans. 2008;36(5):848–52.
Roje S. S-Adenosyl-L-methionine: beyond the universal methyl group donor. Phytochemistry. 2006;67(15):1686–98.
Wolukau JN, Zhang SL, Xu GH, Chen D. The effect of temperature, polyamines and polyamine synthesis inhibitor on invitro pollen germination and pollen tube growth of Prunus mume. Sci Hortic. 2004;99(3–4):289–99.
Martinis DD, Cotti C, Hekker SL. Ethylene response to pollen tube growth in Nicotiana tabacum flowers. Planta. 2002;214(5):806–12.
Tian GW, Chen MH, Zaltsman A, Citovsky V. Pollen-specific pectin methylesterase involved in pollen tube growth. Dev Biol. 2006;294:83–91.
Kim J, Shiu SH, Thoma S, Li WH, Patterson SE. Patterns of expansion and expression divergence in the plant polygalacturonase gene family. Genome Biol. 2006;9:87.
Suen DF, Huang AH. Maize pollen coat xylanase facilitates pollen tube penetration into silk during sexual reproduction. J Biol Chem. 2007;282:625–36.
Swanson R, Clark T, Preuss D. Expression profiling of Arabidopsis stigma tissue identifies stigma-specific genes. Sex Plant Reprod. 2005;18:163–71.
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.
Quiapim AC, Brito MS, Bernardes LA, Dasilva I, Malavazi I, De Paoli HC, Molfetta-Machado JB, Giuliatti S, Goldman GH, Goldman MH. Analysis of the Nicotiana tabacum stigma/style transcriptome reveals gene expression differences between wet and dry stigma species. Plant Physiol. 2009;149:1211–30.
Roach MJ, Mokshina NY, Badhan A, Snegireva AV, Hobson N, Deyholos MK, Gorshkova TA. Development of cellulosic secondary walls in flax fibers requires beta-galactosidase. Plant Physiol. 2011;156:1351–63.
Campbell P, Braam J. Xyloglucan endotransglycosylases: diversity of genes, enzymes and potential wall-modifying functions. Trends Plant Sci. 1999;4(9):361–6.
Huckelhoven R. Cell wall-associated mechanisms of disease resistance and susceptibility. Annu Rev Phytopathol. 2007;45:101–27.
Bernards MA. Demystifying suberin. Can J Bot. 2002;80:227–40.
Dixon DP, Skipsey M, Edwards R. Roles for glutathione transferases in plant secondary metabolism. Phytochemistry. 2010;71:338–50.
Liu HH, Wang LW, Liu X, Ma X, Ning LH, Zhang H, Cui DZ, Jiang C, Chen HB. Proteomic analyses of the early pollen-silk ineraction in maize. Sci Agric Sin. 2010;43(24):5000–8.
Schuler MA, Werck-Reichhart D. Functional genomics of P450s. Annu Rev Plant Biol. 2003;54:629–67.
Group M. Mitogen-activated protein kinase cascades in plants: a new nomenclature. Trends Plant Sci. 2002;7(7):301–8.
Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, Tsai J, Quackenbush J. TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics. 2003;19:651–2.
Conesa A, GÖtz S, García-Gómza JM, Terol J, Talón M, Robles M. Blast2go: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.
Wisniewski JR, Zougman A, Nagraj N, Mann M. Universal sample preparation method for proteome analysis. Nat Methods. 2009;6:359–62.
Cox J, Mann M. MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification. Nat Biotechnol. 2008;26(12):1367–72.
Hubner NC, Mann M. Extracting gene function from protein-protein nteractions using quantitative BAC InteracCtomics (QUBIC). Methods. 2011;53:453–9.
Cox J, Hein MY, Luber CA, Paron I, Nagaraj N, Mann M. Accurate proteome-wide label-free quantification by delayed normalization and maximal peptide ratio extraction, termed MaxLFQ. Mol Cell Proteomics. 2014;13(9):2513–26.
Futcher B, Latter GI, Monardo P, Mclaughlin CS, Garrels JI. A sampling of the yeast proteome. Mol Cell Biol. 1999;19:7357–68.
This study was supported by the National Natural Science Foundation of China (31701948). The funders had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Differential expression genes in 0, 2 and 6 HAP. (XLS 3738 kb)
Identifiers of proteins contained in the protein group. They are sortted by number of identified peptides in descending order.
Quantitative and differential analysis of protein.
The correlation between DEPs and DEGs.
The cor-DEGs-DEPs genes with the same or opposite trend.
GO enrichment analysis of the cor-DEGs-DEPs genes.
KEGG enrichment analysis of the cor-DEGs-DEPs genes.
Detailed protein quantitative information and significant difference analysis results.
About this article
Cite this article
Sun, CQ., Chen, FD., Teng, NJ. et al. Transcriptomic and proteomic analysis reveals mechanisms of low pollen-pistil compatibility during water lily cross breeding. BMC Plant Biol 19, 542 (2019). https://doi.org/10.1186/s12870-019-2166-3