REGENOMICS: A web-based application for plant REGENeration-associated transcriptOMICS analyses.

In plants, differentiated somatic cells exhibit an exceptional ability to regenerate new tissues, organs, or whole plants. Recent studies have unveiled core genetic components and pathways underlying cellular reprogramming and de novo tissue regeneration in plants. Although high-throughput analyses have led to key discoveries in plant regeneration, a comprehensive organization of large-scale data is needed to further enhance our understanding of plant regeneration. Here, we collected all currently available transcriptome datasets related to wounding responses, callus formation, de novo organogenesis, somatic embryogenesis, and protoplast regeneration to construct REGENOMICS, a web-based application for plant REGENeration-associated transcriptOMICS analyses. REGENOMICS supports single- and multi-query analyses of plant regeneration-related gene-expression dynamics, co-expression networks, gene-regulatory networks, and single-cell expression profiles. Furthermore, it enables user-friendly transcriptome-level analysis of REGENOMICS-deposited and user-submitted RNA-seq datasets. Overall, we demonstrate that REGENOMICS can serve as a key hub of plant regeneration transcriptome analysis and greatly enhance our understanding on gene-expression networks, new molecular interactions, and the crosstalk between genetic pathways underlying each mode of plant regeneration. The REGENOMICS web-based application is available at http://plantregeneration.snu.ac.kr.

Introduction

Being sessile organisms, plants are vulnerable to a plethora of environmental stimuli, some of which can cause wounds or the complete loss of organs. To recover from environmental damage, plant cells acquire pluripotency and give rise to new undifferentiated stem cells [49]. The stem cell niche exhibits continuous cell division and subsequent cell differentiation, leading to the formation of various cell types and replenishing the damaged tissues [69], [41]. In addition, owing to the remarkable tissue regeneration potential of plants, several in vitro tissue-culture techniques have been developed [60]. Exogenous hormone treatment with a certain auxin to cytokinin ratio allows pluripotency acquisition and tissue regeneration [60].

Four main modes of plant regeneration have been investigated to date: wound-induced tissue repair, hormone-induced plant regeneration, somatic embryogenesis, and protoplast regeneration [51], [9], [1], [57]. Under natural conditions, wounding causes rapid cell proliferation at injury site to replenish damaged cells and tissues [38]. Upon the excision of Arabidopsis root tips as a model of wound-induced tissue repair, jasmonic acid (JA) and auxin accumulation leads to the activation of CYCLIN D6 (CYCD6;1) and ETHYLENE RESPONSE FACTOR 115 (ERF115), which in turn activates the RETINOBLASTOMA-RELATED (RBR)-SCARECROW (SCR)-SHORTROOT (SHR) module, regenerating the stem cell niche and subsequently the root tip [57], [85]. In some cases, wounding leads to the formation of a dedifferentiated cell mass, called callus [24]. Callus cells are pluripotent and therefore are capable of regenerating tissues, organs, or even plantlets [1]. Wound-responsive WOUND INDUCED DEDIFFERENTIATION (WIND) transcription factors are key players in wound-induced tissue repair and callus development, especially at the injury site [24]. WIND-induced callus formation is dependent on ENHANCER OF SHOOT REGENERATION 1 (ESR1)/DORNRÖSCHEN (DRN) and type-B ARABIDOPSIS RESPONSE REGULATORs (ARRs), which activate cytokinin signaling-associated cell proliferation and tissue regeneration [22].

Plant regeneration can be facilitated through a two-step in vitro tissue-culture method [26], [83], [34]. Pluripotent callus can be obtained from differentiated plant tissues by culturing explants on auxin-rich callus-inducing medium (CIM) [1], [63]. CIM-induced callus formation resembles the lateral root initiation process [1], [63]. Consistently, ABERRANT LATERAL ROOT FORMATION 4 (ALF4), WUSCHEL-RELATED HOMEOBOX 11 (WOX11), and LATERAL ORGAN BOUNDARIES DOMAIN 16/17/18/29 (LBD16/17/18/29) genes, which play crucial roles in lateral root initiation, are also essential for callus formation [63], [35]. In the middle of the callus formation process, pluripotency acquisition is enabled with the activation of root meristem regulators such as PLETHORA 1 (PLT1), PLT2, SCR, and WOX5 [1], [63], [27]. Following pluripotency acquisition, calli are usually transferred to cytokinin-rich shoot-inducing medium (SIM), which promotes de novo shoot organogenesis [26], [34]. The cytokinin-inducible type-B ARR transcription factors directly bind to the WUSCHEL (WUS) promoter, specify the shoot stem cell niche in callus, and promote de novo shoot regeneration [83].

Somatic embryogenesis can be induced by the application of high concentrations of auxin to somatic cells or somatic cell-derived callus [75], [9]. Somatic embryogenesis is stochastically initiated during in vitro tissue culture, and an emerged somatic embryo gives rise to a new plant [9]. Because somatic embryogenesis is similar to normal embryogenesis, both processes share common genetic components [14], including BABY BOOM (BBM), LEAFY COTYLEDON 1 (LEC1), LEC2, ABSCISIC ACID INSENSITIVE 3 (ABI3), FUSCA 3 (FUS3), and AGAMOUS-LIKE 15 (AGL15) [14], [71].

Notably, a single protoplast, which is a cell dissociated from its cell wall, has the ability to proliferate to form microcallus and subsequently regenerate into a whole plant [64]. This mode of plant regeneration involves dedifferentiation from a fully differentiated somatic cell, which is essential for callus formation and tissue regeneration. Protoplast regeneration shares some features of other types of plant regeneration, as crucial molecular factors, such as WINDs, WUS and ESR1/DRN [5], [76], are pervasively involved, although much is unknown about the regulatory mechanisms underlying protoplast regeneration.

Recent studies on plant regeneration have elucidated key genetic players and signaling pathways underlying pluripotency acquisition and plant regeneration, and have generated enormous amounts of high-throughput RNA-seq data relevant to plant regeneration processes [73], [33], [17], [15], [27], [20], [47], [52], [81], [6], [37], [40], [59], [71], [78], [30], [54], [76], [13], [19], [29]. However, despite the numerous RNA-seq datasets available publicly, extracting meaningful information from these data remains challenging. A few web applications, such as ePlant [70] and ATTED [44], allow users to easily explore gene-expression patterns and co-expression networks. However, none of these web applications are solely dedicated to plant regeneration studies or allow customized transcriptome-level analyses. Furthermore, because distinct genetic pathways underlie different modes of plant regeneration, it is important to analyze the gene expression patterns and co-expression networks in the context of each plant regeneration mode. In addition, although studies on plant regeneration have been mainly conducted in Arabidopsis thaliana to date due to its advantages as a model species and a high capacity of cellular reprogramming [18], [39], plant regeneration is an important issue in other plant species [10], [21], [43], [19]. Therefore, efforts to integrate transcriptome data generated from various plant species are required to investigate the universal or specific regulators of plant regeneration across the plant kingdom, highlighting the necessity of a transcriptome analytic hub for regeneration studies on various plant species.

A large volume of single-cell RNA-seq data has also been generated recently from a diverse range of plant species and tissues [72], [58], [76], [80], [82], [36], [45]. Accordingly, tools for single-cell data analysis are being newly developed and improved at an astonishing speed, allowing us to interpret important biological processes at a single-cell level [53], [79]. One of the major developments is trajectory inference, which uses single-cell datasets to computationally predict developmental trajectories within the dataset and order cells according to a pseudotime, an abstract measure of biological progress [66]. High-throughput single-cell datasets also allow the construction of high-confidence gene regulatory networks (GRNs). As single-cell datasets provide massive amounts of transcriptome data with little technical bias, sophisticated gene regulations can be predicted with higher reliability. However, GRNs inferred with plant single-cell RNA-seq data are currently unavailable [67]. Another limitation is that many recently developed tools for single-cell data analysis are not yet readily accessible in the form of web applications.

Here, we developed a web-based database and analytic platform, REGENOMICS, which integrates all currently available transcriptome datasets related to the four main modes of plant regeneration, and enables integrative analysis to gain new insights into the crucial connections between various genetic components and signaling pathways underlying the plant regeneration process. REGENOMICS provides a user-friendly interface and is freely available at http://plantregeneration.snu.ac.kr.

Results and discussion

Overview of REGENOMICS

REGENOMICS is a plant regeneration-dedicated database that contains all currently available transcriptomic data related to regeneration, specifically wound-induced tissue repair, hormone-induced plant regeneration, somatic embryogenesis, and protoplast regeneration (Table 1), along with an analytical pipeline that processes user-submitted queries. The REGENOMICS analytical pipeline is divided into three analysis modes (Single-gene analysis, Multi-gene analysis, and Transcriptome analysis) depending on the type of query (a single gene of interest [GOI], a list of GOIs, or a whole transcriptome) (Fig. 1). In the Single-gene analysis mode, a user can view the expression pattern, co-expression network constructed using bulk RNA-seq, and single-cell RNA-seq-derived GRNs of any GOI involved in plant regeneration. The Multi-gene analysis mode supports the comparison of expression patterns, construction of co-expression networks and GRNs, and differentially expressed gene (DEG) enrichment analysis for a list of GOIs. The Transcriptome analysis mode enables users to select REGENOMICS-deposited and user-submitted transcriptomes for correlation analysis, dimensionality reduction analysis, and differential expression analysis.

Table 1

Datasets Collected and Processed for the REGENOMICS Database.

Category	Tissue (origin)	Plant Materials	Accession Number	Reference
Wound-induced tissue repair	Hypocotyl	Arabidopsis thalianaCol-0	GSE101422	[17]
Wound-induced tissue repair	Leaf	Arabidopsis thalianaCol-0, coi1-2, sdg8-2	GSE120418	[81]
Wound-induced tissue repair	Leaf	Arabidopsis thalianaCol-0, pSPL10:rSPL10-GR	PRJCA001184	[78]
Wound-induced tissue repair	Leaf	Arabidopsis thalianaCol-0	GSE108253	[47]
Wound-induced tissue repair	Root	Arabidopsis thalianaCol-0	GSE168385	[13]
Wound-induced tissue repair	Root	Arabidopsis thalianaCol-0	GSE145488	[40]
Wound-induced tissue repair	Root	Arabidopsis thalianaCol-0	E-MTAB-7609	[52]
Hormone-induced plant regeneration	Hypocotyl	Arabidopsis thalianaCol-0, siz1-2	GSE141188	[6]
Hormone-induced plant regeneration	Hypocotyl	Arabidopsis thalianaCol-0, plt1plt2	GSE178354	[80]
Hormone-induced plant regeneration	Hypocotyl	Arabidopsis thalianaCol-0, wox5-1wox7-1	GSE156990	[80]
Hormone-induced plant regeneration	Hypocotyl	Arabidopsis thalianaCol-0	GSE176161	[29]
Hormone-induced plant regeneration	Root	Arabidopsis thalianaCol-0, ldl3-1	DRP004816	[20]
Hormone-induced plant regeneration	Root	Arabidopsis thalianaCol-0, hag1-6	GSE100966	[27]
Hormone-induced plant regeneration	Root	Arabidopsis thalianaCol-0, atprmt5	E-MTAB-5044	[33]
Hormone-induced plant regeneration	Root	Arabidopsis thalianaCol-0, arr1-4, 35Spro:ARR1	GSE146690	[37]
Protoplast regeneration	Mesophyll cells	Arabidopsis thalianaCol-0	PRJNA739879	[54]
Protoplast regeneration	Mesophyll cells	Arabidopsis thalianaCol-0	PRJNA648028	[76]
Somatic embryogenesis	Immature embryo and seedling	Arabidopsis thalianaCol-0, p35S::LEC2-GR	PRJCA002620	[71]
Somatic embryogenesis	Immature embryo	Arabidopsis thalianaCol-0	E-MTAB-2403	[73]
Root tip cell types and developmental zones	Root tip	Arabidopsis thalianaCol-0	PRJNA323955	[32]
Wound-induced tissue repair	Hypocotyl	Solanum lycopersicumMicro-Tom	PRJNA731333	[30]
Somatic embryogenesis	Seed	Oryza sativaTNG67, IR64	GSE82138	[15]
Somatic embryogenesis	Seed	Oryza sativassp. japonica	PRJNA804674	[59]
Wound-induced tissue repair	Thallus	Marchantia polymorphaTak-1	PRJDB12610	[19]
Wound-induced tissue repair	Thallus	Marchantia polymorphaTak-1, Mperf15ko	GSE196912	Unpublished
Single-cell RNA-sequencing	Hypocotyl explant	Arabidopsis thalianaCol-0	GSE156991	[80]
Single-cell RNA-sequencing	Leaf explant	Arabidopsis thalianaCol-0	GSE147289	[36]
Single-cell RNA-sequencing	Shoot-borne roots	Solanum lycopersicumcv M82	GSE159055	[45]
Single-cell RNA-sequencing	Root tip	Arabidopsis thalianaCol-0	GSE141730	[72]
Single-cell RNA-sequencing	Shoot apex and leaf	Arabidopsis thalianaCol-0	PRJCA003094	[82]
Single-cell RNA-sequencing	Protoplasts	Arabidopsis thalianaCol-0	PRJNA648028	[76]
Single-cell RNA-sequencing	Lateral root primordium	Arabidopsis thalianaCol-0	GSE161970	[58]

Fig. 1

Workflow for REGENOMICS as shown in the REGENOMICS web interface. Transcriptomic datasets related to plant regeneration are collected and processed according to the REGENOMICS analytical pipeline. The pipeline is divided into three sections, Single-gene analysis, Multi-gene analysis, and Transcriptome analysis, each of which accepts a gene of interest (GOI), a list of GOIs, and transcriptomes as the input query, respectively. Then, REGENOMICS outputs the gene expression profile(s), co-expression network and gene-regulatory network for input gene(s). Transcriptome analysis supports correlation, dimensionality reduction, and differential expression analyses for user-defined transcriptome datasets.

The REGENOMICS application has five main features. First, gene expression profiles of any GOI(s) involved in any mode of plant regeneration can be easily accessed by searching for either the gene name or the gene locus code. The mode of plant regeneration can be selected from the dropdown menu. Second, REGENOMICS contains five different co-expression networks: one built with all REGENOMICS-deposited bulk RNA-seq datasets, and four constructed using transcriptome datasets corresponding to wound-induced tissue repair, hormone-induced plant regeneration, somatic embryogenesis, or protoplast regeneration. A sub-network can be drawn with the submitted query gene(s) to visualize its (their) connections with closely co-expressed genes in the context of each mode of plant regeneration. Gene ontology (GO) and DEG enrichment analyses are performed using genes comprising the resulting output sub-network to estimate biological features of the sub-network. Third, DEG enrichment analysis allows the estimation of the biological and functional aspects of a group of genes. This analysis can be performed with genes comprising the co-expression network as well as a user-submitted list of GOIs. Fourth, the REGENOMICS database contains single-cell transcriptome data of tissues undergoing regeneration as well as root and shoot apical meristems. Users can explore cell type-specific expression patterns of GOI(s) and also perform trajectory inference analysis for each dataset. Furthermore, genes that regulate and are regulated by the GOI(s) can be viewed from the sub-networks drawn from the GRNs inferred from each single-cell RNA-seq dataset. Lastly, users can input any Arabidopsis transcriptome profile (RNA-seq data) in comma-separated values (CSV) format to analyze its transcriptome-level correlation using the selected REGENOMICS-deposited transcriptome data and user-deposited RNA-seq data. In addition, users can select two transcriptome datasets and conduct differential expression analysis to obtain a table of DEGs as well as the corresponding GO terms.

It should be noted that REGENOMICS is not limited to Arabidopsis data, but includes data from other species. Users can select the plant species in the dropdown menu located at the top of each page, and perform the same analysis for each species. However, while all analyses shown here are available for use with Arabidopsis studies, some of the analyses provided in REGENOMICS are not yet ready for use for other species, due to the limited availability of RNA-seq data. With future researches, we plan to support all features of REGENOMICS for all plant species to promote regeneration studies covering a wide spectrum of plant species. Overall, REGENOMICS is designed to enable easy access to plant regeneration-related high-throughput transcriptome data to gain a comprehensive insight into the molecular connections between the genetic components and signaling pathways underlying each mode of plant regeneration in various species (Fig. 1).

Single-gene analysis

To demonstrate each analysis supported in REGENOMICS, we used Arabidopsis studies as representative examples. While some features are currently available only for Arabidopsis due to a lack of published data for other species (Fig. 1), all analyses in REGENOMICS will become available for every species in the same format as Arabidopsis through future updates. In the Single-gene analysis mode, the REGENOMICS web application supports three different analyses (Gene expression browser, Co-expression network, and Single-cell analysis) using all currently available transcriptome data related to wound-induced tissue repair, hormone-induced plant regeneration, somatic embryogenesis, and protoplast regeneration. Users can retrieve information for any gene by searching for a gene name (e.g., WOX5) or the gene locus code (e.g., AT3G11260) (Fig. 2A). In the Gene expression browser, the GOI’s expression values normalized in transcripts per million (TPM) are displayed as line plots or bar plots for each study under the selected mode of plant regeneration (Fig. 2B). To validate the gene expression browser of REGENOMICS, we queried several genes involved in wound-induced tissue regeneration, hormone-induced plant regeneration, somatic embryogenesis, or protoplast regeneration. The results obtained using REGENOMICS showed that the wound-responsive ERF109 was transiently up-regulated shortly after wounding (Fig. 2B) as previously reported [81]. Additionally, the expression patterns of LBD16, ESR1/DRN, and FUS3 in REGENOMICS, which were reported to be essential in hormone-induced regeneration, protoplast regeneration, and somatic embryogenesis, respectively, matched the patterns reported previously (Fig. 2B) [14], [34], [71], [76].

Fig. 2

Interface and features of Single-gene analysis. (A) Interface and search page for Single-gene analysis. (B) Gene expression browser. Gene expression browser returns gene expression patterns in line plots for a single query gene. Expression patterns for ERF109, LBD16, ESR1/DRN, and FUS3 are shown in different modes of plant regeneration. (C) Co-expression network analysis. The PLT1 gene (node circled in red) was used as an input query gene. Co-expressed genes include key regulators of the root meristem such as PLT2, BBM, WOX5, RGF8, and RGFR2. Genes with the GO term ‘anatomical structure development’ are colored yellow-green. The network was built with hormone-induced plant regeneration bulk RNA-seq datasets. (D) Uniform Manifold Approximation and Projection (UMAP) for Single-cell analysis. UMAP plot colored according to the expression of WOX5 in single-cell transcriptome of hypocotyl callus is shown. Each dot represents a single cell (see also Fig. S3). (E) Trajectory inference analysis for Single-cell analysis. UMAP plot shows the trajectory inference analysis results performed with the single-cell transcriptome of hypocotyl callus. Each cell is colored according to the calculated pseudotime. (F) Gene regulatory network for Single-cell analysis. The sub-network was drawn using WOX5 gene (node circled in red) as an input query gene. Key regulators of the pluripotency acquisition, such as PLT1 and SCR, were included in the sub-network. Genes with the GO term ‘anatomical structure development’ are colored yellow-green. The network was inferred from single-cell transcriptome of hypocotyl callus. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

REGENOMICS generates a co-expression network that connects functionally-related genes based on correlated expression patterns. Co-expression networks were built by collecting the bulk RNA-seq data of each mode of plant regeneration and then analyzing gene-to-gene expression correlation via the Mutual Rank (MR) [44]. Users can access a sub-network constructed with the user-submitted GOI and its close neighbor genes that have significantly correlated expression levels (Fig. 2C). By selecting a mode of plant regeneration, users can draw a sub-network built with datasets specific for that mode of plant regeneration. Users can interact with the output network by dragging the nodes, which represent the query gene and its co-expressed genes. Edges between nodes are colored, according to the MR gradient, and can be filtered out by changing the MR threshold with the slider (Fig. 2C). A downloadable table including co-expressed genes and their MR values is provided below the sub-network (Table S1). Additionally, the biological processes overrepresented in the output sub-network can be viewed by selecting a GO term, which colors the gene nodes corresponding to the GO term (Fig. 2C). These sub-network genes are also used for DEG enrichment analysis. We generated sets of up- or down-regulated genes (DEGs) from pair-wise comparisons of transcriptomes for each study deposited in the REGENOMICS database. Enrichment analysis allows the identification of DEG sets enriched in the output sub-network (Fig. S1), thus providing an insight into the potential biological impact of the sub-network on plant regeneration processes.

We conducted a case study using PLT1, a gene known to be related to the pluripotency acquisition of callus cells [26], as query. Co-expression network analysis using hormone-induced plant regeneration datasets (Fig. 2C) identified PLT2, another known regulator in establishing pluripotency during callus formation [26]. Similarly, root-meristem regulator genes, such as WOX5 and RGF8 [16], [27], [46], were also shown to be closely co-expressed with PLT1 (Table S1). Consistently, genes with a GO term, ‘anatomical structure development’, were predominantly found in the output sub-network (yellow-green color in Fig. 2C). Furthermore, DEG enrichment analysis showed that the sub-network was enriched with genes up-regulated in the hypocotyl explant-derived calli incubated on CIM (DACs in Fig. S1), relative to samples from wounding or de novo shoot regeneration [6], consistent with the role of PLT1 in callus formation [26].

Moreover, a query gene can be analyzed by the Single-cell analysis, which contains single-cell RNA-seq datasets for different modes of plant regeneration as well as meristem (and meristem-like) tissues (Table 1). Users can view the gene expression dynamics of the query gene occurring during the course of plant regeneration as well as the expression patterns in various cell types found in pluripotent cells and stem-cell derived differentiated cells. We analyzed single-cell RNA-seq datasets according to our unified analysis pipeline and annotated cell types using cell type-specific marker genes. Transcript accumulation of the GOI in each cell type was visualized as a violin plot (Fig. S2) and Uniform Manifold Approximation and Projection (UMAP) plots for each cell (Fig. 2D). To validate the data processing and cell-type annotation, we used PLT1 as query in the single-cell RNA-seq dataset from hypocotyl-derived callus, and confirmed that the expression of each gene was enriched in QC-like-cells (Figs. S2 and 2D) as reported in the original study [80]. Additionally, trajectory inference analysis allows the understanding of the cell identity transition dynamics in terms of pseudotime (Fig. 2E). We were able to confirm that QC-like cells differentiate into different cell types and diverge into epidermis-like cells and vascular-initial-like cells [80]. Furthermore, using single-cell RNA-seq data that provides numerous unbiased samples of transcriptomes at single-cell resolution, we constructed unique GRNs from each single-cell RNA-seq dataset to predict high-confidence gene regulations between transcription factors and their target genes. Upon query input, sub-networks are built with nodes for the query gene and genes having direct regulatory connections with the query (Fig. 2F).

Multi-gene analysis

Multi-gene analysis allows the simultaneous analysis of multiple queries. The Multi-gene analysis mode enables four different analyses (Gene expression heatmap, Co-expression network, DEG enrichment, and Single-cell analysis) for each mode of plant regeneration simply by selecting the menu tab and then submitting a list of genes by their gene locus codes (Fig. 3A). In the Gene expression heatmap analysis, users should submit multi-query genes and select a transcriptome dataset (mode of plant regeneration) for heatmap visualization (Fig. 3B). Rows are ordered by hierarchical clustering based on gene expression patterns, and subsets of the dendrogram can be magnified for a closer look. In this study, as an example, we queried CUC1, CUC2, ESR1/DRN, ESR2/DRNL, PLTs, STM, WOX5, and WUS, key regulators of hormone-induced plant regeneration, and showed that their expression patterns can be simultaneously compared during the process of hormone-induced plant regeneration (Fig. 3B). The results were consistent with previous reports [41], [1], [63], [26], [14], [27], [71].

Fig. 3

Interface and features of Multi-gene analysis. (A) Interface and input page for Multi-gene analysis. (B) Gene expression heatmap. Gene expression heatmap visualizes the expression patterns of the submitted genes in each dataset. Gene expression patterns of CUC1, CUC2, ESR1/DRN, ESR2/DRNL, PLTs, STM, WOX5, and WUS during hormone-induced plant regeneration are shown. (C) Sub-network constructed using key regulators of hormone-induced plant regeneration. The CUC1, CUC2, ESR1/DRN, ESR2/DRNL, PLTs, STM, WOX5, and WUS genes were submitted as query. Query genes and co-expressed genes are each colored violet and magenta, respectively. All plant regeneration-related bulk RNA-seq datasets were used to generate the sub-network. (D) GO enrichment analysis. GO analysis was conducted with all genes shown in (C). (E and F) DEG enrichment analysis. The DEG enrichment analyses were performed with top 200 genes targeted and activated by PLT2 (E; see also Table S2), and top 96 genes targeted and repressed by PLT2 (F; see also Table S2). (G and H) Dot plot and gene regulatory network (GRN) drawn with root stem cell regulators. The PLT1, PLT2, SCR, WOX5, and JKD genes were submitted as input queries to the Single-cell analysis in the Multi-gene analysis. Query genes in the GRN are colored violet, while other genes are colored magenta (H). The single-cell RNA-seq dataset was derived from hypocotyl callus. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Users can also submit multiple query genes to draw sub-networks that show correlations in gene expression between input query genes and other co-expressed genes with the selection of plant regeneration mode (Fig. 3C). Compared with Single-gene analysis, the construction of sub-networks containing multi-query genes employs a higher threshold of Mutual Rank to obtain a reasonable number of close neighbor genes for each query gene, generating a concise output sub-network. Genes consisting this co-expression network are used as input for GO enrichment analysis, and the results of analysis are generated as a table (Fig. 3D) (see below for validation of co-expression networks). Similar to Single-gene analysis, the genes in the network can also be used to perform DEG enrichment analysis.

In addition, users can investigate their GOIs by immediately running the DEG enrichment analysis in Multi-gene analysis. Over-enrichment of DEGs reveals the biological processes overrepresented in the input genes. To validate this approach, we queried top 200 PLT2-activated genes (Table S2) obtained from the PLT2-YFP chromatin immunoprecipitation sequencing (ChIP-seq) analysis [55]. The results showed that genes up-regulated in callus samples relative to their expression in wounding or shoot-regenerated samples (Fig. 3E) [20], [6] were enriched with the query genes, consistent with the pivotal role of PLT2 in callus formation and pluripotency acquisition. By contrast, upon querying 96 genes targeted and repressed by PLT2 (Table S2), the enriched sets of DEGs included genes down-regulated in callus samples compared to other samples (Fig. 3F).

Another feature of Multi-gene analysis is Single-cell analysis. Upon submission of multi-query genes, a dot plot displays the gene expression patterns of query genes in each cell type, the percentage of cells expressing the GOI, and the mean expression values of the GOIs. For instance, when we input several root stem-cell regulators, including PLT1, PLT2, SCR, WOX5, and JKD, the results showed that these genes are primarily expressed in QC-like cells (Fig. 3G), consistent with previous reports [63], [26], [27], [80]. In addition to browsing the expression patterns of query genes in single-cell RNA-seq datasets, users can draw GRNs with these genes (Fig. 3H). GO and DEG enrichment analyses can also be performed with the genes in this output GRN.

Transcriptome-level analysis

REGENOMICS provides transcriptome-level analyses as well. Users can select some or all of the Arabidopsis bulk RNA-seq datasets deposited in the REGENOMICS database (Table 1) to measure and visualize correlations between transcriptomes (Fig. 4A and B). Transcriptome-level correlation analysis is shown as an interactive heatmap generated with Spearman’s rank correlation and hierarchical clustering (Fig. 4B). In addition, different dimensionality reduction methods, such as principal component analysis (PCA; Fig. 4C) or multidimensional scaling (MDS; Fig. 4D), can be performed to allow the visualization of transcriptomes on two-dimensional (2D) or 3D scatter plots.

Fig. 4

Interface and features of Transcriptome analysis. (A) Interface of Transcriptome analysis. Users can select transcriptomes deposited in the REGENOMICS database (Table 1) and/or submit other custom transcriptome data in CSV format. (B) Heatmap of sample-to-sample correlation with hierarchical clustering. Correlation between samples are represented with colors according to the color scale. (C and D) Dimensionality reduction analysis. Principal component analysis (C) and multidimensional scaling analysis (D) were performed with transcriptomes datasets for hormone-induced plant regeneration transcriptomes and somatic embryogenesis transcriptomes. Each replicate is expressed as a dot on the scatter plot and colored according to its dataset. DAC, days after incubation on CIM; DAS, days after incubation on SIM; min, minutes after excision from seedling; HAI, hours after induction on E5 media. (E) Volcano plot of differentially expressed genes (DEGs). RNA-seq data of root samples harvested 0 h (0 h) and 1 h after wounding were used for DE analysis. Each gene is represented as a dot. Red and blue dots represent genes up-regulated and down-regulated, respectively, in 1 h after wounding compared with 0 h after wounding. (F) GO enrichment analysis performed with DEGs. The DEGs were obtained from (E). The analysis was performed separately for up-regulated genes. Top 10 GO terms were shown in the table. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Users can run these analyses with their own transcriptome data (in CSV format) and compare with transcriptomes from the REGENOMICS transcriptome database (Fig. 4A). Although the analysis takes any transcriptome data in the gene-sample matrix format, it is suggested that users follow the analysis pipeline used to process bulk RNA-seq data for REGENOMICS. We conducted correlation analysis using REGENOMICS-deposited bulk RNA-seq datasets of hormone-induced plant regeneration (Fig. 4B). Interestingly, transcriptome datasets were clustered into the stages of two-step plant regeneration, even though the datasets come from different studies (Fig. 4B), indicating that batch effects were resolved for transcriptome-level analysis. We also compared bulk RNA-seq datasets of hormone-induced plant regeneration with somatic embryogenesis transcriptomes using dimensionality reduction methods, PCA (Fig. 4C), and MDS (Fig. 4D). Notably, initial states of samples for hormone-induced plant regeneration and somatic embryogenesis were distant due to the different types of explants (root and hypocotyl explants for hormone-induced plant regeneration versus immature embryo for somatic embryogenesis). However, transcriptome data from later stages of callus formation (4 ∼ 14 DAC) and somatic embryogenesis (48 ∼ 72 HAI) showed high degree of correlation. This result was consistent with the fact that both processes involve high-level auxin-dependent cellular reprogramming. In addition, PCA analysis was performed for the comparison of hormone-induced plant regeneration datasets together with bulk RNA-seq datasets of different cell types and developmental zones of the root tip (Fig. S4). Notably, the transcriptomes from the root meristem (meristematic zone and quiescent center) were closely related to those of hormone-induced hypocotyl- and root-derived calli on CIM (Fig. S4), consistent with the fact that pluripotent calli predominantly express root stem-cell regulators [26], [27]. This analysis showed that transcriptomes of WOX5-expressing QC cells [32] are close to those of hormone-induced hypocotyl- and root-derived calli as well as hormone-induced immature embryos (Fig. S4). In contrast, transcriptomes of root maturation zones and differentiated cell types, such as those expressing S18 (mature xylem pole), COR (mature cortex), and COBL9 (root hair) [32], were the most distantly positioned to the transcriptomes of regenerating tissues (Fig. S4), which is consistent with the pluripotent nature of callus cells undergoing cellular reprogramming.

In addition, REGENOMICS supports the differential expression (DE) analysis using the deposited bulk RNA-seq data. Users can select two transcriptomes from a dataset deposited in REGENOMICS for the comparison. REGENOMICS identified DEGs between two samples, and the results were visualized as a downloadable volcano plot (Fig. 4E). GO analysis was also performed using each up-regulated and down-regulated genes, and the enriched GO terms could be displayed as a table (Fig. 4F).

Validation of co-expression networks in each mode of plant regeneration

One of the main features of REGENOMICS is to present co-expression networks, each of which is constructed based on the bulk RNA-seq transcriptome datasets related to a specific mode or to all modes of plant regeneration (wound-induced tissue repair, hormone-induced plant regeneration, somatic embryogenesis, and/or protoplast regeneration). To confirm whether the co-expression networks produced from different regeneration datasets are valid, we generated sub-networks using query genes that have been previously established as key regulators of each mode of plant regeneration.

ALF4, ERF115, WIND1, ASA1, YUC4, WOX11, ARF7, ARF19, and LBD16, that are known as key regulators of cellular reprogramming and callus formation at wound sites [24], [35], [17], [22], [81], [23], were selected as query genes to build a sub-network using wound-induced tissue repair dataset. Interestingly, the sub-network constructed using only wound-induced tissue repair transcriptome dataset displayed close connections among these query genes along with additional key regulators of wound responses (Fig. 5A). This co-expression sub-network was enriched with GO terms related to wounding responses and de novo root-regeneration processes, such as ‘auxin biosynthetic process’, and ‘lateral root formation’ (Fig. 5B). However, when datasets of other plant regeneration modes were used, the sub-networks were fragmented without higher-order assembly and showed no enrichment of wound-related GO terms (Fig. S5A and S5B).

Fig. 5

Validation of co-expression networks in each mode of plant regeneration. (A, C and E) Sub-networks generated with key regulators of each mode of plant regeneration. The sub-networks were drawn using bulk RNA-seq datasets related to wound-induced tissue repair (A), hormone-induced plant regeneration (C), and somatic embryogenesis (E). Input query genes (A, ALF4, ERF115, WIND1, ASA1, YUC4, WOX11, ARF7, ARF19, and LBD16; C, WUS, STM, CUC1, CUC2, PLT1, PLT2, PLT3, PLT5, PLT7, WOX5, ESR1/DRN, ESR2/DRNL, SCR, and JKD; E, ABI3, BBM, LEC1, LEC2, FUS3, AGL15, SERK1, GA2OX6, GA3OX1, IAA30, MYB118, WUS, YUC2, YUC4 and YUC10) are colored violet and co-expressed genes are colored magenta. (B, D and F) GO enrichment analysis. GO terms enriched in each output sub-network (A, C and E) are shown. A part of the significant GO term list is shown in B, D and F. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Similar results were obtained by submitting key regulators of hormone-induced de novo shoot regeneration, including WOX5, WUS, STM, ESR1/DORNRÖSCHEN (DRN), and ESR2/DRN-LIKE (DRNL) (Fig. 5C, S5C and S5D), as input query genes. The sub-network, which was constructed only with the hormone-induced plant regeneration datasets (Fig. 5C), showed the enrichment of GO terms such as ‘anatomical structure development’, ‘stem cell population maintenance’, and ‘shoot system development’ (Fig. 5D). Additionally, submission of key genes involved in somatic embryogenesis, including ABI3, BBM, LEC1, LEC2, FUS3, and AGL15, resulted in a sub-network enriched with genes associated with GO terms such as ‘developmental process’, ‘regeneration’, and ‘somatic embryogenesis’ (Fig. 5E and F). Sub-networks constructed with the same query genes using different datasets displayed different sub-networks with missing links between the key players (Fig. S5). Taken together, co-expression network analyses indicated that co-expressed genes depend heavily on the samples used as inputs and reflect the connections within a gene regulatory pathway specific to a given mode of plant regeneration.

Currently, several publicly accessible web applications containing extensive amounts of transcriptomes of diverse plant species are available for the examination of co-expression patterns of plant genes. In this regard, REGENOMICS enables plant regeneration-specific co-expression analyses and provides users with functional insights into the plant regeneration process. Notably, the co-expressions derived from the wound-induced tissue repair dataset were distinct from those derived from the co-expression networks publicly available at ATTED (https://atted.jp/) (Fig. S6A). Similarly, the sub-networks generated with hormone-induced plant regeneration and somatic embryogenesis datasets included co-expressions that could not be derived from co-expression network built with other transcriptome data (Fig. S6B and S6C).

As most of the molecular genetic studies conducted on plant regeneration to date have been performed using the model plant Arabidopsis, currently REGENOMICS database is mainly composed of Arabidopsis regeneration-specific RNA-seq datasets. However, we aim to annually expand and update the REGENOMICS web application by including transcriptome datasets obtained from other plant species. Thus, multi-species analyses might also be available in a near future to help finding key genetic factors conserved in plant regeneration. Further, we will attempt to perform multi-omics analyses by including genomics and epigenomics data (e.g. dynamics on DNA methylation and histone modifications during plant regeneration) to provide a comprehensive insight into the plant regeneration process. Through these updates, REGENOMICS will develop as a key resource in exploring the genetic components and pathways underlying plant regeneration.

Conclusions

Determining gene expression profiles and molecular networks provides important insights into gene functions in a given biological context. REGENOMICS is a collection of all available plant regeneration-related transcriptome datasets and implements easily accessible tools for the exploration of novel genetic components and connections underlying plant regeneration process. In addition, this web-based application permits user-friendly transcriptome-level analysis by using the databases constructed within REGENOMICS and customized user-submitted datasets to obtain systematic views on plant regeneration.

Methods

Data collection and processing

To conduct a comprehensive analysis of various plant regeneration processes, all currently available bulk RNA-seq datasets, generated from a total of 391 wound-induced tissue regeneration, hormone-induced plant regeneration, somatic embryogenesis, and protoplast regeneration experiments (Table 1), were downloaded from the NCBI Sequence Read Archive (SRA), DDBJ, EMBL-EBI, and NGDC databases. The raw RNA-seq data were processed using a customized pipeline (Fig. 1). The TAIR10 reference genome sequence and Araport11 annotations (https://www.arabidopsis.org/) for Arabidopsis, the ITAG4.0 genome and annotation (https://solgenomics.net/) for tomato, the MSU7.0 (http://rice.uga.edu/) genome and annotation for rice, and MpTak1_v5.1r2 (https://marchantia.info/) genome and annotation for Marchantia were used as reference for mapping reads with STAR aligner [7] and quantified and normalized with the transcripts per million (TPM) method with RSEM [31].

In addition, single-cell transcriptomes were collected and processed. Standardized workflows, including the 10X Genomics CellRanger pipeline, the Seurat package in R [56], and Scanpy in Python [74], were used for raw data processing and subsequent single-cell analysis. For single-cell transcriptomes generated with the Smart-Seq2 protocol, STAR aligner was used. UMI-tools was used to analyze single-cell transcriptomes generated by the BD Rhapsody system [62]. The TAIR10 and ITAG4.0 reference genomes were used for Arabidopsis and tomato datasets, respectively. Cell type-specific marker gene lists were prepared, based on previous studies [3], [4], [8], [32], [65], [72], [77], [82], [84], and were used with the cluster-specific marker genes given in the original study for each dataset. These marker genes were then used as input to score and annotate each cluster for each cell type using scanpy.tl.score_genes provided in Scanpy.

Gene expression browser and heatmap

Gene counts normalized with the TPM normalization method were displayed on the Gene expression browser in Single-gene analysis. Plotting was done with the ggplot2 package in R. For the Gene expression heatmap in Multi-gene analysis, TPM values were scaled with StandardScaler and MinMaxScaler from the scikit-learn Python module [48], and displayed as Z-scores.

Gene co-expression network analysis and Sub-network generation

Gene counts normalized with the TPM normalization method were log-transformed as log2(x + 1) to create gene expression matrices for each type of plant regeneration. Next, the resulting matrix was used to calculate Pearson correlation coefficients between each pair of gene. The ranks of correlation coefficients were then measured and used to obtain the geometric mean, the Mutual Rank (MR) [44]. Co-expression networks were first generated with all transcriptomes in the REGENOMICS database as input for detecting gene correlations common to a variety of regeneration methods. Then, datasets corresponding to the four plant regeneration modes were used separately to create a co-expression network specific to each mode of regeneration. Sub-networks were generated based on the query gene(s), connecting the top 20 closest neighbor genes (based on the MR values) with significant co-expression values. When constructing sub-networks with multiple query genes, only close neighbor genes with strong co-expression (MR ≤ 30) were used to build the networks in order to avoid creating an overly dense co-expression network. Sub-networks were generated using the Pandas, NetworkX [11], and ForceAtlas2 Python packages, and visualized with d3.js (https://d3js.org/), ObservableHQ (https://observablehq.com), and Sigma.js (https://www.sigmajs.org/).

Gene regulatory network inference and sub-network generation

Gene regulatory network inference for each scRNA-seq dataset was performed using the SCENIC workflow [68]. Pre-processing was performed in Scanpy, and the network was inferred using GRNBoost2 [42]. Transcription factor gene list was obtained from PlantTFDB 4.0 [25]. The directed edges between nodes were ranked by importance, and the top 10 edges connected with the query gene were first taken to build the sub-network in Single-gene analysis. Directly linked neighbor genes were further taken to build an extensive gene regulatory network. For Multi-gene analysis, only gene connections with an importance score of 10 or greater were used to build the sub-networks. Visualization of the sub-networks were done using the same tools and packages as co-expression networks.

Single-cell RNA-seq visualization and trajectory inference analysis

Expression of query gene(s) in Single cell analysis was visualized with Uniform Manifold Approximation and Projection (UMAP) plots, violin plots, and dot plots, which were drawn using the Scanpy toolkit. Trajectory inference was also performed in Scanpy, using scanpy.tl.dpt to compute pseudotime for each gene. The selection of root cell type was done as shown in each original study, and the visualization of the pseudotime was give as a UMAP plot.

GO and DEG enrichment analyses

GO enrichment analysis was conducted using the GOATOOLS Python library [28]. Only protein-coding genes were used as input for analysis. The table of enriched GO terms was filtered, showing terms with a p-value <0.05 and those categorized as biological processes. To perform DEG enrichment analysis, genes with fold-change >2 and p-value <0.05 were identified as DEGs with edgeR [50]. The significance of enrichment was tested with Fisher’s exact test, and p-values were corrected using the Benjamini and Hochberg method [2].

Correlation heatmap and dimensionality reduction

Spearman’s rank correlation coefficients were calculated with the Pandas Python package, and visualized with InCHlib [61] and ObservableHQ. Dimensionality reduction was conducted as PCA and MDS with the Pandas, NumPy [12], and scikit-learn Python packages. Results were visualized as 2D and 3D scatter plots with the Plotly Python package.

Differential expression analysis

Pre-processed RNA-seq data were used for differential expression (DE) analysis. Using edgeR, the reads from the selected RNA-seq data were normalized with CPM, and DEGs were identified if the fold-change >2 and p-value <0.05. DEGs were visualized as a volcano plot with Plotly.

Application architecture

REGENOMICS was built based on the Django web framework (https://www.djangoproject.com) on a Linux platform. Python and R scripts were implemented for the analysis pipeline. STAR, RSEM, CellRanger, UMI-tools, Seurat, and Scanpy were used for pre-processing datasets. The user-interface for REGENOMICS was built with HTML, JavaScript, and Semantic UI (https://semantic-ui.com/). The results were visualized mainly with d3.js, sigma.js, and plotly.js (https://plotly.com/). The co-expression network (https://observablehq.com/@tnsgud996/co-expression-network) and gene regulatory network (https://observablehq.com/@tnsgud996/gene-regulatory-network) for Single-gene analysis, and colored heatmap (https://observablehq.com/@tnsgud996/correlation-heatmap) for Gene expression heatmap in Multi-gene analysis were built with and are available on ObservableHQ.

Funding

This work was supported by the under Project Number SSTF-BA2001-10. S.H.B. is grateful for financial support from .

Author contributions

P.J.S. and Y.S.N. conceived the project; S.H.B. constructed the database and designed the web application; P.J.S. revised the web application; S.H.B. wrote the first draft of the manuscript; P.J.S. and Y.S.N. revised the manuscript. All authors read and approved the final manuscript.