Integrity transcriptome and proteome analyses provide new insights into the mechanisms regulating pericarp cracking in Akebia trifoliata fruit


 Background: Akebia trifoliata (Thunb.) Koidz, a perennial wild woody liana, can be used as biofuel to generate bioenergy, as well as a traditional Chinese medicine plant,and new potential edible fruit crop, due to its high yields in fields, wide adaptability, high economic, medicinal and nutritive values, and tolerance tocultivation conditions. However, the pericarp of A. trifoliata cracks longitudinallyalong the ventral suture during fruit ripening, which is a serious problem that limits its usefulness and causes significant losses in yield and commercial value. Furthermore, there have been no known investigations on fruit cracking and its molecular mechanisms in A. trifoliata . Results: In this study, the dynamic structural changes in fruit pericarps were observed, revealing that the cell wall of fruit pericarp became thinner, and had reduced integrity, and that the cell walls began to degrade in the cracking fruits compared to those observed in non-cracking fruits. Moreover, analyses of the complementary RNA- sequencing-based transcriptomes and tandem mass tag-based proteomes at different development stages during fruit ripening were performed, and the expression of various genes and proteins was found to be changed after cracking, The mRNA levels of 20 differentially expressed genes and 17 differentially abundant proteins (DAPs) involved in cell wall metabolism were further analyzed; 20 DAPs were also validated through parallel reaction monitoring analysis. Among these, pectate lyases and pectinesterase involved in pentose and glucuronate interconversions, β-galactosidases 2 involved in galactose metabolism, were significantly up-regulated in cracking fruits compared to levels in non-cracking fruits, suggesting that they might play crucial roles in A. trifoliata fruit cracking. Conclusions: This study provides new insights into the molecular basis of fruit cracking in A. trifoliata fruits and important clues for further studies on the genetic improvement of A. trifoliata and the breeding of non-cracking varieties.


Changes in pericarp structure
The pericarp of A. trifoliata is known to crack longitudinally at maturity, and the seeds disperse from along the ventral suture with fruit cracking. First, the dynamic structures of the fruit pericarps in different development stages were observed in this study (Fig. 1). In the non-cracking stage (PS), the arrangements of pericarp cells and cuticles were dense, with small intercellular spaces, and were distributed continuously (Fig. 1a, 1d, 1g). However, the cell wall became thinner, the cell volume became larger, the number of cell layers decreased, and the arrangement of cells was loose with poor integrity; further, the spacing between the cells became bigger and the cell of exocarp and mesocarp began to degrade in the initial cracking stage (PM) (Fig. 1b, 1e, 1h). The cells were arranged irregularly and were not compact, and the cell layers continued to reduce, with larger cell spaces, and cells continued to degrade in the total cracking stage (PL) (Fig. 1c, 1f, 1i).

Transcriptomic analysis overview
To obtain an overview of the A. trifoliata transcriptome during fruit development and ripening, nine cDNA libraries (i.e. PS, PM and PL, each with three repeats) were constructed. A total of 47.05, 46.92, and 54.00 million raw sequence reads were produced from the PS, PM and PL libraries, respectively. After removing reads with indeterminate base ratios > 10 %, as well as low-quality reads and adaptor sequences, 46  Most of these were uncharacterized, unknown, and hypothetical proteins in the annotation results (Table  S2-S3).
Among these unigenes, 11 205 were identi ed as differently expressed genes (DEGs) using an absolute log 2 fold-change >1 with p < 0.05 during fruit ripening. There were 779 up-regulated and 1924 downregulated genes in the PM compared to levels in the PS group, which are presented in a volcano plot in Fig. S2a, and 4623 upregulated and 1975 downregulated genes in PL compared to expression in the PM group, and these are presented in a volcano plot in Fig. S2b. There were 1904 DEGs that were co-expressed in PM_PS and PL_PM combinations ( Table 2).

Functional classi cation of the identi ed DEGs
To further understand the function of the identi ed DEGs, bioinformatics analysis was performed based on gene functional classi cation and hierarchical cluster analysis. GO analysis indicated that most of the DEGs in the biological process (BP) category were involved in the cellular amide metabolic processes and amide metabolic processes; structural molecular activity and oxidoreductase activity comprised the highest proportions of DEGs in the molecular function (MF) category of both PM_PS and PL_PM groups, respectively; cytoplasmic parts and intracellular ribonucleoprotein complexes comprised the highest proportion of DEGs in the cell components (CC) in PM_PS and PL_PM, respectively ( Fig. 2a-2b).
Moreover, a KEGG pathway analysis was carried out to further evaluate the DEGs. Many were enriched in metabolic pathways, ribosomes in PM_PS, and biosynthesis of secondary metabolites and metabolic pathways in PL_PM ( Fig. 2e-2f). Additionally, a hierarchical cluster analysis was performed to further understand the expression changes in the cell wall-related DEGs, (Fig. 3). In all, 285 cell wall-related DEGs, including those associated with pentose and glucuronate interconversions, the phenylpropanoid pathway, galactose metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism and transcription factors were clustered closely both in PM_PS and PL_PM group. Notably, most cell wallrelated DEGs were downregulated in the PM_PS group but were upregulated in PL_PM group (Fig. 3a-3f).

Quantitative proteome analysis
To understand the molecular mechanisms of pericarp cracking in A. trifoliata fruits, a quantitative proteomics analysis was performed using the TMT platform and LC-MS/MS analysis during fruit ripening, to complement the transcriptome analysis. Accordingly, totals of 812 625 spectra, 68 151 identi ed spectra, 12 456 peptides, and 10 572 unique peptides were found by proteomic analysis and 2839 proteins were identi ed (Table 1 and Table S3). In terms of protein mass distribution, proteins with molecular weights greater than 9 kDa had a wide range and good coverage, with a maximum distribution area of 10-40 kDa (Fig. S3a). Peptide quantitative analysis of the proteins showed that protein quantity decreased with an increase in the matching peptide (Fig. S3b).
Among of these proteins, 240 were identi ed as differentially abundant proteins (DAPs) using a foldchange >1.2 and < 0.83 with p < 0.05 as the thresholds for upregulated and downregulated, respectively. Further, 84 proteins were more abundant and 106 proteins were less abundant in the PM_PS group and are shown in a volcano plot in Fig. S2c; 20 DAPs were more abundant and 13 DAPs were less abundant in the PL_PM group, and these are shown in a volcano plot in Fig. S2d, whereas 17 were co-expressed in PM_PS and PL_PM.
Functional classi cation of the identi ed DAPs Bioinformatics analysis of DAPs was carried out based on protein functional classi cations and hierarchical cluster analysis. GO analysis showed that most DAPs in the BP category were involved in cellular responses to the chemical stimulus and cellular oxidant detoxi cation processes in the PM_PS group, as well as the metabolic and macromolecule metabolic processes in the PL_PM. The highest portions of DAPs in MF category were the oxidoreductase activity and antioxidant activity in PM_PS, as well as structural constituent of ribosome, and structural molecule activity in PL_PM group. Extracellular region in the PM_PS group and cytoplasmic parts in the PL_PM group comprised the highest portions of DAPs in CC category (Fig. 2c-2d).
Moreover, a KEGG pathway analysis was carried out to further evaluate the DAPs. Many were enriched in two-component system and ribosome pathways in the PM_PS and PL_PM groups, respectively ( Fig. 2g-2h).
Hierarchical cluster analysis was performed to further explore the expression changes in the cell wallrelated DAPs. A total of 40 cell wall-related DAPs, including those associated with pentose and glucuronate interconversions, the phenylpropanoid pathway, galactose metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and cell wall metabolism-related proteins were clustered closely in both PM_PS and PL_PM groups. Notably, most cell wall-related DAPs were upregulated in both the PM_PS and PL_PM groups, whereas those DAPs involved in phenylpropanoid pathways and galactose metabolism were downregulated in the PM_PS group and up-regulated in the PL_PM group ( Fig. 3g-l).
Comparative analysis between protein abundance and gene expression levels There were more DEGs (4607 and 6598) than DAPs (190 and 50; PM_PS and PL_PM, respectively) and many more shared DEGs (1904) than shared DAPs (17) in cracking fruit compared to those in noncracking fruit. Most of these DEGs and DAPs were downregulated in the PM_PS group but upregulated in the PL_PM group, suggesting that greater changes in fruit cracking occurred during fruit ripening. Among the shared DEGs, 1123 were upregulated and 781 were downregulated in the PM_PS group, whereas 808 were upregulated and 1096 were downregulated in PL_PM group. Of the shared DAPs, nine were increased in abundance and eight were decreased in abundance the PM_PS group, whereas eight exhibited increased abundance and nine showed decreased abundance in the PL_PM group (Table 2). To evaluate the relationships between the transcriptomic and proteomic changes during fruit pericarp cracking, the quantitative data for DEGs and DAPs were used for correlation analysis. According to this, 14 and four DAPs and their corresponding DEGs were identi ed in the PM_PS and PL_PM groups, respectively. Of these, 12 DAPs (four with increased abundance and eight with decreased abundance) and four DAPs (two with increased abundance and two with decreased abundance) were regulated in the same direction as their corresponding DEGs in the PM_PS and PL_PM groups, respectively ( Fig. 4a-4b); meanwhile, two were regulated in the opposite direction as their corresponding DEGs in the PM_PS group. There were more DEGs than DAPs in both the PM_PS and PL_PM groups, with signi cant differences in the trends of transcript levels and protein abundance.

Identi cation of DAPs and DEGs associated with candidate pathways
To further clarify the biological functions of the co-regulated DEGs-DAPs genes, an enrichment analysis was conducted based on the GO and KEGG pathways analyses. The largest groups within the BP category were those linked to metabolic process and cellular process; catalytic activity and binding and cell and cell part were predominant in MF and CC categories, respectively, both in PM_PS and PL_PM groups ( Fig. S4a-4b). In the PM_PS group, 14 DAPs were signi cantly enriched in seven pathways, both with respect to DEGs and DAPs, which included fructose and mannose metabolism pathway, phenylpropanoid biosynthesis, glutathione metabolism, ubiquinone and other terpenoid-quinone biosynthesis, pentose and glucuronate interconversions, amino sugar and nucleotide sugar metabolism, and galactose metabolism. In the PL_PM group, four DAPs were signi cantly enriched in the three pathways, both in terms of DEGs and DAPs, which included a calcium signaling pathway, pentose and glucuronate interconversions, and galactose metabolism pathway (Fig S4c-4d). The comparative analysis showed that two pathways, including pentose and glucuronate interconversions and galactose metabolism pathways, were shared between the transcriptome and proteome data, for both PM_PS and PL_PM groups. However, the phenylpropanoid biosynthesis pathway was only shared by DAPs and DEGs in the PM_PS group.
Therefore, those shared metabolic pathways might play potential roles in A. trifoliata fruit cracking.

Validation of data reliability by parallel reaction monitoring (PRM)
The protein expression levels obtained by TMT were con rmed by quantifying their expression by PRM analysis. 20 proteins that exhibited signi cantly different levels by TMT analysis were selected for PRM analysis, of which 18 were successfully quanti ed (Table 4). Here, 14 of the 18 (77.8%) proteins showed the same trend as that observed when the protein levels were quanti ed by TMT, including PE, PL, PG2, F26G, β-GAL2, Auxin e ux carrier, α-HY, PRX2, PG4, PRX5, PRX3, endoglucanase 19, endoglucanase 8, and DIR1. Meanwhile the mean expression levels of PRX, beta-fructofuranosidase, BXL, and BGLU33 proteins were inconsistent with the protein levels quanti ed by TMT. In general, the trends in the expression changes measured by PRM and TMT were basically consistent.

Discussion
Structural changes in the pericarp cell wall might affect A. trifoliata fruit cracking Fruit cracking is a complex phenomenon that is caused by a series of environmental, physiological, biochemical, and genetic changes during fruit ripening. Fruit cracking happens when the stress exerted on the pericarp from the enlarged aril is greater than the strength of the fruit skin, and the mechanical strength of the pericarp depends largely on its cell wall [28]. Studies showed that jujube fruit cracking might be related to the changes of cell wall structure and the rearrangement of the cell wall at the later stages of fruit ripening [29]. Arrangement of the subcutaneous layers of the cells was found to be relatively regular, and cell layers had a closer arrangement in the cracking-resistant tomato genotype [30]. In this study, the pericarp cell structure and ultrastructure in A. trifoliata fruits of different development stages were observed. Compared to those in the non-cracking fruit, the cell wall structures of the fruit pericarp had poor integrity, loose cell wall structures, deformed and reduced cell layers, and larger spaces, and the pericarp cells began to degrade during fruit cracking ( Fig. 1e and 1h); these changes were more pronounced during the total cracking stage ( Fig. 1f and 1i), which was consistent with previous results in grapes and orange [31][32]. These results indicated that the structural changes in the cell wall of A. trifoliata pericarp might play a key role in the occurrence of fruit cracking.
General features of the transcriptomes and proteomes of different A. trifoliata pericarps Fruit cracking is a key factor that affects marketability of fruits, reduces their market acceptability, and causes signi cant losses in eld yield and commercial value [33]. Elucidating the molecular mechanisms that regulate fruit cracking could aid in the utilization of A. trifoliata for biofuels. However, there have been no speci c studies on A. trifoliata fruit cracking and the underlying mechanisms are still largely unknown. In this study, the differences in the transcriptome and proteome were investigated based on RNA-seq and TMT data during different development stages. As the transcriptome database was used for protein identification in this study, the quality of the sequencing and assembly of the transcriptome data was crucial for subsequent analyses. A total of 186 054 unigenes (> 200 bp) were assembled in the transcriptome of A. trifoliata pericarp, which is much more than that previously obtained for this species and those of other Ranunculales such as A. trifoliata (11 749 by Yang et al. [25]; 65 757 by Niu et al. [1]), Dysosma aurantiocaulis (53 929) [34], and Dysosma versipellis (44 855) [35]. The obtained percentages of the Q30 bases (91. 58 Moreover, a total of 812 625 spectra, 10 572 unique peptides, and 2839 proteins were identi ed based on transcriptome data from A. trifoliata pericarp. The RNA-seq and protein sequencing methods identi ed and annotated many genes and proteins, providing the basis for a more precise and detailed description of molecular processes, as well as the elucidation and better understanding of complex physiological processes and their genetic regulation [36]. It can be seen that the A. trifoliata pericarp sequencing and assembly quality presented here were high enough and useful tools for future genetic research on fruit cracking in A. trifoliata and other Lardizabalaceae species. Additionally, there are many 'uncharacterized', 'predicted', or 'putative' transcripts and proteins in the annotation results that might be limited by the lack of genomic information; these were not associated with a de nitive annotation and might also exist in other species [37]. Thus, the roles of these unknown or uncharacterized genes and proteins in A. trifoliata fruit cracking constitute an important issue for future research. Translational and posttranslational regulation of A. trifoliata fruit cracking Some inevitable passive processes will occur in the cracked pericarp following fruit cracking, such as oxidative stress and microbial invasion [37]. Therefore, the differences in the expression of genes and proteins among PS, PM, and PL alone cannot accurately re ect the cause of fruit cracking. In this study, a comparative analysis of protein abundance and gene expression levels was performed, and more DEGs (11205) than DAPs (240), as well as more shared DEGs (1904) than shared DAPs (17), were identi ed in cracking fruits compared to those in non-cracking fruits. A possible explanation for the much lower number of DAPs is the limitation of the proteomic method of choice (MS/MS) in identifying the presence and abundance of the proteins, as compared to that with transcriptomics (RNA-seq) [38]. Moreover, a different trend was observed between protein and transcription levels. The integration of transcriptomic and proteomic data revealed that only 14 and four DAPs corresponding to DEGs and a poor correlation (r = 0.03 and 0.11) between the protein abundance and the corresponding gene expression were identi ed in the PM_PS and PL_PM groups, respectively. These results indicated that there was discordance between the transcript levels and protein abundance, which was similar to that observed in previous reports, suggesting that post-transcriptional and post-translational regulation, reversible phosphorylation, splicing events in cells, and translation e ciency might also play key roles in the regulation of fruit ripening [39][40].
Moreover, the PPI analysis indicated that the interactions of proteins were weak in this study, and only four proteins associated with cell wall metabolism were determined to interact with other proteins. This was associated with the results indicating that the interactions are often weak for many cellular processes, which are regulated by post-translational modi cations that are recognized by speci c domains in protein binding partners [41]. This inconsistency showed that the regulation from mRNA to proteins is a complex process, and the changes in protein abundance were generated after its corresponding transcript was stabilized [42]. Thus, gene translation and post-translation processes might be important regulatory factors during A. trifoliata fruit ripening and cracking, which was consistent with the result that many unigenes were assigned to the category of posttranslational modi cation [35].
Potential regulators and metabolism pathways during fruit cracking GO functional enrichment and KEGG pathways provide prediction information of inner-cell metabolic pathways, as well as the genetic and biologic behaviors of genes. The GO enrichment analysis indicated that most DEGs and DAPs were both involved in oxidoreductase activity and structural molecule activity in MF category in the PM_PS and PL_PM groups, respectively ( Fig. 2a-2d). These results were associated with previous studies and could classify the target genes into different categories based on their functions and amounts, such as peroxidase and cell wall polysaccharide [18,43]. Cell polysaccharides are degraded by the cell wall hydrolases, and the formation of phenolic cross-linking cell wall structural components is catalyzed by cell wall peroxidase. These modi cations reduce the strength of the fruit pericarp, resulting in changes in pericarp development and fruit cracking [44]. In this study, the hierarchical cluster analysis of cell wall-related DEGs and DAPs indicated that most were upregulated in the PL_PM groups both in the transcriptome and proteome (Fig. 3). The pathway enrichment analysis indicated that two cell wall-related pathways, including the pentose and glucuronate interconversions and galactose metabolism pathway, were common pathways shared by the DAPs and DEGs in both the IS_NS and TS_IS groups. Further, the phenylpropanoid biosynthesis pathway was shared by the DAPs and DEGs in the IS_NS group, suggesting that cell well metabolism might have an important role in A. trifoliata fruit ripening and cracking. Moreover, these results revealed that the proteomic data and the transcriptome data were complementary, and that the proteome could con rm the transcriptome data; in addition, genes perform the same function at the transcriptome and proteome levels [42]. The functional classi cation and integrative analysis of the transcriptome and proteome would provide a good basis to better understand the molecular physiology of fruit ripening and cracking.

Candidate DEGs and DAPs might play key roles in fruit pericarp cracking
Fruit cracking is a complex physiological phenomenon that is controlled by many genes working together, rather than a single gene directly controlling the process [16,45]. Researchers have found that several cell wall modi cation genes, including β-GAL, β-GLU, PE, and PG are differentially expressed in cracked fruits compared to levels in non-cracked litchi fruits [46]. The suppression of PpBGAL expression can reduce PpPG transcription and activity, suggesting that the downregulation of PpBGAL expression delays peach fruit softening [47]. Antisense inhibition of PE and PG activity in tomato was also found to reduce fruit cracking [48]. Further, silencing the SIPL gene can enhance fruit rmness and reduce the content of pectin, suggesting that this gene participates in the pericarp cell wall rearrangement during fruit softening [49]. Fruit softening is mainly caused by hemicellulose and pectin degradation proteins, including XYL, BGAL, PE, and PG [50][51]. As a major cellulose degradation enzyme, BGLU is expressed in ripening fruits, and the downregulation of this protein in strawberry can delay fruit maturation [52]. PRXs are well-known cell wall enzymes and involved in the rearrangements of cell wall polysaccharides during plant development [53].
In this study, 20 DEG and 17 DAPs involved in cell wall metabolism were validated using qPCR, and 27 of 37 (13 in DEGs and 14 in DAPs) showed strong correlations with protein expression levels. Moreover, 20 DAPs involved in cell wall metabolism were also validated using PRM and 14 showed the same trend as TMT protein levels, which was associated with the results that cell polysaccharide metabolism plays key roles in fruit ripening and cracking [43,45]. Notably, three proteins, including PE, PL, and β-GAL2, were functionally annotated among thousands of identi ed DEGs and DAPs, and the expression of these three enzymes was signi cantly upregulated in cracking fruit comparing to that in non-cracking fruit. PE, involved in the pentose and glucuronate interconversions pathway, and β-GAL2, involved in galactose metabolism, were upregulated in the PL_PM group; meanwhile, PL, involved in the pentose and glucuronate interconversions pathway was upregulated in the PM_PS group with the same expression trends observed in transcriptome, proteome, qPCR, and PRM data (Fig. 8), suggesting that PE, PL, and β-GAL2 might play important roles in the regulation of A. trifoliata fruit cracking. Moreover, F26G, PG, PG3, XYL, PRX3, and PRX5 also showed higher protein and gene expression in cracking fruit than in noncracking fruit. The signi cantly increased expression of these cell wall metabolism proteins is in accordance with the results indicating that a network of cell wall cellulose, hemicellulose, and pectin might be correlated with fruit cracking, suggesting their dynamic roles in fruit ripening and cracking [52,54].

Conclusions
This study revealed the structural changes in the cell wall during different fruit development stages of A. trifoliata pericarp, suggesting that the structural changes between unripe and ripe fruit might be an important factor in fruit tendency to crack. Therefore, bagging treatment could be performed to strengthen the structure of the pericarp and reduce fruit cracking [55]. Moreover, this study presented comprehensive transcriptome and proteome data to screen the potential effectors of fruit cracking in A. trifoliate, and various genes and proteins were found to be differentially expressed after cracking, providing a global view of the molecular mechanism of A. trifoliate fruit cracking. The co-expressed genes and proteins were found to be implicated in various signaling pathways in fruit development, and those genes involved in cell wall metabolism might play important roles in A. trifoliate fruit ripening and cracking. Therefore, cell wall structure and modi cations might be correlated with the strength of the pericarp and fruit cracking. Moreover, degradation of the cell wall also depends on ethylene during fruit softening. For example, the activities of β-GAL and PG could be suppressed in 1-methylcyclopropene (1-MCP)-treated fruit in response to ethylene treatment [56]. Therefore, 1-MCP treatment could be performed to delay fruit cracking during fruit ripening. In general, these omics data might provide a new perspective to further study the molecular mechanism of fruit cracking and the breeding of non-cracking varieties of Lardizabalaceae and other species. The CRISPR/Cas9 system could be used to determine the molecular functions of these codifferentially expressed proteins/genes in further studies to reveal the mechanisms of fruit cracking in Lardizabalaceae.

Plant materials
The wild germplasm Nong No.8, which had been transplanted for 9 years in the nursery at Hunan Academy of agricultural sciences, Changsha, P. R. China, was used as the research material in this study. According to our observations, in Changsha, Hunan Province, the blooming stage for the germplasm of Nong No.8 was in early April, when 50 % of the A. trifoliata owers were in bloom with a owering period of approximately 30-45 d. The fruit development lasts for 5 months, from the ovary in ation through the fruit setting, and longitudinal and transverse elongation, esh softening, peel cracking, complete the development process [57]. The Nong No.8 variety usually ripens in early October. Since the fruit of A. trifoliata do not develop uniformly, fruits were harvested from different stages of the same Nong No.8 tree at the same time and were then sorted according to their developmental stage, considered to be a sample.
Different fruits, including the non-cracking fruits (PS), the initial cracking fruits (PM), and the total cracking fruits (PL) were randomly taken every 10 days at the ripening stage (September 18, September 28 and October 8, 2018) and three fruits were mixed into one biological replicate for further analysis, and in total, three biological replicates were collected for each stage. Sampled pericarps were rapidly collected and immediately frozen in liquid nitrogen and stored at −80 ºC until use for transcriptome, proteome, qPCR, and PRM analyses.
Anatomical structure of pericarp Anatomy of the pericarp samples taken from the Nong No.8 fruits at different development stages, including PS, PM, and PL, were prepared for para n sections, and scanning electron microscopy (SEM) was carried out according to previous studies [58][59]. Pericarp samples were fixed directly in the field using FAA [70 % ethyl alcohol + 38 % methyl aldehyde + 25 % acetic acid (16:1:1)]. Then the tissues were subsequently dehydrated through an ethanol series with increasing ethanol concentrations and embedded in para n. Subsequently, para n sections were stained with Safranin O and fast green staining and observed with an Axio Imager (Zeiss, Oberkochen, Germany); upright microscopy and images were displayed using Image-pro Plus 6.0 software. RNA isolation, library construction, and sequencing Total RNA used for the RNA-seq assays was isolated from three independent replicates of pericarp in the PS, PM, and PL stages, as described by Tao et al [60]. The RNA samples were detected based on the A260/A280 absorbance ratio with a Nanodrop ND-1000 system (Thermo Scienti c). Paired-end libraries were prepared using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's instructions. The mRNA was puri ed from 3 µg of the total RNA using oligo (dT) magnetic beads followed by fragmentation carried out using divalent cations at elevated temperatures in NEBNext First Strand Synthesis Reaction Buffer. Subsequently, rst strand cDNAs were synthesized with random hexamer primers and Reverse Transcriptase (RNase H-) using mRNA fragments as templates, followed by the second strand cDNA synthesis using DNA polymerase I, RNAseH, buffer, and dNTPs. The synthesized double-stranded cDNA fragments were then puri ed with an AMPure XP system (Beckman Coulter, Beverly, USA). The puri ed double-stranded cDNA was polyadenylated and adapter-ligated for preparation of the paired-end library. Adaptor-ligated cDNA and adaptor primers were used for PCR amplification. PCR products were puri ed (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Finally, sequencing was performed with an Illumina HiSeq2500 instrument by Shanghai Applied Protein technology (Shanghai, China).

Quality control and transcriptome assembly
The raw paired-end reads in fastq format produced from the sequencing were rst processed with the inhouse Perl scripts. Those reads containing adapters, excess "N" nucleotides with more than 10 % of the bases and reads of low-quality (reads with quality values ≤ 10) were removed by lter_fq software. The Q20, Q30, GC-content, and sequence duplication levels of the obtained clean reads were calculated. The assembly of clean reads to unigene collections was performed using the Trinity software package (https://github.com/trinityrnaseq/trinityrnaseq/releases) [61]. The Trinity software consists of three independent software modules, including Inchworm, Chrysalis, and Butterfly, and the transcripts less than 200 bp in length were discarded. Sequences containing the longest cluster transcripts without redundancy extracted from transcripts can be considered unigenes.

Bioinformatics analyses
The de novo assembled unigenes were annotated in ve databases, which include NR, Pfam, the Swiss-Prot, GO, and KEGG pathway databases, based on a BLAST search with an E-value threshold of 1E −5 .
Moreover, to further analyze the annotation results, GO and KEGG results with E-values of 1E −5 were used for functional gene annotation. GO terms could be classi ed into three categories, including BP, MF, and CC. In addition to the GO terms, the pathway maps were determined based on the KEGG database.
The normalized transcript abundances of the genes were estimated using the FPKM based on the length of the gene and read counts mapped to this gene. DESeq2 R package (1.16.1) software was used to identify the differential expression of the genes (DEGs), and the false discovery rate (FDR) was controlled using the Benjamini and Hochberg's approach to adjust the P-value. Genes with an adjusted P-value < 0.05 and absolute fold-change of 2 were deemed to be differentially expressed between the two samples. In addition, GO and KEGG pathway enrichment analysis of DEGs was implemented with the clusterPro ler R package. Transcription factor analysis of DEGs was performed against the PlantTFDB database (http://planttfdb.cbi.pku.edu.cn/). The heat map was visualized using heatmap 2.0 in the gplot R package.

Protein extraction
Protein extraction from A. trifoliata pericarp was performed from each sample as described previously [62]. The samples were frozen in liquid nitrogen and ground into powder. A ve times volume of TCA/acetone (1:9) was added, and the sample was vortexed and mixed and placed at −20 ℃ for 4 h; it was then centrifuged at 6000 × g at 4 ℃ for 40 min. The supernatant was discarded, and the precipitate was washed three times with pre-cooled acetone. After the precipitation and air drying, the precipitate was redissolved in buffer (4 % SDS, 100 mM tris-HCl and 1 mM DTT; pH 7.6). After sonication and boiling for 15 min, the lysate was centrifuged for 40 min, and the supernatant was ltered and quanti ed using the BCA Protein Assay Kit (Bio-Rad, USA).

Trypsin digestion and TMT labeling
For digestion, the samples were added to the buffer (4 % SDS, 100 mM DTT, 150 mM tris-HCl, pH 8.0) and UA buffer (8 M urea, 150 mM tris-HCl pH 8.0) by repeated ultra ltration (Microcon units, 10 kD). Then, iodoacetamide (100 mM IAA in UA buffer) was added into the samples to block reduced cysteine residues, incubating for 30 min in darkness at room temperature. After the lters were washed with UA buffer and triethylamine borane (TEAB) buffer in turn, the suspensions were digested with trypsin (Promega, Madison, WI) in TEAB buffer overnight. After trypsin digestion, the samples (100 μg of protein) were categorized to label them with 129-tag (PS), 130-tag (PM), and 131-tag (PL) (Thermo Fisher Scienti c, Waltham, MA, USA). Finally, TMT-labeled peptide aliquots were pooled for subsequent fractionation using the Pierce high pH reversed-phase fractionation kit (Thermo scienti c).

HPLC fractionation and LC-MS/MS analysis
For the fractionation of labeled peptides, samples were loaded onto a reverse phase trap column (Thermo Scienti c Acclaim PepMap100, 100 μm × 2 cm, nanoViper C18) connected to the C18-reversed phase analytical column (Thermo Scienti c Easy Column, 10 cm long, 75 μm inner diameter, 3 μm resin) in buffer A (0.1 % formic acid) and separated with a linear gradient of buffer B (84 % acetonitrile and 0.1 % formic acid) at a ow rate of 300 nl/min controlled by IntelliFlow technology. The resultant peptides were further processed using a Q Exactive mass spectrometer (Thermo Scienti c) that was coupled to Easy nLC (Thermo Fisher Scienti c). Mass spectrometry analysis was performed in positive ion mode, and MS data were acquired using a data-dependent top10 method dynamically choosing the most abundant precursor ions from the survey scan (300-1800 m/z) through the higher-energy collision dissociation (HCD) fragmentation method. Automatic gain control (AGC) was set as 3E6 and the maximum inject time was 10 ms, with the following parameters: 40.0-s dynamic exclusion duration; 70,000 resolutions with survey scans at m/z 200 and resolution for HCD spectra at 17500 at m/z 200; 2 m/z of isolation width; 30 eV of normalized collision energy and the under ll ratio was de ned as 0.1 %.

Sequence database search and data analysis
The raw data were processed by the MASCOT engine (Matrix Science, London, UK; version 2.2), and Proteome Discoverer 1.4 software was used to process MS/MS spectra. The search was performed using the following settings based on the A. trifoliata database: trypsin for the enzyme and 2 as the maximum missed cleavage allowed; fixed modifications of carbamidomethyl (C), TMT-6plex (N-term), and TMT-6plex (K), variable modification of oxidation (M); mass tolerance for fragment ions of 0.1 Da, and 20 ppm for peptide ions, as well as both peptide and protein levels of FDR less than 0.01, and only unique peptides of the protein were employed for the protein identification and quantification. Proteins with a P value less than 0.05 and fold-change ≥ 1.2 or ≤ 0.83 within a comparison were recognized as DAPs.
Functional categorization was performed using GO and KEGG pathway databases with P values ≤ 0.05. The protein functional network was performed with STRING 9.0 software (http://string-db.org). Clustering analysis of the DEPs was performed using Cluster 3.0 (http://bioservices.capitalbio.com/xzzq/rj/3885.shtml) and the Java Treeview software (http://jtreeview.sourceforge.net). Correlations were analyzed based on the DEGs and DEPs, and Person correlation tests were conducted for each comparison group, including IS vs NS and TS vs IS.

Reverse transcription real-time quantitative PCR
The method of total RNA extraction and synthesis of cDNA were described previously. The Bio-Rad CFX96 Touch detection system (Bio-Rad, Richmond, CA, USA) with SYBR Green PCR master mix (Aidlab Biotechnologies, Co., Ltd) were used for the reactions of each sample. In this study, the EF-1α gene was used as the internal control gene, which was detected by the de novo transcriptome sequencing of A. trifoliata [63]. Primers for qPCR experiments were designed using Primer 5.0 software (Supplementary Table 1) and those gene sequences were blasted against the NCBI database. The ampli cation reactions contained 12.5 μL SYBR Green PCR master mix, 1 μL cDNA, and 0.5 μL of each primer in a nal reaction volume of 25 μL. The thermal cycling program began with 3 min at 95 °C, followed by 40 cycles of 95 °C for 10 s and 55 °C for 30 s, with a melt curve from 65 to 95 °C based on increments of 0.5 °C for 5 s. After PCR ampli cation, the quantitative variation was analyzed using the Delta Ct method, and the analysis of statistically signi cant differences from gene expression was performed by the independent samples t-test analysis at P < 0.05 using GraphPad Prism 8 software. Correlation analysis between cell wall-related gene and protein expression was performed by Pearson's correlation coe cient analysis [64].

PRM analysis
To verify the protein expression levels obtained by TMT analysis, the expression levels of selected proteins were further quanti ed by PRM analysis. The AQUA stable isotope peptide was spiked in each sample as an internal standard reference. The tryptic peptides were directly loaded on C18 STAGE-tips for desalting prior to reversed-phase chromatography on an Easy nLC-1200 system (Thermo Scienti c). Acetonitrile was used in 45 min from 5% to 35% for a 1-h liquid chromatography gradient, and PRM analysis was performed on a Q Exactive TM Plus mass spectrometer (Thermo Scienti c). The mass spectrometer was operated in positive ion mode and the full MS1 scan was acquired with a resolution of 60 000 (at 200 m/z). AGC target values and maximum ion injection times were set at 3e6 and 200 ms, respectively; full MS scans were followed by 20 PRMs (MS2 scans) at 30 000 resolution (at m/z 200) with AGC of 3e6 and the maximum injection time set as 120 ms. The targeted peptides were isolated with a 2 Th window. Ion activation/dissociation was performed at a normalized collision energy of 27 in the HCD collision cell. Skyline version 3.7.0 was used to analyze the MS data where signal intensities for individual peptide sequences for each of the signi cantly altered proteins were quanti ed relative to each sample and normalized to a standard reference [65].  The morphology and structure Changes of 'Nong No.8' pericarp at different development stages. a-c The morphology changes of A. trifoliate pericarp in non-cracking stage (PS) (Fig. 1a), initial cracking stage (PM) (Fig. 1b) and total cracking stage (PL) (Fig. 1c), respectively. d-f Pericarp structure of 'Nong No.8' at different development stages using the semi-thin slices method and stained with Safranin O-staining in PS ( Fig. 1d), PM (Fig. 1e) and PL (Fig. 1f). g-i Pericarp structure of 'Nong No.8' at different development stages using the scanning electron microscope method in PS (Fig. 1g), PM (Fig. 1h) and PL (Fig. 1i)   Red indicates signi cantly upregulated proteins, and blue indicates signi cantly downregulated proteins.

Figure 6
Validation and expression analysis of selected genes using qPCR. The expression levels of the genes revelaed by RNA-seq (Left y-axis) and qPCR (right y-axis). Histograms were gene expression detected by RNA-seq. Line graphs were relative expression validated by qPCR.

Figure 7
Validation and expression analysis of selected proteins using qPCR. The expression levels of the genes revelaed by TMT (Left y-axis) and qPCR (right y-axis). Histograms were protein abundance detected by TMT. Line graphs were relative expression validated by qRT-PCR