Transcriptomic,proteomic,and phosphoproteomic analyses reveal dynamic signaling networks influencing long-grain rice development

Fngyu Chen,Yongsheng Wng,Zesen Zhng,Xiolong Chen,2,Jinpeng Hung,Zhiming Chen,Jingsheng Zheng,Lingrong Jing,Yumin Hung,Houong Wng,Rongyu Hung,*

a Key Laboratory of Ministry of Education for Genetic Improvement and Comprehensive Utilization of Crops,Fujian Provincial Key Laboratory of Crop Breeding by Design,College of Agriculture,Fujian Agriculture and Forestry University,Fuzhou 350002,Fujian,China

b Postdoctoral Station of Biology,School of Life Sciences,Hebei University,Baoding 071000,Hebei,China

c School of Life Sciences,Xiamen University,Xiamen 361102,Fujian,China

Keywords:Proteome Phosphoproteome Transcriptome LGS1/GS2/GL2/OsGRF4 Young panicle Rice (Oryza sativa L.)Grain size

ABSTRACT The LGS1(Large grain size 1)gene,also known as GS2/GL2/OsGRF4,is involved in regulating grain size and quality in rice,but the mechanism governing grain size has not been elucidated.We performed transcriptomic,proteomic,and phosphoproteomic analyses of young rice panicles in Samba (a wild-type cultivar with extra-small grain)and NIL-LGS1(a nearly isogenic line of LGS1 with large grain in the Samba genetic background)at three developmental stages(4–6)to identify internal dynamic functional networks determining grain size that are mediated by LGS1.Differentially expressed proteins formed seven highly functionally correlated clusters.The concordant regulation of multiple functional clusters may be key features of the development of grain length in rice.In stage 5,16 and 24 phosphorylated proteins were significantly up-regulated and down-regulated,and dynamic phosphorylation events may play accessory roles in determining rice grain size by participating in protein–protein interaction networks.Transcriptomic analysis in stage 5 showed that differentially expressed alternative splicing events and dynamic gene regulatory networks based on 39 transcription factors and their highly correlated target genes might contribute to rice grain development.Integrative multilevel omics analysis suggested that the regulatory network at the transcriptional and posttranscriptional levels could be directly manifested at the translational level,and this analysis also suggested a regulatory mechanism,regulation of protein translation levels,in the biological process that extends from transcript to protein to the development of grain.Functional analysis suggested that biological processes including MAPK signaling,calcium signaling,cell proliferation,cell wall,energy metabolism,hormone pathway,and ubiquitin-proteasome pathway might be involved in LGS1-mediated regulation of grain length.Thus,LGS1-mediated regulation of grain size is affected by dynamic transcriptional,posttranscriptional,translational and posttranslational changes.

1.Introduction

Rice (Oryza sativa L.) is one of the most important cereal crops whose seed is the main source of calories for more than half of the global population [1].With the continuous increase of the world’s population and the decrease in arable land,it is a challenge to increase rice production to address the predicted food shortage[2].Rice grain size,which is characterized by grain length,grain width,and grain thickness,affects grain weight and,in turn,yield[3].Grain size also influences the market value of rice grain products.Rice grain with a long and slender shape,such as indica rice,is generally preferred by consumers in southern China,the USA,and south and southeast Asian countries,although the preference for rice grain characteristics varies with consumer groups [4].Increased grain size is a target of rice breeding.

Several major quantitative trait loci (QTL) regulating grain size have been cloned,helping to reveal the molecular basis of the regulation of grain size [5,6].Among these QTL,LGS1 (Large grain size 1),also known as GS2/GL2/OsGRF4 [6–9],was cloned,and the molecular mechanism governing the function of LGS1 was investigated in our previous work[10].LGS1/GS2/GL2 encodes the OsGRF4 transcription factor,which belongs to the growth-regulating factor(GRF) family that is targeted by microRNA miR396 [6–10].LGS1 mutant alleles contain a 2-bp missense mutation in the miR396 targeting sequence,and the associated increased expression leads to larger and longer grain than in the wild-type line [10].The LGS1 mutation also increased cold tolerance at the seedling stage[10],suggesting that LGS1 might function in stress resistance as well as grain size.With its multiple advantages,LGS1/GS2/GL2/OsGRF4 could be useful in rice breeding,especially of indica rice.

LGS1/GS2/GL2/OsGRF4 is involved in brassinolide(BR)responses[7,11]and is preferentially expressed in young panicles,especially in grain husks (hulls),which are the outer cover of the grain and determine grain shape and size.Thus,young panicle development influences the yield and milling quality of the rice grain by determining volume,shape,and size in the LGS1 mutant.However,the mechanism regulating the size of rice grain via LGS1/GS2/GL2/OsGRF4 has not been elucidated.

Novel omics methods have been used to investigate the regulatory mechanism governing the development process of rice seed at various levels[12–14].Transcriptional profiling and network analyses can reveal molecular pathways and signaling networks in rice development and stress response.However,transcript levels are insufficient to predict protein levels in many scenarios,especially during the dynamic developmental stage.Because mRNA expression levels may not translate precisely into protein function,mRNA profiling may not fully characterize the functional proteome.Accordingly,proteomics has been introduced as a powerful approach to study the function and regulation of genes.As a validated,high-throughput,quantitative functional proteomics platform,the isobaric tags for relative and absolute quantification(iTRAQ) method has been developed for the simultaneous quantitative comparison and analysis of the protein expression profiles of multiple samples [15].The iTRAQ-based proteomic technique is considered[16,17]to be one of the most robust techniques for differential quantitative proteomic analyses.With recent advances in mass spectrometry-based analytical technologies [18],deep proteomic profiling with extensive coverage and throughput provides a useful means of comprehensively characterizing proteome dynamics during plant development.As a result of the recent development of novel methods in phosphopeptide enrichment and mass spectrometry,high-throughput identification of phosphopeptides and phosphorylation sites at the proteome level has become available.

In this study,we present a global analysis using proteomics,phosphoproteomics,and RNA sequencing data collected from young panicles in the extra-small-grain cultivar Samba(wild type)and NIL-LGS1(a near-isogenic line of LGS1)with large grains in several developmental stages to investigate the role played by LGS1 in the mechanism governing rice grain size.

2.Materials and methods

2.1.Plant material

NIL-LGS1 carries the LGS1 allele from the extra-large-grain cultivar JF178(donor parent)in the Samba(recurrent parent)genetic background.The extra-small-gain cultivar Samba carries the lgs1 allele and serves as the wild type.

2.2.Growing conditions and sample collection

Rice plants were grown under natural conditions in the experimental field of Xiamen University,Fujian province,China.Seeds were sown in March 2017 and the seedlings were transplanted at an interplant spacing of 20 20 cm.When the plants were close to transitioning from the vegetative to the reproductive stage,panicle initiation was monitored daily by destructive sampling and observation with binoculars.The developmental stages of young panicles were distinguished following Ding [19],and fresh young panicle samples at three developmental stages (stages 4–6) were harvested in triplicate,frozen immediately in liquid nitrogen,and held at -80 °C until use.The lengths of the young panicles in the three stages were 2.0,5.0,and 15.0 cm.Because our previous work[10] indicated that expression differences between NIL-LGS1 and Samba of LGS1 at the transcript level were highly significant at developmental stage 5 in young panicles,young panicle samples at stage 5 and the two adjacent stages 4 and 6 were collected to investigate dynamic changes.

2.3.Total protein extraction,digestion,and iTRAQ labeling

Protein extraction and digestion of each young panicle sample for proteomic analysis followed Han et al.[20]with minor modification.As shown in Fig.1B,desalination of digested protein samples was performed on Sep-Pak C18 columns (Waters,Milford,MA,USA),and the samples were washed with 300 μL Milli-Q water and eluted with 50 μL methyl alcohol.The elution was dried under vacuum and stored at -80 °C for further use.iTRAQ labeling of peptides was performed with an iTRAQ 4-plex kit(Applied Biosystems,Foster City,CA,USA)according to the manufacturer’s instructions,followed by drying by vacuum centrifugation and storage at-80°C until phosphopeptide enrichment or LC-MS/MS analysis.

2.4.Phosphopeptide enrichment

Polymer-based metal ion affinity capture (PolyMAC) was performed to enrich phosphopeptides in the stage-5 samples from Samba and NIL-LGS1 and to remove nonspecifically retained nonphosphopeptides.As shown in Fig.1B,two hundred micrograms of desalted,iTRAQ-labeled,and dry peptide samples (consisting of 100 μg Samba proteins pooled with 100 μg NIL-LGS1 proteins)was completely resuspended in 40 μL Milli-Q water.Next,the detailed procedure for phosphopeptide enrichment was performed using a One-Step Magnetic PolyMAC-Ti Phosphopeptide Enrichment Kit (Tymora Analytical,West Lafayette,IN,USA) according to the manufacturer’s instructions.Eluates were dried completely in a vacuum centrifuge and stored at–80 °C for MS analysis.

2.5.LC-MS/MS

The LC-MS/MS and information acquisition of the peptides and phosphopeptides were performed on a high resolution mass spectrometer TripleTOF 5600 MS (AB SCIEX,Ontario,CA,USA) following Zhong et al.[21].

2.6.Bioinformatics

Fig.1.Summary of phenotype and data types included in this study.(A)Grain phenotype of Samba and NIL-LGS1.Scale bar,2 mm.ns,not significant.*,**,and ***indicate significant at P ＜0.05,P ＜0.01,and P ＜0.001,respectively.(B)Experimental strategy for whole proteome,phosphoproteome,and transcriptome analysis.Scale bars,2.5 cm.(C) Preliminary analysis of the omics data.The left and middle bar plots show the distribution of peptide length from proteins with or without phosphorylation.The right heatmap plot shows consistent clustering across biological replicates from Samba and NIL-LGS1.Samba_1-3 and NIL-LGS1_1-3 indicate biological replicates from Samba and NIL-LGS1,respectively.

Differentially expressed(DE)proteins,phosphorylated proteins,and genes in Samba and NIL-LGS1 were detected using a combination parameter (P-value ＜0.05 and fold change (FC) ＞1.5) calculated with the DESeq R package [22].Principal component analysis (PCA) performed on whole-transcript levels utilizing the pheatmap package in R[23].Functional clusters containing DE proteins with similar expression patterns were identified with Short Time-series Expression Miner (STEM) [24].GO (Gene Ontology)[25] enrichment analysis of DE proteins,phosphorylated proteins,and genes was performed with the BiNGO plugin of Cytoscape[26] with FDR=5%.The interaction network of proteins in rice was downloaded from the STRING database (http://string-db.org/),and functional clusters of the network based on phosphorylated proteins were identified with the MCODE[27]plugin of Cytoscape with default parameters.To detect alternative splicing events identified in a previous study[10],raw reads were first mapped to the genome of indica rice (ASM465 v1) [28],and the software rMATs[29] was employed to identify four alternative splicing events.TFs were detected by homologous sequence alignment using the quality annotated non-redundant protein sequence database of NCBI [30].Subnetworks of complex GRN were selected with the CytoHubba[31]plugin of Cytoscape using the Maximal Clique Centrality(MCC) method.

3.Results

3.1.Phenotyping and profiling of Samba and NIL-LGS1

In the three spikelet-formation stages,Samba and NIL-LGS1 differed in grain length,grain width,and thousand grain weight(Fig.1A).Of 4452 proteins identified,167 were phosphorylated(Table S1).The lengths of peptides,as determined by proteomics,ranged mostly from 4 to 12 aa,whereas the lengths determined by phosphoproteomics ranged mostly from 10 to 20 aa(Fig.1C,left and middle).Principal component analysis (PCA) performed on whole-transcript levels showed that three biological replicates of NIL-LGS1 and Samba were highly correlated and clustered(Fig.1C,right).

3.2.Reprogrammed proteome and multiple functional clusters

The differentially expressed proteins identification showed that 513(289 up-and 224 down-regulated),170(99 up-and 71 downregulated),and 498 (209 up-and 289 down-regulated) proteins were differentially expressed between Samba and NIL-LGS1 during stages 4–6,respectively(Fig.2A;Table S2).Seven highly correlated clusters formed by the majority of DE proteins (842) were distinguished and named clusters 26,39,41,43,44,45,and 47 (Fig.2B;Table S2).Cluster 26 containing 48 DE proteins,exhibited upregulation,unchanged,and down-regulation patterns in NIL-LGS1 during stages 4,5,and 6,respectively,suggesting that the response of these DE proteins to LGS1 might be stage-specific.GO enrichment analysis suggested that LGS1-mediated regulation of grain length in stages 4 and 6 might affect multiple processes such as protein transport,protein localization,and flavonoid metabolic process (Fig.2C).For example,LHCB3,encoding chlorophyll a-b binding protein 3 involved in PSII excitation energy transfer and charge separation during photosynthesis in thylakoids [32],showed increased,unchanged,and decreased patterns during stages 4,5,and 6,suggesting that LHCB3 response to LGS1 might provide energy for grain length increase only in stage 4.Cluster 39,containing the majority (445) of DE proteins showed a down-regulation pattern in stages 4 and 6 but a slight upregulation pattern in stage 5 (Fig.2B) and was enriched in biological processes such as cell-matrix adhesion,cell-substrate adhesion,cell junction assembly,homogalacturonan biosynthetic process and focal adhesion assembly (Fig.2C).These results suggested that the DE proteins in cluster 39 might be involved in grain length determination only in stage 5.Notably,in stage 5,ACS showed higher expression (FC ＞1.5) in NIL-LGS1 than in Samba.This gene encodes acetyl-coenzyme A synthetase,involved in fatty acid synthesis and leaf development in Arabidopsis [33,34],suggesting that ACS response to LGS1 might provide energy for grain length increase only in stage 5.Clusters 41,43,and 47,composed of 19,64,and 111 DE proteins,exhibited a significant decreasing tendency in stages 4 and 6 and only clusters 41 and 47 in NILLGS1 showed increasing trends in stage 5 (Fig.2B).The functional enrichment for these three clusters showed enrichment in glutamine processes such as glutamine family amino acid metabolic processes and glutamine family amino acid biosynthetic processes(cluster 41),pollen processes such as pollen exine formation and sporopollenin biosynthetic processes (cluster 43) and ethylene processes such as ethylene biosynthetic processes and ethylene metabolic processes (cluster 47),respectively (Fig.2C).Clusters 44 and 45,containing respectively 48 and 72 DE proteins,both showed higher expression in stage 5 than in stages 4 or 6 and showed decreased (cluster 44) and increased (cluster 45) expression patterns between Samba and NIL-LGS1 in stage 5 (Fig.2B).For example,GLY1 and ATR3,which generate lipid signals necessary for mediating defense responses by encoding glycerol-3-phosphate dehydrogenase [35] and catalyze the NADP-dependent reduction of cytochrome c by encoding NADPH-dependent diflavin oxidoreductase 1[36],exhibited opposite tendencies in Samba and NIL-LGS1 in stage 5.Taken together,these results indicate that the concordant regulation of multiple functional clusters composed of DE proteins is a key feature of the development of grain length in rice.

3.3.Phosphoproteome profiling and PPI networks

In the phosphoproteome profiling of Samba and NIL-LGS1 in stage 5,167 proteins,accounting for 8.3% of all identified proteins in stage 5(Fig.3A),were found to contain one or more phosphorylation events (Fig.3B).A striking sample is KIN14A,mediating the nuclear avoidance response under high-intensity blue light [37],contained one phosphorylation site in a serine residue.Subsequent differential-expression analysis indicated that respectively 17 and 19 phosphorylated proteins were significantly up-regulated and down-regulated in Samba relative to NIL-LGS1 (Fig.3B;Table S3).For example,UGD,required for the formation of cell wall ingrowths on the outer cell walls of nematode-induced syncytia[38],showed a greater than fourfold increase(FC=4.26,P-value=1.65e-147) in Samba relative to its expression in NIL-LGS1.Clusters 1–6 were composed of respectively 18,15,16,15,13,and 13 phosphorylated proteins (Fig.3C,D;Table S3).Six clusters were enriched in processes such as reproductive developmental process,chromatin modification,regulation of signaling pathway,protein catabolic process,intracellular transport,and developmental process(Fig.3E).For instance,in cluster 6,PWWP3A,encoding PWWP domain-containing DNA repair factor 3A involved in facilitating damage-induced chromatin changes and promoting chromatin relaxation [39],interacted with DLAT and ADA.Both PWWP3A and DLAT were found to be phosphorylated and showed downregulation in NIL-LGS1.In summary,the phosphoproteome profile and PPI networks in Samba and NIL-LGS1 of stage 5 suggest that phosphorylation events might influence rice grain development.

3.4.Genome-wide identification of alternative splicing events and TFs in stage 5

Fig.2.Identification of DE proteins and functional clusters.(A)Volcano plot shows DE proteins in Samba and NIL-LGS1 during the developmental stages 4,5,and 6.W4,5,and 6and M4,5,and 6 indicate developmental stages 4,5,and 6 of rice in Samba and NIL-LGS1,respectively.(B)Number and expression pattern of DE proteins in seven clusters.(C)GO enrichment analysis of DE proteins from seven clusters.The top 10 GO categories of each cluster are displayed.

Fig.3.Protein–protein interaction network of phosphorylated proteins.(A) Number of proteins with phosphorylated modifications in Samba and NIL-LGS1 during stage 5.Yellow,green and blue dashed rings indicate numbers of proteins in stages 4,5,and 6,proteins from stage 5,and phosphorylated proteins from stage 5,respectively.(B)Bar plot shows the number of phosphorylation events in each protein.The right table shows the number of DE phosphorylated proteins.(C) Number of proteins and phosphorylated proteins in six clusters.(D)Six protein–protein interaction clusters.Gray and brown nodes indicate proteins with or without phosphorylation.(E)Functional enrichment analysis of six clusters by Gene Ontology biological process.

To determine whether alternative-splicing (AS) events could regulate grain size in rice,we identified AS events in Samba and NIL-LGS1 of stage 5 from RNA-seq data in our previous study[10].A total of 16,064 typical AS events,including 3350 alternative-acceptor events (AltA),4594 alternative-donor events(AltD),5495 intron-retention events (IntronR),and 2625 exonskipping events(ExonS),were widely but nonrandomly distributed among rice chromosomes (Table S4;Fig.4A).Numbers of AltA,AltD,and IntronR,but not ExonS,were more enriched in Samba than in NIL-LGS1 (Fig.4B).Further investigation of the differentially expressed AS events with parameters P-value ＜0.05 and FC ＞1.5 showed that 2539 AltA,3615 AltD,1938 ExonS,and 4207 IntronR were differentially expressed between Samba and NIL-LGS1(Table S4;Fig.4C),suggesting that these AS events might modulate rice grain size.For example,ABCB1,regulating Arabidopsis photomorphogenesis and root development by mediating polar auxin transport [40],contained IntronR events with a decrease in NIL-LGS1.By GO-BP,as shown in Fig.4D,genes encoding proteins involved in cellular nitrogen compound metabolism,DNA metabolism,nucleic acid metabolism,and nitrogen compound metabolism were enriched in all four AS events.

Fig.4.Posttranslational events and GRN.(A)Schematic of four AS events and their distributionon rice chromosomes.(B)Number of AS events in Samba and NIL-LGS1 during stage 5.Blue and green bars indicate the numbers of AS events in Samba and NIL-LGS1,respectively.(C)Volcano plot shows AS events differentially expressed between Samba and NIL-LGS1.Red and green points indicate respectively up-and down-regulation of AS events in Samba and NIL-LGS1.(D) Heatmap plot shows GO enrichment analysis using biological processes for four AS events.(E)Transcription factors and their targeted genes.Venn plot(top)shows differentially expressed transcription factors.Bar plot indicates Pearson’s correlation coefficient between differentially expressed TFs and their targeted genes.(F)The key subnetwork from GRN.Green and yellow nodes indicate TFs and their targeted genes.(G) GO enrichment analysis using biological processes for genes constituting the subnetwork.

To further investigate the detailed transcriptional regulation of TFs to their targeted genes in Samba and NIL-LGS1 in stage 5,we first identified 39 TFs differentially expressed in Samba and NILLGS1,and then established a high-correlation network based on Pearson’s correlation coefficient (＞0.9 or ＜–0.9) between differentially expressed TFs and their potential target genes(Fig.4E,upper;Table S4).Respectively 5715 and 4633 genes showed positive and negative associations with these 39 TFs (Fig.4E,lower).We next sought to identify hub objects and subnetworks in a complex GRN with the MCC method using the CytoHubba [31] Cytoscape plugin.A network containing 1068 edges was constructed with 100 nodes,including 21 TFs and 79 target genes (Fig.4F).GO enrichment analysis showed that these genes were significantly enriched in biological processes,including regulation of macromolecule metabolism,floral organ morphogenesis,petal development,and ovule development (Fig.4G).

Collectively,these results showed that both alternative splicing events and dynamic GRNs were potential factors contributing to rice grain development.

3.5.Integrative multilevel omics analysis in Samba and NIL-LGS1 in stage 5

Integrative analysis of multidimensional omics data,we first investigated whether gene-generated alternative splicing events could translate into proteins with or without phosphorylation.As shown in Fig.5A and Table S5,respectively 144,196,77 and 176 translated genes also produced AltA,AltD,ExonS,and IntronR.These genes also encoded respectively 2,9,8,and 14 proteins with phosphorylation.An example is GAPCP1,encoding glyceraldehyde-3-phosphate dehydrogenase,which functions in pollen development[41].We next compared the abundance changes of four alternative splicing events and proteins from the same gene between Samba and NIL-LGS1.Distributions of correlations showed that AltA (positive 69 ＜negative 75),AltD (positive 89 ＜negative 107),and IntronR (positive 84 ＜negative 92) were negatively correlated with proteins from the same gene,a result inconsistent with those obtained for ExonS (positive 42 ＞negative 35).These instances suggested a potential relationship between protein levels and alternative splicing events,a speculation consistent with previous findings that diversity of transcripts might affect the process of translation [42].GO enrichment analysis for these genes with alternative splicing events revealed that several biological processes,including DNA metabolism,histone methylation,jasmonic acid metabolism,and response to abiotic stimuli,were enriched for genes that generated AltA,AltD,ExonS,and IntronR (Fig.5B).

We also investigated whether TFs and their targeted genes,which represent GRNs,could exist at the translational level.As shown in Fig.5C,only 5.6% (296/5274) of genes from GRNs were involved in the process of translation.Among these genes,130 showed consistent changes in their transcriptional levels and translational levels,whereas the alteration of 166 genes at the transcriptional level showed a negative association with their translational levels.For example,the transcriptional levels of LYM3 (LysM domain-containing GPI-anchored protein 3) and RH42(DEAD-box ATP-dependent RNA helicase 42)showed respectively positive and negative associations with their translational levels.GO enrichment analysis showed that several key processes including glucose catabolic process and fatty acid beta-oxidation,were significantly enriched in these genes encoding proteins.Interestingly,27 genes from GRNs were translated into proteins that were phosphorylated.Distributions of abundance changes showed that the number of negative correlations between transcriptional levels and translational levels was much greater than the number of positive correlations (23 vs.4) (Fig.5D).These results highlight the importance of specific TFs and their targeted genes and their interplay with grain development at the transcriptional,translational and posttranslational levels.

Overall,integrative multilevel omics analysis in Samba and NILLGS1 in stage 5 yielded several new insights into altering the grain size of rice.The negative correlation of alternative splicing events and proteins translated from their host genes and consistent changes in TFs and their targeted genes at the transcriptional and translational levels indicated that the regulation network at the posttranscriptional and transcriptional levels could be directly manifested at the translational level.Several genes from GRNs showed inconsistent changes at the transcriptional and translational levels,suggesting that a potential regulatory mechanism,such as regulation of translation level,might act in the biological process from extending transcripts to proteins to the development of grain.

3.6.A schematic of key potential mechanisms for regulating grain size mediated by LGS1

We have summarized the multi-omics dataset in a schematic of key mechanisms for grain length increase of rice in Fig.6 and Table S6.The expression of LGS1 showed the influence of multiple biological process including MAPK signaling,calcium signalling,cell proliferation,cell wall,energy metabolism,hormone pathway and ubiquitin-proteasome pathway.In agreement with previous findings [6–9] that LGS1,a major QTL for grain size,increased cell proliferation and cell expansion in the spikelet hull,our study also found that LGS1 could affect three cell proliferation-associated proteins:CDKA-1,POLA3,and NOP2B.Ten cell wall-associated proteins were found to respond to LGS1,providing further support that LGS1-mediated regulation of grain length was involved in cell proliferation,given that the cell wall strongly influences plant cell proliferation[43].Although there is no clear evidence that calcium signaling is involved in grain size control like the mitogenactivated protein kinase (MAPK) and G protein signaling pathway[44–48],we identified six DE proteins involved in calcium signaling:calcium-binding protein CRT3 and calcium-transported proteins HSP81-1,LHCB1.3,CXIP4,ACA4,STR4,and RGLG4,suggesting that calcium signaling might function in LGS1-mediated regulation of grain length.We also found that LGS1 could regulate MAPK signaling proteins such as MAPK6,encoding mitogen-activated protein kinase 6,which have been shown [47]to promote cell proliferation in spikelet hulls.The ubiquitinproteasome pathway has been shown [49–51] to affect seed size.Eight DE proteins were found to be involved in sequential action including activating,conjugating,ligation,deubiquitination,and 26S proteasome.UBC2 and CUL1,SBB,PEX2,and PUB44,affected by LGS1,were involved in conjugation and ligation of ubiquitin to target proteins,suggesting that these proteins might have conserved functions UBP8,encoding ubiquitin carboxyl-terminal hydrolase 8 involved in thiol-dependent deubiquitinase [52],might regulate grain size in a similar way that WTG1,the deubiquitinating enzyme,influenced cell expansion [50,51].ARS5 and PAE1 affected by LGS1,involved in the 26S proteasome by encoding arsenic tolerance 5 and alpha 5 subunit of 20S proteasome,might influence seed size,in agreement with a previous finding[53]that several factors such as RPT2a involved in the 26S proteasome influenced seed size.Recent findings [54–61] suggest that phytohormones such as brassinosteroid and auxin influence grain size.We also identified 3 (TTL1,GRF2 and HERK2) and 6 proteins (PIN4,TTL3,FQR1,ARF6,RPL24B,and CAND1)involved in brassinosteroid and auxin,suggesting that regulation of grain length mediated by LGS1 might involve crosstalk to several known knowledge on grain size control in rice.A striking sample is PIN4,a DE protein response to LGS1,acted as a component of the auxin efflux carrier[62].Other phytohormones such as ethylene,jasmonic acid,and abscisic acid also might be involved in regulation of grain length mediated by LGS1,as some phytohormone related-proteins were responsive to LGS1,suggesting a link between LGS1 and these phytohormones in regulation of grain length.For instance,OPR7 response to LGS1 has been reported[63]to involve the biosynthesis of jasmonic acid(JA).Finally,eight energy metabolism-related genes were identified to be responsive to LGS1 in our study,and might drive energy-requiring reactions such as grain length increase.

Our study demonstrates that the LGS1-mediated regulation of grain length might not only crosstalk some known knowledge on grain size control including MAPK signaling,cell proliferation,hormone pathway and ubiquitin-proteasome pathway,but also affect some unexplained biological processes such as calcium signaling,cell wall,and energy metabolism.

Fig.5.Integrative analysis in Samba and NIL-LGS1 using transcriptome,proteome,and phosphoproteome.(A) Number of genes generating four alternative splicing events and translating proteins.Red,green,blue,purple,light yellow,and deep yellow circles indicate AltA,AltD,ExonS,IntronR,proteins with or without phosphorylation.The plus(+)and minus(-)on the left or right of the slash(/)indicate up-regulation and down-regulation between Samba and NIL-LGS1 at the transcriptional or translational levels.+/+,-/-and +/-,-/+indicate a positive or negative correlation between the transcriptional or translational level.(B) GO enrichment analysis using biological processes for genes generating four alternative splicing events and translated proteins.(C) Genes from gene regulatory networks (GRNs) were translated into proteins.Venn (upper left)plot showing the number of genes in both constructed gene regulatory networks (GRNs) and translated into proteins.The bar plot (lower left) shows the number of genes positively or negatively correlated between the transcriptional and translational levels.The right plot shows the GO enrichment analysis for genes translated into proteins from gene regulatory networks.(D)Heatmap plot showing abundance changes of transcript and protein with phosphorylation differentiating Samba and NIL-LGS1.The right table indicates the number of genes showing positive or negative correlations between transcriptional and translational levels.

4.Discussion

Although LGS1 (Large grain size 1),encoding the OsGRF4 transcription factor,has been reported [6–9] to enlarge cells and increase numbers of cells and thus grain size and yield of rice,the underlying molecular mechanism and biochemical function of LGS1 remain to be elucidated.We performed transcriptomic,proteomic and phosphoproteomic analyses of young rice panicles to investigate internal dynamic functional networks involved in grain size mediated by LGS1 at different developmental stages.

Fig.6.Potential mechanisms for grain size mediated by LGS1. LGS1-mediated regulation of grain length might be coordinated by multiple biological processes including MAPK signaling,calcium signaling,cell proliferation,cell wall,energy metabolic,hormone pathway,and ubiquitin-proteasome pathway.Proteins labeled in red and black are respectively proteins affected by LGS1 and proteins known to be involved in grain size.

At the proteomic level,we identified seven highly correlated clusters grouped by the DE proteins (Fig.2A,B;Table S2).These clusters exhibited different expression patterns and participated in different biological processes (Fig.2B).Clusters 39,43,and 45 were enriched in biological processes of development of leaf,flower and seed,in agreement with previous studies[6,8,10,64,65]in which several GRFs affected leaf development,floral organogenesis,and seed size by increasing cell expansion.Previous studies [54–61] revealed that plant hormones such as brassinosteroid and auxin influenced grain size.Ethylene processes such as ethylene biosynthetic processes and ethylene metabolic processes were enriched in cluster 47,suggesting that ethylene processes similar to those involving brassinosteroid and auxin might be involved in grain length mediation by LGS1.Further investigation of how LGS1 regulates these processes will provide insight into modulating grain length.

At the phosphoproteomic level,the identified 17 up-regulated and 19 down-regulated phosphorylated proteins in Samba and NIL-LGS1(Fig.3B;Table S3),suggesting that phosphorylation modifications might be involved in LGS1-mediated regulation of grain length similarly to other genes controlling grain size such as MAPK6,encoding mitogen-activated protein kinase 6,which has been shown to promote cell proliferation in spikelet hulls [48].For instance,UGD2 including phosphorylated modifications required for cell wall formation[36],showed a more than fourfold increase in Samba compared with that in NIL-LGS1,suggesting that cell proliferation and cell expansion promoted by LGS1[7]might be mediated by phosphorylated proteins such as UGD2.Further integrated analysis suggested that dynamic phosphorylation events affected by LGS1 might form a dynamic PPI network including six clusters containing 18,15,16,15,13,and 13 phosphorylated proteins(Fig.3C,D;Table S3).Consistent with the notion that the QLQ domain from LGS1 is involved in protein interaction with chromatin-remodeling factors [45],we observed that PWWP3A and DLAT in cluster 6 were both phosphorylated and showed down-regulation in NIL-LGS1.Both PWWP3A and DLAT have been reported [37] to facilitate chromatin changes and promote relaxation.How phosphorylated proteins regulated by LGS1 modulate grain length awaits further study.

At the transcriptomic level,the enriched AS events in Samba(Fig.4B),suggesting that AS events might be involved in LGS1-mediated regulation of grain length,a notion supported by the finding that differentially expressed analysis of AS events between Samba and NIL-LGS1(Table S4;Fig.4C).Moreover,hub objects and sub-networks were identified in a complex GRN containing 39 TFs with differential expression and their potential targeted genes(Fig.4E,upper;Table S4) using the MCC method in CytoHubba[31].Surprisingly,the key sub-networks were also enriched in floral processes(Fig.4G),consistent with the result from cluster 43 at the proteomic level,suggesting that integrative multilevel omics analysis need to elaborate their interplay with the development of grain at the transcriptional and translational levels.

Our three-level omics analysis revealed that gene generated alternative splicing events also translate into proteins with or without phosphorylation (Fig.5A).Distributions of correlations between four alternative splicing events and proteins from the same gene showed that ＞50% of all AltA,AltD,and IntronR,but not ExonS,were negatively correlated with proteins from the same gene,consistent with previous findings [42] that diversity of transcripts might affect the process of translation.Only 5.6%(296/5274) of genes from GRNs were involved in the process of translation and 130 and 166 of these genes exhibited consistent and non-consistent changes in their transcriptional levels and translational levels.Overall,integrative multilevel omics analysis indicated that the regulation network at the posttranscriptional and transcriptional levels could not only be directly manifested but also represent a potential regulatory mechanism from extending transcripts to proteins to the development of grain.Further investigation such as ribosome profiling sequencing [66] might reveal the regulatory mechanism linking transcriptional and translational levels involved in LGS1-mediated regulation of grain length.

A recent finding [67] revealed that the final size of rice grain is coordinately controlled by cell proliferation and cell expansion in the spikelet hull.Our study also provided support of the notion that LGS1-mediated regulation of grain length was involved in cell proliferation in the form of 3 cell proliferation-and 10 cell wallassociated proteins responsive to LGS1.Calcium signaling and MAPK pathway were reported [33,47] to modulate cell proliferation.Thus,calcium signaling and MAPK pathway might be identified as key regulators mediating the regulation of grain length by affecting cell proliferation,as LGS1 could regulate several proteins involved in calcium signaling and MAPK pathway.Interestingly,MAPK signaling proteins such as MAPK6 have been shown [47]to increase proliferation in spikelet hulls,suggesting that LGS1 might be involved in a MAPK cascade to control grain length in rice.Our study found that several factors affected by LGS1 including UBC2,CUL1,SBB,PEX2,PUB44,UBP8,ARS5,and PAE1 involved in sequential action including activating,conjugating,ligation,deubiquitination,and 26S proteasome might process conserved function in seed size.Thus,LGS1 might control grain length by modulating the ubiquitin-proteasome pathway in rice.In both Arabidopsis and rice,BR promoted seed growth [54,58,59].For example,GSK2,resulting in typical BR-deficient phenotypes,has been shown to directly interact with the transcriptional activator LGS1 and inhibit its activity [7].Our study found that LGS1 could regulate three brassinosteroid response proteins:TTL1,GRF2,and HERK2,suggesting that LGS1 might be involved in the brassinosteroid pathway in regulation of grain length.Auxin signaling and transport have been reported[55,56,60]to be involved in seed size control.We found that PIN4,affected by LGS1,encoding auxin efflux carrier component 4,acting as a component of the auxin efflux carrier,might be similar to BG1,encoding a membranelocalized protein,increasing auxin basipetal transport and altering cell proliferation and cell expansion of grain size [60],suggesting that LGS1 might modulate auxin transport in regulating grain length.Other proteins involved in the ethylene,jasmonic acid,and abscisic acid pathways,which have been shown [68–70] to be responsive to external plant stresses,were regulated by LGS1,suggesting that LGS1-mediated regulation of grain length might also crosstalk with these phytohormones in regulating stress responses,consistent with our previous finding that LGS1 mutation also enhanced cold tolerance at the seedling stage in rice [48].Taken together,LGS1-mediated regulation of grain length might represent coordination between multiple biological process including MAPK signaling,calcium signaling,cell proliferation,cell wall,energy metabolism,hormone pathway,and ubiquitinproteasome pathway.

5.Conclusions

In this study,we performed transcriptomic,proteomic,and phosphoproteomic analyses of young rice panicles in Samba and NIL-LGS1 at three developmental stages to investigate the internal dynamic functional networks involved in grain size that are mediated by LGS1.Proteomic profiles in stages 4,5 and 6 suggested that concordant regulation of multiple functional clusters are key features of the development of grain length in rice.In stage 5,these differentially abundant phosphorylated proteins participating in protein–protein interaction networks,dynamic alternative splicing events,and gene regulatory networks based on transcription factors might play accessory roles in determining rice grain size.Integrative multilevel omics analysis of young panicles of Samba and NIL-LGS1 in stage 5 revealed that the dynamic regulatory network not only was consistent at the transcriptional,posttranscriptional and translational level but also implicated potential regulation of translation levels in the biological process that extended from transcript to protein to the development of grain.Detailed function analysis revealed that LGS1-mediated regulation of grain length might be coordinated by multiple known biological process such as MAPK signaling,cell proliferation hormone pathway,and ubiquitin-proteasome pathway or several novel biological process including calcium signaling and energy metabolism.Thus,the LGS1-mediated regulation of grain size was affected by dynamic transcriptional,posttranscriptional,translational,and posttranslational changes.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Fangyu Chen:Conceptualization,Methodology,Formal analysis,Writing-original draft,Writing–review &editing;Yongsheng Wang:Visualization,Validation,Data curation,Formal analysis,Writing–original draft,Writing–review&editing;Zesen Zhang:Formal analysis,Visualization;Xiaolong Chen:Formal analysis,Validation;Jinpeng Huang:Formal analysis,Software;Zhiming Chen:Formal analysis,Writing– review &editing;Jingsheng Zheng:Methodology,Investigation;Liangrong Jiang:Visualization,Software;Yumin Huang:Conceptualization,Resources;Houcong Wang:Resources,Supervision;Rongyu Huang:Project administration,Resources,Funding acquisition,Writing– original draft,Writing– review &editing.

Acknowledgments

This work was supported by the National Key Research and Development Program of China (2017YFD0100103),the Seed Industry Innovation and Industrialization Project of Fujian Province (fjzycxny2017004,zycxny2021004),the Program on Technology of Fujian Province (2020NZ08016,2020N0049),and the Open Program of State Key Laboratory of Rice Biology of China(170101).The authors are grateful to Dr.Changchuan Xie and Dr.Yaying Wu (Analysis and Testing Center,School of Life Sciences,Xiamen University) for their help in mass spectrometry analysis.

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2021.11.007.