Identification of Flowering-Related Genes Responsible for Differences in Bolting Time between Two Radish Inbred Lines

Late bolting after cold exposure is an economically important characteristic of radish (Raphanus sativus L.), an important Brassicaceae root vegetable crop. However, little information is available regarding the genes and pathways that govern flowering time in this species. We performed high-throughput RNA sequencing analysis to elucidate the molecular mechanisms that determine the differences in flowering times between two radish lines, NH-JS1 (late bolting) and NH-JS2 (early bolting). In total, 71,188 unigenes were identified by reference-guided assembly, of which 309, 788, and 980 genes were differentially expressed between the two inbred lines after 0, 15, and 35 days of vernalization, respectively. Among these genes, 218 homologs of Arabidopsis flowering-time (Ft) genes were identified in the radish, and 49 of these genes were differentially expressed between the two radish lines in the presence or absence of vernalization treatment. Most of the Ft genes up-regulated in NH-JS1 vs. NH-JS2 were repressors of flowering, such as RsFLC, consistent with the late-bolting phenotype of NH-JS1. Although, the functions of genes down-regulated in NH-JS1 were less consistent with late-bolting characteristics than the up-regulated Ft genes, several Ft enhancer genes, including RsSOC1, a key floral integrator, showed an appropriate expression to the late-bolting phenotype. In addition, the patterns of gene expression related to the vernalization pathway closely corresponded with the different bolting times of the two inbred lines. These results suggest that the vernalization pathway is conserved between radish and Arabidopsis.


INTRODUCTION
Radish (Raphanus sativus L.), an annual plant belonging to the Brassicaceae family, is a familiar root vegetable crop around the world, especially in eastern Asia. The main edible portion of the plant is the fleshy taproot, which has an upper part originating from the hypocotyl and a lower part consisting of true root tissue. The taproot exhibits huge variation in terms of skin color and shape (Tsuro et al., 2008;Zaki Hem and Yokoi, 2012). Some cultivars of radish have been exploited as oilseed crops, silique vegetables, and leafy vegetables. The taproot is an excellent source of carbohydrates and dietary fiber, high levels of glucosinolates and secondary metabolites for human beings: consequently, the properties of the root are the primary determinants of the crop's economic value (Nakamura et al., 2001;Wang et al., 2002). Comparative genomics revealed that the Brassicaceae species exhibit complex genomic syntenies (Panjabi et al., 2008;Li et al., 2011), indicating that a wide range of whole-genome rearrangements occurred during or after species divergence. By contrast, genome syntenies are well conserved in Poaceae (International Brachypodium, 2010) and Solanaceae crops (Tomato Genome, 2012). The genome of the radish is relatively small and diploid (2n = 2x = 18), and analysis of genetic map construction and comparative mapping based on EST-SSR  or EST-SNP markers (Li et al., 2011) revealed that it is closely related to Brassica rapa. Two whole draft genome sequences of radish were recently reported: the 402 Mb sequence of Japanese cultivar "Aokubi" and 510 Mb sequence of Korean cultivar "WK10039" (Kitashiba et al., 2014;Mun et al., 2015). However, establishment of the more accurate chromosome pseudomolecules by clone-end and BACend sequences remains to be accomplished.
In flowering plants, proper timing of the transition from vegetative to reproductive development is important to ensure reproductive success and seed production (Srikanth and Schmid, 2011); accordingly, this trait has considerable economic value. Regulation of Ft is controlled by a diverse range of environmental and internal signals (Koornneef et al., 1998). Environmental signals, which are strongly influenced by season, day length, and temperature adaptation, are integrated with endogenous signals such as developmental stage and age (Zhai et al., 2015). In the model plant Arabidopsis, the majority of key Ft genes have been identified and characterized through genetic and molecular analyses. These genes can be categorized into five main flowering pathways: the photoperiod, vernalization, gibberellin, autonomous, and endogenous pathways (Putterill et al., 2004;Fornara et al., 2010;Srikanth and Schmid, 2011). However, these genetically established pathways are not rigorously distinguished. Instead, a growing body of evidence indicates that extensive crosstalk occurs between the pathways, and that the resultant signals are ultimately integrated by a small number of common target genes (central floral pathway integrator genes) that quantitatively regulate the development of shoot apical meristem. Consequently, control of Ft is plastic and diverse (Simpson and Dean, 2002). Functional genetic analysis in Arabidopsis revealed that the FLOWERING LOCUS T (FT), SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1), and LEAFY (LFY) genes are central floral integrators (Blazquez et al., 1998;Nilsson et al., 1998;Samach et al., 2000). FT has been known for florigen, a mobile signal that transfers the flowering induction signal from leaves into shoot apical meristems through the phloem (Tamaki et al., 2007;Turck et al., 2008), and this protein is highly conserved in most flowering plants (Andrés and Coupland, 2012). SOC1 encodes a MADS box protein that acts as a floral activator to control floral patterning and floral meristem, as well as Ft (Liu et al., 2007(Liu et al., , 2009Melzer et al., 2008). LFY, a unique plant transcription factor that also contains a MADS box, is indispensable for determining the identity of male and female reproductive organs during flower development (Weigel et al., 1992). Together, these key integrator genes interpret signals from a global network of multiple flowering pathways, and their expression levels precisely modulate the expression of floral meristem-specific genes and ultimately determine the exact Ft (Simpson and Dean, 2002;Parcy, 2005).
Throughout the plant life cycle, bolting and flowering are key life-history traits that exercise far-reaching influence on reproductive suitability, mating patterns and opportunities, gene flow, and evolution (Franks, 2015). Currently, due to global warming, a shift toward early flowering is underway in a wide variety of plant species around the world (Parmesan and Yohe, 2003). Premature bolting and flowering pose a serious problem for production because they decrease both yield and economic value of leafy vegetables; therefore, Ft is a vital trait and a target of selection in crop breeding. A number of studies uncovered the genetic basis of Ft through a quantitative trait locus mapping approach in Brassica crops. The BoFLC2 locus controls Ft in Brassica oleracea (Okazaki et al., 2007). The BrFLC1 and BrFLC2 genes are associated with variation in Ft in B. rapa (Yuan et al., 2009;Zhao et al., 2010). BnFRI.A3 is associated with Ft in Brassica napus (Wang et al., 2011). These studies indicate that homologs of FLC and FRI play vital roles in the control of Ft in Brassica species. In a previous study, we characterized two radish inbred lines, generated by a conventional breeding approach, that exhibit different bolting times under the vernalization conditions. NH-JS1, a late-bolting radish inbred line, is less sensitive to cold treatment than the early-bolting line NH-JS2. Expression of genes related to Ft differs between NH-JS1 and NH-JS2: the repressor RsFLC1 is expressed at higher levels in NH-JS1, whereas the positive integrator RsSOC1 is relatively up-regulated in NH-JS2. Although flowering is the most important agricultural trait in cultivation, Ft genes and their associated mechanisms have not been extensively analyzed at a whole-genome scale in radish. Recently, bolting-related miRNAs and their targets were identified in radish using small RNA libraries from leaves at vegetative and reproductive stages (Nie et al., 2015). In addition, differentially expressed genes (DEGs) involved in the transition from vegetative growth to flowering, as well as flowering itself, were identified in radish by RNA-Seq (Nie et al., 2015).
In this study, we performed RNA-Seq and comparative analysis of transcriptomes between two inbred lines with different bolting times, with the goal of identifying regulators of late flowering under vernalization conditions. Our analysis of DEGs related to Ft indicated that most DEGs involved in the vernalization pathway were expressed in patterns that corresponded closely with the different bolting times of the two lines. Negative regulators of the vernalization pathway, such as RsFLC, RsMAF2, RsSPA1, and RsAGL18, were highly expressed in the late-bolting inbred line NH-JS1, whereas positive regulators of vernalization such as RsVRN1, RsVIN3, and RsAGL19 were relatively highly expressed in the early-bolting line NH-JS2. Based on our findings from RNA-Seq and qPCR analysis, we propose a model regulatory network for the flowering pathway.

Bolting Trait Analysis of Two Radish Inbred Lines
NH-JS1 and NH-JS2, two radish inbred lines developed by NongHyup Seed in Korea (Gyeonggi-do, Anseong, Korea), were used as samples of late-bolting and early-bolting radish plants, respectively. Seeds of each line were sown in sterilized soil and grown under normal growth conditions (23 • C, 16 h light/8 h dark) for 2 weeks. For vernalization treatment, 2-week-old plants were grown in the cold room (5 ± 1 • C, 12 h light/12 h dark) for 15 or 35 days. After the vernalization periods, the plants were transferred to a normal growth room and grown for 30 days under the same condition. The percentage of bolting plants was measured by counting the bolting plants when the length of the floral axis was ≥1 cm. Ten plants were used for each bolting test, with two biological replicates.

Plant Materials and Treatments
To generate samples for RNA-Seq analysis, 2-week-old plants of the NH-JS1 and NH-JS2 inbred lines were exposed to the cold (5 ± 1 • C, 12 h light/12 h dark) for 0, 15, or 35 days. Six samples of shoot tissue were collected from each line at the same point in the light/dark cycle. To obtain adequate RNA for each extraction, shoot tissues from three different plants were pooled. Two biological replicates were performed for each vernalization time point. Thus, a total of 12 samples were collected from the two inbred lines. Following harvest, samples were snapfrozen in liquid nitrogen and stored at −80 • C until further processing. Total RNA isolation from shoot tissue was performed as described (Jung et al., 2014).

RNA-Seq Library Construction and Sequencing
Total RNA was prepared for RNA-Seq library construction. mRNA collected from shoot tissues of the two inbred lines was fragmented and used as a template for synthesis of firststrand cDNA using random hexamers and reverse transcriptase. Second-strand cDNA was synthesized using DNA polymerase I (New England BioLabs, Ipswich, MA, USA) and RNase H (Invitrogen, Carlsbad, CA, USA). The resultant cDNA fragments were purified, end-repaired, polyA-tailed, and ligated to index adapters following the Illumina protocol. The ligation products were amplified by PCR and sequenced on an Illumina HiSeq 2000 sequencer system; a 101 bp paired-end sequencing protocol was employed, and two biological replicates were performed for each sample. All raw read data generated in this study were deposited in the GEO of NCBI under accession number GSE89312

Radish Reference-Guided Assembly and Mapping
Raw sequencing data were filtered using standard RNA-Seq parameters (Illumina pipeline). Adapter contamination, lowquality regions, and N-base reads were trimmed from the raw reads, and then reads with a Phred quality score of 31 (Q ≥ 20) or 25 base pairs (bp) were filtered out. These steps were performed using the DynamicTrim and LengthSort programs of the SolexaQA (v.1.13) package (Cox et al., 2010). The cleaned datasets were pooled and mapped to the radish reference gene set. To obtain assembly results, reference-guided assembly was performed using 46,512 genes from the coding sequence regions of the reference genome (Mun et al., 2015). Using the TopHat program, 71,188 assembled unigenes were identified and imposed on the radish reference gene sets. These datasets were pooled and mapped to the radish reference gene set. Mapping was performed using the Bowtie2 (v2.1.0) program (Langmead and Salzberg, 2012) allowing only unique mapping with a maximum of two mismatches; otherwise, the default options were used. The expression levels in each sample were calculated with an in-house script, and read counts for each gene were normalized against library size and rounded to the nearest whole number.

Functional Annotation
Transcripts from RNA-Seq were validated by comparison with gene sequences in the Phytozome database (http://www.phytozome.net/) using BLASTP with E-values of at least 1E-10 (BLAST v.2.2.28+) (Altschul et al., 1997). The Blast2GO software v2.8.0 was further used to compare transcripts (≥200 bp) to the non-redundant (NR) databases with thresholds of E-value ≤ 1E-05 and 20 BLAST hits (Conesa and Götz, 2008). For Gene Ontology (GO) analysis, the GO database (http://www.geneontology.org/) was downloaded, and the transcripts were annotated to the GO database using BLASTP (E-value ≤ 1E-06). GO term annotation was performed using GO classification results from the Map2Slim.pl script. Protein sequences with the highest sequence similarities and cutoffs were retrieved for analysis. Functional enrichment analysis was carried out using DAVID (http://david.abcc.ncifcrf.gov/) and AgriGO (FDR ≤ 0.05) (Du et al., 2010). The transcript lists were further annotated using the TAIR database, and filtered according to default criteria (counts ≥ 2 and EASE score ≤ 0.1). In addition, KEGG pathways were assigned to the sequences by the single-directional best hit method using the KEGG Automatic Annotation Server (Moriya et al., 2007) and DAVID (counts ≥ 10, FDR ≤ 0.05).

DEGs Analysis
Gene expression data were generated from 12 samples from each of two inbred lines. To identify genes differentially expressed during vernalization treatments of each line, raw counts were normalized and analyzed using the DESeq library in R (v3.0) (Anders and Huber, 2010). For a gene to be considered as a DEG, we required |log 2 (fold change)| ≥ 1. In addition, DEGs were filtered by requiring the adjusted p-value (FDR) to be ≤ 0.01. Two types of comparisons were performed: (1) DEGs between the two lines at each time point of vernalization (0, 15, and 35 days); and (2) expression changes of transcripts between vernalization time points: 15 vs. 0 days, and 35 vs. 0 days. Up-and down-regulated transcripts were subjected to a Venn diagram analysis in R.

Identification of Radish Homologs of FT
To discover genes related to Ft in our transcriptome, a set of 174 Ft genes was selected as a reference set based on previous literature and studies in A. thaliana Amasino and Michaels, 2010). Published sequences were obtained from the TAIR database (http://www.arabidopsis.org/index.jsp) based on the Arabidopsis accession numbers of Ft genes. BLASTn was used to query the 174 Ft genes against the assembled 71,188 genes of radish. Top hits were filtered based on the highest percentage of hit coverage and sequence similarity. Cutoffs were E-values ≤ 1E-25 and identity ≥ 65%. Also, the sequences of Ft genes were used to query the protein sequences of radish using BLASTx methods (E-value ≤ 1E-20, identity ≥ 50%).

Validation of a Subset of DEGs by qPCR
Expression data obtained by RNA-Seq were validated by quantitative real-time PCR (qPCR) using cDNA synthesized from the same RNA samples used for RNA-Seq analysis. Total RNA was treated with DNase I (Thermo Fisher Scientific, Waltham, MA, USA) to remove genomic DNA. cDNA was synthesized using reverse transcriptase (RevertAid First-strand cDNA Synthesis Kit; Fermentas, Burlington, Canada). qPCR was performed on a CFX Connect TM Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using SYBR Premix Ex Taq (TaKaRa, Tokyo, Japan). DEGs and Ft genes were determined by qPCR using gene-specific primers (Table S3) and were normalized against the corresponding level of RsACT1. Three technical repeats were performed for each experiment.

Phenotypic Variation in Bolting Time between Two Radish Inbred Lines
First, we recorded the bolting time of the two inbred lines with or without vernalization. After vernalization treatments (0, 15, or 35 days), plants were grown under standard conditions (16 h light/8 h dark, 23 • C) for 30 days in soil, and then we examined variations in flowering phenotypes between the two lines. In the absence of vernalization, the two lines were phenotypically similar, i.e., they did not bolt at all. By contrast, 30 days after vernalization treatment, the two lines exhibited dramatic differences, with NH-JS1 exhibiting a relatively late-bolting phenotype ( Figure 1A). To further investigate the bolting of the two lines under different conditions, we compared the behavior of the two lines following vernalization treatments of various lengths. When grown for 30 days after 15 days of vernalization treatment, 60% of NH-JS2 bolted, whereas no NH-JS1 plants did. Following 35 days of vernalization conditions, all plants of both lines bolted after 30 days of normal growth. At 20 days of growth after 35 days of vernalization, only 63.5% of NH-JS1 plants bolted, whereas almost all (94.5%) of the NH-JS2 plants did. These results confirmed that NH-JS2 bolts earlier than NH-JS1 in response to vernalization ( Figure 1B).
A phenotype is an observable trait caused by differential expression of genes and/or environmental factors. In both inbred lines, vernalization had the most significant effect on the bolting phenotype during growth. Therefore, we predicted that the expression patterns of Ft genes, which we assumed would be closely related to bolting, would differ between the two lines.

RNA-Seq Analysis
To identify transcriptome changes between the two radish inbred lines, we isolated total RNAs from shoot tissues, subjected them to cDNA library construction, and sequenced the libraries on an Illumina HiSeq 2000 Sequencing System ( Figure S1). In total, 587,119,316 paired-end reads (101 bp in length) were generated from 12 libraries generated from the two inbred lines. Raw reads were subjected to quality control, and adapter sequences and low-quality reads were trimmed out. Over 79% of clean reads satisfied the following criteria: quality score Q < 20; minimum read length ≥ 25 bp. In total, 237,331,729 and 228,587,020 clean reads were obtained from NH-JS1 and NH-JS2, respectively (Table S1); the average length of clean reads was 87.85 bp. To verify the similarity between two replicate data, normalized counts were used to created a plot of pairs repeat samples ( Figure S2). Most samples were highly reproducible between pairs. To assess the proportion of unigenes among the 12 transcriptomes, all clean reads were mapped to the reference set of 71,188 unigenes: 94.35% (241,530,265 reads) reads of NH-JS1 and 94.04% (311,512,975 reads) reads of NH-JS2 mapped to the reference unigenes (Table S2). Among all mapped reads, 284,096,662 (48.45%) reads mapped exactly to the reference unigenes, whereas 268,946,578 (45.75%) reads could map to multiple locations in the reference; only 5% of all reads were unmapped. Overall, these results indicated that all of the transcriptome sets had a high proportion of unigenes and were therefore suitable for further DEG and expression profiling analysis.

DEGs between Two Inbred Lines during Vernalization
The total set of expressed genes was subjected to DEG analysis using the DESeq package in R to identify individual characteristics of the transcriptomes of the inbred lines. DEGs were selected using the following criteria: FDR ≤ 0.01 and |log 2 (fold change)|≥1. In total, 3491 DEGs were identified between the two inbred lines at each vernalization time point (0, 15, and 35 days). Of these, 154 genes were shared among the three time points, whereas 309, 788, and 980 genes were specific to samples subjected to 0, 15, and 35 days of vernalization, respectively (Figure 2A). Moreover, in both inbred lines, we investigated genes differentially regulated between two time points (vs. 0 days). Between 15 days and 0 days, 1077 genes (662 up-regulated, 415 down-regulated) were differentially expressed in NH-JS1, and 2552 genes (1429 upregulated, 1123 down-regulated) were differentially expressed in NH-JS2. Between 35 and 0 days, 3253 (1633 up-regulated, 1620 down-regulated) and 2278 (952 up-regulated, 1326 downregulated) genes were identified as DEGs in NH-JS1 and NH-JS2, respectively ( Figure 2B). In the comparison of DEGs between day 15 and day 0, 2.4-fold more DEGs were identified in NH-JS2 than NH-JS1. By contrast, in the comparison between 35 and 0 days, 1.4-fold more DEGs were identified in NH-JS1.
To identify the sets of DEGs that are most relevant to differences in the lines' responses to vernalization, we performed principal component analysis (PCA). The 0 day datasets were similar between the two lines, whereas the 15 and 35 day datasets were dissimilar. The 15 day set from NH-JS2 was the most different from the others. For each inbred line, the results of PCA analysis were consistent with the number of DEGs between the two time points following vernalization treatment, relative to the non-vernalized sample. This suggests that, in both lines, different vernalization times had significantly different effects on the transcription of a subset of transcripts. Thus, these two inbred lines clearly exhibited different expression patterns in response to vernalization treatments of various lengths, as reported in previous studies (Greenup et al., 2011;Villacorta-Martin et al., 2015).
To further investigate the biological pathways active in the two inbred lines, we assigned DEGs from the 15 vs. 0 days comparisons to pathways in the KEGG database. We identified no significantly enriched pathways in NH-JS1, whereas, in NH-JS2, 94 genes were assigned to three pathways: "glycolysis/gluconeogenesis" (45 transcripts), "pyruvate metabolism" (30 transcripts), and "nitrogen metabolism" (19 transcripts). Next, the DEGs from 35 vs. 0 days comparisons were also mapped to the KEGG database. The 352 DEGs from NH-JS1 and the 168 genes from NH-JS2 mapped to 11 and 6 pathways, respectively ( Table 1). Five pathways were shared between the two inbred lines, whereas six and one were specific to NH-JS1 and NH-JS2, respectively. The six specific pathways in the NH-JS1 were "arginine and proline metabolism, " "tryptophan metabolism, " "glutathione metabolism, " "pentose phosphate pathway, " "photosynthesis2, " and "valine, leucine, and isoleucine biosynthesis." The specific pathway in NH-JS2, "circadian rhythm, " included 29 transcripts. Further investigation of the DEGs in the two inbred lines may provide clues about phenotypic variation in response to vernalization.
To identify DEGs related to the flowering pathway, we identified DEGs between the two inbred lines from among the 218 putative Ft genes. In total, 26 and 58 Ft genes were identified as DEGs at the three vernalization time points in the NH-JS1 and NH-JS2, respectively: 11 and 14 genes at 0 days, 5 and 20 genes at 15 days, and 10 and 24 genes at 35 days ( Figure 3A). Interestingly, 12 transcripts were consistently detected as DEGs at all vernalization time points. Furthermore, in both lines, we identified DEGs between vernalization time points; in the comparison between 15 and 0 days, 35 and 45 genes were up-regulated and 20 and 11 were downregulated in NH-JS1 and NH-JS2, respectively. In the comparison between 35 and 0 days, 46 and 55 transcripts were upregulated, and 23 and 11 were down-regulated in NH-JS1 and NH-JS2, respectively ( Figure 3B). These six sets of Ft-related DEGs were analyzed by PCA. This analysis revealed that the datasets from 0 and 15 days of vernalization were similar between the two inbred lines, whereas the sets from 35 days of vernalization were divergent. The 35 day set from NH-JS2 was the most different from the others ( Figure 3C). These results were consistent with the results of DEG analysis of Ft genes between the two lines. Interestingly, most of the upregulated genes (16 of the 19 genes) in NH-JS1 were repressors of flowering, whereas the down-regulated genes in NH-JS1 consisted of 18 enhancer and 12 repressor genes. RsMAF2 (TBIU060196) and RsSPA1 (TBIU023088), which are repressors of the flowering pathway, were always up-regulated, whereas the positive regulators of flowering RsAGL19 (TBIU033101, TBIU033103, and TBIU040817), RsNFYA4 (TBIU046301 and TBIU046298), and RsSOC1 (TBIU065173 and TBIU057467) were consistently down-regulated in NH-JS1, regardless of vernalization time. In addition, most down-regulated genes in NH-JS1 were flowering enhancers (18 of 30 genes), with the exception of RsELF3 (5 of 30 genes) and RsELF4 (3 of 30 genes) ( Table 3).
To identify Ft genes that were differentially regulated between each time point among the Ft DEGs, we investigated the expression profiles of the Ft DEGs based on RNA-Seq data. Fourteen Ft genes were identified based on sequential expression changes between the lines during vernalization (Table 4). RsNFYA4, RsSOC1, and RsSPA1 exhibited > 2-fold expression differences between the two lines at all-time points. RsFUL, RsMAF3, RsPIF3, and RsSPL4 exhibited expression differences between the lines at 0 days: RsELF4, RsLUX/PCL1, and RsSVP     were differentially expressed at 15 days; and RsCDF3, RsELF3, RsMAF2/AGL31, and RsVIN3 genes were differentially expressed at 35 days. The resultant changes in Ft gene transcription are likely to be involved in differences in bolting in response to vernalization.

Validation of DEGs and Expression Profiling of Ft DEGs by qPCR
To evaluate the DEGs identified by RNA-Seq analysis, we performed qPCR. For this purpose, we selected 12 genes that exhibited the highest difference in expression levels between NH-JS1 and NH-JS2. Six of these transcripts (TBIU039665, TBIU054834, TBIU048562, TBIU062462, TBIU047119, and TBIU057835) were up-regulated, and the other six (TBIU057397, TBIU017689, TBIU016676, TBIU056908, TBIU065173, and TBIU057467) were down-regulated, in NH-JS1 vs. NH-JS2. With the exceptions of TBIU048562 and TBIU016676, all of these genes were strongly influenced by vernalization treatments. The qPCR analysis yielded expression patterns consistent with the RNA-Seq data. In particular, the transcript levels of RsELF3 (TBIU017689) and RsSOC1 (TBIU065173, TBIU057467), which are intimately involved in flowering, were, respectively, 43and 8-fold more highly expressed in NH-JS2 than in NH-JS1 (Figure 4). In addition, we performed qPCR analysis to investigate the expression levels of major Ft-regulating genes, which play essential roles in various flowering pathways. A major determinant of Ft, RsFLC1 (FLOWER LOCUS C; suppressor   of flowering), was reduced during vernalization in both lines. Before vernalization treatment, RsFLC1 was 2-fold more highly expressed in NH-JS1; however, the difference in expression level was gradually reduced after 15 and 35 days of vernalization. As in the Ft DEG analysis, most DEGs involved in the vernalization pathway matched with the different bolting phenotype between the lines. Therefore, we also performed qPCR to investigate the expression levels of major Ft genes involved in the vernalization pathway. Interestingly, RsVRN1 and RsVRN2, which are positive regulators of flowering, exhibited different expression levels and opposing expression patterns between the NH-JS1 and NH-JS2 lines. In NH-JS2, the expression levels of these genes increased after vernalization, whereas in NH-JS1 they remained constant or decreased. VRN1, an essential regulator of the floral transition, is gradually up-regulated by vernalization and negatively regulates VRN2 expression (Yan et al., 2003). In particular, RsVRN1 was 3 to 7-fold more highly expressed in NH-JS2 than NH-JS1 during vernalization, but its expression was not induced in NH-JS1 in response to cold exposure or in the absence of vernalization (Figure 5). Expression levels of RsVIN3, a flowering enhancer, as well as the GA pathway-regulated genes RsGID1A, RsAGL19, and RsNFYA4, were highly induced by vernalization treatment in both lines. Moreover, differences in the expression levels between the two lines were increased up to 3-fold during vernalization treatment in NH-JS2. In addition, we checked the expression levels of photoperiod pathway genes (RsELF3, RsCCA1, RsGI, RsLHY, and RsFPA). These genes were expressed in patterns similar to those previously described, and most were more highly expressed in NH-JS2 than NH-JS1. Expression of RsFPA and RsLHY was induced similarly by vernalization treatment in both lines (Figure 5). To address the role of the major FT-regulating genes, RsFLC and RsSOC1, in the GA-regulated flowering pathway, the gene expression levels were analyzed following exogenous GA treatment in two inbred lines under or not under vernalization conditions. The expression levels of the genes differed very little between the two lines in the early GA response, whereas they significantly differed between the two lines when the bolting phenotype was exhibited. NH-JS2 reacted more sensitively to GA treatment than late-bolting NH-JS1 at 20 days after GA application (data not shown). Highthroughput RNA-Seq analysis is now underway to investigate the difference in response to GA treatment between the two lines.

The Transcriptional Regulatory Network Underlying Differences in Bolting Time between the Two Lines
To elucidate the overall flowering pathway governing differences in bolting time, we propose a model regulatory network based on differences in Ft gene expression (as determined by qPCR) between the two inbred lines. We divided the Ft DEGs into three major flowering pathways: vernalization, photoperiod/circadian, and gibberellin. RsFLC was negatively regulated under vernalization and was expressed at higher levels in the late-bolting NH-JS1, whereas another repressor, RsMAF2, was positively regulated in response to vernalization and was also highly expressed in the NH-JS1 line. On the other hand, the enhancers of the vernalization pathway, RsVRN1, RsVIN3, and RsAGL19, were positively regulated upon vernalization and expressed more highly in the early-bolting NH-JS2 line. RsVRN2, a floral repressor, did not exhibit significant changes in expression. Most DEGs involved in the photoperiod pathway (RsCCA1, RsLHY, RsELF3, and RsGI) were remarkably highly expressed in early-bolting NH-JS2; however, RsCO was not highly expressed. RsNFYA4, an enhancer of the photoperiod pathway, also exhibited elevated expression in the early-bolting line. The DELLA domain protein RsGAI, a negative regulator of flowering, exhibited reduced expression in NH-JS1, and a repressor of DELLA proteins, RsGID1A, was also down-regulated in this line (Figure 6).

DISCUSSION
In this study, we investigated the flowering traits of two radish inbred lines under vernalization conditions. Line NH-JS1 was late-bolting following vernalization treatment, whereas NH-JS2 was early-bolting. Our results confirmed that vernalization was necessary for radish bolting and that the vernalization periods required to induce flowering differed between the two inbred lines (Figure 1). RNA-Seq was conducted to identify important genes involved in late bolting and elucidate the molecular network that regulates the flowering pathway in this crop plant. Using radish shoot tissue from young plants, we analyzed individual transcriptomes during vernalization by high-throughput RNA sequencing. Comparative analysis among the radish transcriptomes revealed significant characteristics that were biologically reproducible ( Figure S2). From 40 to 71 million reads were generated per samples (100 bp, paired-end), representing 11 × coverage of the radish genome. Previous studies reported approximately 70,000 unigenes in radish (Wang et al., 2013a,b;Wu et al., 2015). Consistent with this, our transcriptome analysis identified 71,188 genes despite using only shoot tissue for the RNA-Seq analysis.

Differential Expression of Ft-Related Genes Governs Bolting in Response to Vernalization
Genome-wide DEG analysis revealed that vernalization treatments of 15 and 35 days had significantly different effects on the transcription of a subset of genes between the two inbred lines (Figure 2), and KEGG pathway analysis revealed similar trends ( Table 1). The pathway analysis showed that gene expression related to nitrogen and pyruvate metabolism is significantly up-regulated during vernalization in radish. Further investigation of DEGs in the two inbred lines may provide clues about phenotypic variations related to vernalization. By comparison with known Ft-related genes in Arabidopsis, we identified 218 radish Ft genes in our transcriptome analysis  Table 2), of which 49 were differentially expressed between NH-JS1 and NH-JS2. Thus, about 80% of Ft genes were expressed at similar levels in both lines, suggesting that the 49 Ft-related DEGs provide several hints regarding the molecular basis of late-bolting time in NH-JS1. First, expression of repressors of flowering such as RsFLC1 was elevated in NH-JS1, and 16 of 19 up-regulated Ft genes in NH-JS1 were flowering repressors. In addition, 18 of 30 down-regulated Ft genes in NH-JS1 were flowering enhancers. This expression pattern was well correlated with the NH-JS1 phenotype. However, there were some exceptions. For example, the expression levels of RsELF3 and RsELF4, both of which repress flowering, were reduced in NH-JS1. These two genes are components of the circadian clock (Covington et al., 2001;Doyle et al., 2002), suggesting that circadian rhythms may be partly responsible for the differences in Ft between the two lines. Alternatively, these specific circadian clock components may function differently than they do in Arabidopsis. This idea is supported by the observation that RsELF3 was expressed 20-fold more highly in NH-JS2 than in NH-JS1. In Arabidopsis, ELF3 negatively regulates the expression of CO, FT, and GI (Kim et al., 2005). If RsELF3 has similar functions in Ft control, then the expression level of CO, FT, and GI should be lower in NH-JS2 than in NH-JS1. However, RsCO and RsGI were more highly expressed in NH-JS2 than in NH-JS1, suggesting that RsELF3 may function differently than its Arabidopsis counterpart. Second, pathway analysis revealed that several pathways, including photoperiod and vernalization, were altered in NH-JS1. The expression pattern of genes in the vernalization pathway was well correlated with the NH-JS1 phenotype, suggesting that the radish vernalization pathway is similar to that of Arabidopsis. Interestingly, no autonomous pathway genes were found in the DEG list, even though 23 genes in this pathway were identified in radish ( Table 3).
Regardless of vernalization treatment times, two repressor genes and 10 enhancer genes were up-and down-regulated in NH-JS1 vs. NH-JS2, respectively (Figure 3). In radish, these Ft genes are responsible for regulating bolting under vernalization. The expression profile of Ft genes exhibited that most of common Ft DEGs at any time of vernalization (above 12 genes) are overlapped (12 of 14 genes), indicating they expressed quantitatively and temperature-timely to develop different bolting time in the lines. Furthermore, pathway analyses of Ft DEGs suggested that the expression patterns of genes involved in the vernalization pathway perfectly matched the bolting characteristics of the two lines: seven repressors of the vernalization pathway were highly expressed, and five enhancers involved in the vernalization pathway were less expressed, in the late-bolting NH-JS1 line ( Table 4). The 12 common Ft DEGs were biologically confirmed by qPCR analysis (Figures 4, 5).

Possible Divergence of Flowering-Time Regulators between Arabidopsis and Radish
FRI of Arabidopsis is a well-known repressor of Ft that positively regulates FLC expression (Choi et al., 2011). In our data set, we did not detect a radish FRI homolog. This may be because it is not abundantly expressed, or alternatively because there is simply no FRI homolog in radish. To resolve this issue, we performed BLAST analysis of the whole radish genome (Mun et al., 2015) using Arabidopsis FRI as a query sequence. However, we could FIGURE 4 | Validation of RNA-Seq data by qPCR. cDNA synthesized from total RNA extracted from shoots of NH-JS1 and NH-JS2 inbred lines, following vernalization for 0, 15, or 35 days. (A) Up-regulated and (B) down-regulated genes in NH-JS1 vs. NH-JS2 were selected from RNA-Seq data for validation. qPCR values were normalized against the corresponding level of RsACT1 (actin). For each gene in both the qPCR and RNA-Seq analysis, the expression level from NH-JS1 on day 0 was defined as "1." FIGURE 5 | Differential expression between the two inbred lines of major flowering-time genes during vernalization. cDNA synthesized from total RNA extracted from the shoots of NH-JS1 and NH-JS2 inbred lines after vernalization times of 0, 15, or 35 days. qPCR expression level was normalized against the corresponding level of RsACT1. For each, the expression level from NH-JS1 on day 0 was defined as "1." not identify a candidate FRI ortholog with high homology and similar size to the Arabidopsis gene. One gene had 77% nucleotide sequence identity, but the region of homology spanned only half of the Arabidopsis FRI gene (data not shown). At present, we cannot exclude the possibility that this gene acts like Arabidopsis FRI, but it is also possible that yet another gene serves this function in radish. Recently, Nie et al. reported genes associated with bolting in radish based on de novo transcriptome analysis (Nie et al., 2015). They found one FRI gene in radish, but it is not clear whether its sequence is also present in the current FIGURE 6 | Gene regulatory network controlling flowering time in radish. The schema represents the regulatory network of Ft-related genes in late-bolting NH-JS1, in comparison with NH-JS2, under the differential vernalization conditions. The model is based on data obtained by qPCR. Gene expression levels were normalized to the expression levels in non-vernalized NH-JS1 plants. Red indicates higher expression and blue indicates lower expression relative to NH-JS1 on day 0. Arrows indicate transcriptional activation, whereas bars indicate transcriptional repression.
version of the radish genome. Further research is required to resolve this discrepancy. A similar situation also arose regarding LFY, which is an integrator of Ft signaling (Weigel et al., 1992). LFY expression was not detected in our transcriptome data. As with FRI, this could have been either due to low expression or to the absence of such a gene in the genome. When we tried to find the gene using BLAST, we identified two homologous genes, but their cDNAs (5283 and 2490 bp) were much larger than that of Arabidopsis LFY (1263 bp). Most of the LFY sequence was present in the homologous region of the radish genes. Therefore, there may be some errors in gene prediction. Recently, Nie et al. reported three radish LFY genes (Nie et al., 2015). cDNA cloning and detailed analysis may help resolve the apparent discrepancy between our studies.

A Gene Regulatory Network for Control of Bolting Time in Radish
Three major flowering pathways, which differ in their response to vernalization treatment, have been defined, namely, vernalization, photoperiod, and gibberellin, based on Ft DEGs (Figure 6). Vernalization promotes flowering in response to a prolonged period of growth at low temperature, and we detected five Ft DEGs related to the vernalization pathway. FLOWERING LOCUS C (FLC, TBIU004737) is central to the vernalization process (Sung and Amasino, 2005). Two floral repressors, RsFLC and RsMAF2, an RsFLC homolog, were less expressed in the early-bolting NH-JS2, whereas three enhancers of the vernalization pathway were highly expressed in the early-bolting line. We propose that the vernalization pathway is closely associated with the difference in bolting time between the two inbred lines. Ft DEGs in the photoperiod pathway, CCA1, LHY, ELF3, and GI, play a role in facilitating the expression of CONSTANS as floral enhancers (Sawa et al., 2007;Imaizumi, 2010). Although RsCO was weakly expressed in the two lines, most genes related to the photoperiod pathway were highly expressed in the late-bolting line. Similar to the vernalization pathway, the photoperiod pathway showed coincidence with the bolting phenotype between the two lines. In Arabidopsis, the gibberellin pathway promotes flowering by up-regulating the SOC1 genes (Bernier and Périlleux, 2005). GA hormone signaling occurs through proteolytic and non-proteolytic mechanisms when the GA receptor GID1 binds to GA (Murase et al., 2008). In proteolytic GA signaling, GID1 binds to negative regulators of GA responses called DELLA proteins, and the GID1-GA-DELLA complex is formed for destruction via the ubiquitin-proteasome pathway. Both radish DELLA protein and GA receptor genes, RsGAI and RsGID1A, were less expressed in the late-bolting line. However, further DEG analysis related to the GA pathway is needed to understand the relevancy to bolting time. Autonomous pathway-related genes were expressed at very low levels in all samples, suggesting that this pathway may not influence the difference in Ft between NH-JS1 and NH-JS2. These Ft DEGs from the flowering pathways converge on a key floral integrator, RsSOC1, and ultimately regulate LFY to determine the timing of the floral transition (Lee and Lee, 2010). However, we did not detect RsLFY in our transcriptomic analysis. Another floral integrator, RsFT, was also expressed at a low level, possibly because the samples we analyzed were derived from early-stage shoots. RsSOC1 expression differed significantly between the two lines in a manner that depended on vernalization time, and its transcript levels may have an important impact on bolting in these two lines.

CONCLUSIONS
This study demonstrates that differences in flowering traits between two inbred lines were consistent with the expression patterns of flowering-time genes involved in the vernalization pathway. In addition, our comparative transcriptome analysis elucidated the molecular basis of this divergence in bolting time. This is the first genome-wide comparative transcriptome analysis related to flowering traits in the radish.

AUTHOR CONTRIBUTIONS
HC, YK conceived and designed the study and wrote the manuscript. WJ performed the data analysis and wrote the initial manuscript. HP, AL conducted bolting phenotyping and qPCR analysis and wrote the manuscript. SL contributed research materials.

ACKNOWLEDGMENTS
This work was supported by Agricultural Biotechnology Developmental Program (nos. 114061-3 and 115076-2) grants from the Ministry of Agriculture, Food and Rural Affairs and KRIBB Research Initiative Program awarded to HC.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01844/full#supplementary-material Figure S1 | Experimental design and data analysis. Samples of the NH-JS1 and NH-JS2 inbred lines subjected to various vernalization treatments, used for RNA-Seq analysis. Ver, vernalization period; 1st and 2nd, independent biological replicates. The work-flow and algorithm used for detecting DEG and flowering-time genes are shown schematically. Figure S2 | Pairs plot between the two replicate radish transcriptomes.
Table S1 | Summary of RNA-Seq in two radish inbred lines. Table S2 | Results of read mapping in radish transcripts. Table S3 | Gene-specific primers used for qPCR analyses.