Complex Horticultural Quality Traits in Broccoli Are Illuminated by Evaluation of the Immortal BolTBDH Mapping Population

Abstract

Improving horticultural quality in regionally adapted broccoli (Brassica oleracea var. italica) and other B. oleracea crops is challenging due to complex genetic control of traits affecting morphology, development, and yield. Mapping horticultural quality traits to genomic loci is an essential step in these improvement efforts. Understanding the mechanisms underlying horticultural quality enables multi-trait marker-assisted selection for improved, resilient, and regionally adapted B. oleracea germplasm. The publicly-available biparental double-haploid BolTBDH mapping population (Chinese kale × broccoli; N = 175) was evaluated for 25 horticultural traits in six trait classes (architecture, biomass, phenology, leaf morphology, floral morphology, and head quality) by multiple quantitative trait loci mapping using 1,881 genotype-by-sequencing derived single nucleotide polymorphisms. The physical locations of 56 single and 41 epistatic quantitative trait locus (QTL) were identified. Four head quality QTL (OQ_C03@57.0, OQ_C04@33.3, OQ_CC08@25.5, and OQ_C09@49.7) explain a cumulative 81.9% of phenotypic variance in the broccoli heading phenotype, contain the FLOWERING LOCUS C (FLC) homologs Bo9g173400 and Bo9g173370, and exhibit epistatic effects. Three key genomic hotspots associated with pleiotropic control of the broccoli heading phenotype were identified. One phenology hotspot reduces days to flowering by 7.0 days and includes an additional FLC homolog Bo3g024250 that does not exhibit epistatic effects with the three horticultural quality hotspots. Strong candidates for other horticultural traits were identified: BoLMI1 (Bo3g002560) associated with serrated leaf margins and leaf apex shape, BoCCD4 (Bo3g158650) implicated in flower color, and BoAP2 (Bo1g004960) implicated in the hooked sepal horticultural trait. The BolTBDH population provides a framework for B. oleracea improvement by targeting key genomic loci contributing to high horticultural quality broccoli and enabling de novo mapping of currently unexplored traits.

Introduction

Improvement of broccoli and other Brassica oleracea vegetables (cauliflower, cabbage, kale, Gai lan, Brussels sprouts, kohlrabi, and collard) is constrained by complex interactions of many genes affecting plant architecture, developmental processes, and yield. B. oleracea vegetable crop groups have benefited from a number of advances in plant biotechnology, gradually increasing the overall understanding of these quality-based traits. Specifically, diversity and domestication processes (Cheng et al., 2016; Stansell et al., 2018; Yousef et al., 2018; Li et al., 2019) have been clarified, next-generation sequencing and high-quality reference genomes (Liu et al., 2014; Parkin et al., 2014; Golicz et al., 2016) have expedited discovery of molecular markers associated with key traits, and diverse mapping populations segregating for these traits have been characterized (Kianian and Quiros, 1992; Landry et al., 1992; Camargo and Osborn, 1996; Ramsay et al., 1996; Bohuon et al., 1998; Hu et al., 1998; Lan and Paterson, 2000; Sebastian et al., 2000; Lan and Paterson, 2001; Axelsson et al., 2001; Li et al., 2003; Gao et al., 2007; Brown et al., 2014; Lee et al., 2015a; Li et al., 2015). For example, projects integrating these tools such as the Eastern Broccoli Project (SCRI No. 2010-51181-21062) and the USDA Vegetable Brassica Research Project (CRIS No. 6080-21000-019-00D) have developed heat-tolerant broccoli germplasm adapted to novel environments, reducing costs and enabling more sustainable production models (Atallah et al., 2014; Farnham and Björkman, 2011).

A current limitation in B. oleracea vegetable crop improvement is a lack of publicly available mapping populations, constraining information integration across research programs. Attempts to unify existing maps have been limited due to variable germplasm, different marker types, and linkage group nomenclature (Hu et al., 1998). Furthermore, these populations are often difficult to maintain due to self-incompatibility (Farnham, 1998; Bohuon et al., 1998; Sebastian et al., 2000; Pink et al., 2008; Walley et al., 2012).

To address these issues, the double-haploid (DH) BolTBDH population was developed from a cross between morphologically distinct parents (B. oleracea var. alboglabra × B. oleracea var. italica) that segregates for horticultural quality traits specific to broccoli (Iniguez-Luy et al., 2009). BolTBDH offers several distinct advantages over other mapping populations: a large sample size (N∼175), a high degree of self-compatibility, and a short generation time. The rapid-cycling parental taxa ‘TO1000DH3’ (P₁; var. alboglabra) is the reference organism for the B. oleracea v.2 genome (Parkin et al., 2014). ‘Early Big’ (P₂; var. italica) has been evaluated in previous studies (Li et al., 2003; Gao et al., 2005; Tortosa et al., 2018). Moreover, both P₁ and P₂ and included in the B. oleracea pangenome (Golicz et al., 2016). Furthermore, this population has already been used to investigate the genetic control of important traits: glucosinolate content (Sotelo et al., 2014), organ-specific phenylpropanoid metabolism (Francisco et al., 2016), antioxidant content (Sotelo Pérez et al., 2014), and black rot (Xanthomonas campestris pv. campestris) resistance (Iglesias-Bernabé et al., 2019). Under standard greenhouse conditions, BolTBDH lines will typically produce self-seed without the need for hand pollinations. The work presented here increases the value of this population by generating many high-quality genome-wide SNP markers and generating robust phenotypes for 25 horticultural quality traits.

The BolTBDH population provides an unique opportunity to evaluate the genetic basis of the heading broccoli phenotype due to the marked dissimilarity between the parental lines: P₁ is rapid-flowering (∼ 65 days to flowering) and exhibits a leafy, non-head-forming inflorescence, whereas P₂ is relatively late-flowering (∼ 85 days to flowering) and exhibits a heading broccoli phenotype characterized by extensive meristem proliferation during floral bud development and internode elongation (Björkman and Pearson, 1998). While considerable work has investigated head formation under optimal and heat-stressed conditions (Duclos and Björkman, 2008; Farnham and Björkman, 2011; Hasan et al., 2016; Lin et al., 2019), the exact genetic basis of this phenotype remains elusive (Axelsson et al., 2001; Li et al., 2003; Gao et al., 2007; Okazaki et al., 2007; Razi et al., 2008; Matschegewski et al., 2015; Irwin et al., 2016). There is a growing consensus of the central importance of the homologous flowering timing MADS-box transcription factor FLOWERING LOCUS C (BoFLC) in regulating the reproductive transition by inhibiting downstream BoSOC1 and BoFT expression; in turn, delaying a suite of floral-identity genes including BoLFY, BoAP1, and BoCAL (Lin et al., 2018). The genetic basis of the heading broccoli phenotype could be explained by a number of models: simple control by one or several genes, a constrictive-conditional model where multiple genetic factors must be present, or a pleiotropic model, where several key developmental genes would underlie the broccoli heading phenotype and be further modified by additional downstream factors. Under a simple control model, the heading broccoli phenotype would exhibit qualitative control by a limited number of QTL. Under a purely constrictive-conditional model, the broccoli heading phenotype would occur only when a minimum set of independent factors were present. Under a purely pleiotropic model, a small number of developmental loci or genes are implicated in heading quality traits with epistatic interactions with additional loci.

When breeding multiple quality traits in horticultural crops, it is often challenging to determine the degree that individual traits contribute to the overall quality of these crops when predictor traits are correlated. Relative-importance analyses (RIA) allows quantification of the proportional contribution of a predictor variable to the overall quality-model R², considering both unique and joint contributions with other variables (Gromping, 2006) and may be used to establish breeding priorities within a horticultural context (Stansell et al., 2017). Within this study, RIA was used to evaluate the independent contribution of individual traits to overall horticultural quality.

Therefore, our main objectives were to: a) characterize the phenotypic variation of horticulturally important traits within the BolTBDH population; b) produce a reference set of robust and high-quality BolTBDH markers; c) identify optimal QTL models to best explain key horticultural quality traits important to broccoli germplasm; and d) identify which candidate B. oleracea developmental genes collocate with observed QTL.

Materials and Methods

Germplasm

The BolTBDH population was generated via anther culture from the parental lines ‘TO1000DH3’ DH (P₁; B. oleracea var. alboglabra) × ‘Early Big’ DH (P₁; B. oleracea var. italica) (Iniguez-Luy et al., 2009). Seed was provided by the USDA Vegetable Laboratory in Charleston, SC. All initial lines (N = 202; P₁ and P₁ inclusive) were increased in 2016 and closed-bud pollinations were made to verify selfing integrity. Except for P₁ in Y₁ due to inadequate seed quality, all lines were sown into 128 cell trays May 11 in Y₁ and May 9 in Y₂. Seedlings were grown in a greenhouse and transplanted into Lima silt loam fields in Geneva, NY on May 28–29 in Y₁ and June 8–11 in Y₂. All lines were divided into four randomized replications and transplanted onto raised beds with each plot containing 10–12 plants per genotype. Drip irrigation was applied as needed and any additional cultural practices were as previously described (Farnham and Björkman, 2011). Hourly weather data was collected locally at Cornell AgriTech (Figure S1).

Traits Investigated

Plots were evaluated daily and deemed mature when 1/3 of the plants reached a heading or heading-equivalent stage. Traits within six classes considered important to broccoli or other B. oleracea crop groups were chosen: architecture, biomass, bud morphology, leaf morphology, head quality, and phenology (Table 1; Figure 1; Stansell et al. (2017)). LT was measured as the degree of lateral shoot growth. MH was evaluated as flower bud bunching before antithesis. MS was measured as above-ground biomass of a representative central plant. VG was evaluated as overall plant vigor. Leaf color/waxiness (LC) was evaluated visually. LA and LM were evaluated as the leaf-tip angle and degree of leaf margin serration. No intermediate flower color was detected so FC was scored as a binary trait. Other bud morphology traits (SE, SS, SF, ST, and SH) were evaluated visually as unopened buds at head maturity. The traits bud size (BS), bud uniformity (BU), bracting (BR), head compactness (HC), head diameter (HD), head uniformity (HU), head extension (HE), head shape (HS), overall-heading quality, and (OQ) were evaluated following standardized protocols developed by the Eastern Broccoli Project using an ordinal scale (1 = worst; 5 = best) with slight population specific modifications (e.g.: adjusting scale centering to account for smaller DH heads) (Stansell et al., 2017). Days to maturity and flowering (DM and DF) were calculated as days from sowing to head maturity and first flowering respectively. Holding ability HA was defined as DF–DM. Correlation matrices were computed between traits as well as between trial years by invoking the Spearman method with the cor() function in R v3.6.0(R-Core-Team, 2018). RIA of overall heading quality OQ was conducted with the R package ratervar (Stansell et al., 2017) using 1,000 bootstraps under the metric “lmg” by fitting the model:

Figure 1

Marker Development, Map Construction

DNA was extracted from young leaf tissue at the 2–3 true leaf stage, bulked from five plants and extracted according to standard protocols. GBS was accomplished at the University of Wisconsin Biotechnology Center DNA Sequencing Facility following methods of Elshire et al. (2011). Library construction occurred in 96-well plates with ApeKI digestion followed by sequencing on Illumina HiSeq 2500, producing 100-bp single-end reads. SNP production was accomplished using the TASSEL v5.2.35 GBS pipeline (Glaubitz et al., 2014): 214,757,912 raw sequence reads were initially generated and 168,722,056 (78.6%) were good barcoded reads (minimum quality score > 15; minimum K-mer count = 10; min K-mer length = 10) and reads were collapsed into 11,842,938 tags. Alignment of filtered tags were was accomplished with default settings with BWA v0.7.15 (Li, 2013) to the B. oleracea genome v2.1 (Parkin et al., 2014), producing 670,347 mapped tags. Initial filtering removed indels, loci with more than 10% missing data and minor SNP states. Missing data were imputed using the FSFHap plugin (Swarts et al., 2014) invoking the cluster algorithm option. Nucleotide data was recoded as ABH genotypes and heterozygous, missing, or ambiguous calls were removed. Additional quality control steps using the package rqtl (v.1.44-9; Broman et al. (2003)) were performed: taxa with > 5% missing data were removed and taxa pairs exhibiting over 95% pairwise genetic similarity were pruned. Individuals exhibiting three times more crossover events above standard deviation were removed. Ultimately, 175 DH lines were included in the final datasets. Markers exhibiting identical segregation patterns were pruned. Markers with χ² −log[p.adj] > 40 segregation were removed as likely genotyping errors. Linkage disequilibrium analysis identified 87 markers that appeared to be assigned to the wrong linkage groups and these were removed. Markers were imputed using a Viterbi algorithm and genetic maps were constructed using the Kosambi mapping function (Supplementary Data S1–S6).

MQM Mapping

Multiple QTL mapping (MQM) was accomplished using rqtl using the forward and backward search algorithm stepwiseqtl() by searching for QTL models with the highest penalized LOD score (Broman et al., 2003; Arends et al., 2010). Genotypes were first simulated with hidden Markov modeling (N = 1,000) followed by estimation of the true underlying genotype probabilities calculated across a 1 cm fixed stepwidths. Initially, 1,000 permutations of two-dimensional scans per trait were run to establish trait-specific genome-wide significance thresholds and to calculate MQM model penalties. An initial forward scan was used to determine the maximum number of QTL per trait to include in stepwiseqtl(). The normal model was invoked for stepwiseqtl() except for the trait FC which was run under a binary model. Between model selection steps in MQM analysis, QTL positions were refined using iterative maximum likelihood scanning. Additional non-parametric scans using extension of the Kruskal-Wallis tests were run as an additional confirmation step to account for non-normal distributions. Epistatic effects were identified using two-dimensional scans with the scantwo() function and chosen according to the FV₁ model: the log₁₀ likelihood of the full QTL model on chromosomes j and k was compared against a single model QTL on j or k (p_fv1 < 0.05). QTL intervals were called by applying the function find_peaks() by choosing MQM peaks surpassing the genome-wide significance threshold (α = 0.05), estimating interval start and stop locations as 95% Bayes credible intervals (BCIs), and selecting the best adjacent markers. Percent phenotypic variance explained was calculated by

Candidates were determined by subsetting the B. oleracea v2.1 genome annotation by the 95% BCI followed by further extraction of A. thaliana BLASTP hits. Additionally, a review of candidates implicated in developmental control of the horticultural traits under consideration was followed by pan-taxonomic searches using the online tools TAIR (Berardini et al., 2015) and EnsemblePlants (Ruffier et al., 2017) to determine the physical locations of 391 B. oleracea homologs (Supplementary Data S7). These candidates were then cross-referenced against the identified QTL using a custom R script that calculated the difference in physical location of a candidate and LOD peak location, retaining only candidates located within the 95% BCI or within 1 Mbp from the LOD peak.

Results

Phenotyping

For architectural traits, suppression of lateral shoot growth (LT) was considerably higher in Y₁ (Y₁ = 3.1, Y₂ = 1.8; Figure 2A; Table 2) and was moderately correlated between years (r = +0.61; Figure 2B). Presence of floret bunching (MH) was lower in the Y₂ trial (0.98 vs 0.86; r = +0.22). Average above-ground biomass (MS) was 90.5 g lower in Y₂ trials, although plant vigor (VG) was consistent between years. MS and VG were positively correlated (r = +0.57). The leaf morphology traits LA and LM exhibited moderate between-trait correlation (r = +0.45). Leaf color (LC) was moderately correlated between years (r = +0.52). Except for flower color (FC; r = +1.00), other bud morphology traits exhibited low to intermediate year-to-year correlation, ranging from “bud profile:fig” (SF; r = +0.15) to flower bud size (BS; r = +0.55).

Figure 2

All head quality traits were lower in Y₂ when compared to Y₁ except for bud uniformity (BU; −0.03) and head-quality traits were strongly correlated between years (r = +0.82). Head shape (HS; −0.55) and head compactness (HC; −0.47) exhibited the largest changes from Y₁ to Y₂ trials. Overall horticultural quality (OQ) was positively correlated with all other head quality traits such as head uniformity (HU; r = +0.94) and head shape (HS; r = 0.88) except for a negative correlation with head extension HE (r = −0.51). Days to head maturity DM (r = +0.93) and flowering DF (r = +0.88) were correlated between years, although lower in Y₂ trials by 19.7 d and 19.4 d. DM and DF were correlated with each other (r = +0.94). Time from maturity to flowering (HA) was consistent between years (Y₁ = 5.34 d and Y₂ = 5.61 d) but was not strongly correlated with DM (r = +0.13) or DF (r = +0.39).

Relative importance analyses indicated that variation in the traits HU (24.7%), HC (21.3%), HS (16.4%), and BR (15.2%) explained 78.0% of the variability in overall heading quality model (Figure 2C, Table S1). Although OQ was correlated with DM (r = +0.53) and DF (r = 0.46), these traits were not strong predictors of OQ, each explaining < 3% of the variance in the RIA quality model.

Genotyping and Mapping

Genotype-by-sequencing of all initial lines (N = 202) resulted in 168,722,056 quality barcoded reads distributed across 2,529,429 unique tags of which 670,347 were mapped, producing 263,998 SNPs. FSFhap imputation and filtering for minor allele frequencies, missing data, and minor SNPs states reduced this value to 15,774; decreasing percent missing data from 21.0% to 2.7% (Figure S2). Markers assigned to the wrong linkage groups were removed (Figure 3A) and 1881 high-quality, non-duplicated markers were selected and distributed with a mean coverage of 4.25 SNPs/Mbp (Table S2).

Figure 3

Lines identified as hybrids indicated by heterozygous calls, twinned lines, and crossover event outliers were removed, resulting in 175 lines included in multiple QTL mapping (Supplementary Data S1–S3). Deviation from the expected 1:1 segregation pattern (FDR < 0.05) was observed in 61.8% of markers and 50.9% alleles were contributed from P₁. Several chromosomes exhibited strong segregation distortion, notably C03: 73.3% (P₁), C08: 67.6% (P₂), and C07: 63.3% (P₁) (Figure 3B). Although segregation distortion may reduce overall QTL detection power, it does not limit detection in sufficiently dense marker sets, therefore these markers were retained.

A genetic map was constructed from 1,881 markers spanning all linkage groups (Figures 3A–C), with markers per chromosome ranging from 106 (C01) to 283 (C03). Total map distance was 1,060.8 cm with a mean and maximum marker spacing of 0.57 and 16.67 cm/marker. Maximum marker spacing per chromosome ranged from 2.26 cm/marker (C08) to 16.67 cm/marker (C05). Mean crossover events per double haploid were 15.0 ± 9.7; max = 61, min = 4.

MQM Mapping

MQM identified 56 single (Figure 4A; Table 3) and 41 epistatic (Table 5; Figure 5) QTL. QTL per chromosome ranged from 1 (C02) to 12 (C09)(Figure 4B). LOD values ranged from 2.85 (MS_C05@39.5; PVE = 7.3) to 39.9 (FC_C03@55.7; PVE = 65.8). Bayesian confidence intervals (95% CI) ranged from 0.3 to 59.9 Mbp (mean = 16.2 Mbp), and contained on average 1902.2 coding sequences and 330.6 A. thaliana hits (Table 4). QTL per trait class ranged from 3 (biomass) to 20 (head quality). Three architecture trait QTL were identified for LT, and no MH QTL were identified. Four biomass QTL were identified, one MS and three VG. Nine leaf morphology QTL were identified: five LA, three LM, and one LC. Thirteen significant bud morphology QTL were identified: one FC, SE, SF, SH, SS, and ST, five BS, and two BU. Twenty head quality QTL were identified, including two BR, two HC, three HD, three HE, two HS, three HU, and four OQ. Seven phenology QTL were identified: two DM, two DF, and three HA. Epistatic QTL (p < 0.05) were detected in every trait class except FC, MH, and SS.

Figure 4

Figure 5

Optimal MQM models were determined for all traits ranging in complexity from FC: FC ∼ 3@55.7 to LA ∼ 1@3.4 + 3@0.7 + 6@18.7 + 7@37.0 + [3@0.0 × 9@48.5] + [3@0.0 × 7@36.6] + [7@36.9 x 9@24.4] + [1@2.6 x 9@49.5] + [6@18.7 × 7@36.7] + [1@2.6 × 3@0.0] + [2@3.1 × 8.5] (Table 5).

Ten QTL hotspots appearing to harbor QTL across multiple traits were identified: two biomass related (Bio₅ and Bio₇), four morphology related (Bud₁, Lea₃, Bud₄, and Lea₇), three heading quality related (HQ₄, HQ₈, and HQ₉), and one phenology related (Phe₃) (Table 4).

Discussion

Phenotyping

Field evaluations were conducted in growing seasons that differed somewhat in temperature stress (Figure S1). In year two (Y₂), a strong heat wave in the first half of July (mean high = 29.5°C) coincided with the reproductive transition. Horticultural quality scores were lower in Y₂ trails, likely due to less optimal temperatures during the transition to flowering (Table 2) [e.g.: suppression of lateral shoot growth (LT) was evaluated 1.3 points lower in Y₂ trials].

Genotyping and Mapping

The GBS markers (Figure 3A and Figure S2) generated were high quality with less than 5% missing data and distributed across all chromosomes with a mean coverage 1.77 markers/cM (Figure 3C; Table S2; Figure S3), a six-fold improvement in marker coverage from previous BolTBDH maps, which relied upon approximately 300 SSR and RFLP markers (Sotelo Pérez et al., 2014; Francisco et al., 2016). Segregation distortion was prevalent but chromosome-specific (C03, C07, C08; Figure 3B).

MQM Analysis

Architecture

Excessive lateral side-shoot growth is a horticultural defect that increases harvest costs and reduces yield. Although improved broccoli F₁ hybrids typically exhibit strong apical dominance, variability in lateral side-shoot growth occurs among genotypes, typically under environmentally stressed conditions (Potters et al., 2007). In a bulk-segregant analysis of shoot branching in B. juncea, Muntha et al. (2018) identified BjPAT1 and its signal integrator BjBRC1 as branching candidates. The BoPAT1 or BoBRC1 homologs did not collocate with the LT BCI identified in BolTBDH. Tyagi et al. (2019) determined that mutations in BjSOC1 may influence degree of lateral branching. A BoSOC1 ortholog (Bo3g038880) was identified within the 95% LT_C03@5.9 BCI, although this candidate was located >9.5 Mbp from the LOD peak, and is not likely involved in the LT phenotype response observed in BolTBDH. He et al. (2017) conducted GWA and QTL mapping of lateral branching in B. napus and identified BnaC03g63480D as a branching number candidate. The best BnaC03g63480D ortholog, Bo3g159770, was nominally located within the LT_C03@5.9 BCI, but is not considered a likely candidate within this population. Additional candidates associated with meristem identity and fate identified within LT QTL are listed (Table 6).

Biomass

In broccoli, higher yielding genotypes are preferred for commercial production and vegetative biomass is positively correlated with head biomass (Lin et al., 2013). In a F₂ broccoli × broccoli population, Lin et al. (2013) identified biomass QTL on linkage groups C1, C5, C8, and C9, although the authors did not identify candidates associated with these loci. In a GWA study of seedling vigor in B. napus, the vigor candidates Bna.SCO1, Bna.ARR4, and Bna.ATE1 were identified by Hatzig et al. (2015), although no homologous candidates were located within the BolTBDH VG BCIs. A GWA and transcriptome analysis in B. napus by Lu et al. (2017) identified two yield-related candidates: BnaA05g29680D and BnaC04g42030D. These candidates were not identified in BolTBDH biomass QTL, although one ortholog of BnaA05g29680D (Bo5g139830) was identified adjacent to the Bio₅ hotspot.

The homologous candidates BRC1 (Bo1g117490 and Bo5g117410) is closely related to TEOSINTE BRANCHED1, and is a putative transcription factor involved in arresting axillary bud development and limiting axillary bud growth (Muntha et al., 2018) and was identified within Bio₅. Additional homologous VG candidates involved in growth and growth regulation were identified (Table 6).

Leaf Morphology

Variation in leaf morphology is useful to improve and develop novel market classes of B. oleracea leafy greens. Lan and Paterson (2000) identified robust QTL associated with leaf lamina width on C01 and C07, although they did not identify likely candidates. Previous studies have identified leaf-apex QTL on LGO1 and LGO3 (cauliflower × Brussels sprouts) (Sebastian et al., 2002) and C06 and C07 (broccoli × broccoli) (Walley et al., 2012). Li et al. (2009) identified the candidate BrAS1 involved in leaf lamina width. The candidates AtLUG, AtWOX1, and AtAN3 were shown by Zhang et al. (2019) to be involved in leaf blade outgrowth although none of these candidates were identified in LA BCI within BolTBDH. In Arabidopsis, gif1 mutants exhibit a longer, narrow leaf phenotype (Shimano et al., 2018) and a GIF1 homolog (Bo7g093130) collocated within LA_C07@36.6 QTL.

In a QTL-seq analysis of ornamental kale, Ren et al. (2019) identified a lobed-leaf candidate BoLl to C09 (38.82–40.12 Mb), although BoLl did not collocate with BolTBDH LM QTL. Ni et al. (2017) transformed BnLMI1 into Arabidopsis producing serrated leaf margins, similar in appearance to the serrated leaf margin phenotype observed in BolTBDH. In BolTBDH, a BoLMI1 ortholog (Bo3g002560) collocated closely with Lea₃ in agreement with Ni et al. (2017).

Alterations in cuticular wax alters herbivore behavior and may confer broad resistance (Branham and Farnham, 2017). Leaf color in B. oleracea is strongly affected by cuticle wax content and is likely responsible for the variation in blue-green matte and dark-green glossy leaf appearance observed in the BolTBDH population. Xu et al. (2019) conducted fine mapping of the cuticular wax synthesis gene BoWax1 controlling the glossy trait in a F₂ cabbage population although this candidate was not identified in BolTBDH LC BCI. In another cabbage population segregating for the glossy leaves, Zhu et al. (2019) mapped a non-wax glossy NWGL locus to a 99 kb interval in C08. Neither BoWax1 or the NWGL locus collocated with the LC QTL identified in BolTBDH. Lee et al. (2015b) analyzed expression of wax synthesis candidates and determined that the homologous candidates LACS1, KCS1, KCR1, ECR, CER3, and MAH1 were differentially expressed in broccoli lines with elevated cuticular wax levels. Of these candidates, only MAH1 homologs collocated with LC_C09@15.1 and the BCI included five MAH1 copies: Bo9g053360, Bo9g053340, Bo9g053260, Bo9g053220, Bo9g053170. MAH1 encodes CYP96A15, a midchain alkane hydroxylase, involved in cuticular wax biosynthesis (Greer et al., 2007).

Bud Morphology

B. oleracea flower petals are typically white, but a dominant mutation of BoCCD4 implicated in a yellow-flower phenotype via inactivation of a carotenoid-degrading enzyme has been previously described (Zhang et al., 2015; Han et al., 2019). In BolTBDH, a single flower color QTL was identified (FC_C03@55.7) and BoCCD4 (Bo3g158650; C03:56.61Mb) was located < 1 Mbp from this LOD peak, in agreement with (Zhang et al., 2015). Interestingly, 77.2% of BolTBDH lines exhibited a white flower phenotype, an unexpected result given single locus control in the MQM trait model FC ∼ 3@55.7, however the segregation distortion observed at this locus (f = 0.76 P₁) is consistent with this result (Figure 3B).

Certain bud morphology defects contribute to a reduction in head quality, often rendering broccoli heads unmarketable. Hooked sepals are commonly encountered in broccoli hybrids, resulting in a non-uniform crown surface and a reduced ability to shed water. In B. rapa, Zhang et al. (2018) identified loss-of-function BrAP2 alleles with sepal defects similar to the hooked sepal phenotype observed in BolTBDH. In BolTBDH, the optimal MQM model for hooked sepals (SH ∼ 1@2.1 + [1@0.2 × 1@20.2]) co-localized with the bud morphology hotspot Bud₁. Bud₁ QTL include the homologous candidates AP2 (Bo1g004960), in agreement with Zhang et al. (2018). Further evaluation of AP2 (Bo1g004960) homologs and other floral developmental candidates (Table 6) in the Bud₁ hotspot may prove useful for improvement of bud morphology traits.

Small flower buds are preferred in broccoli, but heat-tolerant germplasm typically exhibits larger flower buds. MQM modeling of bead size in BolTHDH resulted in a complex model trait model: BS ∼ Bud₁ + 3@36.8 + Bio₅ + 6@3.8 + [Bud₁ × Bio₅] + [Bio₅ × 6@3.9], suggesting complex genetic control of this trait.

Unequal-sized flower buds are a common horticultural defect in broccoli and bud uniformity requires an arrest of enlargement of older buds until younger buds reach an equivalent size, requiring complex coordination (Roeder et al., 2010). Lin et al. (2018) identified a QTL within CO6 (qCQ-6) associated with a reduction in uneven-sized flower buds and identified PAN and the most probable candidate within this interval, exhibiting strongly differential expression in floral bud at harvest stage. In the BolTBDH model of flower bud uniformity (BU ∼ Lea₃ + 4@51.5), the PAN homolog Bo6g107140 was not harbored within the identified BU QTL.

Head Quality

For heading quality traits, the hotspot HQ₉ was remarkably pronounced within MQM analysis of BolTBDH. HQ₉ is syntenic with the telomeric region of the short arm of Arabidopsis Chr5 (O’neill and Bancroft, 2000), and this region has previously been shown to carry the homologs of the key flowering-time and vernalization-response genes TLF2, COL1, CO, and FLC (Osborn et al., 1997; Lagercrantz et al., 2002; Lin et al., 2005; Okazaki et al., 2007; Uptmoor et al., 2008; Iniguez-Luy et al., 2009; Hasan et al., 2016; Shea et al., 2018). Razi et al. (2008) mapped BoFLC1 to the end of C09, and homology searches in the BOLv.2 genome indicated that two FLC copies (Bo9g173370 and Bo9g173400) appear to be harbored within the HQ₉ hotspot. These putative BoFLC1 copies are located ∼22 Kbp apart and may have been considered as a single copy in previous studies. In broccoli and cauliflower, this region has been implicated in vernalization requirement response (Bohuon et al., 1998; Rae et al., 1999), heat-tolerance (Branham et al., 2017), temperature dependent curd induction (Hasan et al., 2016), variation in heading response to temperature (Lin et al., 2019) and as a target for domestication due to a sharp increase in linkage disequilibrium when comparing improved broccoli and landrace broccoli genotypes (Stansell et al., 2018). Here, we identified HQ₉ as containing disproportionately many key heading-quality broccoli QTL: BU_C04@51.5, BR_C09@49.5, HC_C09@48.8, HU_C09@48.8, HE_C09@47.7, and OQ_C09@49.5, as well as two phenology QTL, DM_C09@50.0 and DF_C09@50.0.

Broccoli is characterized by strong suppression of bract elongation in the inflorescence and the “leaf in curd” bracting phenotype has been previously linked to high temperature stress (>22°C) (Booij and Struik, 1990). Kop et al. (2003) identified a role Boap1-a in bracting suppression within broccoli and cauliflower heads, although the authors suggested that additional candidates may be involved in bract development, e.g., BoFUL. Neither AP1 (Bo2g062650, Bo6g095760, Bo6g095760, Bo6g108600), or FUL (Bo2g161210, Bo7g098190, Bo9g014400) homologs cosegregated with bracting loci identified within the BolTBDH MQM bracting model.

Compact broccoli heads are less susceptible to damage and are more efficient to ship. In broccoli, the degree of head compactness is a consequence of a short rachis at a large angle. In a F₂B. oleracea cross, Lan and Paterson (2000) identified a C04 QTL (EW4D04w+7) likely implicated in curd density. A more complex head compactness model was observed in BolTBDH (HC ∼ HC_6@23.0 + HQ₉ + [HQ₉ × 3@53.6] + [HQ₉ × 5@1.4] + [HQ₈ × 3@53.6]).

Along with head compactness, head diameter at maturity is an important element of high yielding cultivars. In three B. oleracea F₂ populations, Lan and Paterson (2000) identified 14 curd width QTL on chromosomes C01, C03, C4, C5, C7, C8, and C9. Within a broccoli × broccoli mapping population, Walley et al. (2012) identified QTL for head diameter on linkage groups C2, C4, C6, C7, and C9. Head diameter within BolTBDH was captured by the MQM model HD ∼ HD_C05 @2.9 + HD_C07@43.6/46.5 + [HD_C05@2.9 × HD_C07@43.6/46.5] + [Bud₁ × HD_C05@2.9].

Head extension above lead rosette is a useful trait in broccoli by reducing labor during harvest. Alternatively, late stem elongation may be a useful to protect the head during growth. The MQM model of head extension is complex: HE ∼ HQ₉ + HQ₈ + C06@38.5 + [HQ₉ × 6@38.5] + [HQ₉ × 2@5.5].

Convex head shape is an important trait in broccoli by allowing the crown to shed water, thereby reducing disease incidence and exhibits complex control HC ∼ HQ₈ + C02@41.7 + C09@37.1 + [HQ₈ x HQ₉] + [2@41.7 × 9@37.1].

As expected, overall broccoli heading quality is the most complex of the head quality traits, best captured by the MQM model: OQ ∼ 3@57.0 + HQ₄ + HQ₈ + HQ₉ + [HQ₈ × 3@53.6] + [HQ₈ × HQ₉] + [HQ₉ × 5@1.3]. A simple control model of the broccoli heading-phenotype due to quantitative control of heading quality traits may be excluded given the MQM model determined in BolTBDH. Additionally, a constrictive-conditional model explaining the broccoli heading-phenotype seems less probable given the relative independence of heading quality traits. For example, the MQM models for HS and HU share HQ₉, but HS does not share HQ₄ with the overall heading quality model, nor does HU share HQ₈ with the overall heading quality model (Figure S4). The arrested-meristem broccoli heading-phenotype is best explained by a pleiotropic model — where a small number of genes within HQ₉ are implicated in multiple heading quality traits, and these HQ₉ genes exhibit important epistatic interactions with other heading quality loci.

Phenology

In a previous study of days to flowering within a different albogabra × italica mapping population, Bohuon et al. (1998) identified QTL on C02, C03, C05, and C09. In a QTL-seq analysis of broccoli × cabbage, Shu et al. (2018) identified three flowering time regions: C02@0.9–2.9 Mb, C03@1.8–20 Mb, and C06@5.0–5.6 Mb. Only the C03@1.8–20 Mb region collocated with the DF QTL identified in BolTBDH (DF_C03@6.4). Using a relative expression approach, Abuyusuf et al. (2019) report sequence based variations in BoFLC1.C9 (C09:51.0–51.0 Mb) implicated in early and late flowering cabbage genotypes, in agreement with the DM_C09@50.0 and DF_C09@50.0 QTL identified in BolTBDH. In a study of curd initiation in DH cauliflower, Hasan et al. (2016) identified a temperature-dependent time to curd induction QTL on C09 at 49.4 Mb, closely collocating with DM_C09@50.0 and DF_C09@50.0 and identified additional days to flowering QTL on C04, C05, C06, and C07. Within BolTBDH, days to head maturity and flowering appear to be each strongly influenced by two hotspots HQ₉ × PHE₃ given by the QTL pairs (DM C03_6.4 and DF_C03@6.4) and (DM_C09@49.7 and DF_C09@49.7) and these pairs exhibit strong epistatic effects ([DM_3@6.4 × DM_9@50.0]; p_fv1 << 0.01) and [DM_C09@50.0 × DF_C09@50.0] (pfv1 << 0.01). One additional DM epistatic effect was detected ([HQ9 x 7@DM_42.0]; p_fv1 = 0.023). Interestingly, HQ₉ collocates with the region of strongest genome-wide segregation distortion (P₁ allele: f = 0.76). If late-flowering DHs were underrepresented during tissue culture or seed regeneration, the P₁ allele and surrounding region would exhibit segregation enrichment.

Lan and Paterson (2000) identified 15 “days from budding to flowering” QTL, analogous to the holding ability trait measured in BolTHDH population. In BolTBDH, “days from head maturity trait to flowering” was best explained by a complex MQM model: HA ∼ 5@4.0 + 6@25.9 + 7@41.2 + [9@9.5 × 9@10.7] + [5@36.6 × 5@37.3] and no likely homologous candidates were identified.

Conclusions

Evaluation of the BolTBDH population provides new insights into key genomic regions and developmental candidates that define heading broccoli by identifying essential QTL implicated in these phenotypic outcomes. These results support a pleiotropic model of a heading broccoli phenotype. This work demonstrates several key genomic hotspots as essential for the phenotypes observed within this study, and these QTL and markers may prove useful for future marker-assisted breeding efforts. The phenomic and genomic dataset provided herein may be used for additional mapping studies and be integrated with previous work (e.g. metabolic and pathogen resistance studies).

Funding

This work is supported by Specialty Crop Research Initiative grant no. 2016-51181-25402 from the USDA National Institute of Food and Agriculture.

References

• 1

  AbuyusufM.NathU. K.KimH.-T.IslamM. R.ParkJ.-I.NouI.-S. (2019). Molecular markers based on sequence variation in BoFLC1.C9 for characterizing early-and late-flowering cabbage genotypes. BMC Genet.20, 42. doi: 10.1186/s12863-019-0740-1

  • CrossRef
  • Google Scholar

• 2

  ArendsD.PrinsP.JansenR. C.BromanK. W. (2010). R/qtl: high-throughput multiple QTL mapping. Bioinformatics26, 2990–2992. doi: 10.1093/bioinformatics/btq565

  • CrossRef
  • Google Scholar

• 3

  AtallahS. S.GómezM. I.BjörkmanT. (2014). Localization effects for a fresh vegetable product supply chain: broccoli in the eastern united states. Food Policy49, 151–159. doi: 10.1016/j.foodpol.2014.07.005

  • CrossRef
  • Google Scholar

• 4

  AxelssonT.ShavorskayaO.LagercrantzU. (2001). Multiple flowering time QTLs within several brassica species could be the result of duplicated copies of one ancestral gene. Genome44, 856–864. doi: 10.1139/g01-082

  • CrossRef
  • Google Scholar

• 5

  BerardiniT. Z.ReiserL.LiD.MezheritskyY.MullerR.StraitE.et al. (2015). The Arabidopsis information resource: making and mining the “gold standard” annotated reference plant genome. Genesis53, 474–485. doi: 10.1002/dvg.22877

  • CrossRef
  • Google Scholar

• 6

  BjörkmanT.PearsonK. J. (1998). High temperature arrest of inflorescence development in broccoli (Brassica oleracea var. Italica l.). J. Exp. Bot.49, 101–106. doi: 10.1093/jxb/49.318.101

  • CrossRef
  • Google Scholar

• 7

  BohuonE. J. R.RamsayL. D.CraftJ. A.ArthurA. E.MarshallD. F.LydiateD. J.et al. (1998). The association of flowering time quantitative trait loci with duplicated regions and candidate loci in Brassica oleracea. Genetics150, 393–401.

  • Google Scholar

• 8

  BooijR.StruikP. (1990). Effects of temperature on leaf and curd initiation in relation to juvenility of cauliflower. Sci. Hortic.44, 201–214. doi: 10.1016/0304-4238(90)90120-4

  • CrossRef
  • Google Scholar

• 9

  BoseonB.BilichakA.KovalchukI. (2019). Transgenerational response to heat stress in the form of differential expression of noncoding RNA fragments in Brassica rapa plants. Plant Genome12, 180022 doi: 10.3835/plantgenome2018.04.0022

  • CrossRef
  • Google Scholar

• 10

  BranhamS. E.FarnhamM. W. (2017). Genotyping-by-sequencing of waxy and glossy near-isogenic broccoli lines. Euphytica213, 84. doi: 10.1007/s10681-017-1873-9

  • CrossRef
  • Google Scholar

• 11

  BranhamS. E.StansellZ. J.CouillardD. M.FarnhamM. W. (2017). Quantitative trait loci mapping of heat tolerance in broccoli (Brassica oleracea var. Italica) using genotyping-by-sequencing. Theor. Appl. Genet.130, 529–538. doi: 10.1007/s00122-016-2832-x

  • CrossRef
  • Google Scholar

• 12

  BromanK. W.WuH.SenŚ.ChurchillG. A. (2003). R/qtl: QTL mapping in experimental crosses. Bioinformatics19, 889–890. doi: 10.1093/bioinformatics/btg112

  • CrossRef
  • Google Scholar

• 13

  BrownA. F.YousefG. G.ChebroluK. K.ByrdR. W.EverhartK. W.ThomasA.et al. (2014). High-density single nucleotide polymorphism (SNP) array mapping in Brassica oleracea: identification of qtl associated with carotenoid variation in broccoli florets. Theor. Appl. Genet.127, 2051–2064. doi: 10.1007/s00122-014-2360-5

  • CrossRef
  • Google Scholar

• 14

  CamargoL. E. A.OsbornT. C. (1996). Mapping loci controlling flowering time in Brassica oleracea. Theor. Appl. Genet.92, 610–616. doi: 10.1007/BF00224565

  • CrossRef
  • Google Scholar

• 15

  ChengF.SunR.HouX.ZhengH.ZhangF.ZhangY.et al. (2016). Subgenome parallel selection is associated with morphotype diversification and convergent crop domestication in Brassica rapa and Brassica oleracea. Nat. Genet.48, 1218. doi: 10.1038/ng.3634

  • CrossRef
  • Google Scholar

• 16

  DuclosD. V.BjörkmanT. (2008). Meristem identity gene expression during curd proliferation and flower initiation in Brassica oleracea. J. Exp. Bot.59, 421–433. doi: 10.1093/jxb/erm327

  • CrossRef
  • Google Scholar

• 17

  ElshireR. J.GlaubitzJ. C.SunQ.PolandJ. A.KawamotoK.BucklerE. S.et al. (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One6, e19379. doi: 10.1371/journal.pone.0019379

  • CrossRef
  • Google Scholar

• 18

  FarnhamM. W. (1998). Doubled-haploid broccoli production using anther culture: effect of anther source and seed set characteristics of derived lines. J. Am. Soc. Hortic. Sci.123, 73–77. doi: 10.21273/JASHS.123.1.73

  • CrossRef
  • Google Scholar

• 19

  FarnhamM. W.BjörkmanT. (2011). Evaluation of experimental broccoli hybrids developed for summer production in the eastern united states. HortScience46, 858–863. doi: 10.21273/HORTSCI.46.6.858

  • CrossRef
  • Google Scholar

• 20

  FranciscoM.AliM.FerreresF.MorenoD. A.VelascoP.SoengasP. (2016). Organ-specific quantitative genetics and candidate genes of phenylpropanoid metabolism in Brassica oleracea. Front. Plant Sci.6, 1240. doi: 10.3389/fpls.2015.01240

  • CrossRef
  • Google Scholar

• 21

  GaoM.LiG.McCombieW. R.QuirosC. F. (2005). Comparative analysis of a transposon-rich Brassica oleracea bac clone with its corresponding sequence in A. thaliana. Theor. Appl. Genet.111, 949–955. doi: 10.1007/s00122-005-0029-9

  • CrossRef
  • Google Scholar

• 22

  GaoM.LiG.YangB.QiuD.FarnhamM.QuirosC. (2007). High-density Brassica oleracea linkage map: identification of useful new linkages. Theor. Appl. Genet.115, 277–287. doi: 10.1007/s00122-007-0568-3

  • CrossRef
  • Google Scholar

• 23

  GillmorC. S.Silva-OrtegaC. O.WillmannM. R.Buendía-MonrealM.PoethigR. S. (2014). The Arabidopsis mediator CDK8 module genes CCT (MED12) and GCT (MED13) are global regulators of developmental phase transitions. Development141, 4580–4589. doi: 10.1242/dev.111229

  • CrossRef
  • Google Scholar

• 24

  GlaubitzJ. C.CasstevensT. M.LuF.HarrimanJ.ElshireR. J.SunQ.et al. (2014). TASSEL-GBS: a high capacity genotyping by sequencing analysis pipeline. PloS One9, e90346. doi: 10.1371/journal.pone.0090346

  • CrossRef
  • Google Scholar

• 25

  GoliczA. A.BayerP. E.BarkerG. C.EdgerP. P.KimH.MartinezP. A.et al. (2016). The pangenome of an agronomically important crop plant Brassica oleracea. Nat. Commun.7, 13390. doi: 10.1038/ncomms13390

  • CrossRef
  • Google Scholar

• 26

  GreerS.WenM.BirdD.WuX.SamuelsL.KunstL.et al. (2007). The cytochrome p450 enzyme CYP96a15 is the midchain alkane hydroxylase responsible for formation of secondary alcohols and ketones in stem cuticular wax of Arabidopsis. Plant Physiol.145, 653–667. doi: 10.1104/pp.107.107300

  • CrossRef
  • Google Scholar

• 27

  GrompingU. (2006). Relative importance for linear regression in R: the package relaimpo. J. Stat. Software17(1), 1–27. doi: 10.18637/jss.v017.i01

  • CrossRef
  • Google Scholar

• 28

  HanF.CuiH.ZhangB.LiuX.YangL.ZhuangM.et al. (2019). Map-based cloning and characterization of BoCCD4, a gene responsible for white/yellow petal color in B. oleracea. BMC Genomics20, 242. doi: 10.1186/s12864-019-5596-2

  • CrossRef
  • Google Scholar

• 29

  HasanY.BriggsW.MatschegewskiC.OrdonF.StützelH.ZetzscheH.et al. (2016). Quantitative trait loci controlling leaf appearance and curd initiation of cauliflower in relation to temperature. TAG. Theor. Appl. Genet.129, 1273–1288. doi: 10.1007/s00122-016-2702-6

  • CrossRef
  • Google Scholar

• 30

  HatzigS. V.FrischM.BreuerF.NesiN.DucournauS.WagnerM.-H.et al. (2015). Genome-wide association mapping unravels the genetic control of seed germination and vigor in Brassica napus. Front. Plant Sci.6, 221. doi: 10.3389/fpls.2015.00221

  • CrossRef
  • Google Scholar

• 31

  HeY.WuD.WeiD.FuY.CuiY.DongH.et al. (2017). Gwas, qtl mapping and gene expression analyses in Brassica napus reveal genetic control of branching morphogenesis. Sci. Rep.7, 15971. doi: 10.1038/s41598-017-15976-4

  • CrossRef
  • Google Scholar

• 32

  HuJ.SadowskiJ.OsbornT. C.LandryB. S.QuirosC. F. (1998). Linkage group alignment from four independent Brassica oleracea RFLP maps. Genome41, 226–235. doi: 10.1139/g98-007

  • CrossRef
  • Google Scholar

• 33

  Iglesias-BernabéL.MadlooP.RodríguezV. M.FranciscoM.SoengasP. (2019). Dissecting quantitative resistance to Xanthomonas campestris pv. campestris in leaves of Brassica oleracea by QTL analysis. Sci. Rep.9, 2015. doi: 10.1038/s41598-019-38527-5

  • CrossRef
  • Google Scholar

• 34

  Iniguez-LuyF. L.LukensL.FarnhamM. W.AmasinoR. M.OsbornT. C. (2009). Development of public immortal mapping populations, molecular markers and linkage maps for rapid cycling Brassica rapa and Brassica oleracea. Theor. Appl. Genet.120, 31–43. doi: 10.1007/s00122-009-1157-4

  • CrossRef
  • Google Scholar

• 35

  IrwinJ. A.SoumpourouE.ListerC.LigthartJ.KennedyS.DeanC. (2016). Nucleotide polymorphism affecting FLC expression underpins heading date variation in horticultural brassicas. Plant J.87, 597–605. doi: 10.1111/tpj.13221

  • CrossRef
  • Google Scholar

• 36

  KianianS. F.QuirosC. F. (1992). Generation of a brassica oleracea composite RFLP map: linkage arrangements among various populations and evolutionary implications. Theor. Appl. Genet.84, 544–554. doi: 10.1007/BF00224150

  • CrossRef
  • Google Scholar

• 37

  KopE. P.TeakleG. R.McClenaghanE. R.LynnJ. R.KingG. J. (2003). Genetic analysis of the bracting trait in cauliflower and broccoli. Plant Sci.164, 803–808. doi: 10.1016/S0168-9452(03)00068-2

  • CrossRef
  • Google Scholar

• 38

  LagercrantzU.Kruskopf OsterbergM.LascouxM. (2002). Sequence variation and haplotype structure at the putative flowering-time locus COL1 of Brassica nigra. Mol. Biol. Evol.19, 1474–1482. doi: 10.1093/oxfordjournals.molbev.a004210

  • CrossRef
  • Google Scholar

• 39

  LanT.PatersonA. H. (2001). Comparative mapping of QTLs determining the plant size of Brassica oleracea. Theor. Appl. Genet.103, 383–397. doi: 10.1007/s001220100615

  • CrossRef
  • Google Scholar

• 40

  LanT. H.PatersonA. H. (2000). Comparative mapping of quantitative trait loci sculpting the curd of Brassica oleracea. Genetics155, 1927–1954. doi: 10.1007/s001220100615

  • CrossRef
  • Google Scholar

• 41

  LandryB. S.HubertN.CreteR.ChangM. S.LincolnS. E.EtohT. (1992). A genetic map for Brassica oleracea based on RFLP markers detected with expressed DNA sequences and mapping of resistance genes to race 2 of Plasmodiophora brassicae (woronin). Genome35, 409–420. doi: 10.1139/g92-061

  • CrossRef
  • Google Scholar

• 42

  LeeJ.IzzahN. K.ChoiB.-S.JohH. J.LeeS.-C.PerumalS.et al. (2015a). Genotyping-by-sequencing map permits identification of clubroot resistance QTLs and revision of the reference genome assembly in cabbage (Brassica oleracea l.). DNA Res.23, 29–41. doi: 10.1093/dnares/dsv034

  • CrossRef
  • Google Scholar

• 43

  LeeJ.YangK.LeeM.KimS.KimJ.LimS.et al. (2015). Differentiated cuticular wax content and expression patterns of cuticular wax biosynthetic genes in bloomed and bloomless broccoli (Brassica oleracea var. Italica). Process Biochem.50, 456–462. doi: 10.1016/j.procbio.2014.12.012

  • CrossRef
  • Google Scholar

• 44

  LeijtenW.KoesR.RoobeekI.FrugisG. (2018). Translating flowering time from Arabidopsis thaliana to Brassicaceae and Asteraceae crop species. Plants7, 111. doi: 10.3390/plants7040111

  • CrossRef
  • Google Scholar

• 45

  LiF.KitashibaH.InabaK.NishioT. (2009). A Brassica rapa linkage map of EST-based SNP markers for identification of candidate genes controlling flowering time and leaf morphological traits. DNA Res.16, 311–323. doi: 10.1093/dnares/dsp020

  • CrossRef
  • Google Scholar

• 46

  LiG.GaoM.YangB.QuirosC. F. (2003). Gene for gene alignment between the Brassica and Arabidopsis genomes by direct transcriptome mapping. Theor. Appl. Genet.107, 168–180. doi: 10.1007/s00122-003-1236-x

  • CrossRef
  • Google Scholar

• 47

  LiG.ZhangG.ZhangY.LiuK.LiT.ChenH. (2015). Identification of quantitative trait loci for bolting and flowering times in Chinese kale (Brassica oleracea var. Alboglabra) based on ssr and srap markers. J. Hortic. Sci. Biotechnol.90, 728–737. doi: 10.1080/14620316.2015.11668739

  • CrossRef
  • Google Scholar

• 48

  LiH. (2013). “Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM” in arXiv (Cambridge, USA: Broad Institute of Harvard and MIT), 1–3.

  • Google Scholar

• 49

  LiZ.MeiY.LiuY.FangZ.YangL.ZhuangM.et al. (2019). The evolution of genetic diversity of broccoli cultivars in China since 1980. Sci. Hortic.250, 69–80. doi: 10.1016/j.scienta.2019.02.034

  • CrossRef
  • Google Scholar

• 50

  LinC.-W.FuS.-F.LiuY.-J.ChenC.-C.ChangC.-H.YangY.-W.et al. (2019). Analysis of ambient temperature-responsive transcriptome in shoot apical meristem of heat-tolerant and heat-sensitive broccoli inbred lines during floral head formation. BMC Plant Biol.19, 3. doi: 10.1186/s12870-018-1613-x

  • CrossRef
  • Google Scholar

• 51

  LinK.ChangL.LaiC.LoH. (2013). Aflp mapping of quantitative trait loci influencing seven head-related traits in broccoli (Brassica oleracea var. Italica). J. Hortic. Sci. Biotechnol.88, 257–268. doi: 10.1080/14620316.2013.11512964

  • CrossRef
  • Google Scholar

• 52

  LinS.-I.WangJ.-G.PoonS.-Y.WangS.-S.ChiouT.-J. (2005). Differential regulation of FLOWERING LOCUS C expression by vernalization in cabbage and Arabidopsis. Plant Physiol.137, 1037–1048. doi: 10.1104/pp.104.058974

  • CrossRef
  • Google Scholar

• 53

  LinY.-R.LeeJ.-Y.TsengM.-C.LeeC.-Y.ShenC.-H.WangC.-S.et al. (2018). Subtropical adaptation of a temperate plant (Brassica oleracea var. Italica) utilizes non-vernalization-responsive QTLs. Sci. Rep.8, 13609. doi: 10.1038/s41598-018-31987-1

  • CrossRef
  • Google Scholar

• 54

  LiuS.LiuY.YangX.TongC.EdwardsD.ParkinI. A. P.et al. (2014). The Brassica oleracea genome reveals the asymmetrical evolution of polyploid genomes. Nat. Commun.5, 3930. doi: 10.1038/ncomms4930

  • CrossRef
  • Google Scholar

• 55

  LuK.PengL.ZhangC.LuJ.YangB.XiaoZ.et al. (2017). Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits in Brassica napus. Front. Plant Sci.8, 206. doi: 10.3389/fpls.2017.00206

  • CrossRef
  • Google Scholar

• 56

  LuK.WeiL.LiX.WangY.WuJ.LiuM.et al. (2019). Whole-genome resequencing reveals Brassica napus origin and genetic loci involved in its improvement. Nat. Commun.10, 1154. doi: 10.1038/s41467-019-09134-9

  • CrossRef
  • Google Scholar

• 57

  MatschegewskiC.ZetzscheH.HasanY.LeibeguthL.BriggsW.OrdonF.et al. (2015). Genetic variation of temperature-regulated curd induction in cauliflower: elucidation of floral transition by genome-wide association mapping and gene expression analysis. Front. Plant Sci.6, 720. doi: 10.3389/fpls.2015.00720

  • CrossRef
  • Google Scholar

• 58

  MunthaS. T.ZhangL.ZhouY.ZhaoX.HuZ.YangJ.et al. (2018). Phytochrome A signal transduction 1 and CONSTANS-LIKE 13 coordinately orchestrate shoot branching and flowering in leafy Brassica juncea. Plant Biotechnol. J.17, 1333–1343. doi: 10.1111/pbi.13057

  • CrossRef
  • Google Scholar

• 59

  NiX.LiuH.HuangJ.ZhaoJ. (2017). LMI1-like genes involved in leaf margin development of Brassica napus. Genetica145, 269–274. doi: 10.1007/s10709-017-9963-0

  • CrossRef
  • Google Scholar

• 60

  OkazakiK.SakamotoK.KikuchiR.SaitoA.TogashiE.KuginukiY.et al. (2007). Mapping and characterization of FLC homologs and QTL analysis of flowering time in Brassica oleracea. Theor. Appl. Genet.114, 595–608. doi: 10.1007/s00122-006-0460-6

  • CrossRef
  • Google Scholar

• 61

  O’neillC. M.BancroftI. (2000). Comparative physical mapping of segments of the genome of Brassica oleracea var. Alboglabra that are homoeologous to sequenced regions of chromosomes 4 and 5 of. Plant J.23, 233–243. doi: 10.1046/j.1365-313x.2000.00781.x

  • CrossRef
  • Google Scholar

• 62

  OsbornT. C.KoleC.ParkinI. A.SharpeA. G.KuiperM.LydiateD. J.et al. (1997). Comparison of flowering time genes in Brassica rapa, Brassica napus and Arabidopsis thaliana. Genetics146, 1123–1129.

  • Google Scholar

• 63

  ParkinI. A.KohC.TangH.RobinsonS. J.KagaleS.ClarkeW. E.et al. (2014). Transcriptome and methylome profiling reveals relics of genome dominance in the mesopolyploid Brassica oleracea. Genome Biol.15, R77. doi: 10.1186/gb-2014-15-6-r77

  • CrossRef
  • Google Scholar

• 64

  PinkD.BaileyL.McClementS.HandP.MathasE.Buchanan-WollastonV.et al. (2008). Double haploids, markers and QTL analysis in vegetable brassicas. Euphytica164, 509–514. doi: 10.1007/s10681-008-9742-1

  • CrossRef
  • Google Scholar

• 65

  PottersG.PasternakT. P.GuisezY.PalmeK. J.JansenM. A. (2007). Stress-induced morphogenic responses: growing out of trouble? Trends Plant Sci.12, 98–105. doi: 10.1016/j.tplants.2007.01.004

  • CrossRef
  • Google Scholar

• 66

  R-Core-Team. (2018). R: A language and environment for statistical computing. FIX!

  • Google Scholar

• 67

  R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.

  • Google Scholar

• 68

  RaeA. M.HowellE. C.KearseyM. J. (1999). More QTL for flowering time revealed by substitution lines in Brassica oleracea. Heredity83, 586. doi: 10.1038/sj.hdy.6886050

  • CrossRef
  • Google Scholar

• 69

  RamsayL.JenningsD.KearseyM. J.MarshallD. F.BohuonE. J. R.ArthurA. E.et al. (1996). The construction of a substitution library of recombinant backcross lines in Brassica oleracea for the precision mapping of quantitative trait loci. Genome39, 558–567. doi: 10.1139/g96-071

  • CrossRef
  • Google Scholar

• 70

  RaziH.HowellE. C.NewburyH. J.KearseyM. J. (2008). Does sequence polymorphism of FLC paralogues underlie flowering time QTL in Brassica oleracea? TAG. Theor. Appl. Genet.116, 179–192. doi: 10.1007/s00122-007-0657-3

  • CrossRef
  • Google Scholar

• 71

  RenJ.LiuZ.DuJ.FuW.HouA.FengH. (2019). Fine-mapping of a gene for the lobed leaf, boll, in ornamental kale (Brassica oleracea L. var. Acephala). Mol. Breed.39, 40. doi: 10.1007/s11032-019-0944-0

  • CrossRef
  • Google Scholar

• 72

  RidgeS.BrownP. H.HechtV.DriessenR. G.WellerJ. L. (2015). The role of BoFLC2 in cauliflower (Brassica oleracea var. Botrytis l.) reproductive development. J. Exp. Bot.66, 125–135. doi: 10.1093/jxb/eru408

  • CrossRef
  • Google Scholar

• 73

  RoederA. H. K.ChickarmaneV.CunhaA.ObaraB.ManjunathB. S.MeyerowitzE. M. (2010). Variability in the control of cell division underlies sepal epidermal patterning in Arabidopsis thaliana. PLoS Biol.8, e1000367. doi: 10.1371/journal.pbio.1000367

  • CrossRef
  • Google Scholar

• 74

  RuffierM.KähäriA.KomorowskaM.KeenanS.LairdM.LongdenI.et al. (2017). Ensembl core software resources: storage and programmatic access for DNA sequence and genome annotation. (Oxford)2017, (1) bax020. doi: 10.1093/database/bax020

  • CrossRef
  • Google Scholar

• 75

  SchiesslS.Iniguez-LuyF.QianW.SnowdonR. J. (2015). Diverse regulatory factors associate with flowering time and yield responses in winter-type Brassica napus. BMC Genomics16, 737. doi: 10.1186/s12864-015-1950-1

  • CrossRef
  • Google Scholar

• 76

  SebastianR. L.HowellE. C.KingG. J.MarshallD. F.KearseyM. J. (2000). An integrated AFLP and RFLP Brassica oleracea linkage map from two morphologically distinct doubled-haploid mapping populations. Theor. Appl. Genet.100, 75–81. doi: 10.1007/s001220050011

  • CrossRef
  • Google Scholar

• 77

  SebastianR. L.KearseyM. J.KingG. J. (2002). Identification of quantitative trait loci controlling developmental characteristics of Brassica oleracea L. Theor. Appl. Genet.104, 601–609. doi: 10.1007/s001220100743

  • CrossRef
  • Google Scholar

• 78

  ShahS.WeinholdtC.JedrusikN.MolinaC.ZouJ.GroßeI.et al. (2018). Whole-transcriptome analysis reveals genetic factors underlying flowering time regulation in rapeseed (Brassica napus L.). Plant Cell Environ.41, 1935–1947. doi: 10.1111/pce.13353

  • CrossRef
  • Google Scholar

• 79

  SheaD. J.ItabashiE.TakadaS.FukaiE.KakizakiT.FujimotoR.et al. (2018). The role of FLOWERING LOCUS C in vernalization of Brassica: the importance of vernalization research in the face of climate change. Crop Pasture Sci.69, 30–39. doi: 10.1071/CP16468

  • CrossRef
  • Google Scholar

• 80

  ShenY.XiangY.XuE.GeX.LiZ. (2018). Major co-localized QTL for plant height, branch initiation height, stem diameter, and flowering time in an alien introgression derived Brassica napus DH population. Front. Plant Sci.9, 390. doi: 10.3389/fpls.2018.00390

  • CrossRef
  • Google Scholar

• 81

  ShimanoS.HibaraK.-I.FuruyaT.ArimuraS.-I.TsukayaH.ItohJ.-I. (2018). Conserved functional control, but distinct regulation, of cell proliferation in rice and Arabidopsis leaves revealed by comparative analysis of GRF-INTERACTING FACTOR 1 orthologs. Development145, dev159624. doi: 10.1242/dev.159624

  • CrossRef
  • Google Scholar

• 82

  ShuJ.LiuY.ZhangL.LiZ.FangZ.YangL.et al. (2018). Qtl-seq for rapid identification of candidate genes for flowering time in broccoli × cabbage. Theor. Appl. Genet.131, 917–928. doi: 10.1007/s00122-017-3047-5

  • CrossRef
  • Google Scholar

• 83

  SiriwardanaN. S.LambR. S. (2012). The poetry of reproduction: the role of leafy in Arabidopsis thaliana flower formation. Int. J. Dev. Biol.56, 207–221. doi: 10.1387/ijdb.113450ns

  • CrossRef
  • Google Scholar

• 84

  SoteloT.CarteaM. E.VelascoP.SoengasP. (2014). Identification of antioxidant capacity-related QTLs in Brassica oleracea. PLoS One9(9), e107290. doi: 10.1371/journal.pone.0107290

  • CrossRef
  • Google Scholar

• 85

  Sotelo PérezT.Francisco CandeiraM.Cartea GonzálezM. E.Rodríguez GrañaV. M.Soengas FernándezM. D. P. (2014). QTLs controlling antioxidant capacity in leaves and flower buds of Brassica oleracea. PLoS One. doi: 10.1371/journal.pone.0107290

  • CrossRef
  • Google Scholar

• 86

  StansellZ.BjörkmanT.BranhamS.CouillardD.FarnhamM. W. (2017). Use of a quality trait index to increase the reliability of phenotypic evaluations in broccoli. HortScience52, 1490–1495. doi: 10.21273/HORTSCI12202-17

  • CrossRef
  • Google Scholar

• 87

  StansellZ.HymaK.Fresnedo-RamírezJ.SunQ.MitchellS.BjörkmanT.et al. (2018). Genotyping-by-sequencing of Brassica oleracea vegetables reveals unique phylogenetic patterns, population structure and domestication footprints. Hortic. Res.5, 38. doi: 10.1038/s41438-018-0040-3

  • CrossRef
  • Google Scholar

• 88

  SunX.BucherJ.JiY.van DijkA. D. J.ImminkR. G. H.BonnemaG. (2018). Effect of ambient temperature fluctuation on the timing of the transition to the generative stage in cauliflower. Environ. Exp. Bot.155, 742–750. doi: 10.1016/j.envexpbot.2018.06.013

  • CrossRef
  • Google Scholar

• 89

  SwartsK.LiH.NavarroJ. A. R.AnD.RomayM. C.HearneS.et al. (2014). Novel methods to optimize genotypic imputation for low-coverage, next-generation sequence data in crop plants. Plant Genome7, 0. doi: 10.3835/plantgenome2014.05.0023

  • CrossRef
  • Google Scholar

• 90

  TortosaM.CarteaM. E.RodríguezV. M.VelascoP. (2018). Unraveling the metabolic response of Brassica oleracea exposed to Xanthomonas campestris pv. J. Sci. Food Agric.98, 3675–3683. doi: 10.1002/jsfa.8876

  • CrossRef
  • Google Scholar

• 91

  TyagiS.SriT.SinghA.MayeeP.ShivarajS.SharmaP.et al. (2019). Suppressor of overexpression of constans1 influences flowering time, lateral branching, oil quality, and seed yield in Brassica juncea cv. Funct. Integr. Genomics19, 43–60. doi: 10.1007/s10142-018-0626-8

  • CrossRef
  • Google Scholar

• 92

  UptmoorR.SchragT.StützelH.EschE. (2008). Crop model based QTL analysis across environments and QTL based estimation of time to floral induction and flowering in Brassica oleracea. Mol. Breed.21, 205–216. doi: 10.1007/s11032-007-9121-y

  • CrossRef
  • Google Scholar

• 93

  VaraudE.BrioudesF.SzécsiJ.LerouxJ.BrownS.Perrot-RechenmannC.et al. (2011). AUXIN RESPONSE FACTOR8 regulates Arabidopsis petal growth by interacting with the bHLH transcription factor BIGPETALp. Plant Cell23, 973–983. doi: 10.1105/tpc.110.081653

  • CrossRef
  • Google Scholar

• 94

  WalleyP. G.CarderJ.SkipperE.MathasE.LynnJ.PinkD.et al. (2012). A new broccoli × broccoli immortal mapping population and framework genetic map: tools for breeders and complex trait analysis. TAG. Theor. Appl. Genet.124, 467–484. doi: 10.1007/s00122-011-1721-6

  • CrossRef
  • Google Scholar

• 95

  WangB.ZhouX.XuF.GaoJ. (2010). Ectopic expression of a Chinese cabbage brargos gene in Arabidopsis increases organ size. Transgenic Res.19, 461–472. doi: 10.1007/s11248-009-9324-6

  • CrossRef
  • Google Scholar

• 96

  XuX.LuoW.GuoJ.ChenH.AkramW.XieD.et al. (2019). Fine mapping and candidate gene analysis of the yellow petal gene ckpc in Chinese kale (Brassica oleracea L. var. Alboglabra bailey) by whole-genome resequencing. Mol. Breed.39, 96. doi: 10.1007/s11032-019-1011-6

  • CrossRef
  • Google Scholar

• 97

  YousefE. A. A.MüllerT.BörnerA.SchmidK. J. (2018). Comparative analysis of genetic diversity and differentiation of cauliflower (Brassica oleracea var. Botrytis) accessions from two ex situ genebanks. PLoS One13, e0192062. doi: 10.1371/journal.pone.0192062

  • CrossRef
  • Google Scholar

• 98

  ZhangB.LiuC.WangY.YaoX.WangF.WuJ.et al. (2015). Disruption of a CAROTENOID CLEAVAGE DIOXYGENASE 4 gene converts flower colour from white to yellow in Brassica species. New Phytol.206, 1513–1526. doi: 10.1111/nph.13335

  • CrossRef
  • Google Scholar

• 99

  ZhangF.WangH.KalveS.WolabuT. W.NakashimaJ.GolzJ. F.et al. (2019). Control of leaf blade outgrowth and floral organ development by LEUNIG, ANGUSTIFOLIA 3 and WOX transcriptional regulators. New Phytol.223, 2024–2038. doi: 10.1111/nph.15921

  • CrossRef
  • Google Scholar

• 100

  ZhangY.HuangS.WangX.LiuJ.GuoX.MuJ.et al. (2018). Defective APETALA2 genes lead to sepal modification in Brassica crops. Front. Plant Sci.9, 367. doi: 10.3389/fpls.2018.00367

  • CrossRef
  • Google Scholar

• 101

  ZhengM.HuM.YangH.TangM.ZhangL.LiuH.et al. (2019). Three BnaIAA7 homologs are involved in auxin/brassinosteroid-mediated plant morphogenesis in rapeseed (brassica napus L.). Plant Cell Rep.1–15. doi: 10.1007/s00299-019-02410-4

  • CrossRef
  • Google Scholar

• 102

  ZhuX.TaiX.RenY.ChenJ.BoT. (2019). Genome-wide analysis of coding and long non-coding RNAs involved in cuticular wax biosynthesis in cabbage (Brassica oleracea L. var. Capitata). Int. J. Mol. Sci.20, 2820. doi: 10.3390/ijms20112820

  • CrossRef
  • Google Scholar