Identification, expression analysis, and potential roles of microRNAs in the regulation of polysaccharide biosynthesis in Polygonatum cyrtonema Hua

MicroRNAs (miRNAs) are a group of endogenous small non-coding RNAs with important roles in plant growth, development, and metabolic processes. Polygonatum cyrtonema Hua (P. cyrtonema) is an important Chinese traditional medicinal herb with broad pharmacological functions, and polysaccharides are the main biological substance accumulated in the P. cyrtonema rhizome. However, regulation of the process of polysaccharide biosynthesis in P. cyrtonema remains largely unknown.To elucidate the miRNAs and their targets involved in polysaccharide biosynthesis in P. cyrtonema, four small RNA libraries were constructed from flower, leaf, rhizome, and root tissues and sequenced. A total of 69 conserved and 5 novel miRNAs were identified, of which 6 miRNAs (miR156a-5p, miR156f-5p, miR395a-5, miR396a-3p, miR396g-3p, and miR397-5p-1) were down-regulated and 7 miRNAs (miR160, miR160h-1, miR160e-5p, miR319b-1, miR319-1, miR319c-5p-3, and miR319c-1) were up-regulated in rhizomes compared with flower, leaf, and root tissues. Bioinformatics analysis showed that the predicted targets of these miRNAs were mostly transcription factors and functional genes enriched in metabolic and secondary metabolite biosynthetic pathways, and 7 genes and their paired miRNAs were identified in carbohydrate metabolism. qRT-PCR expression analysis demonstrated that 6 miRNAs and their targets involved in carbohydrate metabolism were existed a negative correlation in. P. cyrtonema tissues. MiR396a-3p and one of its target genes, abfA, were possibly involved in polysaccharide biosynthesis pathway. This is the first report on the identification of conserved and novel miRNAs and their potential targets in P. cyrtonema, thus providing molecular evidence for the role of miRNAs in the regulation of polysaccharide biosynthesis in P. cyrtonema.


Introduction
Polygonatum cyrtonema Hua (Liliaceae) is a perennial plant that grows mainly in the temperate regions of the northern hemisphere, including China, Japan, Korea, and Russia (Zhao et al. 2018b). Due to its antioxidant, antiaging, and neuroprotective effects, P. cyrtonema is widely used as a traditional Chinese medicine for the treatment of cardiovascular and other diseases in China (Cui et al. 2018;Zhao and Li 2015). The rich biological activity of P. cyrtonema is attributed to the presence of various compounds such as alkaloids, flavonoids, steroidal saponins, lignin, and polysaccharides (Yelithao et al. 2016;Zhao and Li 2015), and polysaccharides are considered as one of the most important groups of active compounds in P. cyrtonema Zhu et al. 2015). Recently, P. cyrtonema polysaccharides have gained more attention for their various pharmacological applications and biological activities. However, most of the research is focused on the extraction, chemical composition, and molecular structure of polysaccharides, and few studies are conducted on the genetic and epigenetic regulation of polysaccharides in P. cyrtonema. Consequently, molecular mechanisms of polysaccharide biosynthesis and regulation remain unknown.
MicroRNAs (miRNAs) are a group of small non-coding RNAs that play crucial roles in various biological and metabolic processes underlying plant growth and development (Hou et al. 2019). Generally, miRNAs exhibit tissue or stage-specific expression patterns, and participate in the regulation of the accumulation of secondary metabolites such as steroids and ketones by negatively regulating the expression of target genes (Lee and Carroll 2018;Megha et al. 2018;Wei et al. 2015). For instance, miRNAs are involved in the biosynthesis of terpenoid indole alkaloids in Catharanthus roseus (Shen et al. 2017). Additionally, miR156, miR828, miR858, and miR5072 regulate anthocyanin accumulation in apple peel (Qu et al. 2016). In Brazilian rubber tree, miR159b regulates the production of latex (Gebelin et al. 2013). Therefore, the identification of miRNAs and their target genes is essential for understanding the regulation of miRNA-mediated biosynthesis of polysaccharides in P. cyrtonema.
Because miRNAs are important regulators of gene expression, they have been thoroughly studied in many species. However, despite the medicinal significance of P. cyrtonema, miRNAs have not yet been identified in this species. Therefore, the aim of this study was to identify conserved and novel miRNAs and their potential target genes in P. cyrtonema. To achieve this goal, we sequenced small RNAs (sRNAs) from different tissues of P. cyrtonema using the BGISEQ-500 platform, analyzed miRNA expression profiles, and investigated the functions of target genes. Quantitative real-time PCR (qRT-PCR) was used to verify the results of sRNA sequencing and to identify candidate genes involved in the accumulation of polysaccharides in P. cyrtonema. Our results provide valuable information on the miRNAs involved in polysaccharide biosynthesis in P. cyrtonema.

De novo sRNA-seq analysis
To identify miRNAs in P. cyrtonema, four cDNA libraries prepared from flower, leaf, rhizome, and root tissues were sequenced using the BGISEQ-500 sequencing technology. Approximately 162 million raw reads were generated from the cDNA libraries. The raw reads were cleaned to obtain 36,516,720 reads from flower library (PcH_F), 41,537,414 reads from leaf library (PcH_L), 35,905,873 reads from rhizome library (PcH_RH), and 32,754,124 reads from root library (PcH_R) (Fig. 1A, Table S1). The quality score (Q20) of all four libraries exceeded 99%, indicating good quality sequence data. The distribution of base quality in clean reads is shown in Figure S1. Analysis of the length and composition of sRNAs revealed that 24 nt long sRNAs were the most abundant in the four libraries, followed by 21 nt and 22 nt sRNAs (Fig. 1B). In the PcH_F, PcH_L, PcH_RH, and PcH_R libraries, 24 nt sRNAs accounted for 39.42%, 29.28%, 36.79%, and 18.47% of all sRNAs, respectively; similarly, 21 nt sRNAs accounted for 15.81%, 19.88%, 16.03%, and 12.29% of all sRNAs, respectively, while 22 nt sRNAs accounted for 9.21%, 11.77%, 13.42%, and 12.74% of all sRNAs, respectively. The remaining sRNAs accounted for less than 10% of each sample.
A total of 20 miRNA families comprising 69 conserved miRNAs were identified in P. cyrtonema; 65 in PcH_F, 59 in PcH_L, 62 in PcH_RH, and 57 in PcH_R ( Fig. 2A, Table S4). The length of conserved miRNAs varied from 18 to 22 nt, the majority of which were 21 nt long (37 out of 69 conserved miRNAs; 53.62%), followed by 20 nt miRNAs (23 out of 69 conserved miRNAs; 33.33%) ( Table S3). Precursors of conserved miRNAs were also identified, with lengths ranging from 19 to 269 nt. Statistical analysis of the first base of predicted miRNAs showed a strong preference for uridine (U) in 21 nt long miRNAs (Fig. S3A, B). The largest family identified was miR319 with 10 members, followed by miR396 and miR166, with 8 and 7 members, respectively. Of the remaining 17 families, 12 families were represented by 2-13 members while 5 were represented by a single member (Fig. 2B).
To predict potentially novel miRNAs in P. cyrtonema, miRDeep2 and miRA software were used to explore the stem-loop hairpin secondary structures and enzyme cleavage sites, and the minimum free energy of novel miRNAs was measured. Five new hairpin miRNA candidates, ranging from 21 to 30 nt in length, were identified. These miRNA candidates were potentially new miRNAs or new members of conserved miRNA families in P. cyrtonema. (Table S5). The precursors of novel miRNAs ranged from 76 to 248 nt in length (average = 156 nt), which is consistent with the typical length distribution of mature miR-NAs. Nucleotide bias analysis showed that new miRNAs had similar tendency of conserved miRNA in all libraries (Fig. S4). The abundance of novel miRNAs in different tissues varied significantly. For instance, the number of novel_miR4 in flower and leaf tissues was much higher than that in rhizome and root tissues. By contrast, the abundance of novel_miR3 and novel_miR5 was low; these were detected only in flower and leaf tissues, respectively. Similarly, novel_miR2 was detected only in leaves, whereas novel_miR1 was detected in all four tissues, except roots (Table S6).

Differential expression analysis of miRNAs
Based on the number of reads that mapped to the P. cyrtonema genome, the expression level of sRNAs was calculated as Transcripts Per kilobase Million (TPM). Correlation analysis showed that the overall expression pattern of sRNAs in rhizomes was similar to that in roots (Pearson's correlation coefficient = 1.0) but different from that in flowers and leaves (Fig. S5).
Since polysaccharides mainly accumulate in rhizomes, we further investigated the miRNAs expressed in rhizomes. The number of up-regulated and down-regulated miRNAs in rhizomes were 20 and 49 compared with flowers, 34 and 37 compared with leaves, and 45 and 24 compared with roots, respectively (Fig. 3A, Table S7). The total miRNA expression level in rhizome was the same as that in leaf, lower in flower, and higher than root (Fig. 3B). To further investigate miRNAs showing differential expression between the rhizome vs. three other tissues, we performed hierarchical cluster analysis of miRNAs using the heatmap function in R software (Fig. 3C). Specifically, the expression level of miR156a-5p, miR156f-5p, miR395a_5, miR396a-3p, miR396g-3p, and miR397-5p_1 in rhizome was significantly lower than that in flower, leaf, and root, whereas the expression level of miR160, miR160h_1, miR160e-5p, miR319b_1, miR319_1, miR319c-5p_3, and miR319c_1 in rhizome was significantly higher than that in flower, leaf, and root (Table S7).

Prediction of miRNA target genes
To better understand the potential functions of the miRNAs in P. cyrtonema, it is necessary to predict and confirm their target genes. A total of 1,525 target genes were predicted for 81 miRNAs, of which 2218 genes were predicted using psRobot, 1523 using TargetFinder, and 934 were predicted repeatedly. After data filtering, 69 miRNAs had a total of 499 target genes (Table S8), and most of the predicted target genes were involved in a variety of biological processes, including transcriptional regulation, signal transduction, and stress response (Table S9).
It is known that one miRNA may have several target genes. In our study, the number of target genes varied significantly among miRNAs. For example, miR172a_4 was predicted to have 45 target genes, while miR160a-5p and novel_mir4 had only one target. In addition, some miRNAs from the same family or different families were predicted to regulate the same target gene (Table S9). For example, the squamosa promoter-binding protein-like (SPL) transcription factor encoding gene was predicted as the target of miR156a-5p, miR156b_2, miR156f-5p, miR529e_1, and miR529h, whereas genes encoding MYB transcription factors were predicted as targets of miR159a_1, miR319_1, miR319a-3p, miR319b_1, miR319c_1, miR319c_2, miR319f_1, miR319g, miR319_1, and miR858b. Auxin response factors (ARFs) play an important role in plant growth, development, and adaptive responses, and ARF genes were predicted as targets of miR160a-5p, miR160b_1, miR160e-5p, miR160h_1, and miR167d_1 (Table S9).
Although the expression level of novel miRNAs was lower than that of conserved miRNAs, novel miRNAs were found to have candidate target genes, including transcription factor encoding genes and functional genes. For example, gene encoding mannose-specific lectin precursor was predicted as the target of novel_miR1; genes encoding low-temperature-induced 65 kDa protein and ABC transporter C family member as targets of novel_miR2, and those encoding the hypothetical protein, peroxisomal adenine nucleotide carrier 1, and probable protein phosphatase 2C 38 as targets of novel_miR4 (Table S9).

Annotation of miRNA target genes
To gain insights into the possible roles of DEMs, their putative target genes annotated using gene ontology (GO) enrichment analysis. GO analysis grouped the target genes into three main categories: biological process, cellular component, and molecular function. In the biological process category, target genes were associated with cellular process, single organism process, metabolic process, and biological regulation; the number of genes associated with each of these terms was 48, 39, 38, and 18 in the rhizome vs. flower group (Fig. 4A), 37, 35, 31, and 17 in the rhizome vs. leaf group (Fig. 4B), and 31, 16, 23, and 18 in the rhizome vs. root group (Fig. 4C). Thus, the number of differentially expressed genes identified in the rhizome vs. root group was smaller than that identified in the rhizome vs. leaf and rhizome vs. flower groups.
In the cellular component category, target genes were mainly enriched in cell, cell part, and membrane. There was no significant difference in the number of genes enriched in different tissues. However, in the molecular function category, genes were enriched in binding, catalytic activity, and nucleic acid transcription factor activity; the number of genes associated with catalytic activity were 44 in the rhizome vs. flower group, 30 in the rhizome vs. leaf group, and 26 in the rhizome vs. root group.
Genes usually act in networks to regulate biological functions, and pathway-based analysis such as Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis helps to understand the biological functions of these target genes. The results of KEGG enrichment analysis showed that target genes of DEMs were involved in 42 pathways in the rhizome versus flower group (Fig. 5A, Table S10), 39 pathways in the rhizome versus leaf group (Fig. 5B, Table S11), and 38 pathways in the rhizome vs. root group (Fig. 5C, Table S12). In addition to plant hormone signal transduction, most of these pathways were involved in biomaterial metabolism, such as metabolism of sugars, amino acids, lipids, and nucleotides and biosynthesis of secondary metabolites.

Identification of miRNAs and their putative target genes related to polysaccharide biosynthesis
In plants, miRNAs have been associated with secondary metabolism, polysaccharide metabolism, and other biological processes. Here, we screened and analyzed target genes of DEMs and identified seven genes and their paired miRNAs involved in carbohydrate metabolism including those encoding alpha-L-arabinofuranosidase (abfA, miR396a-3p target), L-galactose dehydrogenase (GalDH, miR8175 target), berberine bridge enzyme-like 18 (miR396h target), L-iditol 2-dehydrogenase (SORD, miR396b target), MACPF domain-containing protein CAD1 (miR168b target), succinate dehydrogenase (SDHA, miR397-5p target) and an uncharacterized protein encoding  (Table 1). Among these miRNAs, miR396a-3p, miR396b, and miR168b were up-regulated in rhizomes compared with other tissues, whereas miR390a-5p and miR396h were down-regulated in the rhizome. Additionally, the expression level of miR8175 was the highest in the root and that of miR397-5p which was the lowest in all tissues (Fig. 6A).
The abundance of abfA, GalDH, SORD, and SDHA transcripts was inversely correlated with that of their miRNAs (Fig. 6B), suggesting the possibility of a regulated mechanism. These genes may participate in polysaccharide biosynthesis in P. cyrtonema. . The X-axis shows the number of DEMs, and Y-axis shows the GO terms. GO terms grouped in biological process, cellular component, and molecular function are indicated in blue, yellow, and red, respectively

Validation of miRNA expression by qRT-PCR
To verify sRNA-Seq data, expression levels of seven miRNAs (miR396a-3p, miR396b, miR396g-3p, miR156a-5p, miR828a, miR160a, and miR169d) and their predicted target genes were verified by qRT-PCR using primers listed in Table S13. The results of qRT-PCR analysis of miR396a-3p, miR396g-3p, miR156a-5p, and miR-828a were consistent with those of sRNA-Seq analysis, and expression levels of these miRNAs exhibited a negative correlation with those of their putative target genes (Fig. 7, Tables S14 and S15). For example, compared with flowers and leaves, the expression of miR396a-3p, miR156a-5p and miR828 was significantly down-regulated in rhizomes, whereas that of their target genes were up-regulated.

Discussion
Polygonatum polysaccharides (PSPs) is an important biologically active substance of Polygonatum sibiricum Red, P. cyrtonema Hua, and Polygonatum kingianum Coll. et Hemsl. Previous studies have mainly focused on the transcriptome and pharmacological effects of this genus (Wang et al. 2017). In this study, we identified miRNAs and their target genes involved in the biosynthesis of PSPs. Our results provide valuable information on the molecular mechanisms underlying PSP biosynthesis. Four sRNA libraries from flower, leaf, rhizome, and root tissues were constructed and sequenced, generating more than 32 million clean reads. The 24 nt sRNAs exhibited the highest abundance, which was consistent with the distribution pattern of sRNAs reported in Panax notoginseng (Wei   (Fan et al. 2014). A total of 69 conserved miRNAs and 5 novel miRNAs were identified. Most of the miRNAs were highly conserved in plants, suggesting important regulatory roles of these miRNAs in plant growth and development. Analysis of the spatiotemporal expression patterns of miRNAs provides useful information on the molecular function of these miRNAs. Analysis of conserved miRNAs revealed that miR159a_1, miR168a-5p, and miR166g-3p were highly expressed in all tissues. Among these miR-NAs, miR159a_1 showed the highest expression level in leaves, which was consistent with expression pattern of Eucommia ulmoides (Wang et al. 2016). In rhizomes, expression levels of miR160, miR160e-5p, miR8175, miR319c_2, miR319b_1, miR319a_1, miR319a-3p, and miR319_1 were up-regulated and those of miR167d-5p, miR171b_2, miR397a_3, and miR160a-5p were significantly down-regulated compared with leaves and flowers. A previous study has shown that miR160a promotes leaf blade outgrowth as well as leaf and leaflet initiation and floral organ development through the quantitative regulation of its major target gene, SlARF10A (Damodharan et al. 2018). Additionally, miR160 plays a crucial role in local defense response and systemic acquired resistance (SAR) during the interaction between potato and Phytophthora infestans (Natarajan et al. 2018). The miR319a indirectly regulates the development of pistils by regulating its target gene, Pm-TCP4 . Additionally, miR319 and TCP gene families underlie robust and multilayer control of leaf development (Koyama et al. 2017). Similarly, miR396a-5p is an important conserved miRNA in plants; overexpression of miR396a-5p in tobacco increases tolerance to salt, drought, and cold stresses and susceptibility to Phytophthora nicotianae infection . In this study, novel miRNAs showed lower abundance than conserved miRNAs, suggesting that these novel species-specific miRNAs are either young miRNAs that arose recently through evolution or more volatile than other miRNAs.
The miRNAs exhibit various biological roles by regulating the expression of their target genes. The predicted target genes of DEMs identified in this study showed a wide range of functions, including transcription factor genes and functional gene families related to plant growth and development and metabolite biosynthesis. For instance, gene encoding the SPL transcription factor, a major plant-specific transcription factor family (Li et al. 2018;Zeng et al. 2019), was predicted as a target gene of miR156a-5p. A gene encoding MYB transcription factor was predicted as regulated by miR159a_1, miR858b, and miR319 families, which participate in flavonol biosynthesisand drought or salt stress tolerance Zhai et al. 2019) The ARF genes are a target of miR160 and miR167d, which regulate plant growth and abiotic stress response (Liu et al. 2018;Song et al. 2019). The target gene prediction of novel miRNAs revealed that the mannose-specific lectin precursor is regulated by novel_miR1. Genes encoding the low temperature-induced 65 kDa protein and ABC transporter C family member 5 was predicted as targets of novel_miR2, and those encoding peroxisomal adenine nucleotide carrier 1 and probable protein phosphatase 2C 38 were identified as targets of novel_miR4.
To identify miRNAs and their putative target genes related to polysaccharide biosynthesis, the target genes of DEMs involved in carbohydrate metabolism pathway were screened. Genes encoding abfA, GalDH, berberine bridge enzyme-like 18, SORD, MACPF domain-containing protein CAD1, an uncharacterized protein LOC105157180, and SDHA were predicted as targets of miR396a-3p, miR8175, miR396h, miR396b, miR168b, miR390a-5p, and miR397-5p, respectively. Moreover, PSPs mainly accumulate in the rhizome in P. cyrtonema, and miR396a-3p, miR396b, and miR168b were up-regulated in the rhizome compared with flower, leaf, and root. However, miR390a-5p and miR396h were down-regulated in the rhizome. Both abfA and SORD are key enzymes involved in carbohydrate metabolism, and these enzymes probably participate in polysaccharide biosynthesis in the rhizome of P. cyrtonema. Expression levels of seven tissue-specific miRNAs (miR156a-5p, miR-828, miR396a-3p, miR396g-3p, miR-160a, and miR-169d) and their putative targets were validated by qRT-PCR assays. The expression patterns of four miRNAs (miR156a-5p, miR-828, miR396a-3p, and miR396g-3p) were consistent with sRNA-Seq data and were inversely correlated with the expression pattern of their target genes. Previous studies have shown that miR156a-5p is involved in the regulation of the regulatory network involved in the response to viral infection in watermelon and other cucurbits (Sun et al. 2017)and stamen development in Viburnum macrocephalum f. keteleeri (Li et al. 2017). Additionally, miR-828 is one of the conserved miRNAs in plants, which targets homoeologous MYB2 genes that regulate trichome development in Arabidopsis and fiber development in cotton (Guan et al. 2014) and show strong associations with the color of tuber skin and tomato flesh (Bonar et al. 2018). Moreover, miR396g regulates genes associated with glucosylation and insolubilization of tannin precursors and the regulation of (de)astringency in persimmon fruits under normal development conditions (Luo et al. 2015). Another miRNA, miR396a, regulates growth-regulating factor1 (GRF1), salicylic acid carboxyl methyltransferase (SAMT), glycosyl hydrolases (GHs), and nucleotide-binding site-leucine-rich repeat (NBS-LRR). Additionally, miR396a-3p is critical for abiotic stress tolerance in tomato by affecting salicylic acid (SA) or jasmonate (JA) signaling pathway by the negative regulation of target genes and their downstream genes (Chen et al. 2017). Here, we found that miR396a-3p was down-regulated in roots and rhizomes of P. cyrtonema compared with leaves and flowers, while one of its target genes encoding abfA was up-regulated in rhizomes compared with other tissues. The abfA protein is an enzyme in the pentose and glucuronate interconversion pathway, suggesting its association with the biosynthesis of polysaccharides. Like Tx-Abf, abfA should be a retaining enzyme that catalyzes the hydrolysis of glycosidic bonds with Glu385 as the acid/base and Glu462 as the nucleophilea. Additionally, abfA (retaining enzyme) should be able to catalyze a wide variety of transglycosylation reactions.

Conclusion
This study is the first genome-wide investigation of miR-NAs and their target genes in P. cyrtonema. Four sRNA libraries were generated from flower, leaf, rhizome and root tissues of P. cyrtonema and sequenced using a highthroughput sequencing technology. A total of 69 conserved and 5 novel miRNAs were identified, and their expression patterns, target genes, and functions were characterized. In addition, seven miRNAs and their target genes involved in the accumulation of PSPs were identified.
Although the foundation of the complex miRNA-mediated regulatory networks remains unresolved, this miRNA data set will help elucidate the gene regulatory networks in P. cyrtonema and other species. Therefore, our results act as a valuable resource for studying complex gene regulatory functions of miRNAs, especially during the process of polysaccharide biosynthesis in P. cyrtonema.

Sample collection
P. cyrtonema plants were collected from the Chinese Medicine herb garden of Anhui University of Chinese Medicine with permission of managers and Professionals, and authenticated by Prof. Qingshan Yang. Flowers, leaves, rhizomes and roots were collected from three randomly selected plants.Tissues were cleaned with ultrapure water, dried on filter paper, and immediately frozen in liquid nitrogen for further analyses.

RNA isolation, sRNA library construction, and sequencing
Total RNA was isolated from flowers, leaves, rhizomes, and roots using RNA Plant Kit (Aidlab Biotech, Beijing, China), according to the manufacturer's instructions. RNA quality was evaluated using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Purified RNA was used to construct sRNA libraries, which were sequenced using the BGISEQ-500 platform at BGI Company, Wuhan, China.

Bioinformatics analysis and miRNA prediction
The bio-informatics analysis pipeline for small RNA sequencing was as follows: firstly, the low-quality reads, adapters and other contaminants were eliminated to get clean reads; Then, the length distribution analysis of the clean tags was conducted, and the clean tags were mapped to unigenes assembled by RNA denovo and other sRNA databases (Langmead et al. 2009), the small RNA was mapped and annotated to only one category followed the priority rule: miRbase [ pirnabank [ snoRNA (plant) [ Rfam [ other sRNA. The novel miRNA were predicted via miRDeep2 (Friedlander et al. 2008) and miRA (Evers et al. 2015) through exploring the characteristic hairpin structure of miRNA precursor.

Analysis of DEMs
The expression level of sRNAs was calculated using Transcripts Per Kilobase Million (TPM) (t Hoen et al. 2008), which eliminates the copy number difference caused by inconsistent PCR amplification during database construction. Therefore, the calculated gene expression data could be directly used to compare the difference in gene expression between samples. TPM ¼ C Â 10 6 =N where C represents the number of miRNA reads in the sample, and N represents the total number of mapped reads. The differential expression multiple of the gene between different samples was calculated based on the gene expression level. In our study, genes with fold change (FC) C 2.00 and false discovery rate (FDR) B 0.001 by the Poisson distribution method were considered differentially expressed (Audic and Claverie 1997;Li et al. 2019;Zhao et al. 2018a). Based on the differential miRNA expression results, the pheatmap package of the R software was used for hierarchical clustering analysis.

Target gene prediction, GO function annotation, and KEGG pathway analysis
The target genes of miRNAs were predicted using psRobot (Wu et al. 2012) and TargetFinder (Fahlgren and Carrington 2010). GO enrichment analysis finds all GO terms significantly enriched in a list of differential expression microRNAs (DEMs) target genes that correspond to specific biological functions. All genes were mapped to GO terms in the GO database (http://www.geneontology.org/), which calculates the number of genes associated with each term. The hypergeometric test was used to find significantly enriched GO terms in the input gene list. This test was based on GO-Term-Finder (http://www.yeastgenome. org/help/analyze/go-term-finder). The P-value was corrected using the Bonferroni method, and a corrected Pvalue B 0.05 was considered as a threshold. GO terms fulfilling this condition were defined as significantly enriched. Pathway-based analysis helps to further understand the biological function of DEMs target genes. KEGG was used to perform pathway enrichment analysis using the same formula for calculations as in GO analysis.
Validation of sRNA-Seq data using qRT-PCR To identify candidate genes involved in the regulation of polysaccharide biosynthesis, the expression level of six miRNAs and their potential target genes were selected for qRT-PCR. Total RNA was extracted from roots, rhizomes, leaves, and flowers using the TRIzol Reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions. Primers for qRT-PCR were designed using Primer 5.0 (Premier Biosoft International, Palo Alto, CA, USA) (Table S12). The qRT-PCR was conducted on a PIKOREAL 96 Real-Time PCR System using QuantiNova SYBR Green PCR Kit (Qiagen,USA). The b-actin and U6 with relatively stable expression in different tissues were used as internal references for mRNA and miRNA data normalization, respectively, and the raw data of qRT-PCR were listed in Table S14 and S15. Each PCR was performed in a reaction volume of 20 lL containing 10 ll of 2X SYBR Green mixture, 2 lL diluted cDNA, 0.5 lL forward primer, 0.5 lL reverse primer, and 7 lL double distilled water. The reaction conditions were as follows: 95°C for 2 min, followed by 40 cycles of 95°C for 5 s and 60°C for 10 s. Each sample was replicated at least three times. The relative expression level of each gene was calculated using the 2 -DDCt method. All data were subjected to analysis of variance (ANOVA), and the results were presented as mean ± SD. Differences with P \ 0.05 were considered to be statistically significant.