Induction of Glucoraphasatin Biosynthesis Genes by MYB29 in Radish (Raphanus sativus L.) Roots

Glucoraphasatin (GRH) is a specific aliphatic glucosinolate (GSL) that is only abundant in radish (Raphanus sativus L.). The gene expression regulating GRH biosynthesis in radish is still poorly understood. We employed a total of 59 radish accessions to analyze GSL profiles and showed that GRH was specific and predominant among the aliphatic GSLs in radish roots. We selected five accessions roots with high, moderate and low GSL biosynthesis, respectively, to conduct a comparative transcriptome analysis and the qRT-PCR of the biosynthesis genes for aliphatic GSLs. In this study, among all the accessions tested, roots with the accession RA157-74 had a high GRH content and showed a significant expression of the aliphatic GSL biosynthesis genes. We defined the genes involved in the GRH biosynthesis process and found that they were regulated by a transcription factor (RSG00789) at the MYB29 locus in radish roots. We found 13 aliphatic GSL biosynthesis genes regulated by the RSG00789 gene in the GRH biosynthesis pathway.


Introduction
Glucosinolates (GSLs) are secondary metabolites containing nitrogen and sulfur, mainly found in plants. Approximately 200 GSLs are known to exist naturally in plants [1]. GSLs are classified into several groups by the structure of the side chain [2]. According to the modified amino acid content, three main groups exist: aliphatic, indolic and aromatic. Aliphatic GSLs are derived from methionine, while indolic and aromatic GSLs originate from tryptophan and phenylalanine, respectively [3].

RNA Sequencing and Mapping of Radish Reference Genomes
Among the radish accessions, RA502-92 roots were the low-abundance GRH accessions (LGRHA), and RA157-74 roots were the high-abundance GRH accessions (HGRHA). The GRH level was on average 791 times higher in the HGRHA roots when compared with the LGRHA roots. We selected five accessions for a comparative analysis of GRH biosynthesis. IT119282-15, IT119238-8 and RA280-82 roots were analyzed as the moderate GSL biosynthesis samples (Table 2 and Figure 2).
RNA sequencing was performed with a total of 15 radish root samples (three biological repeats from each of the above five accessions). 25,046,852~59,463,280 reads were produced from each sample after quality trimming. The Q30 values were from 97.38-98.28%, and the overall GC percentages were 44.78-47.16%. The average mapping ratio was 68.31% (Table S3).

RNA Sequencing and Mapping of Radish Reference Genomes
Among the radish accessions, RA502-92 roots were the low-abundance GRH accessions (LGRHA), and RA157-74 roots were the high-abundance GRH accessions (HGRHA). The GRH level was on average 791 times higher in the HGRHA roots when compared with the LGRHA roots. We selected five accessions for a comparative analysis of GRH biosynthesis. IT119282-15, IT119238-8 and RA280-82 roots were analyzed as the moderate GSL biosynthesis samples (Table 2 and Figure 2).
RNA sequencing was performed with a total of 15 radish root samples (three biological repeats from each of the above five accessions). 25,046,852~59,463,280 reads were produced from each sample after quality trimming. The Q30 values were from 97.38-98.28%, and the overall GC percentages were 44.78-47.16%. The average mapping ratio was 68.31% (Table S3).

Analysis of Differentially Expressed Genes between Roots of the Radish Accessions
To investigate the expression level of genes involved in GSL biosynthesis, transcriptome data from the LGRHA and HGRHA roots were used for the analysis of differentially expressed genes (DEGs). A total of 2597 genes were up-regulated in the HGRHA roots when compared with the LGRHA roots, while 2417, 1654 and 1916 genes were up-regulated in the HGRHA roots when compared with the IT119282-15, IT119238-8 and RA280-82 accessions, respectively ( Figure S1).
Gene ontology (GO) enrichment analysis was performed to classify the predicted functions of DEGs that were up-regulated in the HGRHA when compared with the LGRHA roots. The DEGs were classified into three main GO groups, consisting of biological processes, cellular component, and molecular function, including abundant GO terms. The top five GO terms in the biological processes group were the hydrogen sulfide biosynthetic, glucosinolate biosynthetic, sulfur compound biosynthetic, glycosyl compound biosynthetic and amine biosynthetic processes. These highly linked GO terms were related to biosynthesis processes that stimulate biological responses ( Figure 3 and Table S4). Similar results were obtained when comparing the HGRHA roots with the IT119282-15, IT119238-8 and RA280-82 accession roots (Tables S5-S7).
were classified into three main GO groups, consisting of biological processes, cellular component, and molecular function, including abundant GO terms. The top five GO terms in the biological processes group were the hydrogen sulfide biosynthetic, glucosinolate biosynthetic, sulfur compound biosynthetic, glycosyl compound biosynthetic and amine biosynthetic processes. These highly linked GO terms were related to biosynthesis processes that stimulate biological responses ( Figure 3 and Table S4). Similar results were obtained when comparing the HGRHA roots with the IT119282-15, IT119238-8 and RA280-82 accession roots (Tables S5-S7). LGRHA roots. Up-regulated genes in the HGRHA roots, when compared with the LGRHA roots, were divided into three main categories (biological processes, cellular component, and molecular function). The GSL biosynthetic process ranked second in the biological processes category. All GO terms were filtered by the following values; p-value < 0.01, false discovery rate < 0.01, fold-change of selected genes to total genes > 2. The x-axis indicates the fold-change as a% of total assigned transcripts. Detailed information is listed in Table S4.

Identification and Expression Analysis of GSL Biosynthesis Genes in Radish Roots
We selected a total of 64 genes with paralogues involved in GSL biosynthesis by a functional annotation in radish roots (Table S8), and the expression values (FPKM) of these genes are shown in the heat map in Figure 4B. Major genes involved in aliphatic GSL biosynthesis were highly expressed in the HGRHA roots ( Figure 4B).
MYB28 and MYB29 transcription factors are mainly defined as regulators for aliphatic GSL biosynthesis genes in Arabidopsis [6]. Three MYB28 and two MYB29 transcription factors were annotated (Table S8). The MYB29 gene (RSG00789) was significantly differentially expressed in the LGRHA roots. Up-regulated genes in the HGRHA roots, when compared with the LGRHA roots, were divided into three main categories (biological processes, cellular component, and molecular function). The GSL biosynthetic process ranked second in the biological processes category. All GO terms were filtered by the following values; p-value < 0.01, false discovery rate < 0.01, fold-change of selected genes to total genes > 2. The x-axis indicates the fold-change as a% of total assigned transcripts. Detailed information is listed in Table S4.

Identification and Expression Analysis of GSL Biosynthesis Genes in Radish Roots
We selected a total of 64 genes with paralogues involved in GSL biosynthesis by a functional annotation in radish roots (Table S8), and the expression values (FPKM) of these genes are shown in the heat map in Figure 4B. Major genes involved in aliphatic GSL biosynthesis were highly expressed in the HGRHA roots ( Figure 4B).
MYB28 and MYB29 transcription factors are mainly defined as regulators for aliphatic GSL biosynthesis genes in Arabidopsis [6]. Three MYB28 and two MYB29 transcription factors were annotated (Table S8). The MYB29 gene (RSG00789) was significantly differentially expressed in the HGRHA roots only. Interestingly, no significant expression of three paralogous MYB28 genes (RSG53581, RSG23384, RSG16088) was observed in the HGRHA roots when compared with all the tested accessions, despite the fact that MYB28s are the major transcription factors for aliphatic GSL biosynthesis genes [25]. The RSG53581 locus was positively expressed in the HGRHA roots ( Figure 4B and Table S8). Two MYB34 paralogues (RSG13843, RSG52510) were identified in radish roots ( Figure 4B and Table S8), and they are reported to regulate the expression of indolic GSL biosynthesis genes [24].
In the Brassicaceae family, the side chain modification in aliphatic GSL biosynthesis is divided into 3, 4 and 5 carbon stages after the core structure formation, and then the side chains are modified by activating oxygenation, hydroxylation, alkenylation, benzoylation and methoxylation [28]. We verified six FMO GS-OX genes (RSG51376, RSG57085, RSG57158, RSG15645, RSG14855, RSG27761) encoding flavin-monooxygenease glucosinolate S-oxygenases and a GS-OH gene (RSG00041) belonging to the 2-oxoglutarate and Fe (II)-dependent oxygenase superfamily of proteins, which are responsible for the modification of the aliphatic GSLs [28]. Additionally, we found seven CYP81F genes (RSG10180, RSG11506, RSG29212, RSG29215, RSG40296, RSG10178, RSG23271) belonging to CYP81 cytochrome P450 monooxygenases and five IGMT genes (RSG15203, RSG26550, RSG59991, RSG15204, RSG06937) encoding O-methyltransferases for the modification of indolic GSLs. AOP2 genes encoding 2-oxoglutarate-dependent dioxygenases were not detected in this study. However, the expression levels of these genes were very low in all the radish accessions ( Figure 4 and Table S8). It is known that the production of GRH is regulated by a single recessive gene in radish named glucoraphasatin synthase 1 (GRS1), which belongs to the 2-oxoglutarate and Fe (II)-dependent oxygenase superfamily [29,30]. The GRS1 gene (RSG02297) was identified by functional annotation and was strongly expressed in the HGRHA roots (Table 3 and Figure 5A) when compared with all the tested accessions. It has been reported that the accumulation of GRH occurs via the conversion of GER, a process which is catalyzed by GRS1 ( Figure 5B) [29,30]. We examined the expression of 17 DEGs using qRT-PCR. The oligo sequences used for qRT-PCR are listed in Table S1, and RsActin was used as an internal control for the relative expression analysis. The qRT-PCR analysis strongly supported our RNA sequencing data ( Figure 6).

GRH Content and Expression of GRH Biosynthesis Genes in Radish Roots
The expression of genes involved in GRH biosynthesis was significantly high in the HGRHA roots when compared with all the tested accessions (Figures 4-6 and Table S8). We would expect that the GRH content would be high in radish roots with a high expression of GRH biosynthesis genes. However, unexpected expression patterns were revealed among the radish accessions. For example, the GRH content was 274 times higher in the IT119238-8 accession roots when compared with the LGRHA roots (Table 2 and Figure 2). However, the expression levels of GRH biosynthesis genes were low in both accessions, with none being significantly differentially expressed. In fact, several genes were down-regulated in the IT119238-8 accession roots. There was a similar result between the LGRHA and RA280-82 accession roots (Figure 7 and Table S9). RA166-77 accession roots had a high GRH content (98.42 µmol·g −1 dry wt.), which was not significantly different from the HGRHA roots ( Figure 8A). However, the GRH biosynthesis genes were expressed at very low levels in the RA166-77 accession roots when compared with the HGRHA roots. In fact, the expression of these genes in the RA166-77 accession roots was similar to the expression levels in the LGRHA roots ( Figure 8B and Table S10). These results indicated that the GRH content was not dependent on the expression level of GRH biosynthesis genes in the radish roots. Figure 6. qRT-PCR examination of aliphatic GSL biosynthesis genes in radish roots. The GSL biosynthesis genes were analyzed by qRT-PCR, and the value was compared with the FPKM value. The blue lines represent the FPKM (left y-axis) data, and the red bars indicate the qRT-PCR data (right y-axis). The expression level of qRT-PCR was measured by the delta CT method.

GRH Content and Expression of GRH Biosynthesis Genes in Radish Roots
The expression of genes involved in GRH biosynthesis was significantly high in the HGRHA roots when compared with all the tested accessions (Figures 4-6 and Table S8). We would expect that the GRH content would be high in radish roots with a high expression of GRH biosynthesis genes. However, unexpected expression patterns were revealed among the radish accessions. For example, the GRH content was 274 times higher in the IT119238-8 accession roots when compared with the LGRHA roots (Table 2 and Figure 2). However, the expression levels of GRH biosynthesis genes were low in both accessions, with none being significantly differentially expressed. In fact, several genes were down-regulated in the IT119238-8 accession roots. There was a similar result between the LGRHA and RA280-82 accession roots (Figure 7 and Table S9). RA166-77 accession roots had a high GRH content (98.42 µmol·g −1 dry wt.), which was not significantly different from the HGRHA roots Figure 6. qRT-PCR examination of aliphatic GSL biosynthesis genes in radish roots. The GSL biosynthesis genes were analyzed by qRT-PCR, and the value was compared with the FPKM value. The blue lines represent the FPKM (left y-axis) data, and the red bars indicate the qRT-PCR data (right y-axis). The expression level of qRT-PCR was measured by the delta CT method.
( Figure 8A). However, the GRH biosynthesis genes were expressed at very low levels in the RA166-77 accession roots when compared with the HGRHA roots. In fact, the expression of these genes in the RA166-77 accession roots was similar to the expression levels in the LGRHA roots ( Figure 8B and Table S10). These results indicated that the GRH content was not dependent on the expression level of GRH biosynthesis genes in the radish roots.   ( Figure 8A). However, the GRH biosynthesis genes were expressed at very low levels in the RA166-77 accession roots when compared with the HGRHA roots. In fact, the expression of these genes in the RA166-77 accession roots was similar to the expression levels in the LGRHA roots ( Figure 8B and Table S10). These results indicated that the GRH content was not dependent on the expression level of GRH biosynthesis genes in the radish roots.

Abundant GRH Content and Functional Possibility for Radish Breeding
GSLs are sulfur-containing phytochemical compounds mainly found in the Brassicaceae. Radish is a root vegetable belonging to the Brassicaceae family, which has an abundant amount of a specific aliphatic GSL known glucoraphasatin (GRH) [18]. Previous research showed that GRH was the predominant GSL in the roots of eight radish varieties [16]. An analysis of the individual and total GSL contents from 44 Raphanus species revealed that GRH was the predominant GSL in both the leaves and roots of all the tested accessions [17]. The GSL profiling of the roots of 71 radish accessions showed that GRH was the predominant GSL in all accessions, constituting 95.2% of the total GSL content on average [18]. An analysis of the GSL content in Brassicaceae seeds showed that the GRH was only detected in two radish seeds, while it was not detected at all in broccoli (B. oleracea var. italica), rutabaga (B. napus var. napobrassica) and turnip (B. rapa ssp. rapa) seeds [31]. In the seven-day old sprouts of two radishes (Daikon and Sango) and two Brassica varieties (broccoli and Tuscan black kale), GRH was specifically detected in the two radish varieties [32]. Recently, GRH was detected in Chinese cabbage (B. rapa ssp. pekinensis), broccoli (B. oleracea var. italica) and choysum (B. rapa ssp. chinensis), but only at low levels [21,29].
In the present study, the GRH in roots accounted, on average, for 91.3% of the total GSL content in 59 radish accessions. Other GSLs, including GRA, GRE, GNA, GBN, GER, 4HGBS, 4MOGBS, NGBS and GBS, were at very low levels or non-detected in all radish accessions (Figure 1 and Table S2). GRH contributes more to the production of isothiocyanates (ITCs) when compared with GRA, which is the major GSL found in broccoli [32]. Raphasatin (GRH-ITC) is a specific ITC derived from GRH in radish, whereas sulforaphane is derived from GRA in broccoli [22]. These ITCs are beneficial to human health, playing a role in preventing inflammation and carcinogenesis, and aiding in cardiovascular protection [18,33,34]. It is known that GRH-ITC possesses a selective cytotoxic/apoptotic activity on three human colon carcinoma cell lines [22].
GRH was significantly abundant and the predominant GSL in the roots of the radish accessions that were tested. The RA157-74 accession roots contained a significantly high GRH content; this HGRHA will be exploited in an expression study of GSL biosynthesis genes, as well as to improve a new variety of healthy resources.

Induction of GRH Biosynthesis by MYB29 Transcription Factor in Radish Roots
The aliphatic GSL biosynthesis process consists of side-chain elongation, core structure synthesis and side-chain modification in the Brassicaceae family. BCAT4, BCAT3, MAM1, IPMI_LSU1, IPMI_SSU2, IPMDH1, BAT5, CYP79F1, CYP83A1, SUR1, UGT74B1, UGT74C1, SOT17, SOT18, MYB28 and MYB29 genes are known to be the major genes involved in aliphatic GSL biosynthesis [6,11,12,35,36]. Furthermore, the GRS1 gene has been reported to contribute to the accumulation of GRH in radish [29,30]. We measured the expression of these genes by RNA sequencing (Table S8). In previous studies, radish mutants lacking GRH, grs1-1 and 1-2 mutants, were grafted with a cultivar, and the majority of aliphatic GSLs in the radish roots were transported from the leaf to the root [30]. The GRS1 gene, which is responsible for GRH production, is predominantly expressed in the leaf development stages (8 and 12 weeks old) of radish [29]. These studies indicate that the radish root was not the main site of GRH biosynthesis.
The genes involved in aliphatic GSL biosynthesis, including GRS1, were significantly expressed in the HGRHA roots when compared with all the tested accessions (Figures 4-8). It was predicted that these genes would be induced by MYB29 (RSG00789). However, the MYB29 transcription factor is not the dominant regulator, only playing a minor role in the induction of aliphatic GSL biosynthesis genes [25]. In Arabidopsis, a knockout mutant of MYB29 did not influence the expression of aliphatic GSL biosynthesis genes, whereas the MYB28 knockout mutants had significantly decreased levels or lacked these genes [25]. In this study, the expression of MYB28 (RSG53581) was high in the IT119282-15 accession roots, but it had no significant effect on the expression of GRH biosynthesis genes. Meanwhile, one MYB29 paralogue, RSG00789, was strongly amplified in the HGRHA roots (Figures 6-8).
Several biotic elicitors, such as the phytohormones methyl jasmonate (MeJA), jasmonic acid and ethylene, increase the GSL content in the Brassicaceae family [31,[37][38][39]. Aliphatic GSL biosynthesis has been known to be increased by MeJA in Arabidopsis. Under normal conditions, MYB28 is more highly expressed than MYB29. MYB29 probably plays a supplementary role in response to MeJA signaling in Arabidopsis [40]. Our study revealed that some genetic factors influenced the expression of RSG00789, which probably led to the overexpression of GRH biosynthesis genes in HGRHA roots.

Expression Profiling of GRH Biosynthesis Genes in Radish Roots
In the side chain elongation step of aliphatic GSL synthesis, the precursor amino acid is deaminated to 2-oxo acid by BCATs. 2-oxo acid then undergoes a step involving three continuous transformations. Condensation with acetyl-CoA by the MAMs, isomerization by IPMIs and oxidative decarboxylation by IPMDH are known processes involved in the side chain elongation step in the chloroplast [28]. BCAT4, MAM1, IPMI_LSU1, IPMI_SSU2 and IPMDH1 genes were significantly differentially expressed in the HGRHA roots when compared with all tested accessions, but not with BCAT3 ( Figures 4B, 6 and 7). In Arabidopsis, BCAT3 has been shown to change the proportions of aliphatic GSLs in the terminal step of the side chain elongation process in response to amino acid biosynthesis [26]. Arabidopsis BCAT3 mutation did not decrease the total GSL content derived from methionine. The BCAT3 gene mainly plays a role in leucine biosynthesis but can also act as a supporter of BCAT4 by transaminating methionine to 2-oxo acid in the side chain elongation step [26]. Expression of the BCAT3 gene (RSG45676) is not necessary for GRH biosynthesis but seems to participate in leucine biosynthesis in radish roots [26]. All enzymes involved in the side chain elongation step are located in the chloroplast, except for BCAT4. Transporters are necessary for the import of 2-oxo acids into the chloroplast and for the export of chain-elongated amino acids to the cytosol. The BAT5 gene encodes a bile acid transporter [28]. Both BAT5 and MYB29 are commonly induced by MeJA in Arabidopsis [41]. MYB29 (RSG00789) and BAT5 (RSG43366) were significantly expressed in the HGRHA roots when compared with all the tested accessions ( Figure 7 and Figure S5).
During the core structure synthesis process, precursor amino acids are converted into aldoximes by CYP79F1, as seen in Figure 4. Next, the aldoxime is converted to S-alkyl-thiohydroximate by CYP83A1, and then the basic GSL skeleton is formed by the sequential activities of SUR1, UGT74B, UGT74C, SOT17 and SOT19 [6,42]. The two SOTs are responsible for the glucosylation of DS-GSLs [28]. Two genes, SOT17 (RSG03557) and SOT18 (RSG07388), were highly expressed in the HGRHA roots when compared with all the tested accessions ( Figure 6 and Figure S4). The SOTs have different substrate specificities. The chain length of aliphatic DS-GSLs is the important element for substrate affinities of SOTs in Arabidopsis [28,43]. SOT17 and SOT18 prefer long-chain DS-GSLs derived from methionine, such as 7-Methylthioheptyl-GSL (7MTH) and 8-Methylthiooctyl-GSL (8MTO), but they are also active with short-chain aliphatic methionine-derived GSLs, such as 3-Methylthiopropyl-GSL (3MTP) and 4-methylthiobutyl-GSL (4MTB). SOT18 showed a broader substrate specificity than SOT17 in Arabidopsis [43]. GRH, the predominant GSL in radish, is a short-chain GSL derived from 4-methylthiobutyl (4MTB). The expression of RSG03557 was significantly specific for 4MTB when compared with RSG07388 in radish roots.
In the Brassicaceae family, FMO GS-OXs are responsible for the modification of GSL side chains to produce methylsulfinyl types such as GRA, which are derived from methylthio types such as GER by oxygenation [29]. Seven FMO GS-OX genes have been reported in Arabidopsis [29]. We identified a total of six FMO GS-OX paralogous genes in radish. The expression of these genes was very low or not detected in the HGRHA roots, even though GRH is a methylsulfinyl-GSL ( Figure 4 and Table S8). Recently, it was reported that GRS1 is directly responsible for GRH production [30], and we found that GRS1 was highly expressed in the HGRHA roots when compared with all of the tested accessions ( Figure 5).
In mature radishes, most genes involved in aliphatic GSL biosynthesis are found in the leaf. These genes are not expressed much in mature radish roots [11,29]. However, GSL biosynthesis genes were strongly expressed in the HGRHA roots. These results suggest that GRH can be significantly synthesized in the mature radish root under special conditions and that MYB29 (RSG00789) plays an important role in this process.
Based on these results, we outlined a putative GRH biosynthesis process in HGRHA roots in a schematic diagram (Figure 9).
Based on these results, we outlined a putative GRH biosynthesis process in HGRHA roots in a schematic diagram (Figure 9).   [13]. For the analysis of the GSL content, seven seeds from each accession were sown in plastic pots, after which the seedlings were transferred into the field. After eight weeks, taproots from each mature radish were collected for an analysis of the GSL content and RNA-sequencing.

Determination of GSL Content
GSLs were extracted using the method described by Lee et al., 2018 [5]. A total of 100 mg of freeze-dried powder with 1.5 mL of 70% (V/V) MeOH was incubated in a water bath at 70 °C for 5 min. After centrifugation at 4 °C for 10 min (4000 rpm), the supernatants were harvested into a 15 mL tube. The supernatants were used as crude GSL extracts. These extracts were applied to DEAE-Sephadex A-25 (GE Healthcare, USA) in a minicolumn and were then desulfated using 75 µL of aryl sulfatase solution. Desulfo-glucosinolate (DS-GSL) was eluted with 1.5 mL of ultrapure water. The separation of DS-GSLs using a reversed phase Inertsil ODS-3 column (150 mm × 3.0 mm × 3 µm) was carried out with an E type cartridge guard column (10 mm × 2.0 mm × 5 µm) in an Agilent Technologies 1100 series HPLC system (Agilent Technologies, Germany) under the following conditions: the detection wavelength was 227 nm, column oven temperature was 40°C and flow rate was 0.2 mL/min. The solvent A (water) and solvent B (acetonitrile) consisted of the mobile phase. Elution of the DS-GSL samples was performed under the following gradients: 0 min, 100% A/0% B; 2 min, 100% A/0% B; 7 min, 90% A/10% B; 16 min, 69% A/31% B; 19 min, 69% A/31% B; 21 min, 100%   [13]. For the analysis of the GSL content, seven seeds from each accession were sown in plastic pots, after which the seedlings were transferred into the field. After eight weeks, taproots from each mature radish were collected for an analysis of the GSL content and RNA-sequencing.

Determination of GSL Content
GSLs were extracted using the method described by Lee et al., 2018 [5]. A total of 100 mg of freeze-dried powder with 1.5 mL of 70% (V/V) MeOH was incubated in a water bath at 70 • C for 5 min. After centrifugation at 4 • C for 10 min (4000 rpm), the supernatants were harvested into a 15 mL tube. The supernatants were used as crude GSL extracts. These extracts were applied to DEAE-Sephadex A-25 (GE Healthcare, USA) in a minicolumn and were then desulfated using 75 µL of aryl sulfatase solution. Desulfo-glucosinolate (DS-GSL) was eluted with 1.5 mL of ultrapure water. The separation of DS-GSLs using a reversed phase Inertsil ODS-3 column (150 mm × 3.0 mm × 3 µm) was carried out with an E type cartridge guard column (10 mm × 2.0 mm × 5 µm) in an Agilent Technologies 1100 series HPLC system (Agilent Technologies, Germany) under the following conditions: the detection wavelength was 227 nm, column oven temperature was 40 • C and flow rate was 0.2 mL/min. The solvent A (water) and solvent B (acetonitrile) consisted of the mobile phase. Elution of the DS-GSL samples was performed under the following gradients: 0 min, 100% A/0% B; 2 min, 100% A/0% B; 7 min, 90% A/10% B; 16 min, 69% A/31% B; 19 min, 69% A/31% B; 21 min, 100% A/0% B; and 27 min, 100% A/0% B. The classification of GSLs was dependent upon their HPLC retention time with our database. Individual GSLs were quantified using an external standard sinigrin (0.1 mg/mL for its desulfation).

Total RNA Isolation and RNA-Seq Analysis
The total RNAs were isolated from the roots using a Hybrid-R kit (GeneAll, Songpa-gu, Korea) according to the manufacturer's instructions, and their quality and quantity were examined using a Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). About 2 µg of total RNAs were used to construct RNA-Seq libraries with an insert size of 300 bp, using an Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA), according to the manufacturer's instructions. Pooled libraries were sequenced using the Illumina HiSeq X platform with paired-end (PE) reads of 101 bp at Macrogen Co. (Seoul, Korea). Low-quality and duplicated reads and adapter sequences were removed using Trimmomatic (ver. 0.38) [44] with default parameters (removing a read with an average base quality below 20 and dropping a read less than 50 bases long). Trimmed high-quality RNA-Seq reads were mapped on the radish genome sequences of the R. sativus var. hortensis cv. Aokubi doubled haploid (DH) line [12] using HISAT2 (https://ccb.jhu.edu/software/hisat2/index.shtml) [45]. Then, RNA reads mapped on protein coding sequences (64,657 for var. hortensis) were counted using HTSeq-count (https://pypi.org/project/HTSeq) [46]. Fragments per Kilobase of transcript per Million mapped reads (FPKM) values were calculated using the number of RNA-seq reads mapped on protein coding sequences and were used for the expression profiling of genes.

Analysis of Differentially Expressed Genes Involved in GSL Biosynthesis
The bioconductor package DESeq (http://bioconductor.org/packages/release/bioc/html/DE-Seq. html) [47] was used to identify differentially expressed genes (DEGs) between samples with a low GSL content, medium GSL content and high GSL content, respectively. Genes showing over two-fold expression changes with an adjusted p-value of less than 0.05 were considered to be DEGs. A gene ontology (GO) enrichment analysis was performed for the DEGs using Fisher's exact test with an adjusted p-value of 0.05 in BLAST2GO software (ver. 5.2, https://www.blast2go.com) and was then filtered by a p-value of 0.01, false discovery rate (FDR) of 0.01 and fold-change of selected genes to total genes > 2 for an analysis of GO terms. Information on GSL biosynthetic genes was obtained from previous studies [11,12], after which a GO analysis was carried out by functional annotation using nr BLAST, InterProScan and Araport11 database. The expression values (FPKM) of the genes were retrieved from the expression profiles of total gene sets and were used to generate heat maps using pheatmap in the R-studio.

Quantitative Real-Time PCR
A total of 1 µg RNA extracted from radish roots was synthesized using the RNA to cDNA EcoDry™ Premix (TaKaRa). Subsequently, 1 µL cDNA with 0.5 µL (10 pmol) of primers (Table S1) and 10 µL Topreal qPCR premix, SYBR Green with low Rox (Enzynomix) was used for the quantitative real-time PCR (qRT-PCR). A 20 µL reaction mixture was subjected to amplification using a LightCycler ® 480 (Roche) under the following conditions. An initial denaturation at 95 • C for 10 min was followed by 45 cycles of denaturation at 95 • C for 10 s, annealing at 60 • C for 15 s and elongation at 72 • C for 15 s. The radish actin gene was also simultaneously amplified as a reference. The quantification of the gene expression was normalized to that of actin and was then calculated using the delta CT method.