Uncovering transcriptional reprogramming during callus development in soybean: insights and implications

Callus, a valuable tool in plant genetic engineering, originates from dedifferentiated cells. While transcriptional reprogramming during callus formation has been extensively studied in Arabidopsis thaliana, our knowledge of this process in other species, such as Glycine max, remains limited. To bridge this gap, our study focused on conducting a time-series transcriptome analysis of soybean callus cultured for various durations (0, 1, 7, 14, 28, and 42 days) on a callus induction medium following wounding with the attempt of identifying genes that play key roles during callus formation. As the result, we detected a total of 27,639 alterations in gene expression during callus formation, which could be categorized into eight distinct clusters. Gene ontology analysis revealed that genes associated with hormones, cell wall modification, and cell cycle underwent transcriptional reprogramming throughout callus formation. Furthermore, by scrutinizing the expression patterns of genes related to hormones, cell cycle, cell wall, and transcription factors, we discovered that auxin, cytokinin, and brassinosteroid signaling pathways activate genes involved in both root and shoot meristem development during callus formation. In summary, our transcriptome analysis provides significant insights into the molecular mechanisms governing callus formation in soybean. The information obtained from this study contributes to a deeper understanding of this intricate process and paves the way for further investigation in the field.


Introduction
Plant cells exhibit a notably superior capacity for regeneration when compared to animal cells. This remarkable regenerative potential arises from their developmental plasticity, which is evident through the process of dedifferentiation and transdifferentiation. These abilities are prominently shown in phenomena like callus formation and grafting. Particularly, when plants are exposed to stresses such as wounding or pathogen infection, they can generate unorganized masses of cells known as callus (Ikeuchi et al., 2013). Initially, the term "callus" described a substantial accumulation of cells producing callose at a wound site. However, its current definition encompasses disorganized clusters of cells that arise from alterations in cellular state, such as excessive proliferation and the transdifferentiation of specialized cells (Feheŕ, 2019). Callus can be developed from a single somatic cell that has already differentiated and may possess totipotency, which means that each callus cell has the potential capability of regenerating an entire plant via somatic embryogenesis, a process similar to embryonic development in non-zygotic cells (Ikeuchi et al., 2019;Bidabadi and Jain, 2020). During this process, the fate of cells can be reprogrammed depending on the types of imposed stresses and the relative concentrations of phytohormones. The innate regenerative capability of plants has made it possible to introduce functionally useful genes into a callus, usually using Agrobacterium, and subsequently produce beneficial substances in plants, making it a valuable tool in the field of plant-derived biotechnologies (Belide et al., 2017;Hwang et al., 2019;Cordeiro et al., 2023).
Auxin and cytokinin are two phytohormones that play central roles in controlling callus development. The ratios of these hormones are believed to determine the direction of the cell's developmental fate (Skoog and Miller, 1957). Generally, a higher concentration of auxin relative to cytokinin stimulates root formation, while a reversed ratio favors shoot regeneration. In vitro, callus induction is typically conducted on callus-inducing media (CIM) with a proper balance between auxin and cytokinin. The gene expression pattern of CIM-induced callus is known to resemble root meristem development (Sugimoto et al., 2010). Auxin promotes callus formation by activating genes such as lateral organ boundary domain 16 (LBD16), LBD17, and LBD18, which are mediated by auxin response factors (ARFs) like ARF7 and ARF19 (Okushima et al., 2007). LBDs, in turn, activate the expression of genes that promote cell proliferation and modify cell wall structures. Unlike the CIM-induced pathway, wounding can induce callus formation through a pathway that does not follow root initiation. Wounding activates the central regulator wound induced dedifferentiation 1 (WIND1) and its paralogous genes WIND2, WIND3, and WIND4 (Iwase et al., 2011). The wounding signal triggers cytokinin biosynthetic pathway genes such as isopentenyl transferase 3 (IPT3) and lonely guy (LOG1, LOG4, and LOG5) (Ikeuchi et al., 2017). This, in turn, activates type-B Arabidopsis response regulator 1 (ARR1) and ARR12 involved in cytokinin signaling, leading to cell cycle reentry by cyclin D3;1 (CYCD3;1).
Soybean [Glycine max (L.) Merr.], a legume crop, possesses a diploidized paleopolyploid genome that consists of 20 chromosome pairs that underwent the latest whole genome duplication (WGD) about 13 million years ago (MYA). It is widely recognized as a major source of dietary protein and oil for humans. Soybean, along with other legumes, is renowned for its ability to fix nitrogen through symbiosis with Rhizobium bacteria. The whole genome of soybean, estimated to be 1.115 gigabase (Gb) in size, has been fully sequenced, and a recent transcriptome analysis revealed 54,132 confidently expressed genes. The prevalence of paralogous genes and the existence of a substantial proportion of multi-gene families (~75%) in the soybean genome is largely attributed to its WGD nature (Schmutz et al., 2010;Wang et al., 2014).
This study aims to uncover the biological mechanism underlying soybean callus induction and developmental process. To achieve this goal, the transcriptome changes associated with the soybean callus formation process were comprehensively analyzed. This involves dividing the development of soybean callus into seven time points (0, 1, 4, 7, 14, 28, and 42 days) to observe how the transcriptome is affected by morphological changes over time. To confirm the physiological mechanisms related to soybean callus formation, we analyzed the transcriptomic dynamics using various methods such as differentially expressed genes (DEGs), c-mean clustering, and gene ontology (GO) enrichment for genes related to hormones, cell wall, cell cycle, and transcription factors.

Plant materials and callus induction
Callus induction was performed using seeds of soybean (G. max L. Merr, cv. Kwang-an). The surface of seeds was sterilized using chlorine gas, which was generated by combining 5ml of 12N HCl and 45ml of 12% hypochlorite solution in a dessicator for a duration of 16 hours. After that, the seeds were subjected to incubation with 12ml of deionized water (DIW) per 150 seeds. After 24 hours of moist incubation, the seeds were washed with a 1% (v/v) mild bleach solution, followed by three additional washes with autoclaved DIW. The seeds were further incubated for one day in a 50 ml conical tube filled with water. After soaking, the seeds were halved using scalpel (#11 blade), and the embryonic axes were removed. To induce the callus formation, the divided seeds underwent surface scratching in multiple directions, with each seed being scratched 10 times. Additionally, the junction of the hypocotyl and cotyledon was scratched perpendicular to the cotyledon (Jiang et al., 2010). Following that, the divided seeds were implanted with CIM (Gamborg's B5 medium) containing 3 mM MES, 3% (w/v) sucrose, 8 mg/ml of 2,4-D, 1 mg/ml of BA, and 250 mg/ml of cefotaxime. To activate callus induction, the plants were grown in darkness (24 hours dark) at 28°C. To prevent nutrient depletion, the plants were sub-cultured every two weeks. Callus formation, consisting of meristem along with a small portion of cotyledon and hypocotyl, was observed at various time points after wounding; 1, 4, 7, 14, 28 and 42 days after wound (DAW). In contrast, a control sample at 0 DAW remained unwounded. All samples were collected at room temperature, rapidly frozen in liquid nitrogen, and stored at -80°C until RNA extraction.

RNA extraction and sequencing
Total RNA was extracted from the explants at 0, 1, 4, 7, 14, 28, and 42 DAW. For RNA sequencing, the RNA samples were pooled from 10 explants at each time point and isolated using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). The RNA integrity was evaluated through agarose gel electrophoresis, where the ratio of the large subunit and small subunit of rRNA was examined. The analysis revealed consistent intensity across samples, indicating the preservation of RNA integrity. Moreover, the purity of RNA samples was assessed by measuring the A260/A280 ratio using a spectrophotometer (Implen, Munchen, Germany). The obtained results confirmed that both the integrity and purity of the RNA were well-maintained, as the A260/A280 ratio showed minimal deviation from expected range of 1.8~2.0. Sequencing libraries were prepared for each sample using 2 mg of total RNA and the TrueSeq Stranded Total RNA LT Prep Kit (Plant) (Illumina, CA, USA) following the manufacturer's instructions. Paired-end reads for the seven samples were generated using the Illumina NovaSeq 6000 platform (Macrogen, Seoul, ROK) with three biological replicates.

Transcriptome analysis
The sequencing reads generated from RNA sequencing were aligned to the soybean reference genome (Wm82.a2.v1) by utilizing the HISAT2 aligner 1 (Kim et al., 2019) with the default settings. The raw read counts of each gene were determined using HTSeq 2 (Anders et al., 2015) with the options "-mode union -stranded no -format bam". The DESeq2 R package (Love et al., 2014) was used to calculate the normalized count and identify DEGs at each time point compared to day 0, using the criteria of |log2fold change| ≥ 2 and false discovery rate (FDR) adjusted p-value < 0.05. To analyze the relationships among samples, principal component analysis (PCA) was performed using the regularized logarithms (rlog) of read counts calculated by DESeq2. The Venn diagram and heatmap of DEGs at each time point were generated using InteractiVenn 3 (Heberle et al., 2015) and the pHeatmap R package (Kolde, 2015), respectively. Mfuzz R package (Kumar and Futschik, 2007) was used to cluster the expression levels of all significantly changed DEGs, based on c-mean clustering. The optimal fuzzifier m value, which prevented clustering of random gene expression, was estimated as 1.53.

Gene ontology and KEGG enrichment analysis
The enrichment analysis of GO and KEGG pathways was performed using the clusterProfiler R package (Wu et al., 2021). The analysis was conducted with the criterion of p-value < 0.005 to identify the functional categories among the DEGs. The annotation information and GO data were obtained from Phytozome 4 (Goodstein et al., 2012) and PlantRegMap 5 (Tian et al., 2019), respectively. The organism database used in GO enrichment analysis was generated using the AnnotationForge R package. The enzyme information of the induced genes was obtained from the locus mapping table in TGIL platform, as the Phytozome and KEGG gene IDs did not match for KEGG term enrichment. The visualization of enriched terms for each cluster was generated using the dot plot method in clusterProfiler.

Analysis of transcriptional regulation
We obtained information on all soybean transcription factors (TFs) and their interactions from PlantTFDB v4.0 6 (Jin et al., 2017) and PlantRegMap, respectively. To identify significantly enriched TF families, we used the enrichment method in PlantRegMap, which revealed a total of 18 enriched families.

Real-time quantitative PCR analysis
To validate the expression patterns obtained from RNA sequencing, we selected nine genes known to be involved in callus generation in A. thaliana (Ikeuchi et al., 2017). Orthologous genes in soybean were identified using the TGIL platform 7 . One microgram of total RNA was reverse transcribed into cDNA using the SuPrimeScript RT Premix with oligo (dT) 2X kit (GeNet Bio, Korea). qRT-PCR was performed with the following steps: predenaturation at 95°C for 30 seconds, denaturation at 95°C for 5 seconds, and extension at 60°C for 40 seconds. Denaturation and extension were repeated for 40 cycles, and the melt curve was obtained by increasing the temperature incrementally by 0.5°C from 65°C to 95°C at the final step. The delta-delta Ct method (Gao et al., 2017) was used for the relative quantification of gene expression patterns, with tubulin beta 3 (TUBB3) selected as the control gene for normalization (Rao et al., 2013). All primers were designed automatically using an in-house Python program and are listed in Table S1. PCR amplifications were performed with three biological replicates and three technical replicates.

Morphological alterations during the callus development
Significant morphological changes were observed during callus formation (Figure 1), which can be categorized into three stages. The first stage, which lasts from 0 DAW to 4 DAW samples (Figure 1Aa-f), is known as the lag phase. During this phase, the callus is visible (Figures 1Ae, f), and its mass can be measured for the first time ( Figure 1B). The second stage, which begins on day 4, is the log phase. At this stage, the callus grows actively and rapidly, appearing translucent yellow, fragile, and watery. The callus continues to grow fresh until days 7-8 (Figures 1Ag, h, B) as the callus cell divides actively and rapidly. The final stage is the stationary phase, which begins after days 13-14 (Figures 1Ai-n). At this stage, half-seed still grows, but the texture of the callus changes completely. The surface becomes stiff, dark, and opaque (Figures 1Ak-n). The callus appears to be dying, while the cotyledon grows quickly and becomes senescent.

Genome-wide gene expression patterns during callus formation
In order to perform a comprehensive transcriptome analysis of soybean callus formation, we collected a total of 21 samples at seven distinct time points: 0, 1, 4, 7, 14, 28, and 42 days. Each time point consisted of three biological replicate to ensure a robust statistical analysis. This design allowed us to capture the global transcriptome changes occurring during the process of soybean callus formation. The RNA sequencing was conducted using the Illumina NovaSeq 6000 platform, which generated a total of 661 Gb obtained from paired-end reads. After filtering low-quality and adaptor-polluted reads, 84.76% of the clean reads with >Q30 Phred score were mapped to the soybean reference genome (Gmax.v.2.0). The average ratio of concordant pair alignment was 96.31% (ranging from 90.51% to 97.62%) (Table S2). Principal components analysis (PCA) was conducted to evaluate the similarities between individual samples. The results showed that samples from different time points were clearly separated along the first and second principal components, which explained 56% and 28% of the total variation, respectively ( Figure 2). Moreover, the biological replicates at each time point were clustered together, indicating the high reproducibility of the experiment. Additionally, the dot plot revealed a significant difference in sample similarity between 0, 1, 4, 7, and 14 days based on PC1, and the PC2 value gradually increased until day 7 and then decreased. Thus, we hypothesized that genes related to callus formation, cell division, and response to hormone and wounding stress were strongly expressed up to 7 days, while genes involved in plant growth and senescence were activated thereafter. To validate this hypothesis, we conducted GO enrichment analysis using the top 1000 genes associated with PC1 and PC2, respectively. Our results showed that GO terms related to stress response and senescence were overrepresented in PC1, whereas GO terms involved in cell wall modification and response to auxin were highly enriched in PC2 ( Figure S1). Based on these findings, we could conclude that genes associated with stress response were significantly activated by wounding stress on day 1, followed by genes related to callus formation up to day 7, and genes involved in growth and senescence thereafter.
DESeq2 was utilized to compute expression values for each sample. From the total of 42,466 genes that had normalized expression values greater than 3 in at least one sample, DEGs were identified by comparing each time point sample to day 0 (|log2fold change| >= 1 and FDR < 0.05) (Table S3). A total of 27,639 genes exhibited perturbations in their expression pattern during callus formation. Of these,9,455,11,747,13,087,12,163,11,260,and 10,866 genes were upregulated for 1,4,7,14,28,and 42 days,respectively,while 5,232,6,799,7,284,6,795,6,826,and 7,267 genes were down-regulated for the same respective time points (Figure 2). To validate the expression values obtained from RNA-seq, nine genes from the up-regulated genes were selected for qRT-PCR analysis, and the results were comparable to those obtained from RNA-Seq. This finding, presented in Figure S2, supports the robustness of the RNA-Seq experiment. Morphological changes and growth curve during callus induction of soybean. (A) Morphological changes during callus formation. The left half of each image shows seed sections, with samples for RNA-seq (including meristem, part of hypocotyl, and cotyledon) on the right (samples labeled b, d, f, h, j, l, and n). Samples a and b were not wounded and served as control samples, while samples c-n were scratched using a surgical blade and incubated on CIM for varying durations. The time points captured in the images are as follows: 1DAW (days after wounding; c and d), 4DAW (e and f), 7DAW (g and h), 14DAW (i and j), 28DAW (k and l), 42DAW (m and n). Each bar in the images corresponds to 10mm. (B) Growth curve during callus development. The plot shows the callus growth curve, with ten explants pooled at each time point and the weight of each callus sample measured. The blue line represents actual samples used for RNA-sequencing (corresponding to samples b, d, f, h, j, l, and n in panel A), while the orange line represents only callus collected by raking up using a blade, and the gray represents a part of hypocotyl and cotyledon excluding callus.

Clustering and GO enrichment analysis of callus development-related genes
To better understand the functional meaning of transcriptomic changes during callus formation, the Mfuzz R program was used to group the 27,639 induced genes into eight clusters based on their time-series expression patterns ( Figure 3A; Table S4). These clusters were then subjected to downstream analysis, such as GO enrichment, to assign functional criteria to each cluster (Table S5).
Cluster 1, which consisted of 2,599 genes, showed a rapid increase in expression at day 1 followed by a gradual decrease. This cluster mainly contained genes expressed in response to wounding and callus induction. The results of the GO enrichment analysis indicated the involvement of biological processes such as auxin-activated signaling pathway, regulation of hormone levels, and lateral root formation ( Figure 3; Table S5), suggesting that callus initiation process responsive to wounding and auxin was turned on at day 1.
Cluster 2, consisting of 3,630 genes, exhibited strong expression on days 4 and 7, which gradually decreased over time. GO categories related to cell proliferation, including DNA replication, DNA recombination, cell division, cell cycle process, and nuclear division, were found to be overrepresented. This cluster also contained genes associated with organ development, such as multi-organism cellular process, cell wall organization, root development, and regulation of root meristem growth, indicating that cell division and proliferation are critical in generating callus. Furthermore, molecular function analysis identified microtubule binding, structural constituent of cytoskeleton, and auxin transmembrane transporter activity. This result supports that the cell cycle for callus expansion and proliferation was activated from day 4, and that this cluster could be an early stage for cell wall organization.
Cluster 3 comprised 3,053 genes, with highly expressed levels at day 7. The biological processes associated with this cluster were related to secondary metabolites (i.e., secondary metabolite biosynthetic process) and stress response (reactive oxygen species metabolic process, oxalate metabolic process, and cellular oxidant Clustering and GO enrichment analysis of RNA-seq data during callus formation. (A) Clustering of differentially expressed genes. The eight clusters, based on time-series expression patterns at Day 1, Day 4, Day 7, Day 14, Day 28, and Day 42 during callus formation, were grouped by induced genes. The fuzzy c-means (FCM) clustering algorithm, implemented in the Mfuzz R package, was used to cluster genes based on standardized expression change values. A membership value between 0 and 1 is assigned to each gene, and cluster cores are highlighted with gray line. (B) GO and KEGG Pathway enrichment analysis for clusters. The dot plot shows the top 15 overrepresented GO terms of biological process, cell component, and molecular function identified using clusterProfiler with differentially expressed genes for each cluster. detoxification) as confirmed by GO analysis. Notably, this cluster contained genes involved in cell wall biogenesis, cellular polysaccharide metabolic process, lignin biosynthetic process, glucan metabolic process, and cell wall modification, suggesting active participation of genes in cell wall modification for callus formation at day 7. Additionally, KEGG pathway analysis revealed the presence of phenylpropanoid biosynthesis, starch and sucrose metabolism, glycine, serine and threonine metabolism, and biosynthesis of various plant secondary metabolites. Therefore, this finding suggests that the period leading up to day 7 played a crucial role in callus formation by triggering cell cycle activation and cell wall modification.
Cluster 4 consisted of 3,141 genes that exhibited a gradual increase until day 14, followed by rapid expression on day 28, and then a decrease in expression. GO terms related to stress response, such as response to toxic substances, cellular oxidant detoxification, reactive oxygen species metabolic process, and cellular oxidant detoxification, were identified in biological processes. KEGG pathway analysis revealed the presence of phenylpropanoid biosynthesis, MAPK signaling pathway, and plant-pathogen interaction. These findings suggest that genes involved in cell development are down-regulated, while stress-responsive genes are activated after day 14 during callus formation.
Clusters 5 and 6 exhibited gradually increasing expression patterns and contained 3,092 and 2,586 genes, respectively. Cluster 5 was associated with gene silencing by RNA, organophosphate biosynthetic process, carbohydrate derivative biosynthetic process, DNA modification, and DNA methylation or demethylation. Cluster 6 was associated with the recognition of pollen and cell, and pollen-pistil interaction. In contrast, clusters 7 and 8 demonstrated a sharp decline starting from day 1. Cluster 7 consisted of the highest number of genes (6,580). The GO terms associated with response to abscisic acid, regulation of the developmental process, negative regulation of gene expression, and negative regulation of biosynthetic process were identified in the biological process. On the other hand, cluster 8, consisting of 2,958 genes, exhibited a consistent decrease in expression levels starting from day 1. Overrepresented GO terms in this cluster included translational initiation, pyruvate metabolic process, response to cytokinin, and photosynthesis. Furthermore, the metabolic pathway analysis confirmed the involvement of carbon metabolism, glycolysis/gluconeogenesis, and carbon fixation in photosynthetic organisms. These findings suggest that the physiological mechanisms associated with callus formation are suppressed to halt unnecessary cellular events and conserve metabolic energy.
Furthermore, we observed the expressional dynamics of several genes involved in biosynthesis, transport, and signaling pathways of other phytohormones -ethylene, gibberellic acid, abscisic acid, and jasmonic acid -during callus induction. The expression patterns of these genes appeared to be perturbated, indicating the involvement of other hormone signals, in addition to auxin, CK, and BR, in the callus formation.

Analysis of transcriptional regulation related to cell cycle and cell wall development
Because cell cycle and cell wall modification were known as key regulators during callus induction, we examined the transcriptional dynamics to gain insight into their role in callus formation (Table  S7). We identified genes encoding cyclin (CYC), cyclin-dependent kinase (CDK), CDK inhibitor (CKI), retinoblastoma-related (RBR), E2F, and DP (Inzéand De Veylder, 2006;Qi and Zhang, 2020;Shimotohno et al., 2021) (Figure 5A). CYC genes in plants are encoded by a multi-gene superfamily, which is categorized into seven classes: A, B, C, D, H, P, and T. In our analysis, we identified 7, 8, 2, 14, 2, 11, and 3 orthologous genes for each respective class. Among them, several genes encoding D-type, A-type, and B-type CYCs play a role in the G1-to-S, S-to-M, and G2-to-M transition, respectively. Interestingly, these genes exhibited higher expression levels during callus formation compared to other members of the CYC family in our study. The CDK gene, known to interact with CYC genes, is classified into 11 families designated from A to K. Among 21 orthologous CDK genes, CDKB genes responsible for facilitating the transition from G2 to M phase of the cell cycle through interaction with CYCA, CYCB, or CYCD displayed notably high expression levels. Taken together, these findings indicate that G2-to-M phase, which is the process that drives cell division in the cell cycle, is dynamically activated during callus formation. Kiprelated protein (KRP) and WEE1, which are included in the CKI family, act as negative regulators in the G1/S transition and G2/M transition, respectively. Contrary to our expectations, these genes were up-regulated during callus formation. This phenomenon may be attributed to the dynamic expression patterns of cell cycle-related genes, which exhibit changes that occur at a faster rate than observable time intervals between samples. Additionally, the samples used in the study were derived from a mixture of different callus cell types, contributing to the observed variability in gene expression. Therefore, the highly expressed pattern of CKI genes may indicate that the cell cycle is dynamically activated by these genes.
To uncover the transcription regulation in the cell wall coupled with the cell cycle, we investigated several genes involved in cell wall synthesis and modification ( Figure 5B). The cellulose synthase (CesA) and cellulose synthase-like protein (CLS) are involved in cell wall synthesis. We validated that 13 CesA and 18 CLS genes exhibited high expression levels during callus formation. Furthermore, we discovered various genes associated with cell wall modification, such as pectin methylesterase (PME), xyloglucan endotransglycosylase/hydrolase (XTH), polygalacturonase (PG), expansin (EXP), and pectin lyases (PL). This observation suggests that these genes undergo substantial regulation during callus formation, indicating that wounding stimuli influence the process of cell wall synthesis and modification.

Identification of transcription factors and associated genes in relation to callus formation
In order to explore the transcriptional regulation changes during callus formation, we obtained a comprehensive set of 3,866 TFs. Among them, we identified 1,309 differentially expressed TFs belonging to 17 distinct TF families (Table S8). The largest TF families among these genes were the basic/helix loop helix (bHLH), comprising 239 genes, and the ethylene responsive factor (ERF) family, which consisted of 175 genes. These families exhibited significant changes in gene expression following induction. Other TF families also exhibited significant expression changes, including 178 MYBs, 141 NACs, 128 WRKYs, and 99 bZIPs ( Figure 6A). To explore the relationship between clusters and TF families, we conducted an enrichment analysis using hypergeometric distribution in hyperR 8 . This analysis involved assessing the overrepresentation of TF families within each cluster ( Figure 6B).
To ascertain TFs responsible for regulating the transcriptional dynamics of differentially expressed genes within each cluster during callus formation, we performed TF enrichment analysis using genes from each cluster in the PlantRegMap database. Furthermore, we conducted GO enrichment analysis to explore the enriched TFs and their target genes. Additionally, we employed gene regulatory network analysis to unravel the biological significance of these findings (Figures 6C-E). A total of 67 TFs were identified under the following conditions: p-value < 0.05 and inclusion in clusters 1-6. Cluster 1 was governed by 16 TFs associated with biological processes such as secondary cell wall biogenesis, organic cyclic compound, and regulation of DNAtemplated transcription. On the other hand, cluster 2, regulated by 19 TFs, did not exhibited GO terms. However, we verified that ARF TFs and B3 TFs, specifically present in clusters 5 and 3, respectively, exclusively targeted this cluster. Cluster 3 was governed by 35 TFs associated with GO terms related to callus development, including secondary cell wall biogenesis, root development, cell differentiation, and tissue development. In contrast, cluster 4, controlled by 25 TFs, is involved in senescence-related processes such as leaf senescence and plant organ senescence. Cluster 5 consisted of 8 TFs, including 4 Dofs, 1 AP2, 1 E2F/DP, 1 MYB, and 1 WOX. Notably, WOX TF specifically targeted cluster 5 and was not found in any other clusters. On the other hand, cluster 6 was regulated by 27 TFs and exhibited a high representation of GO terms related to cell development, similar to cluster 3.
Our gene regulatory network analysis uncovered the specific biological processes in each cluster that were targeted by enriched TFs. Out of 27 enriched NAC TFs, we discovered the orthologous gene (Glyma.17G154100) corresponding to ANAC011/ANAC071/ ANAC096. This key regulator, known to respond to wounding signals, governs the expression of 126 genes in cluster 3 and 119 genes in cluster 6. AP2 (Glyma.09G248200), an orthologous gene of BBM, was found within cluster 1 and exerted regulatory control over 1299, 1043, and 858 genes in clusters 2, 5, and 6, respectively. These genes are associated with auxin-activated signaling pathway, cell wall organization, cell cycle process, and cell division. ARF8 (Glyma.02G239600) governs the regulation of 54 genes within cluster 2, encompassing processes such as auxin-activated signaling pathway, microtubule-based process, and cytoskeleton organization. BES1 (Glyma.06G034000) exhibits significant upregulation during callus formation and controls the regulation of 121 genes within cluster 1. These genes are involved in various processes such as photosynthesis, light harvesting in photosystem I, rRNA processing, and rRNA metabolic process. OBF-interacting protein 1 (OBP1; Glyma.02G062700), a member of the DOF family, is present in cluster 5 and exerts regulatory control over clusters 3, 5, and 6. It governs genes associated with cell wall biogenesis and cell cycle process. E2F/DP (Glyma.06G024900 and Glyma.17G093600), which plays a role in the cell cycle, is present in cluster 2 and regulates clusters 2 and 5. Additionally, B3 (Glyma.14G079100) exhibits high expression and governs the regulation of 133 genes within cluster 2. Overall, these transcription factors control various biological processes such as DNA replication, cell cycle, nucleosome assembly, and cell division. WUSCHEL (WUS; Glyma.01G166800) governs the regulation of 89 genes in cluster 5, which are involved in processes such as organophosphate biosynthesis, nucleoside phosphate metabolism, and carbohydrate derivative biosynthesis.

Analysis of gene expression changes during the initial 7 days of callus formation
To explore the phenomenon that the day 7 acted as the boundary point of callus morphology from cell division and growth to senescence (Figures 1Ag,h, 2A), we analyzed the expression pattern comparing each sample to day 7 and filtering with the following condition: the first group (log2FC < -1 at 0 DAW, log2FC > 1 from day 0 to 4, and log2FC < -1 from day 14 to 42) and the second group (log2FC < -1.5 from day 0 to 7 and log2FC > 1.5 from day 7 to 42) ( Figure 7B). We identified 76 and 291 genes involved in the two groups, respectively. In order to perform a more detailed analysis of these patterns, we took 13 representative categories from the gene ontology and manually classified them into two groups (Figures 7A, C; Table S9).
Group 1 was highly related to cell homeostasis (GO 01), hormone response (GO 05), cellular metabolism (GO 11), and development (GO 12), including genes related to cell multiplication and enlargement similar to cell growth such as deeper rooting 1-like (DRO1-like), CYCA2;1, WOX5, and EXP. The early stages of callus development seem to involve an increase in cell size during the log phase. EXP (Glyma.19G222900) was expressed in response to wounding and increases cell size by loosening the cell wall (Cho and Cosgrove, 2000). When callus is generated in auxin-rich conditions, it follows a root developmental program (Ikeuchi et al., 2017). One of the representative genes for root formation, DRO1-like (Glyma.07g040800), was highly expressed and contributes to root meristem development, leading to the organized expression of regulators of root tissue (Grimplet et al., 2017;Ikeuchi et al., 2017). In response to hormones, particularly auxin, small auxin up-regulated RNA (SAUR; Glyma.02G049200), type I inositol polyphosphate 5-phosphatase 4 (IP5P4; Glyma.02G145600), and naked pins in YUC mutants 4 (NPY4; Glyma.17G161500) were up-regulated in the early stages and played a role in development ( Figure 7B). Contrary to group 1, group 2 was over-represented genes involved in signal perception and transduction (GO 04), including stress-responsive factors such as leucine-rich receptor-like kinase (LRR), cysteine-rich receptor-like protein kinase (CRR), and G-type lectin S-receptor-like serine/ threonine-protein kinase (GsSRK). These genes are mainly stimulated by stress-responsive phytohormones such as ABA and shown. Each group is filtered based on the specified criteria mentioned below. (C) The percentage of genes associated with each GO category is separated into two groups. The hierarchical diagram for each GO category is displayed on the left side. Genes belonging to group 1 are represented in blue, while genes belonging to group 2 are shown in orange. Park et al. 10.3389/fpls.2023.1239917 Frontiers in Plant Science frontiersin.org ethylene (Sun et al., 2017). These results suggest that the genes involved in root formation, cell wall expansion, and auxin response were highly regulated for callus initiation and proliferation during the initial 7 days, whereas stress response was highly activated after day 7.

Discussion
Studies on callus formation have described two separate mechanisms: the wounding-induced pathway and the CIMinduced pathway (Ikeuchi et al., 2017). However, since all biological processes in plants are interconnected rather than being separate, the pathways related to callus formation are expected to be interlocked as well. When responding to wounding stress, cell reprogramming is activated rapidly and strongly to generate callus, and then the balance between two phytohormones, auxin and cytokinin, determines cell fate (Iwase et al., 2018). Since callus is induced using explants from half-seed or cotyledon in soybean, wounding stress is inevitable (Sairam et al., 2003;Thibaud-Nissen et al., 2003;Jiang et al., 2010). Therefore, it is difficult to divide the callus formation mechanism by wounding and hormones precisely. In this regard, we performed transcriptome analysis on soybean callus induced by wounding and hormones to understand how they affect callus formation. Our data revealed that the complex interaction of pathways by wounding and external hormones activates callus formation (Figure 8).

Recognizing the temporal dynamics:
The distinct stages of callus formation in G. max Based on the analysis of time-series DEGs, which compared each sample to the control (day 0), our research confirmed that 50.31% of all genes revealed transcriptional dynamics at least once during callus formation. Plants tend to activate various metabolic pathways to resist abiotic stress, such as pathogen attack, wounding, UV exposure, high light, and nutrient deficiencies (Akula and Ravishankar, 2011;Yoon et al., 2020). In our study, woundingtreated cotyledon was incubated under abnormal conditions, such as dark and auxin-rich, to induce callus. These processes are necessary for our research to uncover the transcriptional dynamic during callus formation in soybean. So, stress-related TFs such as ERF, MAC, MYB, MYC, and WRKY were highly expressed, as well as hormonal pathways (i.e., auxin, CK, and BR), cell walls, and cell cycles.
Our clustering analysis with DEGs identified genes that change dynamically over time, so we classified dramatic alterations in the expression pattern of genes related to callus formation. Clusters were divided into three classes according to the transcriptional dynamics during callus formation, as shown in Figure 3A: ectopic regulation (clusters 1, 2, 3, and 4), positive regulation (clusters 5 and 6), and negative regulation (clusters 7 and 8). The ectopic regulation class suggests that transcriptional dynamics of the genes involved in specific time points (i.e., days 1, 4, 7, and 28) affect callus formation by responding to hormones and wounding stress. Among them, many genes encoding enzymes for auxin signal transduction and cell proliferation were involved in clusters 1 and 2, which showed to be associated with callus induction. Also, genes related to cell wall modification, which are key regulators in callus development, were differentially expressed before and after day 7, including in cluster 3. In other words, the auxin pathway and cell division for callus induction are initially activated, and cell wall modification, in turn, is rapidly activated for callus proliferation. In contrast, the negative regulation class showed that the genes related to physiological mechanisms are repressed, suggesting that dynamic changes of many cellular events are essential for callus formation, and physiological mechanisms in the normal life cycle are shut down to save the metabolic energy for an abnormal cellular phenomenon.

Influence of endogenous and exogenous hormones on callus formation
For several decades, a multitude of studies have consistently demonstrated that callus development is significantly affected by various phytohormones, notably auxin and cytokinin, which play a pivotal role in triggering the callus formation (Ikeuchi et al., 2013;Ikeuchi et al., 2017;Lee and Seo, 2017;Ikeuchi et al., 2019). Based on our findings, we suggest that the modulation of gene expression during callus formation is governed by the interplay between auxin, cytokinin, and brassinosteroid, working in conjunction with wound signaling mechanisms (Figure 8).
The findings of this study provide confirmation that genes associated with IAM and IPyA were observed to be up-regulated, suggesting the activation of the auxin biosynthesis pathway as a response to wounding. Furthermore, a majority of genes involved in transporter activity exhibited up-regulation, indicating the enhanced expression of genes responsible for auxin transporters, facilitating the influx of external auxin into the cellular environment. Moreover, genes associated with auxin degradation, specifically GH3 and IAMT1, were observed to be up-regulated ( Figure 4A). Based on these findings, it appears highly probable that the level of auxin increased as a result of auxin biosynthesis activation and enhanced transporter activity induced by both wounding and external auxin. Consequently, this leads to the degradation of auxin to uphold auxin homeostasis within the cell.
CK, a crucial hormone along with auxin in callus formation, is synthesized by three biosynthesis-related genes (i.e., IPT, CYP735A, and LOG) in Arabidopsis. In IPT, which is responsible for the initial step in the cytokinin biosynthesis pathway, AtIPT5 and AtIPT7 were induced by auxin, whereas AtIPT1, AtIPT3, AtIPT5, and AtIPT7 were suppressed by CK in shoot meristem (Cheng et al., 2012). In this study, we identified the orthologous genes in soybean, of which only IPT3 exhibited up-regulation ( Figure 4B). Based on our analysis, we deduce that IPT3 is closely associated with the CK response triggered by externally applied CK during callus formation. The up-regulation of CYP73As and LOGs, which play a role in the subsequent step, provides evidence for the activation of CK biosynthesis pathway during callus formation. Within the context of CK signaling, the expression profiles of genes encoding AHK, AHP, ARR, and ESR were observed to be up-regulated. This up-regulation subsequently leads to the activation of CYCD, a key regulator of the cell cycle. These findings strongly suggest that CK synthesis is triggered by wounding signals and, in turn, influences callus formation by promoting cell cycle progression through intricate interactions with signaling factors.
Despite active research, the precise contribution of BR, a phytohormone crucial for development and growth, to callus formation remains uncertain. However, in this study, a significant portion of genes related to BR biosynthesis exhibited up-regulation, shedding light on their potential involvement in the process. Notably, DWF4, a gene already recognized for its robust expression in dividing callus cells (Chung et al., 2011), exhibited an exceptionally elevated level of expression. Moreover, BAK1 and BRLs, which are recognized for their role in regulating various developmental processes, such as cell division, elongation, root growth, vascular differentiation, and resistance to stressors (Oh et al., 2020), exhibited significant up-regulation. This observation Schematic diagram of cellular mechanisms and interactions between wounding and hormones during callus formation. This diagram illustrates the cellular mechanism underlying the interactions between wounding signals and phytohormones, including auxin, cytokinin, and brassinosteroid, during callus formation. The regulatory genes involved in these processes are color-coded: genes regulated by wounding signals are indicated in red, genes related to auxin are shown in blue, genes activated by cytokinin signaling are depicted in green, and genes involved in brassinosteroid signaling are represented in purple. PIN (pin-formed), YUC (YUCCA), AUX (auxin transporter protein 1), ARF (auxin response factor), LBD (lateral organ boundaries domain), WOX (wuschel-like homeobox), PLT (plethora), BRI1 (brassinosteroid insensitive 1), BIN2 (brassinosteroid insensitive 2), BES1 (BRI1-EMS-suppressor 1), ERF (ethylene response factor), ANAC071 (NAC domain-containing protein 71), NAC1 (Nucleus accumbens-associated 1), WIND (Wound induced dedifferentiation), LOG (lonely guy), IPT (isopentenyl transferase gene), AHK (Arabidopsis Histidine Kinase), ARR (Arabidopsis response regulator), WUS (wuschel), ESR (enhancer of shoot regeneration), OPB1 (OBF binding protein 1), DOF2;3 (DNA binding with one finger 2;3).
provides compelling evidence that the accumulation of BR within cells triggers the activation of the cell division process, ultimately leading to callus formation ( Figure 4C). Several studies have explored the signaling pathway, uncovering the role of BIN2 as a positive regulator through the BIN2-ARF-LBD complex, while also noting the activation of BES1/BZR1 during cell division (Wei and Li, 2016;Rane et al., 2017). Consistent with these findings, our research identified distinct expression patterns of these genes throughout callus formation, providing further evidence that injury-related stimuli trigger the synthesis of BR. Consequently, this leads to the stimulation of cell division as a consequential outcome.

Dynamic regulation of cell wall and cell cycle mechanism during callus formation
Several studies have examined transcriptional regulators that play a role in the initiation and progression of callus, such as ARF, LBD, WOX, WIND, and AR2/ERF in different species (Ikeuchi et al., 2013;Ikeuchi et al., 2017;Du et al., 2019;Suo et al., 2021). ERF114, ERF115, ANAC071, and ANAC096 are activated in organ development and regenerative responses via cell wall damage by wounding (Ikeuchi et al., 2017;Zhang et al., 2022). Moreover, ANAC071 is known to directly regulate XTH19 and XTH20, which promote pith cell proliferation (Pitaksaringkarn et al., 2014). The data we obtained demonstrated a substantial up-regulation in the expression pattern of these genes, thereby corroborating previous findings. Consequently, our results strongly indicate that wounding stimuli play a significant regulatory role in cell wall modification during callus formation, mediated by wounding-responsive transcription factors.
Wounding-Induced Networked Developmental Regulator 1 (WIND1), which has homologs WIND2, WIND3, and WIND4, is involved in callus formation and cell division in the apical meristem and is induced by wounding (Iwase et al., 2011). Under wound stress, WIND1 has a direct regulatory effect on ARR1 and ARR12, which participate in cytokinin signaling, as well as on ESR1 and ESR2, which were involved in shoot meristem development (Ikeuchi et al., 2017). In our study, we confirmed that WIND1 and WIND 4 were differentially expressed at day 1, and ARR1/ ARR2, ARR11, ARR12, and ESR1/ESR2 exhibited an increase in the expression profile. Furthermore, ESR2 regulates OBP1, which is the DOF transcriptional regulator. OBP1 plays a role in cell division by targeting the core cell cycle regulator CYCB3;3 and the S phasespecific transcription factor DOF2;3 (Skirycz et al., 2008;Ikeuchi et al., 2013). WUS, a well-known marker gene associated with shoot meristem development, also plays a role in callus formation by governing the fate of meristem cells (Ikeuchi et al., 2013). Analysis of the RNA-seq data unveiled an upsurge in the expression levels of OBP1, DOF2;3, CYCB3;3, and WUS genes. Notably, OBP1 (Glyma.02G062700) and WUS (Glyma.01G166800) exerted significant influence on the transcriptional dynamics of genes within the ectopic regulation cluster. These findings strongly suggest that the cytokinin-dependent pathway involved in shoot meristem development is activated in response to wound stress, potentially playing a role in callus formation by modulating the fate of meristematic cells.
ARFs, TFs that are under the regulation of auxin, have been the subject of extensive research due to their pivotal role in plant development and callus formation (Chen et al., 2011;. ARF7 and ARF19 have been recognized as key players in the process of lateral root development by facilitating the activation of downstream transcription factors, including LBD16, LBD17, LBD18, and LBD29, whose involvement in this pathway is wellestablished (Okushima et al., 2007). Additionally, in the initiation stage of root primordia, the expression of LBD genes was governed by WOX11 and WOX12 (Hu and Xu, 2016). Moreover, it has been observed that LBD18 directly promotes the emergence of lateral roots through the regulation of EXP14 (Lee et al., 2013). Conversely, LBD16 plays a role in root regeneration and lateral root formation, directly activating FAD-BD. This activation occurs through the formation of a heterodimeric complex with bZIP59, thereby facilitating cell wall remodeling (Goh et al., 2012;Liu et al., 2018;Xu et al., 2018). In contrast, LBD29 has been found to regulate PME2 during the process of cell wall modification in callus formation . Our investigation corroborated the up-regulation of LBD and its downstream genes, particularly those implicated in cell wall modification. However, intriguingly, we observed a down-regulation of upstream genes, including ARF7 and ARF19. According to previous research in A. thaliana (Okushima et al., 2007), it was anticipated that the ARF7 and ARF19 genes would have already exhibited expression and initiated the activation of downstream genes prior to 1-day time point examined in our study. Therefore, the noteworthy up-regulation of LBD and its downstream genes observed during callus formation aligns perfectly with this expectation. Consequently, the expression of ARF TFs is induced by auxins, subsequently promoting the expression of LBD genes and facilitating the activation of cell wall modification processes.
WOX5 assumes a critical role as a transcriptional regulator in both lateral root development and callus formation. Its regulation is mediated by LBDs, WOX11, and WOX12, while it acts as a regulator of PLT1 and PLT2 (Ding and Friml, 2010;Sugimoto et al., 2010;Hu and Xu, 2016). Moreover, WOX7 and WOX14, which interact with WOX5, are recognized for their pivotal roles in pluripotency acquisition during callus formation, as well as in lateral root development and vascular cell differentiation (Kong et al., 2016;Denis et al., 2017). Our investigation made a noteworthy observation of substantial up-regulation of WOX5 and WOX7 over a period of seven days. However, it is important to note that orthologous genes of WOX11 and WOX12 were not identified in G. max. Interestingly, PLT3, PLT5, and PLT7 were identified as key regulators of root development, playing a crucial role in de novo shoot regeneration and the acquisition of cell pluripotency through the regulation of PLT1 and PLT2 activity (Kareem et al., 2015;Bustillo-Avendaño et al., 2018). Furthermore, another study reported an up-regulation of PLT family genes in A. thaliana as early as one hour after wounding, which triggers callus formation (Ikeuchi et al., 2017). Through our research, we noted a decline in the expression levels of PLT3, PLT5, and PLT7 genes starting from day 0, while the expression of PLT1 and PLT2 genes showed an increase. These findings indicate that our RNA-Seq analysis was conducted after day 1, preventing us from confirming the rapid upregulation of PLT3, PLT5, and PLT7 within one hour. However, we were able to observe the expression of downstream genes, such as PLT1 and PLT2. These results imply that WOX5 plays a critical role in regulating PLT1 and PLT2, and the interplay between auxin and wounding stimulates the mechanism of root meristem development. Ultimately, this contributes to callus formation by influencing meristematic fate.

Conclusion
In this study, our objective was to explore the transcriptional regulation of callus formation in soybean under the influence of wounding signals and phytohormones. To gain deeper insights into the intricate mechanisms involved in this process, we performed a thorough transcriptome analysis of G. max. Employing a timeseries gene expression profiling approach, we successfully identified key genes associated with transcription factors, biosynthesis, transporters, and signaling pathways related to phytohormones. These findings collectively contribute to a comprehensive understanding of the underlying mechanisms governing callus formation. Notably, this study goes beyond elucidating individual components and explores the intricate interactions among major phytohormones involved in callus formation, including auxin and cytokinin (as previously explained), and the newly implicated brassinosteroid. Such a broader investigation will enhance our understanding of this complex phenomenon. Our results indicate that the coordinated interplay of wounding, auxin, cytokinin, and brassinosteroid signaling pathways triggers the activation of genes involved in determining the fate of meristematic cells. This, in turn, initiates processes such as cell wall modification and cell cycle reentry, ultimately leading to the formation of callus. By unveiling these intricate interactions, our study sheds light on the complex regulatory network orchestrating callus formation in soybean.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/ bioproject/PRJNA949161, PRJNA949161.

Funding
The research described in this paper was supported by multiple sources of grants, including the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (grant numbers; 2020R1A2C201005714 and 2020R1A6A1A03047729) and Green Fusion Technology Program funded by the Ministry of Environment, Republic of Korea.