Transcriptome comparison between pluripotent and non-pluripotent calli derived from mature rice seeds

In vitro plant regeneration involves a two-step practice of callus formation and de novo organogenesis. During callus formation, cellular competence for tissue regeneration is acquired, but it is elusive what molecular processes and genetic factors are involved in establishing cellular pluripotency. To explore the mechanisms underlying pluripotency acquisition during callus formation in monocot plants, we performed a transcriptomic analysis on the pluripotent and non-pluripotent rice calli using RNA-seq. We obtained a dataset of differentially expressed genes (DEGs), which accounts for molecular processes underpinning pluripotency acquisition and maintenance. Core regulators establishing root stem cell niche were implicated in pluripotency acquisition in rice callus, as observed in Arabidopsis. In addition, KEGG analysis showed that photosynthetic process and sugar and amino acid metabolism were substantially suppressed in pluripotent calli, whereas lipid and antioxidant metabolism were overrepresented in up-regulated DEGs. We also constructed a putative coexpression network related to cellular pluripotency in rice and proposed potential candidates conferring pluripotency in rice callus. Overall, our transcriptome-based analysis can be a powerful resource for the elucidation of the molecular mechanisms establishing cellular pluripotency in rice callus.


Results
Transcriptome analysis identifies DEGs in pluripotent rice callus. While several molecular players involved in pluripotency acquisition have been identified in Arabidopsis, we wanted to know whether the gene networks are conserved in monocot rice plants. In addition, since we were able to distinguish pluripotent and non-pluripotent calli which were derived from mature rice seeds, we further tried to elucidate genetic candidates that regulate pluripotency acquisition in rice. Emerging calli could be sorted into pluripotent and non-pluripotent calli at 28 days after CIM incubation (DAC), based on their phenotypes (Fig. 1). Pluripotent calli with a globular-shaped compact cluster of yellowish cells had cellular competence for tissue regeneration, whereas non-pluripotent calli with semitransparent dark yellow color had a limited ability for de novo shoot regeneration (Fig. 1). The separation was reliable and reproducible, and the pluripotent and non-pluripotent calli divided at 28 DAC showed distinguishable regenerative potential during subsequent incubation on SIM (Supplementary Fig. S1 and Supplementary Table S1).
We collected pluripotent and non-pluripotent calli at 28 DAC and extracted mRNAs to perform Illumina sequencing. Using the Illumina Hiseq-2000, we obtained 76,653,249 and 75,047,011 reads from pluripotent and non-pluripotent calli, respectively. Among these reads, 59,747,844 and 57,835,842 reads accounting for 77.9% and 77.1% of the total reads, respectively, were properly paired and mapped to the reference genome of japonica (MSU v.7.0; https ://phyto zome.jgi.doe.gov/pz/porta l.html) ( Table 1). Although transcriptome of the two different types of calli was considerably similar (Fig. 2a), principal component analysis (PCA) showed that the different callus types were also significantly separated along the first principal component, which explains 99.8% of the variability (Fig. 2b).
Based on the adopted cut-off (fold-change > 2 and p value < 0.05) (Fig. 2c), 60 up-regulated genes (Supplementary Table S2) and 184 down-regulated genes (Supplementary Table S3) in pluripotent callus relative   No. of total reads  34,957,214  41,696,035  76,653,249  38,386,915  36,660,096  75,047,011   No. of mapped reads  34,957,214  41,696,035  76,653,249  38,386,915  36,660,096  75,047,011   No. of properly paired and  mapped reads  27,753,148  31,994,696  59,747,844  29,488,630  28,347, Fig. S2 and Supplementary Table S4). Since a limited cell population may be related to pluripotency, observation of the moderate number of DEGs was reasonable. We performed a visual inspection of the hierarchical clustering results to identify major subgroups and clusters (Fig. 2d). Collectively, we successfully distinguished pluripotent and non-pluripotent calli to analyze transcriptome possibly related to cellular pluripotency.

KEGG analysis reveals metabolic processes enriched in DEGs.
Changes in gene expression underlie differences in biological and metabolic functions between pluripotent and non-pluripotent calli. Gene ontology (GO) analysis of DEGs did not enrich a certain category. We thus alternatively explored the Kyoto Encyclopedia of Genes and Genomes (KEGG) terms that were enriched for up-and down-regulated genes in pluripotent calli. Within the group of up-regulated genes, we identified overrepresented KEGG categories for metabolic pathways (map01100), biosynthesis of secondary metabolites (map01110), fatty acid elongation (map00062), α-linolenic acid metabolism (map00592), and ascorbate metabolism (map00053) (Fig. 3a). This might reflect that pluripotent calli possibly promotes accumulation of secondary metabolites, such as ascorbic acid and antioxidants, as well as primary metabolites including fatty acids and lipids, compared with non-pluripotent calli (Fig. 3a), consistent with the fact that reactive oxygen species (ROS) and lipid molecules are closely associated with pluripotency acquisition in callus 22 .
In the case of the down-regulated genes, KEGG terms related to metabolic pathways (map01100), diterpenoid biosynthesis (map00904), phenylpropanoid biosynthesis (map00940), carbohydrate metabolic processes (map01200 and map00500), and hormone metabolism and signaling (map04075 and map04010) were enriched (Fig. 3b). The down regulation of these gene sets may be related to global alterations in plant metabolism and hormone signaling. Furthermore, consistent with the fact that callus tissues resemble to root primordium 7,8 , photosynthesis-related genes (map00195) were also included in this group.
Primary and secondary metabolism and hormone signaling are related to pluripotency acquisition. We next wanted to know molecular factors underlying pluripotency establishment in rice callus. To this end, we first checked whether genes known to regulate cellular pluripotency were included in DEGs of pluripotent rice calli. As a result, the PLT1 gene encoding an auxin-responsive transcription factor responsible for establishing root stem cell niche was up-regulated in pluripotent rice calli ( Fig. 4a; Supplementary Table S2) 11,15 . Furthermore, a HSF transcription factor was also included (Fig. 4a), consistent with a possible role of HSFs in plant regeneration 23 . Other key transcription factors for various developmental processes were found in DEGs of pluripotent calli (Fig. 4a), implying that they might act as potential upstream regulators of pluripotency acquisition. . KEGG pathway enrichment analysis. Enrichment of KEGG terms of up-regulated genes (a) and down-regulated genes (b) in pluripotent callus relative to non-pluripotent callus was shown. X-axis indicates log-transformed p value. Size of circle represents difference between observed and expected ratio for a specific pathway. Statistical significance was tested by binomial test.
Scientific Reports | (2020) 10:21257 | https://doi.org/10.1038/s41598-020-78324-z www.nature.com/scientificreports/ We expanded our analysis and examined more genes implicated in diverse aspects of primary and secondary metabolism in plants. Genes related to photosynthesis and primary metabolism for sugar, lipid, and amino acids were particularly enriched in our DEG lists ( Fig. 4b-f), and they showed dynamic changes in pluripotent callus relative to non-pluripotent callus. Phenylpropanoid and glucosinolate metabolism genes were also differentially expressed (Fig. 4g,h), and most of which were down-regulated in pluripotent calli. Notably, a majority of antioxidant genes were up-regulated in pluripotent calli ( Fig. 4i), as supported by KEGG analysis. In addition, several genes involved in developmental processes and stress responses were also changed in pluripotent calli (Fig. 4j,k), indicating that plant metabolic and developmental processes are massively reprogrammed during pluripotency establishment.
Callus formation and pluripotency acquisition are strongly associated with hormone metabolism and signaling 24 . To further inspect the relevance of hormone signaling in pluripotency acquisition in rice, a list of genes involved in hormone metabolism were investigated in DEGs. Although auxin and cytokinin signaling are known to play crucial functions in plant regeneration 3,21,24 , the genes related to auxin and cytokinin metabolism were negligibly influenced in pluripotent rice calli (see below). Instead, metabolic genes of gibberellic acid (GA) and abscisic acid (ABA) were particularly changed (Fig. 4l,m), suggesting that pluripotency acquisition in a limited cell population of callus might be conferred by intimate interactions between GA and ABA signaling.
Both pluripotent and non-pluripotent calli were generated by incubation on auxin-rich CIM, and they showed a comparable cell proliferation activity (Fig. 1). Since we compared both proliferating calli, genes involved in cell proliferation and cell cycle progression were negligibly influenced in pluripotent calli compared to nonpluripotent calli. It may also explain the observation that auxin and cytokinin signaling, which primarily affect cell division and proliferation, were not significantly changed. Since a limited cell population have distinguishable cellular features between pluripotent and non-pluripotent calli, the number of DEGs was small, and biological processes typically required for plant regeneration [25][26][27][28] were limitedly enriched in the DEGs. Nonetheless, our analysis is still relevant in that the results enable to gain a new insight into the pluripotency acquisition in rice callus.
Antioxidant activity is increased in pluripotent calli. Our RNA-seq analysis showed that many metabolic processes are regulated during callus formation, among which increased antioxidant activity in pluripotent calli was repeatedly suggested (Figs. 3a, 4i). In support, several previous studies have also shown that pluripotent stem cells contain higher levels of antioxidants in the maintenance of stemness 25,[29][30][31][32] . We thus put our focus more on the relevance of antioxidant activity in pluripotent acquisition in rice callus.
We asked whether antioxidant levels are indeed increased in pluripotent calli. To this end, we measured levels of polyphenols, which includes several types of antioxidants such as phenolic acid, hydrolysable tannins and flavonoids with robust ROS-scavenging activity 33 , in pluripotent and non-pluripotent callus extracts. Polyphenol measurement revealed that total polyphenol levels of pluripotent calli were 1.2 times higher than those of nonpluripotent calli independently of callus incubation period (Table 2). www.nature.com/scientificreports/ The radical scavenging activity was further determined by the 2,2-diphenyl-1-picrylhydrazyl (DPPH) test. As a result, antioxidant activity was distinguishable between pluripotent and non-pluripotent calli. Pluripotent callus extract (EC 50 = 1.1) had 12 times greater antioxidant activity compared to extracts from non-pluripotent (EC 50 = 0.09) at 8 weeks after callus incubation (Table 3). These results indicate that enhanced antioxidant activity is critical for pluripotency maintenance in callus.

Network analysis identifies novel genes involved in pluripotency establishment in callus.
Our dataset indicates that 244 genes were differentially expressed between pluripotent and non-pluripotent rice calli. The robustness of our transcriptome prompts us to query the expanded gene set for pluripotency acquisition. To this end, we decided to perform gene regulatory network (GRN) analysis to identify more candidate genes, whose expression is linked to core regulators of pluripotency acquisition. Based on a publicly-available GRN (RiceFREND; http:// https ://ricef rend.dna.affrc .go.jp/), we first extracted the sub-network using QUISCENT CENTER SPECIFIC HOMEOBOX (QHB/OsWOX5) as the seed gene, which is known as a core regulator of pluripotency establishment 16,34,35 . In the subnetwork, PLT1 and PLT4, the key players of plant regeneration, co-expressed with QHB ( Supplementary Fig. S3). We also investigated the DEG genes in this sub-network containing QHB and PLTs and found that LOC_Os05g35070 was included as a potential regulator of pluripotent acquisition ( Supplementary Fig. S3).
To obtain additional candidates regulating cellular pluripotency, we reciprocally extracted subnetworks using our DEGs. All DEGs in pluripotent callus were used as queries and a sub-network for each DEG was extracted. Among the input DEGs, only three genes, LOC_Os04g32300, LOC_Os06g08340, and LOC_Os03g42420, formed GRNs coexpressed with QHB and/or PLTs (Fig. 5a). The three subnetworks obtained could be merged ( Fig. 5a; Supplementary Table S5), constructing a putative GRN related to pluripotency acquisition. Interestingly, the GRN basically contained core regulators of plant regeneration, such as PLT1, PLT4, QHB, WOX9, ESR1, and CUC3, suggesting that at least the three DEGs are potentially important for pluripotency acquisition in rice callus. From the putative pluripotency GRN, the top 10 hub genes were predicted based on the lowest average path length and higher connectivity (Fig. 5b): CHR38, PLT1, ANT, CUC3, WOX9, HB33, OSH6, OBP2, PLT4, and RFL/LFY. Among them, the potential role of CHR38, ANT, HB33, OSH6, and RFL/LFY have not yet been investigated for plant regeneration to date, and thus they can be considered as novel potential regulators of pluripotency acquisition.
To validate our prediction, we chose a candidate gene and analyzed the potential role in plant regeneration using Arabidopsis model system. The OSH6 gene was of particular interest, because it is expressed at a confined pluripotent region of rice embryo, where shoot apical meristem (SAM) will be specified 36,37 . In addition, OSH6 is also highly expressed in callus tissues (Supplementary Figs. S1 and S5). RT-qPCR analysis enabled to measure sensitive gene expression changes, compared with RNA-seq, and revealed that transcript accumulation of OSH6 was significantly increased in pluripotent calli relative to non-pluripotent calli (Supplementary Fig. S4). Hence, we retrieved the Arabidopsis homolog of rice OSH6 gene, KNOTTED1-LIKE HOMEOBOX GENE 6 (AtKNAT6, At1g23380), based on phylogenetic analysis (Supplementary Fig. S6). To convince the AtKNAT6 gene is the functional homolog of rice OSH6 gene, we compared tissue-specific expression of Class I KNOTTED1-LIKE Table 2. Polyphenol content in pluripotent and non-pluripotent calli. Total content of polyphenolic compounds was determined according to Prussian Blue method (Hagerman and Butler 55 ). Quercetin was used as a standard. Three biological replicates were averaged, and the data are represented as the mean ± standard deviation. DW, dry weight.

Callus incubation period on CIM (days)
Pluripotent callus (mg/g DW) Non-pluripotent callus (mg/g DW) 28 48 . Fig. S6) and confirmed that only AtKNAT6, but not KNAT2, is highly expressed in callus tissues similar to OSH6 ( Supplementary Fig. S5), suggesting that the AtKNAT6 gene is a conserved homolog of rice OSH6 gene possibly in the control of pluripotency acquisition during callus formation. We next obtained a knat6 loss-of-function mutant (GK-478F03) and conducted in vitro tissue culture. The knat6 mutant exhibited lower shoot regeneration efficiency with the reduced number of regenerated shoots ( Fig. 5c and 5d), which was attributable to reduced pluripotency acquisition. Overall, we identified potential candidates implicated in pluripotency acquisition in rice callus, which are most likely associated with the stem cell specification and maintenance.

Discussion
In vitro plant regeneration is driven by the exogenous application of phytohormones, auxin and cytokinin, and involves a two-step process: callus formation and de novo shoot regeneration 3 . In Arabidopsis, the hormoneinduced callus is analogous to lateral root primordium [7][8][9] . The callus tissue further establishes root stem cell niche, through the expression of root stem cell regulators, including PLTs and WOXs 9,11,[14][15][16]35 . Since PLT and WOX-induced root stem cell identity is a cellular nature of pluripotency [14][15][16] , they facilitate to form competence for tissue regeneration in callus. In support, only after expression of the pluripotency factors in callus, PLTs and WOXs, de novo shoot organogenesis can be initiated on SIM 11,12 . In this study, we compared the transcriptome of pluripotent and non-pluripotent calli to understand molecular networks involved in pluripotent acquisition in rice. Notably, molecular processes underlying pluripotent callus formation in rice are likely similar to those of dicot plants. PLTs were up-regulated in pluripotent callus relative to non-pluripotent callus in rice. Additionally, based on a putative GRN related to cellular pluripotency in rice (Fig. 5a), WOXs, CUC s and ESR1, which are known to be associated with cellular pluripotency in Arabidopsis 9,11,16,35,[38][39][40] , were also involved in pluripotency acquisition in rice callus, indicating the similarity of hormone-induced callus of rice and Arabidopsis. Together, PLTs and WOXs are major determinants of cellular pluripotency across diverse plant species and can be used as key marker genes for pluripotency acquisition. Several metabolic processes are associated with pluripotency establishment in rice callus. Photosynthetic processes significantly suppressed in pluripotent calli. Since callus resembles root primordium 7-9 , shoot developmental processes are essentially impaired in callus 20 . Consequently, primary metabolism was also globally influenced. In particular, sugar, amino acid, and nitrogen metabolism were suppressed during callus formation and pluripotency acquisition. We suspect that carbon and energy metabolism are substantially reprogrammed to optimize cellular pluripotency.
Notably, genes related to fatty acid metabolism and ROS homeostasis are enriched in upregulated DEGs of pluripotent calli. Several lipid species, such as prostaglandin E2, linoleic acid, and albumin-associated lipids, are demonstrated as key metabolites for pluripotency establishment and maintenance in mammalian stem cells 41 . In addition, ROS acts either signaling molecules or detrimental agents inducing cell death 42,43 . ROS levels are delicately balanced in stem cell niche in plants, which is responsible for pluripotency acquisition 44,45 . Consistently, our study also showed that increased antioxidant activity in pluripotent calli is important for acquiring competence for de novo organogenesis. Overall, acquisition of cellular pluripotency is an active process that accompanies reprogramming of primary and secondary metabolism.
Since the pluripotent and non-pluripotent calli used in our study were incubated for same period of time on CIM, only a limited cell population distinguishes callus samples with different capacity of cellular pluripotency. Consistently, a small number of DEGs were obtained, and several key factors involved in auxin and cytokinin signaling and cell proliferation were unidentified from the DEGs. To understand biological relevance of 244 DEGs of pluripotent callus, coexpression network analysis was employed. Using a public database available for extracting coexpression network in rice, we could propose a putative GRN involved in cellular pluripotency (Fig. 5a). Several DEGs formed a GRN with core regulators of plant regeneration, such as PLTs, WOXs, and CUC s, and interestingly, they could be integrated into a larger coexpression network. Based on the results, several potential regulators can be suggested. The first group includes the DEGs that construct a gene networks with PLTs and WOXs: LOC_Os03g42420, LOC_Os04g32300, LOC_Os05g35070, and LOC_Os06g08340. The second group includes genes, which constitute a hub with high connectivity in the putative pluripotency GRN, but have not studied to date: CHR38, ANT, HB33, OSH6, and RFL. To validate our suggestion, we chose a candidate, OSH6, and confirmed its role in de novo shoot regeneration (Fig. 5c,d). The rice OSH6 gene possibly regulates both cytokinin production and auxin responses [46][47][48] , linking cytokinin and auxin signaling in the control of pluripotency acquisition, although the detailed molecular mechanism should be elucidated in the future. Further, the remaining 8 potential candidates would also be associated with pluripotency acquisition and should be studied in the future to gain a comprehensive view of cellular pluripotency in rice.
Altogether, we propose a hypothesis how plant callus acquires pluripotency in rice. Based on the DEGs in pluripotent callus, we proposed molecular processes as well as putative GRNs involved in pluripotency acquisition. Our study provides candidates for evaluating the involvement of genes in pluripotency acquisition. It also provides a resource for comparative transcriptome analysis of plant regeneration in other species.
RNA isolation and cDNA synthesis. Approximately 0.1 g (fresh weight) of pluripotent and non-pluripotent calli at 28 DAC were harvested respectively, and immediately frozen in liquid nitrogen for RNA isolation. Total RNAs were isolated using the RNeasy plant mini kit (Qiagen). Genomic DNA was eliminated by DNAse I (Invitrogen, Carlsbad, CA, USA) treatment, as recommended by the manufacturer. The purity of RNA samples was examined by use of the Agilent 2100 Bioanalyzer (Agilent Technologies). cDNA was synthesized from Primer design and RT-qPCR conditions. Gene-specific primers for qRT-PCR analysis were designed using the primer 3.0 online tool (http://bioin fo.ut.ee/prime r3/) according to the sequences of reference genes and a target gene. Primers were synthesized by the Bioneer Genomics Institute (Daejeon, Korea) with the following parameters: T m values ranging from 50 to 65 °C, GC percent of 45-55%, primer lengths of 17-25 bp, and product lengths of 50-150 bp. RT-qPCR was carried out with the SYBR Green PCR Master Mix system (Bioneer, Daejeon, Korea) on an Applied Biosystems 7500/7500 Fast Real-time PCR System (ABI, CA, USA). The PCR amplification system and program were performed as described previously 49 . Three biological replicates were analyzed. The RefFinder online software was applied to further estimate the stability of reference genes. Relative gene expression levels were calculated using the 2 −ΔΔCt method 50 . All data were recorded as mean ± standard deviation and statistical significance was determined by SPSS (version 9.0 for Windows 98, SPSS Inc.). Primers used in RT-qPCR analysis were listed in Supplementary Table S6.
Illumina sequencing and data analysis. For RNA-seq library construction, 1 μg of total RNA extracted from pluripotent calli and non-pluripotent calli incubated on CIM for 28 days. Since callus samples are heterogeneous and have significant variation in gene expression, we collected large amounts of samples (> 300 calli) for single biological replicate. Because of the sampling burden, we performed biological duplicates. The total RNA was purified by DNase and mRNA purification kit. Then the purified mRNA samples fragmented in 94 °C for 8 min. Fragmented mRNAs were converted into cDNA using random primer. The cDNA samples were ligated by adapters using TruSeq RNA kit (Illumina, CA, USA) and amplified by PCR. Libraries were qualified using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). To conduct paired end sequencing, libraries were subjected to Illumina HiSeq2000 platform (Illumina, CA, USA) for two biological replicates.
To identify high confidential DEGs, we mapped raw RNA-seq reads and quantified transcript abundance considering distribution of mate-pair inner distance. To this end, briefly, we preliminarily mapped the raw read sequences against primary reference transcriptome sequences of rice (MSU v.7.0; https ://phyto zome.jgi.doe. gov/pz/porta l.html) using Bowtie2 aligner 51 . The mean and standard deviation of mate-pair inner distance were estimated for each sequencing library (Supplementary Table S7). Then, the raw reads were re-aligned against the rice reference genome (MSU v.7.0) considering mean and standard deviation of mate inner distance, and splice junction using Tophat2 and Bowtie2 51,52 . Transcriptome abundance was calculated using Cufflinks and Cuffquant along with application of -frag-bias-correct and -multi-read-correct options for precise quantification 53 . Statistical significance and fold change of each transcriptome was analyzed using Cuffdiff 53 . The differential expression analysis of two samples was performed using criteria, including the absolute value of log 2 fold change ≥ 1 and P < 0.05, to ensure the significance of gene expression difference.
KEGG pathway enrichment analysis. To conduct KEGG pathway enrichment analysis, the number of genes in a specific pathway was counted for gene group of interest and whole reference genome, respectively 54 . The statistical significance was calculated based on the binomial test comparing observed ratio of genes for specific pathway to expected ratio. Antioxidant activity analysis. All samples were homogenized with a mortar, dissolved in three volumes of methanol, and centrifuged at 12,000g for 10 min at room temperature. The supernatant was stored at − 80 °C until use for the analysis. The total content of polyphenolic compounds was measured according to the Prussian Blue method 55 with several modifications. As a standard, Quercetin (Sigma Chemical Co., St. Louis, MO) was used.
Radical scavenging activity was measured according to a spectrophotometric method using an ethanol solution of DPPH 56 . As positive controls, Trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid), BHT (butylated hydroxytoluene), and ascorbic acid were used. The EC 50 values, defined as the amount of antioxidant necessary to decrease the initial DPPH concentration by 50%, were calculated. Three biological replicates were averaged, and statistical significance was determined by SPSS (version 9.0 for Windows 98, SPSS Inc.).

Analysis of GRN for DEGs.
To gain more insights into the acquisition of pluripotency from the DEGs, we searched GRNs carrying DEGs using RiceFREND database containing co-expression gene network which was constructed based on gene expression in various developmental stage and phytohormone treatments in rice 57 . All DEGs were queried in the database. The results were manually inspected whether the DEGs form a subnetwork with key pluripotency regulators, such as QHB, PLTs and CUC3. The sub-networks link DEG with key molecular component of pluripotency were downloaded from the database and merged and visualized using Cytoscape 58 . Degree and average shortest path length of genes in merged network were analyzed using Networ-kAnalyzer embedded in Cytoscape 59 .
Phylogenetic and tissue-specific expression pattern analyses. To identify Arabidopsis orthologous gene of rice OSH6 gene, we constructed phylogenetic tree using all KNOTTED1-LIKE HOMEOBOX (KNOX) genes in the Arabidopsis and rice genomes. The KNOX protein sequences from primary transcripts were retrieved and aligned using MUSCLE algorithm with default parameters 60 . An unrooted phylogenetic tree was constructed using MEGA X software based on the neighbor-joining algorithm with complete deletion for gaps