Transcriptomic dynamics changes in development of carmine radish (Raphanus sativus L.) fleshy roots using RNA-seq method

Radish ( Raphanus sativus L.), belonging to biennial root vegetable crop of Brassicaceae family, is an economically important vegetable crop with an edible taproot. Recently, most of differential expressed genes associating with anthocyanin biosynthesis have been identified in most of important fruit crops. However, transcriptome analysis of anthocyanin biosynthesis and expression of anthocyanin biosynthesis related genes in ‘Hongxin’ radish have not been fully investigated. Here, based on results from HPLC analysis, young fleshy roots obtained from the dynamics development stage of fleshy roots in carmine radish ‘Hongxin 1’ was used for RNA-Seq, including fleshy roots from seedling stage (SS), initial expansion (IE), full-expansion (FE), bolting stage (BS), initial flowering stage (IFS); full-bloom stage (FBS) and podding stage (PS). Subsequently, the putative candidate genes involved in the dynamics development stage of fleshy roots in carmine radish were identified. After that, DGE (differential gene expression) profile analysis was used to identify the pupative transcripts, compared with fleshy roots from seedling stage (SS). In addition, co-modulated DEGs (Common DEGs in the dynamic growing stages of fleshyroot in carmine radish) were also identified, from which most DGEs were more likely to participate in anthocyanin biosynthesis, including two transcription factors RsMYB and Rs RZFP . In addition, some related proteins e.g. RsCHS , RsDFR , RsANS , RsF’3H , RsF3GGT1 , Rs3AT1 , glutathione S-transferase F12, RsUFGT78D2-like and RsUDGT-75C1-like were significantly contributed to the regulatory mechanism during anthocyanin synthesis in the development stage of fleshy roots. Furthermore, GO terms comprised of “anthocyanin-containing compound biosynthetic process” and “anthocyanin-containing compound metabolic process” were commonly overrepresented in the other dynamics growing stages of fleshy roots after initial expansion of fleshy and flavonol biosynthesis, Diterpenoid biosynthesis, Anthocyanin biosynthesis, as well as Benzoxazinoid biosynthesis. These results will expand our understanding of complex molecular mechanism of the putative candidate genes involved in the dynamics development stage of fleshyroot in carmine radish.

growing stages of fleshy roots in carmine radish, including Flavonoid biosynthesis, Flavone and flavonol biosynthesis, Diterpenoid biosynthesis, Anthocyanin biosynthesis, as well as Benzoxazinoid biosynthesis. These results will expand our understanding of complex molecular mechanism of the putative candidate genes involved in the dynamics development stage of fleshyroot in carmine radish.

Background
Anthocyanins, recognized as regulator for red to purple colors in nature, thereby producing water-soluble pigments belonging to the flavonoid group [1]. Most research has demonstrated that anthocyanins, as the benefit food additive worldwide, could release major public health threat (comprised of cardiovascular disease, inflammatory, obesity, and diabetes) caused by chemical-synthesis food additive [2,3]. In addition, most of regulatory genes have been extensively found involved in the anthocyanin biosynthetic pathway, which were largely conserved among flowering plants [4]. Previous studies have demonstrated that anthocyanins are firstly formed from phenylalanine by series enzymes through phenylpropanoid metabolism, such as phenylalanine ammonia-layse (PAL), cinnamic 4-hydroxylase (C4H) and 4-coumarate-CoA ligase (4CL); subsequently followed by chalcone synthase (CHS), and then the product 4, 2′ 4′ 6′-tetrahydrocychalcone are further catalyzed successively by four enzymes [Chalcone Isomerase (CHI), flavanone 3-Hydroxylase (F3H), dihydroflavonol 4-Reductase (DFR), as well as anthocyanidin synthase (ANS/LDOX)] [5,6]. However, the molecular mechanism of anthocyanin biosynthesis regulation is still not fully understood in the dynamic development of fleshy roots in radish.
Recently, based on the global transcriptome technology (such as RNA-seq technology), most of differential expressed genes associating with anthocyanin biosynthesis and expression of anthocyanin biosynthesis have been identified in most of important fruit 4 crops. However, transcriptome analysis of anthocyanin biosynthesis and expression of anthocyanin biosynthesis related genes in 'Hongxin' radish (That is famous for containing natural red pigment (red radish pigment), as produced in Fuling district of Chongqing city, China) have not been fully investigated.

Dynamics Anthocyanidin profiles of fleshy root in development of carmine radish
To demonstrate the dynamics Anthocyanidin profiles of fleshy root in development of carmine radish, Anthocyanidin profiles of fleshy root from seedling stage, initial expansion,full-expansion,Bolting stage, initial flowering stage,full-bloom stage and podding stage in five local cultivars of carmine radish ('Hongxin 1', 'Guanguan','Longquan 1',' Yanzhi 1' and 'Yanzhi 2') were investigated by HPLC analysis. The results showed that Anthocyanidin was significantly increased in 'Hongxin 1' differential fleshy roots types with the development of dynamic growing stages of fleshy roots than other local cultivars of carmine radish, from seedling stage to full-bloom stage, but decresed in podding stage ( Fig. S1).

Illumina sequencing and De novo assembly
Based on results from HPLC analysis, young fleshy roots obtained from the development stage of carmine radish 'Hongxin 1' was used for RNA-Seq in this study. The cDNAs obtained from fleshy roots of seven growth phases (seedling stage, initial expansion, fullexpansion, Bolting stage, initial flowering stage, full-bloom stage and podding stage) were sequenced using Illumina sequencing technology. After filtering out adaptor-only reads, trimming reads and low-quality reads (base quality ≤10), high-quality Reads were aligned to the SSU and LSU rRNA sequences to remove rRNA reads by a home-made perl script.
And the percentage of clean reads counts almost 70% after removing rRNA sequences on average among raw tags in each library, respectively (Table 1). Moreover,198,342 assembled transcipts from the raw sequence reads were constructed with an average length of 411 bp, and 34,927 Unigenes were generated using paired-end reads with an average length of 768 bp through de novo assembly technology (Fig. S2).

DEGs related to the dynamics growing stages of fleshy roots in carmine radish
To identify the putative candidate genes with great changes involved in the dynamics growing stages of fleshy roots in carmine radish, normalized expression levels for all global expressed genes were analyzed, indicating that high distinct gene expression profiles exists in the dynamics growing stages of fleshy roots (Fig. 1, Table S2). More interestingly, we found that the putative candidate genes belong to Cluster 8 was consistently with dynamics anthocyanidin profiles of fleshy root in development of carmine radish, but the putative candidate genes categorized into Cluster 9 were found oppositely. transcripts) ( Fig. 2A). 126 Co-modulated DEGs (Common DEGs in the dynamic growing stages of fleshy roots in carmine radish) were identified based on venny graph (Fig. 2B), and expression changes pattern of co-modulated DEGs were displayed with different colors using heatmap (Fig. 2C). More importantly, we found some of co-modulated DEGs  Table S3). Moreover, these DEGs related to the dynamics growing stages of fleshy roots that involved in different biological processes were validated using qRT-PCR and the results showed higher consistent with expression profiles of RNASeq data (Fig. 3, Table S4).

Functional annotation of DEGs related to the dynamics growing stages of fleshy roots in carmine radish
To explore the regulatory mechanisms of DEGs related to the dynamics growing stages of fleshy roots in carmine radish, GO annotation and KEGG pathway enrichment of those putative DEGs were conducted. The results illustrated that GO terms comprised of "anthocyanin-containing compound biosynthetic process" and "anthocyanin-containing compound metabolic process" were commonly overrepresented in the other dynamics growing stages of fleshy roots after initial expansion of fleshy roots (IE-40 days after planting), "flavoriod biosynthetic process" and "flavoriod metabolic process" were found overrepresented in fleshy roots of IFS, FBS and PS; but for "pigment biosynthetic process" and "pigment metabolic process", which was only found overrepresented in fleshy roots of IFS and FBS; moreover, we found that GO terms comprised of "glucosinolate biosynthetic process" and "glucosinolate metabolic process" were only overrepresented in fleshy roots of IE. (Fig. 4A, Table S5). By conducting pathway enrichment analysis, these results indicated that five significantly enrichment pathways DEG were identified for the dynamics growing stages of fleshy roots in carmine radish, including Flavonoid biosynthesis, Flavone and flavonol biosynthesis, Diterpenoid biosynthesis, Anthocyanin biosynthesis, as well as Benzoxazinoid biosynthesis (Fig. 4B, Table S6). (Cluster_24268), MYB transcription factor (Cluster_28373), as well as Zinc finger, RINGtype protein (Cluster_7186). To date, GSTs were found involved in anthocyanin transport based on genetic and biochemical evidence [17]. Bz2, firstly demonstrated in Zea mays by its mutant bronze-2 as a GST-encoding gene, which was found involved in vacuolar transfer of anthocyanins (bz2) [18]. Anthocyanin accumulation and pigment mislocalization were found reduced in Arabidopsis, caused by Mutations in the GST-encoding genes [19]. In this study, GSTs F12 was significantly up-regulated in differential fleshy roots types with 9 the development of dynamic growing stages of fleshy roots. These findings provided further evidences for the role of GSH in anthocyanin transport mechanisms. Moreover, MYB is a key component of the central regulatory to determine variation of anthocyanin production [20]. More importantly, MYB transcription factors have been extensively studied for their roles in the regulation of pigmentation in plants. It is known that R2R3type MYB proteins and the MYB-bHLH-WD40 complex directly activate the transcription of structural genes in the anthocyanin pathway, such as transcription of the early (CHS, CHI, F3'H and FLS) and late (DFR, ANS and ANR) flavonoid biosynthesis genes, respectively [21].

Discussion
In this study, we also demonstrated that MYB transcription factors was significantly dynamically up-regulated and showed remarkable positive and significant correlation to red pigment content in differential fleshyroots types (Table S3). So we inferred that MYB transcription factors might specifically activate early flavonoid biosynthesis genes comprise of CHS , CHI and F3'H, as well as late flavonoid biosynthesis genes consisted of DFR and ANS in carmine radish, thereby directly playing important roles in anthocyanin biosynthesis. However, their molecular regulation mechanism awaits further investigation. We used a 10 μL injection volume for a VDS C-18 column (4.6 × 250 mm, 5 μm, VDS Optilab, Germany) with 0.8 mL min-1 flow rate. Based on results from HPLC analysis, young fleshy roots obtained from the development stage of carmine radish 'Hongxin 1' was used for RNA-Seq in this study (Fig. S1).

Sample preparation and library construction
Library construction was conducted following NEBNext Ultra RNA Library Prep Kits for Illumina (NEB, USA), the mRNA was isolated through magnetic beads with Oligo (dT) using approximately 5 µg of total RNA and subsequently converted into short fragments by fragmentation buffer. After that, short fragments were converted into the first strands of cDNA used as templates with random hexamers, as well as the second strands of cDNA were also synthesized, and then the desired synthesized cDNA fragments were purified (QiaQuick PCR kit) for PCR amplification, and the quantify and qualify of each sample library were checked by agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System. Ultimately, 200bp paired-end reads were generated from the prepared library with 2 replicates using Illumina HiSeqTM 2000.

Reads processing and differentially expressed genes (DEGs) identification
Using Trimmomatic software, clean reads were obtained through filtering out adaptor-only reads, trimming reads and low-quality reads (base quality ≤10). After that, high-quality Reads were aligned to the SSU and LSU rRNA sequences download from sillva database using bwa with paramers "-n 4 -o 1 -e 1 -i 0 -l 50 -k 2" and the mapped reads were removed by a home-made Perl script as rRNA reads. Subsequently, we used Trinity software to de novo assembled clean reads into transcripts, and all their responding unigenes were annotated through searching in national center for biotechnology information (NCBI) non-redundant protein (Nr) databases and Swiss-Prot protein databases using BLASTx search tool with threshold E-value set as less than 10. After that, the functions of assembled unigenes were annotated through gene ontology (GO, http://www.geneontology.org/) database and Kyoto encyclopedia of genes and genomes (KEGG, http://www.kegg.jp/) database. To assess the abundances of assembled transcript, we firstly mapped the clean reads of seven different fleshy roots libraries to the de novo assembled transcriptome using Bowtie2, and then assessed with RSEM through transcript quantification of the de novo assembly, only transcripts (FPKM ≥ 1) were considered as significant expressed transcripts [7]. At last, DEGs (Differential expressed genes) were then screened by noiseqbio [8] and then identify using a corrected P-value <0.05 between each set of compared samples. (The fold change of gene expression of six cultivars of radish comprised of 'IE_root', 'FE_root', 'BS_root', 'IFS_root', 'FBS_root' and 'PS_root ' were identified by comparing with 'SS_root' respectively, including 'IE_root' Vs 'SS_root', 'FE_root' Vs 'SS_root', 'BS_root' Vs 'SS_root', 'IFS_root' Vs 'SS_root', 'FBS_root' Vs 'SS_root' and 'PS_root' Vs 'SS_root'). Furthermore, DEGs in the dynamic growing stages of carmine radish were analyzed and plotted using Neighbor-Joining cluster through homemade R script.

GO functional annotation and KEGG pathway analysis of co-modulated differently expressed genes (DEGs) in growing stages of carmine radish
Co-modulated DEGs (Common DEGs in the dynamic growing stages of fleshy roots in carmine radish) were identified based on venny graph. Subsequently, we conducted those co-modulated DEGs for GO annotation through Gene Ontology Database (http://www.geneontology.org/) and KEGG pathway enrichment analysis using KOBAS software, respectively [9]. In addition, the degree of KEGG enrichment was evaluated as the rich factor, q-value, and the number of genes in the enriched pathway. The rich factor refers to the ratio of the number of DEGs to the number of total annotated genes in a certain pathway. The q-value is a multiple hypothesis-corrected P value. The q-value can take on values between 0 and 1; values closer to 0 indicate more significant enrichment.
After that, R script was used to construct their relative graphs.

Validated of candidate DEGs invovled in the growing stages of carmine radish using real-time qRT-PCR
To confirm the results obtained from the RNA-Seq assay, 11 DEGs with great alteration that related to growing stages of carmine radish were chosen and validated by qRT-PCR.
The primers are designed by Primer 5.0 software for qRT-PCR experiments and radish gene (Actin) is used as a standard control ( Table S1). The amplification programs were performed according to the standard protocol of the ABI7500 system, and conducted in Table   Table 1 Table S1: List of primers for qRT-PCR analysis of 11 candidate DEGs involved in the dynamics growing stages of fleshy roots in carmine radish        Validate of candidate Co-modulated DEGs involved in the dynamics growing stages of fleshy roots in carmine radish using qRT-PCR and then correlation between RNA-seq and qPCR data were conducted. Each RNA-seq expression data was plotted against that from quantitative real-time PCR and fit into a linear regression. Both x-and y-axes were shown in log2 scale and each color represented a different gene.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.