Transcriptome Profiling Reveals Differential Gene Expression during the Process of Microtuber Formation in Pinellia ternata

Using petiole material as explants and directly inducing the formation of microtubers without going through the callus stage is an essential way to rapidly expand scarce medical plants such as Pinellia ternata. However, the early molecular mechanism underlying the formation of the microtuber is largely elusive. Here, we conducted cytology and dynamic transcriptome analyses of inchoate microtubers in Pinellia explants and identified 1092 differentially expressed genes after their cultivation in vitro for 0, 5, and 15 days. Compared with 0 day, the number and size of the microtuber cells were larger at 5 and 15 days of culture. Detailed categorization revealed that the differentially expressed genes were mainly related to responses to stimulus, biological regulation, organelles, membranes, transcription factor activity, and protein binding. Further analysis revealed that the microtuber at different incubation days exhibited quite a difference in both hormone signaling pathway transduction and the regulation pattern of transcription factors. Therefore, this study contributes to a better understanding of the early molecular regulation during the formation of the microtuber and provides new insights for the study of the rapid expansion of P. ternata and other medical plants.


Introduction
Pinellia ternata (Thunb.) Breit is a valuable medical plant that belongs to the Araceae family. It is native to the eastern part of Asia and now mainly distributed in China [1,2]. The tuber of P. ternata is the central medicinal part, which has been frequently used in traditional Chinese medicine for thousands of years [3][4][5]. Modern pharmacological research has demonstrated that P. ternata tubers have anti-depressant, wound-healing, anti-coughing, and anti-vomiting functions [6,7]. In addition, alkaloid compounds, including guanosine, inosine and ephedrine, have been recognized as the main active ingredients and are believed to exert anticancer effects [8].
In recent years, with climatic and environmental changes and indiscriminate logging all intensifying, the natural sources of P. ternata have gradually decreased, and artificial cultivation has become mainstream. The natural seed-setting rate of P. ternata is very limited as it is mostly cultivated using bulbs and explants for asexual reproduction [9]. However, the former is less efficient at reproduction, and annual production is far too inadequate to meet the market demand [10]. With the advancement of plant tissue culture technology, direct induction of microtuber formation using petiole material as explants without going through the callus stage has become an important way to ensure the rapid expansion of P. ternata [11].
Microtubers can be produced in vitro from explants as a single node cutting or from petioles in tuber-inducing medium in small stationary-container bioreactors and fermenters, 2 of 14 and different concentrations of plant growth regulators have different effects on the formation of microtubers in P. ternata [12,13]. It has been found that in vitro, P. ternata cell cultures produce similar amounts of alkaloid compounds to field-grown tubers [14]. Effective microtuber induction could be achieved by subculturing the morphologically expanded lower portion of suspension-cultured petioles into the solid MS medium without the addition of additional plant growth regulators. And an additional 100 µM salicylic acid (SA) increased the concentration of total alkaloid compounds, guanosine, inosine, and ephedrine by 2.5-, 2.1-, 2.8-, and 3.1-fold, respectively, suggesting that the induction of in vitro microtubers in such cultures is a promising approach for large-scale alkaloid production [15]. After making microtubers induced in the test tube into artificial seeds for cultivation, the tuber yield was twice that of cultivated seedlings, and the content of the active ingredient alkaloids had increased by 1.6-fold [16]. Hence, microtuber technology has significant advantages in tuber seed production, germplasm conservation, breeding programs, and research.
With advanced sequencing technology, a large amount of data resources and research methods have been applied to study P. ternata. Lu et al. (2013) [17] used EST library sequencing to identify the expression of candidate genes related to tuber development and high temperature response. Another study identified a number of potential genes associated with delayed sprout-tumble formation and tuber growth by full-length transcriptome sequencing of shade-induced P. ternata [18]. However, the early molecular mechanisms of P. ternata microtuber formation as well as the hormones and the transcription factors that are involved in the formation process are largely unknown. In this study, we performed dynamic transcriptome analysis of P. ternata microtuber cultured for different time intervals in order to analyze differentially expressed genes (DEGs) and to mine and screen candidate genes associated with microtuber formation. Our primary aim was to gain an overall understanding of the genes that regulate microtuber formation, which could benefit the rapid expansion and varietal improvement of P. ternata.

Phenotypic Comparison of P. ternata Microtubers with Different Cultivation Times
Early formation of the microtuber occurs in P. ternata explants, with gradual expansion occurring at one end of the petiole as time progresses. As shown in Figure 1a, obvious signs of expansion were observed after 5 d of culture in Murashige and Skoog (MS) medium, and when cultured for 15 d, the average diameter of the expansion site reached 3 mm (Figure 1a). To further analyze the differences in microtubers, fresh samples (microtubers cultured for 0, 5, and 15 d) were stained with Safranin O-Fast Green, and then sectioning was performed. The changes in microtubers at different cultivation times were observed by light microscopy. Compared with 0 d, microtuber cells cultured for 5 d showed dense protoplasm. The number and size of cells increased when cultured for 15 d, and black starch granules were observed inside the cells (Figure 1b). These results suggest that the microtubers in P. ternata explants undergo progressive expansion and cellular changes over time. The increase in cell number, size, and the presence of starch granules indicate the development and maturation of the microtubers.

Transcriptome Sequencing of P. ternata Microtubers with Different Cultivation Times
To further elucidate the early molecular mechanisms of microtuber formation in P. ternata, transcriptome analysis was performed on microtubers cultured for 0 d, 5 d, and 15 d. There were nine samples in total, BX0 (culture for 0 d), BX5 (culture for 5 d), and BX15 (culture for 15 d), with three replicates for each group. A total of 996,835 transcripts was extracted from the cDNA libraries, which were further filtered and assembled into 363,624 unigene clusters with an impressive N50 value of 610 bp (Table S1). After removing low-quality reads, the total reads and total bases produced by each library ranged from 41.37 million to 58.22 million and 59.17 million to 87.39 million. The contents of Q30 and GC were 93.48-94.88% and 52.46-53.82%, respectively, indicating the high quality of the transcriptome sequencing data (Table S2). The expression level of each gene was normalized by the transcripts-per-million-reads (TPM) method. A clustering heatmap ( Figure 2a) and principal component analysis (PCA) analysis plot (Figure 2b) revealed that the gene expression data of the samples were highly reproducible between the three biological replicates of each sample, and the trends of transcriptome expression levels were similar in microtubers cultured for 5 d and 15 d, while the trends were the opposite for the 0 d samples.

Transcriptome Sequencing of P. ternata Microtubers with Different Cultivation Times
To further elucidate the early molecular mechanisms of microtuber formation in P. ternata, transcriptome analysis was performed on microtubers cultured for 0 d, 5 d, and 15 d. There were nine samples in total, BX0 (culture for 0 d), BX5 (culture for 5 d), and BX15 (culture for 15 d), with three replicates for each group. A total of 996,835 transcripts was extracted from the cDNA libraries, which were further filtered and assembled into 363,624 unigene clusters with an impressive N50 value of 610 bp (Table S1). After removing low-quality reads, the total reads and total bases produced by each library ranged from 41.37 million to 58.22 million and 59.17 million to 87.39 million. The contents of Q30 and GC were 93.48-94.88% and 52.46-53.82%, respectively, indicating the high quality of the transcriptome sequencing data (Table S2). The expression level of each gene was normalized by the transcripts-per-million-reads (TPM) method. A clustering heatmap ( Figure  2a) and principal component analysis (PCA) analysis plot (Figure 2b) revealed that the gene expression data of the samples were highly reproducible between the three biological replicates of each sample, and the trends of transcriptome expression levels were similar in microtubers cultured for 5 d and 15 d, while the trends were the opposite for the 0 d samples.

Transcriptome Sequencing of P. ternata Microtubers with Different Cultivation Times
To further elucidate the early molecular mechanisms of microtuber formation in P. ternata, transcriptome analysis was performed on microtubers cultured for 0 d, 5 d, and 15 d. There were nine samples in total, BX0 (culture for 0 d), BX5 (culture for 5 d), and BX15 (culture for 15 d), with three replicates for each group. A total of 996,835 transcripts was extracted from the cDNA libraries, which were further filtered and assembled into 363,624 unigene clusters with an impressive N50 value of 610 bp (Table S1). After removing low-quality reads, the total reads and total bases produced by each library ranged from 41.37 million to 58.22 million and 59.17 million to 87.39 million. The contents of Q30 and GC were 93.48-94.88% and 52.46-53.82%, respectively, indicating the high quality of the transcriptome sequencing data (Table S2). The expression level of each gene was normalized by the transcripts-per-million-reads (TPM) method. A clustering heatmap ( Figure  2a) and principal component analysis (PCA) analysis plot (Figure 2b) revealed that the gene expression data of the samples were highly reproducible between the three biological replicates of each sample, and the trends of transcriptome expression levels were similar in microtubers cultured for 5 d and 15 d, while the trends were the opposite for the 0 d samples.

Analysis of DEGs in BX5 vs. BX0, BX15 vs. BX0, and BX15 vs. BX5 Comparisons
Each sample was compared with the other two to identify the DEGs (q-value ≤ 0.05 and |FoldChange| ≥ 2), and 37,990 DEGs were obtained from these three comparisons (BX5 vs. BX0, BX15 vs. BX0, and BX15 vs. BX5). From these pairwise comparisons, we  The Venn statistics of the comparative analysis between the groups identified 1092 DEGs that were common to microtuber formation ( Figure 3b). This further indicates that DEGs associated with microtuber formation are abundantly expressed throughout the culture period and that these genes are significantly up-or down-regulated in expression during the microtuber formation phase, suggesting that the mechanisms underlying their early formation may be complex and variable. The specific expression levels and details of these DEGs are shown in Table S3.

Functional Classification of DEGs
A gene ontology (GO) enrichment analysis was performed in order to determine the function of the identified DEGs. The 17,028 DEGs were enriched in 3627 GO terms in the biological process categories (Table S4). Eleven GO terms included more than 1000 DEGs. Biological process (GO:0008150) and cellular process (GO:0009987), with 2514 and 2019 DEGs, respectively, were the two most dominant terms in the biological process category.
Among the DEGs identified among the pairwise comparisons, 15 common GO terms with a q-value of <0.05 were identified in the biological process category ( Figure 4). As depicted in Figure 4, the significantly enriched biological processes (BPs) encompassed response to stimulus (GO:0050896), biological regulation (GO:0065007), and regulation of biological processes (GO:0050789), while the most significantly enriched cellular components (CCs) belonged to the organelle (GO:0043226), membrane (GO:0016020), and membrane parts (GO:0044425). In addition, transporter activity (GO:0005215), transcription factor activity, protein binding (GO:0000988), and signal transducer activity (GO:0004871) were significantly enriched within the molecular function (MF) categories. The Venn statistics of the comparative analysis between the groups identified 1092 DEGs that were common to microtuber formation ( Figure 3b). This further indicates that DEGs associated with microtuber formation are abundantly expressed throughout the culture period and that these genes are significantly up-or down-regulated in expression during the microtuber formation phase, suggesting that the mechanisms underlying their early formation may be complex and variable. The specific expression levels and details of these DEGs are shown in Table S3.

Functional Classification of DEGs
A gene ontology (GO) enrichment analysis was performed in order to determine the function of the identified DEGs. The 17,028 DEGs were enriched in 3627 GO terms in the biological process categories (Table S4). Eleven GO terms included more than 1000 DEGs. Biological process (GO:0008150) and cellular process (GO:0009987), with 2514 and 2019 DEGs, respectively, were the two most dominant terms in the biological process category.
Among the DEGs identified among the pairwise comparisons, 15 common GO terms with a q-value of <0.05 were identified in the biological process category ( Figure 4). As depicted in Figure 4, the significantly enriched biological processes (BPs) encompassed response to stimulus (GO:0050896), biological regulation (GO:0065007), and regulation of biological processes (GO:0050789), while the most significantly enriched cellular components (CCs) belonged to the organelle (GO:0043226), membrane (GO:0016020), and membrane parts (GO:0044425). In addition, transporter activity (GO:0005215), transcription factor activity, protein binding (GO:0000988), and signal transducer activity (GO:0004871) were significantly enriched within the molecular function (MF) categories. To further investigate the biological functions of these DEGs, pathway-based analysis was conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Cell growth and death, transport and catabolism, signal transduction, replication and repair, and translation were the topmost enriched pathways in the microtuber formation of P. ternata ( Figure 5).

Analysis of DEGs Encoding Transcription Factors (TFs)
Next, the DEGs encoding TFs were analyzed. A total of 761 DEGs encoding TFs were identified in P. ternata tissues that responded to microtuber formation, and these TFs belonged to 31   To further investigate the biological functions of these DEGs, pathway-based analysis was conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Cell growth and death, transport and catabolism, signal transduction, replication and repair, and translation were the topmost enriched pathways in the microtuber formation of P. ternata ( Figure 5).  To further investigate the biological functions of these DEGs, pathway-based analysis was conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Cell growth and death, transport and catabolism, signal transduction, replication and repair, and translation were the topmost enriched pathways in the microtuber formation of P. ternata ( Figure 5).

Analysis of DEGs Encoding Transcription Factors (TFs)
Next, the DEGs encoding TFs were analyzed. A total of 761 DEGs encoding TFs were identified in P. ternata tissues that responded to microtuber formation, and these TFs belonged to 31 TF families. Most of the identified DEGs encoded members of the NAC, ERF, WRKY, bHLH, and MYB-related TF families (Figure 6a), and 22 TF families included more

Analysis of DEGs Encoding Transcription Factors (TFs)
Next, the DEGs encoding TFs were analyzed. A total of 761 DEGs encoding TFs were identified in P. ternata tissues that responded to microtuber formation, and these TFs belonged to 31 TF families. Most of the identified DEGs encoded members of the NAC, ERF, WRKY, bHLH, and MYB-related TF families (Figure 6a), and 22 TF families included more than 100 differentially expressed TFs. The NAC family, with 89 DEGs, was the largest TF family responding to microtuber formation in P. ternata, including 38, 40, and 11 DEGs in the three comparisons, respectively. A total of 60 DEGs belonging to the ERF family were identified, including 20, 29, and 11 DEGs in the three comparisons, respectively. Furthermore, a total of 56 WRKY TFs, 51 bHLH TFs, and 50 MYB-related TFs were found.  We analyzed the TFs common to the three comparison groups, and a total of 27 genes were identified. Compared with BX0, all of the common differentially expressed NACs and ERFs were down-regulated in BX5 and BX15; conversely, MYB-related and HB-other were typically up-regulated (Figure 6b). Transcription factor genes have different expression patterns in the P. ternata microtuber at different days of culture, suggesting that their formation processes are controlled by a combination of factors.

Analysis of DEGs Associated with Plant Hormones
To identify crucial genes involved in important pathways, DEGs enriched in the auxin (IAA), cytokinine (CK), gibberellin (GA), abscisic acid (ABA), ethylene (ET), and jasmonic acid (JA) signaling pathways were further identified from the KEGG analysis, and about 70 DEGs involved in hormone biosynthesis and signaling transduction pathways were identified. In the IAA signaling pathway, a total of five DEGs were associated with the AUX/IAA signaling pathway, of which two were up-regulated and three were down-regulated ( Figure  7a). This suggests that the transmission of growth hormone signals during microtuber formation may be multi-pathway. For DEGs associated with the CK signaling pathway, most genes were down-regulated; however, for DEGs associated with the GA signaling pathway, most were up-regulated (Figure 7b,c). These results indicate that GA and CK play distinct roles in microtuber formation. In the ABA signaling pathway, DEGs associated with the PYR/PYL, SnRK2, and ABF signaling pathways were mostly up-regulated, while two DEGs associated with PP2C were significantly down-regulated (Figure 7d). This implies that ABA has a crucial role in microtuber formation. Similarly, we also found that DEGs associated with the ET and JA signaling pathways were mostly down-regulated, suggesting they may play similar roles in microtuber formation (Figure 7e,f). We analyzed the TFs common to the three comparison groups, and a total of 27 genes were identified. Compared with BX0, all of the common differentially expressed NACs and ERFs were down-regulated in BX5 and BX15; conversely, MYB-related and HB-other were typically up-regulated (Figure 6b). Transcription factor genes have different expression patterns in the P. ternata microtuber at different days of culture, suggesting that their formation processes are controlled by a combination of factors.

Analysis of DEGs Associated with Plant Hormones
To identify crucial genes involved in important pathways, DEGs enriched in the auxin (IAA), cytokinine (CK), gibberellin (GA), abscisic acid (ABA), ethylene (ET), and jasmonic acid (JA) signaling pathways were further identified from the KEGG analysis, and about 70 DEGs involved in hormone biosynthesis and signaling transduction pathways were identified. In the IAA signaling pathway, a total of five DEGs were associated with the AUX/IAA signaling pathway, of which two were up-regulated and three were downregulated (Figure 7a). This suggests that the transmission of growth hormone signals during microtuber formation may be multi-pathway. For DEGs associated with the CK signaling pathway, most genes were down-regulated; however, for DEGs associated with the GA signaling pathway, most were up-regulated (Figure 7b,c). These results indicate that GA and CK play distinct roles in microtuber formation. In the ABA signaling pathway, DEGs associated with the PYR/PYL, SnRK2, and ABF signaling pathways were mostly up-regulated, while two DEGs associated with PP2C were significantly down-regulated (Figure 7d). This implies that ABA has a crucial role in microtuber formation. Similarly, we also found that DEGs associated with the ET and JA signaling pathways were mostly downregulated, suggesting they may play similar roles in microtuber formation (Figure 7e,f)

Validation of RNA-Seq Analysis by Quantitative Real-Time PCR (qRT-PCR)
To validate the reliability of the gene expression data obtained by the RNA-seq analysis in the microtubers of P. ternata, eight DEGs that play different roles in the formation of the microtuber were selected from the three different culture-day samples for qRT-PCR analysis (Table S5). The ratio of expression levels found between the samples using qRT-PCR was compared with the ratio of expression as measured by RNA-Seq. A significant correlation (r 2 = 0.8279, n = 24, Figure 8) was observed between the RNA-Seq and qRT-PCR data, which confirmed the authenticity of the DEGs in this study. Thus, these comparisons of data from qRT-PCR and RNA-seq analyses of microtubers fully validated the findings from our transcriptome study.

Validation of RNA-Seq Analysis by Quantitative Real-Time PCR (qRT-PCR)
To validate the reliability of the gene expression data obtained by the RNA-seq analysis in the microtubers of P. ternata, eight DEGs that play different roles in the formation of the microtuber were selected from the three different culture-day samples for qRT-PCR analysis (Table S5). The ratio of expression levels found between the samples using qRT-PCR was compared with the ratio of expression as measured by RNA-Seq. A significant correlation (r 2 = 0.8279, n = 24, Figure 8) was observed between the RNA-Seq and qRT-PCR data, which confirmed the authenticity of the DEGs in this study. Thus, these comparisons of data from qRT-PCR and RNA-seq analyses of microtubers fully validated the findings from our transcriptome study.

Discussion
As a valuable medicinal plant, the tuber of P. ternata serves as a crucial medicinal herb resource [19]. Recently, due to climate change and increased indiscriminate logging, wild P. ternata resources have declined, making artificial cultivation the primary method for obtaining P. ternata herbs. An effective approach for the rapid propagation of P. ternata involves using the petiole as an explant, inducing callus tissue through dedifferentiation, and then differentiating and culturing to form microtubers as propagules [20,21]. However, the early molecular regulatory mechanisms of microtuber formation remain unclear, limiting genetic improvements.

Transcriptomic Analysis Yields 1092 DEGs Associated with Microtuber Formation
The identification of genes expressed during microtuber formation at the transcriptome level is an essential approach to understanding the regulatory circuits that control the microtuber formation process. This study provided cDNA sequence data from three microtuber formation stages. A total of 1092 overlapping genes that were differentially expressed in P. ternata explants from the three different comparisons over 15 days were focused on, and the early molecular mechanisms of microtuber formation were summarized. GO enrichment analysis revealed various biological processes, cellular components, and molecular activities relevant to P. ternata microtuber formation. It is reported that low auxin and cytokinin concentrations stimulate P. ternata microtuber polar growth, utilizing the petiole as an explant [12]. These include the response to stimulus, biological regulation, organelle function, membrane processes, as well as transcription factor activity and

Discussion
As a valuable medicinal plant, the tuber of P. ternata serves as a crucial medicinal herb resource [19]. Recently, due to climate change and increased indiscriminate logging, wild P. ternata resources have declined, making artificial cultivation the primary method for obtaining P. ternata herbs. An effective approach for the rapid propagation of P. ternata involves using the petiole as an explant, inducing callus tissue through dedifferentiation, and then differentiating and culturing to form microtubers as propagules [20,21]. However, the early molecular regulatory mechanisms of microtuber formation remain unclear, limiting genetic improvements.

Transcriptomic Analysis Yields 1092 DEGs Associated with Microtuber Formation
The identification of genes expressed during microtuber formation at the transcriptome level is an essential approach to understanding the regulatory circuits that control the microtuber formation process. This study provided cDNA sequence data from three microtuber formation stages. A total of 1092 overlapping genes that were differentially expressed in P. ternata explants from the three different comparisons over 15 days were focused on, and the early molecular mechanisms of microtuber formation were summarized. GO enrichment analysis revealed various biological processes, cellular components, and molecular activities relevant to P. ternata microtuber formation. It is reported that low auxin and cytokinin concentrations stimulate P. ternata microtuber polar growth, utilizing the petiole as an explant [12]. These include the response to stimulus, biological regulation, organelle function, membrane processes, as well as transcription factor activity and protein binding. Our findings provide insights into the molecular events underlying microtuber formation and shed light on the functional roles of specific genes in this process. KEGG pathway analysis showed that KEGG-enriched pathways mainly included cell growth and death, transport and catabolism, signal transduction, replication and repair, and translation. These results show that microtuber formation is a complex phenomenon involving many biological processes.

DEGs Encoding Transcription Factors
The highly complex biological process of tuberization can be influenced by both environmental and endogenous factors, such as photoperiod, plant growth regulators, and plantlet growth stage [22]. Understanding the regulatory mechanism of tuberization is essential for developing strategies to improve tuber yield and quality. As a result, it is necessary to conduct a more comprehensive molecular biological analysis of microtuber formation from a new perspective. Transcription factors regulate almost all aspects of plant growth and development and can orchestrate regulatory networks to improve resistance to abiotic and biotic stresses in plants [23]. Plant TF families such WRKY, NAC, AP2/ERF, bZIP, and MYB regulate numerous key biological processes throughout growth and development [24][25][26][27][28]. In this study, a total of 31 TF families containing 761 differentially expressed TFs were identified by paired comparison. Most of the differentially expressed TFs belonged to the NAC, ERF, WRKY, bHLH, and MYB TF families ( Figure 6). The NAC TF family, with 89 DEGs, is the largest class that plays the largest role in the formation of the P. ternata microtuber, including 38, 40, and 11 DEGs in the three comparisons, respectively. A previous study indicated that ANAC043, ANAC066, and ANAC012 play key regulatory roles in secondary cell wall biosynthesis in Arabidopsis [29]. The identification of a large number of NAC TFs as DEGs suggests that NAC proteins may play a widespread role in the early formation of microtubers. ERF TFs are essential for floral growth, spikelet meristem destiny, root initiation, and seed/fruit development [30]. We found 20, 29, and 11 differentially expressed ERFs in each of the three comparisons during the early formation of the microtuber. Notably, AtERF38, an AP2/ERF family TF in Arabidopsis, was found to be involved in the regulation of several genes related to secondary wall thickening [31]. AtWRKY44 (also known as TTG2), bZIP29, and AtMYB59 were all shown to be involved in plant cell development and root growth [32][33][34]. In this study, we identified 56, 51, and 50 WRKY, bZIP, and MYB TFs involved in microtuber formation. The findings contribute to a better understanding of the early formation of microtubers and the roles of specific TF families in this process.

Effect of Plant Hormones on Microtuber Formation
IAA, CK, GA, ABA, ET, and JA are small-molecule plant hormones that affect many phases of plant development [35]. In this study, we observed a significant increase in cell number and size with the expansion of P. ternata microtubers (Figure 1), likely due to a combination of hormones. IAA is a major regulator of tuber formation and increases during early tuberization to drive cell proliferation and storage tuber formation [36]. Our transcriptome analysis revealed a total of five DEGs associated with the AUX/IAA signaling pathway, of which two were up-regulated and three were down-regulated. In Arabidopsis, IAA5 and IAA19 were rapidly and highly induced during callus induction, and they both belong to the AUX/IAA signaling pathway [37], suggesting that expression of similar genes may be highly induced during microtuber formation in P. ternata. The ratio of auxin to cytokinin has been suggested to be critical in determining cell fate in in vitro systems in typical plant tissue cultures [38]. Our analysis showed six DEGs involved in the CK signaling pathway, suggesting that CK may play a coordinating role for microtuber formation. GA can regulate cell elongation, which is essential for biological yield [39]. We found six transcription factor DEGs in the GA signaling pathway in this work, demonstrating that GA is essential to microtuber formation. PYR/PYL, a family of ABA receptor proteins, mediates plant immunity, seed germination, and seedling establishment [40], and our study has shown that four genes are involved in encoding PYR/PYL-type proteins, suggesting that ABA also plays a regulatory role in microtuber formation. ETH is essential for plant growth and development, including seed germination, leaf senescence, fruit ripening, and responses to environmental stresses [41]. We found genes encoding the relevant proteins on many branches of the ETH pathway, showing that the ETH pathway is crucial for microtuber formation. The JAZ protein and MYC2 are key components in the crosstalk between JA signaling and other plant hormone signaling pathways that regulates plant growth and stress responses [42]. Our investigation identified six DEGs involved in JAC and MUC2 protein coding, indicating that JA also contributes to the formation of microtubers. Although this study did not demonstrate the existence of multiple hormone signaling pathways under these experimental conditions, based on previous research and our transcriptome analysis, it can be inferred that P. ternata has numerous hormonal pathways associated with microtuber formation.
Taken together, our study used transcriptomic analysis to identify genes expressed during the microtuber formation of P. ternata, identifying 1092 DEGs and 31 transcription factor families involved in the process. Further analysis revealed DEGs associated with various hormone signaling pathways, highlighting the importance of plant hormones in the early formation of microtubers and providing insights into their complex regulatory mechanisms.

Plant Materials and Growth
Potted P. ternata plants were cultivated in a greenhouse with a 16 h/8 h light/dark photoperiod at 25 • C. The mature bulbil of healthy plants was then cut and rinsed under running water, disinfected on the surface of an ultra-clean bench with 70% alcohol for 30 s, disinfected with 0.1% mercury solution for 8-10 min, rinsed 6-7 times with sterile water, and then suspension-cultured in MS medium containing 0.2 mg/L 2,4-dichlorophenoxyacetic acid (2,4-D) and 1 mg/L kinetin (KT) until complete plants were formed. When the plants grew to a height of 4-5 cm, they were cut off at about 1.5 cm of the petiole and placed in MS medium containing 0.5 mg/L 6-benzylaminopurine (6-BA), 0.5 mg/L indolacetic acid (IAA), 5% sucrose to induce the formation of the microtuber.

Microscopic Structure Analysis of P. ternata Microtubers
Histological analysis was performed on microtubers of P. ternata sampled from three different cultivation times. Following an 8-hour fixation in FAA fixative (composed of formalin, acetic acid, and 70% alcohol in a ratio of 1:1:18), the sections of samples were rehydrated in BioDewax and put into the 0.1% safranin O staining solution for 3-8 s. Afterward, the sections were stained with 0.1% Fast Green staining solution for 6-20 s after washing away the excess dye with 1% acetic acid [43]. Imaging was performed under a Motic PA53 FS6 microscope (Motic, Hong Kong SAR, China).

cDNA Library Construction and Sequencing
For total RNA extraction, samples from the swollen part of the petiole were collected and frozen in liquid nitrogen at 0, 5, and 15 d after inducing microtuber formation (0.5 cm length of petiole taken close to the side of the MS medium at 0 d); samples from each cultivation time contained three biological replicates. The methods for total RNA extraction and RNA-seq library construction were the same as previously described [44]. Total RNA from the spike tissues was extracted using the TRIzol kit (Invitrogen, Carlsbad, CA, USA). The RNA was examined for quality and quantified using an agarose gel and the NanoDrop 2500 (Thermo Fisher Scientific, Wilmington, DE, USA). The mRNA was subjected to enrichment through the use of magnetic beads featuring Oligo (dT), with subsequent fragmentation through the application of fragmentation buffer to random fragments. The initial cDNA strand was produced via the utilization of small fragments as templates and arbitrary hexameric primers, which was further followed by the creation of the secondary cDNA strand using dNTPs, RNase H, and DNA polymerase I. The double-stranded complementary DNA (cDNA) was purified through the use of magnetic beads, followed by an end-repair process that involved the addition of a single nucleotide A (adenine) to the 3 ends and subsequent ligation with sequencing adaptors. The appropriate fragments were selectively chosen as templates and enriched through PCR amplification. The quality and concentration of the RNA were evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Equal amounts of RNA from each sample were sequenced on an Illumina HiSeq 2500 sequencer (Illumina, San Diego, CA, USA) at the Sangon Biotech company (Shanghai, China) to construct the strand-specific RNA-seq libraries.

Data Processing and De Novo Assembly
The original data produced by the sequencer were defined as raw reads. Raw reads from each sample were filtered by removing reads containing more than 10% of unknown nucleotides, low-quality reads that contained more than 40% of low-quality (Q value ≤ 20) bases, and adapter sequences according to the Illumina adapter list. The remaining highquality reads were assembled de novo using the Trinity software (version 2.14.0) [45].

Functional Annotation of Unigenes and Differential Expression Genes Analysis
The assembled unigenes were aligned with a non-redundant protein sequence (Nr), cluster of orthologous groups of proteins (COG), gene ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG), and SwissProt databases, as well as highly similar proteins to annotate the unigenes IDs. Gene expression level read-count data was analyzed using DESeq2 for differential genes [46]. Genes with a false discovery rate (FDR) of <0.05 and |log2FC| > 1 were considered significant DEGs. Using Blast2go software (version 2.3.5, http://www.blast2go.de/, 16 July 2023), the significantly different unigenes were compared with the KEGG and GO databases for annotation and enrichment analysis. The results of functional annotation and pathway enrichment were combined to screen for candidate genes related to microtuber formation.

Validation of RNA-seq Data by qRT-PCR
Eight DEGs were selected for validation by qRT-PCR. Single-stranded cDNA was obtained using the Evo MMLV RT Premix for qRT-PCR (Accurate Biotechnology Co., Ltd., Changsha, China) according to the manufacturer's protocol. Gene-specific primers were designed using Primer3Plus (http://www.primer3plus.com/, accessed on 6 December 2022). FastStart Essential DNA Green Master (Roche, Basel, Switzerland) was used for the PCRs. The reaction system and procedure are described in a previous paper [47]. The Pt18SrRNA gene was used as an internal control for normalization [48], and three technical replicates of each cDNA sample were analyzed. The primer sequences for qRT-PCR are listed in Table S4. The 2 −∆∆CT method was used to calculate the relative expression level of each gene [49].

Conclusions
We analyzed RNA-seq to provide a global overview of the differences in the transcriptomes of microtubers from P. ternata at 0 d, 5 d, and 15 d of formation. A total of 1092 common DEGs were identified in microtubers at the three different culture times. The analysis of DEGs not only revealed that a large number of genes encoding TFs were highly induced during this process, but also identified some DEGs associated with hormones during microtuber formation. This will provide new insights for studying microtuber formation and rapid propagation in P. ternata.