Comparative Transcriptome Analysis Reveals the Effect of Lignin on Storage Roots Formation in Two Sweetpotato (Ipomoea batatas (L.) Lam.) Cultivars

Sweet potato (Ipomoea batatas (L.) Lam.) is one of the most important crops with high storage roots yield. The formation and expansion rate of storage root (SR) plays a crucial role in the production of sweet potato. Lignin affects the SR formation; however, the molecular mechanisms of lignin in SR development have been lacking. To reveal the problem, we performed transcriptome sequencing of SR harvested at 32, 46, and 67 days after planting (DAP) to analyze two sweet potato lines, Jishu25 and Jishu29, in which SR expansion of Jishu29 was early and had a higher yield. A total of 52,137 transcripts and 21,148 unigenes were obtained after corrected with Hiseq2500 sequencing. Through the comparative analysis, 9577 unigenes were found to be differently expressed in the different stages in two cultivars. In addition, phenotypic analysis of two cultivars, combined with analysis of GO, KEGG, and WGCNA showed the regulation of lignin synthesis and related transcription factors play a crucial role in the early expansion of SR. The four key genes swbp1, swpa7, IbERF061, and IbERF109 were proved as potential candidates for regulating lignin synthesis and SR expansion in sweet potato. The data from this study provides new insights into the molecular mechanisms underlying the impact of lignin synthesis on the formation and expansion of SR in sweet potatoes and proposes several candidate genes that may affect sweet potato yield.


Introduction
Sweet potato is one of the seventh most important food crops in the world [1] and produces approximately 88.9 million tons of storage root (SR) from an area of 7.4 million ha [2]. China is the largest sweet potato-producing country with 47.8 million tons, which accounts for 53.8% of the total production of sweet potato worldwide. SR is the main edible tissue, its formation and development is the most important agronomic trait in sweet potato producing, which is accompanied by complex biological processes such as adventitious root morphogenesis and accumulation of carbohydrates, storage proteins, and secondary metabolites [3]. SR yield varied in sweet potato cultivars and is prone to environmental change [4][5][6]. It was reported that SR yield is not only dependent upon the rate and duration of SR expansion but also on the beginning of SR formation, the leaves' longevity, and the growth stage [7,8].
Studies have shown that lignin biosynthesis is connected with root formation and development in sweet potatoes. Firon compared the expression profiles of initiating SRs and fibrous roots and identified the lignin biosynthesis down-regulated at an early stage of SR formation [9]. Ectopic expression of maize Lc regulatory gene in sweet potato induced the expression of lignin biosynthesis genes and affected SR development [11]. Kim discovered many differently expressed genes related to phenylpropanoid biosynthesis in adventitious root formation through RNA-seq analysis [20]. However, the molecular mechanism of lignin in SR development keeps it obscure.
In this study, we selected SR at three developmental stages of two sweet potato varieties, compared the full-length transcriptome data and investigated the gene expression profiles by using full-length and second-generation transcriptome. The primary objective is to reveal the molecular mechanism by which lignin synthesis affects the formation and expansion of sweet potato SR.

Plant Materials
Sweet potato cultivars, "Jishu25 (J25)" and "Jishu29 (J29)", were planted in the experimental station of the Crops Research Institute, Shandong Academy of Agricultural Sciences, Jinan, China. SRs at the three stages were collected from sweet potato plants 32 days after planting (DAP) (D1), 46 DAP (D2), and 67 DAP (D3) [21], respectively. Three independent biological replicates were taken from each stage of each variety, and each biological replicate came from three independent sweet potato SRs. The samples were immediately frozen in liquid nitrogen and stored at −80 • C for hormone and RNA isolation. At the root developmental stages, the numbers and yield of SR were counted and weighed. The root/shoot (R/S) ratio was calculated according to the root weight divided by the shoot biomass of per plant.

Determination of Lignin Content
Analysis of the lignin content was performed with slight modification according to the methods described by [22]. In short, 10 mg dried roots (W) were digested in 2.5 mL HoAc solution containing 25% (v/v) acetylbromide and 0.1 mL 70% perchloric acid. The sealed vessel was mixed fully and incubated for 40 min at 80 • C while shaking. After incubation and cooling, the slurry was centrifuged at 23,477× g for 15 min. The supernatant was added to 2.5 mL of 2 M NaOH and 1 mL acetic acid. After 20 min, the absorbance (A) was measured at 280 nm. The lignin content was calculated using the following formula: lignin (mg/g) = 0.147 × (∆A − 0.0068) ÷ W × 50.

RNA Extraction, Full-Length cDNA Library Construction and Sequencing
Total RNAs were extracted using RNA prep Pure Plant plus Kit (Tiangen Biotech (Beijing) Co., Ltd., Beijing, China) and purified with the RNA easy Plant Mini Kit (Qiagen, Valencia, CA, USA). RNA quality was verified using a 2100 Bioanalyzer RNA Nanochip For full-length cDNA library, cDNA was synthesized using SMARTer PCR cDNA Synthesis Kit, and optimized for PCR amplification. The fragments for large-scale PCR were performed using magnetic beads to obtain sufficient total cDNA. The complete SMRT bell library was constructed with using a SMARTer PCR cDNA Synthesis Kit and assembly was performed on the PacBio Sequel platform, the second-generation sequencing and assembly was implemented on the Hiseq 2500 sequencing platform (Illumina) with PE150 by Novogene Co., Ltd. (Beijing, China).

Real-Time Quantitative PCR Validation
To confirm the expression of unigenes, 12 unigenes were selected for qRT-PCR. The analysis was performed using samples with tuberous roots at the three stages from two cultivars. The qRT-PCR and data analyses were performed as described by [28]. A total of 12 unigenes, including 4 plant hormone biosynthesis-related genes, 5 lignin biosynthesisrelated genes, and 3 transcription factors, from the RNA-Seq were validated, and the primers used for the validation were listed in Supplementary Table S4. Sweet potato Ibactin gene was used as the reference gene for normalizing quantities of gene expression.

Statistical Analysis
Statistical analysis was performed using SPSS version 13.0 with ANOVA and Duncan's test for dry matter, lignin, hormone, and qPCR results. Data are means with three biological replicates, with each error bar representing standard error. The statistical significance difference was calculated with Duncan's new multiple ranges test and marked with asterisks at p < 0.05.

Characteristic of SR Development at Different Stages in Two Cultivars
We analyzed the phenotypes of two varieties (cv. Jishu25 (J25) and cv. Jishu29 (J29)) at three different time periods (D1, D2, and D3) ( Figure 1A), the two varieties had a significant difference in the characteristics of SR. The expanding level, number of expansion SRs, and yield of J29 were significantly higher than J25 at D1, D2, and D3 ( Figure 1B Statistical analysis was performed using SPSS version 13.0 with ANOVA and Duncan's test for dry matter, lignin, hormone, and qPCR results. Data are means with three biological replicates, with each error bar representing standard error. The statistical significance difference was calculated with Duncan's new multiple ranges test and marked with asterisks at p < 0.05.

Characteristic of SR Development at Different Stages in Two Cultivars
We analyzed the phenotypes of two varieties (cv. Jishu25 (J25) and cv. Jishu29 (J29)) at three different time periods (D1, D2, and D3) ( Figure 1A), the two varieties had a significant difference in the characteristics of SR. The expanding level, number of expansion SRs, and yield of J29 were significantly higher than J25 at D1, D2, and D3 ( Figure 1B-D).
In the D2 stage, the SR of J25 had just begun to expand, but there were already obviously expanded SRs in J29, and the SR number of J29 during the D2 stage was significantly higher than J29. The yield and R/S value also indicate that J29 starts root expanded earlier than J25, and the SR expansion of J29 storage was faster during the D1-D2 stage. Data are means ± SE of three biological repeats, error bars indicate error standard. Means denoted by the same letter were not significantly different at P > 0.05, and different letters indicate statistically significantly differences (p < 0.05).

Determination of Lignin Content
Lignin affects SR development in sweet potato. To understate the dynamic changes in lignin accumulation in sweet potato, the lignin content in SR was detected during the root development stages in two cultivars. The results showed that lignin content in SR decreased gradually from D1 to D3 stages in the two cultivars, and the lignin content in J29 was significantly higher than those of J25 at the D1 and D2 stages ( Figure 1E). In the D2 stage, the SR of J25 had just begun to expand, but there were already obviously expanded SRs in J29, and the SR number of J29 during the D2 stage was significantly higher than J29. The yield and R/S value also indicate that J29 starts root expanded earlier than J25, and the SR expansion of J29 storage was faster during the D1-D2 stage.

Determination of Lignin Content
Lignin affects SR development in sweet potato. To understate the dynamic changes in lignin accumulation in sweet potato, the lignin content in SR was detected during the root development stages in two cultivars. The results showed that lignin content in SR decreased gradually from D1 to D3 stages in the two cultivars, and the lignin content in J29 was significantly higher than those of J25 at the D1 and D2 stages ( Figure 1E).

Global Analysis for RNA-Seq Data
To compare the molecular mechanisms of tuberous root development of both cultivars, the samples were sampled in three stages (D1, D2, and D3). Three independent biological replicates were taken from each stage of each variety, and each biological replicate came from three independent sweet potato SRs (Table S1). Meanwhile, the RNA of each sample was mixed for SMRT sequencing. Transcriptome sequencing analysis yielded a total of 142.931 GB of clean data, with 7.94 GB of data per sample on average. In addition, Q20 ≥ 96.54% and Q30 ≥ 90.95% were identified in all samples (Supplementary Table S1). We obtained 741,084 circular consensus sequencing (CCS) reads ranging from 51 bp to 14,896 bp with an average length of 1322 bp, which included 77.1% of full-length non-chimeric (FLNC) and 20.7% of non-full-length (NFL) reads. After correcting with Hiseq2500 sequencing, 52,137 transcripts and 21,148 unigenes were obtained with an average length of 1253 bp and 1457 bp, which indicates that the data are of high quality (Supplementary Table S2). The sequence length after redundancy varied from 53 bp to 6868 bp with the mean length of 1457 bp. Pearson analysis of the transcriptome data found that the three replicates of each line had good consistency and met the requirements of subsequent analysis (excepted for J29_D2_3, which has been deleted) (Supplementary Figure S1).

Annotation and Classification of Sweetpotato Unigenes
To identify the predictive functions of unigenes in sweet potato SRs, all of the assembled unigenes were matched against the seven databases. Based on the sequence similarity, a total of 96.51% of the unigenes were annotated by alignment in at least one database. Of them, 19 Table S3).

Analysis of Differential Expression Genes (DEGs)
In this study, we analyzed the global gene expression profiles of sweet potato SR in different developmental stages. According to the false discovery rate and fold change, 14,598 DEGs were screened through cluster analysis in the SRs during the three different developmental stages of two sweet potato cultivars. To examine the gene expression difference during the root development stages in the two genotypes, the DEGs were identified by the comparisons of the nine DEG libraries, i.e., J25-D2 vs. J25-D1 and J25-D3 vs. J25-D2 (Supplementary Figure S5). The largest number of DEGs occurred between J29-D3 vs. J25-D3 with 3063 up-regulated and 4304 down-regulated unigenes. Furthermore, 6675 and 6571 unigenes were significantly differentially expressed between J25-D3 vs. J25-D2, and J29-D1 vs. J25-D1, respectively. It is interesting that the number of DEGs marked in the two genotypes increased from D1-D2 comparison to D2-D3 ( Figure 2). Meanwhile, the number of up-or down-regulated DEGs in Jishu25 was much bigger than those in J29, which revealed that the transcriptome of J25 changed drastically compared to J29 ( Figure 2). Moreover, cluster analysis of 14,598 DEGs also showed this result (Figure 3). and 6571 unigenes were significantly differentially expressed between J25-D3 vs. J25-D2, and J29-D1 vs. J25-D1, respectively. It is interesting that the number of DEGs marked in the two genotypes increased from D1-D2 comparison to D2-D3 ( Figure 2). Meanwhile, the number of up-or down-regulated DEGs in Jishu25 was much bigger than those in J29, which revealed that the transcriptome of J25 changed drastically compared to J29 ( Figure  2). Moreover, cluster analysis of 14,598 DEGs also showed this result (Figure 3).

Pathway Analysis of DEGs
To determine the involvement of these differentially expressed genes in SRs, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment of DEGs was performed in J29-D2 vs. J29-D1 and J25-D2 vs. J25-D1. The upregulated genes in J29-D2 vs. J29-D1 were identified to be involved in 29 distinct metabolic pathways. Of them, the top and J29-D1 vs. J25-D1, respectively. It is interesting that the number of DEGs marked in the two genotypes increased from D1-D2 comparison to D2-D3 ( Figure 2). Meanwhile, the number of up-or down-regulated DEGs in Jishu25 was much bigger than those in J29, which revealed that the transcriptome of J25 changed drastically compared to J29 ( Figure  2). Moreover, cluster analysis of 14,598 DEGs also showed this result (Figure 3).

Pathway Analysis of DEGs
To determine the involvement of these differentially expressed genes in SRs, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment of DEGs was performed in J29-D2 vs. J29-D1 and J25-D2 vs. J25-D1. The upregulated genes in J29-D2 vs. J29-D1 were identified to be involved in 29 distinct metabolic pathways. Of them, the top

DEGs in Starch and Sucrose Metabolism
Starch is the main component of dry matter in sweet potato SR, while sucrose is the main form of long-distance transportation of assimilated carbon in sweet potato

DEGs in Phenylpropanoid Biosynthesis
Lignin is an important secondary metabolite in plants and plays an important biological role in plant growth and development. In two cultivars D1 and D2, a total of 41 DEGs were involved in phenylpropanoid biosynthesis.

Differentially Expressed Transcription Factors in Root Development
The development of SRs is regulated by transcription factors (TFs) that control various key gene expressions. A total of 1063 transcription factors were identified in three developmental stages of two varieties. Among them, the largest types were AP2/ERF-ERF (103), bHLH (80), C3H (66), MYB (64), bZIP (61), NAC (56), and WRKY (52) (Figure 6), which were reported to be related to root formation and development [30][31][32]. NAC and WRKY TF are the key regulators of the lignifications of vessel cell differentiation [33,34]. In addition, 47 GRAS and 47 AUX/IAA transcription factors which play a crucial role in gibberellins and auxin signal transduction [35,36]. The development of SRs is regulated by transcription factors (TFs) that control various key gene expressions. A total of 1063 transcription factors were identified in three developmental stages of two varieties. Among them, the largest types were AP2/ERF-ERF (103), bHLH (80), C3H (66), MYB (64), bZIP (61), NAC (56), and WRKY (52) (Figure 6), which were reported to be related to root formation and development [30][31][32]. NAC and WRKY TF are the key regulators of the lignifications of vessel cell differentiation [33,34]. In addition, 47 GRAS and 47 AUX/IAA transcription factors which play a crucial role in gibberellins and auxin signal transduction [35,36]. Figure 6. Transcription factors family analysis. The horizontal axis represents the names of different transcription factor families, while the vertical axis represents the number of detected transcription factor families. To identify the DEGs correlating with tuberous root formation and development, WGCNA was implemented to construct a gene network from 27 sweet potato root samples in two cultivars by R package [37]. In the analysis, 17 stable co-expressed modules were obtained through WGCNA ( Figure 7A). In the brown module, GO enrichment analysis showed that the unigenes were mainly involved in transferase activity and transcription factor activity in molecular function category, S-adenosylmethionine biosynthetic and metabolic process in biological process category, and apoplast in cellular component ( Figure 7B and Supplementary Table S6). KEGG enrichment analysis showed that the unigenes were significantly enriched in phenylpropanoid biosynthesis, and plantpathogen interaction ( Figure 7C). To identify the DEGs correlating with tuberous root formation and development, WGCNA was implemented to construct a gene network from 27 sweet potato root samples in two cultivars by R package [37]. In the analysis, 17 stable co-expressed modules were obtained through WGCNA ( Figure 7A). In the brown module, GO enrichment analysis showed that the unigenes were mainly involved in transferase activity and transcription factor activity in molecular function category, S-adenosylmethionine biosynthetic and metabolic process in biological process category, and apoplast in cellular component ( Figure 7B and Supplementary Table S6). KEGG enrichment analysis showed that the unigenes were significantly enriched in phenylpropanoid biosynthesis, and plant-pathogen interaction ( Figure 7C).

Validation of RNA-Seq Results
To validate our transcriptome data, qRT-PCR analyses were performed to determine the expression of 12 random DEGs in three root developmental stages of two cultivars. The qRT-PCR results showed that the expression patterns of the 12 DEGs were in good agreement with their RNA-seq results in the root development stages in every cultivar (Figure 8 and Supplementary Table S7), and the positive correlation coefficient (R2) was 0.9195. Therefore, the transcriptome data was highly reliable.

Validation of RNA-Seq Results
To validate our transcriptome data, qRT-PCR analyses were performed to determine the expression of 12 random DEGs in three root developmental stages of two cultivars. The qRT-PCR results showed that the expression patterns of the 12 DEGs were in good agreement with their RNA-seq results in the root development stages in every cultivar (Figure 8 and Supplementary Table S7), and the positive correlation coefficient (R2) was 0.9195. Therefore, the transcriptome data was highly reliable.

Discussion
SR formation and development are important for sweet potato production. SR expansion is affected by many factors. In this study, we compared two sweet potato cultivar, J25 and J29, and found that J29 started to expand earlier in the early SR period (32-45 days), and its yield was higher (67 days). In order to clarify the differences in SR extension mechanisms between J25 and J29, we carried out comparative transcriptome analysis of two varieties in three stages (D1; 32d, D2; 45d, and D3; 67d), and analyzed the changes of gene expression in three stages. A total of 21,148 genes were identified and annotated into NR, KOG, Pfam, Swiss Prot, and GO databases.
Previous studies have extensively understood the anatomical structure of sweet potato roots. In fibrous roots, the degree of lignification of columnar cells is high, and the activity of vascular cambium is weak, on the contrary, in thick roots, the degree of lignification of columnar cells is high, and the activity of vascular cambium is strong [38].The

Discussion
SR formation and development are important for sweet potato production. SR expansion is affected by many factors. In this study, we compared two sweet potato cultivar, J25 and J29, and found that J29 started to expand earlier in the early SR period (32-45 days), and its yield was higher (67 days). In order to clarify the differences in SR extension mechanisms between J25 and J29, we carried out comparative transcriptome analysis of two varieties in three stages (D1; 32d, D2; 45d, and D3; 67d), and analyzed the changes of gene expression in three stages. A total of 21,148 genes were identified and annotated into NR, KOG, Pfam, Swiss Prot, and GO databases.
Previous studies have extensively understood the anatomical structure of sweet potato roots. In fibrous roots, the degree of lignification of columnar cells is high, and the activity of vascular cambium is weak, on the contrary, in thick roots, the degree of lignification of columnar cells is high, and the activity of vascular cambium is strong [38].The carbon flux distribution in the starch and lignin metabolic pathways can affect the development of fibrous roots towards pencil roots or SRs [11]. Down-regulation of lignin biosynthesis and up-regulation of starch biosynthesis are the main events leading to the initiation of SRs [9].
In this study, we found that the overall number and yield of SR in J29 were higher than those in J25. In addition, at D2, the root-shoot ratio of J29 was significantly higher than that of J25. The measurement of lignin content showed that After D1 to D3 of transplantation, the lignin content of both varieties decreased, but J29 had a more significant reduction in lignin at D1 and D2 (the lignin content in J29 was significantly higher than that in J25 at D1, and it was slightly lower than that in J25 at D3). Research shows GA promotes lignification and secondary wall formation; SR formation is accompanied by marked reductions in GA signaling [39]. High levels of lignification appear to be detrimental to storage organ formation [40], the application of exogenous GA on sweet potato branches leads to a down-regulation of starch biosynthesis genes, while an up-regulation of lignin biosynthesis genes enhances root lignification, leading to a decrease in SR expanding [41]. Overexpression of the maize LC gene in sweet potatoes stimulates lignin biosynthesis, leading to enhanced lignification of vascular cells in early SRs, severely reducing the expansion of SRs [11]. Therefore, we selected swbp1 and swpa7 (peroxidase; transcript24226/f8p0/1239 and transcript26346/f3p0/1164), two DEGs involved in phenylpropanoid biosynthesis, peroxidase is widely involved in plant physiological processes including growth, maturation, seed germination, and regulation and crosstalk of plant hormone signals. It acts downstream of the phenylpropanoid pathway and aggregates lignin monomers to form lignin A [42]; Increased lignin and phenolic content in transgenic plants overexpressing the sweet potato peroxidase gene swpa4 [43]; Lee found that overexpression of the sweet potato peroxidase gene IbLfp increased lignin content in SRs [44]. The quantitative RT-PCR results indicated that swbp1 and swpa7 were strongly down-regulated during SR expansion, especially the expression level in J29 was significantly lower than that in J25 ( Figure 9A,B and Supplementary Table S8). These results indicate that swbp1 and swpa7 might play a critical role in SR expansion by reduction of lignin content.
The development of SR and lignin synthesis in plants is regulated by multiple transcription factors. Transcription factors that regulate lignin synthesis have been identified in various plants instantaneously. The AP2/ERF transcription factors (TFs) regulate various processes of plant growth, development, and response to environmental stimuli [45]. Expression of CsERF1B in citrus peel can enhance the activity of enzymes related to lignin synthesis, such as POD and COMT, and promote lignin accumulation [46], overexpression of ERF139 significantly increases the total lignin accumulation in hybrid poplar [47], overexpression of sweet potato ERF transcription factor IbRAP2.4 inhibits SR formation by activating the expression of genes involved in lignin biosynthesis pathways [48]. We identified two AP2/ERF transcription factors, IbERF061 and IbERF109 (transcript22489/f2p0/1329 and transcript28975/f5p0/1069) which down-regulated during SR development, and our qPCR results revealed that the expression of two ERF transcription factors in J29 is significantly higher than those in J25 ( Figure 9C,D and Supplementary Table S8). This is consistent with the trend of lignin synthesis regulatory genes described above. We speculate that during the SR development of J29, the POD activity in the SRs may be affected by reducing the expression level of IbERF061 and IbERF109 transcription factors, resulting in a decrease in lignin content and earlier expansion than J25. that overexpression of the sweet potato peroxidase gene IbLfp increased lignin content in SRs [44]. The quantitative RT-PCR results indicated that swbp1 and swpa7 were strongly down-regulated during SR expansion, especially the expression level in J29 was significantly lower than that in J25 ( Figure 9A,B and Supplementary Table S8). These results indicate that swbp1 and swpa7 might play a critical role in SR expansion by reduction of lignin content.

Conclusions
The rate of SR formation and expansion affects sweet potato yield. In this study, we compared the transcriptomes of SR at the three different development stages in the two cultivars. The results of de novo assembly identified and annotated 52,137 transcripts and 21,148 unigenes were obtained after corrected with Hiseq2500 sequencing. Through the comparative analysis, 9577 unigenes were found to be differently expressed in the different stages of the two cultivars. Phenotypic analysis of two cultivars, combined with analysis of GO, KEGG, and WGCNA showed the regulation of lignin synthesis and related transcription factors plays a crucial role in the early expansion of SR. The four key genes swbp1, swpa7, IbERF061, and IbERF109 were proved as potential candidates for regulating lignin synthesis and SR expansion in sweet potato. However specific experiments are needed to further verify the function of these genes.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes14061263/s1, Figure S1: Pearson analysis of the transcriptome data in the independent biological replicates; Figure S2: Gene ontology (GO) classifications in sweet potato (Jishu25 and Jishu29 cultivars); Figure S3: Gene ontology (GO) classifications in sweetpotato (Jishu25 and Jishu29 cultivars); Figure S4: Clusters of orthologous groups (COG) classification in Jishu25 and Jishu29 sweetpotato cultivar. Genes from the same Orthologous have the same function, so that direct functional annotations to other members of the same KOG cluster; Figure S5: Volcano plot of differentially expressed unigenes in different developmental stages in Jishu25 and Jishu29; Table S1: The data quality of raw sequencing; Table S2: Statistics of transcripts length distribution before and after correction; Table S3: Unigenes annotation information; Table S4: The primer of qRT-PCR; Table S5: Expression analysis of key DEGs involved in different pathways; Table S6: Enriched GO terms genes name; Table S7: qRT-PCR genes name; Table S8: qRT-PCR analysis of four key genes.

Conflicts of Interest:
The authors declare no conflict of interest.