RNA-Seq Analysis of the Arabidopsis Transcriptome in Pluripotent Calli.

Plant cells have a remarkable ability to induce pluripotent cell masses and regenerate whole plant organs under the appropriate culture conditions. Although the in vitro regeneration system is widely applied to manipulate agronomic traits, an understanding of the molecular mechanisms underlying callus formation is starting to emerge. Here, we performed genome-wide transcriptome profiling of wild-type leaves and leaf explant-derived calli for comparison and identified 10,405 differentially expressed genes (> two-fold change). In addition to the well-defined signaling pathways involved in callus formation, we uncovered additional biological processes that may contribute to robust cellular dedifferentiation. Particular emphasis is placed on molecular components involved in leaf development, circadian clock, stress and hormone signaling, carbohydrate metabolism, and chromatin organization. Genetic and pharmacological analyses further supported that homeostasis of clock activity and stress signaling is crucial for proper callus induction. In addition, gibberellic acid (GA) and brassinosteroid (BR) signaling also participates in intricate cellular reprogramming. Collectively, our findings indicate that multiple signaling pathways are intertwined to allow reversible transition of cellular differentiation and dedifferentiation.

INTRODUCTION

Plant somatic cells can undergo dedifferentiation processes to give rise to pluripotent cell masses, called calli. The dedifferentiated cells are able to regenerate new organs or whole plants (Sugimoto et al., 2010). The remarkable plasticity of cellular differentiation allows plants to optimize their growth and development in response to environmental stimuli and thus overcome their sessile nature (Grafi and Barak, 2015).

Several genetic pathways are associated with the cellular dedifferentiation process. Consistent with the observations that two growth promoting hormones, auxin and cytokinin, stimulate callus formation (Ikeuchi et al., 2013), hormone signaling pathways are implicated in cellular dedifferentiation. In particular, the AUXIN RESPONSE FACTOR7 (ARF7) and ARF19 proteins regulate several members of the LATERAL ORGAN BOUNDARIES DOMAIN (LBD)/ASYMMETRIC LEAVES2-LIKE proteins, including LBD16, LBD17, LBD18, and LBD29, to promote auxin-induced callus formation (Fan et al., 2012; Okushima et al., 2007). In parallel with auxin signaling, cytokin-in-responsive type-B ARABIDOPSIS RESPONSE REGULATORs (ARRs) also positively regulate callus formation (Mason et al., 2005).

Reacquisition of cell proliferative activity is an important feature of callus formation, and the hormone signaling pathways are consistently integrated into the cell cycle program. LBD proteins control the E2 PROMOTER BINDING FACTOR a (E2Fa) transcription factor that promotes DNA replication together with DIMERIZATION PARTNER (DP) (Berckmans et al., 2011). Accordingly, overexpression of E2Fa contributes to callus formation in some plant species (Kosugi and Ohashi, 2003). Furthermore, auxin-regulated PROPORZ1 (PRZ1) suppresses KIP-RELATED PROTEIN (KRP) genes encoding cyclin-dependent kinase (CDK) inhibitors (Sieberer et al., 2003), influencing callus formation. In addition to the auxin signaling pathways, type-B ARRs may also promote the cell cycle by activating cyclins, such as CYCD3s (Argyros et al., 2008). Two APETALA2/ETHYLENE-RESPONSIVE FACTOR (AP2/ERF) transcription factors, ENHANCED SHOOT REGENERATION 1 (ESR1) and ESR2, are involved in cytokinin-dependent callus formation by activating CYCD1;1 and CYCD3;1 (Ikeda et al., 2006; Ikeuchi et al., 2013).

Embryonic or meristematic regulators are also involved in callus formation. Ectopic expression of a transcriptional activator of embryogenesis, such as LEAFY COTYLEDON 1 (LEC1), LEC2, or AGAMOUS-LIKE 15 (AGL15), results in embryonic callus formation from somatic cells even on a hormone-free medium (Gaj et al., 2005; Guo et al., 2013; Harding et al., 2003; Thakare et al., 2008). AP2/ERF transcription factors involved in embryo development, including BABY BOOM (BBM) and EMBRYOMAKER (EMK)/AINTEGUMENTA-LIKE5 (AIL5)/PLETHORA5 (PLT5), also participate in cellular dedifferentiation (Boutilier et al., 2002; Tsuwamoto et al., 2010). In addition, the homeodomain-containing transcription factor WUSCHEL (WUS), a key regulator of meristem homeostasis, plays a role in callus induction (Zuo et al., 2002).

Notably, the callus has organized structures that resemble lateral root primordia (Sugimoto et al., 2010). Transcriptome analysis indicates that gene expression profiles of the callus tissues are similar to those of lateral root meristem (Sugimoto et al., 2010). Although root developmental pathways are unequivocally associated with cellular dedifferentiation, leaf explants can be used for in vitro callus formation. Therefore, genome-wide massive reprogramming of gene expression is required for callus formation from leaf explants (He et al., 2012), although the molecular basis of the acquisition of pluripotency in leaves remains to be fully elucidated. Here, we performed transcriptional profiling using high-throughput next-generation sequencing (RNA-Seq) to highlight genes that are differentially regulated in leaves and leaf explant-derived calli. Multiple pathways were massively interconnected in the cellular dedifferentiation process, and we confirmed the biological relevance of selected biological processes. Our study provides biological insight into the intricate molecular signaling networks underlying the cellular dedifferentiation process in Arabidopsis.

MATERIALS AND METHODS

Plant materials and growth conditions

Arabidopsis thaliana (Columbia-0) seeds were germinated on Murashige and Skoog (MS) medium at 22–23°C with a 16-h light/8-h dark photoperiod. The toc1-3 (SALK_203853) (Lee et al., 2016) and Qhai1-1 and Qabi2-2 (Rodrigues et al., 2013) mutants were previously reported.

For callus induction, leaf explants of two-week-old plants were placed on callus-inducing medium (CIM) (B5 medium supplemented with 0.5 μg/ml 2,4-dichlorophenoxyacetic acid [2,4-D] and 0.05 μg/ml kinetin), followed by incubation at 22°C in the dark for 2 weeks (Fan et al., 2012). To determine the effects of paclobutrazol (PAC) and brassinolide (BL) on callus formation, 1 μM PAC (MB-P5699, MB cell, USA) and 0.1 nM and 1 nM BL (E1641, Sigma-Aldrich, USA) were added to CIM.

Quantitative real-time RT-PCR analysis

Total RNA was extracted using TRI agent (TAKARA Bio, Japan) according to the manufacturer’s recommendations. Reverse transcription (RT) was performed using Moloney Murine Leukemia Virus (M-MLV) reverse transcriptase (Dr. Protein, Korea) with oligo (dT18) to synthesize first-strand cDNA from 2 μg of total RNA. Total RNA samples were pretreated with an RNAse-free DNAse. cDNAs were diluted to 100 μl with TE buffer, and 1 μl of diluted cDNA was used for PCR amplification.

Quantitative RT-PCR reactions were performed in 96-well blocks using the Step-One Plus Real-Time PCR System (Applied Biosystems). The PCR primers used are listed in Supplementary Table 1. The values for each set of primers were normalized relative to EUKARYOTIC TRANSLATION INITIATION FACTOR 4A1 (eIF4A) (At3g13920) levels. All RT-qPCR reactions were performed in triplicate using total RNA samples extracted from three independent biological replicates. A comparative ΔΔCT method was employed to evaluate relative quantities of each amplified product in the samples. The threshold cycle (CT) was automatically determined for each reaction using the instrument software with default parameters. The specificity of each RT-qPCR reaction was determined by melting curve analysis of the amplified products using the standard method installed in the system.

mRNA-Seq data

To construct RNA libraries with the TruSeq RNA library kit, 1 μg of total RNA was used. The procedure included polyA-selected RNA extraction, RNA fragmentation, random hexamer primed reverse transcription, and 100-nt paired-end sequencing by Illumina HiSeq2000. Libraries were quantified using qPCR according to the qPCR Quantification Protocol Guide and qualified using an Agilent Technologies 2100 Bioanalyzer.

To estimate expression levels, the RNA-Seq reads were mapped to the Arabidopsis reference genome (ftp://ftp.arabidopsis.org/home/tair) using TopHat (Trapnell et al., 2009), which is capable of reporting split-read alignments across splice junctions. Transcript counts were calculated, and the relative transcript abundances were measured in FPKM (Fragments Per Kilobase of exon per Million fragments mapped) using Cufflinks.

Statistical analysis of gene expression levels

Raw data were calculated as the FPKM of each transcript in each sample using cufflinks software. We excluded transcripts with zeroed FPKM values of more than one for total samples. We added 1 to FPKM values of filtered transcripts to facilitate log2 transformation. Filtered data were transformed logarithmically and normalized using a quantile normalization method. For each transcript, we calculated fold change between case and control. Differentially expressed transcripts were determined by adjusting |fold change ≥2 of more than at least one of total comparisons.

Gene ontology (GO) term enrichment analysis

Gene functional annotation analysis for DEG list was performed using the DAVID tool (http://david.abcc.ncifcrf.gov/) (Huang et al., 2009) to understand the biological meanings behind a large list of genes. The DAVID tool provides functional annotation in over 40 annotation categories, including GO terms, protein-protein interactions, protein functional domains, disease associations, bio-pathways, sequence general features, homologies, gene functional summaries, gene tissue expressions, and literature. In the DAVID annotation system, a modified Fisher Exact p value (EASE score) is adopted to measure gene-enrichment in annotation terms. If the EASE Score is lower than 0.05 for the specific GO-term, we interpret that the given gene list is specifically associated with the GO term rather than being due to random chance. All data analysis and visualization of differentially expressed genes were conducted using R 3.1.2 software (www.r-project.org).

Immunoblot analysis

Harvested plant materials were ground in liquid nitrogen, and total cellular extracts were suspended in sodium dodecyl sulphate-polyacrylamide gel electrophoresis (SDS-PAGE) sample loading buffer. Protein samples were analyzed using SDS-PAGE (10% gels) and blotted onto Hybond-P+ membranes (Amersham-Pharmacia). Proteins were immunologically detected using anti-H3K4me3 or anti- H3K36me3 antibodies (Millipore, USA).

RESULTS AND DISCUSSION

Validation of RNA-Seq analysis

To identify the transcriptional networks controlled during the cellular dedifferentiation process, we compared the transcriptomic profiles of 3rd leaves of two-week-old wild-type plants and leaf explant-derived calli using genome-wide RNA sequencing (RNA-Seq) (Fig. 1A). RNA was isolated from a pooled callus sample rather than from biological replicates, which is similar to the averaged samples used for other studies (Holmes-Davis et al., 2005; Niederhuth et al., 2013). We sequenced the libraries on Illumina HiSeq 2000 and generated a total of 24,261,279 high quality reads. The reads were aligned onto the Arabidopsis reference genome assembly (TAIR10) (Supplementary Table 2).

Fig. 1.

RNA-Seq analysis to identify genes differentially expressed in leaves and calli. (A) Leaf explant-derived callus formation. Representative samples used for RNA-Seq analysis were photographed. Leaf explants from third-leaves of two-week-old wild-type Columbia-0 (Col-0) plants were used to induce calli on callus-inducing medium (CIM). (B) Two-way hierarchical clustering heatmap. In total, 24,497 transcripts were normalized using Z-score calculation. Differentially expressed transcripts in leaves (yellow bar) and calli (red bar) were clustered. The color key in the top left-hand corner is for colors in the heat map.

We found 5,708 up-regulated and 4,697 down-regulated genes in calli relative to leave tissues (> two-fold change in expression) (Supplementary Tables 3 and 4). These genes were divided into two groups using a hierarchical clustering algorithm (Fig. 1B). Among the genes differentially expressed, well-known positive regulators of callus formation, such as AGL15, ARF19, EMK/AIL5/PLT5, LBD18 and RESPONSE REGULATOR 1 (RR1) (Fan et al., 2012; Sakai et al., 2001; Thakare et al., 2008; Tsuwamoto et al., 2010), were included (Supplementary Table 5), supporting the reliability of our analysis.

To validate the gene expression profiles revealed by RNA-Seq analysis, we performed quantitative real-time RT-PCR (RT-qPCR) analysis and examined transcript accumulation of 10 randomly selected genes: AT1G01060, AT2G43010, AT2 G45420, AT2G46830, AT3G20810, AT4G20400, AT4G25470, AT4G33470, AT4G37650, and AT5G13790 (Supplementary Tables 3 and 4). As expected, all genes examined were significantly and differentially expressed in calli (Supplementary Fig. 1), which is highly consistent with our estimates using RNA-Seq data.

Various biological processes involved in callus formation

The global transcriptome changes in callus samples were further categorized based on their GO that suggest predicted or experimentally defined biological processes, molecular functions, and cellular components. Functional categorization of differentially regulated genes revealed that a wide variety of biological processes are associated with cellular reprogramming. In particular, functional annotations were highly enriched for functions related to protein metabolism and plant responses to biotic and abiotic stresses (Fig. 2A). Molecular functions were also widely distributed with particular enrichment related to biological macromolecule binding and transferase activity (Fig. 2B). While cellular component is largely a prediction, enriched cellular component terms included categories related to nucleus and chloroplast (Fig. 2C). Since leaf characteristics disappear during callus formation, it seems supportive of our analysis. Because biological process categories are typically derived empirically and thus tend to be more stringent, we focused on biological functions enriched in calli to understand the molecular mechanism underlying cellular dedifferentiation.

Fig. 2.

GO plant term enrichment. Genes differentially expressed in calli were categorized and annotated based on biological processes (A), molecular functions (B), and cellular components (C).

Leaf development

Although callus formation is known to follow lateral root developmental programs (Sugimoto et al., 2010), pluripotent cells can also be produced from aerial explants. Reactivation of root developmental pathways and/or genetic elimination of aerial features may accompany cellular dedifferentiation from leaf somatic cells. In support of this, the core components of Polycomb Repressive Complex 2 (PRC2), including CURLY LEAF (CLF) and SWINGER (SWN) that catalyze histone H3 lysine 27 trimethylation (H3K27me3) during callus formation to suppress leaf characteristics, contribute to stimulating callus formation from leaf explants (He et al., 2012).

We found 149 genes involved in leaf development were down-regulated in calli compared with leaves (Supplementary Table 6). In particular, predominant members of key transcription factor families involved in leaf development were included, including TEOSINTE BRANCHED1, CYCLOIDEA, AND PCF (TCP), BEL1-LIKE HOMEODOMAIN (BLH), and NGATHA (NGA) (Kumar et al., 2007; Nag et al., 2009; Trigueros et al., 2009). The Arabidopsis genome contains 13 class I and 11 class II TCP genes (Koyama et al., 2010; Tao et al., 2013). Consistent with the fact that the class II members have been characterized as major regulators of leaf development (Tao et al., 2013), expression levels of TCP2, TCP3, TCP4, and TCP10 genes were reduced in calli (Table 1 and Supplementary Table 6). In addition, BLH2/SAWTOOTH1 (SAW1) and BLH4/SAW2, which establish shoot architecture together with SHOOTMERISTEMLESS (STM) and BREVIPEDICELLUS (BP) (Kumar et al., 2007), were also down-regulated (Table 1 and Supplementary Table 6). NGA transcription factors are also implicated in leaf development, as shown in quadruple nga1 nga2 nga3 nga4 mutants that display abnormal leaf development (Trigueros et al., 2009), and accordingly, the NGA2 and NGA3 genes were significantly repressed in calli (Table 1 and Supplementary Table 6).

Table 1.

Transcript profiles of key genes involved in selected biological processes. Genes differentially expressed in calli were functionally categorized according to gene ontology (GO) at Arabidopsis Information Resource. The key components of selected biological processes were shown to estimate their biological relevance. Genes were rank ordered in each category by fold change (FC) in expression.

Transcript	Description	Calli/Leaves (FC)
Leaf development
AT4G18390.1	TEOSINTE BRANCHED 1, CYCLOIDEA AND PCF TRANSCRIPTION FACTOR 2 (TCP2)	−2.39
AT3G61970.1	NGATHA2 (NGA2)	−3.88
AT1G53230.1	TCP FAMILY TRANSCRIPTION FACTOR 3 (TCP3)	−4.28
AT1G01030.1	NGATHA3 (NGA3)	−5.29
AT3G15030.1	TCP FAMILY TRANSCRIPTION FACTOR 4 (TCP4)	−7.81
AT2G23760.1	BEL1-LIKE HOMEODOMAIN 4 (BLH4)	−10.58
AT4G36870.1	BEL1-LIKE HOMEODOMAIN 2 (BLH2)	−14.27
AT2G31070.1	TCP DOMAIN PROTEIN 10 (TCP10)	−30.50
Circadian clock
AT5G61380.1	TIMING OF CAB EXPRESSION 1 (TOC1)	9.16
AT2G25930.1	EARLY FLOWERING 3 (ELF3)	3.64
AT2G40080.1	EARLY FLOWERING 4 (ELF4)	3.23
AT2G46830.1	CIRCADIAN CLOCK-ASSOCIATED 1 (CCA1)	−18.11
AT1G01060.1	LATE ELONGATED HYPOCOTYL (LHY)	−25.42
Response to virus
AT3G03300.2	DICER-LIKE 2 (DCL2)	2.67
AT2G27040.1	ARGONAUTE 4 (AGO4)	2.13
AT1G48410.2	ARGONAUTE 1 (AGO1)	2.12
Response to bacterium
AT3G52430.1	PHYTOALEXIN DEFICIENT 4 (PAD4)	11.79
AT3G48090.1	ENHANCED DISEASE SUSCEPTIBILITY 1 (EDS1)	8.15
AT4G31800.1	WRKY DNA-BINDING PROTEIN 18 (WRKY18)	7.14
AT5G45110.1	NPR1-LIKE PROTEIN 3 (NPR3)	7.03
AT3G56400.1	WRKY DNA-BINDING PROTEIN 70 (WRKY70)	4.23
AT5G65210.1	TGACG SEQUENCE-SPECIFIC BINDING PROTEIN 1 (TGA1)	3.05
AT2G25000.1	WRKY DNA-BINDING PROTEIN 60 (WRKY60)	2.55
AT4G31550.2	WRKY DNA-BINDING PROTEIN 11 (WRKY11)	2.29
AT2G24570.1	WRKY DNA-BINDING PROTEIN 17 (WRKY17)	2.24
AT4G23810.1	WRKY DNA-BINDING PROTEIN 53 (WRKY53)	2.20
AT1G22070.1	TGA1A-RELATED GENE 3 (TGA3)	2.02
Response to osmotic stress
AT5G63650.1	SNF1-RELATED PROTEIN KINASE 2.5 (SNRK2.5)	37.84
AT5G66880.1	SNF1-RELATED PROTEIN KINASE 2.3 (SNRK2.3)	11.29
AT3G50500.1	SNF1-RELATED PROTEIN KINASE 2.2 (SNRK2.2)	7.84
AT1G49720.1	ABSCISIC ACID RESPONSIVE ELEMENT-BINDING FACTOR 1 (ABF1)	5.19
AT5G59220.1	HIGHLY ABA-INDUCED PP2C GENE 1 (HAI1)	3.58
AT4G33950.2	OPEN STOMATA 1 (OST1)	3.32
AT4G26080.1	ABA INSENSITIVE 1 (ABI1)	2.05
Response to oxidative stress
AT1G08830.1	COPPER/ZINC SUPEROXIDE DISMUTASE 1 (CSD1)	4.97
AT2G48150.1	GLUTATHIONE PEROXIDASE 4 (GPX4)	4.74
AT2G43350.1	GLUTATHIONE PEROXIDASE 3 (GPX3)	3.88
AT1G20630.1	CATALASE 1 (CAT1)	3.72
AT1G63460.1	GLUTATHIONE PEROXIDASE 8 (GPX8)	3.57
AT4G11600.1	GLUTATHIONE PEROXIDASE 6 (GPX6)	2.38
AT3G10920.2	MANGANESE SUPEROXIDE DISMUTASE 1 (MSD1)	−2.00
AT4G35000.1	ASCORBATE PEROXIDASE 3 (APX3)	−2.17
AT5G23310.1	FE SUPEROXIDE DISMUTASE 3 (FSD3)	−2.20
AT4G32320.1	ASCORBATE PEROXIDASE 6 (APX6)	−3.13
AT1G07890.6	ASCORBATE PEROXIDASE 1 (APX1)	−3.19
AT5G18100.2	COPPER/ZINC SUPEROXIDE DISMUTASE 3 (CSD3)	−5.36
AT2G31570.1	GLUTATHIONE PEROXIDASE 2 (GPX2)	−6.78
AT4G25100.1	FE SUPEROXIDE DISMUTASE 1 (FSD1)	−7.89
AT2G25080.1	GLUTATHIONE PEROXIDASE 1 (GPX1)	−16.83
AT4G31870.1	GLUTATHIONE PEROXIDASE 7 (GPX7)	−17.40
AT4G35090.2	CATALASE 2 (CAT2)	−42.58
AT4G09010.1	ASCORBATE PEROXIDASE 4 (APX4)	−177.58
Polysaccharide metabolic process
AT4G09020.1	ISOAMYLASE 3 (ISA3)	11.23
AT1G10760.1	A-GLUCAN WATER, DIKINASE/STARCH EXCESS 1 (GWD1/SEX1)	2.63
AT5G26570.2	PHOSPHOGLUCAN WATER DIKINASE (PWD)	2.02
AT5G64740.1	CELLULOSE SYNTHASE 6 (CESA6)	−2.11
AT4G39350.1	CELLULOSE SYNTHASE 2 (CESA2)	−5.31
AT5G09870.1	CELLULOSE SYNTHASE 5 (CESA5)	−8.84
Gibberellin signaling
AT5G17490.1	RGA-LIKE PROTEIN 3 (RGL3)	10.63
AT3G63010.1	GA INSENSITIVE DWARF1B (GID1B)	6.56
AT3G03450.1	RGA-LIKE 2 (RGL2)	3.43
AT3G05120.1	GA INSENSITIVE DWARF1A (GID1A)	3.42
AT1G66350.1	RGA-LIKE 1 (RGL1)	2.86
AT5G27320.1	GA INSENSITIVE DWARF1C (GID1C)	2.58
Brassinosteroid signaling
AT4G39400.1	BRASSINOSTEROID INSENSITIVE 1 (BRI1)	2.34
AT4G33430.1	BRI1-ASSOCIATED RECEPTOR KINASE (BAK1)	2.20
Chromatin organization
AT1G26760.1	SET DOMAIN PROTEIN 35 (SDG35)	4.73
AT5G24330.1	ARABIDOPSIS TRITHORAX-RELATED PROTEIN 6 (ATXR6)	3.20
AT1G76710.1	SET DOMAIN GROUP 26 (SDG26)	2.66
AT5G09790.2	ARABIDOPSIS TRITHORAX-RELATED PROTEIN 5 (ATXR5)	2.61

In addition, a significant portion of genes related to leaf identity, lateral organ meristem establishment, chloroplast development, and photosynthesis was substantially repressed in calli (Supplementary Table 6). These results indicate that genetic components establishing leaf identity are suppressed during callus formation to turn off differentiated cell identity.

Circadian clock

The circadian clock generates the biological rhythm with a period of ∼24 h, synchronizing internal processes with environmental cycles. Correct matching of the circadian period with external cues ensures plant fitness, and consistently, nearly all stages of plant development are under the control of the circadian clock (Yakir et al., 2007). Previous findings have shown that a biological rhythm retains in dedifferentiated Arabidopsis protoplast cells (Nakamichi et al., 2004), suggesting a possible role of the circadian clock in cellular differentiation and dedifferentiation.

We found that at least 138 circadian-related genes were differentially regulated in calli (Supplementary Table 7). Notably, many core clock oscillator components in multiple-interlocked circadian feedback loops, including CIRCADIAN CLOCK-ASSOCIATED 1 (CCA1), LATE ELONGATED HYPOCOTYL (LHY), TIMING OF CAB EXPRESSION 1/ PSEUDO-RESPONSE REGULATOR 1 (TOC1/PRR1), EARLY FLOWERING 3 (ELF3), and ELF4, were included in the lists (Table 1 and Supplementary Table 7), supporting the requirements of endogenous circadian networks in the establishment and maintenance of pluripotency.

To estimate the relevance of circadian functions in cellular dedifferentiation, we employed a genetic mutant that harbors a defect in the core clock oscillator TOC1, the toc1-3 mutant, and compared its callus formation rate with wild-type. Callus formation was reduced when derived from toc1-3 mutant leaves (Fig. 3). Fresh weight measurements revealed that the toc1-3 mutant showed a 1.2-fold reduction in callus formation capability (Fig. 3). Considering that leaf explants were placed in darkness for 2 weeks to induce calli, it is unlikely that circadian oscillation persists over the period of cellular reprogramming. Instead of circadian fluctuations, clock-controlled signaling networks might be associated with robust callus formation.

Fig. 3.

Callus formation capability of the toc1-3 mutant. (A) Callus formation. Leaf explants from third-leaves of two-week-old plants were used to induce calli on CIM (n > 30). Plates were incubated for 2 weeks under continuous dark conditions and photographed. Scale bars = 5 mm. (B) Fresh weight measurement. Thirty calli were collected to measure fresh weight. Bars indicate the standard error of the mean. Statistically significant differences between wildtype and toc1-3 mutant are indicated by asterisks (Student’s t-test, *P < 0.05).

Plant responses to biotic stress

Since plants are sessile organisms, they have evolved elaborate dedifferentiation abilities to ensure developmental plasticity. Biotic stress occasionally stimulates cell fate reprogramming. For instance, Agrobacterium tumefaciens causes crown gall tumor in a wide-range of host plants (Gohlke and Deeken, 2014). Infection with Rhizobium meliloti also leads to cellular dedifferentiation in Medicago root cortical cells (Charon et al., 1997). A few genetic components responsible for biotic stress-induced cellular dedifferentiation have been identified. The wounding-responsive AP2/ERF transcription factor WOUND-INDUCED DEDIFFERENTIATION 1 (WIND1) and its homologs, WIND2, WIND3 and WIND4, are possible candidates for the regulation of callus formation upon exposure to pathogens or herbivory (Iwase et al., 2011).

In addition to the biotic stress-stimulated cellular dedifferentiation process (Grafi, 2004), calli are known to retain a high protective ability against various microbial invaders, suggesting intimate bidirectional regulation. Our RNA-Seq analysis showed that pivotal defense genes were up-regulated in calli. RNA silencing is a representative antiviral defense mechanism in plants. DICER-LIKE (DCL) RNase III enzymes are primarily responsible for catalyzing RNA processing (Henderson et al., 2006). Arabidopsis has four DICER-like proteins, DCL1–DCL4, with functional specificity. Among them, DCL2 is mainly responsible for the generation of virus-derived siRNAs (Gasciolli et al., 2005). Small RNA duplexes are then incorporated into RNA-induced silencing complex (RISC) to target viral RNA for processing (Pumplin and Voinnet, 2013). A central component of RISC complex is the ARGONAUTE (AGO) protein, and in Arabidopsis, ten AGO proteins (AGO1-AGO10) are possible candidates for RISC formation (Garcia-Ruiz et al., 2015). Our RNA-Seq analysis revealed that 60 virus-responsive genes, including DCL2, AGO1 and AGO4, showed higher expression in calli than in leaves (Table 1 and Supplementary Table 8).

Plants have also developed strong innate immunity against bacterial pathogens. As a first layer of defense, plants activate pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) upon recognition of PAMPs by the plant pattern recognition receptor (PRR) (Kim et al., 2008; Nicaise et al., 2009). Signaling cascades built-up with MITOGEN-ACTIVATED PROTEIN KINASEs (MAPKs), CALCIUM-DEPENDENT PROTEIN KINASES (CDPKs), and WRKY transcription factors are subsequently activated (Eulgem and Somssich, 2007; Rasmussen et al., 2012). In addition, plants have a second barrier to pathogen invasion, which stimulates effector-triggered immunity (ETI). Pathogen effector molecules are sensed by various plant resistance (R) proteins, and the recognition triggers hypersensitive response (HR) and systemic acquired resistance (SAR) mediated by SA biosynthesis (Dempsey et al., 2011). Notably, 76 bacterial-responsive genes that mediate PTI and/or ETI, such as WRKY11, WRKY17, WRKY18, WRKY53, WRKY60, WRKY70, ENHANCED DISEASE SUSCEPTIBILITY 1 (EDS1), and PHYTOALEXIN DEFICIENT 4 (PAD4), were up-regulated in calli (Table 1 and Supplementary Table 9). Accordingly, SA signaling genes, such as TGACG SEQUENCE-SPECIFIC BINDING PROTEIN 1 (TGA1), TGA3, and NPR1-LIKE PROTEIN 3 (NPR3), were also influenced (Table 1 and Supplementary Table 9), possibly contributing to maintaining disease-free dedifferentiated cells.

Plant responses to abiotic stress

Environmental stress stimulates cellular dedifferentiation (Grafi et al., 2011). Transcriptome profiling has revealed that differentiated organs experiencing environmental challenges, such as drought and continuous dark, have gene expression profiles similar to those of dedifferentiated protoplast cells (Grafi et al., 2011). Consistently, stress-related genes form the molecular signature for pluripotent animal cells (Ramalho-Santos et al., 2002). These observations suggest that abiotic stress-responsive transcription factors may play a key role in cellular reprogramming. For example, the NAC DOMAIN-CONTAINING PROTEIN 2 (NAC2) transcription factor is not only responsive to abiotic stresses (Nuruzzaman et al., 2013) but also regulates meristem proliferation activity (Florentin et al., 2013). Thus, NAC2 acts with potential molecular crosstalk of abiotic stress and stem cell signaling and contributes to the reacquisition of stem cell fate in Arabidopsis somatic cells under unfavorable environmental conditions (Florentin et al., 2013).

We found that a significant number of abiotic stress- and/or abscisic acid (ABA)-responsive genes are differentially expressed in calli (Supplementary Table 10). Osmotic stress is largely dependent on ABA, and consistently, osmotic stress primarily influences genes encoding ABA sensory and up stream components, such as REGULATORY COMPONENT OF ABA RECEPTORs (RCARs), Type 2C protein phosphatases (PP2Cs), SNF1-related protein kinase 2s (SnRK2s), and abscisic acid-responsive element-binding factors (ABFs) (Chan, 2012). Notably, 152 osmotic stress-responsive genes were up-regulated in calli (Supplementary Table 10), and most notably ABA-INSENSITIVE 1 (ABI1), ABF1, HIGHLY ABA-INDUCED1 (HAI1), SnRK2.2, SnRK2.3, SnRK2.5, and SnRK2.6/OPEN STOMATA 1 (OST1) (Table 1 and Supplementary Table 10). Given that osmotic stress induces somatic embryogenesis in Arabidopsis (Ikeda-Iwai et al., 2003), these genes might play a role in switching cell fates under environmentally unfavorable conditions.

To address whether the ABA signaling genes are indeed involved in switching cell fates, we employed two quadruple mutants, Qhai1-1 (lacking HAI1, PP2CA, HOMOLOGY TO ABI1 (HAB1), and ABI1; Rodrigues et al., 2013) and Qabi2-2 (lacking ABI2, PP2CA, HAB1, and ABI1; Rodrigues et al., 2013), and induced callus formation on CIM. Both quadruple mutants exhibited decreased callus formation capability relative to wild-type plants (Fig. 4), indicating that ABA signaling homeostasis is necessary for proper callus formation. While further detailed analysis is required to unravel the underlying molecular mechanisms, ABA upstream regulators might be a molecular linker of stress adaptation and cellular dedifferentiation.

Fig. 4.

Callus formation of ABA signaling mutants. (A) Callus formation. Leaf explants from the third leaves of two-week-old plants were used to induce calli on CIM (n > 30). Plates were incubated for 2 weeks under continuous dark conditions and photographed. Scale bars = 5 mm. (B) Fresh weight measurement. Thirty calli were collected to measure fresh weight. Bars indicate the standard error of the mean. Statistically significant differences between wild-type and mutants are indicated by asterisks (Student’s t-test, *P < 0.05).

Reactive oxygen species (ROS) are important signaling molecules in stress responses and cellular differentiation. In addition to their stress-responsive defense functions (Tsukagoshi et al., 2010), ROS also regulate dedifferentiation as exemplified by the promotive effects of hydrogen peroxide (H2O2) on cellular dedifferentiation in Lycium barbarum (Cui et al., 2002). Optimal balances of ROS accumulation are achieved by multiple antioxidative systems. Enzymatic detoxification is catalyzed by superoxide dismutase (SOD), ascorbate peroxidase (APX), catalase (CAT), glutathione peroxidase (GPX), and peroxiredoxin (PRX) (Caverzan et al., 2012). The Arabidopsis genome encodes seven SOD proteins classified depending on their metal cofactors: three Cu/Zn SOD proteins [COPPER/ZINC SUPEROXIDE DISMUTASE 1 (CSD1), CSD2, CSD3], three Fe SOD proteins [Fe SUPEROXIDE DISMUTASE 1 (FSD1), FSD2 and FSD3], and one Mn SOD protein [MANGANESE SUPEROXIDE DISMUTASE 1 (MSD1)] (Xing et al., 2013). There are also nine APX (APX1-9), three CAT (CAT1-3) and eight GPX (GPX1-8) genes (Caverzan et al., 2012; Passaia et al., 2014; Xing et al., 2007). In calli, CSD1, CAT1, GPX3, GPX4, GPX6 and GPX8 were up-regulated, whereas APX1, APX3, APX4, APX6, CAT2, CSD3, FSD1, FSD3, GPX1, GPX2, GPX7 and MSD1 were repressed (Table 1 and Supplementary Table 11). Furthermore, many peroxidase and reductase genes were also differentially regulated (Supplementary Table 11), supporting the importance of ROS homeostasis during cellular dedifferentiation.

While it seems clear that abiotic stress triggers cell fate changes, the molecular components responsible for the process are largely unknown. The genes that we proposed in this study might play a role as molecular linkers of stress-induced cellular dedifferentiation. Future studies will also provide biological insights as to why plants dedifferentiate in suboptimal conditions and how pluripotent cells improve plant adaptation ability to environmental challenges.

Carbohydrate metabolism

The plant cell wall is an important cellular structure that provides structural support and maintains cellular differentiation states (Keegstra, 2010). In support of this, treatments with cell wall-degrading enzymes produce dedifferentiated protoplast cells (Grafi, 2004). Moreover, endo-1,4-β-d-glucanase KORRIGAN 1/TUMOROUS SHOOT DEVELOPMENT 1 (KOR1/TSD1), which is involved in cellulose synthesis, also participates in the maintenance of somatic cell identity. The tsd1 loss-of-function mutants display callus-like structures with reduced cell wall formation on a hormone-free medium (Frank et al., 2002).

The plant cell wall is composed of cellulose, hemicellulose, lignin, and pectin. In particular, cellulose, a linear chain of β (1–4)-linked D-glucose, is the main component of the cell wall and comprises approximately 50% of the total cell wall mass (Endler and Persson, 2011). In Arabidopsis, 10 CELLULOSE SYNTHASE (CESA) proteins (CESA1-CESA10) are supposed to catalyze cellulose synthesis (Endler and Persson, 2011). Here, we found that 29 genes involved in polysaccharide metabolism, including CESA2, CESA5, and CESA6, were down-regulated in calli (Table 1 and Supplementary Table 12). The genes encoding hemicellulose-metabolizing enzymes, such as xyloglucan endotransglucosylases/hydrolases (XTHs), were also included (Supplementary Table 12).

The storage polysaccharide starch is the energy source in plants (Kötting et al., 2005). Starch is converted into simple soluble sugars by the action of hydrolytic enzymes, such as amylases (Smith et al., 2005; Streb et al., 2012). In Arabidopsis, additional proteins that stimulate starch metabolism have been identified; α-GLUCAN WATER, DIKINASE/STARCH EXCESS 1 (GWD1/SEX1) and PHOSPHOGLUCAN, WATER DIKINASE (PWD), which phosphorylate C6 and C3-position of glucosyl residues, respectively, catalyze starch hydrolysis (Santelia et al., 2011). Notably, ISOAMYLASE 3 (ISA3), GWD1, and PWD were up-regulated in calli (Table 1 and Supplementary Table 12), whereas many starch synthase genes were repressed (Supplementary Table 12). These results suggest that active cell proliferation requires a sufficient energy supply (Narbonne and Roy, 2006; Skylar et al., 2011), and simple sugar production might drive the cellular dedifferentiation process.

Hormone signaling

Two phytohormones, auxin and cytokinin, play a crucial role in cellular dedifferentiation in a variety of plant species. Since the coordinated actions of phytohormones are important for various aspects of plant growth and development (Vanstraelen and Benková, 2012), we hypothesized that additional hormones, such as ABA, gibberellin (GA), ethylene (ET), brassinosteroid (BR), jasmonic acid (JA), and salicylic acid (SA), may be implicated in callus formation.

Our RNA-Seq analysis revealed that differential expression of GA and BR signaling components in calli were noticeable. In particular, 22 GA-related genes, which encode GA metabolic enzymes, GA receptors, and DELLA proteins, were differentially regulated in calli (Table 1 and Supplementary Table 13). Considering that GA shows a functional overlap with auxin in plant growth and development (Fu and Harberd, 2003), crosstalk between GA and auxin may underlie the cellular dedifferentiation process. To prove this hypothesis, we employed a potent chemical inhibitor of GA biosynthesis, PAC, and examined its effects on callus formation capability. Treatment with PAC resulted in significant reduction in callus formation ability, which reaches up to approximately 70%, supporting a positive function of GA in cellular dedifferentiation (Fig. 5A).

Fig. 5.

Effects of PAC and BL on callus formation. Leaf explants from third-leaves of two-week-old wild-type plants were used to induce calli on CIM in the presence of 1 μM paclobutrazol (PAC) (A) and 0.1 nM and 1 nM brassinolide (BL) (B). Plates were incubated for 2 weeks and photographed (left panels). Thirty calli were collected to measure fresh weight. Bars indicate the standard error of the mean. Statistically significant differences compared with the mock-treated sample are indicated by asterisks (Student’s t-test, *P < 0.05) (right panels).

The plant steroid hormone BR is closely associated with the processes of cell proliferation and differentiation (González-García et al., 2011). The BR biosynthetic gene DWARF7-deficient dwarf7-1 mutant consistently displays reduced callus formation, possibly due to reduced cell division (Cheon et al., 2010). Our analysis further supported the previous reports. BR signalling genes, including BRASSINOSTEROID INSENSITIVE 1 (BRI1) and BRI1-ASSOCIATED RECEPTOR KINASE1 (BAK1), were up-regulated in calli (Table 1 and Supplementary Table 14). In addition, we also noticed that several key BR signaling genes, including BRI1-EMS-SUPPRESSOR 1 (BES1) and BRI1 SUPPRESSOR 1 (BRS1), were repressed in calli (Supplementary Table 14). Consistent with the controversial results, opposite effects of BR on callus formation were observed depending on BL concentrations applied. Callus formation was stimulated by low concentration of BL, but inhibited at higher concentration (Fig. 5B). Different modes of BR actions would underlie the balanced cellular reprogramming. Collectively, these results indicate that additional hormones play a supplemental role in inducing callus formation.

Chromatin organization

Chromatin structure determines gene expression by regulating the accessibility of transcriptional machinery. In higher eukaryotes, global gene expression changes facilitated by epigenetic regulation are associated with cellular differentiation status. For instance, in mammals, differentiated cells have a closed chromatin state, whereas pluripotent dedifferentiated cells have an open chromatin state (Gaspar-Maia et al., 2011).

Genome-wide reorganization of chromatin structure is also crucial for cellular dedifferentiation in plants. The PRC2 and PRC1 complexes coordinate establishment and maintenance of repressive histone mark H3K27me3 to control cellular differentiation (Lafos et al., 2011). Mutations in PRC2 or PRC1 components lead to spontaneous callus formation from root explants with ectopic expression of LEC1, LEC2, AGL15, BBM, WUS, and WUSCHEL RELATED HOMEOBOX 5 (WOX5) (Bouyer et al., 2011; Bratzel et al., 2010; Chen et al., 2010).

While H3K27 methylation is involved in callus formation, the relevance of other histone marks in this process remains to be fully elucidated. The activation histone marks, H3K4me2/3 and H3K36me2/3, are associated with higher gene expression (Guo et al., 2010; Yang et al., 2012), and the TRITHORAX GROUP (TrxG) and ABSENT, SMALL OR HOMEOTIC DISCS 1 (ASH1) proteins are responsible for the establishment of the activation marks (He et al., 2013). The Arabidopsis genome has multiple homologs: five ARABIDOPSIS HOMOLOG OF TRITHORAX proteins (ATX1-5), seven ATX-RELATED proteins (ATXR1-7), four ASH1 HOMOLOG proteins (ASHH1-4), and three ASH1-RELATED proteins (ASHR1-3) (Berr et al., 2009). We found that at least 115 genes functionally associated with chromatin remodeling were differentially regulated in calli (Supplementary Table 15). Notably, expression of ATXR1/SDG35, ATXR5/SDG15, ATXR6/SDG34, and ASHH1/SDG26 were elevated (Table 1 and Supplementary Table 15). Accordingly, immunoblot analysis showed that H3K4me3 and H3K36me3 levels were significantly increased in calli relative to vegetative tissues (Fig. 6).

Fig. 6.

Total H3K4me3 and H3K36me3 levels in calli. Leaf explants from third-leaves of two-week-old wild-type plants were used to induce calli on CIM. H3K4me3 (A) and H3K36me3 (B) levels (arrowheads in each) were detected immunologically using the corresponding antibodies. A part of a Coomassie blue-stained gel (C) is shown as a loading control (left panels). Bands from three blots were quantified and averaged using Image J software. Bars indicate the standard error of the mean. Statistically significant differences between values of leaf and other vegetative tissue samples are indicated by asterisks (Student’s t-test, *P < 0.05) (right panels).

The activation marks established by TrxG and ASH1 counteract, in general, polycomb group (PcG) repression. Given the roles of PcG in cellular dedifferentiation in plants (Bouyer et al., 2011; Bratzel et al., 2010; Chen et al., 2010), it seems meaningful that TrxG functions also underlie the process. In addition, considering the large number of chromatin remodeling genes differentially regulated in calli (Supplementary Table 15), epigenetic regulation is a key regulatory layer for cellular dedifferentiation.

In conclusion, by comparing RNA-Seq read count, we assessed the relative accumulation of RNAs present in the leaf and callus and studied the molecular processes controlling cellular dedifferentiation. We found and confirmed that a variety of biological processes are engaged during callus formation. In particular, molecular components related to shoot development, circadian clock regulation, biotic and abiotic stress responses, nutrient metabolism, and chromatin modification are possibly associated. These results demonstrate the complexity of cellular reprogramming and that the multiple pathways are integrated to facilitate callus formation. Here, we provide valuable initial resources for investigating cellular dedifferentiation in many crop plants and pave the way for improving crop transformation efficiency for genetic engineering.