Genome-wide sRNA and mRNA transcriptomic profiling insights into dynamic regulation of taproot thickening in radish (Raphanus sativus L.)

Background Taproot is the main edible organ and ultimately determines radish yield and quality. However, the precise molecular mechanism underlying taproot thickening awaits further investigation in radish. Here, RNA-seq was performed to identify critical genes involved in radish taproot thickening from three advanced inbred lines with different root size. Results A total of 2606 differentially expressed genes (DEGs) were shared between ‘NAU-DY’ (large acicular) and ‘NAU-YB’ (medium obovate), which were significantly enriched in ‘phenylpropanoid biosynthesis’, ‘glucosinolate biosynthesis’, and ‘starch and sucrose metabolism’ pathway. Meanwhile, a total of 16 differentially expressed miRNAs (DEMs) were shared between ‘NAU-DY’ and ‘NAU-YH’ (small circular), whereas 12 miRNAs exhibited specific differential expression in ‘NAU-DY’. Association analysis indicated that miR393a-bHLH77, miR167c-ARF8, and miR5658-APL might be key factors to biological phenomenon of taproot type variation, and a putative regulatory model of taproot thickening and development was proposed. Furthermore, several critical genes including SUS1, EXPB3, and CDC5 were characterized and profiled by RT-qPCR analysis. Conclusion This integrated study on the transcriptional and post-transcriptional profiles could provide new insights into comprehensive understanding of the molecular regulatory mechanism underlying taproot thickening in root vegetable crops.


Background
Radish (Raphanus sativus L., 2n = 2x = 18) is an important worldwide root vegetable crops belonging to Brassicaceae family. The fleshy taproot is the main product organ and determines the final yield and quality of radish. Abundant nutrient substances exist in fleshy taproot including carbohydrates, crude fiber, vitamin C and protein. Currently, extensive researches on the molecular mechanism of root development had been conducted in a range of plant species such as Arabidopsis [1], tobacco [2], maize [3][4][5], and rice [6,7]. However, the molecular mechanism underlying taproot thickening is still far away from being fully clarified in root vegetable crops such as radish.
In the past decades, with the rapid development of 'omics' methodology, RNA-seq has become a valuable strategy for systematical identification of differentially regulated genes, miRNAs, and regulation pathways in different tissues, organs, and developmental stages in several plant species, such as Rosa chinensis [8], Glycine max [9], Citrus sinensis [10], Myrica rubra [11], Solanum lycopersicum [12], and R. sativus [13]. Furthermore, the available genome database of radish provided a useful genome information platform for the investigation of molecular mechanism underlying radish taproot thickening [14][15][16].
In recent years, several studies on identification and dissection of gene expression and complex regulatory network during taproot thickening in radish have been performed. Using RNA-Seq technology, many miRNAs and transcripts were identified to be differentially expressed during radish taproot thickening [13,17], and carbohydrate metabolism pathway was significantly activated during taproot thickening, particularly in cell proliferating tissues [15]. Nevertheless, previous studies were mainly focused on a single radish cultivar, and the key genes involved in radish taproot thickening among different root-type radish genotypes remains to be accurately identified. Increasing evidences indicated that root shape and size (root-type) are significantly different among different cultivars, and genes related to root-type difference are considered to be candidates that promote or repress taproot thickening in radish [15,18,19]. So far, there is no report on identification of genes involved in taproot development in different cultivars, which limits the genetic improvement and germplasm innovation of radish cultivars.
To identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) involved in taproot thickening, three advanced inbred lines with different root shape and size consisting of radish 'NAU-DY' (large acicular), 'NAU-YB' (medium obovate) and 'NAU-YH' (small circular) were used in this study. An integrated mRNA-seq and sRNA-seq analysis was performed at three development stages of taproot thickening from three advanced inbred lines. Furthermore, RT-qPCR analysis was carried out to validate the expression patterns of several important candidate genes. The outcome of this study could reveal critical genes and miRNAs underlying the taproot thickening, and provide new insights into dynamic regulation of taproot thickening in radish.

Plant materials
Three advanced inbred lines of radish, 'NAU-DY' (large acicular), 'NAU-YB' (medium obovate), and 'NAU-YH' (small circular), were used in this study, and the seeds were developed from college of Horticulture, Nanjing Agricultural University, Nanjing, China. (Fig. 1). The seeds were germinated on a wet filter paper in darkness at room temperature for 3 days, and then planted in plastic pots and cultured in the growth chamber with 16 h light (25°C) and 8 h dark (18°C). For each advanced inbred line, the characteristics of radish growth and development were shown in Fig. 1, and the time point with 'two leaves and one heart', cortex splitting, and the highest rate of taproot thickening were as selection criteria of pre-cortex splitting stage (PSS), cortex splitting stage (CSS) and expanding stage (ES), respectively. The taproot samples S1, S2, and S3 were harvested at PSS
Small RNA library was generated from total RNA according to the instruction of NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) in 'NAU-DY'. The detailed experimental operation procedures of sRNA library construction and sequencing followed the previous method [22], and the corresponding sRNA-seq libraries were named as DS1, DS2, and DS3 library, respectively. Among these, mRNA-seq and sRNA-seq data of 'NAU-YH' were cited from our previous studies [Sequence Read Archive (SRA) with the GenBank accession No.: SRX707630] [13,17]. The technical workflow of the integrated mRNA-seq and sRNAseq analysis was shown in Fig. 2.

Genome mapping and differential expression analysis
Reference genome sequences were downloaded from the radish genome website (ftp://ftp.kazusa.or.jp/pub/radish), and two transcriptome sequences of radish from 'NAU-YH' and 'CKA' were download from our two published papers, and all these sequences were acted as reference sequences in this study [23,24]. Prior to the alignment of RNA-seq reads to reference sequences, the raw reads were screened to achieve high quality clean reads by removing adaptor reads and contaminants. For mRNA-seq data analysis, the clean reads were mapped to reference sequences of radish with no more than two mismatchs using TopHat2 software [25]. The FPKM method was used to calculate gene expression level. Read counts data were normalized using TMM method, and p-value was calculated with Poisson distribution model. The FDR (false discovery rate) is determined by p-value ranges in multiple tests. In this study, the threshold of |log 2 FC (fold change)| > 1 with q-value < 0.005 was selected as simulated biological variation to determine whether a gene is significantly differential expression in the DEGseq analysis [26,27] For sRNA-seq data analysis, the unique small RNAs were mapped to radish genomic sequence by Bowtie software [28]. Sequences matching non-coding RNAs included rRNAs, tRNAs, snRNAs, and snoRNAs were removed. The remaining unique sequences were searched against with known miRNA sequence by miRBase 21 software for known miRNA identification. Then the remaining unknown sRNAs was used to predict novel miRNAs by miREvo [29] and mirdeep2 software [30]. Differentially expressed miRNAs (DEMs) from different developmental stages of taproot thickening were identified using DESeq software [31]. For all comparisons, miRNAs with |log 2 FC | > 1 and an adjusted q-value < 0.01 were assigned as DEMs.

miRNA target prediction and annotation
Prediction of miRNA target gene was performed by psRobot_tar in psRobot [32]. GO classification (GO database, http://www.geneontology.org.) and KEGG pathway (KEGG database, http://www.genome.jp/ kegg/) methods were carried out for allocating genes to different functional categories and predicting their biological functions, respectively. The KEGG pathway enrichment and GO enrichment analysis were performed with the condition of corrected p value < 0.05.

Reverse transcription quantitative PCR (RT-qPCR) analysis
Total RNA and miRNA extraction (Tiangen) and reverse transcription (Takara) were conducted according to the manufacturer's instructions. RT-qPCR was performed using a SYBR Primix Ex Taq kit (TaKaRa), and the amplification reactions were conducted on ROCHE LightCycler 480 instruments [33]. The RsActin and 5.8S ribosomal RNA (rRNA) were used as the reference genes for normalization, respectively. The relative expression level of each gene was calculated by 2 -△△CT method. Three replicates and Duncan's test (P < 0.05) were conducted, and the Pearson correlation coefficient was calculated by DPS software to evaluate the correlation of gene expression patterns from RNA-Seq and RT-qPCR. Primers were designed by Beacon Designer 7.0 (Additional file 1: Table S1).  reference sequences corresponding to DS1, DS2, and DS3 library, respectively (Additional file 1: Table S2).

Analysis of sRNA and mRNA sequencing data
To identify regulated genes during taproot thickening in radish, six cDNA libraries from two advanced inbred lines were constructed from plants at three developmental stages of taproot thickening (S1, S2, and S3).  Table S3).

DEMs identification during three developmental stages of taproot thickening
A total of 77, 85, and 56 DEMs were identified from DS2 vs DS1, DS3 vs DS1, and DS3 vs DS2, respectively ( Fig. 3a, b, c and Additional file 1: Table S4). In detail, compared with pre-cortex splitting stage, 28 upregulated miRNAs and 49 down-regulated miRNAs were identified at cortex splitting stage. Meanwhile, 85 DEMs including 26 up-regulated and 59 down-regulated miR-NAs were identified in DS3 vs DS1. However, only 56 miRNAs were identified to be differentially expressed in DS3 vs DS2, with 19 up-and 37 down-regulated miR-NAs in radish. Interestingly, only eight, seven and three miRNAs were specifically expressed in DS2 vs DS1, DS3 vs DS1, and DS3 vs DS2, respectively; whereas 33, 17 and eight DEMs were shared with each pairwise comparison, respectively; but 28 DEMs were shared among three comparisons (Fig. 3d). Heatmap clustering of DEMs was shown in Additional file 1: Figure S1, and the results indicated that the expression levels of those miR-NAs exhibited characteristics of dynamic change during radish taproot thickening.
A total of 16 DEMs were shared between 'NAU-YH' and 'NAU-DY' from sRNA-seq data ( Table 1). Among these DEMs, 14 DEMs were conserved miRNAs, and two DEMs were non-conserved miRNAs. Interestingly, the expression patterns of miR394a, miR408-5p, and miR828 were similar between 'NAU-YH' and 'NAU-DY'; whereas that of miR395a was changed between 'NAU-YH' and 'NAU-DY' in two lines of three corresponding comparisons (Table 1). Meanwhile, compared with the small-size genotype 'NAU-YH', a total of 12 DEMs were specifically differential expressed during taproot thickening in large-size radish genotype 'NAU-DY' ( Table 2). Among these, miR165a-3p and miR165a-5p were downregulated in DS3 vs DS1 and DS3 vs DS2 pairs, and miR167c-5p and miR167d were down-regulated in DS2 vs DS1 and DS3 vs DS1 pairs; whereas miR167a-3p were up-regulated in each comparison pairs. These results suggested that miR394a, miR408-5p, and miR828 might be involved in taproot thickening, while miR395a and DEMs specific to 'NAU-DY' might contribute to roottype differences.
GO enrichment analysis was performed among 2606 DEGs that shared according to the same comparison pairs from 'NAU-YB' and 'NAU-DY'. 'cellulose microfibril organization' (GO:0010215), 'cell growth' (GO: 0016049), 'carbohydrate biosynthetic process' (GO: 0016051), and 'S-adenosylmethionine biosynthetic process' (GO:0006556) were specifically enriched in S3 vs S1, whereas 'carbohydrate metabolic process' (GO: 0005975) was shared with each comparison pairs in taproot thickening of radish (Table 3). Meanwhile, KEGG enrichment analysis was used to identify the critical pathway that genes involved in taproot thickening. 'starch and sucrose metabolism' (ath00500) was shared among three comparison pairs, 'phenylpropanoid biosynthesis' (ath00940), 'glucosinolate biosynthesis' (ath00966) was shared between two comparison pairs (S2 vs S1 and S3 vs S1), whereas 'thiamine metabolism' (ath00730) was specifically enriched in S3 vs S2 ( Table 4). The results suggested that it was a process of substances and energy metabolism, which promotes cell growth and organ enlargement during taproot thickening in radish.
A total of 12 DEMs were specifically differential expressed during taproot thickening in large-size radish genotype 'NAU-DY' when compared with the small-size genotype 'NAU-YH' ( Table 2). To explore miRNA-mRNA regulatory network in taproot thickening, the association analysis between miRNA and mRNA of 12 DEMs that specific to 'NAU-DY' was performed. The results showed that miR167c-5p (targeted by ARF8), miR393a-5p (targeted by bHLH77), miR5658 (targeted by APL) were identified to be specifically differentially expressed during radish taproot thickening of 'NAU-DY' (Additional file 1: Table S6).
To verify the expression patterns of miRNAs and their corresponding targets in radish, a total of 14 DEMs and three differentially expressed target genes were selected for RT-qPCR analysis. As shown in Figure S4, the expression patterns of RT-qPCR analysis were in line with Table 1 Overview of DEMs shared between 'NAU-YH' and 'NAU-DY' DEMs NAU-YH NAU-DY log 2 (S2/S1) log 2 (S3/S1) log 2 (S3/S2) log 2 (S2/S1) log 2 (S3/S1) log 2 (S3/S2)  that obtained from sRNA-seq, whereas they had differences on the degrees of differential expression between the two methods. This inconsistency might be caused by two different calculation methods. Moreover, the negative correlations between miRNAs and their corresponding targets (miR167-ARF8, miR393-bHLH77, miR5658-APL) could be found at the expression levels, suggesting the reliability of the sRNA-seq data and miRNA-mediated gene silencing involved in radish taproot thickening (Fig. 5).

Discussion
Integrative mRNA-seq and sRNA-seq approach provided valuable tool for exploiting the potential critical genes and uncovering complex regulatory networks for the traits [48][49][50]. Radish taproot thickening is a complex biological process that consisted of a series of material accumulation and signal transduction pathway. To date, the molecular mechanism of taproot formation was still not fully uncovered in radish. In this study, an integrated mRNA-seq and sRNA-seq analysis was performed during the taproot thickening in three radish advanced S2 VS S1-S3 VS S1-S3 VS S2 alpha-Linolenic acid metabolism ath00592 3 0.005396308 Starch and sucrose metabolism ath00500 5 0.006062768 inbred lines to further understanding the molecular mechanism underlying the taproot formation. As size varied in fleshy root, three representative radish advanced inbred lines were used in this study. A total of 2606 DEGs were shared between 'NAU-DY' and 'NAU-YB', whereas 16 DEMs were shared between 'NAU-YH' and 'NAU-DY' and 12 DEMs were specifically differentially expressed in 'NAU-DY'. Moreover, several critical genes including SUS1, EXPB3, and CDC5 were characterized and profiled by RT-qPCR analysis. This study represents a systematical report on characterization of the potential critical genes involved in taproot formation by integrative mRNA-seq and sRNA-seq approach in three radish advanced lines. Critical genes related to root -type differences Genes related to root-type difference are considered to be candidates that promote or repress taproot thickening in radish [51]. A total of 140 and 70 specifically expressed genes were identified in skinny-type roots and thick-type roots using suppression subtractive hybridization (SSH), respectively [18]. Several genes involved in phenylpropanoid metabolism were overexpressed in the skinny-root-type cultivar, and phenylpropanoid was the precursor in lignin synthesis, suggesting lignin biosynthetic pathway was involved in radish taproot thickening [52]. Furthermore, genes related to ethylene production and root hair elongation were also overexpressed in skinny-type roots, and the reason may be that ethylene and lateral root development inhibit cell expansion and elongation of main root [53]. Therefore, the DEGs between skinny-and thick-root cultivars were probably involved in root-type variation.

SUS1 and CDC5 might contribute to taproot thickening in radish
Carbohydrate metabolism was prominently activated in thickening roots, particularly in cell proliferating tissues, among which the expression level of SUS1 was associated with root thickening rates [15]. In this study, the GO term 'carbohydrate metabolism process' (GO: 0005975) was significantly enriched and shared by any comparison pair for taproot thickening of radish, irrespective of the size of the root, which was in agreement with a previous report in radish [13,15], suggesting that carbohydrate metabolism would probably be vital for taproot thickening.
Interestingly, among these genes involved in carbohydrate metabolism process, SUS1 gene showed differential expression patterns among three advanced inbred lines of radish. For large acicular radish, the expression level of SUS1 (Rsa1.0_00483.1_g00003.1) reached peak at cortex splitting stage and much higher than those in medium obovate and small circular radish (FD536105). Furthermore, SUS gene was of importance on the development of potato tuberization [54], tomato fruits setting [55], carrot root formation [56], maize grain formation [57], wheat grain formation [58], and rice grain formation [59]. In this study, RT-qPCR validation result indicated that SUS1 displayed high expression patterns at CSS, particularity in 'NAU-DY', indicating SUS1 might play a vital role on taproot thickening process, particularly for the thickening of large acicular radish rather than the small circular radish (Fig. 4).
CDC5 is a cell cycle regulator that encoding a MYBrelated protein. Recently, increasing reports of AtCDC5 plays a positive regulation for miRNAs accumulation, which control plant growth and development [60]. Meanwhile, CDC5 gene was crucial to cell cycle during the G2 period, and the phase transition of G2 to M (G2/ M) was affected in the AtCDC5-RNAi transformants [37], and AtCDC5-VIGS transformants died before bolting and accelerated cell death [38]. Interestingly, CDC5 was identified to be up-regulated during taproot thickening, whatever RNA-seq or iTRAQ-seq method, which would be play important role in cell division of taproot thickening in radish [13,35]. In this study, RT-qPCR validation result was approximately consistent with those of previous studies, and CDC5 showed high abundance expression at ES (Fig. 4). These results preliminarily suggested that CDC5 might play an important role on the growth and development of radish taproot thickening.
COBRA might be required for cell elongation of radish taproot Plant growth and development are promoted by a series of targeted cell division and cell expansion. The plant cell wall provides fundamental mechanical support for the plant body and determinant of cell size and shape. As the main component of cell wall, cellulose microfibril organization is one of determinant of cell expansion. In this study, the GO term 'cellulose microfibril organization' (GO: 0010215) was significantly enriched in S3 vs S1 pair comparison between 'NAU-DY' and 'NAU-YB', and totally two COBRA and three COBRA-like genes including COBL2, COBL5 and COBL8 were identified to be differentially expressed (Additional file 2). Previous studies showed that COBRA involved in the cellulose synthesis, controlling the content of cellulose in plant cell wall and the function of cell directional elongation [61][62][63][64]. Interestingly, in this study, all COBRA genes were up-regulated, whereas all COBRA-like genes were down-regulated in S3 vs S1 pair comparison between 'NAU-DY' and 'NAU-YB', suggesting that COBRA might play a critical role on cell elongation for taproot thickening in radish.

Conclusions
This is the first report on integrative analysis of transcriptome and miRNA in three radish advanced inbred lines during taproot thickening. A total of 2606 DEGs were shared between 'NAU-DY' and 'NAU-YB', which significantly enriched in 'phenylpropanoid biosynthesis', 'glucosinolate biosynthesis', and 'starch and sucrose metabolism' pathway. Meanwhile, a total of 16 DEMs were shared between 'NAU-DY' and 'NAU-YH', whereas 12 miRNAs showed specifically differential expression during taproot thickening in 'NAU-DY' with large acicular root when compared with 'NAU-YH' with small circular root. Association analysis between DEMs and DEGs indicated that miR393-bHLH77, miR167-ARF8, and miR5658-APL might be related to root-type variation in radish. Furthermore, RT-qPCR validation results indicated that the DEGs/ DEMs evaluated were highly in agreement with the RNA-Seq data. These finding would provide valuable information on comprehensive understanding of the molecular regulatory mechanism underlying taproot thickening in radish, and facilitate further genetic manipulation and quality improvement of root vegetable crops.
Additional file 1 Figure S1. The sequence length distribution of small RNAs during radish taproot thickening. Figure S2. The heatmap of DEGs and DEMs in radish taproot thickening. a. Heatmap of DEMs in 'NAU-DY'; b, c represents heatmap of DEGs in 'NAU-DY' and 'NAU-YB', respectively. Figure S3. Comparative analysis of gene expression profiles from RNAseq and RT-qPCR. 20 randomly selected genes by RNA-seq and RT-qPCR showing different expression patterns in three comparative groups (DS2 vs DS1, DS3 vs DS1, and DS3 vs DS2) during taproot thickening. Each data point represents the log 2 normalized expression level obtained from RNA-seq (x axis) and RT-qPCR (y axis) analyses. Figure S4. RT-qPCR validation of 14 DEMs during radish taproot thickening. The relative expression of DEMs between DS2 and DS1 libraries (a), DS3 and DS1 libraries (b) and DS3 and DS2 libraries (c) were analyzed by the 2 −ΔΔCT method. Table S1. Primer sequences for RT-qPCR assay. Table S2. Summary of small RNA sequencing data. Table S3. Summary of mRNA sequencing data. Table S4. Detailed information of DEMs during radish taproot thickening in 'NAU-DY'. Table S6. DEMs and their corresponding targets extracted from RNA-seq.