Integrative Analysis of mRNA and miRNA Expression Profiles of the Tuberous Root Development at Seedling Stages in Turnips

The tuberous root of Brassica rapa L. (turnip) is an important modified organ for nutrition storage. A better understanding of the molecular mechanisms involved in the process of tuberous root development is of great value in both economic and biological context. In this study, we analyzed the expression profiles of both mRNAs and miRNAs in tuberous roots at an early stage before cortex splitting (ES), cortex splitting stage (CSS), and secondary root thickening stage (RTS) in turnip based on high-throughput sequencing technology. A large number of differentially expressed genes (DEGs) and several differentially expressed miRNAs (DEMs) were identified. Based on the DEG analysis, we propose that metabolism is the dominant pathway in both tuberous root initiation and secondary thickening process. The plant hormone signal transduction pathway may play a predominant role in regulating tuberous root initiation, while the starch and sucrose metabolism may be more important for the secondary thickening process. These hypotheses were partially supported by sequential DEM analyses. Of all DEMs, miR156a, miR157a, and miR172a exhibited relatively high expression levels, and were differentially expressed in both tuberous root initiation and the secondary thickening process with the expression profiles negatively correlated with those of their target genes. Our results suggest that these miRNAs play important roles in tuberous root development in turnips.


Introduction
The modified tuberous root in many vegetables, such as radishes, carrots, turnips, and sweet potatoes provides a great food source for human consumption. Tuberous roots are rich in nutrition and dietary fiber [1], and convenient for storage. Some are also of medicinal importance [2][3]. Thus, it is important, in the context of economic value, to study the molecular mechanisms underlying the tuberous root formation and development. In additions, since tuberous roots are important modified organs, studies on tuberous root formation and development is also of high biological value.
greenhouse at 20 ± 2°C under a 16 h light/8 h dark cycle. The plants were watered once a week with the Hoagland solution. All the plants were grown in the same room until collected for sample preparation. Samples consisting of hypocotyl and main root tissues were collected on day 18 (the early stage before cortex splitting, ES), day 28 (the stage of cortex splitting, CSS) and day 64 (the stage of root thickening, RTS) after sowing, respectively. The samples collected from at least 10 seedlings at each stage were pooled together and stored in liquid nitrogen for RNA extraction. All samples were harvested at the same time during the day. To reduce error rates, two independent biological replicates were used for each stage. The samples were labeled as ES1, CSS1, RTS1 for the first biological replicate, and ES2, CSS2, RTS2 for the second biological replicate, respectively.

RNA isolation
Total RNA was extracted from each sample using Trizol reagent (Invitrogen, Carlsbad CA, USA). The quality of RNA was verified using a NanoDrop ND-1000 Spectrophotometer (Nano-Drop, Wilmington, DE, USA) and a 2100 Bioanalyzer RNA Nanochip (Agilent, Santa Clara, CA, USA), and the quantity of total RNA was determined using the NanoDrop ND-1000 Spectrophotometer. RNA samples were selected for further cDNA library and small RNA (sRNA) library construction according to the following selection criteria: OD260/280!1.8, OD260/230!1.8, RIN!9.0, 28S/18S!1.6, a concentration of !200 ng/μl, and total mass of !30 μg. For mRNA sequencing, the total RNA was cleaned up using the RNeasy Plant Mini Kit (74904; Qiagen, Valencia, CA, USA). For sRNA sequencing, RNA fragments of 18-30 nt were isolated and purified from total RNA samples extracted using Trizol reagent after electrophoresing on a denaturing 15% polyacrylamide gel.
Preparation, sequencing and data processing of mRNA libraries A total of six mRNA libraries were constructed in this study. The preparation and sequencing of the libraries were performed at Beijing Genomics Institute (BGI, Shenzhen, China) according to the manufacturer's protocol (Illumina). The sample libraries were single end sequenced via an IlluminaHiSeq TM 2000 (Illumina Inc., San Diego, CA, USA). Raw sequence reads were generated and filtered using the Illumina pipeline. Adapters, empty reads, and low quality reads (reads in which the percentage of unknown bases (N) greater than 10%, and reads in which the percentage of the low quality bases with a quality value 5 greater than 50%) were removed. The high-quality, clean sequences were mapped against the B. rapa (Chiifu-401) reference genome (http://brassicadb.org/brad/) using SOAPaligner/soap2 [34]. No more than two mismatches were allowed in the sequence alignment. The clean tags mapped to multiple reference genes were filtered, and the remaining tags were designed as unambiguous clean tags. Quality assessment of reads, statistics of alignment, sequencing saturation analysis and randomness assessments were carried out subsequently to assess the quality of sequencing.
Identification of DEGs among different seedling stages in the tuberous root of turnips [36]. The threshold with a FDR (False Discovery Rate) of 0.001 and the absolute value of log2Ratio ! 1 was used to assess the significant difference of expression levels of a given gene between two seedling stages. To reduce the false positive rates, only DEGs detected in both biological replicates were considered. To understand the biological functions of DEGs involved in tuberous root development of turnips, Gene Ontology (GO) enrichment analysis using the GO database (http://www.geneontology.org/), and pathway enrichment analysis using the KEGG database (http://www.genome.jp/kegg/) were conducted subsequently.
Preparation, sequencing and data processing of sRNA libraries A total of six sRNA libraries were constructed and sequenced at Beijing Genomics Institute (BGI, Shenzhen, China). The sRNA Libraries were prepared following the procedure by   [33], and were sequenced on an Illumina HiSeq TM 2000. Raw samples of the 49 nt RNA tags obtained from the HiSeq sequencing were cleaned by removing the low quality tags (reads with more than one base with SQ value 10 or more than two bases with SQ values 13), 5' adaptor contaminants, oversized insertions, no insert tags, poly A tags and small tags (reads shorter than 18 nt). Length distribution of the remaining clean tags was analyzed to determine the compositions of sRNA samples. Subsequently, we mapped the final clean reads against the Chiifu Chinese cabbage genome (http://brassicadb.org/brad/) using SOAP v1.11 to determine gene expression levels and distributions in the genome [34]. The annotation of the clean sRNA tags and prediction of novel miRNAs were performed following the procedure by   [33]. Putative target genes of the miRNAs were predicted using the rules developed by Allen et al. (2005) [37] and Schwab et al. (2005) [38].
The reads annotated as rRNA, scRNA, snoRNA, snRNA, and tRNA from GenBank and Rfam were removed, while those annotated as conserved miRNAs from the miRBase database and those predicted as novel miRNAs by Mireap (http://sourceforge.net/projects/mireap/) were used for further analysis in their differential expression profiles among different seedling stages in tuberous roots of turnips.

Identification of DEMs among different seedling stages in the tuberous root of turnips
The expression level of miRNA in all samples was normalized as transcripts per million (TPM, TPM = actual miRNA count / total count of clean reads × 10 6 ). The fold-change and P-value from the normalized expression level were calculated among different samples. The fold change was calculated as log2 Ratio, and the P-value was calculated using the formula described in   [33]. The threshold with the P-value 0.05 and an absolute value of log2Ratio ! 1 was used to evaluate the significance of the difference of miRNA expression levels. To reduce the false positive rates, only DEMs were detected in both biological replicates were considered.

Quantitative RT-PCR analysis
To confirm the identified DEGs and DEMs, DEMs and their differentially expressed target genes were subjected to the quantitative RT-PCR (qRT-PCR) analysis. For DEM validation, the first strand cDNA was synthesized using the miRcutemiRNAcDNA Synthesis Kit (DP501, Tiangen, Beijing, China), and qRT-PCR analyses were performed using the miRcutemiRNAqPCR Detection Kit (SYBR Green) (FP401, Tiangen, Beijing, China). Spliceosomal RNA U6 was used as a constitutive expression control. PCR reactions were performed in a 96-well plate at 94°C for 2 min, followed by 45 cycles of reaction (94°C for 20 s, followed by 60°C for 34 s). For DEGs validation, the first strand cDNA was synthesized using M-MLV Reverse transcriptase (M1701, Promega, Madison, USA) and qRT-PCR analyses were performed using a PrimeScript™ RT reagent Kit (Perfect Real Time) (DRR037A, Takara, Dalian, China). PCR reactions were performed in a 96-well plate at 95°C for 10 min, followed by 40 cycles of reaction (95°C for 10 s, followed by 60°C for 60 s). Actin gene was used as a constitutive expression control. All qRT-PCR analyses were performed using an IQ5 Real-Time PCR System (BIO-RAD, Hercules, CA, USA). The specific primers for qRT-PCR are described in S1 Table. Results mRNA Library construction, sequencing and mapping against the B. rapa (Chiifu-401) reference genome To study the global transcriptome during the tuberous root development of turnips, a total of six mRNA libraries (three seedling stages during tuberous root development, two biological replicates for each stage) were constructed and sequenced using an IlluminaHiSeq TM 2000. Approximately 12 million raw tags were generated for each library, and more than 99.0% of the raw tags were identified as clean tags (Table 1). We then mapped all the clean tags against the B. rapa reference genome (http://brassicadb.org/brad/), and found that over 71% of the clean tags for each library could be mapped onto the reference genes with 2 bp mismatch. Less than 29% of the clean tags for each library could not be aligned to any reference genes due to incomplete sequences, and these tags were designated as unknown tags. More than 67.3% (about 8 million) of the clean tags in each library were mapped to a single gene, while less than 4.7% (about 0.5 million) were mapped to multiple reference genes (Table 1). Unknown tags and tags mapped to multi-genes were filtered out, and the unique clean tags mapped to a single gene were reserved for further DEG analysis.
The saturation for the six libraries was analyzed, and the results showed that the number of detected genes increased with the number of reads until the number of reads reached up to 2 million (S1 Fig), which indicates that the sequencing depth is sufficient for the transcriptome coverage in this study. We also assessed the randomness of RNA fragmentation in the six libraries to estimate whether it was sufficient for a subsequent bioinformatics analysis of gene expression, and the results showed that the reads in each position of the reference gene were distributed evenly and demonstrated highly similar tendencies in all libraries (S2 Fig).

Differentially regulated genes in tuberous root development of turnips
Based on deep sequencing of the six mRNA libraries, approximately thirty thousand genes (about 70% of reference genes in Chinese cabbage) were detected in each library (Table 1, S2  Table). To compare gene expression differences in different stages of tuberous root development in turnips, the gene expression level was calculated using the RPKM method [35]. Gene expression correlations based on RPKM between the two biological replicates were analyzed, and the correlation coefficients (R 2 ) were high (0.97, 0.95, and 0.96 for ES, CSS, and RTS, respectively) ( Fig 1A). Genes expressed in both replicates in the same stage were screened. 28,798, 28,315, and 27,288 genes were expressed in both replicates in ES, CSS, and RTS, respectively. Most of them (25,729) were expressed in all three stages, while only a small portion (460-1,365) was expressed at two stages or specifically at one stage ( Fig 1B).
We compared the gene expression profiles between every two stages in tuberous root development in both biological replicates (ES1 vs. CSS1, ES1 vs. RTS1, CSS1 vs. RTS1, ES2 vs. CSS2, ES2 vs. RTS2, and CSS2 vs. RTS2, where the former is the control and the latter is the experimental group). The results are presented in S2 and S3 Tables. DEGs detected in both biological replicates were screened for subsequent analysis. A total of 2693 (1176 up-and 1517 down-regulated), 4655 (1094 up-and 3561 down-regulated), and 2133 (436 up-and 1677 down-regulated) genes were differentially expressed in both replicates in the stage pairs, ES vs. CSS, ES vs. RTS, and CSS vs. RTS, respectively ( Fig 1E). More genes, which were down-regulated, were detected than the up-regulated ones in all pairs ( Fig 1E).
Gene Ontology (GO) enrichment analysis of the DEGs was performed to gain insight into the possible mechanisms involved in tuberous root development. All the well-annotated DEGs were functionally clustered into three main groups (cellular component, molecular function, and biological process). For cellular component ontology, a large number of the DEGs in the stage pair, ES vs. CSS, involved in cell, cell part, and organelle, accounted for 63.73%, 63.73%, and 44.15% of all the well-annotated DEGs, respectively. Under the molecular function category, catalytic activity (50.26%) and binding (48.62%) were the main groups, followed by nucleic acid binding transcription factor activity (10.54%) and transporter activity (6.88%) in ES vs. CSS. For biological process category, the main groups of the DEGs in ES vs. CSS were involved in metabolic process (53.74%), cellular process (51.42%), response to stimulus (39.89%), and single-organism process (32.31%). Similar results were obtained for GO enrichment analysis of the DEGs in ES vs. RTS, and CSS vs. RTS (Fig 2).
We performed pathway enrichment analysis by mapping all the annotated genes to terms in the KEGG database to identify significantly altered biological pathways associated with the tuberous root development in turnips. The results are shown in Table 2. A total of 17, 20, and 21 pathways were significantly enriched with Q value 0.05 in ES vs. CSS, ES vs. RTS, and CSS vs. RTS, respectively. Of these, 11 pathways were significantly enriched in all comparison pairs, indicating that these pathways may be involved in both tuberous root initiation and the secondary thickening process. Four pathways, including plant hormone signal transduction, amino sugar and nucleotide sugar metabolism, other glycan degradation, and cutin, suberine and wax biosynthesis, were only enriched in ES vs. CSS, implying the great importance of these pathways in tuberous root initiation. The plant hormone signal transduction pathways were significantly enriched in ES vs. CSS. Among these seven phytohormones, namely auxin, CTK, brassinosteroid (BR), gibberellin (GA), JA, ethylene, and ABA, signal transduction pathway of auxin had the most DEGs (S4 Table). Metabolic pathways, biosynthesis of secondary metabolites, and plant hormone signal transduction were the top three enriched pathways in ES vs. CSS, while metabolic pathways, biosynthesis of secondary metabolites, and starch and sucrose metabolism were the top three enriched pathways in both ES vs. RTS and CSS vs. RTS (Table 2).    Table 3). The percentages of mapped total clean sRNAs and unique sRNAs were similar in all of the six sRNA libraries.

Identification of conserved and novel miRNAs in the tuberous root of turnips
The clean sRNAs were annotated following the procedure described in Wang et al. (2013) [33], and the distribution of sRNAs among different categories is shown in S5  Table). We aligned all these miRNA candidates to the miRNA precursor/mature miRNA of all plants deposited in the miRBase database to identify conserved miRNAs. Only perfectly matched or closely related ( 2 mismatches) sequences and precursors that could form hairpins were defined as conserved. A total of 173, 187, 162, 165, 179, and 157 conserved miRNAs were identified from the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively. In all, 279 conserved miRNAs were identified in this study, among which, 98 were expressed in every library (S6 Table). Conserved miRNAs expressed in both replicates in the same stage were screened for further analysis. It was found that 139, 149, and 125 were expressed in both replicates in ES, CSS, and RTS, respectively ( Fig 1C). One hundred of the total detected conserved miRNAs were expressed in all three stages, while 2-18 were only expressed at two stages or specifically at one stage ( Fig 1C). Currently, 43 miRNAs (35 miRNAs, 4 miRNA-5p and 4 miRNA-3p) in B. rapa have been deposited in the miRBase database (miRNA Registry, Release 20.0, 2013). In this study, 38, 36, 31, 39, 35, and 31 bra-miRNAs were identified from the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (S7 Table). Collectively, 42 of these 43 bra-miRNAs deposited in the miRBase database, except bra-miR5723, were identified in at least one of the six libraries; 27 bra-miRNAs were found in each of the six libraries. However, over half (23-26) of these 43 bra-miRNAs (e.g. bra-miR5720, bra-miR5717, and bra-miR5724) showed a very low level (less than 100 reads) or no expression in all libraries, while only 3-6 (e.g. bra-miR157a and bra-miR164a) were expressed dramatically high (more than 5000 reads) in all or most of the libraries (S7 Table). This finding indicates that the expression of conserved miRNAs varies significantly in the tuberous root of turnips. Bra-miRNAs expressed in both replicates in the same stage were screened, and 36, 34, and 29 were expressed in both replicates in ES, CSS, and RTS, respectively (S7 Table).
Novel miRNAs were predicted by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated sRNA tags which could be mapped to the genome using the Mireap program developed by BGI (http://sourceforge.net/projects/mireap/ ). Accordingly, 292, 292, 279, 266, 287, and 266 novel miRNAs were identified in the ES1, CSS1, RTS1, ES2, CSS2, and RTS2 libraries, respectively (S6 Table). Most of the novel miRNAs (more than 85%) showed very low expression level (less than 100 reads), and only three (novel_mir_14, novel_mir_176, and novel_mir_26) showed an expression level higher than 1000 reads in all six libraries. A total of 912 novel miRNAs were identified, but only 48 were expressed in all six libraries (S6 Table), indicating significant variation of novel miRNA expression in different samples. A total of 119, 128, and 121 were expressed in both replicates in ES, CSS, and RTS, respectively (Fig 1D). Venn diagram analyses of stage specific novel miRNAs in tuberous root of turnips are shown in Fig 1D. Differentially regulated miRNAs in tuberous root development of turnips miRNA expression correlations between the two biological replicates were analyzed according to TPM. The correlation coefficients (R 2 ) were 0.98, 0.98, and 0.98 for ES, CSS, and RTS, respectively (Fig 1A). Based on the criteria described in the Methods section, DEMs between different stages in each replicate were identified, and the detailed results are listed in S6 Table. DEMs detected in both the biological replicates were screened for subsequent analysis. In all, 46 (29 up-and 17 down-regulated), 50 (22 up-and 28 down-regulated) and 36 (11 up-and 25 down-regulated) conserved miRNAs (Fig 1F), and 53 (26 up-and 27 down-regulated), 62 (26 up-and 36 down-regulated) and 40 (12 up-and 28 down-regulated) novel miRNAs (Fig 1G) were differentially expressed in both replicates in ES vs. CSS, ES vs. RTS and CSS vs. RTS, respectively. Eight miRNAs (miR156a, miR157a, miR157a-3p, miR5718, novel_mir_4, novel_-mir_85, novel_mir_142, and novel_mir_236) were differentially expressed in both replicates of all comparison pairs, indicating that these miRNAs may be involved in both tuberous root initiation and the secondary thickening process. Twenty miRNAs were differentially expressed in ES vs. CSS, implying a potential role of these miRNAs in tuberous root initiation (S6 Table).
We  Table). GO enrichment analyses of the target genes of DEMs (both conserved and novel DEMs) detected in both replicates were conducted subsequently. For cellular component ontology, cell, cell part, and organelle were the main groups of the target genes; under the molecular function category, binding, nucleic acid binding transcription factor activity, and catalytic activity were the main groups; for the biological process category, cellular process, metabolic process, and single-organism process were the top three groups ( Fig  2B). The results of GO analysis in the target genes of DEMs were similar to those of DEGs in this study (Fig 2).
Pathway enrichment analyses of the target genes of DEMs detected in both replicates were also performed using the KEGG pathway database. No pathway was significantly enriched with Q value of 0.05 in all the three comparisons, which may be because of the small numbers of the target genes. More genes were enriched in plant hormone signal transduction pathway in ES vs. CSS and ES vs. RTS, while more were enriched in metabolic pathways in CSS vs. RTS (S9 Table).
The expression profiles of the target genes of DEMs were analyzed according to the highthroughput sequencing of mRNA. The results showed that a number of DEMs had target gene (s) differentially expressed in the DEG analysis. Some differentially expressed target genes had their expression profiles matched up with those of their corresponding DEMs (S8 Table). These DEGs and DEMs are the most important candidates for further study of their functions in tuberous root development of turnips.

Validation of DEGs and DEMs by qRT-PCR
The expression profiles of six DEMs between different stages and their seven differentially expressed target genes were verified by qRT-PCR. The results showed that all seven genes and four miRNAs (miR156a, miR157a, miR172a, and novel-mir-208) exhibited the same  (Fig 3, S2 and S6 Tables). The expression profile of miR395b did not match with that of the original results (Fig 3B, S6 Table), probably due to its low expression level. Similarly the expression profile of novel-mir-108 did not match either. Of the four miRNAs exhibiting the same expression tendency as the original results, miR156a, miR157a, and miR172a exhibited relatively high expression level in all samples, while novelmir-208 expressed extremely low in all samples in the sRNA sequencing.

Quality assessment of sequencing
Tuberous root development is a complex biological process controlled by interactions of genetic, physiological and environmental factors. High-throughput sequencing, owing to its advantages of high coverage and low cost, has been widely used to study complex biological processes of plants in recent years. Here, in order to gain insight into the underlying molecular mechanisms regulating tuberous root development, large-scale analyses of both mRNA and miRNA expression profiles of tuberous roots at different seedling stages in turnips were performed by high-throughput sequencing technology. Good repeatability was achieved, according to the high correlation coefficients (R 2 >0.95, Fig 1A) between the two replicates on both mRNA and miRNA expression levels. The number of raw reads, the percentages of clean tags, and the percentages of clean tags mapped onto the reference genome in both mRNA and miRNA expression profile analyses were greater than, or similar to, those in previous studies in the brassica species [33,[39][40][41][42][43][44]. According to the results of sequencing saturation analysis and randomness assessment, the clean tags obtained in this study were sufficient for subsequent bioinformatics analyses of gene expression. The size distribution of sRNAs in each sRNA

Genes and biological pathways involved in tuberous root development of turnips
Genome-wide gene expression profiles analysis, based on high-throughput RNA sequencing, could facilitate the revelation of quantitative changes in transcript abundance and the identification of candidate genes. We identified a large number of DEGs among different developing stages of tuberous roots of turnips (Fig 1E, S2 Table), implying a high complexity of molecular regulations in tuberous root development. Tuberous root is a nutrition storage organ. Sucrose and starch, proteins and amino acid, dietary fiber, vitamin, and other relevant molecules are accumulated in the tuberization process [1,11,13]. Many metabolic pathways may be actively involved in this process. Based on our genome-wide DEGs analysis, metabolic and biosynthesis of secondary metabolite pathways were the top two pathways with a large number of DEGs significantly enriched, which is consistent with previous studies on the tuberous root of sweet potatoes [16,[21][22][23], radishes [14], and cassavas [20]. Of these metabolism pathways, starch and sucrose metabolism pathway and phenylpropanoid biosynthesis pathway are the top two evident pathways, which may be active in both tuberous root initiation and the secondary thickening process. In cassavas, glycolysis/ gluconeogenesis was identified to be the most evident pathway involved in tuberous root development [20]. In turnips, this pathway was probably not so actively involved, because we did not see a significant enrichment.
Phytohormones play important roles in the regulation of tuberous formation and development [4][5][6][7]. Auxin [4][5], CTK [4][5]7], and JA [7] were reported to be crucial regulators for tuberous initiation, while auxin [4], CTK [6], and ABA [6] probably regulate the secondary thickening growth of tuberous roots. In our study, plant hormone signal transduction pathways were significantly enriched in ES vs. CSS, but not in ES vs. RTS or CSS vs. RTS (Table 2). This suggests that DEGs related to plant hormone signal transduction may play a more important role in tuberous root initiation than in the secondary thickening growth of tuberous roots in turnips. Besides DEGs involved in auxin, CTK and JA signal transduction pathways, DEGs involved in BR, GA, and ethylene signal transduction pathways were also identified in ES vs. CSS (S4 Table), indicating potential roles of these plant hormones in regulating tuberous root initiation in turnips. Furthermore, over half of the DEGs involved in plant hormone signal transduction pathways in ES vs. CSS belong to auxin and CTK pathways, implying the predominant roles of these two phytohormones in regulating tuberous root initiation in turnips.
Several genes related to tuberous development have been identified in previous studies. MADS box genes, such as IbMADS1, SRD1, and IbAGL17, were identified to be preferentially expressed in tuberous root of sweet potatoes [4,8,48]. In this study, four homologous of MADS box genes (Bra017340, Bra004928, Bra006051, and Bra009913) were identified to be possibly involved in tuberous root development in turnips. Bra017340 was more predominantly expressed in CSS libraries than in ES and RTS libraries, indicating a possible positive role of Bra017340 in tuberous initiation process. Bra004928 and Bra006051 were up-regulated in ES vs. RTS, which are likely to be positively involved in the secondary thickness process of tuberous root development in turnips. Bra009913 were down-regulated in both ES vs. RTS and CSS vs. RTS, implying its negative role in regulating tuberous root development. Expansin (EXP1) has been reported to play a negative role in the formation of storage root in sweet potatoes [10]. In our study, the homologous gene EXP1 (Bra007864) and five other Expansin genes (Bra012218, Bra036225, Bra028878, Bra025372, and Bra036499) were identified to be downregulated in both ES vs. CSS and ES vs. RTS. Previous studies reported SuSy plays a key role in regulating tuberous root development [12], and homologous of the enzyme are up-regulated in tuberous roots of radishes [14]. In this study, one homologous of SuSy 3 was up-regulated, and two homologous of SuSy 6 were down-regulated in ES vs. RTS. No SuSy homologous genes were identified differentially expressed in ES vs. CSS. This indicates that SuSy genes probably are involved in tuberous root thickness, but not in tuberous root initiation. New members of MADS box genes, expansin genes and SuSy were identified for tuberous root development, indicating that high-throughput RNA sequencing is effective in identifying DEGs. miRNAs involved in tuberous root development of turnips miRNAs are important regulators of many important biological processes, such as development, hormone signaling, abiotic stress responses, and pathogen responses [49][50][51]. They negatively regulate their target gene expression at the post-transcriptional level. Large-scale miRNA expression profiles were analyzed on the modified leafy head development of Chinese cabbage, providing a relatively complete view of the miRNAs involved in the modified leafy head development of B. rapa [33]. Whether miRNAs are involved in the modified tuberous root development in B. rapa is still unclear. According to our results, many miRNAs were found differentially expressed in different development stages of tuberous roots in turnips [S6 Table]. In previous studies, it was reported that miR156 and miR172 are involved in tuberous stem development process in potatoes. Overexpression of miR156 reduces the yield of tuberous stem [29], while miR172 overexpression accelerates the tuberization in potatoes [30]. miR156 down-regulates the expression of miR172 by targeting one of the squamosa promoter binding protein-like (SPL) genes, StSPL9, a transranscription factor of miR172 [29]. In this study, we found miR156a was down-regulated in the root tuberization process, while miR172a was upregulated in CSS and RTS compared with ES, which is coincide with that in stem tuberization process in potatoes. miR157 families are predicted to target the SPL genes at the same binding sites as miR156 families [52]. In our study, similar to miR156a, the expression of miR157a was down-regulated during the tuberization process. They are likely to play the same roles in regulating tuberous root development in turnips. Furthermore, miR156a, miR157a and miR172a exhibited a relatively high expression level in all samples, and their expression patterns were negatively correlated with all or some of their target DEGs. It indicated that these three miRNAs may play important roles in tuberous root development of turnips. Further studies need to be done to fully understand their functions in regulating tuberous root development.
Gene Ontology analyses of the target genes of the DEMs were similar with those of the DEGs (Fig 2). Although no KEGG pathway was significantly enriched due to the small numbers of the target genes, more genes were clustered in the plant hormone signal transduction pathway in ES vs. CSS and ES vs. RTS, and more were enriched in metabolic pathways in CSS vs. RTS (S9 Table). It supports our hypothesis that plant hormone signal transduction pathway plays more important roles in tuberous root initiation.
Since miRNAs negatively regulate gene expression at the post-transcriptional level, the expression profiles of DEMs are expected to be negatively correlated with the target gene expression profiles. Previous studies have shown that some miRNAs were co-expressed with their target genes [53][54][55]. We analyzed the expression profiles of the putative target genes of DEMs based on the DEGs sequencing, and found that many miRNA expression profiles were not correlated with their target gene expressions. Similar results were obtained in Chinese cabbage for the modified leafy head development [33]. The results indicate that target gene expression is not regulated by miRNAs alone, but many other factors.
qRT-PCR validation of the six DEMs showed that four DEMs (miR156a, miR157a, miR172a, and novel-mir-208) exhibited the same expression tendency as the original results, while two (miR395b, and novel-mir-108) did not match with those of the original results from the high-throughput sRNA sequencing (S6 Table, Fig 3B). This is probably due to the extremely low expression levels (less than 300 reads were detected in each sample) of the two miRNAs detected in the sRNA sequencing.

Conclusions
We analyzed both mRNA and miRNA expression profiles of the tuberous root of turnips at different stages based on high-throughput sequencing technology. Our results indicate that a large number of DEGs and DEMs are involved in tuberous root development. According to the DEG analysis, we propose that metabolism is the dominant pathway in both tuberous root initiation and secondary thickening process; the plant hormone signal transduction pathway may play a more important role in regulating tuberous root initiation; and the starch and sucrose metabolism may be more important for the secondary thickening process. Our analyses for DEMs partially supported these hypotheses. Based on the DEM analyses, miR156a, miR157a, and miR172a may play important roles in tuberous root development in turnips. These findings will contribute to our further understanding of the molecular mechanisms of tuberous root development in turnips.