Global Quantitative Proteomics Studies Revealed Tissue-Preferential Expression and Phosphorylation of Regulatory Proteins in Arabidopsis

Organogenesis in plants occurs across all stages of the life cycle. Although previous studies have identified many genes as important for either vegetative or reproductive development at the RNA level, global information on translational and post-translational levels remains limited. In this study, six Arabidopsis stages/organs were analyzed using quantitative proteomics and phosphoproteomics, identifying 2187 non-redundant proteins and evidence for 1194 phosphoproteins. Compared to the expression observed in cauline leaves, the expression of 1445, 1644, and 1377 proteins showed greater than 1.5-fold alterations in stage 1–9 flowers, stage 10–12 flowers, and open flowers, respectively. Among these, 294 phosphoproteins with 472 phosphorylation sites were newly uncovered, including 275 phosphoproteins showing differential expression patterns, providing molecular markers and possible candidates for functional studies. Proteins encoded by genes preferentially expressed in anther (15), meiocyte (4), or pollen (15) were enriched in reproductive organs, and mutants of two anther-preferentially expressed proteins, acos5 and mee48, showed obviously reduced male fertility with abnormally organized pollen exine. In addition, more phosphorylated proteins were identified in reproductive stages (1149) than in the vegetative organs (995). The floral organ-preferential phosphorylation of GRP17, CDC2/CDKA.1, and ATSK11 was confirmed with western blot analysis. Moreover, phosphorylation levels of CDPK6 and MAPK6 and their interacting proteins were elevated in reproductive tissues. Overall, our study yielded extensive data on protein expression and phosphorylation at six stages/organs and provides an important resource for future studies investigating the regulatory mechanisms governing plant development.


Introduction
Unlike animals, which complete most of their organ differentiation during embryonic development, plants undergo a continuous developmental process, with organogenesis occurring throughout the life cycle. With the two apical meristems that form during embryogenesis, plants develop the

Quantitative Identification of 2187 Proteins and 1194 Phosphoproteins in Six Stages/Organs
To identify proteins that likely function at various Arabidopsis developmental stages, we collected six developmental stages/organs, including roots and above-ground parts (AGPs containing hypocotyl, cotyledons, and leaves) of 16-day-old plants, cauline leaves (CL) from 6-week-old plants, stage 1-9 floral buds (F1-9) including floral meristem, stage 10-12 floral buds (F10- 12), and open flowers (OF) ( Figure 1A-C). All plants for tissue sampling were grown in the same growth chamber with the same growth conditions, and for each tissue/organ, the sample for either proteomics or phosphor-proteomics was collected from more than 100 plants to reduce the variation due to individual plant differences. In addition, three replicates for each tissue were obtained, and global profiling of the quantitative proteome and phospho-proteome of each sample was performed following the workflow indicated in Figure 1D. between different plant stages/organs, enabling the formulation of hypotheses regarding protein functions and regulatory networks during plant development.

Quantitative Identification of 2187 Proteins and 1194 Phosphoproteins in Six Stages/Organs
To identify proteins that likely function at various Arabidopsis developmental stages, we collected six developmental stages/organs, including roots and above-ground parts (AGPs containing hypocotyl, cotyledons, and leaves) of 16-day-old plants, cauline leaves (CL) from 6-weekold plants, stage 1-9 floral buds (F1-9) including floral meristem, stage 10-12 floral buds (F10- 12), and open flowers (OF) ( Figure 1A-C). All plants for tissue sampling were grown in the same growth chamber with the same growth conditions, and for each tissue/organ, the sample for either proteomics or phosphor-proteomics was collected from more than 100 plants to reduce the variation due to individual plant differences. In addition, three replicates for each tissue were obtained, and global profiling of the quantitative proteome and phospho-proteome of each sample was performed following the workflow indicated in Figure 1D.  In total, 11,341 unique peptides corresponding to 2187 non-redundant proteins were identified with an estimated FDR (false discovery rate) of 1% (Figure 2A, Dataset 1), of which 2186 proteins with quantitative information were obtained. Pearson correlation was performed to determine the analytical reproducibility of the three replicates of each tissue ( Figure S1), and the results showed that all values of the correlation coefficients were approximately 0.8, indicating good correlation among different replicates. Int  In addition, quantitative phospho-protein analyses were performed for each of the six developmental stages/organs following the workflow indicated in Figure 1D. A total of 2497 unique phospho-peptides with 2893 detected phosphorylation sites in 1194 proteins were identified (Figure 2A, Dataset 2). Among these proteins, 926 phospho-proteins with 1475 phosphorylation sites have been reported in the Plant Protein Phosphorylation Database (P3DB, http://www.p3db.org) ( Figure 2B,C), suggesting that our data are relatively reliable. In addition, 268 phospho-proteins with 1022 phosphorylation sites were not included in the P3DB ( Figure 2B,C), indicating the phosphorylation modification of these proteins are newly identified and further suggesting a great depth of our detection.
We also analyzed the times of detection for each identified peptide in total three replicates for both proteomic and phospho-proteomic analysis. In the proteomic data, 47%, 25%, and 28% of peptides were identified three times, two times, and once, respectively ( Figure 2E). In addition, 40%, 28%, and 32% peptides were identified three times, two times, and once, respectively, in phospho-proteomic analysis ( Figure 2F). Furthermore, we analyzed the Peptide-Spectrum Matches (PSMs) in total proteome peptides and phospho-peptides. The PSMs of most peptides varied from 1 to 50, and the PSMs of phospho-peptides ranged from 1 to 80 ( Figure 2H).
To present the expression differences between the AGPs and roots and between the reproductive organs (including F1-9, F10-12, and OF) and the cauline leaves, we generated volcano plots for each pair in a comparison for proteomic and phospho-proteomic data. The differences in the expression levels of most peptides were more than 1.5-fold, with a p-value cutoff of 0.05 ( Figure 3). Compared to those in roots, most of the detected peptides and phospho-peptides were increased significantly in AGPs ( Figure 3A,B). Comparing to cauline leaves, approximately half of the detected peptides were up-regulated significantly in F1-9 ( Figure 3C), F10-12 ( Figure 3E), and OF ( Figure 3G), and almost all phospho-peptides were up-regulated significantly in F1-9 ( Figure 3D), F10-12 ( Figure 3F), and OF ( Figure 3H). These results suggest that the AGPs might require more proteins and phosphorylation modifications than roots, as do the reproductive organs/stages in comparison with the vegetative cauline leaves. Int  Volcano plot between OF and CL phosphopeptide samples. R, 14-day-old root; AGP, 14-day-old tissues above ground part; CL, cauline leaves; F1-9, stage 1-9 flowers; F10-12, stage 10-12 flowers; OF, opened flowers. The red denotes log2FC > 1 and p-value < 0.05. The blue denotes log2FC < −1 and p-value < 0.05. The black denotes −1 < log2FC < 1 or p-value > 0.05.

Sixteen-day AGPs and Roots have Distinct Protein Profiles
To investigate proteins that potentially have broad or distinctive functions in roots and/or AGPs, we first carried out protein expression comparisons between these two tissues. Interestingly, 1985 of the 2002 AGP-root differentially expressed proteins showed obviously increased expression in AGPs compared with roots (fold change of AGP/R, |FC| ≥ 1.5-fold), and only 17 proteins showed obviously reduced expression in AGPs ( Figure 4A,B). These results suggest that many more proteins are required for development and physiology of AGP than that for roots. Among the 1985 AGPpreferential proteins, 369 showed a 1.5-to 2-fold increase, 940 showed a 2-to 4-fold increase, and 676 showed a greater than 4-fold increase ( Figure 4B). Among the 17 root-preferential proteins, the fold change of seven proteins was over 4-fold ( Figure 4B). (H) Volcano plot between OF and CL phosphopeptide samples. R, 14-day-old root; AGP, 14-day-old tissues above ground part; CL, cauline leaves; F1-9, stage 1-9 flowers; F10-12, stage 10-12 flowers; OF, opened flowers. The red denotes log2FC > 1 and p-value < 0.05. The blue denotes log2FC < −1 and p-value < 0.05. The black denotes −1 < log2FC < 1 or p-value > 0.05.

Sixteen-Day AGPs and Roots Have Distinct Protein Profiles
To investigate proteins that potentially have broad or distinctive functions in roots and/or AGPs, we first carried out protein expression comparisons between these two tissues. Interestingly, 1985 of the 2002 AGP-root differentially expressed proteins showed obviously increased expression in AGPs compared with roots (fold change of AGP/R, |FC| ≥ 1.5-fold), and only 17 proteins showed obviously reduced expression in AGPs ( Figure 4A,B). These results suggest that many more proteins are required for development and physiology of AGP than that for roots. Among the 1985 AGP-preferential proteins, 369 showed a 1.5-to 2-fold increase, 940 showed a 2-to 4-fold increase, and 676 showed a greater than 4-fold increase ( Figure 4B). Among the 17 root-preferential proteins, the fold change of seven proteins was over 4-fold ( Figure 4B). roots, for short, 62/84), glucose metabolic processes (12/59), and the response to cadmium ions (38/182). In addition, proteins in abiotic stimulus processes (104/284), salt stress responses (32/113), and responses to cold (38/89) were enriched as well ( Figure 4C,D). These results are consistent with previous findings that light intensity, CO2 concentration, water availability, and temperature have combined effects on proliferation, cell expansion, and endoreduplication and subsequently affect shoot organ growth [62]. The molecular function (MF) analysis also showed that many of the 1985 proteins were involved in GTPase activity, translation factor activity, rRNA binding, and similar processes ( Figure 4E). However, the 17 root-preferential proteins were clearly enriched in response to oxidative stress and hydrogen peroxide ( Figure 4D). These 17 proteins included ( Figure S2): ZCE2, which is reported to function in plant bolting [63]; a class III peroxidase member RCI3, which has shown peroxidase activity in root cells [64] and is involved in reactive oxygen species (ROS) production [65]; IAA-CONJUGATE-RESISTANT4 (IAR4), which is critical in root development, root hair formation, and auxin response [66]; and JAL34 and PLAT2.  To investigate the functions of these differentially expressed proteins, we performed a gene ontology (GO) analysis of the 1985 AGP-preferential proteins and the 17 root-preferential proteins ( Figure 4C-E). Most of the 1985 proteins were enriched in basic cellular processes, including photosynthesis (62 of the 84 enriched proteins showed ≥4-fold increased expression compared with roots, for short, 62/84), glucose metabolic processes (12/59), and the response to cadmium ions (38/182). In addition, proteins in abiotic stimulus processes (104/284), salt stress responses (32/113), and responses to cold (38/89) were enriched as well ( Figure 4C,D). These results are consistent with previous findings that light intensity, CO 2 concentration, water availability, and temperature have combined effects on proliferation, cell expansion, and endoreduplication and subsequently affect shoot organ growth [62]. The molecular function (MF) analysis also showed that many of the 1985 proteins were involved in GTPase activity, translation factor activity, rRNA binding, and similar processes ( Figure 4E).
However, the 17 root-preferential proteins were clearly enriched in response to oxidative stress and hydrogen peroxide ( Figure 4D). These 17 proteins included ( Figure S2): ZCE2, which is reported to function in plant bolting [63]; a class III peroxidase member RCI3, which has shown peroxidase activity in root cells [64] and is involved in reactive oxygen species (ROS) production [65]; IAA-CONJUGATE-RESISTANT4 (IAR4), which is critical in root development, root hair formation, and auxin response [66]; and JAL34 and PLAT2.

Hierarchical Clustering Analysis of the Quantification of the Total Proteome
To identify possible correlations between protein expression level and plant development and putative organ/stage preferential functional proteins, we performed a cluster analysis (K_means = 15) on the 2186 proteins using RStudio and Euclidean distance classification. A total of 15 clusters with different expression patterns in the evaluated stages/organs were generated ( Figure S3). Among them, 89 proteins were expressed at similar levels in all six stages/organs ( Figure S3A). Gene ontology (GO) analysis of these 89 proteins showed that the enriched GO categories are response to cadmium ions, oxidation-reduction processes, response to salt stress, tricarboxylic acid cycle, malate metabolic process, calcium ion transmembrane transport, and similar processes ( Figure S4A). In total, 17 proteins showed significantly preferential expression in 16-day-old roots ( Figure S3) and were the same as the 17 root-preferential protein detected in the AGP/R comparison ( Figure 4 and Figure S2). In addition, 12, 40, 72, and 134 proteins were found to show significantly preferential expression in 16-day-old AGP, 6-week-old CL, F1-9, and F10-12, respectively ( Figures S3 and S4).
Interestingly, a total of 729 proteins showed obviously higher expression level in the three floral stages than in the other three tissues ( Figure S3H-J). The result that more proteins are relatively highly expressed in the reproductive organs/stages suggests that the normal reproductive process might need more complex functions of many proteins. In addition, proteins that were found to be preferentially expressed in certain organs/stages might be particularly important for the normal development and the physiology of the corresponding organ/stage.
To confirm the reliability of the quantitative proteomic data, especially for those protein clusters that showed some pattern(s) of stages/organ-preferential expression, we compared our proteomic data from six stages/organs to the transcriptomic data of the same or highly similar stages that are available in the eFP database (http://bar.utoronto.ca/efp/cgi-bin/efpWeb.cgi). As eFP has no transcriptomic data for stages exactly matching the 16-day AGP, stage 1-9 floral buds (F1-9), stage 10-12 floral buds (F10- 12), and open flowers (OF) as those in our proteomic study, we used transcriptomic data from rosette leaves, flower 9, flower 10/11 and 12, and flower 15, respectively. The results of the comparison showed a relatively high degree of correlation, as illustrated by the following clusters: (1) for the 17 proteins that showed root-preferential proteomic expression pattern, 15 (88%) showed a consistent pattern at the transcript level ( Figure 5A); (2) for the 40 cauline leaf (CL)-preferentially expressed proteins, transcript level of 40 (100%) in CL was considerably higher than that in roots and floral buds ( Figure 5B); (3) for the 72 early floral buds (EF, flower 1-9)-preferentially expressed proteins, 37 (51%) exhibited the highest transcript level in stage 9 flowers ( Figure 5C); (4) 43 of the genes encoding proteins that were identified as preferentially expressed in the late flower (LF, flower 10-12) were found to be preferentially expressed in flower 10/11 and/or flower 12 ( Figure 5D). Overall, the stage/organ preferential expression in transcript level further supported the proteomic results in our study.

Distinct Expression Profiles of Proteins in Three Floral Stages Compared to Cauline Leaf
To further investigate proteins preferentially expressed in reproductive development, we first compared the expression levels of proteins among three different floral stages and the cauline leaf.  Figure  6A). Next, we compared the protein expression data of these three flower developmental stages to uncover common and distinctive changes. Expression of 1212 of the 1766 proteins was altered in all F1-9, F10-12, and OF stages, 125 were found in both F1-9 and F10-12 stages, 34 were found in both F1-9 and OF stages, and 117 were found in both F10-12 and OF ( Figure 6B). In addition, 74,190, and 14 proteins, respectively, were found higher in F1-9, F10-12, and OF than cauline leaves ( Figure 6B). The 1212 proteins shared by all three floral stages were enriched in GO categories of both cellular processes and environmental responses. For cellular processes, ribosome biogenesis, ubiquitindependent protein catabolic process, formation of translation initiation complex, response to cytokinins, photosynthesis, fatty acid biosynthetic process, and microtubule-based processes were highly enriched. In addition, response to abiotic stresses including salt stress, cold, heat, and cadmium ion were strongly enriched ( Figure 6C).

Distinct Expression Profiles of Proteins in Three Floral Stages Compared to Cauline Leaf
To further investigate proteins preferentially expressed in reproductive development, we first compared the expression levels of proteins among three different floral stages and the cauline leaf.  Figure 6A). Next, we compared the protein expression data of these three flower developmental stages to uncover common and distinctive changes. Expression of 1212 of the 1766 proteins was altered in all F1-9, F10-12, and OF stages, 125 were found in both F1-9 and F10-12 stages, 34 were found in both F1-9 and OF stages, and 117 were found in both F10-12 and OF ( Figure 6B). In addition, 74,190, and 14 proteins, respectively, were found higher in F1-9, F10-12, and OF than cauline leaves ( Figure 6B). The 1212 proteins shared by all three floral stages were enriched in GO categories of both cellular processes and environmental responses. For cellular processes, ribosome biogenesis, ubiquitin-dependent protein catabolic process, formation of translation initiation complex, response to cytokinins, photosynthesis, fatty acid biosynthetic process, and microtubule-based processes were highly enriched. In addition, response to abiotic stresses including salt stress, cold, heat, and cadmium ion were strongly enriched ( Figure 6C). The p-value is the enrichment score (E-score) from DAVID. Higher −log10p-values indicate more significant enrichment. The dark and light red denote the higher and lower degree of GO enrichment respectively.

Anther-, Meiosis-, and Pollen-Preferentially Expressed Proteins were Identified
Angiosperm male reproductive development includes anther morphogenesis, meiosis, and pollen development. Previous studies have identified many genes preferentially expressed during these processes, such as anther-preferential expressed genes [67], meiosis-preferential genes [68], and pollen genes [67]; however, expression and function of proteins encoded by most of these genes have not been investigated.
To identify proteins with putative functions in anther development, meiosis, and pollen development, we compared our quantitative proteomic data with previously published antherpreferentially [67], meiosis-related [68], and pollen-preferentially expressed [67] gene clusters from transcriptomic profiling ( Figure 7A). We identified 15 anther-preferentially expressed genes that encode proteins among our proteomics results. Among the 15 proteins, 14 were more highly expressed in the three reproductive stages than in the cauline leaves, including 10 with the highest expression level in F10-12 (QRT3, KIN1, TKPR2, SCPL49, ACP1, AT3G23770, ACOS5, ATA27, MEE48/A6, and AT1G66850) ( Figure 7B). The other four proteins, RGP3, AGO9, NUDT3, and ACLB-2, showed similar expression levels in all three flower stages ( Figure 7B). Previous transcriptomic data showed that these genes were preferentially expressed in stage 4-7 anthers, which are part of stage 9 flowers [67]. In this study, our quantitative proteomics data indicated that their protein products might have been present for a considerably longer time and played roles even in the open flower when pollen development had completed and pollen grains were released from the anther. The high expression levels of these proteins during flower development strongly suggests their

Anther-, Meiosis-, and Pollen-Preferentially Expressed Proteins Were Identified
Angiosperm male reproductive development includes anther morphogenesis, meiosis, and pollen development. Previous studies have identified many genes preferentially expressed during these processes, such as anther-preferential expressed genes [67], meiosis-preferential genes [68], and pollen genes [67]; however, expression and function of proteins encoded by most of these genes have not been investigated.
To identify proteins with putative functions in anther development, meiosis, and pollen development, we compared our quantitative proteomic data with previously published anther-preferentially [67], meiosis-related [68], and pollen-preferentially expressed [67] gene clusters from transcriptomic profiling ( Figure 7A). We identified 15 anther-preferentially expressed genes that encode proteins among our proteomics results. Among the 15 proteins, 14 were more highly expressed in the three reproductive stages than in the cauline leaves, including 10 with the highest expression level in F10-12 (QRT3, KIN1, TKPR2, SCPL49, ACP1, AT3G23770, ACOS5, ATA27, MEE48/A6, and AT1G66850) ( Figure 7B). The other four proteins, RGP3, AGO9, NUDT3, and ACLB-2, showed similar expression levels in all three flower stages ( Figure 7B). Previous transcriptomic data showed that these genes were preferentially expressed in stage 4-7 anthers, which are part of stage 9 flowers [67]. In this study, our quantitative proteomics data indicated that their protein products might have been present for a considerably longer time and played roles even in the open flower when pollen development had completed and pollen grains were released from the anther. The high expression levels of these proteins during flower development strongly suggests their involvement in these processes. Previous studies based on phenotypic analyses suggest that QRT3 plays a direct role in degrading the pollen mother cell wall during microspore development [69], ACOS5 participates in a conserved biochemical pathway of sporopollenin monomer biosynthesis and is required for primexine formation [70,71], and MEE48/A6 is important for the dissolution of callose wall around the microspore tetrads [72,73]. In addition, our observation on another two T-DNA insertional mutant alleles, one ACOS5 allele acos5-4, and one MEE48 allele mee48-2 ( Figure 7C), also showed that, unlike the wild type plants that produced normal pollen grains with well-organized pollen exine structure, the acos5-4 mutant had shriveled pollens with abnormal pollen exine structure. Pollen defects of mee48-2 were stronger, most pollen grains were shrunken, and the pollen exine components were not even connected ( Figure 7D).

11
involvement in these processes. Previous studies based on phenotypic analyses suggest that QRT3 plays a direct role in degrading the pollen mother cell wall during microspore development [69], ACOS5 participates in a conserved biochemical pathway of sporopollenin monomer biosynthesis and is required for primexine formation [70,71], and MEE48/A6 is important for the dissolution of callose wall around the microspore tetrads [72,73]. In addition, our observation on another two T-DNA insertional mutant alleles, one ACOS5 allele acos5-4, and one MEE48 allele mee48-2 ( Figure 7C), also showed that, unlike the wild type plants that produced normal pollen grains with well-organized pollen exine structure, the acos5-4 mutant had shriveled pollens with abnormal pollen exine structure. Pollen defects of mee48-2 were stronger, most pollen grains were shrunken, and the pollen exine components were not even connected ( Figure 7D).
On the other hand, one of the 15 anther-preferentially expressed genes, BGLU44, was shown to be preferentially expressed in stage 4-7 anthers, but the BGLU44 protein showed nearly equal levels in all examined organ/tissues, except for its lower expression in roots ( Figure 5B), suggesting possible post-transcriptional regulation of BGLU44. In addition, proteins of four meiosis-and 15 pollenrelated genes were identified, and most of them were more strongly expressed in the reproductive organs than in the vegetative tissues, except for MPA1 ( Figure 7E,F).  On the other hand, one of the 15 anther-preferentially expressed genes, BGLU44, was shown to be preferentially expressed in stage 4-7 anthers, but the BGLU44 protein showed nearly equal levels in all examined organ/tissues, except for its lower expression in roots ( Figure 5B), suggesting possible post-transcriptional regulation of BGLU44. In addition, proteins of four meiosis-and 15 pollen-related genes were identified, and most of them were more strongly expressed in the reproductive organs than in the vegetative tissues, except for MPA1 ( Figure 7E,F).

Different Expression of Defense-Related Proteins in Different Organs
To identify the organ/stage-preferential expression of defense proteins, the 275 proteins that show clearly sage/organ-preferential expression pattern ( Figure S3B-H) were searched with keyword "defense". Nine defense related proteins were identified ( Figure S5 and Dataset 3). Among them, three proteins showed strong cauline leaf-preferential expression. These three proteins include EPITHIOSPECIFIER MODIFIER 1 (ESM1) that mediates indol-3-acetonitrile production from indol-3-ylmethyl glucosinolate for defenses against insect herbivores and various pathogens [74], CARBONIC ANHYDRASE 2 (CA2) that is structurally required for the assembly of Summary Complex I mitochondrial electron transport chain (mETC) and plays important role in reproductive development [75], and PHOTOSYSTEM II SUBUNIT P-1 (PSBP-1) that participates in the regulation of oxygen evolution and is involved in defense response to temperature and bacterium [76,77].
On the other hand, compared to their expression in vagetative tissues, the expressions of six proteins were clearly higher in the three flower developmental stages, especially F10-12 ( Figure S5 Figure S5 and Dataset 3). Among these proteins, ESP is suggested to be involved in pathogen resistance and leaf senescence [78][79][80]; TSA1 is suggested to be Jasmonic acid (JA) inducible and facilitates ERECTA (ER) body formation [81] and is involved in nuclear architecture [82] and seedling development in darkness [83].

Quantitative Phosphoproteomics Overview
Protein phosphorylation is one of the most important post-translational modifications regulating plant cellular functions. To investigate global protein phosphorylation patterns in major Arabidopsis organs/stages, we compared the phosphorylation levels of proteins among the six examined plant organs/stages, using quantitative proteomic analyses. To identify certain tissue-specific protein phosphorylation, we searched the phosphor-proteomic data for peptides being fully absent in one or more organs/stages. Among the 2497 detected phospho-peptides with 2893 detected phosphorylation sites in 1194 proteins, 51 peptides with 83 detected phosphorylation sites in 45 proteins were absent in root (R) but present in the other five organs/stages (AGP, CL, F1-9, F10-12, OF); two peptides with two detected phosphorylation sites in two proteins were absent in above-ground-parts (AGP) but present in the other five organs/stages (R, CL, F1-9, F10-12, OF); and two peptides with two detected phosphorylation sites in two proteins were absent in cauline leaves (CL) but present in the other five organs/stages (R, AGP, F1-9, F10-12, OF); another 2442 (97.8%) phospho-peptides with 2806 (97.0%) phosphorylation sites in 1145 (95.9%) proteins were present in all the six examined organs/stages (Supplemental Figure S6 and Dataset 4). These results suggest that most of the phosphorylation regulation likely exhibited a more-or-less but not yes-or-no differences among different organs/stages. Therefore, we further investigated the organ/stage-preferential phosphorylation. Between AGPs and root, we identified 995 phospho-proteins with 1956 phosphorylation sites that were differentially phosphorylated ( Figure S7A,B and Dataset 5). Similar to the previously mentioned greater number of proteins showing increased expression in AGPs, most of the phosphoproteins and the phosphorylation sites exhibited higher levels in AGPs than that in roots. In addition, 1149 proteins with 2313 phosphorylation sites showed different phosphorylation levels between OF and CL, with most of them showing obviously increased levels in OF ( Figure S7A,B and Dataset 6). These results further demonstrated regulation of protein function via phosphorylation likely plays a greater role in the AGPs than that in roots, and the normal development and physiology of reproductive stages/organs need functions of more protein and phosphorylation modification.
To uncover the potential biological processes in which these phosphoproteins are involved, we then examined predicted subcellular localizations of the differentially phosphorylated proteins between flower and leaf (OF/CL;1149) and between AGP and root (AGP/R; 995) using a public database (WoLF PSORT, http://www.genscript.com/wolf-psort.html). As shown in Figure S7C, 701 (61.0%) of the OF/CL differentially phosphorylated proteins and 603 (60.6%) of the AGP/root differentially phosphorylated proteins were predicted to localize to the nucleus ( Figure S7C), suggesting an important role of phosphorylation for nuclear-localized proteins, including transcription factors (TFs), in both vegetative and reproductive tissues. In addition, the numbers of phospho-proteins were similar in both the OF/CL and the AGP/R comparisons for localization to several subcellular locations or structures, including cytosol, cytoskeleton, extracellular space, mitochondria, chloroplast, and various membranes.
To obtain clues about regulation of TFs, we compared protein phosphorylation levels between the three floral developmental stages (F1-9 including floral meristem, F10-12, open flower), and the cauline leaves and examined the identified phosphoproteins for members of known TS families and their difference in phosphorylation levels. Among the 2296 Arabidopsis TFs (1717 loci) in 58 families in the PlantTFDB v4.0 (http://planttfdb.cbi.pku.edu.cn/), 62 TFs belonging to 22 families were identified here as phosphoproteins, including two ARF, three bHLH, six bZIP, eight C2H2, five GeBP, six MYB, two WRKY, and five ZF-HD family transcription factors ( Figure S7D). Furthermore, the numbers of identified phosphoproteins in each TF family were the same across all three floral stages, except those in the Trihelix family ( Figure S7D), suggesting that the extent of TF phosphorylation in the three floral developmental stages was similar.
We further compared the phosphorylation levels of the members of each TF family between F10-12 and F1-9 (F10-12/F1-9) and between OF and F10-12 (OF/F10-12). The results showed that only 11 TFs were differentially phosphorylated in F10-12/F1-9, including one bZIP, two C 2 H 2 , one C 3 H, two HB-other, one WRKY, and five ZF-HD family members, whereas 28 TFs belonging to 17 TF families were differentially phosphorylated in OF/F10-12, including the two ARF, two bZIP, three C 2 H 2 , one WRKY, and four ZF-HD family members ( Figure S7E). These results suggest that the phosphorylation modifications of TFs involved in developmental regulation between F1-9 and F10-12 are more similar than those between OF and F10-12. Moreover, 27 of the TFs showed decreased phosphorylation levels in F10-12/F1-9, but only one, ARF family member ARF8, showed an increased phosphorylation level. In OF/F10-12, 10 TFs showed decreased phosphorylation, and one ZF-HD family member, MIF2, increased ( Figure S7F).

Phosphorylation of Anther-, Meiosis-, and Pollen-Preferential Proteins
To uncover putative regulation by phosphorylation of protein encoded by previously reported anther-, meiosis-, and pollen-preferential genes, we compared our quantitative phosphorylation data to the anther-preferentially [67], meiosis-related [68], and pollen-preferentially expressed [67] gene clusters as well. Among all the detected phosphoproteins, one anther-preferential, three meiosis-related, and nine pollen-related gene encoded proteins were identified.
PolyA binding protein 5 (PAB5) was suggested to be anther-preferentially expressed, and in this study, we found a novel phosphorylation site in the second serine of the LASDLALpSPDK peptide of the PAB5 protein ( Table 1). The phosphorylation level of this serine was higher in stage 1-9 floral buds than that in cauline leaves (Log2FC(F1-9/CL) = 1.58). Later in flower development, its phosphorylation became higher in F10-12 and significantly decreased in OF ( Table 1). The phosphorylation level of the peptides from the meiosis-preferential protein MEI2-LIKE PROTEIN 5 (AML5) in F1-9 was higher than that in cauline leaves as well, then its phosphorylation decreased in F10-12 and OF. In comparison, two other meiosis-preferential proteins, SISTER-CROMATED COHESION PROTEIN 3 (SCC3) and ATM, showed dramatic increase in their phosphorylation levels in F1-9 compared with CL, and their phosphorylation levels remained high at later floral developmental stages ( Table 1). As meiosis occurs at anther stage 6 in flower stage 9, the relatively higher levels of phosphorylation modification level of the meiosis-related genes in F1-9 were consistent with their transcriptomic expression pattern and functional stage.
Nine phospho-proteins related to pollen were also identified, and their phosphorylation levels were up-regulated in F1-9 related to CL, with six out of these nine proteins having only one detected phosphorylation site and three proteins with 3-to-5 phosphorylation sites ( Figure 8A and Table 1). For instance, five, three, and four phospho-sites were identified in GLYCINE RICH PROTEIN 17 (GRP17), At4G28000, and At1G30470, respectively ( Figure 8A and Table 1). GRP17 encodes a glycine-rich protein and is expressed specifically during flower stages 10 to 12 [84], At4G28000 is involved in nematode-induced syncytium development and abiotic stress responses [85], and At1G30470 is a phosphoatase-associated family protein [86]. Remarkably, the phosphorylation levels of the peptide HpTpSGNDLHSR of At4G28000.1 protein were extremely high in the three floral stages with similiar phosphorylation levels in F1-9 and F10-12 (Log2FC(F1-9/CL) = 3.45, Log2FC(F10-12/F1-9) = 0.00) and even higher phosphorylation level in OF (Log2FC(OF/F10-12) = 3.45), suggesting potentially important function of this phosphorylation in flower and pollen development. In addition, the phosphorylation level of the detected serine site in the CMSGGMpSGSEGGMSR peptide of GRP17 showed progressive increases in F1-9 and F10-12 but dramatic reduction at OF stage, although its phosphorylation levels in these three floral stages were all obviously higher than those in the CL ( Figure 8A and Table 1).
In addition, several proteins encoded by pollen-related genes based on transcriptomic data [67] showed extremely high phosphorylated level in at least one floral organ. One of these proteins is AT1G52680.1, which is a late embryogenesis abundant protein-related/ late embriogenesis abundant (LEA) protein-like protein that is reported to be involved in flower development, pollen germination, and tube growth in Arabidopsis [87,88]. The phosphorylation level of the identified threonine (T) site in the NTLGMSPATNSPSSPAGpTTR peptide was much higher in the three floral stages than that in CL (Log2FC(F1-9/CL) = 2.41), remained the same level in F1-9 and F10-12, and increased to its highest level in OF (Log2FC(OF/ F10-12) = 1.37), suggesting a possible role for phosphorylation of this protein in flower developmental stages, especially in the pollen ( Figure 8A and Table 1). Another one (AT4G28000.1) is a member of the P-loop containing nucleoside triphosphate hydrolases superfamily. Three phosphorylation sites were detected, and two of them showed progressive increase in phosphorylation level during flower development and reached its highest phosphorylation level in OF ( Figure 8A and Table 1). Moreover, a member (AT1G30470.3) of the SIT4 phosphotase-associated family had four identified phosphorylation sites, and each of them showed highly reproductive-preferential phosphorylation pattern (Log2FC(F1-9/CL) = 3.63, 2.74, 2.7, 2.4, respectively), with three sites showing equal phosphorylation level in the three floral stages and one site showing obviously F1-9 abundant phosphorylation pattern ( Figure 8A and Table 1).  The asterisk represents the phosphorylated protein.

Protein Kinases and Phosphatases are More Likely to Be Phosphorylated in Reproductive Tissues
Protein kinases and phosphatases are often themselves regulated by phosphorylation. To obtain clues about possible functions of kinases and phosphatases in floral development, we compared protein and phosphorylation levels of kinases and phosphatases in the three floral stages, F1-9, F10-12, and OF, to those of the CL (Tables 2 and 3). Five kinases of four families were identified from the quantitative proteomics analysis, among which four kinases of CDPK, LRR_3, and SnRK2 families, CDPK6, RLK1, RLK902, and SNRK2-10, were up-regulated in floral stages ( Table 2).

Protein Kinases and Phosphatases Are More Likely to Be Phosphorylated in Reproductive Tissues
Protein kinases and phosphatases are often themselves regulated by phosphorylation. To obtain clues about possible functions of kinases and phosphatases in floral development, we compared protein and phosphorylation levels of kinases and phosphatases in the three floral stages, F1-9, F10-12, and OF, to those of the CL (Tables 2 and 3). Five kinases of four families were identified from the quantitative proteomics analysis, among which four kinases of CDPK, LRR_3, and SnRK2 families, CDPK6, RLK1, RLK902, and SNRK2-10, were up-regulated in floral stages (Table 2).   Note: Values mean Log2fold change/P-value. The fold change represents the ratio of the phosphorylation level of the detected peptide in a flower development stage to its phosphorylation level in cauline leaves. F1-9 means stage 1-9 flower buds, F10-12 means stage 10-12 floral buds, OF means opened flowers, and CL means cauline leaves. The domains in which each peptide is located are shown. PAS indicates the Per-ARNT-Sim (PAS) domain; S-TKc indicates the serine/threonine-protein-kinases catalytic domain; "-" indicates that the peptide was not in a specific domain.
In addition, 27 kinases/phosphatases belonging to 15 families showed obviously enhanced phosphorylation levels in at least one of the floral stages, as indicated by the quantitative phosphoproteomics results (Table 3). These results suggest that the functions of kinases and phosphatases in the floral stages are likely regulated by phosphorylation. The phosphorylated kinases include three MAP3K proteins, one MAP2K protein, three CDK proteins, four CDPK proteins, and four LRR proteins ( Figure 8B and Table 3). Among the phosphorylated kinases are the receptor kinase ERECTA (ER) and MPK6; previously ER and its closely related ER-like 1 and ER-like 2 (ERL1 and ERL2) together were shown to be important for early anther development, as were MPK6 and its close paralog MPK3 [22]. These results suggest that the functions of ER and MPK6 in anther development are regulated by phosphorylation.  Furthermore, we compared their phosphorylation levels of the kinases and the phosphatases among the three floral stages with F10-12/F1-9 and OF/F10-12 (Table 4). In total, 12 proteins showed changed phosphorylation level during flower development, including nine kinases with reduced levels in F10-12 as compared to those in F1-9 and three proteins (two kinases and one phosphatase) showing reduced phosphorylation in the OF related to F10-12. Notably, two phosphorylation sites were detected in the peptides of MPK6, ATSK11, STY8 (STY8-3), and MPK16 (MPK16-2), respectively (Table 3), and both the phosphorylation level and the change pattern of the two sites in each peptide during floral development were found to be the same. These findings of the 27 floral-preferential phospho-proteins suggest important functions and regulation by phosphorylation in flower development. Moreover, almost half (12/27) of them showed much higher phosphorylation levels in F1-9, suggesting that early floral development likely involves more protein phosphorylation compared with later stages.
To further test for changes of the level of phosphoproteins, we scanned the entire list of detected phospho-proteins here for ones that have available antibodies; we found such antibodies for three detected phosphoproteins, GRP17, CDC2/CDKA.1, and ATSK11. We obtained the commercially available antibodies for western blot analysis. As shown in Figure 8C, only one GRP17 protein band was detected by the anti-GRP17 antibody in the protein gel in each examined organ/stage, whereas in the phos-tag gel, another band was present in both F1-9 and F10-12 stages, with relatively weak and much stronger signals, respectively ( Figure 8C). This result confirmed the phosphorylation of GRP17 and its F10-12-preferential phosphorylation pattern ( Figure 8A and Table 1). In addition, similar analyses of CDC2/CDKA.1 and ATSK11 also uncovered bands of reduced mobility of phospho-proteins in flower organs but not in cauline leaves and rosette leaves, further validating the quantitative phospho-proteomics results.

Phosphorylation Regulatory Network during Floral Development
To investigate the possible regulatory network mediated by protein phosphorylation, we predicted kinase-substrate relationships according to the results from motif analysis and the PhosPhAt database [42]. The analysis supported 71 predicted kinase-substrate pairs and their interaction types from preferentially phosphorylated proteins from F1-9/CL, 22 in F10-12/F1-9, and two in OF/F10-12 comparisons (Dataset 7). We found that the putative substrates of MAPK6 were highly enriched in both F1-9/CL and F10-12/F1-9 comparisons, as part of the interaction network illustrated using Cytoscape 3.4.0 ( Figure 9). Our results are consistent with the previous conclusion that MAPK6 plays an important role in flowering, especially in anther development [22]. Among the putative substrates of MAPK, one transcription factor, bZIP16, was identified as differentially expressed during floral development. bZIP16 has demonstrated that it can integrate light and hormone signaling pathways to regulate early seedling development and can heterodimerize with other G group members [89]. Our phosphorylation results on bZIP16 suggest an important role of phosphorylation in regulating bZIP16 activity during floral development.
Among other possible interactive kinases and substrates, components of the JA and the Abscisic Acid (ABA) signal regulatory network were up-regulated in flowers 1-9 when compared with the leaf. We also identified BSK1, which is one of three homologous Brassinosteroid (BR)-signaling kinases belonging to the RLCK_2 family and is autophosphorylated as part of BR signaling and ROS pathways [90]. Another identified kinase is KIN10, an SnRK1 family kinase, which is autophosphorylated during calcium signaling of plant growth and phosphorylates TPS5 and NIA2 during sugar signaling. A member of the Arabidopsis CDPK gene family, CDPK6, also had target proteins overrepresented in our network during floral stages 1-9 and 10-12, suggesting that CDPK6 may play important roles in floral development. The substrates of YDA and the LRR STY8 were also significantly identified in comparisons of both floral stages 10-12 with 1-9 and floral 1-9 to leaf.

22
pathways [90]. Another identified kinase is KIN10, an SnRK1 family kinase, which is autophosphorylated during calcium signaling of plant growth and phosphorylates TPS5 and NIA2 during sugar signaling. A member of the Arabidopsis CDPK gene family, CDPK6, also had target proteins overrepresented in our network during floral stages 1-9 and 10-12, suggesting that CDPK6 may play important roles in floral development. The substrates of YDA and the LRR STY8 were also significantly identified in comparisons of both floral stages 10-12 with 1-9 and floral 1-9 to leaf.

Discussion
In flowering plants, vegetative processes and reproductive development both involve crucial pathways that directly determine the growth and other functions of the plant. Using iTRAQ-based quantitative proteomic and phosphoproteomic analyses, this study provides an overview of dynamic protein expression and phosphorylation during six Arabidopsis developmental stages at a whole-

Discussion
In flowering plants, vegetative processes and reproductive development both involve crucial pathways that directly determine the growth and other functions of the plant. Using iTRAQ-based quantitative proteomic and phosphoproteomic analyses, this study provides an overview of dynamic protein expression and phosphorylation during six Arabidopsis developmental stages at a whole-genome scale. In addition, comparative analyses of protein expression and phosphorylation modification reveal the distinct regulation of gene function at both translational and post-translational levels between vegetative and reproductive developmental processes and between different floral developmental stages.
Our data reveal that the photosynthetic process was enriched in rapidly growing organs, AGP in vegetative processes, and flowers in reproductive processes. More detailed analyses showed that most of the proteins related to photosynthesis were involved in light reactions, harvesting, and stimulus responses. Chlorophyll A/B binding proteins such as CAB3, LHCA1/2/3, and LHCB2/3/6 were identified in our differentially expressed proteins (DEPs) with stronger gene expression levels in leaves than in any other organs or tissues [91][92][93]. However, there are many differences between the vegetative and the reproductive phases. Translation and metabolic processes are very active in reproductive organs, whereas responses to different stresses and oxidation-reduction are more significant in seedlings. From the perspective of phospho-proteomics analysis, there are more phosphorylation events during reproductive development than in the vegetative process. The phospho-proteins with changes in phosphorylation levels showed more distribution in the nucleus during the reproductive process. Transcription factors (TFs) are among the most important elements in the nucleus. There is only one difference in the distribution of TF families between the two processes, namely, the expression of the Trihelix families.
Our proteomics and phospho-proteomics analyses focused on comparisons between vegetative and reproductive processes and among the reproductive phases from floral stages 1-9 to open flowers. Hundreds of phospho-proteins were found to be phosphorylated and dephosphorylated during flowering, especially transcription factors and kinases, which suggests their possible molecular mechanisms in regulating floral development. There were 27 transcription factors phosphorylated in the early floral stages compared to vegetative organs but dephosphorylated in floral stages 10-12. Another 10 TFs were dephosphorylated after flowers opened. In addition to transcription factors, kinases play significant roles in these processes. Nearly one third of the kinases showing high phosphorylation levels in floral stages 1-9 were dephosphorylated in stages 10-12. The phosphorylation networks between kinases and their targets also indicated that there are multiple hormone signaling pathways involved in the flowering process, including JA, ABA, BR, and ROS.
Proteomics has been widely used to generate comprehensive and quantitative map. Meanwhile, phospho-proteomics provides insights into the dynamic regulation of proteins by post-translational modifications (PTMs), which are critical for protein function. However, some questions resulting from iTRAQ-based proteomics remain to be resolved. Although our iTRAQ-based data provided protein abundance information for different organs, the proportions of specific proteins in various tissues are still unknown. Other questions are how the threshold value of differentially expressed proteins can be standardized, and when it is appropriate to use 1.2, 1.5, or 2. The use of a lower ratio such as 1.2 would include more quantification variations; however, there would be more false positive variations. In comparison, a higher ratio such as 2 might enhance the reliability of detected variations, but it also could possibly overlook useful information with small changes. Therefore, the key issue of selection of candidates should base on two aspects; one is the confidence level of the quantitative data, the correlation of biological triplicate, for example, and the other is how willing the scientist is to spend extra effort to validate the selected candidates with an alternative method including immunoblot, multiple reaction monitoring (MRM), or selected reaction monitoring (SRM) of samples.
Thus, iTRAQ-based proteomics generates relative expression levels, but it has a number of disadvantages. For instance, the iTRAQ method provides both a lower number and a lower percentage of differentially abundant proteins than the label-free method. At present, label-free quantitative methods are becoming more popular in proteomics and other biological studies. However, when time on the instrument or cost are important considerations, iTRAQ may be the method of choice. Increasingly, efforts have been made to improve proteome sample preparation and mass spectrometry instruments. In the future, a greater number of scientific issues will be investigated using proteomics and phospho-proteomics, allowing the identification of key genes, signaling pathways, or regulatory mechanisms.

Plant Materials and Experiment Design
Arabidopsis thaliana plants of the Columbia ecotype (Col-0) were used in this study. Seeds were surface-sterilized with 75% EtOH for 3 min, rinsed 3 times with sterile water, plated on plates containing half-strengthen Murashige and Skoog (1/2 MS) medium (Wako) with 1.5% agar and 1% sucrose (pH 5.8), and then stratified at 4 • C in the dark for 2 days. Plants were germinated and grown under long day conditions (16 h light/8 h dark; light intensity: 120 µmol m −2 s −1 ; humidity: 60-70%) at 22 • C until the roots (R) and the above-ground parts (AGP) of 16-day plants were harvested. The 16-day-old Col-0 plants were transferred to the soil and grown for another 4 weeks under the same light and temperature conditions, then cauline leaves (CL), stage 1 to 9 flowers (F1-9), stage 10 to 12 flowers (F10- 12), and open flowers (OF) were collected and immediately frozen in liquid nitrogen, respectively, and stored at −80 • C until further analyses. For each tissue/organ in this study, the sample for either proteomics or phospho-proteomics was collected from more than 100 plants.

Protein Extraction and Peptide Preparation
Two grams of tissues from each stage/organ were ground to fine powder in liquid nitrogen and suspended with buffer containing 6 M Urea, 2 M Thiourea, and 100 mM NH4HCO3 (final concentration), respectively. Suspensions were sonicated for 60 min (2 s sonication with 5 s intervals) on ice, and the extracted proteins were collected by centrifugation at 20,000× g for 20 min at 4 • C. Then, proteins were reduced with 10 mM DL-Dithiothreitol (DTT) for 1 h at 37 • C and alkylated with 30 mM iodoacetamide for 1 h at room temperature (25 • C) in the dark. The excess iodoacetamide was quenched by the addition of 10 mM DTT. After another centrifugation at 30,000× g for 20 min at 4 • C, the protein concentration was determined using Quick Start Bradford Protein assay (Biorad, Hercules, CA, USA).
An aliquot of total protein (200 µg) was taken from each sample and digested with endoproteinase Lys-C at the ratio of 1:100 (Lys-C: protein) for 4 h at room temperature, followed by digestion with Trypsin (Promega, Madison, WI, USA) at the ratio of 20:1 (protein: trypsin) at 37 • C for 12 h.

iTRAQ Labeling and Desalting
After trypsin digestion, the peptides were dried by vacuum centrifugation, and 100 µg peptides for each tissue were processed according to the manufacturer's protocol for 8-plex iTRAQ (Applied Biosystems). The six samples were labeled with the six iTRAQ-tags: 115-iTRAQ tag for F1-9, 116 for F10-12, 117 for OF, 118 for R, 119 for AGP, and 121 for CL.
Then, the iTRAQ-labeled peptides were mixed together, and 100 µg peptide mixtures were acidified with Formic Acid (FA) and loaded onto a pre-equilibrated homemade Poros R3 microcolumn. After washing the R3 resin twice with 5% FA, the peptides were eluted by 20 µL of 30% Acetonitrile (ACN), followed by 20 µL of 60% ACN.

Nanoflow LC-ESI-MS/MS Analysis Based on Linear Trap Quadrupole (LTQ)-Orbitrap Elite
Liquid chromatography (LC)-electrospray ionization (ESI) tandem MS (MS/MS) analysis was performed using a nanoflow EASY-nLC 1000 system coupled to an LTQ-Orbitrap Elite mass spectrometer (Thermo Fisher Scientific) following previously reported procedures (Zhang et al., 2015). A two-column system was adopted for all analyses. Each peptide fraction was first loaded onto an Acclaim prepMap100 C18 Nano Trap Column and then analyzed on an Acclaim Prep Map RSLC C18 Column. The peptides were eluted using the following gradients: 5-35% solution B (0.1% FA in 95% ACN) in 58 min, 35-90% B in 5 min, 90% B for 13 min at a flow rate of 200 nL/min.
The MS/MS mode was set as follows: activation type, higher collision energy dissociation (HCD); collision energy, 35 eV; resolution, 60,000 FWHM; scan range, 300-1800 m/z; The peptide ions were detected in the Orbitrap mass spectrometer, and up to 15 of the most intense peptide ions (>5000 counts) were selected and fragmented by MS/MS. Information on peptides and peptide fragments m/z was collected as follows: 15 fragment spectra were collected after every full scan (MS2 scan), HCD fragmentation, full scan at a resolution of 60,000, and normalized collision energy of 35 eV. Three technical replicates were conducted to ensure reliable statistical consistency.

Proteomic Data Analysis
Raw data files acquired from the LTQ-Orbitrap Elite were converted into MGF files using Proteome Discoverer 1.4 (Thermo Fisher Scientific, Germany), and the MGF files were queried against the Arabidopsis database (ftp://ftp.arabidopsis.org/home/tair/Sequences/blast_datasets/TAIR10_blastsets/ TAIR10_pep_20110103_representative_gene_model_updated) using an in-house Mascot search engine 2.3.02 (Matrix Science, London, UK).
Data were searched using the following parameters: 10 ppm mass tolerance for MS and 0.05 Da for MS/MS fragmented ions, with an allowance for only tryptic peptides with up to two missed cleavage sites allowed. Only peptides with a significance score greater than 99% confidence interval by a Mascot probability analysis were counted as identified. The quantitative protein ratios were weighted and normalized by the Mascot search engine (Matrix Science), and only ratios with p < 0.01, as determined according to a Student's t-test, were employed.
The summation of six labeled sample mixtures was used as background (BC) based on the weighted average of the intensity of report ions in each identified peptide. The relative expression of protein was normalized by the ratio of value in certain samples to that of the summation (F1-9/BC, F10-12/BC, OF/BC, R/BC, AGP/BC, CL/BC). Proteins that were once identified in three technical replicates were all considered. In the vegetative process, the final ratio of protein was obtained from the comparison AGP/BC to CL/BC. At the reproductive stage, the ratios of CL/BC were used as reference (REF), and all the three flower organs, F1-9, F10-12, and OF, were compared to REF, among which ratios greater than ±1.5-fold were considered to be differentially expressed proteins (DEPs).

Phosphoproteomic Data Analysis
As with the proteome data, raw LC-MS/MS files from phoshoproteomics experiments were queried against the Arabidopsis database using an in-house processed using Mascot search engine 2.3.02 (Matrix Science, London, UK), with the same parameters except for four significant points: 1. phospho_STY (serine, threonine, and tyrosine) was added in the specified parameters in the protein database searches; 2. identified peptides were further validated with a Target Decoy PSM validator, and phosphorylation sites were evaluated with phosphoRS3.0; 3. a false discovery rate (FDR) of 1.0% was used for the identification of phosphorylation residues (phosphosites); 4. regarding the phosphorylation levels of phosphoproteins, we determined the ratios of each phosphosites manually instead of using the calculated results from the Proteome Discover.

Bioinformatic Analysis
Gene ontology (GO) analysis of differentially expression genes was performed using online DAVID analysis (http://david.abcc.ncifcrf.gov/). Pathway enrichment analysis was performed using the list of genes that encode for proteins containing the peptide/phosphopeptides from each of the inferred cluster. Pathway enrichment with a set of genes was evaluated by comparing that set of genes against genes within known pathways using Mapman software. The subcellular localization of proteins and phosphoproteins was predicted with WoLF PSORT (http://above-ground.genscript. com/wolf-psort.html). R language (x64 3.3.1) was used to cluster the DEPs and to draw heatmaps. The relationship of detected kinases and potential substrates was generated through PhosPhATDB (http://phosphat.uni-hohenheim.de/) [42].