Normalization and clustering of RNA-Seq-based transcriptomes
A comparative statistical analysis of gene expression was performed with the edgeR program (v3.36.0) (Robinson et al. 2009). A total of 16386 differentially expressed genes (DEGs) were detected (Supplementary Table 1/ All_EdgeR.LogData). Of this total, 5083 were down-regulated (log2 <= -1) and 3872 were up-regulated (log2 > = 1) in the protoplasts. The remaining 7431 genes did not fluctuate.
Figure 1.1 shows a heat map of the DEGs of leaves, leaves 3 hours after the protoplast isolation procedure, and the protoplasts. The heat map analysis shows that the expression data of DEGs were almost identical to each duplicate sample. The DEG expression data were analyzed by principal component analysis (PCA). The PCA (Fig. 1.2) shows that the variability of gene expression data in leaves increased three hours after the protoplast isolation procedure (Leaf3h).
MultiExperiment Viewer (MeV, v4.8.1) is based on the k-means method and produced the hierarchical clustering analysis of the DEGs. The clustering analysis of co-regulated gene expression patterns was performed with MeV. The 16386 DEGs were clustered by MeV (v4.8.1) based on the k-means method (k = 9). Clustering of the DEG expression data by k-means cluster analysis (k = 9) showed a temporal change of the gene expression from leaves to protoplasts (Supplementary Fig. 1). The data for the six sampling points were clustered into nine groups. The cluster number and corresponding expression values for each DEG in the six sampling point clusters are shown in Supplementary Table 2.
Comparison between two different time points during the process of protoplast isolation
The number of genes with a CPM greater than 0.5 between the two samples was 20,286 (Supplementary Table 1/Count.all.AT1). Between the leaves three hours after the onset of the protoplast isolation procedure (Leaf3h) and the untreated leaves as controls (Leaf0h), 8223 DEGs with an FDR less than 0.01 were detected from 20286 genes (Supplementary Table 1/Leaf0h_Leaf3h_significant); 3741 genes were up-regulated and 4482 genes were down-regulated (Supplementary Table 1/Leaf0h_Leaf3h_induction/Leaf0h_Leaf3h_repression). A total of 4537 DEGs were found (Supplementary Table 1/Leaf3h_Protoplast_significant) between Protoplasts and Leaf3h; 2132 were up-regulated and 2405 were down-regulated (Supplementary Table 1/Leaf3h_Protoplast_induction/Leaf3h_Protoplast_repression). Moreover, 8734 DEGs were found between Protoplasts and controls (Leaf0h) (Supplementary Table 1/ Leaf0h__Protoplast_significant ); 3846 were up-regulated and 4888 were down-regulated (Supplementary Table 1/Leaf0h__Protoplast_induction/ Leaf0h__Protoplast_repression). These results indicated that there were more down-regulated DEGs than up-regulated ones during the protoplast isolation.
Figures 2.1a and b show Venn diagrams of the up- and down-regulated genes between the two treatments. The number of down-regulated genes was higher than the number of up-regulated genes in all treatments; 391 were up-regulated and 1164 were down-regulated. The most up-regulated gene in the protoplasts was the serine protease inhibitor gene (AT5G43570: 11792.48-fold), followed by a peroxidase (AT2G38390: 9538.36-fold). The expression levels of the two genes in Leaf3h decreased 63.257-fold and 122.31-fold. A gene with significant down-regulation was an MYB transcription factor (AT5G07700, MYB76) that was down-regulated 0.0002595-fold in protoplasts and 0.009194-fold in Leaf3h.
GO over-representation analysis
The up- or down-regulated DEGs overlapping in the sample treatment of Leaf0h, Leaf3h, and Protoplast were analyzed through GO over-representation analysis by using a web application, PANTHER (http://geneontology.org) (Fig. 2.2). Supplementary Table 3 includes the data. The protoplast isolation procedures activated genes associated with catalytic activity (GO:0003824, molecular function category) and metabolic process (GO:0008152, biological process category), while the procedures also inactivated genes involved in metabolic process (GO:0008152, biological process category).
A Venn diagram was created from the GO analysis data. The four GOs overlapped in the three treatments, and one of the four GOs was related to ethylene. Transcription factor ATG5G39610 (ANAC092, ATNAC2) was found in the ethylene-related GO. DEGs were assigned to biological process (BP), cellular component (CC), and molecular function (MF) categories based on the enrichGO package of the R program.
As a result of the GO over-representation analysis by enrichGO, various transport systems in the BP category were enhanced in both Leaf3h and Protoplasts compared with Leaf0h (Fig. 3a). STP13 (AT5G26340) encodes a protein with high affinity, i.e., hexose-specific/H + symporter activity. The study of the expression of this gene in PCD mutants demonstrated a tight correlation between this gene's expression and PCD (Norholm et al. 2006). The expression levels of this gene in Leaf3h and Protoplasts exhibited a 160.65-fold and 39.840-fold increase compared to Leaf0h, but decreased to 0.24799-fold in Protoplasts compared to Leaf3h. The protoplast isolation procedure led to enhanced expression of proteasome-related genes, as there were many DEGs involved in the proteasome assembly in the BP category in Protoplasts compared with Leaf0h and Leaf3h (Fig. 3a).
The gene expression associated with oxidative stress and fungus GO increased significantly (Fig. 3a). GO terms for endoplasmic reticulum (CC category) were enhanced in both Leaf3h and Protoplasts. Mitochondrial-related GO terms were also enhanced in Protoplasts (Fig. 3b). The results of the GO over-representation analysis for the MF category are presented in Supplementary Fig. 2. Gene groups associated with chloroplast and photosynthesis in BP and CC categories were repressed in both Leaf3h and Protoplasts compared with Leaf0h (Figs. 3c and d). When DEGs in Protoplasts were compared with DEGs in Leaf3h, GO terms of various transport systems to the nucleus were inhibited (Fig. 3c). The genes encoding protein kinases in the MF category were repressed in Protoplasts compared with Leaf0h (Supplementary Fig. 2).
DEGs Related to Senescence
Plant cell death is closely related to cellular senescence and photosynthesis. A total of 182 genes were listed by PODC (http://bioinf.mind.meiji.ac.jp/podc/) as senescence-associated genes (Supplementary Table 4). A total 64 genes were up-regulated and 40 were down-regulated in Leaf3h, whereas 68 and 38 genes were up- and down-regulated in protoplasts, respectively. When the gene expression of protoplasts was compared with Leaf3h, 52 genes were up- and down-regulated. These results suggest that the expression of senescence-associated genes was up-regulated relatively soon after the protoplast isolation treatment.
Temporal distribution of primary metabolites during protoplast isolation
Water-soluble primary metabolites measured by GC–MS (Supplementary Table 5/GC–MS normalized data) were analyzed by principal component analysis (Fig. 4.1). The analysis shows that the primary component separated metabolites of protoplasts from other leaf metabolites. Figure 4.2 shows a heat map and clustering analysis of log2 fold ratios between the two samples. The content of most metabolites decreased during the first 6 h of protoplast isolation and slightly increased from the intermediate point to the end of protoplast isolation. However, the metabolite content in protoplasts hardly recovered, except for proline.
Figure 4.3 shows correlation matrices of water-soluble primary metabolites extracted from leaves 3 and 6 h after the onset of protoplast isolation and from protoplasts. The metabolites with a high correlation coefficient are shown in red in the panels. There was a drastic metabolical conversion from 3 to 6 h, which then partially returned to the control (Leaf0h) level at the end of the protoplast isolation.
Sugar and organic acid contents were measured via LC–MS/MS to obtain more quantitative data (Table 1). The galactose, arabinose, fructose, glucose, and sucrose contents decreased in leaves 3 h after protoplast isolation but increased more than 10-fold in protoplasts. Organic acids measured in the LC–MS/MS experiments showed almost the same result as GC–MS data except for galacturonate, whose content increased in protoplasts.
These figures indicated that sucrose and amino acids increased in protoplasts compared with Leaf3h or Leaf6h. Considering these results and the increase in catabolic activity shown by the GO over-representation analysis (Fig. 2.2), the protoplast isolation procedures induced protein degradation.
Carbon and oxidative phosphorylation metabolic pathways based on KEGG pathway analysis
The expression profiles of all DEGs and the primary metabolites were color-coded in each pathway with expression data from the three time points in the following order: Leaf3h vs. Leaf0h, Protoplasts vs. Leaf3h, and Protoplasts vs. Leaf0h.
Figure 5.1 shows changes in sugar metabolic pathways during the process of protoplast isolation. DEGs involved in photosynthesis were down-regulated as in the GO over-representation analysis, but the process of protoplast isolation hardly affected the DEGs involved in oxidative phosphorylation (Supplementary Fig. 3). Most of the DEGs involved in the TCA cycle were up-regulated, as predicted by the results of the GO over-representation analysis. The DEGs involved in the glycolysis and pentose phosphate pathways were up- or down-regulated during the process of protoplast isolation. The sucrose and glucose ratio of Leaf3h and Leaf0h decreased to 0.33 and 0.00, respectively, but increased to 13.74 and 24.28 at the end of the protoplast isolation.
The content of organic acids for TCA cycle components decreased during the first 3 h of the protoplast isolation procedure and slightly increased until the end of the protoplast isolation, but the content in protoplasts remained at a lower level than in the leaves.
DEGs involved in protein degradation machinery, plant–pathogen interaction, and plant hormone signal transduction
The GO over-representation analysis and primary metabolite analysis by GC–MS revealed that catabolic activity was enhanced throughout the protoplast isolation (Fig. 2.2, Fig. 4.2). Therefore, the KEGG pathway analyzed DEGs related to autophagy and the proteasome for protein degradation (Fig. 5.2). These results showed that most of the DEGs related to autophagy and the proteasome were significantly up-regulated, suggesting that a stimulus generated by the process of protoplast isolation triggered proteolysis. AT4G01610 encodes a cathepsin B-like protease (ATCATHB3). The cathepsin gene expression increased throughout the protoplast isolation, i.e., 1.04-fold in Leaf3h and 2.30-fold in Protoplasts compared to Leaf0h as a control (Supplementary Table 1/All_EdgeR.LogData).
Since the GO over-representation analysis also revealed that the immune system process in the biological function category was up-regulated (Fig. 2.2), DEGs were analyzed by the KEGG pathway of plant–pathogen interaction (Fig. 5.3). The results shown in Fig. 2.2 and Fig. 5.3 indicate that elicitors initiated a cascade of defense response events and senescence-related syndromes, such as activation of catabolism and proteolysis.
Gene expressions of the chitin elicitor receptor, CERK1, bacterial flagellar protein (flg22) receptor, FLS2, bacterial translational elongation factor (elf18) receptor, EFR, and another pattern-recognition receptor, SERK4, were up-regulated during the first 3 h of protoplast isolation and decreased or remained constant thereafter. These elicitor receptor-mediated signal transductions stimulated the immunity triggered by pathogen-associated molecular patterns (PAMPs). Interestingly, hypersensitive response genes were down-regulated or remained constant except for the CR88 gene. Up-regulation of WRKY22 and WRKY29 accounted for defense-related gene induction.
The KEGG pathway analysis assigned DEGs for plant hormone signal transduction (Fig. 5.4). These results indicated ethylene, jasmonic acid, and salicylic acid-activated ESR2, JAZ5, and TGA9 gene expression during the first 3 h of protoplast isolation and at the end of protoplast isolation. These genes were the initial steps for endoplasmic reticulum stress responses, proteolysis, and disease resistance responses. The KEGG pathway analysis suggested that some unknown elicitors and plant growth substances promoted senescence-related syndromes and pathogen defense immune systems.
Transcriptome profiles of transcription factors
Transcription factors play essential roles in many biological aspects of life cycles. Supplementary Fig. 4 shows a heat map for the temporal distribution profile of differentially expressed transcription factor genes during the protoplast isolation. There were more down-regulated transcription factors than up-regulated ones. These results indicated that the protoplast isolation procedure triggered transcriptional reprogramming during the protoplast isolation.
The Plant TFDB (plant transcription factor database v4.0) lists 1726 transcription factor genes (Supplementary Table 6). The number of DEGs on the list was 1079. Leaf3h vs Leaf0h had genes with significantly up-regulated (p < 0.05) ERF, HFS, Myb, Trihelix, and WRKY families. Many other transcription genes showed no significant difference between Protoplasts and Leaf0h and suggested relatively early responses. On the other hand, the late-response TF families that significantly increased in Protoplast vs Leaf0h were ARF, bHLH, C2H2, CO-like, ERF, GRAS, Myb-related, MFY-A, TCP, YABYY, and the ZFHD family.
Two TF families (HSF and Trihelix) contained significantly up-regulated genes in Leaf3h, but not in the comparison between Protoplasts and Leaf0h. On the other hand, several TF families (C2H2, CO-like, GATA, GRAS, SBP, TCP, YABBY, and ZF-HD), containing genes whose expression was suppressed, were found in Leaf3h but not in the comparison between Protoplasts and controls. The TF families containing genes that were significantly suppressed in Leaf3h were b-HLH, C2H2, CO-like, G2-like, GATA, GRAS, HD-ZIP, SBP, TCP, YABBY, and the ZFHD family. Among these, there were many TF families that were highly expressed in Protoplasts. In other words, the gene expression of some TF families was once suppressed by Leaf3h but then restored in Protoplasts, but the expression was higher than that under Leaf0h.
A total of 391 up-regulated DEGs in the Venn diagram overlapped among Leaf0h, Leaf3h, and Protoplasts (Fig. 2.1). There were 26 transcription factors in these DEGs (Supplementary Table 1); five WRKY, six MYB, four MAC, and four ERF transcription factors accounted for more than 70% of the 26 up-regulated transcription factor genes that overlapped in the Venn diagram of Leaf0h, Leaf3h, and Protoplasts. Figures 6a, b, and c show the co-expression of these up-regulated transcription factor genes.
AT5G39610 (ANAC092, ATNAC2) encodes an NAC-domain transcription factor and positively regulates senescence-induced cell death and senescence in leaves. ANAC092 gene expression was significantly activated throughout the entire period of protoplast isolation. There was a 277.97-fold increase in the leaves three hours after protoplast isolation (Leaf3h) and a 2442.40-fold increase in Protoplasts, compared with the control leaves (Leaf0h). Therefore, to explore the regulatory network administered by ANAC092, co-expression of ANAC092 and other up-regulated genes extracted from the overlapping Venn diagram in Leaf0h, Leaf3h, and Protoplasts was analyzed by Cytoscape (cytoscape.org), which is an open-source software platform for visualizing complex networks (Figs. 6d, e and f). These figures indicate that the NAC and WRKY transcription factors were closely co-expressed during the process of protoplast isolation. Supplementary Fig. 5 shows the co-expression of STP13 genes and transcription factor genes found as DEGs. ANAC092 was co-expressed with STP13 and was significantly up-regulated throughout the protoplast isolation.
AT3G11020 (DREB2B) encodes a member of the DREB subfamily in A-2 of the ERF/AP2 transcription factor family, which includes four major subfamilies: the AP2, RAV, ERF, and dehydration-responsive element-binding protein (DREB) subfamilies (Sakuma et al. 2002). Plant abiotic stress-induced DREB subfamily members regulate gene expression via the cis-acting dehydration-responsive element/C-repeat (DRE/CRT) element. ERF subfamily members that bind to ethylene response elements (EREs) respond to abiotic stress. Co-expression analysis indicated that 120, 122, and 96 transcription factor genes that were either up- or down-regulated during the protoplast isolation were identified as the first neighbors of DREB2B (Supplementary Fig. 6).
The three overlapped treatments resulted in 1163 down-regulated DEGs (Fig. 2.1). There were 48 transcription factors, which amounted to about 4% in these overlapped treatments. In contrast to up-regulated transcription factors, only one ethylene-responsive factor (ERF), two MYB, and four bHLH transcription factors were found in the three overlapped treatments. Up-regulated DEGs in the three overlapped treatments did not indicate bHLH transcription factors.
The AT1G50030-encoded TOR (target of rapamycin) protein belongs to the family of phosphatidylinositol 3-kinases and is the target of the antiproliferative drug rapamycin. Figure 9 shows the co-expression of TOR and other first neighbor DEGs analyzed by Cytoscape and indicates that the expression of TOR decreased during the protoplast isolation. There was no transcription factors co-expressed with TOR. Almost all of these co-expressed genes were down-regulated. However, AT5G54620 encoding an ankyrin repeat-containing repeat protein was up-regulated during the early stage of protoplast isolation.