For transcriptome analysis, the Chinese cabbage (B. rapa ssp. pekinensis) ‘Chiifu’ inbred line and yellow sarson (B. rapa ssp. trilocularis) ‘LP08’ line were used in this study. For the mapping population, 151 individuals in the F5 generation of a Brassica rapa segregating population were obtained by crossing the ‘Chiifu’ and ‘LP08’ lines in the greenhouse of Rural Development Administration of Korea as described previously (Kim et al., 2017; Kim et al., 2022). This population was used for QTL mapping of GSL contents.

Considering that GSL content is highly affected by a diversity of environmental stresses, in this study, we wanted to exclude environmental stress factors as much as possible and identify genetic factor(s) contributing to different GSL amounts between Chiifu and LP08. Thus, we grew seedling plants in a plant culture dish containing MS media in an environment-controlled growth chamber (22°C, 16h light/8h dark photoperiod). For quantification of GSL contents, 97 F₅ seedlings of recombinant inbred lines (RILs), as well as the parental lines ‘Chiifu’ and ‘LP08’, were used for extraction of GSLs. Seeds were surface-sterilized, spread on half-strength Murashige and Skoog agar medium and stored in the dark at 4°C for 3 days for stratification. Seedlings were grown for 1 week in a growth chamber at 22°C under long-day (LD) (16-h light/8-h dark) conditions. At least 5 seedling plants per line were harvested for high-performance liquid chromatography (HPLC) analysis.

All 151 F₅ lines, as well as the two parental lines, were subjected to GBS using the Illumina NextSeq500 sequencing platform, as described previously (Kim et al., 2022). GBS libraries were sequenced on the NextSeq500 system (Illumina, USA) and single-end reads 150 bp in length were obtained. After sequencing of GBS libraries, the raw reads were de-multiplexed to sort samples using the GBSX tool (Herten et al., 2015). The Brassica rapa reference genome (Brassica_rapa.Brapa_1.0.dna.toplevel.fa) was obtained from the EnsemblPlants genome database (https://plants.ensembl.org/index.html). After de-multiplexing, single-end reads were mapped to the B. rapa reference genome using Bowtie2 (Langmead and Salzberg, 2012). For calling of single nucleotide polymorphism (SNP) variants, the Genome Analysis Toolkit (GATK) and Picard tools (McKenna et al., 2010) packages were used. Local realignment of reads was performed to correct misalignment resulting from the presence of indels, using the GATK ‘RealignerTargetCreator’ and ‘Indoleamine’ sequence data processing tools. Lastly, the GATK ‘HaplotypeCaller’ and ‘SelectVariants’ tools were used for SNP variant calling.

The variant call format SNP data were transformed into the input format using customized code for the R/Qtl package (eGenome, Republic of Korea). Markers with a duplicated pattern or distorted segregation ratio, as estimated using the Chi-square test with Bonferroni correction, were removed. The est.map function of R/Qtl was used to construct the linkage map, with the Kosambi map function converting genetic distances into recombination fractions (p-value threshold of 1e-06, and EM iterations in 1000 times). Subsequently, the composite interval mapping function of the R/Qtl package was used for QTL mapping with the Kosambi function. Regions inferred from peak positions with local maximum limit of detection (LOD) values exceeding the threshold determined using 1,000 permutation tests were considered significant QTLs. To detect QTLs responsible for the differing GSL profiles of the two parental lines, the GSL contents of 97 F₅ lines and the parental lines were measured and applied to the SNP linkage map. The LOD score threshold value for significance (ɑ = 0.05) was estimated based on 1,000 permutation tests. Peaks exceeding the estimated LOD threshold value were selected for further analysis.

Plants were grown at 22°C under a 16-h light/8-h dark photoperiod. One-week-old plants were harvested at ZT4 (4 h after light on) and then immediately ground in liquid nitrogen. GSLs were extracted as desulfo-glucosinolates (DS-GSLs), as reported previously (Han et al., 2019). Approximately 500 mg of fresh sample was incubated with 70% methanol at 70°C for 25 min to deactivate myrosinases. Contents of DS-GSLs were analyzed through ultra-high-performance liquid chromatography (UHPLC; 3000 U-HPLC system; ThermoFisher Scientific, USA). As a standard, sinigrin (0.5mg/ml) injection was used in all analyses (Sigma-Aldrich, USA). The DS-GSLs were resolved with a C18 reverse phase column (Zorbax XDB-C18, 4.6 × 250 mm, 5-µm particle size; Agilent, USA) with a water and acetonitrile gradient system. Samples were injected and maintained at a flow rate of 0.5 ml/min. Peaks were identified using standard compounds (Phytoplan, Germany). The samples were analyzed independently (three replicates) and the results are presented in nmol/g based on fresh weight (FW).

DS-GSLs were analyzed using the Accela UHPLC system (ThermoFisher Scientific) fitted with an ion trap mass spectrometer (LTQ Velos Pro; ThermoFisher Scientific). The samples were resolved using a C18 reverse phase column (Zorbax XDB-C18, 4.6 × 250 mm, 5-μm particle size; Agilent) with water and acetonitrile as the mobile phase, and measured in negative ion mode ([M-H]-). Mass spectrometry was conducted with the following settings: capillary temperature, 275°C; capillary voltage, 5kV; source heater temperature, 250°C; sheath gas flow, 35 arb; auxiliary gas flow, 5 arb; and spectral scanning range, m/z 100–1,500.

Total RNA was extracted from 1-week-old seedlings grown at 22°C under light (16-h light/8-h dark photoperiod) or dark conditions. Three biological replicates were harvested at each time point and frozen in liquid nitrogen. Total RNA was extracted using the RNeasy Plant Mini Kit (QIAGEN, Germany) and subsequently treated with DNase I (NEB, USA) to remove contaminating genomic DNA. Purified total RNA was used for construction of RNA-seq libraries using the TruSeq RNA Sample Preparation Kit (Illumina) according to the manufacturer’s instructions. Paired-end sequencing was performed on the HiSeq 2500 system (Illumina).

Prior to alignment of RNA-seq reads to the Brassica rapa reference genome, FastQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) was employed to evaluate the quality of the RNA-seq reads. Reads with > 90% Q values > 30 were filtered and used only for genome alignment. The B. rapa FASTA genome file (Brassica_rapa.Brapa_1.0.dna.toplevel.fa) and gff3 file (Brassica_rapa.Brapa_1.0.54.gff3) were downloaded from the EnsemblPlants genome database. TopHat2 mapping software was employed with the default parameters for alignment of reads to the reference genome (Kim et al., 2013). Digital read counts were obtained using featureCounts (Liao et al., 2014) and subsequently analyzed for differentially expressed genes (DEGs) using edgeR. The cut-off for DEGs was a two-fold difference (p < 0.05). A multi-dimensional scaling (MDS) plot, correlation map, and PlotSmear results were produced using R software (ver. 3.6.1; https://www.rstudio.com/products/rpackages/) packages. The web-based tool VENNY was used to generate a Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). The web-based tool ShinyGO (ver. 0.61) was used for Gene Ontology (GO) enrichment analysis (http://bioinformatics.sdstate.edu/go/). Hierarchical clustering heatmap analysis was performed using Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm) and the JAVA TreeView program (Saldanha, 2004). Mapping results were visualized with the Integrative Genomics Viewer (IGV) program of the Broad Institute (Thorvaldsdóttir et al., 2013).

Total RNA was extracted from 1-week-old seedlings grown at 22°C under light (16-h light/8-h dark photoperiod) or dark conditions using the RNeasy Plant Mini Kit (Qiagen). Extracted RNA was treated with DNase I (NEB) to remove contaminating genomic DNA. Approximately 5 µg of total RNA was used for cDNA synthesis with EasyScript reverse transcriptase (TransGen Biotech, China). qRT-PCR was performed using Solg 2× Real-Time PCR Smart Mix (SolGent, Republic of Korea) on a LineGene 9600 Plus Real-Time PCR system (BIOER, China), according to the manufacturer’s instructions. qRT-PCR was conducted under the following conditions: denaturation at 95°C for 12 min, followed by 50 cycles of amplification (95°C for 15 s, 60°C for 25 s, 72°C for 35 s) and sampling at 72°C. The relative transcript level of each gene was determined through comparison with that of BrPP2Aa (Bra012474), a housekeeping gene consistently expressed in our RNA-seq analysis. The primers were designed based on sequences obtained from the B. rapa genome database (BRAD; http://brassicadb.cn). The primers used for qPCR analysis are presented in Supplementary Table S6. Student’s t-test was used for statistical analysis (*P < 0.05, **P < 0.01, ***P < 0.001).

PCR amplification was performed to clone two major QTL candidate genes, BrMYB28.1 and BrGSL-OH.1, using gene-specific primers (Supplementary Table S6). PCR reaction was conducted with the following conditions: denaturation at 94°C for 5 min, followed by three initial cycles of amplification (94°C for 30 s, 55°C for 35 s, 72°C for 5 min), 33 additional cycles of amplification (94°C for 30 s, 64°C for 35 s, 72°C for 4 min 30 s) and a final extension step at 72°C for 10 min. The PCR product was extracted from the agarose gel after 1% agarose gel electrophoresis and then purified with the Dyne Power Gel Extraction Kit (Dynebio, Republic of Korea). The purified PCR products were inserted into the pPZP211 plant expression vector using the In-fusion HD cloning kit (Takara Bio, Japan). Entire genomic sequences of BrMYB28.1 and BrGSL-OH.1 were obtained through Sanger sequencing (Bionics, Republic of Korea) using the series of primers listed in Supplementary Table S6.

The vegetable-type ‘Chiifu’ and oilseed-type ‘LP08’ inbred lines grown under LD conditions were harvested to quantify GSL compounds using HPLC. Twelve GSL compounds representing all three GSL groups (aliphatic, indolic, and aromatic) were successfully identified using our detection system (Supplementary Figure S1). Among these 12 GSLs, 7 aliphatic, 4 indolic, and 1 aromatic GSL compound were identified (Supplementary Figure S1 and Table S1). The aromatic GSL compound GNT was detected at a very low level, and was therefore neglected in further analyses. Total GSL levels were apparently higher in ‘LP08’ than ‘Chiifu’ (Figure 1A and Supplementary Table S1). For example, the total GSL levels were 5,434.3 nmol/g FW in ‘LP08’ and 2,869.9 nmol/g FW in ‘Chiifu’; thus, ‘LP08’ had around double the GSL level of ‘Chiifu’ (Figure 1B and Supplementary Table S1). In the LD sample, 94.74% of GSLs in ‘LP08’ and 91.37% in ‘Chiifu’ were aliphatic, while the remaining GSLs were indolic (Figure 1B). While the total amount of aliphatic GSLs in ‘LP08’ was almost double that in ‘Chiifu’, the total amounts of indolic GSLs did not differ dramatically between the two inbred lines, although the difference was found to be statistically significant (Figure 1A and Supplementary Table S1). Taken together, these results indicate that, in B. rapa seedlings, the majority of total GSLs are aliphatic GSLs.

Among aliphatic GSLs, two compounds (GNP and GER: glucoerucin) were present at significantly higher levels in ‘LP08’ than ‘Chiifu’ (Figure 1C). In particular, in ‘LP08’, GNP accounted for an overwhelming proportion of aliphatic GSLs (97.5% of total aliphatic GSLs) (Figure 1D). These results indicate that the higher level of total GSLs in ‘LP08’ is attributable to high abundance of GNP. Meanwhile, among aliphatic GSLs in ‘Chiifu’, five compounds (PGT, GRA, GAS, GNA: gluconapoleiferin, and GBN: glucobrassicanapin) had higher levels in ‘Chiifu’ than ‘LP08’ (Figure 1C). Among these compounds, GRA accounted for 50.0% (1,310.5 nmol/g FW) of total aliphatic GSLs, while GNP and GBN accounted for 30.3% (793.5 nmol/g FW) and 13.2% (347.2 nmol/g FW) of total GSLs, respectively (Figure 1D).

Two indolic GSLs (4-methoxyglucobrassicin: 4-MTGB and NGB) had higher levels in ‘LP08’ than ‘Chiifu’, while two other GSL compounds (glucobrassicin [GBC] and 4-hydroxyglucobrassicin [4-HGB]) were present at higher levels in ‘Chiifu’ than ‘LP08’, demonstrating a dynamic compositional difference of indolic GSLs between the two lines (Figure 1E). In terms of total indolic GSL levels, ‘Chiifu’ and ‘LP08’ had similar amounts (241.0 and 280.5 nmol/g FW, respectively; Figure 1F and Supplementary Table S1). However, detailed analysis of the composition of indolic GSL compounds revealed differing profiles (Figure 1F). For example, in ‘Chiifu’ 4-HGB was dominant (49.2%) among the four measured indolic GSL compounds. Meanwhile, NGB was dominant (55.7%) among the four indolic GSLs in ‘LP08’. This indicates that ‘Chiifu’ and ‘LP08’ have differing compositional profiles for both aliphatic and indolic GSLs. Taken together, these results show that ‘LP08’ had a higher abundance of total GSLs than ‘Chiifu’, which can be attributed to greater accumulation of GNP, whereas several aliphatic GSL compounds were present in large proportions in ‘Chiifu’.

To identify the candidate genes responsible for the observed differences in amounts and profiles of GSLs between ‘Chiifu’ and ‘LP08’, RNA-seq was performed using ‘Chiifu’ and ‘LP08’ grown under light and dark conditions. MDS and correlation heatmap analyses of RNA-seq samples showed close clustering within each sample group, indicating that RNA-seq libraries were properly generated for ‘Chiifu’ and ‘LP08’ (Supplementary Figures S2A, B). We isolated DEGs between the two parental lines based on comparison of pairwise samples (Supplementary Figure S2C). Under light conditions, 3,054 and 5,222 genes were up- and down-regulated, respectively, in ‘LP08’ compared to ‘Chiifu’. Under dark conditions, 4,659 and 5,219 genes were up- and down-regulated, respectively, in ‘LP08’ relative to ‘Chiifu’ (Figure 2A and Supplementary Table S3).

GO analysis allows for functional annotation by classifying individual genes based on their biological, cellular, and molecular functions. Thus, lists of up- and down-regulated genes in ‘LP08’ compared to ‘Chiifu’ grown under light and dark conditions were assessed with the ShinyGO analysis tool. We set 30% as an arbitrary threshold for significance. Interestingly, GO analysis of up-regulated genes in ‘LP08’ grown under light conditions revealed four categories related to GSL metabolism among the top 10 enriched categories (Supplementary Figure S3A). This result indicates that GSL-related metabolism is affected significantly more in ‘LP08’ compared to ‘Chiifu’. We also obtained the top 10 categories of down-regulated genes in ‘LP08’ compared to ‘Chiifu’. However, we detected no significant enrichment of GO categories above the threshold (Supplementary Figure S3B). These results suggest that GSL metabolism is an actively enhanced metabolic process in ‘LP08’ relative to ‘Chiifu’ grown in the light.

We performed a search using the Basic Local Alignment Search Tool (BLAST) to identify GSL metabolic genes in the B. rapa genome, using sequence information for Arabidopsis GSL pathway genes collected from The Arabidopsis Information Resource (TAIR) website. In total, we collected 162 GSL pathway genes (19 TFs and 143 metabolic pathway genes) from the BRAD website (Supplementary Table S2). Among 162 B. rapa GSL pathway genes, 79 were DEGs between the ‘Chiifu’ and ‘LP08’ lines under light conditions (Supplementary Table S4). Among 79 GSL-related genes, 55 and 24 were up- and down-regulated, respectively, in ‘LP08’ compared to ‘Chiifu’ (Figure 2B and Supplementary Table S4). Interestingly, a majority of the up-regulated genes in ‘LP08’ (49 of 55 genes; 89%) were involved in the aliphatic GSL biosynthetic pathway (Figure 2C and Supplementary Table S4). This finding is in an agreement with the result showing significantly higher amounts of aliphatic GSLs in ‘LP08’ than ‘Chiifu’ (Figures 1A, B). Among the 24 down-regulated genes in ‘LP08’ compared to ‘Chiifu’, both aliphatic and indolic GSL pathway genes were present in almost equal numbers (Figure 2C).

A subgroup of R2R3-type MYB TFs was reported to regulate GSL pathway genes. In total, 14 B. rapa MYB homologs (3 BrMYB28 homologs, 1 BrMYB29, 3 BrMYB34, 3 BrMYB51, 2 BrMYB118, and 2 BrMYB122) were identified in BRAD (Supplementary Table S2). We examined the expression profiles of 14 BrMYB TFs between ‘Chiifu’ and ‘LP08’ grown under light and dark conditions using IGV. First, we analyzed six BrMYBs (BrMYB28.1–3, BrMYB29, and BrMYB118.1–2) involved in the aliphatic GSL pathway (Figure 3A). Expression levels of BrMYB28.1 and BrMYB28.3 were significantly higher in ‘LP08’ than ‘Chiifu’, whereas transcript levels of BrMYB28.2 were similar between the two lines. Transcripts of BrMYB29 and the two BrMYB118 genes (BrMYB118.1–2) were not detected in young seedling plants grown under light or dark conditions. To validate the RNA-seq data, we conducted qRT-PCR analysis of aliphatic GSL pathway BrMYB genes (Supplementary Figure S4A). The results were similar to the RNA-seq data, confirming higher expression of BrMYB28s genes in ‘LP08’ than ‘Chiifu’. Therefore, higher expression of BrMYB28.1 and BrMYB28.3 in ‘LP08’ likely contributes to the greater accumulation of aliphatic GSLs seen in ‘LP08’ than ‘Chiifu’.

Among indolic pathway BrMYB TFs, the expression levels of eight BrMYB TFs (BrMYB34.1–3, BrMYB51.1–3, and BrMYB122.1–2) were compared between ‘Chiifu’ and ‘LP08’ under light and dark conditions (Figure 3B). Among these genes, the expression of BrMYB34.1, BrMYB51.1, and BrMYB51.3 was strongly enhanced under dark conditions in both ‘Chiifu’ and ‘LP08’, whereas BrMYB34.2 and BrMYB34.3 were highly expressed under light conditions in the ‘Chiifu’ line but not the ‘LP08’ line (Figure 3B). While BrMYB51.2 was expressed consistently regardless of genotype or light conditions, the two BrMYB122 genes (BrMYB122.1 and BrMYB122.2) and BrMYB34.1 were expressed only in young seedlings under both light and dark conditions (Figure 3B). This suggested that individual BrMYB genes related to the indolic GSL pathway might regulate indolic GSL biosynthesis under certain environmental conditions or genotypes via diverse processes. Validation of RNA-seq data using qRT-PCR provided similar results (Supplementary Figure S4B). For example, BrMYB34.1, BrMYB34.2, and BrMYB51.3 in the indolic GSL pathway were highly expressed in the dark, while BrMYB34.3 was expressed in ‘Chiifu’ but not in ‘LP08’ (Figure 3B and Supplementary Figure S4B). Taken together, these results show that indolic GSL pathway BrMYB genes exhibit dynamic expression patterns in ‘Chiifu’ and ‘LP08’ under light and dark conditions, whereas aliphatic GSL pathway BrMYB genes are abundantly expressed in the light, particularly in the ‘LP08’ line.

Notably, the expression levels of aliphatic BrMYB genes were substantially higher than those of BrMYB genes in the indolic GSL pathway. The maximum intensity (y-axis of each IGV track view) of BrMYB28.2 and BrMYB28.3 was set to 5,000 and 2,500 read counts, respectively, whereas for BrMYB51.2 and BrMYB34.3 the values were 300 and 200, respectively. This difference indicates that aliphatic GSL pathway BrMYBs had relatively high expression compared to BrMYB genes in the indolic GSL pathway (Figure 3). This finding is in accordance with the level of aliphatic GSLs being substantially higher than the level of indolic GSLs in both the ‘LP08’ and ‘Chiifu’ lines.

In general, genes involved in aliphatic and indolic GSL metabolism have been identified using Arabidopsis as a model plant (Supplementary Table S2). The GSL metabolic process can be grouped into five stages: ‘side-chain elongation’, ‘core structure formation’, ‘secondary modification’, ‘co-substrate’, and ‘breakdown’ (Supplementary Table S2). Among 143 B. rapa GSL metabolic genes (excluding TF genes), 72 (50%) were differentially expressed between the two lines grown in light conditions, with 55 being up-regulated and 24 down-regulated in ‘LP08’ compared to ‘Chiifu’ (Figure 2C and Supplementary Figure S4). A majority of up-regulated genes (46 of 55 in total; 84%) were involved in the aliphatic GSL pathway. Thus, it is plausible that higher expression of upstream BrMYB TF genes (BrMYB28.1 and BrMYB28.3) in ‘LP08’ relative to ‘Chiifu’ results in up-regulation of downstream metabolic genes involved in aliphatic GSL biosynthesis (Figure 4).

For the 24 down-regulated genes in ‘LP08’ (up-regulated in ‘Chiifu’), 11, 11 and 2 genes were related to the aliphatic, indolic, and both pathways, respectively (Figure 5 and Supplementary Table S4). This result indicates that indolic GSL metabolism is more diverse and complicated in these two lines than aliphatic GSL metabolism. We reasoned that the diverse composition of aliphatic and indolic GSLs in ‘Chiifu’ and ‘LP08’ might be derived in part from the dynamic expression of GSL metabolic genes. Taken together, our findings indicate that oilseed-type B. rapa ‘LP08’ had higher expression of BrMYB28 TFs and aliphatic GSL metabolic genes compared to ‘Chiifu’, resulting in greater accumulation of aliphatic GSLs in ‘LP08’. Among aliphatic GSL compounds, GNP was extraordinarily dominant. How the upregulation of BrMYB28s and downstream GSL metabolic genes drives specific accumulation of GNP among aliphatic GSL compounds remains unclear and requires further investigation.

Validation of the RNA-seq results was performed using qRT-PCR for 36 selected genes (24 aliphatic and 12 indolic pathway genes) involved in various stages of GSL metabolism. Among the 24 metabolic genes in the aliphatic GSL pathway, 23 were more abundantly expressed in ‘LP08’ than ‘Chiifu’; the 1 exception was BrGSL-OH.1 (Figure 5). Aliphatic GSL genes involved in the ‘secondary modification’ stage, such as BrFMOGS-OX2, BrFMOGS-OX4, BrAOP2.1, BrAOP2.2, BrST5b.1, BrST5b.4, and BrST5c, were more abundant in ‘LP08’ than ‘Chiifu’. For the GSL-OH.1 gene, the transcript level was much higher in ‘Chiifu’ than ‘LP08’, and was also higher in the dark than light condition. GSL-OH.1 had low expression in ‘LP08’ under both light and dark conditions, in accordance with the RNA-seq results. In addition, between the light and dark conditions, most of the tested aliphatic GSL genes were more dominant in the light than dark, confirming the results of RNA-seq. In total, 12 genes related to the indolic GSL pathway were subjected to qRT-PCR and the results were similar to the RNA-seq data, confirming that transcript patterns of indolic pathway genes differ between ‘Chiifu’ and ‘LP08’, under both light and dark conditions, more than aliphatic GSL genes (Figures 6, 7). Taken together, these results indicate that the oilseed-type B. rapa ‘LP08’ has higher expression of GSL metabolic genes, particularly aliphatic GSL pathway genes, than the vegetable-type ‘Chiifu’ line.

Previously, we developed a mapping population comprised of F₅ 151 RILs through the crossing of ‘Chiifu’ and ‘LP08’ (Kim et al., 2022). Among the 151 F₅ RILs, the seeds of 97 germinated successfully and grew well. These F₅ lines, and the parental lines, were used for extraction of GSLs (Supplementary Table S5). The amounts of individual GSL compounds from individual F₅ and parental lines were applied to linkage mapping based on 8,707 SNPs (Supplementary Figure S7). Although there are no peaks exceeding the LOD for total GSL, aliphatic GSL, and indole GSL, three QTL peaks exceeding the LOD (threshold = 5.906, 5.454, and 6.745, respectively) for GRA, GNP, and GBN were detected in chromosome A03 (Supplementary Figure S7; Figure 8A). these three peaks were corresponding to two peaks in upper arm regions (757 kb–787 kb) and the lower arm regions (approximately 1,950–2,210 kb) of total GSL and total Aliphatic GSL (Supplementary Figure S7). Within those regions, only two GSL pathway genes, BrGSL-OH.1 (Bra022920) and BrMYB28.1 (Bra012961), were located in the upper and lower arm QTL peak regions of chromosome A03, respectively. While BrGSL-OH.1 was found at about 776 kb, BrMYB28.1 was located at about 2,132 kb of chromosome 3. These two genes are both involved in the aliphatic GSL pathway, and no QTL candidate genes related to the indolic GSL pathway were detected in this study.

Two QTL genes (BrMYB28.1 and BrGSL-OH.1) were included in the list of DEGs between the two lines (Figure 2C and Supplementary Table S4). BrMYB28.1, which was located in the highest QTL peak, was more abundant in ‘LP08’ than ‘Chiifu’, whereas BrGSL-OH.1 showed the opposite pattern, with substantially higher expression in ‘Chiifu’ than ‘LP08’ (Figure 6). To elucidate the molecular details underlying the differential expression of these two genes, we cloned the genomic sequences of BrMYB28.1 and BrGSL-OH.1 from both the ‘Chiifu’ and ‘LP08’ lines. The results confirmed that both genes underwent significant genomic “context changes” from the promoter region to the 3′ downstream region between the two inbred lines (Figures 8B, C). For BrGSL-OH.1, several mutations including large deletions and insertions were identified in the promoter region of ‘LP08’ (Figure 8B). In the coding region, two point mutations were identified in the first exon (G to A) and second exon (C to T). In the 3′ downstream region, seven point mutations and two single-base insertions were identified.

Multiple mutations were also detected in the BrMYB28.1 genomic region in ‘LP08’ compared to ‘Chiifu’. Specifically, three point mutations, two deletions, and one 13-base insertion were identified in the promoter region of BrMYB28.1 of ‘LP08’ (Figure 8C). In the coding region, two point mutations (A to T and A to C) were identified in the first intron and one point mutation (T to A) was found in the second exon (Figure 8C). In the 3′ downstream region, eight point mutations and three deletions were identified in ‘LP08’ compared to ‘Chiifu’. As shown in Figure 3A and Supplementary Figure S4A, the expression of BrMYB28.1 was elevated in ‘LP08’ compared to ‘Chiifu’, suggesting that genomic context changes in ‘LP08’ may enhance its transcriptional activity.