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-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-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.

Methods

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 (‘NAU-DY’, 20 days after sowing, DAS; ‘NAU-YB’, 20 DAS; ‘NAU-YH’, 10 DAS), CSS (‘NAU-DY’, 30 DAS; ‘NAU-YB’, 25 DAS; ‘NAU-YH’, 20 DAS) and ES (‘NAU-DY’, 55 DAS; ‘NAU-YB’, 45 DAS; ‘NAU-YH’, 40 DAS) from five randomly selected individual plants, respectively. Equal amount of samples from five individuals was pooled for library preparation and sequencing. All harvested taproot samples were immediately frozen in liquid nitrogen and stored at − 80 °C for RNA extraction.

Fig. 1

The characteristics of radish growth and development in three advanced lines. a. NAU-DY, b. NAU-YB, c. NAU-YH, d. The growth and development curve of radish. The asterisk represents the plant materials of RNA-Seq, and a bar represents 5 cm. 1 represents pre-cortex splitting stage (PSS), 2 represents cortex splitting stage (CSS) and 3–6 represents expanding stage (ES)

mRNA-seq and sRNA-seq library construction and sequencing

Prior to mRNA library construction, total RNAs were extracted at three different developmental stages of taproot thickening with the TRIzol reagent (Invitrogen) according to the instruction manuals, respectively. The detailed experimental procedures of cDNA library construction and sequencing were performed according to the reported approaches [20, 21]. The corresponding mRNA-seq libraries were named as S1 (‘NAU-DY’, DY_S1; ‘NAU-YB’, YB_S1), S2 (‘NAU-DY’, DY_S2; ‘NAU-YB’, YB_S2), and S3 (‘NAU-DY’, DY_S3; ‘NAU-YB’, YB_S3) library, respectively.

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 sRNA-seq analysis was shown in Fig. 2.

Fig. 2

The basic workflow of the integrative transcriptome and miRNA experiment in radish

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 |log2FC (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 |log2 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-△△ 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).

Results

Analysis of sRNA and mRNA sequencing data

A total of 14.81 M, 15.65 M, and 16.95 M clean reads were obtained from three development stages during taproot thickening in ‘NAU-DY’. For the length of 18 to 30 nt reads, a total of 10.86 M, 12.20 M, and 13.06 M clean reads were generated from DS1, DS2, and DS3 library, respectively. All the three libraries in length distribution showed a similar size characteristic ranging from 18 to 30 nt, especially for the majority of sRNA reads enriched at the length of 21 and 24 nt (Additional file 1: Figure S1). Among them, 7.56 M (69.59%) sRNA, 9.58 M (78.52%) sRNA, and 10.55 M (80.79%) sRNA reads were successfully mapped to reference sequences corresponding to DS1, DS2, and DS3 library, respectively (Additional file 1: Table S2).

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). In the line of ‘NAU-DY’, a total of 34.84 M, 25.42 M, and 28.28 M clean reads were obtained from DY_S1, DY_S2, and DY_S3, respectively, from which 22.62 M (64.91%), 16.19 M (63.69%), and 18.10 M (64.01%) were correspondingly uniquely mapped (Additional file 1: Table S3). While in the line of ‘NAU-YB’, totally 32.54 M, 33.10 M, and 31.62 M clean reads were generated from YB_S1, YB_S2 and YB_S3, respectively, among which 21.14 M (64.95%), 22.36 M (67.56%), and 21.21 M (67.08%) were uniquely mapped from YB_S1, YB_S2 and YB_S3, respectively (Additional file 1: 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 up-regulated miRNAs and 49 down-regulated miRNAs were identified at cortex splitting stage. Meanwhile, 85 DEMs including 26 up-regulated and 59 down-regulated miRNAs 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 miRNAs 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 miRNAs exhibited characteristics of dynamic change during radish taproot thickening.

Fig. 3

Identification of DEMs during radish taproot thickening. a. Volcano plot of miRNAs; b. Venny chart of miRNAs

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 down-regulated 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 root-type differences.

Table 1

Overview of DEMs shared between ‘NAU-YH’ and ‘NAU-DY’

DEMs	NAU-YH			NAU-DY
	log2 (S2/S1)	log2 (S3/S1)	log2 (S3/S2)	log2 (S2/S1)	log2 (S3/S1)	log2 (S3/S2)
miR156a-3p	9.52	N	−9.52	−5.3591	− 9.405	− 4.3581
miR156c-3p	−13.12	− 13.12	N	− 3.3905	− 5.0363	−1.9579
miR157a-3p	N	−1.72	− 2.44	− 2.2267	− 6.8916	− 4.977
miR171b-3p	10.26	N	−10.26	−2.627	− 3.4505	− 1.1357
miR171b-5p	−9.86	−9.86	N	−2.3005	−5.8566	−3.6757
miR172c	3.37	2.74	N	5.5808	8.3388	2.4459
miR172e-3p	−4.22	−4.28	N	6.0129	7.5173	1.1922
miR394a	7.73	7.25	N	1.6853	1.9103	N
miR395a	−7.57	−7.57	N	3.3711	3.9826	N
miR395b	7.41	6.36	−1.05	2.4611	1.8126	N
miR397a	−3.92	−6.29	−2.38	−5.7315	−4.3851	1.0342
miR408-5p	−2.01	− 4.93	− 2.91	−1.9648	−2.6612	− 1.0086
miR824-3p	−3.36	−3.53	N	−1.6303	N	1.3092
miR828	−1.80	−1.66	N	−6.2059	−5.7862	N
miR400	−1.19	−1.86	N	−1.2587	N	N
miR858a	N	1.38	1.19	−3.3814	−1.035	2.0342

Table 2

Identification of DEMs specific to ‘NAU-DY’ with large acicular root

sRNA	DS2 reads	DS1 reads	log2FC	DS3 reads	DS1 reads	Log2FC	DS3 reads	DS2 reads	log2FC
miR165a-3p	N	N	N	4219	12,833	−2	3692	11,787	−2
miR165a-5p	N	N	N	83	295	−2	73	315	−2
miR167a-3p	63	13	2	113	10	4	1646	760	1
miR167c-3p	39	183	−2	1881	546	2	N	N	N
miR167c-5p	42	332	−3	42	137	−2	N	N	N
miR167d	677	1452	−1	58	248	−2	N	N	N
miR319a	162,378	373,420	−1	N	N	N	475,335	131,886	2
miR5658	N	N	N	7	0	4	6	0	4
miR8175	N	N	N	407	141	2	356	102	2
miR857	3	199	−6	8	149	−4	7	2	2
miR170-5p	0	57	−7	0	43	−6	N	N	N
miR393a-5p	N	N	N	5	14	−1	N	N	N

DEGs identification and functional enrichment analysis

A total of 4131, 4979, and 1635 genes were differentially expressed in YB_S2 vs YB_S1, YB_S3 vs YB_S1, and YB_S3 vs YB_S2, respectively. Meanwhile, 2499, 3970, and 1924 DEGs were identified from DY_S2 vs DY_S1, DY_S3 vs DY_S1, and DY_S3 vs DY_S2, respectively. Of these DEGs, a total of 2606 DEGs were shared between ‘NAU-YB’ and ‘NAU-DY’ (Additional file 2: Table S5). Heatmaps of DEGs indicated that the expression patterns of DEGs were similar between ‘NAU-YB’ and ‘NAU-DY’ (Additional file 1: Figure S2); whereas more DEGs existed and their expression patterns were more various in ‘NAU-YH’ than those in ‘NAU-YB’ and ‘NAU-DY’ [13].

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.

Table 3

GO significantly enrichment of DEGs shared between ‘NAU-DY’ and ‘NAU-YB’

GO_accession	Description	Corrected_p-Value	DEG number
S2 VS S1
GO:0005975	carbohydrate metabolic process	0.0036482	32
S3 VS S1
GO:0010215	cellulose microfibril organization	0.024918	5
GO:0016049	cell growth	0.024918	5
GO:0016051	carbohydrate biosynthetic process	0.024918	29
GO:0006556	S-adenosylmethionine biosynthetic process	0.024918	4
S2 VS S1-S3 VS S1
GO:0005975	carbohydrate metabolic process	1.40E-06	72
GO:0055114	oxidation-reduction process	2.49E-06	90
GO:0044283	small molecule biosynthetic process	5.82E-06	34
S2 VS S1-S3 VS S2
GO:0005975	carbohydrate metabolic process	4.67E-09	22
GO:0008152	metabolic process	0.0048322	56
S2 VS S1-S3 VS S1-S3 VS S2
GO:0005975	carbohydrate metabolic process	0.010675	11

Table 4

KEGG significantly enrichment of DEGs shared between ‘NAU-DY’ and ‘NAU-YB’

#Term	ID	Input number	Corrected P-Value
S3 VS S1
Phenylalanine metabolism	ath00360	16	0.049828146
Cysteine and methionine metabolism	ath00270	14	0.051167179
Phenylpropanoid biosynthesis	ath00940	18	0.0527419
Ubiquinone and other terpenoid-quinone biosynthesis	ath00130	6	0.218324112
ABC transporters	ath02010	5	0.223834934
Biosynthesis of secondary metabolites	ath01110	64	0.414801288
Isoquinoline alkaloid biosynthesis	ath00950	4	0.414801288
Cutin, suberine and wax biosynthesis	ath00073	4	0.414801288
S3 VS S2
Thiamine metabolism	ath00730	3	0.020228102
S2 VS S1-S3 VS S1
Glucosinolate biosynthesis	ath00966	13	3.79E-07
Starch and sucrose metabolism	ath00500	27	0.001822699
2-Oxocarboxylic acid metabolism	ath01210	15	0.002400877
Biosynthesis of secondary metabolites	ath01110	87	0.002498601
Phenylpropanoid biosynthesis	ath00940	22	0.003925529
Phenylalanine metabolism	ath00360	18	0.004088511
Sulfur metabolism	ath00920	10	0.004445186
Selenocompound metabolism	ath00450	5	0.067434074
Tryptophan metabolism	ath00380	8	0.067434074
Cysteine and methionine metabolism	ath00270	13	0.067434074
Nitrogen metabolism	ath00910	7	0.126626602
Phenylalanine, tyrosine and tryptophan biosynthesis	ath00400	8	0.159801465
Glycine, serine and threonine metabolism	ath00260	9	0.159801465
Glutathione metabolism	ath00480	11	0.159801465
Tyrosine metabolism	ath00350	6	0.215571217
Flavonoid biosynthesis	ath00941	4	0.239223498
S2 VS S1-S3 VS S2
Fructose and mannose metabolism	ath00051	5	0.007626905
Amino sugar and nucleotide sugar metabolism	ath00520	7	0.007626905
Starch and sucrose metabolism	ath00500	8	0.013150643
S3 VS S2-S3 VS S1
Flavonoid biosynthesis	ath00941	4	0.035061213
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

DEM target prediction and association analysis of DEGs and DEMs

A total of 1280 targets were predicted with psRobot software from sRNA data of ‘NAU-DY’ (Additional file 1: Table S6). Among these targets, totally 849 and 431 targets belonged to 68 known miRNAs and 36 novel miRNAs, respectively (Additional file 1: Table S6). Of these, even though these miRNAs targeted different genes, the annotations of several targets belonged to the same families, for example, several known miRNAs including miR156a-3p, miR159b-3p, miR159c, miR5658, miR827, miR828, miR858a, miR858b, and novel miRNAs including novel_119, novel_136, novel_1, novel_37, and novel_96 targeted Myb domain genes (Additional file 1: Table S6).

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).

RT-qPCR validation

To explore genes expression patterns during radish taproot thickening, 20 genes were performed for RT-qPCR analysis from RNA-seq data of ‘NAU-DY’ and ‘NAU-YB’. RT-qPCR results showed that the expression levels of auxin-induced in root cultures protein 12 (RsAIR12), auxin transporter protein 1 (RsAUX1), auxin efflux carrier component 3 (RsPIN3), dihomomethionine N-hydroxylase (CYP79F1), glutamate synthase 1 (GLT1), glutathione S-transferase F3 (GSTF3), and acid beta-fructofuranosidase 4 (BFRUCT4) were higher at PSS than that at CSS and ES, whereas the expression level of auxin-responsive protein (IAA26) and sucrose synthase 1 (SUS1) were higher at CSS than that at PSS and ES (Fig. 4). Meanwhile, expansin B3 (EXPB3), cell division cycle 5-like protein (CDC5), sucrose-phosphate synthase 1 (SPS1), jacalin-related lectin 34 (JAL34), glutamine synthetase cytosolic isozyme 1–2 (JLN1–2), aspartic protease in guard cell 2 (ASPG2), aspartic proteinase A1 (APA1), beta-amylase 5 (BAM5), phosphoenolpyruvate carboxylase 2 (PPC2), and glucose-1-phosphate adenylyltransferase (AGPS1) were higher at ES than that at CSS and PSS (Fig. 4). In addition, RT-qPCR expression profiles showed that 18 of 20 randomly selected genes were in agreement with the RNA-Seq data, and two genes including JAL34 and AUR showed difference only at ES and CSS, respectively. Overall, the Pearson correlation coefficient results displayed a positive correlation between RNA-seq data and RT-qPCR analysis at the mRNA level (r = 0.78, p = 0.0006), indicating the reliability of the transcriptomic data (Additional file 1: Figure S3).

Fig. 4

The validation of expression levels of selected DEGs related to radish taproot thickening

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 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).

Fig. 5

A putative regulatory model of radish taproot thickening and development. a. Heatmap of shared DEGs among three advance inbred lines; b. The expression profile of DEMs and their corresponding target genes; c. A proposed regulatory model of radish taproot thickening and development

Molecular dynamic regulation network of taproot formation and development

Radish taproot thickening has a complex molecular regulation mechanism [13, 15, 17]. In general, the substances involved in morphogenesis including sucrose, starch, and protein are indispensable for taproot formation and development in radish, and a series of genes (e.g. SUS1, SUS3, SUC1, SUC2) were involved in these processes. Among them, SUS1 might contribute to taproot thickening in radish with long and thick root [15]. Genes involved in plant hormone biosynthesis and signal transduction including auxin (e.g. SAUR40, IAA26, IAA3, IAA7), CTK (e.g. LOG8, CKX6), and GA (e.g. GA2ox1, GA3ox2, RGA2) are likely to be important regulators in signal transduction process for radish taproot thickening [34, 35]. Cell cycle is regulated by F-box protein (e.g. CYCU4–1, CYCD3–1, CDC5) and cyclin-dependent kinase (e.g. KRP5) [36-38], which might determine cell number of radish taproot thickening. Cell expansion is regulated by expansin proteins (e.g. EXPA9, EXPB1/3), xyloglucan endotransglucosylase/ hydrolase protein (e.g. XTH9, XTH7), and root cell elongation factor (e.g. COBRA) [35, 39], which might determine cell size during taproot thickening. In addition, several transcription factors (TFs) were found to be involved in root development in other plants including rice, and Arabidopsis, and some TFs including bHLHs (e.g. miR393-bHLH77, bHLH96, bHLH148), MYBs (MYB4, MYB28, MYB44), and WRKYs (e.g. WRKY12, WRKY19, WRKY33) were identified in this study, indicating that they may play critical roles in taproot thickening in radish (Additional file 2: Table S5) [40-43].

Totally, 45 shared genes were differentially expressed in three advanced inbred lines with different root size, and their corresponding expression patterns were shown in Fig. 5a. The results indicated that most DEGs shared the similar expression patterns, whereas several DEGs (e.g. SUS1, SUC1, EXPB3, EXLA2, and ERF109) exhibited differential expression patterns. Interestingly, gibberellin 3 oxidase 2 (GA3ox2), cytokinin oxidase 6 (CKX6), miR167c-5p (targeted by ARF8), and miR393a-5p (targeted by bHLH77) were specifically differentially expressed in ‘NAU-DY’. Based on these results, a putative regulatory model of radish taproot thickening and development was put forward (Fig. 5b, c). In short, sucrose (the key gene of sucrose metabolism: SUS1) as a signal molecular could induce miR393-bHLH77 specific expression [44, 45], and then regulation of downstream genes (CYCs) played important roles in cell division during radish taproot thickening [40, 43, 46]. Some genes involved in the biosynthesis and signal transduction of auxin (miR167-ARF8), CTK (CKX6) and GA (GA3ox2) also regulated several functional genes (COBRA, EXPB1/3, CDC5), indicating that they might play roles in initiating radish taproot thickening and development (Fig. 5) [47].

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-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 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.

In this study, a total of 12 DEMs belonging to miR165, miR167, miR319, miR5658, miR8175, miR857, miR170, and miR393 family miRNAs were specifically differentially expressed during taproot thickening in large acicular radish ‘NAU-DY’, which might be related to root-type variation. On the other hand, many DEGs were identified to be specifically or commonly expressed between ‘NAU-YH’ and ‘NAU-DY’ (Fig. 5). For example, genes involved in starch and sucrose metabolism including sucrose synthase 4 (SUS4), sucrose-phosphate synthase 2 (SPS2), and 6-phosphofructokinase (FK1, FK7) were specifically expressed in ‘NAU-DY’, whereas SUS1, sucrose synthase 3 (SUS3), sucrose transport protein (SUC1, SUC2), sucrose-phosphate synthase 2 (SPS1), and 6-phosphofructokinase (FK3, FK6) were co-expressed in ‘NAU-DY’ and ‘NAU-YH’. Association analysis of mRNA and sRNA results showed that miR167c-5p (targeted by ARF8), miR393a-5p (targeted by bHLH77), miR5658 (targeted by APL) were identified to be significantly differentially expressed during radish taproot thickening of ‘NAU-DY’, indicating that miR393-bHLH77, miR167-ARF8, and miR5658-APL (a member of MYB) might play critical roles in taproot thickening of radish.

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 MYB-related 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-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 RNA-seq 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 log2 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−ΔΔ 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.

Additional file 2 Table S5. Overview of DEGs shared between ‘NAU-DY’ and ‘NAU-YB’.