Comparative proteomic analysis of eggplant (Solanum melongena L.) heterostylous pistil development

Heterostyly is a common floral polymorphism, but the proteomic basis of this trait is still largely unexplored. In this study, self- and cross-pollination of L-morph and S-morph flowers and comparison of embryo sac development in eggplant (Solanum melongena L.) suggested that lower fruit set from S-morph flowers results from stigma-pollen incompatibility. To explore the molecular mechanism underlying heterostyly development, we conducted isobaric tags for relative and absolute quantification (iTRAQ) proteomic analysis of eggplant pistils for L- and S-morph flowers. A total of 5,259 distinct proteins were identified during heterostyly development. Compared S-morph flowers with L-morph, we discovered 57 and 184 differentially expressed proteins (DEPs) during flower development and maturity, respectively. Quantitative real time polymerase chain reactions were used for nine genes to verify DEPs from the iTRAQ approach. During flower development, DEPs were mainly involved in morphogenesis, biosynthetic processes, and metabolic pathways. At flower maturity, DEPs primarily participated in biosynthetic processes, metabolic pathways, and the formation of ribosomes and proteasomes. Additionally, some proteins associated with senescence and programmed cell death were found to be upregulated in S-morph pistils, which may lead to the lower fruit set in S-morph flowers. Although the exact roles of these related proteins are not yet known, this was the first attempt to use an iTRAQ approach to analyze proteomes of heterostylous eggplant flowers, and these results will provide insights into biochemical events taking place during the development of heterostyly.


Introduction
In flowering plants, different strategies have evolved to avoid selfing and promote outcrossing, of which heterostyly is one of the most effective mechanisms. Heterostyly, a complex floral polymorphism, can aid in environmental adaptations of plants and accelerate species diversification [1,2]. Heterostyly has arisen independently in at least 20 lineages and is present in 199 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 genera, distributed among 28 families in 15 orders [1,3]. Heterostylous plants usually include two (distyly) or three (tristyly) genetic morphs with reciprocal displacement of sexual organs (stigmas and anthers) within an individual [4]. For example, in eggplant (Solanum melongena L.), plants produce two types of flowers (distyly): either long-styled flowers with anthers attached midway along the floral tube (L-morph or pin), or short-styled flowers with anthers attached at the top of the floral tube (S-morph or thrum). This character promotes outcrossing between morphs via delivery and uptake of pollen by pollinators [5].
Although many angiosperms are heterostylous, only a few differentially expressed genes (DEGs) have been detected for the condition, and the regulatory molecular mechanisms are not well understood. Ushijima et al. [6] elaborated molecular differences by comparing transcripts and proteins in the thrum and pin flowers of Linum grandiflorum Desf. These floral phenotypes were known to be regulated by the S locus and differed in style length, pollen size, and anther length [7]. Four genes, TSS1, AP1, MYB21, and SKS1, were predicted to be related to heterostyly development [6]. However, there was no difference in messenger ribonucleic acid accumulations of these four genes, indicating that they were controlled by floral morphspecific post-transcriptional regulation [6]. Transcriptome analysis for both Primula veris L. and the closely related species P. vulgaris Huds. demonstrated that 113 candidate heterostyly genes showed significant floral morph-specific differential expression [8].
Since the development of RNA sequencing (RNA-seq) technology, it has been easy to document changes in gene expression at the transcriptional level. However, the transcriptional changes are not always directly related to expression of the corresponding proteins because of posttranslational regulatory mechanisms [9,10], alternative splicing, and protein degradation [11]. Proteomic approaches provide valuable tools for monitoring developmental profiles directly at the protein level and therefore have been widely used [12][13][14]; approaches include two dimensional gel electrophoresis (2-DE), differential in gel electrophoresis (DIGE), tagbased labeling of proteins (isotope-coded affinity tags (ICAT)), stable isotope labeling with amino acids in cell culture (SILAC), isobaric tags for relative and absolute quantification (iTRAQ), protein-protein interaction, and protein modifications [15]. iTRAQ is a powerful technology that allows identification of numerous proteins between different samples [13]. Furthermore, if iTRAQ returns a sufficient number of differentially expressed proteins (DEPs), pathway and protein-protein interaction analyses can be conducted [13,14,16].
iTRAQ provides advantages in labeling of complex samples, with comparatively high throughput and identification of low-abundance proteins in complex samples [17,18]. Since Ross et al. [19] first published an approach using iTRAQ to examine the global protein expression of a wild-type yeast strain, this technique has been widely used to document quantum changes in DEPs in plants and animals [12][13][14]. Meng et al. [20] then applied iTRAQ in a proteomic study of blood cells infected with Spiroplasma eriocheiris. Yang et al. [14] discovered genes related to grain development in an analysis of wheat grain protein expression at different stages. The study of regulated protein expression levels in different eggplant flower morphs should provide insights into heterostyly developmental mechanisms.
Eggplant is a heterostylous plant that is widely cultivated [21]. The S-morph flowers generally possess a small and highly reduced gynoecium, and are often functionally staminate, limiting production. Therefore, understanding the molecular genetics that regulate heterostyly could lead to improved selection for better production. In this study, we performed an iTRAQ-based quantitative proteome analysis of pistils of two flower morphs from the budding to the blooming stage. The general workflow is shown in S1 Fig. Our results provide information about differences in proteins during heterostylous development and highlight the value of proteomics in characterizing complex biochemical processes.

Floral measurements
Eggplant B3-3 lines were grown at the Guangxi Academy of Agricultural Sciences and cultivated in a field using conventional methods. Plants produced two morphologically different types of flowers: long morph (L-morph) and short morph (S-morph). We randomly selected 30 S-morph and 40 L-morph flowers from budding to blooming to measure the pistil length and bud length with a Vernier caliper and conducted a correlation analysis using Microsoft Excel 2016.

Stigma-pollen interactions
The flowers were emasculated in bud and bagged to exclude pollinators before and after hand pollination. The stigmas were collected at 1, 2, 4, 6, 8, 12, 24, or 48 h after pollination and fixed for more than 24 h in a mixture of Formalin-acetic acid-ethanol. The samples were washed with distilled water and immersed in 2 molÁL −1 NaOH in a 60˚C water bath for 12 h. Specimens were washed again with distilled water and were stained with 0.1% aniline blue for 4 h, and then mounted on glass slides. The samples were covered with 80% glycerin and observed under an inverted fluorescence microscope (Olympus BX51, Japan).

Female gametophyte development
To investigate whether differences in embryo sac development reduced successful fertilization, we sectioned flowers embedded in paraffin. From budding to blooming, we collected S-morph and L-morph pistils every 2 days and fixed them in a Formalin-acetic acid-ethanol mixture. The samples were dehydrated in an ethanol series (10 min each in 35%, 55%, 75%, 85%, 95%, and 100% [v/v]), cleared in a xylene series (10 min each in 35%, 55%, 75%, 85%, 95%, and 100% [v/v]), and embedded in paraffin (melting point: 54-56˚C) for 48 h at 56˚C. Embedded specimens were serially sectioned at a thickness of 6 μm and mounted on glass slides. They were prepared with a 2% ferrovanadium mordant for 30 min, stained with 5% hematoxylin for 1.5 h, and then destained with saturated picric acid for 1.5 h. Finally, coverslips were mounted with neutral balsam and the slides were observed under an OLYMPUS BX-51 light microscope (Olympus Co. Ltd., Japan).

Protein extraction
We divided flower development into five stages from budding to blooming: 0, 3, 6, 10, or 13 days after budding (DAB). The flowers at 0, 3, and 6 DAB were considered development samples. We collected mixed samples of pistils at 0, 3, and 6 DAB from S-morph and L-morph flowers. Pistils collected at 13 DAB were kept separately as mature flower samples. All tissue samples were stored at −80˚C liquid nitrogen until protein extraction.
The pistil samples were ground into powder under liquid nitrogen and extracted with lysis buffer A (7 M urea, 2 M thiourea, 4% 3-[(3-cholamidopropyl)dimethylammonio]-1-propanesulfonate, 40 mM Tris-HCl, pH 8.5) containing 1 mM phenylmethylsulfonyl fluoride and 2 mM ethylenediaminetetraacetic acid. After 5 min, 10 mM dithiothreitol was added. After sonication and centrifugation, the suspension was mixed well with a 5-fold volume of chilled acetone containing 10% trichloroacetic acid and incubated overnight at −20˚C. After centrifugation (4˚C, 30,000 ×g), the precipitate was washed three times with chilled acetone. The pellet was air-dried and dissolved in lysis buffer B (7 M urea, 2 M thiourea, 4% nonylphenol ethoxylate (NP-40), and 20 mM Tris-HCl, pH 8.5). The suspension was sonicated for 15 min and centrifuged at 4˚C and 30,000 ×g for 15 min. Subsequently, 10 mM dithiothreitol was added to reduce disulfide bonds in proteins in the supernatant, and the solution was incubated at 56˚C for 1 h. Next, 55 mM iodoacetamide was added to bind to cysteines and the solution was incubated for 1 h in the dark. The supernatant was mixed well with a 5-fold volume of chilled acetone for 2 h at −20˚C. After centrifugation (30,000 ×g for 20 min), the pellet was air-dried for 5 min, then dissolved in 500 μL of 0.5 M triethylammonium bicarbonate and sonicated for 15 min. Finally, after centrifugation at 4˚C and 30,000 ×g for 15 min, the supernatant was transferred into a new tube and quantified by Bradford's method [22]. The proteins in the supernatant were stored at −80˚C for further analysis.
iTRAQ labeling and chromatography fractionation Total protein (100 μg) taken from each sample solution was digested with Trypsin Gold (Promega, USA) with a 20:1 ratio of protein:trypsin at 37˚C for 4 h. After trypsinization and drying by vacuum centrifugation, peptides were redissolved using 0.5 M triethylammonium bicarbonate and iTRAQ reagent (Applied Biosystems, USA) according to the manufacturer's instructions. Each group of peptides was marked by different iTRAQ tags and incubated at room temperature for 2 h. We mixed all groups of tagged peptides, purified them using a strong cation exchange chromatography column (Phenomenex, USA), and separated them by liquid chromatography (LC) using a LC-20AB high pressure LC pump system (Shimadzu, Japan). Then we redissolved tagged mixed peptides with 4 mL of buffer A (25 mM NaH 2 PO 4 in 25% acetonitrile (ACN), pH 2.7) and loaded them onto a 4.6 × 250 mm Ultremex strong cation exchange column containing 5 mm particles (Phenomenex). Gradient elution was applied to peptides at a flow rate of 1 mL min −1 , in which we initially used buffer A for 10 min elution and then progressively interfused 5-35% buffer B (25 mM NaH 2 PO 4 , 1 M KCl in 25% ACN, pH 2.7) for 11 min elution. Finally, we conducted 1 min elution with 35-80% buffer B. The entire elution process was monitored by measuring the absorbance at 214 nm, and each component was desalted with a Strata X C18 column and vacuum dried.

LC-electrospray ionization-tandem mass spectrometry analysis
A nanoACQuity (Waters, USA) rapid separation LC system connected with the mass spectrometer, including a Symmetry C18 column (5 μm, 180 um × 20 mm), was used for peptide absorption and desalting, and a BEH130 C18 column (1.7 μm, 100 um × 100 mm) was used for separation. Both mobile phase buffer A (98:2:0.1 H 2 O:ACN:HCOOH) and buffer B (2:98: 0.1 H 2 O:ACN:HCOOH) were added with a certain ratio of correction fluid (Thermo Fisher Scientific, USA). A 2.25 μg (9 μl) amount was loaded each time. Peptide absorption and desalting were carried out with buffer A at a flow rate of 2 μL min −1 for 15 min elution. The samples were loaded with 5% buffer B at 300 nL min −1 for 1 min, and then a 40 min gradient was run starting with 5-35% buffer B, followed by 5 min of linear gradient to 80%, followed by 5 min of maintenance at 80%, and a final 2 min at 5%.
Data were acquired with a TripleTOF 5600 System (AB SCIEX, Concord, ON) fitted with a Nanospray III source (AB SCIEX), with a pulled quartz tip as the emitter (New Objectives, Woburn, MA), controlled by the software program Analyst 1.6 (AB SCIEX). The following mass spectrometry conditions were used: 2.5 kV ion spray voltage, 30 psi curtain gas, 15 psi nebulizer gas, and 150˚C interface heater temperature. The resolution was approximately 30,000. For independent data acquisition, survey scans were acquired in 250 ms and as many as 30 product ion scans were collected if they exceeded a threshold of 120 counts s −1 and had a 2+ to 5+ charge state. The total cycle time was fixed at 3.3 s. The second quadrupole transmission window was 100 Da for 100%. Four time bins were summed for each scan at a pulsed frequency value of 11 kHz through monitoring of a 40 GHz multichannel time-to-digital detector with four anode channels. An adjusted iTRAQ rolling collision energy was applied to all precursor ions for collision-induced dissociation. Dynamic exclusion was set for 1/2 of peak width (15 s), and then the precursor was refreshed off the exclusion list.

Protein identification and bioinformatics analysis
Raw mass spectrum files were converted into Mascot generic files and protein identification was conducted using Mascot (version 2.3.02) to search for predicted proteins in the Eggplant Genome DataBase (http://eggplant.kazusa.or.jp/). For further functional analysis, differential expression of proteins was analyzed for significant downregulation or upregulation. A change in expression was determined by comparing the S-morph and L-morph flower pistils during two developmental stages, and t-tests were used to identify significant (p < 0.05) differences. The proteins with an average fold-change ! 1.5 or 0.667, and unique proteins with at least two peptide matches, were confidently defined as differentially expressed proteins (DEPs).
Proteins were classified by Gene Ontology (GO) analysis with AmiGO 2 (http://amigo. geneontology.org/amigo) based on three categories: biological process, cellular component, and molecular function. The Kyoto Encyclopedia of Genes and Genomes (KEGG) was used to annotate pathways in the KEGG pathway database (http://www.genome.jp/kegg/pathway. html) [23]. The statistical significance of GO terms or KEGG pathway enrichment was determined by a hypergeometric test following Yu et al. [24]. Additionally, protein-protein interaction networks for DEPs were explored using the publicly available Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database [25].

Quantitative real-time PCR analysis
We selected genes from the DEPs of developing and mature flowers to validate the high throughput data at a transcriptional level by qRT-PCR. GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was used as a reference gene and the primers used in this study are listed in S1 Table. We added 10 μL 2× Master Mix, 0.6 μL 10 μM forward primer, 0.6 μL 10 μM reverse primer, 0.8 μL complementary DNA (cDNA) template reverse transcripted to a PCR reaction system (Fermentas), and enough diethyl pyrocarbonate water for a total volume of 20 μL. A Stepone Plus thermocycler was used (Applied Biosystems, Foster City, CA). The amplification protocol consisted of 95˚C denaturation for 10 min; 40 cycles of 95˚C for 15 s, 60˚C for 30 s, and 72˚C for 15 s, and 95˚C for 15 s, 60˚C 1 min, and 95˚C 15 s for a melting curve and PCR specificity test. The analysis was performed using two independent cDNA preparations in triplicate PCR reactions. The relative expression ratio was calculated using the 2 -ΔΔCt method [26] with eggplant GAPDH as the internal reference gene.

Phenotypic differences between L-morph and S-morph flowers in eggplant
Heterostyly is a type of flower polymorphism that leads to separation of the stigma and anthers (hercogamy), preventing self-pollination, and is widely distributed in angiosperms [1,27,28]. In this study, two morphologically different types of flowers, L-morph and S-morph, were observed at different developmental stages in eggplant. As shown in Fig 1A, the flowers are clearly different, especially in their pistil length. During flower development, the pistil length increased linearly with flower bud length in L-morph flowers (R 2 = 0.974) (Fig 1B). Although the S-morph flowers experienced similar growth when the pistil length was less than 10 mm (R 2 = 0.968), the pistil generally did not elongate when the buds were more than 10 mm long  (Fig 2B). With time, pollen tube elongation and germination increased ( Fig  2C). However, after self-pollination of S-morph flowers, the pollen failed to germinate, and not as much pollen adhered to the stigma (Fig 2D, 2E and 2F), indicating incompatibility in self-pollination. Hence, compared with the L-morph flowers, pollen does not adhere as well to the stigma of S-morph flowers.
To determine whether the germination failure of pollen in S-morph flowers originated with the pollen or pistil, we reciprocally cross-pollinated L-morph and S-morph flowers. Germination occurred 24 h after cross-pollination of S-morph pollen to L-morph stigmas ( Fig 2G). Smorph pollen germinated and pollen tubes grew normally on the stigmas of L-morph flowers, which suggested that of pollen from S-morph flowers can function normally. However, no obvious germination occurred 24 h after cross-pollination of L-morph pollen to S-morph stigmas ( Fig 2H). No penetration of pollen tubes into the short style was detected, nor did any pollen grains germinate on the stigma. Given that the pollen of L-morph flowers can effectively germinate on L-morph flowers, this may indicate that the failure of pollen germination stems from the S-morph stigmas. Therefore, we inferred that the structure of the stigmatic surface in S-morph flowers may inhibit pollen germination, and the lower fruit set of S-morph flowers may reflect a shift toward functionally staminate flowers.

Structural observations of pistils in L-morph and S-morph flowers
We examined eggplant embryo sacs of both S-morph and L-morph flowers from budding to blooming to better understand early embryonic development and compare the structural integrity of egg cells, central cells, antipodal cells, and synergids (Fig 3). The egg cell and two synergids were located at the micropylar end of the nucellus. The central cell occupied the majority of the embryo sac. Three smaller antipodal cells were located at the chalazal end. The structure of embryo sacs from S-morph and L-morph flowers did not differ; we therefore inferred that low fertilization rates of S-morph flowers probably resulted from stigmatic incompatibility, rather than defects in the embryo sac.

Protein identification overview
We conducted a comparative proteome survey using the iTRAQ technique to detect molecular expression differences between S-morph and L-morph flowers during heterostyly development. We identified 16,460 high-quality unique peptides from 346,625 secondary spectra. We searched the Eggplant Genome DataBase (http://eggplant.kazusa.or.jp/) [29] and compared the peptides with predicted proteins using Mascot 2.3.02, identifying 4,728 proteins (S2A Fig). in L-morph and S-morph flowers. The green dots indicate the relationship between bud length and pistil length in Lmorph flowers. In S-morph flowers, the relationship is indicated in red when bud length < 10 mm and blue when bud length > 10 mm. https://doi.org/10.1371/journal.pone.0179018.g001 We classified 16,460 peptides in 42 categories based on peptide length; the most frequent category (~10%) was 9-11 amino acids (S2B Fig). The majority of the 4,728 identified proteins included fewer than 10 peptides, and as the peptide number increased, the number of corresponding proteins decreased (S2C Fig). All the proteins with a false discovery rate (FDR) less than 1% were included in downstream analyses including (GO, Cluster of Orthologous Groups of proteins (COG) and KEGG Pathway).
COG is a database that is generated by comparing predicted and known proteins in all completely sequenced genomes to infer sets of orthologs. We compared proteins in our study with the COG database to predict possible functions and conduct a functional classification analysis (S3A Fig). The 4,728 identified proteins were categorized as involved in post-translational modification, protein turnover, chaperone functions, translation, ribosomal structure and biogenesis, carbohydrate metabolism and transport, energy production and conversion, and amino acid metabolism and transport. Additionally, many proteins (! 100) were involved in lipid transport and metabolism, transcription, replication, recombination and repair, signal transduction mechanisms, cell wall/membrane/envelope biogenesis, secondary metabolite biosynthesis, transport, and catabolism.
GO is a community-based bioinformatics resource that standardizes descriptions of functions and classifies gene product functions through the use of structured, controlled vocabularies for associated biological processes, cellular components, and molecular functions. As shown in S3B Fig, the highest percentages of GO terms related to biological processes were cellular progress, metabolism progress, and single-organism process, while those related to cellular components were cell, cell part, and organelle. The highest percentages of GO terms related to molecular functions were binding and catalytic activity.
The KEGG pathway database is a collection of pathway maps representing our knowledge of molecular interaction and reaction networks [23]. A total of 135 pathways were annotated for 3,380 of the 4,728 proteins. More than 5% of the identified proteins belonged to pathways of metabolism, biosynthesis of secondary metabolites, carbon metabolism, and biosynthesis of amino acids (S3C Fig).

Analysis of protein expression patterns in L-morph and S-morph flowers during development and maturity
Here, we defined proteins with expression level fold changes > 1.5 and p-values < 0.05 as DEPs. We analyzed protein expression levels in the pistils of both L-morph and S-morph flowers during development (0-6 DAB) and discovered 24 downregulated and 33 upregulated proteins in S-morph flowers compared with L-morph flowers (Tables 1 and 2). Alanine-glyoxylate aminotransferase 2 (AGT2) homolog 3 (mitochondrial-like) had the greatest downregulation. AGT2 is a peroxisomal photorespiratory enzyme that catalyzes transamination reactions with multiple substrates [30]. Agamous-like MADS-box protein AGL61-like had the greatest upregulation. During maturity (13 DAB), 83 proteins were upregulated and 101 were downregulated in S-morph flowers compared with L-morph flowers (Tables 3 and 4). Polyubiquitin-like had the greatest downregulation. Interestingly, pistil extension-like protein (Sme2.5_02187.1_g00002.1) was downregulated in S-morph flowers during maturity, which may indicate that the short pistil of S-morph flowers is related to its lower expression level. Pectinesterase 2-like was the top differentially upregulated protein. Additionally, during flower maturity, the expression levels of cysteine proteinase and peroxidase, which are involved in senescence and programmed cell death [31], were upregulated in S-morph flowers, which might contribute to the failure of S-morph flowers to set fruit. We found that methionine sulfoxide reductase A3 in S-morph flowers, which is involved in the response to oxidative stress [32], was downregulated during both flower development and maturity compared with Lmorph flowers. Moreover, the cysteine proteinase expression level increased in S-morph pistils and was probably related to cell apoptosis.
We compared expression differences of L-morph and S-morph pistils between development and maturity. For L-morph flowers, there were 322 upregulated and 322 downregulated Heterostyly in eggplant genes during maturity compared with development (S2 and S3 Tables). Among upregulated proteins, olee1-like protein-like, fasciclin-like arabinogalactan protein 14-like, profilin-1-like, profilin-1, and cysteine proteinase 3-like, had the highest expression differences (S2 Table). Pectinesterase 2-like was the top differentially upregulated protein.
Other proteins, such as nonspecific lipid-transfer protein-like protein At2g13820-like were differentially downregulated (S3 Table). In S-morph flowers, there were 446 upregulated and 550 downregulated genes during maturity compared with development (S4 and S5 Tables). Olee1-like proteinlike, nucleoside diphosphate kinase IV, pectinesterase 2-like, profilin-1, and other proteins were the top differentially upregulated proteins (S4 Table). Many downregulated proteins, such as histone H1, ubiquitin extension protein, and cyclic nucleotide-gated ion channel 1-like, had significant expression changes (S5 Table). Pistil extension-like protein was upregulated in L-morph flower pistils from development to maturity but showed no expression difference in S-morph flowers. In S-morph and L-morph flowers, suberization-associated anionic peroxidase 2, which may play an important role in cell wall suberization, superoxide elimination, and oxidative stress response [33,34], was upregulated. We also found that cinnamoyl-CoA reductase, which participates in lignin biosynthesis, was downregulated in both S-morph and L-morph flowers [35]. To test pistil protein expression levels of messenger RNA (mRNA) transcript differences between L-morph and S-morph flowers and validate the iTRAQ results, we selected nine genes for mRNA qRT-PCR analysis from DAB 0, 3, 6, 10, and 13 of both L-morph and S-morph pistils. As shown in Fig 4, except for Sme2.5_06391.1_g00003.1 and Sme2.5_13401.1_g00002.1, the expression of the other seven genes agreed well with our iTRAQ results from developing and mature flowers, suggesting that most proteins were regulated directly at the transcriptional level. The expression patterns of the two exceptions at the protein and mRNA levels were not consistent with the iTRAQ results during development but were consistent during maturity.  This may indicate that the abundance of these proteins depends not only on transcript levels but also on posttranslational modification.

GO annotation and KEGG pathway analysis of DEPs
We conducted GO enrichment and KEGG pathway analyses based on a hypergeometric test of all DEPs, with all proteins as the background. As shown in Fig 5, we itemized the respective top significant GO terms in Fig 5 during flower development and maturity. Fig 5A shows that significant biological processes during flower development concerned morphogenesis and metabolic processes. Significant biological processes varied greatly during flower maturity and included translation, gene expression, biosynthetic processes, and metabolic progress (Fig 5B). These findings are in agreement with current understanding of mature plant biochemical and physiological activities. For cellular components, terms related to the ribosome were enriched, suggesting that mature flowers had extensive and abundant translation, transcription, and regulation, more than in the developmental stage. We also conducted a GO enrichment analysis for DEPs of L-morph and S-morph flowers between developing and mature stages. In the L-morph flowers, significant biological processes included primary metabolic processes and carbohydrate metabolic processes (S4A Fig). However, for the S-morph flowers, the significant enriched biological processes were primarily metabolic processes, translation, and gene expression (S4B Fig).
The top 20 significant KEGG pathways are illustrated in Figs 6 and S4. Phenylpropanoid biosynthesis and metabolic pathways were enriched during flower development (Fig 6A). The enriched pathways during flower maturity were ribosome formation, starch and sucrose metabolism, cyanoamino acid metabolism, pentose and glucuronate interconversions, phenylpropanoid biosynthesis, flavonoid biosynthesis, and proteasome formation (Fig 6B). Enriched DEP pathways in L-morph flowers during maturity relative to development included metabolic pathways, molecular biosynthesis, ribosome formation, carbon fixation, and ribonucleic acid (RNA) polymerase biosynthesis (S5A Fig). Enriched DEP pathways in S-morph flowers during maturity relative to development included carbon fixation, ribosome formation, metabolic pathways, proteasome biosynthesis, molecular biosynthesis, and aminoacyl-transfer RNA biosynthesis (S5B Fig).

Protein-protein interaction analysis
Proteins in organisms do not act as single entities but rather form a variety of functional connections with each other, and these connections are fundamental in cellular processes. To uncover functional aspects associated with proteins in eggplant flowers, 225 proteins (Tables  1-4, with 16 overlapping DEPs between flower development and maturity) with significantly changed expression were analyzed by searching the STRING database. Fig 7 illustrates that six separate interaction networks were predicted. In each network, proteins that increased interacted with proteins that decreased to constitute a gene regulation network. These proteins were directly or indirectly related to pistil development differences of S-morph and L-morph flowers and were associated primarily with ribosome function, metabolic pathways, phenylpropanoid biosynthesis, starch and sucrose metabolism, and biosynthesis of secondary metabolites (Fig 7). Although these predicted interaction networks still need to be verified and further analyzed in future studies, the interactions between these proteins suggest that they have important roles in heterostylous development in eggplant.

Discussion
Heterostyly, which produces a distinctive flower polymorphism and results in hercogamy, is widely distributed in angiosperms [1,27,28]. Heterostyly had important implications in plant adaptability and yield because it affects breeding systems. Phylogenetic analysis has shown that this phenotype has originated independently among the heterostylous species and arisen via distinct evolutionary developmental pathways [3]. However, eggplant heterostyly has not been well studied. The study of regulated gene expression levels in producing different flower morphs should give insights into mechanisms involved in heterostyly development. We found that buds of both L-and S-morph lengthened linearly during development. However, around 10 DAB, the L-morph pistils continued to increase but the S-morph pistils generally did not elongate. The self-and cross-pollination of L-morph and S-morph flowers indicated that the structure of the stigmatic surface in S-morph flowers may inhibit pollen germination. A previous study [20] found that the pollination efficiency of L-morph flowers was the greatest, with a pollen germination percentage up to 100%. However, the percentage of pollen germination on the stigmas of S-morph flowers was much lower (< 5%). Another study found that the viability of pollen from L-morph and S-morph flowers exceeded 98% [36]. A third study indicated that there were differences in physiological and biochemical properties of stigmas between S-morph and L-morph flowers, and potassium levels, which may be influenced by auxin, were lower in S-morph flowers [37]. Rylski et al. [38] showed that the stigma of a short style differed from that of a long style in its smaller size, underdeveloped papillae, and lower sugar content. Structural observations of pistils may indicate that the low Protein-protein interaction network analyzed by STRING software. Network analysis results for significantly changed proteins between S-morph and L-morph flowers. The confidence score was set to ! 0.4 (medium). Different line colors represent the types of evidence for association. Known interactions: magenta = experimental evidence; light blue = database evidence. Predicted interactions: green = neighborhood evidence; red = fusion evidence; blue = co-occurrence evidence. Other: black = coexpression evidence; yellow = text-mining evidence; purple = protein homology.
https://doi.org/10.1371/journal.pone.0179018.g007 fertilization rates of S-morph flowers probably resulting from stigmatic incompatibility, rather than defects in the embryo sac. However, it is difficult to determine the mechanism underlying this difference in successful fruit set in S-morph and L-morph flowers.
iTRAQ is an advantageous technology used in comparative proteomic analysis that has provided greater reliability and accuracy for analysis in numerous studies [19,39]. In the present study, we used an iTRAQ approach to investigate the proteomic differences underlying heterostyly in eggplant. We identified 255 DEPs of S-morph and L-morph pistils during floral development and maturity. There were 33 upregulated and 24 downregulated genes during development, and 83 upregulated and 101 downregulated genes during maturity in S-morph compared with L-morph flowers (Tables 1-4). The results of qRT-PCR were highly consistent with our iTRAQ data. We characterized the differentially expressed proteins through GO, KEGG, and protein-protein interaction analyses. The protein expression differences were mainly involved in biosynthesis and metabolism, ribosomes, and expression regulation.

Proteins involved in biosynthesis and metabolism
During development, molecular biosynthesis and metabolism, including amino acid and fatty acid biosynthesis, are generally upregulated to meet needs for growth. In this proteomic study, 86 proteins involved in biosynthesis and metabolism were differentially expressed between Smorph and L-morph flowers (S6 and S7 Tables), providing evidence that differences in storage and energy availability can play key roles in pistil development. DEPs were involved in the biosynthesis of amino acids including valine, leucine, and isoleucine. Other DEPs were involved in metabolic pathways including glycometabolism, glycerolipid metabolism, and phenylalanine metabolism. Pectin methyl esterase was downregulated in S-morph flowers during maturity; this enzyme is involved in structural modifications of the cell wall during growth and development and is involved in plant-pathogen interactions [40]. Most DEPs (15/19, S7 Table) involved in starch and sucrose metabolism were upregulated in S-morph flowers during maturity. We may infer here that differences in metabolic pathways may play important roles in the development of pistil length. In addition, many proteins were involved in other biosynthetic pathways, such as those for phenylpropanoids, flavonoids, isoquinoline alkaloids, terpenoids, anthocyanins, and other secondary metabolites.

Proteins involved in ribosomes
The ribosome serves as the site of biological protein synthesis (translation), which links amino acids together in the order specified by messenger RNA molecules. Eukaryotes have 80S ribosomes, each consisting of a small (40S) and large (60S) subunit. In our study, the majority of both 60S and 40S ribosomal proteins was downregulated in S-morph flowers compared with L-morph flowers during maturity (Table 4). There were few differences in expression between the morphs during development (Tables 1 and 2), but the expression of ribosomal proteins was downregulated during maturity in S-morph flowers (Table 4). This may suggest that many proteins need to be translated during maturity as pollination and fruit development proceed.

Proteins involved in expression regulation
Many binding proteins play important roles in complex and intricate regulatory networks. Here, we identified 14 differentially expressed proteins related to nucleotide binding (RNA-, deoxyribonucleic acid-, and protein-binding factors) between S-morph and L-morph flowers during maturity (Table 2) based on their putative molecular functions; these binding factors might regulate protein activity in complex biological processes. For instance, the expression level of glycine-rich RNA-binding protein (GRP) was upregulated in S-morph flowers during maturity (1.56-fold, Table 3); this protein produces robust circadian oscillations and can suppress expression of the FLOWERING LOCUS C (FLC) protein, thereby promoting flowering [41,42]. FLC expression is directly promoted by histone H2A and leads to delayed flowering during vegetative growth [43]. In S-morph flowers at maturity, the histone H2A expression level was downregulated (1.58-fold, Table 4). Hence, we may infer that GRP and histone H2A tend to downregulate the expression of FLC to promote flowering in S-morph flowers. Additionally, the expression level of profilin-1 tended to be upregulated in S-morph flowers during development and maturity and was also upregulated in L-morph flowers during maturity. Pandey et al. [34] demonstrated that constitutive overexpression of cotton profilin-1 in tobacco induced early flowering [44]. Our results suggest that the expression changes of these binding proteins could affect blooming in S-morph flowers. Therefore, proteins involved in RNA-, deoxyribonucleic acid-, or protein-binding may be connected with heterostyly and regulate the pistil length of eggplant flowers.
Our study presents the first proteomic analysis of heterostyly development in eggplant. We identified approximately 225 DEPs in both developing and mature S-morph and L-morph flowers, whose expression levels might be closely related to heterostyly. Differences in flower development stages of the morphs were apparent in morphogenesis and metabolic processes. We discovered that some proteins associated with senescence and programmed cell death were upregulated in S-morph pistils, which may prevent them from developing into fruit. Our results provide valuable information on heterostyly in eggplant, and future studies of other heterostylous species may find similar mechanisms.

Conclusions
We performed self-and cross-pollination of L-morph and S-morph flowers and compared embryo sac development in eggplant, which suggested that S-morph flowers have stigma incompatibility features that inhibit pollen germination and subsequent fertilization. To explore the molecular mechanisms underlying heterostylous development, we conducted an iTRAQ-based proteomic analysis of eggplant pistils for L-morph and S-morph flowers. There were 225 DEPs in both developing and mature stages of S-morph and L-morph flowers whose expression levels might be closely related to heterostyly. We also conducted qRT-PCR for nine genes to verify DEPs from the iTRAQ data. Differences during flower development between the morphs were primarily observed in proteins related to morphogenesis and metabolic processes. Important biological processes, including translation, gene expression, biosynthetic processes, and metabolic progress, varied greatly during flower maturity. Additionally, we discovered that some proteins associated with programmed cell death were upregulated in Smorph pistils; these may be associated with low fruit set in these flowers. Our research provides important information for understanding eggplant heterostyly and establishes a foundation for the study of relevant mechanisms. It also highlights the importance of complex characters for understanding relationships between proteomic and phenotypic variation.