Integrated transcriptomic and proteomic analysis reveal new insights into stone cell formation in Korla pear

Background Korla fragrant pear( P•sinkiangensis Yü)is a famous local variety of Xinjiang China. One difficulty is the high stone cell content of these pears, which causes the formation of rough skins on the fruit. To elucidate the underlying mechanisms of stone cell formation, parallel analyses of the transcriptome and proteome was performed to identify important regulators and pathways involved in stone cell formation. Results Fruit samples were collected at three important time points depending on the stages of stone cell formation (20, 50 and 80 days after flowering). A total of 24268 differentially expressed genes (DEGs) and 1011 differentially accumulated proteins (DAPs) were identified from all the time points. Function analysis of the differential genes/proteins revealed that a set of candidates was associated with stone cell formation. These candidates mainly enriched in pathways involved in lignin biosynthesis, cellulose and xylan biosynthesis, S-adenosylmethionine (SAM) metabolic process, Reactive oxygen species (ROS) production, and cell death. We mined a total of 253 DEGs, and 100 DAPs, 63 of which were significantly changed at both the transcript and protein levels during fruit development. Conclusions Our findings reveal that some intriguing genes/proteins were previously unrecognized related with the sclereid formation, which provided new insights into molecular processes regulating sclereid accumulation in pear pulp. of parenchyma cells [ 19 – 21 ] . Our results showed that stone cell formation mainly occurred during the early stage of fruit development. Similar results have been reported in other pear varieties [ 1 , 2 , 11 ] . In plant secondary walls, i.e. the proportion of cellulose, hemicelluloses and lignin, may vary among different plant species and even in different cell types of the same species. Some specialized secondary walls may be composed predominantly of one or two polymers. For example, cotton fibers contain > 90% cellulose, and phloem fibers, on the other hand, contains no lignin [ 6 ] . We found that pear stone cells mainly composed of cellulose (52%), hemicellulose (23%), lignin (20%) as well as a little amount of polysaccharide (3%). Our results are supported by several early reports in pear [ 4 , 22 ] , indicating that the stone cell formation is more likely to result from the differentiation of distinct cell types.

inference, in this study, parallel analyses of the transcriptome and proteome of Korla pear fruits from the three developmental stages were carried out using RNA-seq and TMT technologies, aiming to identify differential genes/proteins that may be involved in the formation of stone cells. Our integrated omics data (transcriptprotein) revealed several new and interesting insights regarding pear fruit stone cell formation. Knowledge gained through this study could be helpful for understanding the molecular mechanisms of stone cell formation in pear and could be further used to develop strategies to reduce or eliminate the probability of stone cell formation.

Physiological and morphological changes
Stone cell staining and cell apoptosis detection results as well as cell wall polymers (cellulose, hemicellulose, lignin and polysaccharide) composition in the stone cell is shown (Fig. 1) to illustrate the difference between experimental time stages. Greater differences for stone cell, cell apoptosis and cell wall polymers composition were observed in Day20 compared to fruits sampled in later stage. It is evident from the figure that stone cell formation mainly occurred during the early stage of fruit development, and the apoptotic period basically overlapped with the period of stone cell formation (Fig A,B). On the other hand, the main components for pear stone cells were cellulose (52%) representing about half of stone cell biomass, hemicellulose (23%), lignin (20%) as well as a little amount of polysaccharide (3%). Therefore, formation of stone cells cannot be fully explained by the deposition of lignin in the secondary cell walls. In the later section, we will pay close attention on the enrichment of genes and proteins related to secondary cell wall polymer synthesis.

Overview of transcriptomic and proteomic data
In order to identify differential genes/proteins that may be involved in stone cell formation in Korla pear, combined transcriptome and proteome analysis were performed by the RNA-seq and the TMT techniques. A total of 42893 transcripts and 7904 proteins were obtained, respectively, in which 24268 genes and 1011 proteins were showed differential expression during the developmental time course. Correlation analysis showed that approximately 13% (1010/7904) of the proteins correlated with DEGs, including 700 DAPs. Of the 1011 DAPs, 311 DAPs (30.8%) were not correlated with DEGs ( Fig. 2A). The analysis of DEGs and DAPs in our study focused on pathways that involved in lignin biosynthesis (94 genes and 31 proteins), cellulose and xylan biosynthesis (46 genes and 18 proteins), SAM metabolic process(10 genes and 3 proteins), apoplastic ROS production (16 genes and 2 proteins), and cell death (14 genes and 6 proteins) (Fig. 2B). Information on all DEGs and DAPs on these pathways can be found in Data S1-S3. However, information provided in the text included genes and proteins with relatively high expression that were differentially expressed during the developmental time course . To validate the RNA-seq results, 15 genes were randomly selected from the DEGs for qPCR analysis in the transcript levels of three time courses. As shown in Figure S1, a high agreement was observed between the RNA-seq and qPCR data, demonstrating the validity of RNA-seq data of genes with different transcript abundance.

DEGs and DAPs involved in lignin biosynthesis
Omics data revealed that 94 DEGs and 31 correlated DAPs encoding several monolignol biosynthetic enzymes were identified in Korla pear pulp (Data S1). These genes and proteins encodes upstream enzymes of phenylpropanoid pathway (PAL, C4H, C3'H,4CL, HCT, CSE and CCoAOMT ), enzymes in the monolignol-specific pathway ( CCR, COMT, F5H and CAD) and enzymes involved in lignin polymerization (PRX, LAC). Majority of these genes showed peak expression in the early stage of fruit development, and gradually down-regulated from Day20 to Day80, showed similar trends with lignin content in pear during fruit development. Some transcripts, however, showed no differential expression from Day 50 to Day 80, indicating that lignin accumulation mainly occured in the early stage of fruit development. Interestingly, expression patterns of some transcripts encoding PAL (one transcript), C4H (two transcripts), 4CL (two transcripts), HCT (one transcript) and CCR (five transcripts), CAD (two transcripts), CSE (two transcripts) and PRX (five transcripts) showed the opposite trend, indicating that they may be involved in the biosynthesis of other secondary metabolites in the phenylalanine pathway.
However, the expression levels of these genes were almost negligible, compared with expression levels of the genes positively regulated with lignin content.

DEGs and DAPs associated with cellulose and xylan biosynthesis
Cellulose, consisting of linear chains of β-1, 4-linked glucosyl residues, is the predominant constituent of secondary cell walls. 29 DEGs in transcriptome profile were assigned to cellulose biosynthesis during the experimental time course, while 7 DAPs in proteome data (Fig. 4). UDP-glucose is the direct precursor for cellulose. Therefore pathways that lead to UDP-glucose production can potentially have an impact on cellulose biosynthesis UDP-glucose pyrophosphorylase (UGPs) and sucrose synthase (SuSy) are the main enzymes that have been shown to produce UDP-glucose in plants [16] . But we didn't find any differential expressed genes encoding UGPs. We found 4 DEGs and 1 DEP that encodes sucrose synthase (SuSy) and their expression pattern were the same trend as in cellulose content, indicating that SuSy tightly associated with the processes of cellulose biosynthesis. Sucrose can be directly converted to UDP-Glcose and fructose by SuSy, Hexose phosphates pool can also have an effect on UDP-glucose pool size [17] . Fructose is a feedback inhibitor of sucrose synthase, and fructokinase(FRK)phosphorylate fructose to produce fructose-6-phosphate. In our omics data, 4DEGs and 2 DEP that encode FRK were selected. All of these genes and proteins showed similar expression as SuSy during the experimental period, indicating that FRK might indirectly affect cellulose production through modulation of sucrose synthase activity and UDP-glucose production by reducing intracellular fructose pools. Cellulose synthase (CesA) use UDP-glucose as substrate for the biosynthesis of cellulose chains. In arabidopsis CESA1, CESA3 and CESA6 (or CESA6-like) are required for primary wall cellulose synthesis whereas CesA4, CesA7, CesA8 constitute the secondary wall CesA complexes. In our study, we obtained a total of 21 functional genes and 4 proteins encoding different CesAs enzymes and cellulose synthase like proteins (CSLs), among which homologous genes of arabidopsis CesA4 (2 transcripts), CesA7 (2 transcripts), CesA8 (2 transcripts) showed the same expression patterns as the genes related to lignin synthesis, in accordance with trends of lignin and stone cell contents (Fig. 1), suggesting that these genes and proteins might be involved in the development of stone cells in Korla pear pulp. However, homologous genes of arabidopsis CesA1 (2 transcripts), CesA2 (2 transcripts), CesA3 (2 transcripts) and most of the CSLs (5 transcripts) showed the opposite trend, indicating that these genes mainly participate in the growth of primary walls of parenchyma cells in Korla pear pulp. In addition we screen for 2 DEGs and 1DEP encodes COBRA-like4 proteins (COBL4). COBRA, encodes a glycosylphosphatidylinositol-anchored protein, is involved in the orientation of deposition of cellulose microfibrils. In arabidopsis COBL4 is required for cellulose synthesis in the secondary cell wall. Expression trends of these genes were similar to that of CesA genes associated with secondary wall synthesis, implying that this protein likely be involved in the assembly of cellulose synthase complexes.
Xylans are mainly composed of a linear backbone of β-(1, 4)-linked D-xylosyl residues. UDP-xylose is a nucleotide sugar required for xylan backbone elongation. Similarly, 25 DEGs and 13 DAPs were identified as xylan biosynthetic genes and proteins that showed the same expression trends as the genes related to cellulose and lignin synthesis (Fig. 4). These genes and proteins encoding a suite of enzymes that is required for the UDPxylose formation, xylan backbone elongation, side chain addition and modification as well as the reducing end sequence. Among these, 1 gene with no correlated protein encodes UDP-Glucose dehydrogenase(UGD), which is an enzyme converting UDP-glucose to UDP-glucuronic acid in an irreversible manner; 6 genes and 4 proteins encodes UDP-Xylose synthase (UXS), which irreversibly convert UDP glucuronic acid to UDP-xylose; 4 genes and 4 correlated proteins encodes xylosyltransferases (XylT) required for elongation of xylan backbone; 4 genes with no correlated proteins encodes glucuronosyltransferase (GUXs ) required for the addition of both glucuronic acid and 4-O-methylglucuronic acid branches to xylan backbone; 6 genes and 3 correlated proteins encode glucuronoxylan 4-O-methyltransferase (GXMT ) which is responsible for the glucuronic acid (GlcA) side chain methylation in xylan; 3 genes (IRX7) with one correlated proteins required for formation of the xylan reducing end sequence (XRES) which has been speculated to act as either a primer or terminator for glucuronoxylan biosynthesis [18] (Fig. 4 and Data S2). Based on the data in Fig. 4, expression levels of genes related to UDP-Xylose formation (UXSs), elongation of xylan backbone (IRXs) and side chain methylation (GXMs) were significantly higher than other xylan related genes, implying their important roles in xylan or glucuronoxylan synthesis. Most of secondary wall genes and xylan related genes were preferentially accumulated on Day20, while primary wall genes were relatively abundant on Day80 (Data S2). These results suggested that greater secondary thickening in parenchymal cells was occurred in the early stage of fruit development.

Other genes and proteins involved in secondary cell wall biosynthesis
Genes involved in apoplastic ROS production (16 genes and 2 proteins), SAM metabolic process (10 genes and 3 proteins), and cell death (12 genes and 6 proteins) were well represented in the candidate gene lists ( Fig. 5 and Data S3). Apoplastic ROS can be used as an oxidant in cell wall cross-linking, and are required for lignin polymerization. ROS being a major stress factor also play an important role in apoptosis induction. Apoplastic ROS related genes and proteins were preferentially accumulated in Day20, with significantly lower expression in later stages. For example, 3 genes and 1correlated protein annotated as respiratory burst oxidase homolog (RBOH) were expressed at high levels in the early stage then were down-regulated to varying extents. More accumulation of transcripts of genes encoding superoxide dismutase (SOD), polyamine oxidase (PAO), PRXs and germin-like proteins (GLPs), was also observed in the early stage of fruit development.
Genes encoding enzymes that related to S-adenosylmethionine (SAM) metabolic process were differentially expressed during fruit development in a similar manner as ROS related genes ( Fig. 5 and Data S3). SAM acts as a methyl group donor for various biologically important processes including the synthesis of cell wall polymers. Moreover, SAM acts as amino carboxylpropyl group donor in polyamine generation, which is an important substrate for apoplastic H 2 O 2 production. Numerous genes related to programmed cell death were also expressed in pear pulps, and all showed their highest expression in Day20 consistent with the trends of ROS generating genes during fruit development ( Fig. 5 and Data S3).

Discussion
Stone cells are a type of sclerenchyma cell formed by the secondary thickening of cell walls followed by the deposition of lignin on the primary walls of parenchyma cells [19][20][21] . Our results showed that stone cell formation mainly occurred during the early stage of fruit development. Similar results have been reported in other pear varieties [1,2,11] . In plant secondary walls, i.e. the proportion of cellulose, hemicelluloses and lignin, may vary among different plant species and even in different cell types of the same species. Some specialized secondary walls may be composed predominantly of one or two polymers. For example, cotton fibers contain > 90% cellulose, and phloem fibers, on the other hand, contains no lignin [6] . We found that pear stone cells mainly composed of cellulose (52%), hemicellulose (23%), lignin (20%) as well as a little amount of polysaccharide (3%). Our results are supported by several early reports in pear [4,22] , indicating that the stone cell formation is more likely to result from the differentiation of distinct cell types.
Lignification is an integral part of the differentiation process of several specialized cell types to fulfil their physiological function [13] . Lignification process is vary among different cell type, for instance, differentiation process in tracheary elements (TEs) require cell death to achieve their lignification, in which enzymes and/or substrates required for cell wall lignification was provided by neighbouring cells [23][24] . In contrast to TEs, lignification in fibre cells is carried out by cell autonomous processes [23][24] . Lignin composition varies among these cell types, in TEs, mainly composed of G-units, and the sclerenchyma fibres have lignified secondary cell walls mainly composed of S-units [25][26] . In our previous research, we found that sinapyl alcohol, which gives rise to S-unit, was the dominant lignin monomer in Korla pear pulp, accounting for 77.2% of total alcohol monomer [21] . These observations have led to the idea that stone cell differentiation in pear pulp might proceed through in a cell-autonomous process. Since lignification in sclerenchyma fibre cells begins with a microtubuledependent deposition of xylan and cellulose in the secondary cell wall. Therefore being a type of sclerenchyma, differentiation of stone cells must begin with a cellulose and xylan secondary cell wall deposition. Supporting this theory, we found that many secondary cell wall genes and proteins (53 genes and 20 proteins) showed peak expression in the early stage of fruit development, consistent with the changing trends of stone cell and lignin.
Cellulose is the load-bearing unit, in which cellulose microfibrils cross-link with hemicelluloses to form the framework of secondary walls. We found that cellulose content accounted for more than a half of stone cell biomass, indicating that it plays a decisive role in the differentiation of stone cells. Cellulose synthase (CesA) use UDP-glucose as substrate for the biosynthesis of cellulose chain [27][28] . Therefore, pathways that lead to UDPglucose production can potentially have an impact on cellulose biosynthesis. SUSs and UGPs are the main enzymes that produce UDP-glucose in plants [16] . In our study, we found 4 DEGs and 1 DEP that encodes SuSy. It is reported that in hybrid aspen inhibition of the SuSy activity reduced carbon allocation to all the secondary cell wall polymers, however, its over expression only increased cellulose deposition in hybrid poplar [29][30] . In accordance with these findings, in our present work several genes encoding SuSy in pear pulp showed peak expression in Day20, then down-regulated more than 20-fold in Day80 in parallel with cellulose content, indicating that SuSy tightly associated with the processes of cellulose biosynthesis. In our early research, we found that SuSy expressed higher enzymatic activity in the early stage of fruit development in Korla pear, SuSy activity decreased first and then increased gradually 60 days after flowering. Accordingly, sucrose accumulations rose slowly in the later stage,began to increase 60 days after flowering [31] , suggesting that in the early stage of fruit development sucrose are mainly used to cellulose synthesis. Fructose is a feedback inhibitor of SuSy [17] . FRK phosphorylate fructose to produce fructose-6-phosphate.We identified 4DEGs and 2 DEP that encode FRK showed similar expression as SuSy during the experimental period, indicating that FRK might indirectly affect cellulose production through modulation of SuSy and UDP-glulcose production by reducing intracellular fructose pools [16][17]32] .
In our study, CESA genes that encode cellulose synthases were differentially expressed during the experimental time course. As expected, genes encoding CESAs that are involved in secondary cell wall cellulose synthesis [33][34] (CESA4, CESA7 and CESA8) showed peak expression in the early stage of fruit development and gradually decline in expression in the later stage. In contrast, genes required for primary wall cellulose synthesis [33][34] (CESA1, CESA2 and CESA3) expressed the opposite trend. In addition correlated proteins of AtCESA7 and AtCESA8 were also found to be differentially expressed during the experimental time course, and showed the same trend as their corresponding genes. It has been demonstrated that presence of AtCESA4, AtCESA7 and AtCESA8s is required for their assembly into complexes and their correct targeting to the plasma membrane [35][36][37][38][39][40][41][42][43][44][45][46] . Mutation of any of these three genes causes a severe reduction in cellulose content and secondary wall thickening [33][34] . The same results were reported in rice, brachypodium and populus [37][38][39] . These findings suggest that CESA genes required for secondary cell wall cellulose synthesis participate in the process of stone cell differentiation.
Hemicellulose is another major component of pear stone cells, accounting for 23% of stone cell biomass. Xylan, a major hemicellulosic component, is essential for secondary wall strength [40] . UDP-xylose is a nucleotide sugar required for xylan backbone elongation. UDP-xylose was formed by serial actions of UGD and UXS. Genes encoding these two enzymes showed higher expression in the early stage of development. Expression levels of UXSs were more than 100-fold higher than UGDs, indicating that UXS is the major enzyme responsible for UDPxylose. This is supported by the evidence that antisense downregulation of UXS leads to high glucose-to-xylose ratios in xylem walls due to fewer xylose-containing polymers. Such plants also have altered vascular organization and reduced xylans in their secondary walls [41][42] .
A number of enzymes have been identified as potentially involved in the elongation of the xylan backbone and its side chains addition. Coincident with expression of CESA genes involved in secondary cell wall biosynthesis, genes or proteins encoding secondary cell wall xylan synthesis enzymes [43][44][45][46][47] (IRX10, IRX15, IRX15-L) showed peak expression in the early stage of fruit development and down-regulated in the later stage. In Arabidopsis, all of irx mutants exhibit similar phenotypes, including severe reductions in xylan content and chain length and secondary wall thickness [43][44][45][46][47] . Interestingly, there was no significant reduction of xylan:XylT activity in irx15 and irx15-L mutants, and the irx15 irx15-L double mutant has replacement of GlcA with meGlcA xylan side chain. Therefore, IRX15 and IRX15-L define a new class of genes involved in xylan biosynthesis, whose function remains elusive [48][49][50] .
Another group of genes involved in addition of GlcA side chains and side chain methylation such as GUXs and GXMTs showed peak expression in the early stage of fruit development, and expressed at low levels in the later stages. Our results were supported by the evidences that these genes are preferentially expressed in secondary wall-forming cells, and mutations in these genes lead to a complete loss of GlcA or methylated GlcA side chains in xylan [51][52][53] . Furthermore, gux1 gux2 gux3 triple mutants also has found with an irregular xylem phenotype with reduced secondary wall thickness [52] , indicating that GlcA substitution of xylan is important in maintaining wall integrity. Genes such as FRA8/IRX7, F8H, IRX8 and PARVUS have been identified to be required for the biosynthesis of XRES, all of which are expressed in secondary wall-forming cells [54][55][56] . Consistent with this notion, we screened three DEGs and one DAP that annotated as probable IRX7, which were substantially upregulated in the early stage, and down-regulated in the later stage, indicating that these genes are required for xylan synthesis. What is noticeable is that mutation in XRES related genes can't affect either XylT or GlcAT enzyme activity, despite the loss of XRES [48,50] .
Lignin comprises approximately 20% of stone cell biomass (Fig. 1). Previous studies have shown that lignin related genes are closely associated with the development of lignified xylem and interfascicular fibers in stems [57][58][59] . Similarly, in present study more than 90 genes and 31 proteins involved in lignin biosynthesis were differentially expressed during pear fruit development (Data S1). Majority of these genes showed peak expression in the early stage of fruit development, positively correlated with lignin content, but some transcripts showed the opposite trend. However, the expression levels of genes positively regulated with lignin content were universally much higher than that of negatively regulated genes, indicating that they are the core genes responsible for lignin synthesis during stone cell formation in pear fruits. Genes and enzymes involved in monolignol biosynthesis have been studied extensively in different pear varieties [1,2,8,21,[60][61] . Therefore roles of these genes during lignin accumulation in pear pulp will not be discussed in detail. Instead, we focus mainly on the genes involved in lignin polymerization.
After their synthesis monolignols are transported into cell walls, where they are oxidized by laccases and peroxidases to form monolignol radicals, which are then polymerized into lignin by radical coupling [62][63][64][65] . In our study all of LAC genes and majority of PRX genes showed similar expression profile as monolignol biosynthetic genes. Although some PRX genes showed the opposite trend, their expression levels almost negligible during the entire period, compare to up regulated PRXs in the early stage. LACs are reported to be coregulated with monolignol synthetic genes and secondary cell wall-forming genes [50,66] . We found that 6 of 17 LAC genes, three transcripts of LAC4-L, two transcripts of LAC17 -L and one transcripts of LAC2-L, expressed 10-150-fold higher than other transcripts, indicating that these are the core LACs for the formation of monolignol radicals. Mutant studies have shown that simultaneous mutation of LAC4 and LAC17 results in a reduction in lignin content and lac4lac11lac17 triple mutant results in a phenotype of severe growth defect and a loss of lignification in root vessels [63,[67][68] . Like laccases, are expressed in lignifying vascular tissue.In model plant down-regulation or silencing of genes (PRX2, PRX3, PRX25, PRX60, PRX71, PRX72, etc) encoding peroxidases resulted in reduced lignin accumulation and altered lignin composition [65,69] . Consistent with these findings, our study showed that six transcripts of PRXs (two PRX12-L, two PRX42-L, one PRX25 and one PRX72-L) were expressed remarkably higher (up to several hundred times) than other transcripts throughout the experiment. Correlated proteins of PRX12-L and PRX72-L also up regulated 10-100-fold higher than other correlated proteins, suggesting their importance for lignin polymerization in Korla pear.
Beside the lignin monomers, peroxidases and laccases require H 2 O 2 and O 2 , respectively to form monolignol radicals. Our study showed that ROS related genes and their encoding proteins were more abundantly expressed in the early stage, and down-regulated in the later stage. Certain ROS, like H2O2, can be used as an oxidant in cell wall cross-linking, or as a signalling molecule controlling various biological processes [70][71] . It is reported that casparian strip and the sclerenchyma cells are found with an increase of H 2 O 2 . Scavenging of H 2 O 2 reduces the lignification of the casparian strip, TEs and the extracellular lignin accumulation in cell cultures [24,[72][73] . Together these findings indicate that apoplastic ROS play important role in the lignification process (stone cell formation). H 2 O 2 is also a key modulator of programmed cell death (PCD) which is related with a number of developmental processes including tissue (tracheary elements, fibers, ect) differentiation [74] . The present study found several highly transcribed cell death related genes and proteins (AED3-like, MC1-like, ACD11-like), expression pattern of which were similar to that of ROS related genes.

Conclusions
On the basis of these multi-level results, we conclude that stone cells are mainly composed of cellulose, hemicellulose and lignin as well as a little amount of polysaccharide. In addition to lignin synthetic genes, many genes related to secondary wall development, cell apoptosis and ROS production were highly expressed during the critical period of cell differentiation. Therefore, formation of stone cells cannot be fully explained by the deposition of lignin in the secondary cell walls. Together with lignin related genes, other genes/proteins encoding enzymes responsible for the synthesis of other secondary cell wall components (cellulose, xylan) co-regulate the differentiation process of stone cells. These previously unrecognized genes/proteins provided new insights into molecular processes regulating stone cell accumulation in pear pulp.

Plant materials and sampling periods
Fruits of 'Korla pear' were obtained from 20-year-old pear trees grown on a farm (Korla, Xinjiang, China). The first author is responsible for the entire field sampling process. In the last 5 years, all of our field experiments on Korla pear have been carried out in this experimental park to ensure the consistency of plant materials and the reliability of pear variety. There is no voucher specimen of this material has been deposited in a publicly available herbarium. Three important time points were sampled depending on the stages of stone cell formation in 2018, Day20 (20 days after flowering) represents the early development stage of stone cells; Day50 (50 days after flowering) represents the later stage of stone cell formation, and Day80 (80 days after flowering) represents the stationary phase of stone cells (with lower stone cell contents ). In the early stage (Day20 and Day50), each biological replicate consisted of fifteen single fruit samples collected from five independent trees, and ten fruits were collected for each biological replicate in the later stage. The fruit skins were peeled with a peeler, and cylindrical samples, 15 mm in diameter, were sampled from the pulp, frozen in liquid nitrogen, ground into mixed samples and stored at −80°C until further analysis. The same batch of sampling materials was used for transcriptome, proteome, physiological and biological analyses.

Stone cell separation and structural carbohydrate detection
Stone cell separation was carried out as described by Brahem et al (2017) [4] . Briefly 2 kg of frozen pear flesh was used. The separation of stone cells from parenchymatous tissues was based on their higher density. After washing with acetone: water, the pulp was stored overnight at 4℃in the same acetone:water (60:40, v:v) solution. Stone cells sedimented and formed a distinctive yellow layer at the bottom of the container. The stone cells were collected, oven-dried, ground into a uniform powder. The lignin content was measured using the Klason method [75] . A small amount (0.2 g) of this powder was extracted with 15 ml of 72% H 2 SO 4 at 30 • C for 1 h, combined with 115 ml of distilled water and boiled for 1 h. The volume was kept constant during boiling. The liquid mixture was filtered and the residue was rinsed with 500 ml of hot water, air-dried and weighed. The hemicellulose content was determined by 2% hydrochloric acid hydrolysis method combined with DNS method [76] . The cellulose content was determined by Anthrone-Sulfuric acid colorimetry method [76] . The polysaccharide content was tested by phenol-sulfuric acid method [77] .

Stone cell staining and apoptosis detection
Stone cells were sectioned using a procedure reported by Cai et al. (2010) [2] . Pulp, from 2.0mm under the peel was collected, cut into appropriate sizes, and fixed with FAA fixative solution (50 ml 50% ethanol, 5ml glacial acetic acid, 5ml formalin). Pulp sections were cut manually, stained with 0.1% Safranin, and studied under an optical microscope. Apoptosis was detected by TUNEL assay kits from Shanghai beyotime (www.beyotime.com). The kits were carried out according to the manufacturer's instructions.

Transcriptomic sequencing
Transcriptomic sequencing was performed by Novogene Bioinformatics Technology Co., Ltd (Beijing, China). In brief, total RNA was isolated by use of TRIzol reagent and the purity, concentration and integrity of RNA were determined. Once RNA samples were qualified, Oligo-dT beads were used for the isolation of mRNA. After that, the isolated mRNA was fragmented into short sections and used for the synthesis of first-strand cDNA. The synthesis of second-strand cDNA was conducted by use of dNTPs, DNA polymerase I and RNase H. The purification of double-stranded cDNA was performed using AMPure XP beads, and then the purified doublestranded cDNA were used for end reparation, "A" base addition and ligation of sequencing adapters. After size selection with AMPure XP beads, PCR amplification was performed and the production was purified using AMPure XP beads for library construction. Finally, the library was sequenced as 2 × 150 bp paired-end reads on Illumina HiSeq TM P150 sequencer (Illumina, San Diego, CA, USA). Raw data were cleaned by removing reads with adaptors, low quality (> 50 %) or high-proportion unknown bases (> 10 %) and then were assembled by use of Trinity methodology. The annotation of transcriptome sequences were performed using the following seven databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI nucleotide sequences), Pfam (Protein family), KOG/COG (eukaryotic Ortholog Groups and Clusters of Orthologous Groups of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology). The read numbers were transformed to FPKM (fragments per kilobase of transcript sequence per millions base pairs sequenced) value for gene expression quantification, and differentially expressed genes (DEGs) between day 50 and day20, day80 and day 20, and day 80 and day50 were selected based on padj< 0.05. KEGG pathway analysis was conducted and corrected Pvalue (FDR) cut-off was set at 0.05. The sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) database under the accession number PRJNA588520.

Proteomic sequencing
Proteomic sequencing was performed by Novogene Bioinformatics Technology Co., Ltd (Beijing, China). Briefly, the concentrations of proteins were determined using Bradford assay and 100 µg of protein from each replicate was digested with Trypsin Gold (Promega, Madison, WI) at 37°C for 16 h. After that, peptides were dried by vacuum centrifugation, and Desalted peptides were labeled with TMT6/10-plex reagents (TMT6/10plex™ Isobaric Label Reagent Set, Thermofisher), following the manufacturer's instructions. After labeling, the samples were mixed and lyophilized. TMT-labeled peptide mix was fractionated using a C18 column (Waters BEH C18 4.6×250 mm, 5 5µm,) on a Rigol L3000 HPLC. Eluent was collected every minute and merged to generate 15 fractions.
Shotgun proteomics analyses were performed using an EASY-nLC TM 1200 UHPLC system (Thermo Fisher) coupled to an Orbitrap Q Exactive HF-X mass spectrometer (Thermo Fisher) operating in the data-dependent acquisition (DDA) mode. Peptides were separated on a Reprosil-Pur 120 C18-AQ analytical column (15 cm×150 μm, 1.9 μm). Full MS scans from 350 to 1500 m/z were acquired at a resolution of 60000 (at 200 m/z) with an AGC target value of 3×10 6 and a maximum ion injection time of 20 ms. From the full MS scan, a maximum number of 40 of the most abundant precursor ions were selected for higher energy collisional dissociation (HCD) fragment analysis at a resolution of 15000 for TMT6-plex (at 200 m/z) or 45000 for TMT10-plex (at 200 m/z) with an automatic gain control (AGC) target value of 1×10 5 , a maximum ion injection time of 45ms, a normalized collision energy of 32%, an intensity threshold of 8.3×10 3 , and the dynamic exclusion parameter set at 60 s. Differentially accumulated proteins (DAPs) were selected based on p< 0.05 and |log2FC| > 0.585. To identify the most important biochemical metabolic pathways and signal transduction pathways involved in differentially expressed proteins, KEGG analysis were performed. GO enrichment analysis of differentially expressed proteins was performed to identify the major biological functions.

Quantitative real-time PCR validation
Several genes which were significantly enriched in lignin, cellulose and xylan biosynthetic pathways were randomly selected for quantitative real-time PCR validation. Total RNA was isolated with TRIzol reagent (Takara, Dalian, Liaoning, China). First-strand cDNA synthesis was conducted by use of Prime Script TM RT reagent kits (Takara, Dalian, Liaoning, China). Quantitative real-time PCR was performed using SYBR Green kits (Takara, Dalian, Liaoning, China). Primer sequences were designed with the online designing tool in NCBI (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The reaction conditions of PCR were 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Expressions of genes were concluded by using the 2 −∆∆CT method.

Statistical analysis
Statistical significance of differences in this study was analyzed using one-way ANOVA followed by Tukey's post hoc test with a significance level of 0.05 (p < 0.05) (GraphPad Prism 7.0 software for Windows).   Expression of genes and proteins involved in lignin biosynthetic process. Heat maps were produced using standardized figures that were transformed to a value between 0.0 and 1.0 by Min-Max normalization method.

Figure 5
Other genes and proteins involved in apoplastic ROS productionSAM metabolic process and programmed cell death. Heat maps were produced using standardized figures that were transformed to a value between 0.0 and 1.0 by Min-Max normalization method. In the proteome data, the 0 value for all three periods representing that no corresponding protein was detected.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.