Four broccoli genotypes were used within the research project: (1) Ironman, a widely grown commercial F₁ hybrid broccoli that is homozygous for the MYB28 C2 broccoli allele (B/B). (2) 1199 (tradename Beneforte), an F₁ hybrid that is heterozygous for the MYB28 C2 B. villosa/broccoli allele (V/B). (3) HG Inbred, an inbred genotype that is homozygous for the B. villosa MYB28 C2 allele (V/V). (4) 1086, an elite broccoli F₁ hybrid that is homozygous for the MYB28 C2 B. villosa allele (V/V). Whilst both 1086 and HG Inbred have the same MYB28 C2 genotype (V/V) but the genetic background is different with 1086 exhibiting an elite broccoli agronomic phenotype in contrast to HG Inbred. These four genotypes were grown under normal agronomic conditions in an experimental field plot in Norwich in 2017. Leaves were sampled from three plants of each genotype for genome sequencing and gene expression analyses, and florets were sampled for metabolite analyses. All leaf material was harvested whilst the plants were in a vegetative stage prior to the emergence of inflorescences. Additional leaves were sampled from a single glasshouse-grown B. villosa plant for genomic and gene expression analyses.

DNA was isolated from leaf material using a phenol-chloroform extraction protocol (Kang and Yang, 2004). Genomic DNA libraries were prepped and sequenced using Illumina HiSeq paired-end sequencing with 30X coverage by Novogene©¹ producing an average of 125 million total reads per sample with 35X average depth coverage (Supplementary Table 1). The total length of assemblies varied between 396,236,782 and 427,100,513 bp, with the number of contigs ranging from 108,592 to 195,078 (Supplementary Table 2). Raw data from each line were trimmed using Trimmomatic version 0.33 (Bolger et al., 2014) with parameters LEADING:20 TRAILING:20 SLIDINGWINDOW:10:20 MINLEN:50. Reads were assembled into contigs using CLC Assembly Cell² with default parameters.

Contigs from the CLC draft assemblies were anchored using gene models of the reference genome of B. oleracea TO1000 (Assembly BOL, INSDC Assembly GCA_000695525.1 version 97.1) (Parkin et al., 2014). Coding sequences (CDS) of each gene model were aligned to an assembly using NCBI BLASTN (Chen et al., 2015) with default parameters. For each gene model, the contig with the highest alignment score was selected and the genomic position of the gene model was assigned to the contig. Plots were visualised and generated using RStudio v.1.0.143 (Racine, 2012).

K-mers were counted in raw-data of B. villosa (donor) and Ironman (background) using jellyfish³ version 2.1.4. (Marçais and Kingsford, 2011). The parameter “-C” was used for jellyfish to handle a k-mer along with its reverse complement as one item. Sequence data from the draft assembly of the introgression line and k-mer presence or absence within the donor and background were used to determine the ratio of k-mers that are shared between donor and introgression line vs. the ratio of k-mer that is shared between background and introgression line. A custom java programme calculated the ratio of k-mers between the “reference” high glucoraphanin assembly and the “wildtype” (Ironman) and “donor” (B. villosa) per contig. For each contig that was able to be anchored to the TO1000 B. oleracea assembly (Parkin et al., 2014), this ratio was visualised. For each contig of a CLC assembly, the set of k-mers (k = 31) was recorded. Subsequently, it was determined how many of those k-mers were present only in the raw data of the “donor” B. villosa and how many were present only in the raw data of the commercial broccoli background, Ironman. The resulting plot displays the ratio of k-mer alignment between these lines in which k-mer alignment with the B. villosa appears in black.

Transcriptome analysis was undertaken on vegetative leaves from field-grown individual plants of field-grown Ironman and HG Inbred in addition to a glasshouse-grown B. villosa leaf of similar developmental stage. In Arabidopsis, leaves are the site of the biosynthesis of aliphatic glucosinolate that is subsequently transported to roots and reproductive tissues (Andersen and Halkier, 2014; Madsen et al., 2014; Jorgensen et al., 2015). Likewise, in broccoli, leaf expression of MYB28 correlates with the expression of glucosinolates in florets as opposed to the expression of MYB28 in florets themselves (Traka et al., 2013). RNA extraction was performed using E.Z.N.A.^® Plant RNA Kit provided by Omega Bio-tek Inc. A total of seven individual plant RNA samples (three Ironman, three HG Inbred, and one B. villosa) were used to generate TruSeq non-directional RNA libraries at Earlham Institute (United Kingdom). Libraries were sequenced in one pool of seven (7-plex) and run on two lanes of the Illumina HiSeq2500 with a 125 bp paired-end read metric that generated an average of 160,170 million reads per sample. Sequencing was performed by the Earlham Institute. Individual plant sample RNA was used to generate 7 TruSeq non-directional RNA libraries. These libraries were sequenced through pooling in one pool of 7 (7-plex) and run on 2 lanes of the Illumina HiSeq2500 with a 125 bp paired-end read metric with an average of 160–170 million reads per sample. This data are available on NCBI under Project ID PRJNA623495.

Data analysis of RNA-seq raw data was conducted following the protocol for the Tuxedo suite for short reads with default parameters (HISAT2, v.2.0.4, StringTie, v.1.2.2, Ballgown, v.2.8.4) (Pertea et al., 2016). Differential expression analysis was performed using the Ballgown package in RStudio (Racine, 2012; R Core Team, 2013). An enrichment analysis was conducted on the list of differentially expressed (p < 0.05) genes that was generated using the Ballgown software included in the “Tuxedo Suite” for RNA-seq data analysis (R Core Team, 2013), which provided the significantly differentially expressed genes in the HG Inbred MYB28 V/V broccoli when compared to standard broccoli cultivar Ironman MYB28 B/B (p < 0.05). This data set was processed using TopGo analysis to gain a list of enriched gene ontologies in the dataset (Alexa and Rahnenführer, 2009). The gene set analysis statistically compared the representation of Gene Ontology (GO) terms, which were first assigned to the B. oleracea gene models using the Arabidopsis thaliana org.At.tair.db library (version 3.2.3) (Carlson, 2016), in this gene set to that of the “expected” value to determine those to be considered “over-represented” using the BP parameter with a node size of 10, then confirming the statistical significance of this using Fisher’s exact test. This was performed in RStudio v.1.0.143 (Racine, 2012). A Multidimensional Scaling (MDS) of Fragments Per Kilobase Million (FPKMs) displaying Euclidean distances between these samples was analysed as a “principal component” analysis to address the comparison of clustering between replicates of the different genotypes. This was performed using gene expression data, following the removal of low abundance transcripts. Independent comparative analyses included t-tests of Transcripts per Kilobase Million (TPM) in the RNA TPM in the RNA-seq analysis when comparing expression in the standard broccoli with the HG Inbred (Figure 3).

Gene expression analyses of three genes, APS3 (Bo5g021810), MAM C7 (Bo7g098000), and CYP79F1/F2 (Bo5g021810), were additionally quantified by quantitative reverse transcription PCR (RT-qPCR). MAM C7 has previously been referred to as MAM3 (Yin et al., 2017). Primer sequences can be found in Supplementary Table 3. RT-qPCR amplification was from the RNA extractions that had generated cDNA using the REVERSE TRANSCRIPTASE M-MLV Kit Supplied by Life Technologies Ltd. Reactions included 5 μl RNA, 1 μl oligo dT, 1 μl deoxynucleotide 5′-triphosphates (dNTPs), and 5 μl distilled H2O and were run at 65°C for 5 min. This was followed by the addition of 4 μl of 5 × buffer and 2 μl 0.1 M Dithiothreitol (DTT) before being kept at 37°C for 2 min. Finally, 1 μl Moloney Murine Leukemia Virus (MMLV) reverse transcriptase was added to the mixture before being incubated at 37°C for 50 min followed by 70°C for 15 min. Concentrations of 150 ng/μl of cDNA were used for RT-qPCR reactions. Gene expression was quantified using the QuantiNova SYBR Green PCR Kit from Qiagen. PCRs were carried out in a Bio-Rad CFX96 machine (C1000 Touch). The PCR cycling conditions were 95°C for 15 min, 40 cycles of 95°C for 15 s, and 60°C for 60 s.

Sulphate and S-methylcysteine sulfoxide (SMCSO) analyses were performed as described previously with minor modifications (Koprivova et al., 2008; Kubec and Dadáková, 2009). Glucosinolate analysis was performed as previously described (Saha et al., 2012). Freeze-dried inflorescences resembling commercial edible broccoli floret samples were sent to the Eurofins Food and Feed Testing Laboratories for quantification of total sulphur, cysteine, and methionine by high-performance liquid chromatography (HPLC) analysis. Phenotypic analyses of metabolites included a one-way ANOVA comparing metabolite content of each of the broccoli cultivars for a single year along with Tukey’s multiple comparison tests. All statistical analyses were carried out on GraphPad Prism (version 8.2.0).

Field-grown broccoli florets of the four genotypes described above were analysed for their glucosinolate phenotype. Broccoli genotypes that had either one or two MYB28 V alleles had significantly greater total content of the aliphatic methionine-derived glucosinolates glucoraphanin and glucoiberin in their florets (Figure 1A and Supplementary Table 4), but there was no difference in glucosinolate content between genotypes that were homozygous (HG Inbred and 1086) or heterozygous (1199) for the MYB28 C2 V allele. There were no significant differences in tryptophan-derived glucosinolates or the other major sulphur-containing broccoli metabolites, S-methylcysteine sulfoxide, methionine, cysteine, and sulphate (Figure 1B).

Three broccoli genotypes (Ironman, HG Inbred, and 1086) and one accession of B. villosa were sequenced to ≈35-fold coverage using Illumina PE sequencing. Two of the broccoli genotypes, 1086 and HG Inbred, contained introgressions from B. villosa and had high levels of glucoraphanin. The third genotype, the commercial broccoli Ironman, lacked any introgression from B. villosa and was used as a partial surrogate for the broccoli genetic background of the genotypes with B. villosa introgressions.

Draft gene-space assemblies of 1086 and HG Inbred were generated. Gene-containing contigs were anchored to the B. oleracea reference assembly of TO1000 (Parkin et al., 2014). For each anchored contig of an assembly, its set of k-mers (k = 31) was determined. We additionally obtained all k-mers from the raw data of Ironman and B. villosa (Table 1). Testing for each contig of an introgression genotypes assembly involved quantifying how many of its contigs k-mers were unique to B. villosa and how many were unique to Ironman. This enabled to ascertain which contigs represented sequence from introgressed B. villosa segments. The anchored positions for contigs that likely resembled a signal for true were identified (Figure 2). As was expected, the most pronounced region of introgression was towards the end of C2 where MYB28 has previously been positioned (Traka et al., 2013), but there was evidence of additional introgression across all other linkage groups. The distribution of introgressed segments of B. villosa genome was similar in the elite F₁ hybrid 1086 to the HG Inbred (Figure 2 and Supplementary Figure 1). It is notable that despite the widespread introgression of B. villosa into 1086, this genotype has an elite agronomic phenotype with no apparent negative traits.

The overall percentage alignment of RNAseq reads from the three samples, each of Ironman and HG Inbred to the TO1000 B. oleracea reference genome (Parkin et al., 2014), was 90.0 ± 0.31 and 90.7 ± 0.48, respectively. The alignment of B. villosa and TO1000 was significantly lower at 82.3% (p < 0.001), indicative of its greater divergence from TO1000 when compared to broccoli. 3286 genes were differentially expressed between Ironman and HG Inbred (uncorrected p < 0.05). Of the three MYB28 paralogues on C2, C7, and C9, only the MYB28 C2 allele was differentially expressed between Ironman and HG Inbred, with levels of expression in HG Inbred similar to that in B. villosa (Figure 3A). Examples of relative expression of genes involved in sulphate reduction, methionine metabolism, and glucosinolate biosynthesis are shown in Figures 3B–D, respectively. RT-qPCR independently confirmed the enhanced expression of three genes—APS3, MAM C7, and CYP79F1/2—in HG Inbred when compared to Ironman. Gene ontology analyses confirmed the significant difference of several metabolic pathways associated with sulphur metabolism and glucosinolate biosynthesis (Table 2).

The enhanced level of aliphatic methionine-derived glucosinolates in the high glucoraphanin genotypes is likely to be due to the enhanced expression of the MYB28 C2 allele (Figure 3A). To provide further insight as to the possible cause of the enhanced expression, the CDS and 5′ upstream regions of MYB28 C2 of B. villosa, HG Inbred, and 1086 were compared with that of the commercial broccoli, Ironman. Only one single nucleotide polymorphism (SNP) was found in the sequence coding of the R2R3 DNA-binding domain (Figure 4). This SNP is synonymous and was previously used to confirm the introgression of the B. villosa allele into the high glucoraphanin broccoli genotypes (Traka et al., 2013; Supplementary Figure 2). The most significant difference between the sequences was a polymorphism in an AT microsatellite or short tandem repeat (STR) upstream of the ATG start codon in which the length of the A/T sequence was longer in B. villosa, HG Inbred and 1086 compared to Ironman in which it was largely absent (Figure 4). To further explore this variable region, RNA-Seq data from B. villosa and the HG Inbred genotypes were mapped against TO1000. It was apparent that the AT-STR (AT STR) was situated within the 5′ untranslated region (5′UTR), and that the significant difference in sequence prevented the mapping of the broccoli and B. villosa sequence data to TO1000. It was also found that the AT STR was absent from the domesticated broccoli B. oleracea reference sequence DH1012, in a similar manner to Ironman (Chen et al., 2021).⁴