Catabolic reprogramming of Brassica rapa leaf mesophyll protoplasts during the isolation procedure

The use of leaf mesophyll protoplasts for transformation and genome editing in plants is expected due to the distinctive features of protoplasts, such as cell-wall-free single cells. However, the application of leaf mesophyll protoplasts for molecular breeding is limited because of the recalcitrance of protoplasts to regenerate. We speculated that the primary reason for this recalcitrance is senescence during protoplast isolation before initialization. In this study, we performed profiling, clustering, co-expression analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of primary metabolites and transcriptomes of Brassica rapa leaves and leaf mesophyll protoplasts to reveal the reason for senescence. Primary metabolite profiling indicated that the metabolism of B. rapa protoplasts was converted to catabolism. The number of downregulated differentially expressed genes (DEGs) increased during protoplast isolation. An analysis of GO overexpression revealed the activation of genes involved in catabolic and immune system processes and the inactivation of chloroplasts and photosynthetic genes. KEGG pathway analysis showed that the activation of autophagy and proteasomes accounted for senescence-associated proteolysis during protoplast isolation. It also revealed activation of the genes involved in endoplasmic reticulum stress responses and disease resistance responses during the first 3 h of isolation. These results suggest that elicitor receptor-mediated signal transduction stimulates the pathogen-associated molecular patterns (PAMPs), which triggers immunity. There were more downregulated than upregulated transcription factors during protoplast isolation. AT5G39610 (ANAC092, ATNAC2) gene expression was significantly activated throughout the entire period of protoplast isolation. ANAC092 was co-expressed with upregulated AT5G26340 (STP13), which encodes a protein with high affinity, hexose-specific/H + symporter. AT1G50030 encodes target of rapamycin (TOR) proteins, and the expression of TOR decreased during protoplast isolation. The presence of ethylene and the inhibition of photosynthesis-related genes by glucose and sucrose promote senescence-associated gene expression of protoplasts. Since a decrease in glucose downregulates glucose TOR signaling, the inactivated TOR signaling will promote catabolic reprogramming, senescence, and autophagy in protoplasts during protoplast isolation. We concluded that the initialization of protoplasts requires the blocking of these complicated crosstalk signaling pathways to prevent the catabolic reprogramming of protoplasts, which will lead the way to new plant breeding techniques.


Introduction
Plant protoplasts can provide solutions for improving plants by either DNA transfer or genome editing (Yue et al. 2021). The removal of plant cell walls allows for the introduction of foreign DNA, RNA, or proteins into protoplasts through either polyethylene glycol (PEG) treatment or electroporation (Reed and Bargmann 2021). Due to the single cell, protoplasts can potentially avoid chimerism in transformation. These characteristics represent advantages over leaf-or callus-based plant transformation and pave the way to new plant breeding techniques.
Since Cocking developed the method of protoplast isolation by enzymatic digestion, protoplasts from many higher plant species and organs have been isolated (Cocking 1960;Roest and Gillissen 1989;Yue et al. 2021). However, most procedures in protoplast isolation and regeneration have been developed using an empirical approach (Davey et al. 2005). Previous studies attempted to optimize the physiological state of donor plants through foliar absorption of plant growth regulators, nutrients, and cycloheximide (Cassells and Tamma 1987;Kaur-Sawhney et al. 1976); preconditioned plants under short periods of low-intensity illumination (Shepard and Totten 1977); excised leaves placed in the dark (Johnson et al. 1981); and plants grown with calcium nitrate and supplementary light (Cassells and Cocker 1982). It is believed that protoplasts derived from developed organs retain the properties of the original organs (Yue et al. 2021). However, protoplasts are not simply cell-wall-free replicates of original cells, but are subjected to abiotic stress and complex metabolic modifications during their isolation procedure (De Marco and Roubelakis-Angelakis 1996). Several factors affect the isolation, culture, and regeneration of plants from protoplasts (Roest and Gillissen 1989;Davey et al. 2005). Kartha et al. (1974) isolated and cultured Brassica napus leaf mesophyll protoplasts. Since then, many studies have developed reliable and efficient protoplast culture procedures for plant regeneration of different Brassica species. B. rapa contains many cultured varieties within the genus. Leaf mesophyll protoplasts of B. rapa respond very differently when cultured. B. rapa protoplasts appear to be more recalcitrant in culture than those of other Brassica species and only rarely regenerate from protoplasts into plants. More details can be found in a review by Uddin (2013).
There is little understanding of the nutritional and physiological states of leaf mesophyll protoplasts. An increase in RNase activity was found to be accompanied by morphological instability and deterioration of oat leaf protoplasts (Altman et al. 1977). Leaf excision and wounding are prerequisites for leaf protoplast isolation and enhance the activity of senescence-inducing enzymes such as proteases, nucleases, and other hydrolases (Udvardy et al. 1969;Wyen et al. 1971). An osmotic stress-induced increase in putrescine and a blockage of its conversion to higher polyamines partly account for the failure of cereal mesophyll protoplasts to divide consistently (Tiburcio et al. 1986). When B. napus leaves are wounded before the enzymatic digestion-induced senescence of protoplasts prior to culture, the free endogenous putrescine, ethylene, and NADH-GDH activity of protoplasts increased (Watanabe et al. 1994). A previous study also showed that leaf protoplasts senesced during culture, which may account for the recalcitrance of B. napus leaf protoplasts (Watanabe et al. 1998). Propidium iodide (PI) and DAPI staining, transmission electron microscopy, and flow cytometry analysis detected the cleavage of nuclear DNA into oligonucleosomal fragments during the culture of B. napus leaf protoplasts. Further, TUNEL detected the cleavage of nuclear DNA into oligonucleosomal fragments on a southern hybridization membrane; nuclear DNA fragments labeled with DIG-dUTP at the free 3′-OH ends by terminal deoxynucleotide transferase (TdT) were used as a probe for Southern hybridization (Watanabe et al. 2002a). These results suggest that senescence of B. napus with culture time is accompanied by apoptosis-like cell death.
However, Petunia hybrida leaf protoplasts increased their spermidine and spermine contents with little increase in ethylene emission during culture (Watanabe et al. 2002b). Petunia leaf protoplasts can initiate cell division and form calluses. Their plating efficiency was 0.7% (Watanabe et al. 1992). The reasons why recalcitrant leaf protoplasts are unable to divide remain elusive. Leaf senescence and the associated cell death are developmentally programmed processes. Considering the previous studies, a key to overcoming the recalcitrance of protoplast initialization is to elucidate the underlying molecular mechanism of protoplast isolationinduced senescence. In the present study, we investigated changes in primary metabolites and transcriptome expression during protoplast isolation.
These data indicate that protoplast metabolism was converted to catabolism during protoplast isolation. This catabolic reprogramming agrees with the results of previous studies, in which protoplasts were already senesced at the time of isolation. Transcriptome expression analysis revealed the upregulation of genes involved in the pathogen defense immune system process and TOR signaling, which account for catabolic reprogramming and senescence during protoplast isolation.

Plant materials
Seeds of B. rapa cv. Himehassen-chingensai were purchased from Nakahara Seed Product Co. Ltd. (Hukuoka, Japan). Plants were grown under a 14 h photoperiod at a photon fluence density of 50 µmol m −2 s −1 beneath LED tubes and under a day/night temperature cycle of 24/20 °C. The upper developing leaves of plants were harvested 2-3 months after germination and used for water-soluble metabolite and RNA extraction.
Water-soluble metabolites and mRNA were extracted from B. rapa leaves as control and leaf mesophyll protoplasts, then 3 h after treatment of cell wall digestion enzymes, water-soluble metabolites and mRNA were extracted; 6 h after treatment, water-soluble metabolites were extracted.

Preparation of RNA sequencing samples
The B. rapa leaves were immersed in 0.5% sodium hypochlorite and rinsed several times with sterilized distilled water. The sterilized B. rapa leaves were cut into narrow strips and incubated in a filter-sterilized enzyme solution containing 1.0% Cellulase Y-C (Kyowa Kasei, Osaka, Japan), 0.05% Pectolyase Y-23 (Kyowa Kasei, Osaka, Japan), 5 mM MES (pH 5.8), and 0.6 M sorbitol. Then, 3 h after incubation, excised leaves were filtered through nylon mesh and rinsed several times with sterilized distilled water. The leaves were placed on filter paper and dried. Then, 14 h after incubation, protoplasts were isolated as previously described (Watanabe et al. 2002a). The dried leaves and isolated protoplasts (5 × 10 5 ) were placed in Eppendorf tubes, frozen by liquid nitrogen, and stored in a freezer at − 80 ℃ until RNA extraction.

RNA extraction and sequencing
Two independent biological replicates were performed for the transcriptome analysis. Total RNA extraction from leaves and protoplasts was performed using the NucleoSpin RNA Plant (TAKARA BIO INC, Japan) according to the manufacturer's instructions. The quantity and quality of extracted RNA were assayed by a 2100 Bioanalyzer and RNA6000 Nano kit (Agilent Technologies Japan, Ltd). For this assay, 1 µg total RNA per sample was used as input material to prepare RNA sequencing libraries. Sequencing libraries were generated using a TruSeq Stranded mRNA library kit (Illumina K.K.) following the instruction manual. Library quality was assessed with a 2100 Bioanalyzer DNA1000 kit and KAPA Library Quantification Kit (Roche). The libraries were subjected to 50-base single read sequencing using the Illumina HiSeq 2500 System, and 22-26 million reads were obtained for each sample. Clean reads were generated by removing reads containing adapter sequences and/or lowquality sequences from the raw data using CLC Genomics Workbench v10.0.0 (QIAGEN). STAR 2.5.3a (https:// github. com/ alexd obin/ STAR) was used to align the reads against the B. rapa reference genome downloaded from the Ensembl genome database (ftp:// ftp. ensem blgen omes. org/ pub/ plants/ relea se-32/ fasta/ brass ica_ rapa/ dna/ Brass ica_ rapa. IVFCA ASv1. dna. tople vel. fa). The transcriptome abundance count data were generated by htseq-count (Anders et al. 2014), and the TAIR ID of the annotation file was added (Brassica_rapa.IVFCAASv1.34. gff3). Because many B. rapa genes have the same TAIR ID, duplicated B. rapa counts were summed as the TAIR ID counts shown in Supplementary Table 1, Count.all.AT1. The read count data were normalized, and differentially expressed genes (DEGs) were identified with edgeR. Genes with a CPM value higher than 0.5 in at least two libraries are listed in Supplementary Table 1, All_EdgeR.LogData.

GO, KEGG pathway, and co-expression analyses
The genes whose expression ratio increased more than twofold or decreased to less than half with FDR less than 0.01 were analyzed. GO and KEGG pathway analyses were performed by the cluster Profiler R package (Yu et al. 2012) using the above data. Cytoscape 3.9.0 software (cytoscape. org) with the GeneMANIA application was used to visualize co-expression networks.

Primary metabolite extraction and GC-TOF-MS analysis
Leaves and protoplasts (n = 6, biological replicates) were immediately frozen in liquid nitrogen to quench enzymatic activity. After lyophilization, each sample was extracted with a concentration of 2.5 mg dry weight of tissue per 1 mL of extraction medium (methanol/chloroform/water, 3:1:1 v/v/v) containing 10 stable isotope reference compounds using a Retsch MM310 mixer mill at a frequency of 15 Hz for 10 min at 4 °C. The extracts were evaporated to dryness in a Savant SPD2010 SpeedVac Concentrator (Thermo Electron Corporation, USA). For methoximation, 30 μL of methoxyamine hydrochloride (20 mg/mL in pyridine) was added to the sample. After 15 h of derivatization at room temperature, the sample was trimethylsilylated for 1 h using 30 μL of MSTFA with 1% TMCS at 37 °C with shaking; 30 μL of n-heptane was added following silylation. All the derivatization steps were performed in a VSC-100 vacuum glove box (Sanplatec, Japan) filled with 99.9995% (G3 grade) dry nitrogen.
Metabolites were analyzed according to previously established methods (Kusano et al. 2007a, b) using a GC-TOF-MS system, which consisted of an Agilent 7890 N gas chromatograph (Agilent Technologies, Wilmington, DE, USA) connected to a Pegasus IV TOF mass spectrometer (LECO, USA). One microliter of each sample was injected in splitless mode by a CTC CombiPAL autosampler (CTC Analytics, Zwingen, Switzerland) into the Agilent 7890 N GC equipped with a fused-silica capillary column with an inner diameter of 30 m × 0.25 mm, with a chemically bound 0.25 μL film Rtx-5 Sil MS stationary phase (RESTEK, USA) for metabolite profiling. Helium gas (G1 grade) was used as the carrier gas at a constant flow rate of 1 mL/min. The temperature program for metabolome analysis started with a 2 min isothermal step at 80 °C, followed by temperature ramping at 30 °C steps to a final temperature of 320 °C, which was maintained for 3.5 min. Data acquisition was performed on a Pegasus IV TOF MS with an acquisition rate of 30 spectra/s in a mass-to-charge ratio range of m/z = 60-800. Alkane standard mixtures (C8-C20 and C21-C40) purchased from Sigma-Aldrich were used to calculate the retention index (RI).
The signal intensity of each metabolite was calculated by a method described previously (Kusano et al. 2007b). The extracted MS spectra were finally identified or annotated according to their RI and by comparison with the reference mass spectra in our libraries in RIKEN CSRS. The data were normalized using the CCMN algorithm (Redestig et al. 2009). All raw data in the NetCDF format were preprocessed by MATLAB (Mathworks, Natick, MA, USA). Data processing was described by Kusano et al. (2020).

Primary metabolite extraction for LC-MS-MS and determination of content
Six biological replicates of leaves and protoplasts were used for extraction and analysis. The metabolites were extracted with 70% ethanol at 60 °C twice and rinsed with cold ethanol. The extract was evaporated by a vacuum evaporator (EYELA CVE-100, Japan) and dissolved in 500 µL ultrapure water filtered through a 0.22 µm membrane. Determination of sugars (arabinose, fructose, galactose, and trehalose) and organic acids (3PGA, galacturonate, lactate, malonate, pyruvate, malate, OAA, and salicylate) was performed by LC-MS-MS (LCMS-8030, Shimadzu, Japan) with a normal phase TSKgel NH 2 -100 column (2.0 mm × 15 cm, 3 µm particles; TOSOH, Japan) at a flow rate of 0.5 mL/min at 40 °C. The sample was analyzed in positive mode. Solvent A was 20 mM ammonium acetate in 95:5 (V/V) water/acetonitrile, and solvent B was acetonitrile. The gradient was t = 0, 85% B; t = 12, 0% B; t = 24, 0% B; t = 26, 85% B; t = 35, 85% B. The injection volume was set to 1 µL.
Mass spectrometry was set to default parameters determined through autotuning. Multiple reaction monitoring (MRM) mode was used for metabolite analysis. The Q3 m/z of each targeted metabolite was determined based on the standard compounds. Supplementary Table 5, Target ions, lists the precursor and Q3 m/z of production.
The contents of ammonium ion, ascorbic acid, ethylene, dehydroascorbic acid, GSH, GSSG, salicylic acid, and spermine are referenced in our previous reports (Watanabe et al. 2011(Watanabe et al. , 2002b. The combined data were used for principal component analysis (PCA) and correlation analysis as log2 transformed data with R packages.

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 16,386 differentially expressed genes (DEGs) were detected (Supplementary Table 1, All_Edge R.Log Data). Among them, 5083 DEGs were downregulated (log2 ≤ − 1) and 3872 were upregulated (log2 ≥ 1) in the protoplasts. The remaining 7431 genes did not fluctuate. Figure 1a shows a heat map of the DEGs of leaves (Leaf0h), leaves 3 h after the protoplast isolation procedure (Leaf3h), and protoplasts. The heat map analysis shows that the expression data of DEGs were almost identical for each duplicate sample. The DEG expression data were analyzed by PCA. PCA (Fig. 1b) shows that the variability of gene expression data in leaves increased 3 h after the protoplast isolation procedure (Leaf3h).
MultiExperiment Viewer (MeV, v4.8.1), which is based on the k-means method, was used to produce the hierarchical clustering analysis of DEGs. Clustering analysis of co-regulated gene expression patterns was performed with MeV. The 16,386 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 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 value for each DEG in the six sampling point clusters are shown in Supplementary Table 2.

Comparison between time points during the process of protoplast isolation
The number of genes with CPM greater than 0.5 between the two samples was 20,286 (Supplementary Table 1, Count.all. AT1). Between the leaves 3 h after the onset of protoplast isolation (Leaf3h) and the untreated control leaves (Leaf0h), 8223 DEGs with an FDR less than 0.01 were detected from 20,286 genes (Supplementary Table 1, Leaf3h_Leaf0h_significant); 3741 genes were upregulated and 4482 genes were  These results indicate that there were more down-regulated than upregulated DEGs during protoplast isolation. Figure 2A (a,b) shows Venn diagrams of the up-and downregulated genes between the two treatments. There were more downregulated than upregulated genes in all a b Fig. 2 Venn diagrams and pie charts of overlapping DEGs of B. rapa leaves, leaves 3 h after protoplast isolation procedure, and leaf mesophyll protoplasts. A Venn diagrams showing overlapping DEGs in B. rapa leaf mesophyll protoplasts during isolation procedure. B Pie charts showing overlapping DEGs categorized by GO overrepresentation analysis: (a) upregulated and (b) down-regulated DEGs treatments: 391 were upregulated and 1164 were downregulated. The most upregulated gene in the protoplasts was the serine protease inhibitor gene (AT5G43570, 11792.48fold), followed by a peroxidase (AT2G38390, 9538.36-fold). The expression level of the two genes in Leaf3h decreased 63.257-fold and 122.31-fold, respectively. A gene with significant downregulation was an MYB transcription factor (AT5G07700, MYB76), which was downregulated 0.0002595-fold in protoplasts and 0.009194-fold in Leaf3h.

GO overrepresentation analysis
The up-and downregulated DEGs overlapping in the sample treatment of Leaf0h, Leaf3h, and protoplasts were analyzed through GO overrepresentation analysis by using the PANTHER web application (http:// geneo ntolo gy. org) ( Figure 2B). Supplementary Table 3 includes the data. The protoplast isolation procedure activated genes associated with catalytic activity (GO:0003824, molecular function category) and metabolic process (GO:0008152, biological process category), while it 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 was related to ethylene. Transcription factor ATG5G39610 (ANAC092, ATNAC2) was found in the ethylene-related GO. DEGs were assigned to the biological process (BP), cellular component (CC), or molecular function (MF) category based on the enrichGO package in R.
As a result of GO overrepresentation 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. A 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 level of this gene in Leaf3h and protoplasts exhibited a 160.65-fold and 39.840-fold increase compared to Leaf0h, but decreased 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 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 terms 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 overrepresentation analysis for the MF category are presented in Supplementary Fig. 2. Gene groups associated with chloroplasts and photosynthesis in BP and CC categories were repressed in both Leaf3h and protoplasts compared with Leaf0h (Fig. 3c, d). When DEGs in protoplasts and Leaf3h were compared, the GO terms of various systems for transport to the nucleus were inhibited (Fig. 3c). 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 are listed by the PODC (http:// bioinf. mind. meiji. ac. jp/ podc/) as senescence-associated genes (Supplementary Table 4). A total of 64 genes were upregulated and 40 were downregulated in Leaf3h, whereas 68 were upregulated and 38 were downregulated in protoplast. When the gene expression of protoplasts was compared with Leaf3h, 52 genes were up-and downregulated. These results suggest that the expression of senescence-associated genes was upregulated relatively soon after 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 (Figure 4a). The analysis shows that the primary component (PC1) separated metabolites of protoplasts from other leaf metabolites. Figure 4b 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 isolation. However, the metabolite content in protoplasts hardly recovered, except for proline. Figure 4c shows correlation matrices of water-soluble primary metabolites extracted from leaves 3 and 6 h after the onset of protoplast isolation and from protoplasts. Metabolites with a high correlation coefficient are shown in red in the panels. There was a drastic metabolic conversion from 3 to 6 h, which then partially returned to the control (Leaf0h) level at the end of protoplast isolation.
Sugar and organic acid contents were measured via LC-MS/MS to obtain more quantitative data (Fig. 6d). The galactose, arabinose, fructose, glucose, and sucrose contents decreased in leaves 3 h after protoplast isolation but increased more than tenfold 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 indicate that sucrose and amino acids increased in protoplasts compared with Leaf3h and Leaf6h. Considering these results and the increased a b catabolic activity shown by GO overrepresentation analysis (Fig. 2B), it can be concluded that the protoplast isolation procedure induced protein degradation. The expression profiles of all DEGs and 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 5A (a,b) and Supplementary Fig. 3a show changes in sugar metabolic pathways during the process of protoplast isolation. DEGs involved in photosynthesis were downregulated [ Fig. 5A (c,d)] as in the GO overrepresentation analysis, but the process of protoplast isolation hardly affected the DEGs involved in oxidative phosphorylation. Most of the DEGs involved in the TCA cycle were upregulated, as predicted by the results of the GO overrepresentation analysis. The DEGs involved in the glycolysis and pentose phosphate pathways were up-or downregulated during the process of protoplast isolation. The sucrose and glucose ratios of Leaf3h and Leaf0h decreased to 0.33 and 0.00, respectively, but increased to 13.74 and 24.28 at the end of isolation.

Carbon and oxidative phosphorylation metabolic pathways based on KEGG pathway analysis
The content of organic acids for TCA cycle components decreased during the first 3 h of protoplast isolation (Leaf3h) and slightly increased until the end of isolation (protoplasts), but the content in protoplasts remained at a lower level than in leaves (Leaf0h).

DEGs involved in protein degradation machinery, plant-pathogen interaction, and plant hormone signal transduction
The GO overrepresentation analysis and primary metabolite analysis by GC-MS revealed that catabolic activity was enhanced throughout the protoplast isolation (Figures 2A  and 4b). Therefore, the KEGG pathway analyzed DEGs related to autophagy and the proteasome for protein degradation ( Figure 5 1 (e,f)). These results showed that most of the DEGs related to autophagy and the proteasome were significantly upregulated, 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 by 1.04-fold in Leaf3h and 2.30-fold in protoplasts compared to Leaf0h as control (Supplementary Table 1, All_EdgeR.LogData).
The KEGG pathway analysis assigned DEGs for plant hormone signal transduction (Fig. 5 2a). 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 isolation. These genes were the initial steps for endoplasmic reticulum stress responses, proteolysis, and disease resistance responses.
Since the GO overrepresentation analysis also revealed that the immune process in the biological function category was upregulated (Fig. 2B), DEGs were analyzed by the KEGG pathway of plant-pathogen interaction (Fig. 5 2 (b)). These results indicated that elicitors initiated a cascade of defense response events and senescence-related syndromes, such as activation of catabolism and proteolysis.
The gene expression 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) was upregulated 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 downregulated or remained constant, except for the CR88 gene. The increased expression of WRKY22 and WRKY29 implies the induction of defense-related genes.
The KEGG pathway analysis suggested that some unknown elicitors and plant growth substances promote senescence-related syndromes and pathogen immune defense.

Transcriptome profiles of transcription factors
Transcription factors play essential roles in many biological aspects of life cycles. Supplementary Fig. 4 shows a heat map of the temporal distribution profile of differentially expressed transcription factor genes during protoplast isolation. There were more downregulated than upregulated transcription factors. These results indicate that the protoplast isolation procedure triggered transcriptional reprogramming.
The Plant Transcription Factor Database v4.0 (Plant TFDB) lists 1726 transcription factor genes (Supplementary Table 6, Gene number of TF family). The number of DEGs on the list was 1079. Leaf3h vs. Leaf0h had genes with significantly upregulated (p < 0.05) ERF, HFS, Myb, Trihelix, and WRKY families. Many other transcription genes showed no significant differences between protoplasts and Leaf0h, suggesting relatively early responses. On the other hand, the late-responding TF families that significantly increased in protoplasts vs. Leaf0h were ARF, bHLH, C2H2, CO-like, ERF, GRAS, Myb-related, MFY-A, TCP, YABYY, and ZFHD.
Two TF families (HSF and Trihelix) contained significantly upregulated genes in Leaf3h, but not in protoplasts vs. Leaf0h. On the other hand, several TF families containing genes whose expression was suppressed (C2H2, CO-like, GATA, GRAS, SBP, TCP, YABBY, and ZF-HD) were found in Leaf3h but not in protoplasts vs. Leaf0h. 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 ZFHD. Among these, many TF families were highly expressed in protoplasts. In other words, the gene expression of some TF families was suppressed by Leaf3h but then restored in protoplasts, but the expression level was higher than that in Leaf0h.
A total of 391 upregulated DEGs among Leaf0h, Leaf3h, and protoplasts overlap in the Venn diagram ( Fig. 2A). There are 26 transcription factors in these DEGs (Supplementary Table 1, Prot_Leaf3h_Leaf0h_EdgR_UP_TFac): five WRKY, six MYB, four MAC, and four ERF, accounting for more than 70% of the 26 upregulated transcription factor genes that overlap in the Venn diagram. Supplementary  Fig. 5a-c show the co-expression of these genes.
AT5G39610 (ANAC092, ATNAC2) encodes an NACdomain 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.97fold increase in leaves 3 h after protoplast isolation (Leaf3h) and a 2442.40-fold increase in protoplasts, compared with control leaves (Leaf0h). To explore the regulatory network administered by ANAC092, co-expression of ANAC092 and other upregulated genes in Leaf0h, Leaf3h, and protoplasts from the overlapping Venn diagram was analyzed by Cytoscape (cytoscape.org), an open-source software platform for visualizing complex networks (Fig. 6a-c). Figure 6 indicates that the NAC and WRKY transcription factors were closely co-expressed during isolation. Supplementary  Fig. 6 shows the co-expression of STP13 and transcription factor genes found as DEGs. ANAC092 was co-expressed with STP13 and was significantly upregulated throughout 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: AP2, RAV, ERF, and dehydration-responsive element-binding protein (DREB) (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) factor. ERF subfamily members that bind to ethylene response elements (EREs) respond to abiotic stress. Coexpression analysis indicated that 120, 122, and 96 transcription factor genes that were either up-or downregulated during protoplast isolation were identified as the first neighbors of DREB2B ( Supplementary Fig. 7).
The three overlapping treatments resulted in 1163 downregulated DEGs ( Fig. 2A). There were 48 transcription factors, amounting to about 4% of these overlapping treatments. In contrast to upregulated transcription factors, only one ethylene-responsive factor (ERF), two MYB, and four bHLH transcription factors were found in the three overlapping treatments. Upregulated DEGs in the three treatments did not indicate bHLH transcription factors.
The AT1G50030-encoded target of rapamycin (TOR) protein belongs to the family of phosphatidylinositol 3-kinases and is the target of the antiproliferative drug rapamycin. Figure 7a-c shows the co-expression of TOR and other first neighbor DEGs analyzed by Cytoscape, indicating that TOR expression decreased during protoplast isolation. There were no transcription factors co-expressed with TOR. Almost all of the co-expressed genes were downregulated. However, AT5G54620 encoding an ankyrin repeat-containing repeat protein was upregulated during the early stage of protoplast isolation. Fig. 6 Co-expression of ANAC092 (AT5G39610) transcription factor gene and first direct neighbor genes analyzed by Cytoscape (cytoscape.org). Gene expression is compared with (a) Leaf3h and Leaf0h, (b) Protoplast and Leaf0h, and (c) Protoplast and Leaf3h. Red and blue indicate genes induced and suppressed, respectively

Discussion
In order to use cell-wall-free protoplasts in the production of genetically modified plants, the protoplasts must be initialized before transferring useful genes and undergoing genome editing. The cell wall must be digested by cellulase and pectinase in a relatively hyperosmotic medium. Therefore, protoplasts are exposed to a great deal of stress during their isolation procedure.
Up-and downregulated genes were counted during the first 3 h of the protoplast isolation procedure. There were 3741 upregulated genes and 4482 downregulated genes (Supplementary Table 1, Leaf0h_Leaf3h_induction/Leaf0h_ Leaf3h_repression). The abundance of downregulated genes means decreased viability. Therefore, the stresses generated in the early stage of protoplast isolation play an essential role in the fate of subsequent protoplasts. There were more than 1000 downregulated than upregulated genes at the end of protoplast isolation (Supplementary Table 1, Leaf0h_Pro-toplast_induction/Leaf0h_Protoplast_repression), which indicates that protoplast viability is reduced throughout the isolation procedure.
GO overrepresentation analysis using PANTHER (Fig. 2B) showed that protoplast isolation activated genes involved in catabolic activity and metabolic processes, while it also inactivated genes involved in the metabolic process. Figure 3A (a,b) shows that many chloroplast-and photosynthesis-related genes were repressed during protoplast isolation. On the other hand, the procedure activated the expression of many mitochondrial electron transport-related genes ( Supplementary Fig. 3.1c). A decline in mitochondrial inner membrane potential triggers apoptosis in animals (Gottlieb et al. 2003;Hengartner 2000). However, the mitochondrial inner membrane potential of protoplasts was high immediately after isolation and decreased with culture time (Rajesh Kumar Tewari et al. 2012). These data suggest that chloroplasts, rather than mitochondria, are preferentially involved in the signal transduction of senescence and the apoptoticlike cell death of protoplasts.
Numerous physiological and biochemical studies have suggested the end-product inhibition of photosynthesis (Azcón-Bieto 1983;Blechschmidt-Schneider et al. 1989;Jeannette et al. 2000;Krapp et al. 1991;Quick et al. 1989). The inhibition of Calvin cycle enzymes and activation of glycolytic enzymes were attributed to photosynthesis inhibition . Photosynthetic end products of sucrose and starch in maize mesophyll protoplasts specifically and coordinately repressed photosynthetic gene expression via their promoters (Sheen 1990). Jang et al. proposed a novel role of hexokinase (HK) as a putative sugar sensor and signal transmitter in higher plants (Jang and Sheen 1994;Jang et al. 1997;Xiao et al. 2000).
Leaf senescence seems to be induced by an accumulation of carbohydrates in leaves (Quirino et al. 2000;Rolland et al. 2002;Yoshida 2003). Metabolic profiles showed that the metabolism of leaves and protoplasts differed (Fig. 4). The protoplast isolation stress triggered the induction of catabolic metabolism at 3-6 h of isolation. The sucrose content in protoplasts increased 13.78-fold compared to leaves (Leaf0h), while the glucose content decreased 5.20 × 10 -fivefold in the first 3 h of protoplast isolation (Figs. 5A and 7d; Supplementary Table 5, LC-MS-MS data summary). The gene expression of hexokinase (ATHKL1) in protoplasts, known as a sugar signal (Jang et al. 1997), increased 2.68-fold compared to that in leaves (Leaf0h) (Supplementary Table 1, All_ EdgeR.LogData). These data suggest that the photosynthetic product of sucrose promotes protoplast senescence.
The NADPH content of B. napus protoplast chloroplasts increased more than sixfold (Rajesh Kumar Tewari 2014). These data suggest that NADPH accumulation in protoplast chloroplasts results from inactivation of the Calvin-Benson cycle in protoplasts. Because thiol enzymes in the Calvin-Benson cycle are susceptible to ROS (Foyer 2018;Wise 1995), it is quite possible to inhibit the cycle under the stress of the protoplast isolation procedure.
Consequently, glucose was reduced to nearly the detection limit of LC-MS/MS during the first 3 h of protoplast isolation (Leaf3h) (Fig. 7d; Supplementary Table 5, LC-MS-MS data summary). Glucose that is not used for catabolism in protoplasts is synthesized into sucrose. Since sucrose is not loaded onto phloem, protoplasts will accumulate it. STP13 (AT5G26340) encodes a protein with high affinity, hexose-specific/H + symporter. A study of this gene's expression in programmed cell death (PCD) mutants demonstrated a tight correlation between STP13 expression and PCD (Norholm et al. 2006). Compared to Leaf0h, this gene's expression significantly increased during the early stages of protoplast isolation, by 160.65-fold in Leaf3h and 39.84-fold in protoplasts. STP13 was upregulated in cell suspension cultures undergoing PCD (Swidzinski et al. 2002), senescing leaves (Buchanan-Wollaston et al. 2005), and the two-to-sixleaf stage in the accelerated cell death 11 (acd11) mutant that develops cotyledons normally but subsequently undergoes spontaneous PCD (Brodersen et al. 2002).
The upregulation of STP13 was also observed in Arabidopsis plants treated with a fungal toxin, fumonisin B1 (FB1), and a pathogen, Pseudomonas syringae (Norholm et al. 2006). Our previous study reported that ROS were Fig. 7 Co-expression analysis of target of rapamycin (TOR; AT1G50030) gene and first direct neighbor differentially expressed genes analyzed by Cytoscape (cytoscape.org) and heat map of primary metabolites. Gene expression compared with (a) Leaf0h and Leaf3h, (b) Protoplast and Leaf0h, and (c) Protoplast and Leaf3h, and (d) heat map of primary metabolites measured by LC-MS/MS. Red and blue indicate genes induced and suppressed, respectively ◂ released outside of protoplasts and accumulated in the chloroplasts of protoplasts (Rajesh Kumar Tewari et al. 2012;Rajesh Kumar Tewari 2014;Yasuda et al. 2007). The previous results also indicated that unknown elicitors accounted for the upregulation of STP13 and accumulation of ROS in protoplasts during isolation.
Ethylene is known as a senescence hormone of plants. In order for cell wall digestion enzymes to penetrate into leaf mesophyll tissues, leaves were cut into narrow strips. This manipulation induces ethylene in wounded leaves (Watanabe et al. 2002b). Wounded leaves release jasmonic acid and oligosaccharide before ethylene production ). Salicylic acid is also produced in wounded leaves and protoplasts (Watanabe et al. 2011). These plant hormone molecules seem to promote senescence during protoplast isolation. These results suggest that the ethylene released from wounded leaves during protoplast isolation induces senescence in a different regulatory network than the ANAC029/EIN2/miR164 network. Figure 5B shows that the gene expression of acceptors for chitin and a bacterial flagellar protein increased during the early stages of protoplast isolation. Chitin is a kind of oligosaccharide. Since cellulase and pectinase in this study are crude enzymes, oligosaccharides derived from cell wall digestion and unknown proteins in the crude enzymes seemed to be candidates for eliciting ROS, subsequently inducing a hypersensitive response and immune defense. The ROS-inducing defense response finally leads to hypersensitive response (HR)-associated cell death. These results imply that sucrose, oligosaccharides, and protein elicitors contained in the cellulase/pectinase crude enzymes promote senescence and immune defense responses during the early stage of protoplast isolation. In preparation for protoplasts, cells were disrupted and formed cell debris, which was removed by a sucrose density gradient. The cell debris resulted from sucrose trigger senescence or HR-associated cell death. Therefore, the surviving protoplasts already induced senescence, which seems to be the difficulty of initializing leaf protoplasts in many plant species.
TOR-related proteins from yeast and mammals are cell growth regulators in response to nutrient availability. While AtTOR is involved in embryogenesis and is expressed in embryos, endosperm, and meristems, it negatively modulates ethylene signals. The molecular mechanism of AtTOR is likely involved in the regulation of ethylene biosynthesis by affecting 1-aminocyclopropane-1-carboxylic acid synthase (ACS) at the transcription and protein level. The regulation of autophagy by TOR and sucrose non-fermenting related protein kinase 1(SnRK1) is conserved in the plant kingdom. Tor RNAi, raptor1b, and tap46 mutants, as well as inactivation of TOR by the inhibitor AZD8055, induce autophagy, as indicated by the increased numbers of autophagosomes and higher ATG8e expression levels (Shi et al. 2018). TOR inhibition upregulated senescence-related gene expression, accompanied by upregulated expression of ethylene biosynthetic and response genes (Zhuo et al. 2020).
Our results show that the gene expression of TOR was downregulated during the isolation of protoplasts (0.367 Leaf3h/Leaf0h, 0.651 Protoplast/Leaf0h, and 1.775 Protoplast/Leaf3h) (Supplementary Table 1, All_EdgeR.LogData) and strongly support previous reports showing that TOR negatively regulated autophagy in yeast and animals (Liu and Bassham 2010). This TOR downregulation would cause ethylene-induced senescence and autophagy in protoplasts.
Cathepsin B is strongly implicated in PCD induced by abiotic stress such as UV-C, oxidative stress, and ER stress, as well as biotic stress, such as HR caused by Pseudomonas infection (Ge et al. 2016). The gene expression of cathepsin B-like protein (ATCATHB3, AT4G01610) increased throughout protoplast isolation (1.045 Leaf3h/Leaf0h, 2.304 Protoplast/Leaf0h, and 2.206 Protoplast/Leaf3h) (Supplementary Table 1, All_EdgeR.Log.Data). The cathepsin B-like protease is a candidate of proteolysis involved in cellular protein catabolic processes, especially autophagy. AT5G54620 encodes ankyrin repeat-containing proteins universally present in eukaryotes, prokaryotes, and some viruses. They function in protein-protein interactions. A total of 509 ankyrin repeats are present in 105 proteins in Arabidopsis, and the most abundant group contains proteins with ankyrin repeats and transmembrane domains (AtANKTM, cluster A) (Becerra et al. 2004). The AT5G54620-encoded protein is classified as AtANKTM. Gene expression of AT5G54620 was significantly activated during the early process of protoplast isolation (28.525 Leaf3h/Leaf0h, 3.535 Protoplast/Leaf0h, and 0.124 Protoplast/Leaf3h) (Supplementary Table 1, All_EdgeR.Log. Data, Fig. 7d). Although there are no reports showing that AtANKTM accepts protein elicitors, it is a transmembrane protein. Therefore, it likely binds unknown protein elicitors contaminated in the crude cellulase/pectinase enzyme. The AT5G54620-encoded protein may downregulate TOR.
Recently, glucose-TOR signaling has been shown to be involved in the transcriptional reprogramming of a broad range of genes (Shi et al. 2018;Xiong and Sheen 2015). In this study, the glucose content decreased to nearly the detection limit of LC-MS/MS during the first 3 h of protoplast isolation (Figs. 5A and 7d; Supplementary Table 5, LC-MS-MS data summary). Decreased glucose will inactivate glucose-TOR signaling, which consequently activates catabolic reprogramming signals and leads to ethyleneinduced senescence and autophagy in protoplasts. Following Altman et al.'s (1977) report that protoplasts isolated from Avena sativa L. leaves underwent progressive senescence, our study provides the molecular basis for protoplast senescence via TOR signaling during isolation. However, protoplast senescence seems to involve complex crosstalk signaling pathways, such as TOR signaling, ROS, elicitor signaling, wound-induced ethylene signaling, etc. Therefore, preventing these senescence signaling pathways individually may be necessary for protoplast initialization. If single-cell protoplasts can be efficiently initialized, Crisper/Cas9 can be used for genome editing by microinjection or liposomes.
Funding This research was funded by the Ichimura Foundation for New Technology.
Data availability All data are available in the manuscript and its supplementary files.