Transcriptomic Analysis Identifies Differentially Expressed Genes (DEGs) Associated with Bolting and Flowering in Radish (Raphanus sativus L.)

The transition of vegetative growth to bolting and flowering is an important process in the life cycle of plants, which is determined by numerous genes forming an intricate network of bolting and flowering. However, no comprehensive identification and profiling of bolting and flowering-related genes have been carried out in radish. In this study, RNA-Seq technology was applied to analyze the differential gene expressions during the transition from vegetative stage to reproductive stage in radish. A total of 5922 differentially expressed genes (DEGs) including 779 up-regulated and 5143 down-regulated genes were isolated. Functional enrichment analysis suggested that some DEGs were involved in hormone signaling pathways and the transcriptional regulation of bolting and flowering. KEGG-based analysis identified 37 DEGs being involved in phytohormone signaling pathways. Moreover, 95 DEGs related to bolting and flowering were identified and integrated into various flowering pathways. Several critical genes including FT, CO, SOC1, FLC, and LFY were characterized and profiled by RT-qPCR analysis. Correlation analysis indicated that 24 miRNA-DEG pairs were involved in radish bolting and flowering. Finally, a miRNA-DEG-based schematic model of bolting and flowering regulatory network was proposed in radish. These outcomes provided significant insights into genetic control of radish bolting and flowering, and would facilitate unraveling molecular regulatory mechanism underlying bolting and flowering in root vegetable crops.

Transcriptomic Analysis Identifies Differentially Expressed Genes (DEGs) Associated with Bolting and Flowering in Radish (Raphanus sativus L.)

INTRODUCTION
The developmental transition from vegetative growth to bolting and flowering is one of the most important traits in plant life cycle. Bolting and flowering time must be appropriately determined to ensure reproductive success under most favorable conditions (Amasino and Michaels, 2010;Srikanth and Schmid, 2011). Plants have evolved an intricate bolting and flowering genetic circuitry in response to various endogenous and environmental signals including development, age, plant hormones, photoperiod, and temperature (Fornara et al., 2010;Capovilla et al., 2015;Kazan and Lyons, 2015). Molecular and genetic regulation of flowering has been extensively studied in the model plant Arabidopsis thaliana. Five major flowering pathways including vernalization, photoperiod, autonomous, aging and gibberellin (GA) pathways have been identified to govern bolting and flowering time (Amasino and Michaels, 2010;Fornara et al., 2010;Srikanth and Schmid, 2011), and a number of flowering-related genes involved in these pathways have been isolated and characterized in Arabidopsis (Fornara et al., 2010;Srikanth and Schmid, 2011).
The signals from flowering pathways converge on several floral pathway integrators such as FLOWERING LOCUS T (FT), SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1) and LEAFY (LFY), which are integrated into the genetic networks of flowering (Moon et al., 2005;Parcy, 2005). Among these integrators, the florigen gene FT is a central node of floral transition, whose transcriptional expression is positively regulated by CONSTANS (CO) encoding a putative zinc finger transcription factor (Suárez-López et al., 2001), while it is negatively regulated by FLOWERING LOCUS C (FLC), a flowering repressor encoding a MADS-box transcription factor (Michaels and Amasino, 1999). Different environmental factors affect plant flowering by modulating the expression of floral integrators and stimulating changes in plant hormone levels (Yaish et al., 2011;Riboni et al., 2014;Kazan and Lyons, 2015). Increasing evidences have revealed the connections between flowering time and plant hormones including salicylic acid (SA), jasmonic acid (JA), GA, abscisic acid (ABA) and auxin (Davis, 2009;Kazan and Lyons, 2015). The effects of phytohormone signaling on flowering, particularly GA pathway, have been extensively described in Arabidopsis (Mutasa- Göttgens and Hedden, 2009). Therefore, understanding the roles of floweringrelated genes and crosstalk between diverse genetic pathways is fundamental for elucidating the regulatory mechanisms underlying bolting and flowering in plants.
RNA sequencing (RNA-Seq), a powerful strategy for global discovery of functional genes, has provided a better qualitative and quantitative description of gene expressions under certain conditions in many plant species (Lister et al., 2009;Wang et al., 2009a). Digital gene expression (DGE) tag profiling is a revolutionary approach for identifying differentially expressed genes (DEGs) in diverse plant tissues, organs and developmental stages (Bai et al., 2013;Zhang et al., 2014a;Zhu et al., 2015). Moreover, RNA-Seq combined with DGE profiling has been employed for flowering-related gene discovery and expression analysis in some species such as bamboo (Gao et al., 2014), Lagerstroemia indica (Zhang et al., 2014b), sweetpotato (Tao et al., 2013) and litchi (Zhang et al., 2014c). However, to our knowledge, there are no studies on global expression profile analysis of bolting and flowering-related genes in radish (Raphanus sativus L.).
Radish (2n = 2x = 18), belonging to Brassicaceae family, is an important annual or biennial root vegetable crop worldwide. Premature bolting is a seriously destructive problem and results in poor root growth and the reduced harvest during radish production, especially in spring. Appropriate timing of bolting and flowering is significant for reproductive success at suitable conditions, as well as preventing the premature bolting in radish. Progress on bolting and flowering time control (Fornara et al., 2010;Srikanth and Schmid, 2011), especially in Arabidopsis, has provided a solid foundation and reference for identifying numerous functional genes during radish bolting and flowering. Recently, the transcriptomes from radish roots and leaves have been assembled and analyzed Zhang et al., 2013;Xu et al., 2015). Moreover, a list of microRNAs (miRNAs) and functional genes related to bolting and flowering were successfully isolated from late-bolting radish based on transcriptomic datasets (Nie et al., 2015). Therefore, to further identify the DEGs involved in bolting and flowering regulation is of importance for understanding the genetic regulatory network of bolting and flowering in radish.
In this study, to investigate the gene expression patterns during the transition of vegetative growth to bolting stage in radish, using the late bolting radish advanced inbred line as material, two DGE libraries were constructed and sequenced with RNA-Seq technology. The aims were to comprehensively identify DEGs involved in bolting and flowering regulatory network and to explore their roles in determining radish bolting and flowering time. Expression patterns of several critical DEGs related to bolting and flowering were validated by quantitative real-time PCR (RT-qPCR) analysis. Finally, to characterize the bolting and flowering-related genes and miRNAs in flowering pathways, a putative miRNA-DEG-based model of bolting and flowering regulatory network was put forward in radish. These results could provide significant insights into the molecular mechanism underlying bolting and flowering regulation in radish and other root vegetable crops.

Plant Materials
The late bolting radish advanced inbred line 'NAU-LU127, ' which was self-pollinated for more than 20 generations, was used in this study. The genetic background and structure of this line are stable and highly homozygous. After surface-sterilization, the seeds were sowed and grew in a growth chamber with 16 h light at 25 • C and 8 h darkness at 16 • C. The radish leaves used for DGE sequencing and RT-qPCR analysis were separately collected at two different developmental stages: vegetative stage (VS) and reproductive stage (RS), with three biological replicates. Each sample was collected at two developmental stages from three randomly selected individual plants, respectively. All the samples were immediately frozen in liquid nitrogen and stored at −80 • C until use.

DGE Library Construction and Illumina Sequencing
Total RNA from radish leaves at vegetative stage and reproductive stage was individually extracted using Trizol R Reagent (Invitrogen) following the manufacturer's instructions. The equivalent quantity of total RNA from three replicates was pooled and used for library preparation and sequencing. Two cDNA libraries named NAU-VS and NAU-RS were constructed and sequenced according to the previously reported method . The library construction and Illumina sequencing were performed using HiSeq TM 2500 platform at Beijing Genomics Institute (BGI, Shenzhen, China). The RNA-Seq data were deposited in NCBI Sequence Read Archive (SRA, http:// www.ncbi.nlm.nih.gov/Traces/sra/) with accession numbers of SRX1671036 (NAU-VS) and SRX1671054 (NAU-RS).

Data Processing and Expression Analysis of DEGs
The raw reads were primarily produced for data processing. After filtering low quality reads, adaptor sequences and reads containing ploy-N, the clean reads were obtained. These clean reads were then matched to the radish reference sequences which contained the public radish genomic survey sequences (GSS) and expressed sequence tag (EST) sequences and leaf transcriptome sequences from 'NAU-LU127' (Nie et al., in press) with no more than two mismatches. These sequences from radish leaf transcriptome have been deposited in NCBI Transcriptome Shotgun Assembly (TSA, http://www.ncbi.nlm.nih.gov/genbank/ tsa/) database under the accession number GEMG00000000.
To screen the DEGs between two DGE libraries, the expression level of each transcript is calculated using RPKM (Reads Per kb per Million reads) method (Mortazavi et al., 2008). Prior to differential gene expression analysis, the read counts of each transcript were adjusted by edgeR program package  through one scaling normalized factor. Trimmed Mean of M values (TMM), an appropriate normalization method implemented in the edgeR package , was employed to obtain the normalized read counts. The differential expression analysis of two libraries was performed using the DEGSeq R package 1.20.0 (Wang et al., 2010). Subsequently, the false discovery rate (FDR) was used to determine P-value threshold in multiple testing (Benjamini et al., 2001). A strict algorithm was used to further perform DEG identification according to the previous reports (Audic and Claverie, 1997). The absolute value of log 2 Ratio (NAU-RS/NAU-VS) ≥ 1, P < 0.05 and FDR ≤ 0.001 were used as threshold for judging the significance of gene expression difference. The cluster analysis of gene expression patterns was performed with cluster software and Java Treeview software (Saldanha, 2004).

Functional Annotation and Enrichment Analysis of DEGs
To investigate the biological function and involvement in functional pathways, all the identified transcripts were mapped to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. For GO annotation, the unique transcripts were subjected to BLASTX searching against the NCBI Nr database using the E < 10 −5 . Then the Blast2GO (Conesa et al., 2005) and WEGO software (Ye et al., 2006) were used to obtain GO annotations and functional classifications. GO enrichment analysis of DEGs was implemented by the GOseq R package . KOBAS software (Xie et al., 2011) was used to test the statistical enrichment of DEGs in KEGG pathways. The significantly enriched functional terms and pathways were identified using the criterion of a Bonferronicorrected P ≤ 0.05.

RT-qPCR Analysis
Total RNAs from radish leaves were isolated and obtained as described above. RT-qPCR was performed on an iCycler Real-Time PCR Detection System (Bio-Rad, USA) with three replications according to previous reports (Nie et al., 2015;Xu et al., 2015). All PCR reactions were carried out in a total volume of 20 µL with RsActin gene as the internal control (Xu et al., 2012). The relative gene expression levels were calculated using 2 − C T method (Livak and Schmittgen, 2001). The specific PCR primers were designed using Beacon Designer 7.0 (Premier Biosoft International, USA) and listed in Table S7.

DGE Library Sequencing and Data Analysis
To obtain global unique sequences from radish leaves, de novo assembly and analysis of transcriptome prepared from radish leaves of 'NAU-LU127' were carried out using Illumina RNA-Seq technology. Totally 111,167 contigs and 53,642 unigenes were generated from the radish leaf transcriptome (Nie et al., in press). The available dataset of radish leaf transcriptome integrating with the available radish GSS and EST sequences released in NCBI database enriched the radish reference sequences for DEG identification during radish bolting and flowering.
In this study, two DGE libraries from leaves of radish advanced inbred line 'NAU-LU127' at vegetative and reproductive stages were constructed and sequenced by Illumina HiSeq TM 2500 platform, respectively. As a result, 8,825,790 and 12,382,793 raw reads were obtained in NAU-VS and NAU-RS libraries, respectively (Table 1; Figure 1A). After removing adaptor sequences and low quality reads, 8,541,912 and 10,154,256 clean reads were generated in the two libraries ( Table 1). These clean reads were then mapped to the radish reference sequences, resulting in the generation of 87.83% (7,502,455 reads) and 68.50% (6,955,386 reads) matched reads in NAU-VS and NAU-RS libraries, respectively. For the variation of clean reads mapping percentage between two libraries, it may arise from the sample differences and the specific pre-processing of obtained reads . The more mapped reads from NAU-VS library implied that some stage-specific genes may be expressed only at vegetative stage and differentially expressed between these two DGE libraries. Further analysis revealed that 4,124,632 reads (48.29 %) in NAU-VS library and 5,200,880 reads (51.22 %) in NAU-RS library were uniquely matched ( Table 1).

Identification and Functional Enrichment Analysis of DEGs
The transcript abundance of each gene from two DGE libraries was calculated and analyzed by RPKM method. The threshold of |log 2 Ratio| ≥ 1 and FDR ≤ 0.001 were further used to determine the significantly DEGs. A total of 5922 significantly DEGs  including 779 up-regulated and 5143 down-regulated genes were obtained from NAU-VS and NAU-RS libraries (Table S1; Figure 1B).
To better classify the functions of these identified DEGs, GO classification and enrichment analysis were carried out in this study. These DEGs were categorized into three main GO categories including 23 biological processes, 14 cellular components and 13 molecular functions (Figure 2). Functional enrichment analysis revealed that 140 GO terms were significantly enriched with a Bonferroni-corrected P ≤ 0.05 ( Table S2). The terms of "metabolic process" (GO: 0008152) and "organic substance metabolic process" (GO: 0071704) were the dominant groups in biological processes; "cell" (GO: 0005623) and "cell part" (GO: 0044464) were the highly represented groups in the cellular components. For the molecular functions, a large proportion of genes were significantly enriched in "organic cyclic compound binding" (GO: 0097159) and "heterocyclic compound binding" (GO: 1901363) categories. Moreover, some enriched GO terms were related to plant flowering and meristem development, including "regulation of photoperiodism, flowering" (GO: 2000028), "regulation of timing of meristematic phase transition" (GO: 0048506), "meristem maintenance" (GO: 0010073), "meristem growth" (GO: 0035266), "meristem development" (GO: 0048507) and "flower development" (GO: 0009908) ( Table S2).

DEGs Involved in the Transition of Vegetative Growth to Bolting in Radish
In this study, to identify DEGs during radish bolting and flowering, BLAST searching was performed and the putative functions of DEGs were assessed. A total of 95 DEGs representing 128 unique sequences related to bolting and flowering were identified ( Table S5). The analysis of flowering pathways revealed that these genes were involved in five different flowering pathways including photoperiod, vernalization, autonomous, GA and aging pathways.

Expression Profile Analysis by RT-qPCR
To validate the differential expression patterns of DEGs during radish bolting and flowering, totally 21 functional genes were randomly selected and subjected to RT-qPCR analysis. These selected genes included six genes related to hormone signaling and 15 genes related to bolting and flowering regulation. The relative expression levels of these genes between vegetative growth and reproductive stage were analyzed and compared (Figure 4). Further, comparative analysis revealed that these gene expression trends except MYC2-CL4584.Contig1 were in agreement with the transcript abundance changes by RNA-Seq (Figure 4), indicating the highly accuracy and quality of DGE sequencing.

The Regulatory Network Underlying Bolting and Flowering in Radish
Considerable studies have revealed that some miRNAs regulating corresponding target genes played important roles in the transition from vegetative growth to bolting and flowering (Spanudakis and Jackson, 2014). In our recent study, several bolting and flowering-related miRNA-target gene pairs were identified and characterized in late-bolting radish (Nie et al., 2015). To better understand the genetic regulatory network of radish bolting and flowering, correlation analysis between the DEGs identified in the present study and bolting and floweringrelated miRNAs previously reported (Nie et al., 2015) was performed. As expected, 24 miRNA-mRNA pairs including 16 miRNAs and 27 target DEGs were identified (Table S6). Among them, 19 miRNA-mRNA pairs showed negative correlations in expression patterns. Several DEGs including AP2 (targeted by miR172), VRN1 (targeted by miR5227), PRP39 (targeted by miR6273), and NF-YB3 (targeted by miR860), were found to be involved in bolting and flowering regulation (Wang et al., 2007;Kumimoto et al., 2008;Zhu and Helliwell, 2010). To gain insights into the bolting and flowering regulatory network in radish, a putative model for summarizing the bolting and flowering-related DEGs and miRNAs was proposed (Figure 5). The critical genes involved in various flowering pathways and phytohormone signaling pathways were displayed in the schematic regulatory network of radish bolting and flowering. According to the known Arabidopsis flowering regulatory network (Fornara et al., 2010;Srikanth and Schmid, 2011), we speculated that the transcriptional regulations of several floral integrators including FT, CO, SOC1, FLC, and LFY, could integrate the signals from various pathways and modulate the radish bolting and flowering ( Figure 5). Moreover, the models of miR172-AP2 and miR5227-VRN1 have been shown to be important participants in the regulatory network of bolting and flowering (Wang et al., 2007;Zhu and Helliwell, 2010;Nie et al., 2015).

DISCUSSION
Radish bolting and flowering are integral stages in its complete life cycle. The timing of bolting and flowering is coordinately regulated by various endogenous and environmental signals integrating into a complexity of flowering regulation (Amasino and Michaels, 2010;Srikanth and Schmid, 2011). Recent advances in flowering genes and regulatory networks have greatly enhanced our knowledge of molecular basis underlying bolting and flowering-time control in Brassicaceae crops. However, no studies on comprehensive identification of DEGs related to radish bolting and flowering have been reported, and the regulatory mechanism of bolting and floweringtime control remains largely unexplored in radish. In this study, two cDNA libraries from leaves of radish advanced inbred line 'NAU-LU127' at vegetative and reproductive stages were constructed, respectively. A list of DEGs related to phytohormone signaling and transition from vegetative growth to bolting and flowering were identified and comprehensively profiled.

The Roles of Plant Hormone Signaling in Bolting and Flowering
Plant hormones are endogenously occurring compounds that regulate multiple aspects of plant growth and development including flowering time (Davis, 2009;Santner and Estelle, 2009). Various phytohormones have been implicated in the developmental transition of flowering (Davis, 2009;Domagalska et al., 2010). The pathways of several hormones including auxin, GA, ABA, SA, and JA signaling were significantly enriched by pathway-based analysis in our study (Table 3). GA pathway is one of the genetic flowering pathways, which could interact with several pathways and is integrated into the flowering regulatory complexity (Srikanth and Schmid, 2011). The role of GA pathway in flowering time has been thoroughly investigated in Arabidopsis and several fruit trees (Wilkie et al., 2008;Mutasa-Göttgens and Hedden, 2009). Many genes related to GA metabolism and signaling were involved in GA-mediated regulatory process of flowering (Mutasa-Göttgens and Hedden, 2009;Domagalska et al., 2010). GA exerting its biological functions on floral transition and development is mainly dependent on the growth inhibitor DELLA proteins (Mutasa-Göttgens and Hedden, 2009). GA signaling promotes flowering through initiating the degradation of transcriptional regulator DELLA and activating the expression of SOC1, AGL24 and LFY (Davis, 2009;Mutasa-Göttgens and Hedden, 2009). As expected, both the decreased transcript abundance and expression level of DELLA (CL7599.Contig1) were detected in radish reproductive stage compared with vegetative phase (Figure 4). The ABA pathway, which is antagonistic to GA, has been demonstrated to delay flowering through modulating DELLA activity and affecting the transcriptional expression of floral repressor FLC (Achard et al., 2006;Domagalska et al., 2010). In the current study, unique transcripts annotated as PP2C and ABF, ABA signaling components, were identified and differentially regulated during radish bolting and flowering (Figure 3), which is consistent with the results in litchi (Zhang et al., 2014c) and soybean (Wong et al., 2013). These findings suggested that the differential expressions of ABA signalingrelated genes may be associated with the timing of radish transition to bolting and flowering.
Function of SA in accelerating transition to flowering is pronounced by SA-deficient mutants of Arabidopsis (Martínez et al., 2004). SA could negatively regulate the floral repressor FLC and activate the flowering promoter FT which strongly highlights the positive role of SA in flowering transition (Martínez et al., 2004). SA promotes the activation of NON-EXPRESSOR OF PR-1 (NPR1) proteins, whose interaction with TGA transcription factors could induce the expression of PR genes (Wu et al., 2012). Moreover, JA is also implicated in flowering regulatory process and delays flowering in Arabidopsis (Krajnčič et al., 2006;Riboni et al., 2014). JA signaling pathway has been involved in three molecular elements including JA receptor gene COI1, transcriptional repressor JAZ protein and some transcription factors, e.g., the bHLH family (Krajnčič et al., 2006). Notably, recent studies have demonstrated the regulatory role of COI1 in delaying flowering mediating the repressed expression of FT (Zhai et al., 2015). In this study, some transcripts belonging to the main components of SA and JA signaling were found, including NPR1, TGA, PR, JAZ, COI1, and MYC2 (Table 3). In addition, previous studies reveal that auxin is necessary for flower initiation and floral organ identity (Cheng and Zhao, 2007). We also found the critical genes related to auxin signaling such as AUX1, SAUR, TIR1, and ARFs (Table 3; Table S4). Overall, these results reveal that phytohormone-mediated transcriptional reprogramming are crucial to the transition of bolting and flowering and participate in its regulatory network of radish. The characterization of critical genes in plant hormone signaling pathways would greatly help to illuminate the complex genetic network of bolting and flowering in radish.

The Complex Bolting and Flowering Regulatory Network in Radish
Multiple genetic flowering pathways integrating endogenous and environmental signals determine the transition from vegetative growth to reproductive development. Studies in Arabidopsis have revealed the participation of more than 200 flowering-related genes in the intricate regulatory network (Fornara et al., 2010;Srikanth and Schmid, 2011). In this study, 95 candidate genes related to bolting and flowering were isolated and involved in five major flowering pathways within genetic regulatory network (Table S5; Figure 5). It is inferred that known genetic pathways and critical flowering genes may conservatively present in radish, being consistent with the reports in maize (Dong et al., 2012), soybean (Jung et al., 2012), and citrus (Zhang et al., 2011). Gene expression profiling revealed that these genes were differentially expressed between NAU-VS and NAU-RS libraries, suggesting their putative important roles in radish bolting and flowering.
The complex regulatory network of Arabidopsis is composed of five major converging pathways (Fornara et al., 2010;Srikanth and Schmid, 2011). It is believed that endogenous developmental signals such as developmental stages of plants and phytohormones monitor flowering time through age, autonomous and GA pathways, while environmental cues regulate flowering time through the photoperiod and vernalization pathways in response to day length or temperature (Srikanth and Schmid, 2011;Capovilla et al., 2015). The signals from photoperiodic process are converted into the transcriptional regulation of key genes such as FT, CO, AP1, and AP2 to affect flowering time (Kikuchi and Handa, 2009;Amasino, 2010;Srikanth and Schmid, 2011). The florigen gene FT as a floral integrator is central for the photoperiodic flowering pathway of long-day plant Arabidopsis, which is perceived in leaves and transported to the shoot apex initiating floral transition (Huang et al., 2005;Parcy, 2005). The role of FT in promoting flowering has been proven by mutants and overexpressed transgenic analysis in Arabidopsis (Amasino, 2010;Srikanth and Schmid, 2011). As expected, the homolog of FT (FD571044) was up-regulated in reproductive stage of radish (Figure 4), indicating that the RsFT gene could positively regulate the development transition of bolting and flowering ( Figure 5). Under long-day condition, the FT expression is activated by CO, which is a floral activator and modulated by the circadian clock and day length (Suárez-López et al., 2001;Amasino, 2010;Johansson and Staiger, 2015). The link between circadian clock and flowering control may be mainly mediated by the transcriptional expression of CO (Fujiwara et al., 2008;Johansson and Staiger, 2015). In Arabidopsis, two essential circadian clock components LATE ELONGATED HYPOCOTYL (LHY) and CIRCADIAN CLOCK ASSOCIATED1 (CCA1) function in photoperiodic flowering and regulate flowering pathway by controlling the rhythmic expression of CO and FT (Fujiwara et al., 2008). In this study, some transcripts belonging to CO, CCA1 and LHY homologs were found to be up-expressed in reproductive phase with DGE sequencing and RT-qPCR analysis (Figure 4), suggesting the critical roles of these genes in the transition of radish bolting and flowering.
The vernalization and autonomous pathways converge on the flowering repressor FLC, and many genes involved in these two pathways could control flowering time through affecting FLC expression (Amasino, 2010). The high level of FLC delays flowering and requires its activator FRI (Michaels and Amasino, 1999;Choi et al., 2011). Recently, several naturally occurring spliced transcripts of FLC were found and isolated from B. rapa (Yuan et al., 2009) and orange , which were proven to be associated with variations in flowering time. The transcriptional co-expression analysis in B. rapa indicated that BrFLC2 may be the major regulator of flowering time in genetic flowering network (Xiao et al., 2013). In this study, we found putative homologs of FLC (CL1584.Contig3) from late-bolting radish, which was down-regulated in reproductive stage compared with vegetative stage, with similar patterns being detected in FRI and FRL (Figure 4). In addition, similar results were found in other homologous genes in vernalization pathway, including FIP2, EMF2, VRN1, and VRN2 (Table S5). These results indicate that the genetic elements of the vernalization pathway may be of importance for the manipulation of radish bolting and flowering time.
Furthermore, miRNAs as central regulators of gene expression have been shown to be implicated in multiple genetic pathways governing flowering time (Spanudakis and Jackson, 2014;Wang, 2014). The newly defined age pathway of flowering, which is controlled by miR156 and its target SPL transcription factors (Wang et al., 2009b), regulates flowering time and interacts with vernalization, photoperiodic and GA pathways (Zhou et al., 2013;Spanudakis and Jackson, 2014;Wang, 2014). Several members of SPL family were identified in this study, including SPL1, SPL2, SPL3, SPL9, SPL13, and SPL15 (Table S5). It was known that miR172 is down-regulated by the age-dependent expression of SPL9 Spanudakis and Jackson, 2014). The target genes of miR172 are a class of AP2-like transcription factors including AP2, TARGET OF EAT 1-3 (TOE1-3), SCHLAFMÜTZE (SMZ), and SCHNARCHZAPFEN (SNZ), which act as floral repressors (Zhu and Helliwell, 2010). The levels of these AP2-like genes are relatively high during plant seedling stage and decline with plant development, ultimately relieving the repression of flowering to trigger flowering (Aukerman and Sakai, 2003;Zhu and Helliwell, 2010). Consistent with these evidences, the down-expressed patterns of AP2 (CL1275.Contig1), SMZ (Rsa#S42015352), and RAP2.7 (CL2600.Contig3) were detected at reproductive stage in this study (Figure 4). Moreover, correlation analysis revealed that some bolting and flowering-related DEGs were targeted by specific miRNAs forming the transcriptional model of miRNA-mRNA pairs (Table S6). These findings reveal that some miRNA-DEG models including miR5227-VRN1, miR6273-PRP39, and miR860-NF-YB3 are crucial participators and integrated into the intricate genetic networks of bolting and flowering in radish ( Figure 5).

CONCLUSIONS
In summary, RNA-Seq technology was employed to systematically identify DEGs at transcriptome-wide level during radish transition from vegetative growth to bolting and flowering in this study. To our knowledge, this is the first investigation to illustrate the expression profiles of bolting-related genes and dissect the bolting and flowering regulatory network in radish. In this study, a total of 5922 DEGs were identified from late-bolting radish leaves. Several candidate genes related to plant hormone signal and bolting and flowering regulatory pathways were characterized and implicated in the complex networks of bolting and flowering regulation. Correlation analysis suggested that the miRNA-mRNA regulatory models played pivotal roles in determining bolting and flowering time. Moreover, a schematic regulatory network of radish bolting and flowering was put forward for characterization of DEGs and miRNAs. These results provided essential information for genetic control of radish bolting and flowering, and would facilitate unraveling the molecular regulatory mechanism underlying bolting and flowering in radish and other root vegetable crops.

AUTHOR CONTRIBUTIONS
SN, CL, and LL designed the research. SN, XS, and MT conducted experiments. SN, LX, and YW participated in the design of the study and performed the statistical analysis. SN analyzed data and wrote the manuscript. LL and EM helped with the revision of manuscript. All authors read and approved the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00682 Table S1 | All the identified DEGs in NAU-VS and NAU-RS libraries.