De novo and comparative transcriptomic analysis explain morphological differences in Panax notoginseng taproots

Panax notoginseng (Burk.) F. H. Chen (PN) belonging to the genus Panax of family Araliaceae is widely used in traditional Chinese medicine to treat various diseases. PN taproot, as the most vital organ for the accumulation of bioactive components, presents a variable morphology (oval or long), even within the same environment. However, no related studies have yet explained the molecular mechanism of phenotypic differences. To investigate the cause of differences in the taproot phenotype, de novo and comparative transcriptomic analysis on PN taproot was performed. A total of 133,730,886 and 114,761,595 paired-end clean reads were obtained based on high-throughput sequencing from oval and long taproot samples, respectively. 121,955 unigenes with contig N50 = 1,774 bp were generated by using the de novo assembly transcriptome, 63,133 annotations were obtained with the BLAST. And then, 42 genes belong to class III peroxidase (PRX) gene family, 8 genes belong to L-Ascorbate peroxidase (APX) gene family, and 55 genes belong to a series of mitogen-activated protein kinase (MAPK) gene family were identified based on integrated annotation results. Differentially expressed genes analysis indicated substantial up-regulation of PnAPX3 and PnPRX45, which are related to reactive oxygen species metabolism, and the PnMPK3 gene, which is related to cell proliferation and plant root development, in long taproots compared with that in oval taproots. Furthermore, the determination results of real-time quantitative PCR, enzyme activity, and H2O2 content verified transcriptomic analysis results. These results collectively demonstrate that reactive oxygen species (ROS) metabolism and the PnMPK3 gene may play vital roles in regulating the taproot phenotype of PN. This study provides further insights into the genetic mechanisms of phenotypic differences in other species of the genus Panax.


Introduction
Panax notoginseng (Burk.) F. H. Chen (PN) is a perennial herb belonging to the Araliaceae family. PN taproot, as the most vital organ for the accumulation of bioactive components, is a raw material of many famous patented Chinese medicines, including Yunnan Baiyao, Xuesaitong, and Compound Danshen Dripping Pills [1,2]. The main bioactive components isolated from PN taproots include ginseng saponins, notoginsenoside, and dencichine [2,3]. Extensive pharmacology studies have shown that PN prevents cardio-and cerebrovascular disease, has anti-inflammatory effects, and aids immune regulation [4,5]. Xuesaitong, produced from the total saponins extracted from PN Open Access *Correspondence: liuyuan513@kust.edu.cn; sanqi37@vip.sina.com 7 Sanqi Research Institute of Yunnan Province, Kunming 650000, China Full list of author information is available at the end of the article taproots, is in broad clinical use for the prevention and treatment of hyperlipidemia, coronary heart disease, stroke sequelae, and other chronic diseases, and is particularly favored by patients. In general, the growth rate of PN is relatively slow in the wild, and even with good field management it takes at least three years to grow and accumulate sufficient bioactive components in its taproots before they can be harvested and used as a medicinal material. However, in the process of PN cultivation, there has been a perplexing phenomenon wherein there are obvious differences in the root phenotype during growth, even when many PN plants are grown in the same environment. Specifically, some of them present long strip-like ginseng (long taproot of PN: LPN), whereas others are oval-like (oval taproot of PN: OPN); they are called "Luobo qi" or "Chang qi" and "Geda qi" or "Tuan qi, " respectively.
To date, many studies have focused on the biosynthetic pathways, whole-genome expression profile, and pharmacological effects of the active ingredients from PN taproots [6][7][8][9][10][11][12]. To the best of our knowledge, no studies have yet investigated the cause of morphological differences in PN taproots, and the metabolic regulation pathways involved in taproot morphology differences remain unknown. Plant traits are shaped by two vital factors, external environment condition and genetic regulation. To a large extent, visible changes in plants are caused by a series of invisible physiological changes and molecular regulation based on transcriptional expression levels, which reveal that gene expression changes are closely correlated with wide variations in plant development characteristics [13]. In recent years, many researchers have sought to understand the regulation mechanism of plant root anatomy and architecture by using the model plant Arabidopsis thaliana and the diverse metabolic pathways involved in root development regulation [14]. Consistent with the growth of plant aerial parts, the formation and development of roots follow the internal laws of gene regulation [15]. In other words, changes in root morphology are inevitably accompanied by changes in the internal metabolic pathways. Therefore, research based on transcriptional regulation reveals the causes of plant morphological changes from the perspective of molecular mechanisms. Moreover, the rapid advances in next-generation sequencing (NGS) can provide abundant transcriptome data, which offers insights into the molecular regulation mechanism of plant growth and development, such as rhizome formation in Nelumbo nucifera, taproot thickening in PN, fruit morphology of pumpkin cultivars, secondary metabolite accumulation in tuberous roots of Aconitum heterophyllum, and improvement of disease resistance in strawberries [16][17][18][19][20]. Collectively, analysis of these data can reveal the internal molecular interaction mechanisms of various plant types, including morphology, yield, and disease resistance.
Both reactive oxygen species (ROS) and class III peroxidase (PRX) play pivotal roles in modulating plant root development. Either ROS balance or PRX activity affects plant root growth and elongation, and this pathway is completely independent of the signaling pathways of plant hormones, such as cytokinin and auxin; this has long been confirmed in A. thaliana [21][22][23]. Additionally, numerous studies have confirmed that ROS are important plant growth regulators and widely involved in various processes of plant root development, such as meristematic expansion and root elongation [24][25][26][27][28][29][30][31]. ROS homeostasis maintains a delicate balance under the coordination of ROS production and scavenging, maintaining an appropriate threshold boundary between redox potential and cytotoxicity [32]. Two important enzymatic families, PRX and NADPH oxidase family, contribute to ROS production [26,33]. ROS scavenging is a two-armed regulation mechanism. One arm comprises antioxidative enzymes, including catalase (CAT), L-Ascorbate peroxidase (APX), PRX, superoxide dismutase (SOD), and glutathione peroxidase (GPx) [34][35][36][37]. The other arm consists of nonenzymatic antioxidants, such as reduced glutathione, ascorbic acid, and flavonoids [38]. It has already been reported that oxidative stress in plant cells is regulated by the mitogen-activated protein kinase (MAPK) family cascade [39]. The MAPK family widely participates in biological processes in plants, such as cell proliferation and differentiation, as well as responding to and tolerating diverse stresses and environmental stimuli [39][40][41].
Therefore, to investigate the causes of phenotypic differences between LPN and OPN at the molecular regulation level, we performed an in-depth study of the transcriptome data of PN taproot and then verified the results. Three genes were very likely involved in regulating the taproot morphogenesis of PN, which provides a reliable explanation for the phenotypic variation of PN taproot.

Plant materials
PN (three years old) was collected at the vigorous growth stage from three sampling sites: Shiping (SP) county (N23°73′, E103°48′), Shilin (SL) county (N24°77′, E103°27′), and Luxi (LX) county (N24°52′, E103°76′), Yunnan Province, People's Republic of China. The specimens undertook the formal identification by Xiuming Cui that have been preserved in the greenhouse of the Faculty of Life Science and Technology, Kunming University of Science and Technology. Taproots with obviously different phenotypes (LPN and OPN) were collected separately from individual plants and washed, followed by measurement of their length and width. The samples were immediately stored in liquid nitrogen for further processing. Ten samples were collected from each sampling site, including five LPN and five OPN samples; thus, a total of 30 samples were collected from three sampling sites.

RNA isolation, library construction, and sequencing
Five samples each of LPN and OPN were collected from each sampling site, for a total of three mixed samples each of LPN and OPN (marked as LPN1, LPN2, and LPN3 for LPN and as OPN1, OPN2, and OPN3 for OPN) for RNA extraction and further experimentation. Total RNA per mixed sample was extracted from the taproot tissue using the RNeasy Plus Mini kit (Qiagen, Valencia, CA, USA), and purified using oligo (dT) magnetic beads. After completing both cDNA library synthesis and PCR amplification for each mixed sample, the PCR product was purified on an AMPure XP system (Beckman Coulter, Beverly, MA, USA). RNA integrity was assessed using an Agilent 2100 Bioanalyzer Instrument (Agilent Technologies, Santa Clara, CA, USA), and RNA degradation and contamination were monitored on 1% agarose gels.
One cDNA library was constructed from each mixed sample, and a total of six libraries were generated. Clustering of the index-coded samples was performed using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) based on the cBot Cluster Generation System, following the manufacturer's instructions. After cluster generation, 150 bp paired-end sequencing were performed using NGS based on the HiSeq 2500 system (Illumina). The clean reads were generated after removing the adapter, ploy-N, and low-quality of raw reads by using fastp software v0.19.4 [42].

Quality control and de novo assembly
The quality of clean data was evaluated using FASTQC software v0.11.9 (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). Sequenced reads with a per base average quality score below 28 were filtered, and the first 15 bases of each read were removed using Trimmomatic software v0.39 [43]. To create a better reference unigenes catalog, the LPN1 sample was used as the test object, and then three assemblers including Trinity v2.8.4 [44], SPAdes v3.14.1 [45], and SOAPdenove-trans v1.03 [46] were used to perform de novo transcriptome assembly of the LPN1 dataset. At the same time, BUSCO v5.1.0 [47] was used to assess the assembly quality of the three assemblers, as an indication of selecting an optimal assembler. Subsequently, based on the BUSCO assessment result and previous published reports on the comparative analysis of de novo transcriptome assemblers [48,49], we combined all trimmed reads into one sample using Trinity v2.8.4 [44] to complete the assembly process of the reference transcriptome of PN taproot for further analysis. Then, the assembled transcriptome was clustered, and redundancy was removed using CD-HIT software v4.8.1 [50] with a similarity parameter of 0.95. The longest transcripts were extracted as a reference transcriptome for subsequent functional annotation and DEG analysis between LPN and OPN.

Functional annotation and analysis
Trinotate pipeline v3.2.0 [44] was used to annotate both the reference transcriptome and protein-coding genes predicted by TransDecoder v5.5.0. Trinotate pipeline can integrate several databases including SwissProt, RNAmmer v1.2 (predicting ribosomal RNA), SignalP v5.0 (predicting signal peptides), TMHMM v2.0 (predicting transmembrane helices), and HMMER v3.3.1 (identifying protein domains) to populate its own database. In addition to the databases mentioned above, Trinotate can map BLAST results of the unigenes to GO and KEGG databases to obtain corresponding GO and KO numbers, respectively [51][52][53]. In brief, Trinotate merged Swis-sProt, Pfam, and other related databases into the SQLite database. BLASTX with an E-value cut-off of 1.0 × 10 −5 , BLASTP, RNAmmer, SignalP, TMHMM, HMMER, egg-NOG, KEGG, and GO homology searches were performed against the SQLite database. Finally, the search results were compiled in a report file as a table. The subcategories of GO terms of the annotated unigenes were visualized using Panther GO-slim, and the related metabolic pathways found from KEGG annotation were visualized using KEGG Mapper [52,54].

DEG analysis
The gene expression level of each sample from LPN and OPN was calculated by mapping clean reads back to the reference transcriptome using RNA-seq by expectationmaximization (RSEM) software v1.3.3 [55] and Bowtie 2 v2.4.1 alignment [56]. The expression matrix generated by RSEM was imported into DESeq2 for DEG analysis between LPN and OPN [57]. Genes with adjusted p < 0.05 and |log 2 fold change| > 1 were selected as significant DEGs. The PPI network diagram based on DEGs was constructed using the STRING database (https:// string-db. org) with default parameters, and the PPI network was further optimized using Cytoscape software v3.8.2 [58].

Measurement of H 2 O 2 content and enzyme activities
The H 2 O 2 content in OPN and LPN was determined using the H 2 O 2 content assay kit (Solarbio, Beijing, People's Republic of China) following the manufacturer's instructions. The outcome was expressed as H 2 O 2 content per gram of fresh taproot weight (μmol/g fresh weight). The enzyme activities of APX and PRX in OPN and LPN were measured using APX and PRX assay kits (Solarbio) following the manufacturer's instructions. The enzyme activity was calculated in terms of enzyme activity units per g of fresh taproot weight (U/g fresh weight). Three biological replicates were used for each experiment.

RT-qPCR validation of target genes
Total RNA was extracted from six PN taproot samples using the TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was reverse-transcribed using the RT6 cDNA Synthesis Kit v2 (Beijing Qingke Biotechnology, Beijing, People's Republic of China). Specific primers for each gene were designed for RT-qPCR amplification using Primer-BLAST (https:// www. ncbi. nlm. nih. gov/ tools/ primer-blast/). Actin was used as an internal reference gene to normalize the mRNA expression levels of target genes in each sample [59]. The specific primers for the target genes and actin are listed in Supplementary Table 1. Quantitative reactions were performed using the RT-qPCR Detection System (FQD-96A, Hangzhou Bioer Technology, Hangzhou, People's Republic of China). The reaction mixture (20 μL) contained 10 μL 2× T5 Fast qPCR Mix (SYBR Green I), 0.8 μL each of the forward and reverse primers, and 1 μL of template cDNA. Finally, qPCR amplification was conducted under the following conditions: 95 °C for 60 s, followed by 40 cycles of 95 °C for 15 s, 60 °C for 15 s, and 72 °C for 30 s. The relative expression of internal reference and target genes was calculated using the 2 −ΔΔCT method. Three biological replicates were used for the validation study.

Statistics of phenotype data
The typical phenotypes of OPN and LPN are shown in Fig. 1a. Length-width ratios of 30 PN taproot samples (OPN and LPN) collected from Shipin (SP), Shilin (SL), and Luxi (LX) county were calculated. Length-width ratios for OPN and LPN samples collected from the SP county were 1.42 ± 0.23 and 2.96 ± 0.90, those from SL county were 1.47 ± 0.26 and 2.37 ± 0.45, and those from LX county were 1.27 ± 0.20 and 1.95 ± 0.38, respectively. From the calculated results, the length-width ratios of OPN samples were approximately 1.00 to 1.50, whereas those of the LPN samples were approximately 2.00 to 3.00; thus, the length-width ratio of LPN was approximately twice that of OPN. These data were statistically analyzed using the PASW software v18.0.0. Statistical results (Fig. 1b) indicated that there was a significant difference in the length-width ratio between OPN and LPN taproots for each sampling site (p < 0.05).

Illumina sequencing results
After removing adapter, ploy-N, and low-quality reads, Illumina sequencing produced a total of 133,730,886 (39.11 Gb) and 114,761,595 (33.31 Gb) paired-end clean reads from three samples each of OPN and LPN, respectively (Table 1). After assessing read quality using FASTQC software and trimming the first 15 bases of each read using Trimmomatic software, the length of the retained reads was 135 bp, and all of them had a highquality score (> 28) for subsequent de novo assembly.  Table 2). The reference transcriptome generated from Trinity (after clustering, moving redundancy, and extracting longest transcripts) consisted of 121,955 unigenes with an average length of 914 bp, contig N50 = 1,774 bp, and GC content percentage = 38.27%; the length distribution statistics of all unigenes are shown in Fig. 2a.

De novo transcriptome assembly and function annotation
In total, 42,207 protein-coding genes were predicted from the reference transcriptome using the TransDecoder software. The final annotation results showed 38,313 Basic Local Alignment Search Tool (BLASTX) hits and 24,820 BLASTP hits. Moreover, 41,448, 25,706, 31,269, and 26,563 unigenes were mapped to RefSeq non-redundant proteins (Nr), the Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), and eggNOG databases, respectively, and 20,285 unigenes were common among the Nr, KEGG, GO, and eggNOG databases (Fig. 2b).

Mapping of KEGG pathways
The analysis of metabolic pathway assignment from the annotated unigenes was performed using the KEGG database.

Analysis of differentially expressed genes (DEGs)
The volcano plot of gene expression levels illustrates the distribution of up-and down-regulated genes in LPN as compared with in OPN (Fig. 4a). A total of 127 DEGs were identified between LPN and OPN, and the heatmap clearly presents the clustering results of these DEGs in OPN and LPN samples, showing that genes that were significantly up-regulated in LPN were significantly down-regulated in OPN (Fig. 4b) (Fig. 5b and c). Thus, up-regulated APX and PRX enzyme activities led to lower H 2 O 2 content in LPN taproots.

RT-qPCR analysis results
The Interestingly, the relative expression levels of the three genes in LPN were approximately 3-, 4-, and 3-fold higher than those in OPN, respectively. Thus, the overall relative expression levels of PnAPX3, PnPRX45, and PnMPK3 genes were significantly higher in LPN than in OPN (p < 0.01) (Fig. 5e). Both RT-qPCR and statistical analyses results of FPKM values confirmed the analysis of DEGs; that is, the relative expression levels of the PnMPK3, PnAPX3, and PnPRX45 genes were up-regulated in LPN.   molecule [61,62]. Depending on the type and concentration of ROS produced in the cells, different physiological reactions may occur. ROS can induce the expression of stress-responsive genes at low concentrations but can lead to oxidative damage to important biomacromolecules such as lipids, nucleic acids, and proteins, eventually causing cell death at high concentrations [63,64]. To date, many studies have reported important roles of ROS, including regulation of growth throughout plant root development, maintenance of apical meristem, and promotion of root elongation; the related internal molecular mechanisms have also been explored [24]. Specifically, the molecular mechanism of ROS-mediated control of the cell cycle, including cellular proliferation, elongation, and differentiation, has been widely demonstrated in A. thaliana and other plants [25][26][27][28][29][30]60]. Plant cell expansion and elongation are mostly determined by the plasticity and structure of the cell  (n = 3). The asterisk and double asterisks represent significant differences determined by the independent-sample t-test in PASW at p < 0.05 (*) and p < 0.01 (**), respectively wall (CW). The cellular expansion rate, namely cellular dynamics, is closely associated with the balance between CW stiffening and loosening, both of which are controlled by ROS metabolic processes [26]. Thus, ROS production and scavenging are directly related to the elongation and development of plant roots. ROS is generated mainly by the PRX and NADPH oxidase families [26,33], whereas its scavenging is mainly regulated by antioxidant enzymes such as CAT, APX, PRX, SOD, and some nonenzymatic antioxidants, such as reduced glutathione and ascorbic acid [34,38]. A number of early studies on A. thaliana have confirmed that accumulation of H 2 O 2 in plant root cells can enhance CW rigidity to hinder cellular elongation and proliferation by repressing the expression of cell cycle-related genes CyclinB and CyclinD, thus leading to shorter roots. In contrast, scavenging an excess of H 2 O 2 in plant root cells contributes to maintenance of intracellular redox homeostasis in favor of cellular proliferation rather than differentiation, thus leading to the development of longer roots [21,23,26,28]. In the present study, it was found that the H 2 O 2 concentration in LPN was significantly lower than that in OPN (p < 0.05) (Fig. 5a). As shown by green arrows in Fig. 6, low H 2 O 2 concentration in root cell, which is conducive to the expression of cell cycle-related genes CyclinB and CyclinD, to promote cell division and root growth. Therefore, it is possible that the emergence of the LPN phenotype is tightly correlated with lower H 2 O 2 concentration caused by the high expression of some antienzymes in root cells.

Discussion
APXs, bifunctional enzymes with CAT and broadspectrum peroxidase activity, are the core enzymes that scavenge ROS in plants [65,66]. In plant cells, APXs scavenge H 2 O 2 , which mainly participates in the first step of the ascorbate-glutathione cycle, catalyzing the reaction of L-ascorbate + H 2 O 2 + 2 H + = L-ascorbate + L-dehydroascorbate + 2 H 2 O to convert excess H 2 O 2 into H 2 O [67]. During plant development, an increasing number of APX-encoding genes, such as APX1-6, were detected in A. thaliana, and their activity and function have been further confirmed [68]. In addition, other APX-encoding genes have been found in rice, wheat, and potato tubers. A number of studies on the effects of APX-knockout genes on plant physiological functions, growth processes, and antioxidant metabolism suggest that APXs can influence plant growth and development by regulating cellular redox signaling pathways involved in plant growth [69][70][71][72][73]. The up-regulated expression of gene PnAPX3 in LPN was validated by RT-qPCR results (Fig. 5e) and APX enzyme activity analysis (Fig. 5b) in this study. Based on the known function of APX3, up-regulated PnAPX3 also confirmed the lower content of H 2 O 2 in LPN compared with in OPN. Therefore, it is reasonable to infer that the up-regulated PnAPX3 in LPN can maintain the content of H 2 O 2 in cells at a lower level compared with that in OPN, increase the rate of cell proliferation, and promote root elongation (green arrows in Fig. 6)   Fig. 6 A Putative network diagram of molecular regulatory mechanism leading to LPN formation. The blue, green, and purple arrows represent the metabolic pathways involving PnPRX45, PnAPX3 and PnMPK3, respectively PRXs are plant CW-localized proteins that have been shown to be involved in plant CW dynamics [74]. PRXs can generate ROS such as •OH and HOO• to promote cell elongation in the hydroxylic cycle and oxidize various substrates to polymerize typical CW lignin compounds to stiffen the CW in the peroxidative cycle; both cycles are capable of regulating the H 2 O 2 level [26,75]. To be specific, during the peroxidative cycle, some PRXs can oxidize molecules such as monolignols, suberin units, and ferulic acids linked to diverse polymers to form impregnable cross-links between CW polymers and proteins in CW network, leading to a tight CW, which directly reduces cellular elongation and expansion capacity [23,76]. However, during the hydroxylic cycle, other PRXs with extremely strong activity can produce •OH and HOO•, which can nonspecifically break polysaccharide covalent bonds in various types of organic molecules to enhance nonenzymatic CW loosening and consequently promote cellular elongation [26]. A previous study on A. thaliana confirmed that the expression inhibition of genes encoding PRXs results in an increase in H 2 O 2 content, whereas up-regulated expression of genes encoding PRXs can continually scavenge excess H 2 O 2 in plant root cells and maintain intracellular redox homeostasis in favor of cellular proliferation rather than differentiation, leading to the development of longer roots [21]. Briefly, the function of PRXs can be classified into two main categories: stiffening and loosening CW. Additionally, based on previous studies of multiple phenotypes of the corresponding AtPrx-deficient mutants of A. thaliana (such as mutants deficient in AtPrx02/25/71, AtPrx33/34, and AtPrx53), one PRX category, including AtPrx37 [77], AtPrx02/25/71 [78], AtPrx64 [79], and AtPrx72 [80], can promote the hardening of the plant CW by taking H 2 O 2 as an electron acceptor to catalyze lignin formation and polymerization [81][82][83]. In other words, they can catalyze lignin formation and polymerization to stiffen the plant CW in the presence of H 2 O 2 , thus leading to shorter roots. Other categories of PRXs, including AtPrx36, AtPrx39, AtPrx40, AtPrx57, AtPrx33/34, and AtPrx53, were found to be closely associated with CW loosening and root elongation [22,84,85]. Specifically, A. thaliana with overexpression of gene encoding PRX34 had longer roots than the wild type (WT), whereas the plants with double knockdown of genes encoding PRX33 and PRX34 had shorter roots than the WT. The gene identified in this study, which PnPRX45 (named AtPrx45 in A. thaliana), was significantly up-regulated in LPN compared with in OPN (Fig. 5c). PRX45, like other PRXs, mainly exists in roots and hypocotyl tissues and acts in response to oxidative stress according to description in A. thaliana database. However, the specific molecular regulation mechanism of PRX45 in the plant root CW has not been reported. The up-regulated expression of the gene PnPRX45 in LPN was validated by RT-qPCR results (Fig. 5e) and PRX enzyme activity analysis (Fig. 5c). It is reasonable to speculate that PnPRX45 may participate in the hydroxylic cycle to produce •OH and HOO• and split glycosidic bonds or certain covalent bonds in CW components, resulting in cell elongation and expansion, and eventually root elongation (blue arrow in Fig. 6). The specific molecular function of PRX45 may be comprehensively confirmed in A. thaliana roots, similar to other AtPrxencoding genes in the foreseeable future.
MAPKs, also known as MPKs in plants, are a highly conserved enzyme family with essential functions; they are widely found in plant and animal cells. The MAPK cascade is composed of diverse proteinases involved in various biological processes, including cell proliferation, differentiation, response to diverse stresses, and tolerance to environmental stimuli [39][40][41]. For instance, it has been shown that the activity of PsMPK2 in pea and ZmMPK5 in maize can be elevated by H 2 O 2 stress [86,87]. TaMPK4 plays a critical role in mediating plant tolerance to various stresses by inducing root growth and regulating cellular ROS metabolism [88]. StMAPK11 upregulation can enhance CAT and PRX activity to increase the antioxidant activity in potato, tobacco, and A. thaliana [89]. The activities of MPK3 and MPK6 in A. thaliana are positively correlated with plant defense against oxidative stress triggered by salt stress [90]. Importantly, the regulatory functions of the MAPK cascade in plant shoot apical meristem (SAM) have long been proposed, as MPK3 and MPK6 are crucial regulators of stem cell homeostasis in A. thaliana by participating in CLAVATA peptide receptor-WUSCHEL transcription factor (CLV-WUS) signaling pathways of SAM development [91,92]. In addition, two recently reported studies indicated that MPK3 can positively regulate root meristem growth factor 1 (RGF1)-mediated root growth and development and promote cell division in the root apical meristem, leading to plant root elongation [93,94]. In the present study, RT-qPCR was conducted to validate the expression of the PnMPK3 gene in LPN and OPN. The validation results were consistent with the analysis of DEGs; the PnMPK3 gene was up-regulated in LPN as compared with in OPN (Fig. 5e). Based on recent studies on MPK3 function in A. thaliana, it is rational to assume that up-regulated PnMPK3 gene is more favorable for stem cell maintenance of plant root meristem and can mediate RGF1 expression, leading to promotion of cell proliferation in the PN taproot. Therefore, up-regulated PnMPK3 activity is likely to be one of the driving forces for the formation of the LPN phenotype (purple arrows in Fig. 6).
Additionally, previous studies have confirmed that plant root development is jointly determined by the rates of cell proliferation and the extent of cell elongation [30,95]. Several published studies on different sweet potato root types suggest that up-regulated antioxidant enzyme levels could improve plant root growth and development, and under certain circumstances, may increase yield [26,33,96]. In the present study, DEGs analysis showed that the PnAPX3, PnPRX45, and PnMPK3 genes were significantly up-regulated in LPN as compared with in OPN. Further, H 2 O 2 content, APX and PRX enzyme activity, and RT-qPCR analyses in LPN and OPN further verified the transcriptome analysis. This illustrates that PnAPX3, PnPRX45, and PnMPK3 may be directly or indirectly involved in the process of promoting PN taproot development and elongation. A hypothetical molecular regulatory network leading to LPN formation may be based on the joint interference of H 2 O 2 , PnAPX3, PnPRX45, and PnMPK3 (Fig. 6).

Conclusions
In this study, we performed de novo transcriptome assembly and functional annotation from six PN taproot samples (three each of LPN and OPN) and determined the causes of phenotypic differences in the development process of PN by analyzing DEGs of LPN and OPN. DEGs analysis showed that PnAPX3, PnPRX45, and PnMPK3 genes were significantly up-regulated in LPN compared with in OPN. These three enzymes (APX, PRX and MPK) play pivotal roles in ROS metabolism and oxidative stress. Many previous studies have demonstrated that the process of ROS metabolism is closely related to the proliferation, elongation, and differentiation of plant root cells. It has recently been reported that MPK3 can positively regulate RGF1-mediated root growth and is indispensable for stem cell maintenance in the shoot apical stem of A. thaliana.
In summary, we confirmed that the PN taproot phenotype is influenced by a network controlling ROS metabolism during the taproot development process. Based on the results of this study and those of previously published studies on the relationship between ROS metabolism and plant root development, it can be concluded that the taproot phenotype of LPN is due to the up-regulated expression of PnAPX3, PnPRX45, and PnMPK3. The results of this study provide a reliable explanation for the phenotypic differences between OPN and LPN, and they offer further insights into the genetic mechanism of phenotypic differences for other species of the Panax genus. Our results will be useful for future molecular breeding of PN.