Changes in pericarp structure
The pericarp of A. trifoliata is known to crack longitudinally as it matures. The crack in the fruit and subsequent dispersal of the seeds from along the ventral suture, which is similar to the cracking seen in other legume species, so the dynamic structures of the fruit peels in different development stages were observed (Fig. 1). In the non-cracking stage (PS), the arrangement of pericarp cell and cuticle were dense, small intercellular space and distributed continuously. The fruit cracking however, was accompanied by the cells becoming thinner and bigger, a reduction in the number of cell layers, and furthermore, the arrangement of the cells was loose and poor integrity, and they began to degrade in the initial cracking stage (PM). The cells were arranged irregularly and continue to reduce cell layers, and there was atrophy degradation in the total cracking stage (PL).
Transcriptomic analysis overview
To obtain an overview of the A. trifoliata transcriptome during fruit development and ripening, three cDNA libraries (ie., PS, PM and PL) 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 of > 10%, low-quality reads and adaptor sequences, 46.45, 46.35, 53.46 million clean reads with the percentage of Q30 bases and GC contents of 91.58 − 93.75% and 45.94%−48.48%, were obtained respectively (Table S1). The resultant A. trifoliata transcriptome contains 241 376 transcripts, ranging from 201 to 2000 bp, and a total of 186 054 unigenes were identified (Table 1), and the details of size distribution of the transcripts and unigenes are shown in Fig. S1.
Table 1
Summary of the transcriotome and proteome data in Akebia trifoliata fruits.
RNA-seq data | MS data based on transcriptome |
Total number of transcripts | 241,376 | Total spectra | 812625 |
Mean length of transcripts (bp) | 515 | Identified spectra | 68151 |
Total number of unigenes | 186,054 | Identified peptides | 12456 |
Mean length of unigenes (bp) | 447 | unique peptides | 10572 |
N50 length of transcripts (bp) | 713 | Identified proteins | 2839 |
N50 length of unigenes (bp) | 518 | | |
To determine the putative functions of the assembled transcripts, all unigenes were annotated using Basic Local Alignment Search Tool (BLAST) searches against the five databases, including National center for biotechnology information non-redundant protein sequences database (NR) (100 329, 100%), SwissProt (56346, 56.16%), Protein families database (Pfam) (34428, 34.32%), Gene Ontology database (GO) (44558, 44.41%), and Kyoto Encyclopedia of Genes and Genomes pathway database (KEGG) (30298, 30.20%); indicating that the NR database provided the largest number of annotations. 100 329 unigenes were shown to correspond with sequences from at least one of the public databases, and 7283 unigenes annotated in all databases, resulting in the annotation of 100,329 unigenes (41.57% of all the assembled transcripts) in A. trifoliata pericarp (Table S2).
Among of these unigenes, 11 205 were identified as differently expressed genes (DEGs) using absolute log2 fold change > 1 with p < 0.05 during the fruit ripening. There were 779 up-regulated and 1924 down-regulated unigenes in the PM compared with the PS group, which are presented in a volcano plot in Fig. S2a; 4623 up-regulated and 1975 down-regulated in PL, compared with the PM group, and are presented in a volcano plot in Fig. S2b. There were 1904 DEGs that were co-expressed at PM_PS and PL_PM (Table 2).
Table 2
Summary of proteins and transcripts detected from TMT and RNA sequence data.
| Protein | Transcriptome |
PM_PS | PL_PM | PM_PS | PL_PM |
Unique proteins/genes detected | 2839 | 2839 | 100329 | 100329 |
Significantly DAPs /DEGs | 190 | 50 | 4607 | 6598 |
Up-regulated | 84 | 28 | 1945 | 4623 |
Down-regulated | 106 | 22 | 2662 | 1975 |
Shared proteins/genes | 17 | 1904 |
Shared proteins/genes (up-regulated) | 9 | 8 | 1123 | 808 |
Shared proteins/genes (down-regulated) | 8 | 9 | 781 | 1096 |
Functional classification of the identified DEGs
To further understand the function of the identified DEGs, bioinformatics analysis was performed on the basis of gene functional classification and hierarchical cluster analysis. GO analysis indicated that most of the DEGs in the biological processes were involved in the cellular amide metabolic processes and amide metabolic processes; structural molecular activity and oxidoreductase activity were the highest portion of the DEGs in the category of molecular functions, both in PM_PS and PL_PM; cytoplasmic parts and intracellular ribonucleoprotein complexes in PM_PS cells and PL_PM cell parts, were the highest portion of the DEGs in category of cell components respectively (Fig. 2a-2b).
Moreover, a KEGG pathway analysis was carried out to further evaluate the DEGs. In the comparison between the PM_PS and PL_PM, many DEGs were enriched in metabolic pathways, ribosomes, and biosynthesis of the secondary metabolites. Notably, most of the DEGs involved in the cell wall-related DEGs, were downregulated in the PM, compared with the PS group, but however they were upregulated in the PL compared with the PM group (Fig. 2e-2f).
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 were clustered closely both in PM_PS and PL_PM group. Most of these were involved in pentose and glucuronate interconversions, phenylpropanoid pathway, galactose metabolism, starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism and transcription factors (Fig. 3a-3f).
Confirmation of fruit cell-wall related genes in A. trifoliata by Reverse transcription real-time quantitative PCR (qPCR)
To validate the results of the RNA-Seq data, an expression and correlation analysis between the qPCR and the fragments per kilobase per million reads mapped (FPKM) values obtained from the RNA-seq were performed. The 20 selected genes had shown differential expression patterns in the PM_PS and PL_PM groups, and the results of the qPCR are shown in Fig. 4. Specifically, in the PM_PS group, most of these DEGs were down-regulated, including the phenylpropanoid pathway related genes- 4-coumarate-COA-ligase (4CL), peroxidase (PRX), and PRX2; galactose metabolism related genes-β-galactosidases (β-GAL1 and β-GAL2); amino sugar and nucleotide sugar metabolism related gene-beta-D-xylosidase (BXL), starch and sucrose metabolism related genes- cellulase (CEL), cellulose synthase-like protein (CSLG) and Glucan endo-1,3-beta-D-glucosidase (ENDOB), and the transcription factor and cell wall metabolism genes-NAC, NAC like and EXP1. While the phenylpropanoid pathway related genes- cinnamyl-alcohol dehydrogenase (CAD) and shikimate O-hydroxycinnamoyltransferase (HCT), transcription factor and cell wall metabolism genes-BHLH and dirigent protein (DIR2), and pentose and glucuronate interconversion related genes (PL, PG and PE) were significantly up-regulated. In the PL_PM group, most of these genes were significantly up-regulated in the PL_PM group, except for 4CL, CAD, β-GAL, and EXP1. Moreover, the expression of 13 candidate genes, including DIR2 (r = 0.7977, p < 0.05), NAC like (r = 0.9464, p < 0.01), EXP1 (r = 0.8582, p < 0.01), CAD (r = 0.8015, p < 0.01), β-GAL1 (r = 0.8440, p < 0.05), β-GAL2 (r = 0.6675, p < 0.05), 4CL (r = 0.7466, p < 0.05), ENDOB (r = 0.7052, p < 0.05), PE (r = 0.8042, p < 0.05), BHLH (r = 0.7485, p < 0.05), PG3 (r = 0.6819, p < 0.05), CEL (r = 0.8732, p < 0.05), and PRX2 (r = 0.8325, p < 0.01), showed strong correlations with the RNA-Seq data, and 7 genes showed poor correlation with the corresponding protein expression (Fig. 4; Table 3).
Table 3
Correlation analysis between cell wall related gene and protein expression levels
Gene/Protein name | Gene ID | Pearson correlation efficient | P-value | Numbers |
Gene | | | | |
DIR2 | TRINITY_DN135342_c3_g3 | 0.7977 | 0.0100 | 18 |
HCT | TRINITY_DN138969_c0_g6 | 0.2848 | 0.4577 | 18 |
CSLG | TRINITY_DN136551_c1_g11 | -0.06589 | 0.8663 | 18 |
NAClike | TRINITY_DN136342_c2_g5 | 0.9464 | 0.0001 | 18 |
EXP1 | TRINITY_DN141308_c1_g4 | 0.8582 | 0.0031 | 18 |
CAD | TRINITY_DN141875_c0_g4 | 0.8015 | 0.0094 | 18 |
β-GAL1 | TRINITY_DN138388_c1_g1 | 0.8440 | 0.0169 | 18 |
β-GAL2 | TRINITY_DN142386_c5_g1 | 0.6675 | 0.0495 | 18 |
PG | TRINITY_DN196976_c0_g1 | 0.7863 | 0.0636 | 18 |
4CL | TRINITY_DN141686_c0_g5 | 0.7466 | 0.0208 | 18 |
ENDOB | TRINITY_DN138197_c1_g5 | 0.7052 | 0.0338 | 18 |
NAC | TRINITY_DN131789_c1_g4 | 0.1755 | 0.6515 | 18 |
PE | TRINITY_DN143028_c0_g1 | 0.8042 | 0.0292 | 18 |
PRX | TRINITY_DN142336_c1_g1 | 0.4798 | 0.1912 | 18 |
BXL | TRINITY_DN141432_c1_g2 | 0.1761 | 0.6505 | 18 |
PL | TRINITY_DN143250_c1_g6 | 0.6021 | 0.0862 | 18 |
BHLH | TRINITY_DN138043_c3_g7 | 0.7485 | 0.0327 | 18 |
PG3 | TRINITY_DN142042_c0_g2 | 0.6819 | 0.0430 | 18 |
CEL | TRINITY_DN76417_c0_g1 | 0.8732 | 0.0103 | 18 |
PRX2 | TRINITY_DN141264_c1_g1 | 0.8325 | 0.0054 | 18 |
Protein | | | | |
DIR1 | TRINITY_DN135837_c0_g4 | 0.8316 | 0.0402 | 18 |
PG2 | TRINITY_DN142943_c1_g1 | 0.8336 | 0.0101 | 18 |
EXP1 | TRINITY_DN141308_c1_g4 | -0.7482 | 0.0204 | 18 |
PRX3 | TRINITY_DN139379_c0_g3 | -0.1040 | 0.7900 | 18 |
F26G | TRINITY_DN142424_c1_g1 | 0.8907 | 0.0013 | 18 |
BGLU33 | TRINITY_DN137437_c3_g1 | 0.7248 | 0.0272 | 18 |
PE | TRINITY_DN143028_c0_g1 | 0.7596 | 0.0176 | 18 |
PRX | TRINITY_DN142336_c1_g1 | 0.7489 | 0.0202 | 18 |
BXL | TRINITY_DN141432_c1_g2 | 0.6766 | 0.0453 | 18 |
α-HY | TRINITY_DN141662_c1_g4 | 0.8270 | 0.0060 | 18 |
PL | TRINITY_DN143250_c1_g6 | 0.4460 | 0.2289 | 18 |
PRX2 | TRINITY_DN141264_c1_g1 | 0.5422 | 0.1315 | 18 |
PRX4 | TRINITY_DN137008_c1_g6 | 0.7481 | 0.0204 | 18 |
PRX5 | TRINITY_DN139660_c0_g1 | 0.6748 | 0.0462 | 18 |
ENBG | TRINITY_DN141880_c0_g1 | 0.8917 | 0.0029 | 18 |
β-GAL1 | TRINITY_DN138388_c1_g1 | -0.8215 | 0.0234 | 18 |
β-GAL2 | TRINITY_DN142386_c5_g1 | 0.7797 | 0.0132 | 18 |
Quantitative proteome analysis
To understand the molecular mechanisms of peel cracking in A. trifoliata fruits, a quantitative proteomics analysis was also performed using the TMT platform and LC-MS/MS analysis during fruit development, to complement the transcriptome analysis. Accordingly, a total of 812 625 spectra, 68 151 identified spectra, 12 456 peptides, and 10 572 unique peptides, were found by proteomic analysis and 2839 proteins were identified (Table 1). In terms of protein mass distribution, proteins with molecular weight greater than 9 kDa have a wide range and good coverage, with the maximum distribution area of 10–40 kDa (Fig. S3a). The peptide quantitative analysis of the proteins showed that the protein quantity decreased with the increase of the matching peptide (Fig. S3b).
Among of these proteins, 240 were identified as differentially abundant proteins (DAPs) using a fold-change > 1.2 and < 0.83 with p < 0.05 as the up-regulated and down-regulated threshold, respectively. In the comparison between PM and PS, 84 proteins were more abundant and 106 proteins were less abundant in the PM than PS and are shown in a volcano plot in Fig. S2c; 20 DAPs were more abundant and 13 DAPs were less abundant in the PL than PM, and are shown in a volcano plot in Fig. S2d, and 17 were co-expressed at PM_PS and PL_PM.
Functional classification of the identified DAPs
Bioinformatics analysis of DAPs was carried out based on protein functional classifications and hierarchical cluster analysis. GO analysis showed that most of the DAPs in the biological processes (BP) were involved in cellular responses to the chemical stimulus and cellular oxidant detoxification processes in the PM_PS, and the metabolic and macromolecule metabolic processes in the PL_PM, respectively. The oxidoreductase and antioxidant activity in PM_PS and the cytoplasmic part and protein-containing complex in PL_PM, were the highest portions of DAPs in the category of molecular functions (MF). Extracellular regions in the PM_PS, and cytoplasmic parts of the PL_PM, were the highest portions of DAPs in the category of cell components (CC) (Fig. 2c-2d).
Moreover, a KEGG pathway analysis was carried out to further evaluate the DAPs. In the comparison between the PM_PS and PL_PM, many DAPs were enriched in two-component systems and ribosome pathways. Notably, most of the DAPs involved in cell wall-related DAPs were upregulated in both the PM_PS and PL_PM groups. While those DAPs involved in phenylpropanoid pathways and peroxisomes were downregulated 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 wall-related DAPs. A total of 40 cell wall-related DAPs were clustered closely, in both the PM_PS and PL_PM groups. Most of these were involved in 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 (Fig. 3g-l).
Confirmation of fruit cell-wall related proteins in A. trifoliata by qPCR
To validate the results of the TMT data, an expression and correlation analysis between the qPCR and the FPKM values was obtained when the TMT was performed. There were 17 genes selected that had shown differential expression patterns in the PM_PS and PL_PM groups, and the results of the qPCR of these genes are shown in Fig. 6. Specifically, in the PM_PS group, the phenylpropanoid pathway related genes (PRX, PRX3, PRX4 and PRX5), galactose metabolism related genes (β-GAL1 and β-GAL2), amino sugar and nucleotide sugar metabolism related gene (BXL), and pentose and glucuronate interconversion related genes (PG and PG2), were significantly down-regulated. While starch and sucrose metabolism related genes- glucan endo-1,3-beta-glucosidase (ENBG), Furostanol glycoside 26-O-beta-glucosidase (F26G), pentose and glucuronate interconversion related genes (PL and PE), and cell wall metabolism related genes (EXP1, DIR1) were significantly up-regulated. In the PL_PM group, most of the genes were up-regulated, except for DIR1, EXP1, BGLU33, PG, PG2, β-GAL1, and alpha/beta hydrolase (α-HY). Moreover, the expression of the 14 candidate genes, including DIR1 (r = 0.8316, p < 0.05), PG2 (r = 0.8336, p < 0.05), EXP1 (r = -0.7482, p < 0.05), F26G (r = 0.8907, p < 0.01), BGLU33 (r = 0.7248, p < 0.05), PE (r = 0.7596, p < 0.05), PRX (r = 0.7489, p < 0.05), BXL (r = 0.6766, p < 0.05), α-HY (r = 0.8270, p < 0.01), PRX5 (r = 0.7481, p < 0.05), PRX6 (r = 0.6748, p < 0.05), β-GAL2 (r = 0.7797, p < 0.05), ENBG (r = 0.8917, p < 0.01) and β-GAL1 (r = -0.8215, p < 0.05), showed strong correlations with the RNA-Seq data, 3 genes showed poor correlation with their corresponding protein expression (Fig. 5; Table 3).
Comparative analysis between protein abundance and gene expression levels
To evaluate the relationships between the transcriptomic and proteomic changes during fruit pericarp cracking, the quantitative data for DEGs and DAPs was used for correlation analysis. According to the association analysis, 1904 shared DEGs and 17 shared DAPs were identified in the comparison between the PM_PS and PL_PM. Among the shared DEGs, 1123 were up-regulated and 781 were down-regulated in the PM_PS, whereas 808 were up-regulated and 1096 were down-regulated in PL_PM. Of the common DAPs, 9 were more abundant and 8 were less abundant in PM_PS than in PL_PM, while 8 were more abundant and 9 were less abundant in PL_PM than in the PM_PS. (Table 3). Furthermore, 14 and 4 DAPs and their corresponding DEGs were identified in the PM_PS and PL_PM groups, respectively. Of these, 12 DAPs (4 with increased abundance and 8 with decreased abundance) and 4 DAPs (2 with increased abundance and 2 with decreased abundance) were regulated in the same direction as their corresponding DEGs in the PM_PS and PL_PM groups, respectively (Fig. 6a-6b). There were more DEGs than DAPs in both the PM_PS and PL_PM groups, with significant differences of the trends in transcript levels and protein abundance.
Furthermore, the fold-changes of the DAPs were significantly negatively correlated with their corresponding DEGs by Pearson’s correlation tests (r = 0.03 and 0.11, p < 0.01, in PM_PS and PL_PM, respectively), indicating a poor correlation between transcript levels and protein abundance (Fig. 6c-6d).
Identification 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 analysis. The largest groups within the biological processes category were those linked with metabolic and cellular processes; catalytic activity and binding mainly included membranes, cell and cell parts were predominant in the category of molecular functions, both in PM_PS and PL_PM (Fig. S4a-b).
Moreover, a pathway enrichment analysis was conducted in the PM_PS and PL_PM groups, based on the KEGG database. In the PM_PS, 14 DAPs were significantly enriched in 7 pathways, both in DEGs and DAPs, which including 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, 4 DAPs were significantly enriched in the 3 pathways, both in the 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 2 pathways, including pentose and glucuronate interconversions (2 DEGs and DAPs), and galactose metabolism pathways (2 DEGs and DAPs) were shared in the transcriptome and proteome data, for both the PM_PS and PL_PM groups. While the 2 DEGs and DAPs involved in the phenylpropanoid biosynthesis pathways were only shared by DAPs and DEGs in the PM_PS group. Therefore, involved in the shared DAPs and DEGs of the two pathways were further investigated as candidate genes related to the A. trifoliata fruit peel cracking.
Analysis of proteins expressed in A. trifoliata fruits identifies genes that might play relevant roles in fruit peel cracking
Of the two candidate pathways, a total of 13 and 3 DAPs and 28 and 46 DEGs were detected in the PM_PS and PL_PM groups, respectively, where indicated that more DEGs than DAPs were involved in the peel cracking of A. trifoliata fruit (Fig. 3a-3H). Moreover, most of these DEGs and DAPs were down-regulated in the PM_PS group, and up-regulated in the PL_PM group, both in the transcriptome and proteome data. Notably, PL was up-regulated in the PM_PS group, and the PE and β-GAL2 genes were up-regulated in the PL_PM both in the transcriptome and proteome data and showed strong positive correlations. While β-GAL1 had a negative correlation between the transcriptome and proteome data, as it was down-regulated in the transcriptome and up-regulated in the proteome. Furthermore, the results indicated that most of the genes encoding DAPs were not included in the DEGs, which was accorded with the observed differences between the proteome and transcriptome data.
Additionally, the signal transmissions underlying the A. trifoliata fruit peel cracking were studied by analyzing the protein–protein interactions using the STRING database. In the interaction network of these DEGs and DAPs (Fig. 7), there were two pathways identified, PE interacted with PL and β-fructofuranosidase was interacted with raffinose synyhase. These results indicated that PE and PL should be further investigated as candidate proteins in the development of A. trifoliata fruit peel cracking.