Upregulation of tandem duplicated BoFLC1 genes is associated with the non‑flowering trait in Brassica oleracea var. capitata

Key message Tandem duplicated BoFLC1 genes ( BoFLC1a and BoFLC1b ), which were identified as the candidate causal genes for the non-flowering trait in the cabbage mutant ‘nfc’, were upregulated during winter in ‘nfc’. Abstract The non-flowering natural cabbage mutant ‘nfc’ was discovered from the breeding line ‘T15’ with normal flowering characteristics. In this study, we investigated the molecular basis underlying the non-flowering trait of ‘nfc’. First, ‘nfc’ was induced to flower using the grafting floral induction method, and three F 2 populations were generated. The flowering phenotype of each F 2 population was widely distributed with non-flowering individuals appearing in two populations. QTL-seq analysis detected a genomic region associated with flowering date at approximately 51 Mb on chromosome 9 in two of the three F 2 populations. Subsequent validation and fine mapping of the candidate genomic region using QTL analysis identified the quantitative trait loci (QTL) at 50,177,696–51,474,818 bp on chromosome 9 covering 241 genes. Additionally, RNA-seq analysis in leaves and shoot apices of ‘nfc’ and ‘T15’ plants identified 19 and 15 differentially expressed genes related to flowering time, respectively. Based on these results, we identified tandem duplicated BoFLC1 genes, which are homologs of the floral repressor FLOWERING LOCUS C , as the candidate genes responsible for the non-flowering trait of ‘nfc’. We designated the tandem duplicated BoFLC1 genes as BoFLC1a and BoFLC1b . Expression analysis revealed that the expression levels of BoFLC1a and BoFLC1b were downregulated during winter in ‘T15’ but were upregulated and maintained during winter in ‘nfc’. Additionally, the expression level of the floral integrator BoFT was upregulated in the spring in ‘T15’ but hardly upregulated in ‘nfc’. These results suggest that the upregulated levels of BoFLC1a and BoFLC1b contributed to the non-flowering trait of ‘nfc’.


Introduction
Cabbage (Brassica oleracea var. capitata), a leafy vegetable belonging to the family Brassicaceae, is cultivated worldwide. Many plants belonging to the family Brassicaceae are promoted to flower after exposure to low temperatures, and this process is known as vernalization. Flower bud development and stem elongation (bolting) significantly reduce the yield and commercial value of Brassicaceae crops whose Communicated by Lixi Jiang.
41 Page 2 of 23 vegetative organs such as leaves and roots are harvested. Therefore, breeding cultivars with a late bolting trait is required for the stable production of these crops in the spring (Jung and Müller 2009). Meanwhile, the extremely late bolting trait is a trade-off for seed production efficiency. Hence, the breeding material and knowledge of a late bolting trait are limited.
Flowering is the most important event in the plant life cycle. In the model plant Arabidopsis thaliana, the following six flowering pathways have been elucidated: vernalization, photoperiod, autonomous, gibberellin, age, and ambient temperature pathways (Kim et al. 2009;Whittaker and Dean 2017). In particular, the vernalization requirement plays an important role in local adaptability and has been well studied (Stinchcombe et al. 2004(Stinchcombe et al. , 2005Shindo et al. 2005Shindo et al. , 2006. Vernalization response can be classified into the following two types according to the age at which the plants can respond to low temperatures (Friend 1985): seed-vernalization-type, in which plants, including A. thaliana and Brassica rapa, can sense low temperatures at the imbibed seed stage; plant-vernalization-type, in which plants, including Arabis alpina and B. oleracea, can sense low temperatures only after they reach a certain plant size (Miller 1929;Ito and Saito 1961;Hyun et al. 2019). In Arabidopsis, the floral repressor FLOWERING LOCUS C (FLC), encoding a MADS-box transcription factor, plays a central role in vernalization and autonomous pathway (Sheldon et al. 1999(Sheldon et al. , 2000Amasino 1999, 2001). FLC represses the expression of the floral integrator FLOWERING LOCUS T (FT), SUPPRESSOR OF OVEREXPRESSION OF CON-STANS1 (SOC1), and FD by binding directly to the first intron of FT and the promoter regions of SOC1 and FD (Lee et al. 2000;Helliwell et al. 2006;Searle et al. 2006). Vernalization epigenetically represses the expression of FLC, enabling the activation of these flowering-related genes and the subsequent induction of flower bud differentiation (Bastow et al. 2004;Sung and Amasino 2004).
Four FLC homologs [BoFLC1, BoFLC2 (also known as BoFLC4), BoFLC3, and BoFLC5] have been reported in B. oleracea (Schranz et al. 2002;Lin et al. 2005;Okazaki et al. 2007). BoFLC1, BoFLC2, and BoFLC3 have conserved functions as floral repressors (Lin et al. 2005;Irwin et al. 2016;Itabashi et al. 2019). B. oleracea includes many varieties with diverse vernalization requirements. The correlation between each BoFLC homolog and flowering time has been previously reported (Lin et al. 2005(Lin et al. , 2018Okazaki et al. 2007;Razi et al. 2008;Ridge et al. 2015;Irwin et al. 2016;Abuyusuf et al. 2019;Guo et al. 2021;Li et al. 2022). Comparative analysis of early-flowering and late-flowering cabbage lines revealed that a 67 bp indel in the second intron of BoFLC1 correlates with flowering time (Abuyusuf et al. 2019). Quantitative trait loci (QTL) analysis revealed that BoFLC2 (also known as BoFLC4) is strongly associated with the vernalization requirement in broccoli (B. oleracea var. italica) and cabbage and that the non-vernalization type had a non-functional allele with a frameshift mutation in the fourth exon (Okazaki et al. 2007). These functional and non-functional alleles of BoFLC2 were correlated with curd induction and flowering time in cauliflower (B. oleracea var. botrytis) (Ridge et al. 2015). Other polymorphisms in BoFLC2 have been reported to be involved in flowering time in broccoli and cabbage (Irwin et al. 2016;Li et al. 2022). BoFLC3 is involved in the timing of curd induction in subtropical broccoli breeding lines (Lin et al. 2018). Comparative analysis of cabbage and cauliflower revealed the presence of a 263 bp deletion and a 49 bp insertion at the first intron of BoFLC3 cauliflower allele (Guo et al. 2021). BoFLC5 is considered a pseudogene based on the frameshift mutation in the second exon or the lack of expression (Pires et al. 2004;Okazaki et al. 2007;Razi et al. 2008). FT, a target gene of FLC, has the following two homologs in B. oleracea: BoFT.C2 and BoFT.C6. BoFT.C2 is not transcribed or exhibits extremely low expression levels (Wang et al. 2012;Irwin et al. 2016). The overexpression of BoFT. C6 induced flowering in Arabidopsis (Itabashi et al. 2019). The expression level of FT homologs correlates with that of FLC homologs in several B. oleracea species (Ridge et al. 2015;Irwin et al. 2016).
The non-flowering natural cabbage mutant 'nfc' was discovered from the breeding line 'T15' in the spring of 1978 (Kinoshita et al. 2021). Since its discovery, 'nfc' plants have been propagated using cutting. 'nfc' plants rarely flower even under floral inductive conditions where the flowering rate of 'T15' and other commercial cabbage cultivars is 100% (Kinoshita et al. 2021). 'nfc' has potential value as a breeding material for late-bolting breeding. Additionally, the analysis of 'nfc' plants can provide useful insights into the factors required for flowering in B. oleracea. However, the molecular basis underlying the non-flowering mechanism of the 'nfc' has not been elucidated. QTL mapping methods with segregating populations are often used to identify the gene responsible for a certain phenotype in seed-propagating plants. This genetically based QTL mapping method has the advantage of identifying genomic regions associated with a phenotype regardless of the expression site or pattern of the causal gene, the mutation type, or the expression changes in the genes affected by the causal gene. Flowering, in particular, is a dynamic and systemic process that accompanies the expression changes of various gene clusters depending on the developmental stage of the plant and the surrounding environment. Therefore, the QTL mapping method is thought to be an effective means to identify the causal gene responsible for the nonflowering trait of 'nfc'. However, the application of QTL mapping methods in 'nfc' plants is challenging as they Page 3 of 23 41 rarely flower. Recently, a method has been developed to induce the flowering of cabbage without vernalization by grafting cabbage onto radishes (Raphanus sativus) (Motoki et al. 2019;Motoki et al. 2022;Motoki et al. 2023). This method enables the induction of flowering in the non-flowering cabbage mutant 'nfc' (unpublished data). Subsequently, a mapping population for gene identification can be generated.
This study aimed to identify the candidate gene responsible for the non-flowering trait of 'nfc'. We generated F 2 segregating populations by crossing the flowering 'nfc' scion obtained via grafting with other cultivars. Phenotypic evaluation of the F 2 population and subsequent QTLseq analysis and QTL analysis revealed the QTL involved in the non-flowering trait of 'nfc'. Additionally, we identified differentially expressed genes between 'nfc' and 'T15' using RNA sequencing (RNA-seq) analysis. Based on these results, we identified the candidate genes responsible for the non-flowering trait of 'nfc'.

Plant materials and population construction
The 'nfc' was subjected to the non-vernalization grafting floral induction method to artificially induce flowering, following the previously reported methods (Motoki et al. 2019;Motoki et al. 2022;Motoki et al. 2023). To generate segregating populations, flowering 'nfc' scions via grafting were crossed with three B. oleracea cultivars exhibiting different degrees of relatedness to 'nfc'. The most closely related cultivar was 'T15' (B. oleracea var. capitata), which is the parental line of 'nfc'. One individual of 'T15'(#43) was used for crossing. The second most closely related cultivar was a commercial cabbage cultivar 'Watanabe-seiko No. 1 (W1)' (B. oleracea var. capitata) (Genebank Project, NARO, Japan, accession No. 25974), and one individual was used for crossing. The most distantly related cultivar was a commercial Chinese kale cultivar 'kairan' (B. oleracea var. alboglabra) (Tsurushin Seed Ltd., Japan), and one individual was used for crossing. Each F 1 was cultivated in the open field as described below, and each F 2 seed was obtained after the self-pollination of six individuals of 'nfc' × 'T15' F 1 , one individual of 'nfc' × 'W1' F 1 , and three individuals of 'kairan' × 'nfc' F 1 . Self-pollination was performed using bud pollination or CO 2 treatment to overcome self-incompatibility, following the methods of previous studies (Nakanishi and Hinata 1975). To prevent cross-pollination, the inflorescences were covered with non-woven bags before flowering until pod formation.

Growth conditions in the field
For 'nfc', plants propagated in vitro via cutting were used. For 'T15', plants propagated in vitro via cutting or cultivated from seeds were used. For 'W1', 'kairan', F 1 , and F 2 , plants cultivated from seeds were used. Plants were cultivated between 2018 and 2021. Cultivation year represents the year the plants were transplanted into the open field. Seedlings were raised as described previously (Kinoshita et al. 2021). The seedlings were transplanted into the open field of the Kyoto Farm or the Kizu Farm, Kyoto University in September or October. Cultivation, fertilization, and pest management were performed following commercial practices. The maximum leaf length of each plant was measured on January 19, , December 15-18, 2019, December 19-20, 2020, and December 18-19, 2021 of each cultivation year. Plants with defective growth were excluded from the analysis. Since March of the following year after transplanting, cabbage heads were divided vertically with a knife to enable the main stem to bolt and the lateral shoots to elongate. Cultivation was continued until July of the following year.

Quantitative evaluation of flowering characteristics
'nfc' is not a mutant that has lost the flowering ability but exhibits rare flowering in a field. However, the terminal bud does not flower even when flowering occurs, although the lateral buds elongate and flower at the end of spring (Kinoshita et al. 2021). We consider these flowering characteristics as important indicators to characterize 'nfc'. Therefore, we evaluated flowering characteristics quantitatively using the following three criteria: the flowering date of the plant, the flowering date of the terminal bud, and the number of flowering shoots (Fig. S1). The flowering date of the plant was defined as the day when the first flower opened irrespective of whether the flowering site was a lateral bud or the terminal bud. Additionally, the flowering date of the terminal bud was examined. In plants exhibiting terminal bud flowering at first, the flowering date of the terminal bud coincides with that of the plant. Meanwhile, in plants exhibiting lateral shoot flowering at first, the flowering date of the terminal bud does not coincide with that of the plant. The terminal bud was examined to check if it developed into a flower bud or continued vegetative growth in plants whose terminal bud did not flower. The terminal bud that developed into flower buds but ceased development before the flower opening was considered an aborted bud. The number of flowering shoots was defined as the number of terminal and primary lateral shoots with flowers. This is an indicator that allows non-flowering and flowering plants to be evaluated with a continuous number. For 'kairan' and 'kairan' × 'nfc' F 1 and F 2 , only the flowering date of the plant was examined. The data for 'nfc' and 'T15' cultivated in 2018 and 2019 were obtained from a previous study (Kinoshita et al. 2021).

Library preparation and sequencing
For whole genome sequencing, the genomic DNA was extracted from leaves using the DNeasy Plant Mini Kit (Qiagen, Valencia, USA) or the method described by Mizuno et al. (2020). The DNA samples of 'nfc', 'W1', 'T15'(#43), and 'kairan' were extracted from the same clone individuals that were used for crossing. 'nfc' × 'T15' F 2 comprised the progeny of six F 1 individuals (n = 125, 82, 68, 68, 60, and 8), and the distribution of flowering dates differed slightly depending on the F 1 parent. Therefore, from each F 2 population derived from five F 1 individuals (n = 125, 82, 68, 68, 60), except for one F 2 population (n = 8), 10% of plants (n = 13, 8, 7, 7, 6) with the earliest flowering date were selected for Early-flowering bulk and 10% of plants with the latest flowering date including non-flowering plants were selected for Late-flowering bulk. In 'nfc' × 'W1' F 2 , 15% of plants (n = 16) with the earliest flowering date were selected for Early-flowering bulk and 15% of plants with the latest flowering date including non-flowering plants were selected for Late-flowering bulk. Chinese kale does not require vernalization for flowering because of the loss of BoFLC2 (Tang et al. 2021). Preliminary experiments revealed that 'kairan' used in this study lacks BoFLC2, while 'nfc', 'T15', and 'W1' harbors BoFLC2. Additionally, the expression level of BoFLC2 was comparable between 'nfc' and 'T15'. The genotypes for BoFLC2 in 'kairan' × 'nfc' F 2 should segregate into three types. Hence, each genotype should be considered separately to prevent the detection of BoFLC2 as a QTL in QTL-seq analysis. We developed a codominant DNA marker for BoFLC2 to distinguish the genotypes of 'kairan' × 'nfc' F 2 individuals using multiplex polymerase chain reaction (PCR) (Fig. S2). Consequently, 'kairan' × 'nfc' F 2 individuals were divided into 'nfc' genotype with homozygous BoFLC2 (n = 93), 'kairan' genotype without BoFLC2 (n = 77), or a heterozygous genotype (n = 174). From each BoFLC2 genotype, 20% of plants (n = 18, 16, 34) with the earliest flowering date were selected for Early-flowering bulk and 20% of plants with the latest flowering date were selected for Late-flowering bulk. Before library preparation, the genomic DNA of the F 2 plants selected for each bulk was quantified using a Qubit 2.0 Fluorimeter (Invitrogen, Carlsbad, CA, USA) and was mixed in equal amounts. Library preparation was performed using the TruSeq DNA PCR-Free Sample Prep Kit (Illumina, San Diego, CA, USA). The libraries of 'nfc' sample were sequenced using Illumina NextSeq 500 in a 150 bp paired-end mode. Meanwhile, the libraries of all other samples except 'nfc' were sequenced using Illumina NovaSeq 6000 in a 150 bp paired-end mode (Table S1).

QTL-seq analysis for the flowering date
The data in FASTQ files were filtered and trimmed using Trimmomatic (Bolger et al. 2014) to obtain a quality score of at least 20. The remaining 150 bp paired-end reads were aligned to the reference genomes of B. oleracea TO1000  using Burrows-Wheeler Aligner (BWA) with the default parameters ). Next, the SAM files were converted into BAM files, which were sorted using SAM tools ). Variant calling including both single nucleotide polymorphisms (SNPs) and insertions/deletions (indels) was performed using Genome Analysis Toolkit (GATK) (version 4) (McKenna et al. 2010). The genomic variant call format (GVCF) files of both parents and both bulks from the same F 2 population were merged using GenomicsDBImport function and GenotypeGVCFs function from GATK. The final merged variant call format (VCF) file was converted into a tab-delimited table using the VariantsToTable function from GATK.
QTL-seq analyses were performed following the methods of a previous study (Abe et al. 2012;Takagi et al. 2013;Itoh et al. 2019). Polymorphisms containing both SNPs and indels were filtered under the following three conditions; (1) genotypes (GT) of both parents were homozygous; (2) presence of genotypic (GT) polymorphism between parents; (3) read depth (DP) > 9 for both parents and both bulks. By dividing the allele depth (AD) of both bulks and comparing it to the genotypes of the TO1000 reference genome and the parental genome, we calculated the AD of 'nfc' type and the other parental type for each bulk. The frequency of 'nfc' type AD was calculated as the SNP-index of each bulk. The SNP-index in this analysis was calculated not only for SNPs but also for indels. The average ∆SNP-index (SNP-index. Late-bulk-SNP-index.Early-bulk) was calculated using the R package 'QTLseqR' (Mansfeld and Grumet 2018). Polymorphisms with SNP-index of both bulks less than 0.3 were excluded because they were considered false positive polymorphisms. Additionally, the window size was set to 10 Mb for 'nfc' × 'T15' F 2 and 2 Mb for 'nfc' × 'W1' F 2 or 'kairan' × 'nfc' F 2 to account for the number of polymorphisms. Simultaneously, 95%, 99%, and 99.9% confidence intervals were calculated for 10,000 simulations.

Validation of the genomic region associated with the flowering date using QTL analysis
The significant ∆SNP-index peak for the flowering date identified using QTL-seq was confirmed using QTL analysis with the 'nfc' × 'W1' F 2 population. Cleaved amplified polymorphic sequence (CAPS) markers were designed for candidate QTL regions on chromosome 9 (Chr. 9) using 'nfc' and 'W1' genome sequences mapped to the TO1000 reference genome (Table S2). Genomic DNA of all 'nfc' × 'W1' Page 5 of 23 41 F 2 individuals cultivated in 2020 (n = 109) was extracted from leaves, following the methods described by Mizuno et al. (2020). Specific regions were amplified with primers of each CAPS marker using KOD One ® PCR Master Mix (Toyobo Co., Ltd., Osaka, Japan), following the manufacturer's instructions. The PCR conditions were as follows: 30 cycles of 98 °C for 10 s, 65 °C for 5 s, and 68 °C for 1 s. The PCR products (1 µL) were digested with restriction enzyme (0.5-1.5 U) for at least 3 h at 37 °C and analyzed using 1% agarose gel at 100 V for 30 min. The genotype was determined based on the lengths of the detected bands.
Days to flowering of the plant, days to flowering of the terminal bud, and the number of flowering shoots were used for phenotypic characterization. Days to flowering were calculated based on the flowering date of the earliest flowering individual (for example, 1 day for the earliest flowering individual and 6 days for the individual that flowered 5 days later). Non-flowering individuals were not included in QTL analysis for days to flowering. The number of flowering shoots of non-flowering individuals was considered as zero in QTL analysis for flowering shoots.
QTL analysis was performed with R/qtl (Broman et al. 2003). The Haldane map function was used to estimate genetic distances and map order. The logarithm of odds (LOD) scores were estimated using both interval mapping (IM) and composite interval mapping (CIM). In IM analysis, the LOD scores were estimated using nonparametric interval mapping (scanone function; model = 'np') and the EM algorithm. Meanwhile, in CIM analysis, the LOD scores were estimated using Haley-Knott regression. To identify the presence of a significant QTL, 95%, 99%, and 99.9% significance thresholds were determined using a permutation test with 10,000 permutations in IM analysis and 1000 permutations in CIM analysis.

Fine mapping of QTL region using the QTL analysis
To narrow down the identified QTL regulating the nonflowering trait, QTL analysis was performed using 'kairan' × 'nfc' F 2 population cultivated in 2021. CAPS markers were designed for candidate QTL regions on Chr. 9 using 'nfc' and 'kairan' genome sequences mapped to the TO1000 reference genome (Table S2). DNA purificationfree PCR was performed using a previously reported method ) with some modifications. Approximately 1 cm 2 of leaf sampled from each seedling was rubbed onto filter papers that were air-dried after submerging in the soaking buffer. A disk with leaf samples and filter papers punched with a 1 mm diameter biopsy punch (KAI Corporation and KAI Industries Co., Ltd, Tokyo) was used as a template. Specific regions were amplified with primers of each CAPS marker using KOD One ® PCR Master Mix (Toyobo Co., Ltd., Osaka, Japan) with 4% Tween 20 (Wako Pure Chemical Industries, Ltd., Osaka, Japan). PCR and restriction digestion were performed as mentioned above. QTL analysis using only recombinant plants was associated with a high power of detection and fine mapping of the loci (Topcu et al. 2021). Therefore, genotypes of the CAPS markers P350-351 (Chr. 9: 46,064,092) and P407-408 (Chr. 9: 53,311,640) designed at both ends of the candidate region for all seedlings were determined. Only recombinant plants with different genotypes of the two CAPS markers were transplanted into the open field and genotyped using the additional seven CAPS markers. Subsequently, the flowering dates of these plants were examined. For phenotypic characterization, days to flowering of the plant calculated based on the flowering date of the earliest flowered individual was used. QTL analysis was performed as described above.

Expression analysis
For gene expression analysis, 'nfc' (n = 4), 'T15' (plants cultivated from seed, n = 3), and 'W1' (n = 3), which were transplanted into the open field on October 15, 2021, were used. The upper leaves and the maximum leaves were sampled 15-30 min before sunset at the following time points: one week after planting (October), in the middle of each month from November to March, early April, and mid-April. Meanwhile, the plants grew gradually and encountered low temperatures during winter. The upper leaf was defined as the most recently developed leaf with a length of approximately 2 cm. As the head formed after December, new leaves were sampled by pushing aside the surrounding leaves. The maximum leaf was defined as the intact leaf with the longest length. Leaf disks of the maximum leaf were sampled with a 6 mm diameter biopsy punch (KAI Corporation and KAI Industries Co., Ltd, Tokyo).
Total RNA was extracted from the leaf samples using Sepasol RNA I Super G (Nacalai Tesque, Inc., Kyoto, Japan), purified using an Econospin ™ for RNA (Epoch Life Sciences, Missouri, TX, USA), and reverse transcribed into complementary DNA (cDNA) using ReverTra Ace ® (Toyobo Co., Ltd., Osaka, Japan) and RNase Inhibitor Recombinant (Toyobo Co., Ltd., Osaka, Japan), following the manufacturer's instructions. Subsequently, 1 μL of tenfold diluted cDNA was used as a template for quantitative real-time PCR (qPCR). qPCR analysis was performed using a THUNDERBIRD ® SYBR ® qPCR Mix kit (Toyobo Co., Ltd.), following the manufacturer's instructions. PCR was performed using a LightCycler ® 480 system (Roche Diagnostics K.K., Tokyo, Japan) under the following conditions: 95 °C for 5 min, followed by 40 cycles of 95 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s. Amplification of singletarget products was evaluated using a melting curve. The expression levels of target genes were normalized to the expression level of the housekeeping gene BoActin. Gene 41 Page 6 of 23 expression was quantified using the comparative Ct method (2 −∆Ct method). The experiments were repeated two times as technical replicates. All qPCR primers used in this study are listed in Table S3. The amplification efficiency of each primer set that was calculated using a standard curve of serially diluted samples (five-fold dilution), was 100 ± 5%.

RNA-seq analysis
RNA-seq analysis of leaves was performed using 'nfc' (n = 3) and 'T15' (plant cultivated from seed, n = 3). Total RNA extracted from the upper leaf sampled in mid-November as above was used for RNA-seq analysis. RNA-seq libraries were prepared with NEBNext Ultra Directional RNA LP kit (New England Biolabs Ltd., UK) and subjected to 150 bp pair-end sequencing using an Illumina NovaSeq 6000 (Table S4). High-quality clean data were obtained by discarding the paired reads containing adapter sequences, more than 10% uncertain nucleotides, and more than 50% low-quality nucleotides (base quality less than 5), following the standard protocol provided by Novogene (Novogene Bioinformatics Technology Co. Ltd., Beijing, China). The reads were mapped against the TO1000 reference genome, and the expression levels were calculated using RSEM (Li and Dewey 2011) with Bow-tie2 alignment program. Differential expression analysis was performed using edgeR (ver 3.38.4) with glmQLFTest (Robinson et al. 2010) in R (R Core Team 2022).
RNA-seq analysis of shoot apical meristem (SAM) was performed using 'nfc' and 'T15' plants with a stem diameter of 7.5 mm or more which were cultivated in the plant growth room maintained at 22 ± 2 °C and 16 h/8 h (light/dark) photoperiods with a light intensity of 120 µmol m −2 s −1 . Due to the small amount of SAM, apices from several plants (12 plants for 'nfc' and 14 plants for 'T15') were collected, pooled and used for RNA extraction and RNA-seq (n = 1). SAM containing two or three leaf primordia was sampled under the microscope and immediately immersed in RNAlater (Invitrogen, Camarillo, USA). After centrifugation at 3000 × g for 2 min, the RNAlater was removed and the samples were homogenized using a BioMasher (Funakoshi Co., Ltd., Tokyo, Japan). Total RNA was extracted using RNeasy Plant Mini kit (Qiagen, Valencia, USA), following the manufacturer's instructions. RNA-seq libraries were prepared with TruSeq stranded mRNA Library (Illumina Inc., San Diego, USA) and subjected to 100 bp pair-end sequencing using an Illumina NovaSeq 6000 (Table S4). Reads were filtered using TRIMMOMATIC (Bolger et al. 2014) with the following options: SLIDINGWINDOW:4:15, LEAD-ING:20, TRAILING:20, MINLEN:30. Read mapping and calculation of the expression levels were performed in the same way as above. Genes with TPM > 1 in either 'nfc' or 'T15' and log 2 |Fold-Change| > 1 were considered as differentially expressed genes.

Sanger sequencing analysis
Total RNA extraction and RT reactions were performed as described above. cDNA was amplified using a gene-specific primer with KOD One ® PCR Master Mix (Toyobo Co., Ltd.), following the manufacturer's instructions, and sequences were determined using Sanger sequencing.

Phenotypic evaluation of flowering characteristics in F 2 populations
Flowering characteristics were quantified based on the following three phenotypes: the flowering date of the plant, the flowering date of the terminal bud, and the number of flowering shoots (Fig. S1). The non-flowering mutant 'nfc' rarely flowered (Table 1). In the 'nfc' flowering plants, only the lateral shoots generated from the base of the main stem flowered, while the terminal bud continued to grow vegetatively. The number of flowering shoots in 'nfc' ranged from 0 to 10. The flowering rates of 'T15' (the parental cabbage line of 'nfc'), a cabbage cultivar 'W1', and a Chinese kale cultivar 'kairan' were 100%. Additionally, the flowering rates of the terminal bud were 100% in these cultivars. 'kairan' flowered from November to  December owing to its non-vernalization requirement due to the loss of BoFLC2. The number of flowering shoots in 'T15' and 'W1' were 10-43 and 18-32, respectively, which were significantly higher than those in 'nfc'. In the 'nfc' × 'T15' population, the flowering rate of the plant of the F 1 plants was 97%, while the median flowering date of the plant was May 1 (Fig. 1a; Table 1). The flowering rate of the terminal bud was 0%, which was similar to that of 'nfc' (Fig. 1b; Table 1). The average number of flowering shoots was 7.9 ( Fig. 1c; Table 1). The flowering characteristics of 411 F 2 plants derived from six F 1 plants were examined. The flowering dates of the plant of F 2 plants were widely distributed, ranging from April 9 to non-flowering ( Fig. 1d; Table 1). The overall flowering rate of the plant was 97%, and only 12 plants were non-flowering. The flowering rate of the terminal bud was 30%. Meanwhile, the flowering dates of the terminal bud were distributed from April 11 to vegetative growth with a bias toward 'nfc' (Fig. 1e; Table 1). The number of flowering shoots was widely distributed, ranging from 0 to 27 ( Fig. 1f; Table 1).
In the 'nfc' × 'W1' population, three F 1 plants were cultivated in 2019, and only one of which flowered (Table 1). The cause of the flowering/non-flowering segregation in F 1 was unknown with no apparent difference between them. F 2 seeds were harvested from one flowering plant, and 109 F 2 plants were cultivated in 2020. The flowering dates of the plant of F 2 plants were widely distributed, ranging from March 30 to non-flowering, and tended to be biased toward 'W1' (Fig. 1g; Table 1). The flowering rate of the terminal bud was 89%. Additionally, the flowering dates of the terminal bud also tended to be biased toward 'W1' (Fig. 1h; Table 1). The number of flowering shoots was widely distributed from 0 to 61, and approximately half of the F 2 had more flowering shoots than 'W1' (Fig. 1i; Table 1).
In the 'kairan' × 'nfc' population, three F 1 plants were cultivated in 2018, and all of which flowered vigorously from the terminal bud in mid-March (Table 1). F 2 seeds were harvested from three plants, and 344 F 2 plants were cultivated in 2019. All F 2 plants flowered, and the flowering dates of the plant were widely distributed with the earliest plants flowering in early December and the latest flowering in mid-April ( Fig. 1j; Table 1). Additionally, the terminal bud flowered in all plants ( Table 1). As 'kairan' (B. oleracea var. alboglabra) and 'nfc' (B. oleracea var. capitata) are distantly related, multiple loci may be responsible for the difference in flowering time between them. This may be one of the reasons for the lack of non-flowering plants in F 2 . BoFLC2 genotyping in F 2 revealed that the flowering date markedly varied within the same BoFLC2 genotype, as well as among BoFLC2 genotypes (Fig. 1j).
Phenotypic evaluation of the three F 2 populations revealed that the phenotypic distribution in F 2 , namely heritability of the non-flowering trait of 'nfc', varied depending on the parents. Meanwhile, the distribution of the flowering dates of the plant for all F 2 populations showed continuous and wide variation. We hypothesized that the gene responsible for the non-flowering trait of 'nfc' contributes to this variation in flowering dates. Therefore, we selected plants that exhibited extreme phenotypes for the flowering date as Early-bulk and Late-bulk and subjected them to QTL-seq analysis.

Identification of a genomic region associated with flowering date using QTL-seq analysis
We performed genome sequencing of 'nfc', the other parents, and the Early-bulk and Late-bulk from each F 2 population (Table S1). After filtering, 252-336 and 117-162 million Illumina short reads were obtained for each parent and bulk, respectively. These Illumina short reads were mapped to the TO1000 reference genome, and polymorphisms were called. After filtering the merged VCF table with polymorphisms of parents and bulks, 3332, 1,620,903, and 2,563,771 polymorphisms were identified in 'nfc' × 'T15' F 2 , 'nfc' × 'W1' F 2 , and 'kairan' × 'nfc' F 2 , respectively, which were used to calculate the SNP-index.
QTL-seq analysis did not reveal any significant peaks in 'nfc' × 'T15' F 2 (Fig. 2a). The reduced power of QTL detection may be attributed to the similar genetic background of the parental lines, resulting in a low number of polymorphisms. In both 'nfc' × 'W1' F 2 and 'kairan' × 'nfc' F 2 , significant peaks were detected on Chr. 9 with the confidence interval exceeding 99%, ranging approximately from 45 to 54 Mb centered near 51 Mb (Fig. 2b, c). Therefore, the region Chr. 9: 45-54 Mb was most likely the candidate genomic region associated with the non-flowering trait of 'nfc'.

Validation of the candidate genomic region using QTL analysis
To validate the candidate genomic region identified using QTL-seq, traditional QTL analysis was performed with 'nfc' × 'W1' F 2 using 10 CAPS markers designed to span the Chr. 9: 45-54 Mb region. IM and CIM analyses identified significant QTL for days to flowering of the plant (n = 101, flowering plants), days to flowering of the terminal bud (n = 97, plants with flowering terminal bud), and the number of flowering shoots (n = 109, all plants) (Figs. 3, S3). The QTL for days to flowering of the plant exhibited a 1.5-LOD support interval (SI) extending from 8.6 to 22.2 cM (between P350-351 and P360-361 markers, corresponding physical interval 46,064,092-52,925,397 bp) in IM and exhibited a significant LOD at a 5% level extending from 12.8 to 22.7 cM (between P354-355 and P360-361 markers, corresponding physical interval 48,153,540-52,925,397 bp) in CIM (Table 2). Meanwhile, the QTL for days to flowering of the terminal bud and the number of flowering shoots revealed similar results (Fig. S3). This QTL explain 27.6%, 29.4% and 32.6% of the phenotypic variation, respectively. The values of each phenotype differed significantly according to the genotype of the closest CAPS marker to the peak LOD interval (Figs. 3d, S3c, f). Therefore, the candidate genomic region identified using QTL-seq analysis was validated, and this QTL was considered to harbor the candidate gene responsible for the non-flowering trait of 'nfc'.

Fine mapping of the QTL region using QTL analysis
To further delineate this QTL interval, traditional QTL analysis was performed with 'kairan' × 'nfc' F 2 cultivated in 2021. 'kairan' × 'nfc' F 2 seedlings (n = 1075) were genotyped with the CAPS markers P350-351 (Chr. 9: 46,064,092) and P407-408 (Chr. 9: 53,311,640), and only recombinant plants were transplanted into the field in 2021 (n = 575). Additionally, we selected only plants with a maximum leaf length of ≥ 40 cm in December as well-growing plants (n = 432). These plants were genotyped with additional seven CAPS markers. The frequency distribution of this population for the flowering date in 2021 exhibited continuous variation and was similar to that in 2019 ( Fig. 4a; Table 1). IM analysis revealed a 41 Page 10 of 23 1.5-LOD SI extending from 25.8 to 31.4 cM (between P501-502 and P505-506 markers, corresponding physical interval 50,649,301-51,474,818 bp). Meanwhile, CIM analysis revealed a significant LOD at a 5% level extending from 23.4 to 33.3 cM (between P358-359 and P505-506 markers, corresponding physical interval 50,177,696-51,474,818 bp) (Fig. 4b-d; Table. 2). This QTL explained 15.0% of the phenotypic variation. Based on the results of IM and CIM analyses, the QTL region harboring the candidate gene for the nonflowering trait of 'nfc' was narrowed down to the region 50,177,696-51,474,818 bp of Chr. 9. Additionally, the genotype of the closest CAPS marker to the peak LOD, as well as the genotype of BoFLC2 marker, had a significant effect on days to flowering of 'kairan' × 'nfc' F 2 population ( Fig. 4e-f). The genotype of BoFLC2 significantly interacted with that of the closest CAPS marker to the peak LOD, suggesting that they exerted additive effects on the flowering date (Fig. 4g).

Estimation of the candidate genes by gene function prediction
T h e n a r r o w e d -d o w n c a n d i d a t e r e g i o n 50,177,696-51,474,818 bp of Chr. 9 covered 241 genes in the TO1000 reference genome. Basic local alignment search tool (BLAST) analysis of these 241 genes was performed against the Arabidopsis thaliana cDNA database. The  (Sheldon et al. 1999(Sheldon et al. , 2000Amasino 1999, 2001). CDKC;2, a component of positive transcription elongation factor b (P-TEFb), influences global RNA polymerase II phosphorylation and transcription of COOLAIR, which is a long noncoding RNA transcribed from the antisense direction in the FLC locus (Cui et al. 2007;Liu et al. 2010;Hornyik et al. 2010;Wang et al. 2014). HAM1/2 influence global histone 4 acetylation and H4K5ace at FLC and MADS AFFECTING FLOWERING 3/4 (MAF3/4) (Latrasse et al. 2008;Xiao et al. 2013). Then, we compared the coding sequences (CDS) and amino acid sequences of these four genes between 'T15' and 'nfc'. LOC106318712 (FLC homolog) had three isoforms with different start positions for the eighth exon, one of which also had a different termination codon position (Fig.  S4a). The sequences of all four genes of 'nfc', including isoforms, were identical to those of 'T15' (Fig. S5).

RNA-seq analysis and identification of the candidate genes
Although QTL-seq analysis did not reveal any peaks in 'nfc' × 'T15' F 2 (Fig. 2a) and there were no differences in the amino acid sequences of the four flowering time-related gene homologs within the candidate region (Fig. S5), we hypothesized the presence of at least some differences in gene expression between 'nfc' and 'T15' that contribute to the different flowering characteristics. The transcriptomes in upper leaves in November of 'nfc' and 'T15' were comparatively analyzed. We detected 535 upregulated genes [false discovery rate (FDR) < 0.05, log 2 Fold-change (FC) > 1] and 267 downregulated genes (FDR < 0.05, log 2 FC < 1) in 'nfc' (Fig. 5a). Among the differentially expressed genes (DEGs), 19 were homologs of flowering time-related genes in Arabidopsis (Fig. 5b, blue and red heatmaps indicate genes with negative and positive effects on flowering time, respectively). Meanwhile, among the DEGs with negative effects on flowering time, LOC106318713 (FLC homolog) and LOC106318712 (FLC homolog) had the first and eighth smallest FDRs, respectively, and were upregulated in 'nfc'. According to these results, the flowering time-related gene homologs and DEGs within the candidate region associated with the flowering date were LOC106318712 and LOC106318713 (Fig. 5c). The expression levels of the other two flowering time-related gene homologs located within the candidate region were comparable between 'nfc' and 'T15' (Fig. S6a, LOC106318968, 'T15' transcripts per million [TPM] vs. 'nfc' TPM = 19.0 vs. 18.5, FDR = 1.0; LOC106315498, 'T15' TPM vs. 'nfc' TPM = 11.5 vs. 13.2, FDR = 0.8). Additionally, the transcriptomes in SAM of 'nfc' and 'T15' cultivated in the growth room maintained at 22 ± 2 °C were comparatively analyzed. LOC106318713 and LOC106318712 were detected as the negative-effective flowering time-related genes with the largest log 2 FC, while LOC106318968 and LOC106315498 had comparable relative expression levels (Figs. 5d, S6b, LOC106318968, 'T15' TPM vs. 'nfc' TPM = 28.7 vs. 28.8; LOC106315498, 'T15' TPM vs. 'nfc' TPM = 43.4 vs. 47.7). Therefore, these two FLC homologs (LOC106318712 and LOC106318713) were considered as the candidate genes responsible for the nonflowering trait of 'nfc'.
The two FLC homologs are located close to each other (17.5 kb apart) in the TO1000 reference genome. BoFLC1 tandem duplication event is thought to have occurred independently after a whole genome triplication in B. oleracea, and a corresponding tandem duplication also exists in Brassica napus (Lysak et al. 2005;Schiessl et al. 2019;Calderwood et al. 2021). Based on the names of the corresponding B. napus FLC homologs, LOC106318712 and LOC106318713 were named BoFLC1a and BoFLC1b, respectively (Fig. S7a). LOC106318714 (annotated as protein DOWNSTREAM OF FLC-like) was located in the 17.5 kb region between BoFLC1a and BoFLC1b. The gene was not a homolog of flowering time-related genes and was barely expressed in both leaves and SAMs of 'nfc' and 'T15' (Fig. S8). BoFLC1a had three isoforms in the Table 2 Results of QTL analysis for days to flowering of the plant in 'nfc' × 'W1' F 2 population and 'kairan' × 'nfc' F 2 population a QTL position in interval mapping (IM) was defined as a 1.5-logarithm of odds (LOD) interval. QTL position in composite interval mapping (CIM) was defined as the region above the 5% significance level b The marker interval was defined as the region flanked by neighboring cleaved amplified polymorphic sequence ( sequence analysis (from shortest to longest: BoFLC1a-1, BoFLC1a-2, BoFLC1a-3) (Fig. S4a). The ratios of these isoforms (BoFLC1a-1: BoFLC1a-2: BoFLC1a-3) calculated from the RNA-seq short reads were approximately 3:6:1 for the upper leaves of 'T15' and 'nfc', approximately 7:3:0 for the SAM of 'T15', and approximately 6:4:0 for the SAM of 'nfc' (Fig. S4b), with no significant differences between 'nfc' and 'T15'.

Estimation of protein function of three BoFLC1a isoforms and BoFLC1b
Since BoFLC1a and BoFLC1b were identified as the candidate genes, we performed the protein function prediction to estimate whether both of the two genes contribute to floral repression or only one of them contributes to floral repression. The amino acid sequence of BoFLC1b of  Table S2. formation, respectively, while the C domain is a variable primary sequence with few conserved structural motifs (Lai et al. 2019). Conformational prediction by Alpha-fold2 (Akdel et al. 2022) showed that the conformations of the monomer and homodimer of AtFLC, BoFLC1a-1, and BoFLC1b were similar (Fig. S10). On the other hand, BoFLC1a-2 and BoFLC1a-3 retained the MADS-box domain, but the K-box domain of the homodimer was very different due to the extra 126 or 134 aa (Fig. S10). Therefore, BoFLC1a-1 is likely to function as a floral repressor like BoFLC1b, but it was unclear whether BoFLC1a-2 and BoFLC1a-3 retained the floral repressor function.

Expression analysis of BoFLCs and BoFT
The relative expression levels of BoFLCs and BoFT in 'T15', 'nfc', and 'W1' in the upper and maximum leaves sampled over time from October to April were quantified (Fig. 6).
In order to accurately distinguish between BoFLC1a and BoFLC1b, primers were designed based on the sequences of the de novo assembled transcripts of BoFLC1a and BoFLC1b using the RNA-seq reads of 'nfc' and 'T15' (Fig. S7b). No contigs belonging to the BoFLC5 clade were detected in the phylogenetic tree, suggesting that BoFLC5 is a pseudogene in 'nfc' and 'T15' (Fig. S7a). The expression levels of BoFLC2 and BoFLC3 exhibited similar trends in all three cultivars and were downregulated upon exposure to low temperatures during winter (Fig. 6). Meanwhile, the expression patterns of BoFLC1a and BoFLC1b differed among the cultivars. In the upper and maximum leaves of 'T15' and 'W1', the expression levels of BoFLC1a and BoFLC1b decreased gradually from October to January and remained low thereafter. However, the expression level of BoFLC1a in the upper leaves of 'nfc' gradually decreased from October to January and subsequently increased. The expression level of BoFLC1a in April was higher than that in October. In mid-April, the BoFLC1a expression level in 'nfc' was approximately 10 times higher than that in 'T15'. The BoFLC1b expression levels in the upper leaves of Fig. 6 Expression analysis of BoFLCs and BoFT using quantitative real-time PCR (qPCR) in 'T15' (n = 3), 'nfc' (n = 4), and 'W1' (n = 3). Different letters or asterisks indicate significant differences at a 5% significance level using the Tukey-Kramer test. An asterisk indicates that the expression level in the cultivar is significantly higher than that in the other two cultivars. The upper and lower panels for each gene indicate the upper or maximum leaf samples, respectively. The x-axis indicates the month of sampling (from October to March, leaves were sampled in the middle of each month); 4e and 4m indicate early April and mid-April, respectively 41 Page 16 of 23 'nfc' increased gradually after October. In mid-April, the expression level of BoFLC1b in 'nfc' was approximately 20 times higher than that in 'T15'. The expression levels of BoFLC1a and BoFLC1b in the maximum leaves of 'nfc' decreased gradually, similar to those in 'T15' and 'W1', but were significantly higher than those in the maximum leaves of 'T15' and 'W1'. 'T15' and 'W1' that flowered on April 10-11 and April 13, respectively, were used for this expression analysis. Cauline leaves with a length of approximately 2 cm were sampled as the upper leaves in early April and mid-April. 'nfc' used for the expression analysis did not flower. In 'T15' and 'W1', the expression level of BoFT began to increase in early April and further increased in mid-April in both the upper and maximum leaves. However, the expression level of BoFT in the maximum leaves of 'nfc' slightly increased in April but remained at a low level, whereas that in the upper leaves did not increase. This suggested that BoFT expression was repressed in 'nfc' owing to the upregulation of both BoFLC1a and BoFLC1b. The similar expression patterns of BoFLC1a and BoFLC1b in 'nfc' suggest that a cis-regulatory mutation in the vicinity of these genes may have contributed to their upregulation.

Search for a causal mutation responsible for expression changes of BoFLC1a and BoFLC1b in 'nfc'
Prediction of protein function and expression analysis suggested that upregulation of both BoFLC1a and BoFLC1b in 'nfc' contributes to the non-flowering trait and the presence of a cis-regulatory mutation affecting the expression of both genes. From a comparative analysis of cabbage inbred lines with different flowering times in a previous study, it was reported that early-flowering cabbage lines harbor the BoFLC1 allele with a 67 bp insertion in the second intron, and BoFLC1 expression in early-flowering lines was lower than in late-flowering lines (Abuyusuf et al. 2019). Sanger sequencing suggested that 'nfc', 'T15', and 'W1' harbored the late-flowering type BoFLC1 allele without a 67 bp insertion in the second intron (data not shown). This suggested that a novel BoFLC1 allele with a cis-regulatory mutation may be associated with the non-flowering trait of 'nfc'. Therefore, we attempted to identify the causal mutation in 'nfc' using several strategies. Comparative analysis of a genomic sequence from 100 kb upstream of the transcription start site (TSS) of BoFLC1a to 100 kb downstream of the termination codon of BoFLC1b using short reads revealed no polymorphism with an absolute value of ∆SNP-index (SNP-index. nfc-SNP-index.T15) greater than 0.75 (Fig. S11). PCR of the genomic regions of 7290 bp from 1752 bp upstream of the TSS of BoFLC1a to 481 bp downstream of the termination codon and 6495 bp from 1490 bp upstream of the TSS of BoFLC1b to 780 bp downstream of the termination codon revealed no differences in band length among 'nfc', 'T15', 'W1', and 'kairan' (Fig. S12). Therefore, a causal mutation may be present in a region more than 100 kb distant from BoFLC1a and BoFLC1b or may be a few hundred nucleotides indel mutation that is not large enough to be mapped to the reference genome but not small enough to be detected using genomic PCR in this study, or may be an epigenetic mutation.

Discussion
'nfc' is a non-flowering natural cabbage mutant, which was discovered from the breeding line 'T15' with normal flowering characteristics (Kinoshita et al. 2021). Flowering and seed set are important events for the survival of offspring in nature or maintaining seed production efficiency of crops cultivated from seed. Therefore, plants that cannot flower are eliminated in nature or during the crossbreeding process. Additionally, rapid-cycle lines tend to be preferred as research materials for conducting experiments easily and quickly. Therefore, the non-flowering trait of 'nfc' (namely the trait that does not flower even under floral inductive conditions) is rare and interesting. The elucidation of the non-flowering mechanism of 'nfc' may lead to the discovery of a novel floral repression mechanism or a novel factor necessary for cabbage to flower. In this study, 'nfc' was induced to flower using the grafting floral induction method (Motoki et al. 2019;Motoki et al. 2022;Motoki et al. 2023), and segregating populations were generated, which enabled the identification of QTL responsible for the non-flowering trait.

Identification of the candidate genes responsible for the non-flowering trait of 'nfc'
The QTL was detected using QTL-seq and QTL analysis, while DEGs were identified using RNA-seq analysis (Figs. 2,3,4,5). These analyses revealed tandem duplicated BoFLC1 genes near 51 Mb on Chr. 9 as the candidate genes responsible for the non-flowering trait of 'nfc' (Fig. 5c). The tandem duplicated BoFLC1 genes were named BoFLC1a and BoFLC1b. BoFLC1b of 'T15'/'nfc' was considered to function as a floral repressor and BoFLC1a-1 among the three BoFLC1a isoforms was inferred to function as a floral repressor based on its amino acid sequence and predicted protein conformation (Figs. S9, S10). In addition, the expression levels of both BoFLC1a and BoFLC1b were gradually decreased Page 17 of 23 41 in 'T15', whereas both were upregulated in 'nfc' (Fig. 6). Therefore, both BoFLC1a and BoFLC1b were thought to contribute to the non-flowering trait of 'nfc'. Considering that FLC acts on floral repression in a dosage-dependent manner (Calderwood et al. 2021), we speculate that BoFLC1b has a greater effect on floral repression because the TPM of BoFLC1b in leaves and SAM of 'nfc' were about 8 and 11 times higher than that of BoFLC1a, respectively (Fig. 5b, d).
In Arabidopsis, FLC expression is activated by FRIGIDA (Johanson et al. 2000;Geraldo et al. 2009). On the other hand, FLC expression is repressed by genes in autonomous pathways [FCA, FPA, FY, FLOWERING LOCUS D (FLD), FLK, LUMINIDEPENDENS, and FVE] whose activities appear to be largely independent of the environment (Simpson and Dean 2002). In addition, FLC expression is also repressed under low-temperature conditions by genes in the vernalization pathway (VIN3, VRN1, VRN2, VRN5) and the non-coding RNA COOLAIR (Whittaker and Dean 2017). FLC expression is epigenetically silenced by elevated histone 3 at Lys 9 trimethylation (H3K9me3) and H3K27me3 levels according to vernalization (Bastow et al. 2004;Sung and Amasino 2004). In Arabidopsis, the flowering behavior of each ecotype is determined by the pre-vernalization expression level and the reduction ratio during vernalization and these differences are caused by genetic variation in FLC (Shindo et al. 2005(Shindo et al. , 2006Stinchcombe et al. 2005;. For example, in Lov-1 and Var2-6, two of the most vernalization-requiring accessions in Arabidopsis, four non-coding SNPs in the 5' region or one non-coding SNP in the first intron of FLC are responsible for decreasing the silencing ratio or increasing the expression, respectively (Coustham et al. 2012;Li et al. 2015). In B. rapa, a large insertion of approximately 5 kb into the first intron resulted in weak repression of BrFLC2 and BrFLC3 during cold exposure (Kitamoto et al. 2014). In B. oleracea, a 67 bp insertion into the second intron of BoFLC1b significantly reduced the pre-vernalization expression level (Abuyusuf et al. 2019).
In this study, it was considered that a single cis-regulatory mutation in 'nfc' is responsible for the expression changes of the two genes because 'nfc' is a natural mutant arising from 'T15', BoFLC1a and BoFLC1b are located tandemly (Fig. S12a), and the expression patterns of the two genes are similar (Fig. 6). However, we could not find the causal mutation in 'nfc' in this study (Figs. S11, S12). Future studies must focus on the de novo assembly of 'nfc' and 'T15' genomes using long reads and their comparative analysis to identify the causal mutation in 'nfc'. In the previous study, it was suggested that epigenetic regulation is involved in the flowering characteristics of 'nfc' because 'nfc' flowered spontaneously at a low frequency and it was not due to tissue chimerism or differences in environmental conditions (Kinoshita et al. 2021). Therefore, we need to search for the cis-regulatory mutation in 'nfc' also taking into account the possibility that it is an epigenetic mutation. The expression levels of BoFLC1a and BoFLC1b in 'nfc' were significantly higher before vernalization and did not decrease under floral induction conditions (Fig. 6). Especially, there was a difference in the expression levels between 'T15' and 'nfc' even in the SAM which was not exposed to low temperatures at all (Fig. 5d). Therefore, we speculate that cis-regulatory mutations in 'nfc' prevent the repression of the expression of BoFLC1a and BoFLC1b by the genes of the autonomous and/or vernalization pathways. Additionally, the expression patterns of BoFLC1a and BoFLC1b in 'T15' and 'W1' were similar in the upper and maximum leaves, whereas their patterns in 'nfc' were different between the upper and maximum leaves (Fig. 6). The causal cis-regulatory mutation in 'nfc' may be responsible for the differences in expression patterns between the upper and maximum leaves. Future identification of the causal cisregulatory mutation in 'nfc' will provide new insights into the regulation of FLC expression.

Estimation of gene networks downstream of BoFLC1a and BoFLC1b
FLC plays a central role in the floral regulatory network, and in Arabidopsis, it has been reported that more than 500 genes were differentially expressed in flc mutant leaves compared to that of wild type and that more than 300 genomewide binding sites of FLC were present by ChIP-seq analysis (Mateos et al. 2015). Therefore, flowering time-related DEGs other than BoFLC1a and BoFLC1b identified by RNA-seq analysis may be directly or indirectly regulated by BoFLC1a and BoFLC1b (Fig. 5b, d). RNA-seq analysis of 'nfc' and 'T15' leaves revealed that two BRANCHED1 (BRC1) homologs and three TEMPRANILLO1/2 (TEM1/2) homologs, as well as BoFLC1a and BoFLC1b, were upregulated in 'nfc' (Fig. 5b). One additional BRC1 homolog and three additional TEM1/2 homologs were expressed at higher level in 'nfc' than in 'T15', although not significantly (data not shown). In Arabidopsis, BRC1, which encodes a TCP transcription factor, is an integrator gene that determines the activity of axillary buds (Aguilar-Martínez et al. 2007;González-Grandío et al. 2013). BRC1 protein interacts with FT and TWIN SISTER OF FT (TSF) to repress the early floral transition of axillary buds, and the ectopic expression of BRC1 in the terminal bud delays the floral transition in the terminal shoot (Niwa et al. 2013). The upregulated expression of BRC1 homologs in the upper leaves of 'nfc' may be associated with the continued vegetative growth of the terminal bud. In winter oilseed rape, warming during winter delays flowering through an ABA-related bud dormancy, 41 Page 18 of 23 which may be regulated by BRC1 (Lu et al. 2022). Differences in flowering response to warming during winter among cultivars were significantly associated with the haplotypes of BnaFLC.A03b and BnaFLC.C02 (Lu et al. 2022). TEM1 and TEM2, which belong to the RAV (related to ABI3/ VP1) subfamily of transcription factors, have a redundant function to repress the expression of FT, GA4 biosynthesis genes GA3-oxidase1/2 (GA3ox1/2), and miRNA172 (Castillejo and Pelaz 2008;Osnato et al. 2012;Aguilar-Jaramillo et al. 2019). Previous studies have suggested that TEM1 and TEM2 are the direct targets of the FLC-SHORT VEGETA-TIVE PHASE (SVP) complex (Mateos et al. 2015). Each of the three BRC1 homologs and the six TEM1/2 homologs identified using RNA-seq may be coordinately regulated by a common transcription factor, which may be BoFLC1a, BoFLC1b, or their downstream factors. BoFT expression or BoFT activity in 'nfc' may be additively repressed not only by the upregulated expression of BoFLC1a and BoFLC1b but also by the upregulated expression of BRC1 and TEM1/2 homologs. RNA-seq analysis of SAM revealed that a SQUAMOSA PROMOTER BINDING PROTEIN-LIKE 15 (SPL15) homolog was downregulated in 'nfc' (Fig. 5d). SPLs, which are involved in the phase transition from juvenile to adult phase, are directly repressed by microRNA156 (miR156) and FLC/PERPETUAL FLOWERING 1 (PEP1, FLC homolog in Arabis alpina) and activate miR172 and floral meristem identity genes such as LEAFY, APETALA1, FRUITFULL independently of FT (Wang et al. 2009;Bergonzi et al. 2013;Zhou et al. 2013;Hyun et al. 2016Hyun et al. , 2019Xu et al. 2016). Therefore, the high expression of BoFLC1a and BoFLC1b in 'nfc' may repress the SPL15 expression and the transition to the adult phase at the SAM.

Role of BoFLC1a and BoFLC1b in B. oleracea
The following B. oleracea FLC homologs have been previously reported: BoFLC1, BoFLC2 (also known as BoFLC4), BoFLC3, and BoFLC5 (Lin et al. 2005;Okazaki et al. 2007;Razi et al. 2008). However, in the B. oleracea reference genome sequences of TO1000 and JZS (also known as 02-12) (Liu et al. 2014;Parkin et al. 2014;Cai et al. 2020), two tandem duplicated BoFLC1 genes are located on Chr. 9. The transcripts derived from each locus are registered in the National Center for Biotechnological Information (NCBI) Reference Sequence Database of TO1000. In B. napus (AC genome), BnaFLC.C9a and BnaFLC.C9b correspond to tandem duplicated BoFLC1 genes. However, the corresponding gene does not exhibit tandem duplication in B. rapa (A genome) and B. nigra (B genome), suggesting that the tandem duplication event occurred in B. oleracea (C genome) independently after whole genome triplication (Lysak et al. 2005;Schiessl et al. 2019;Calderwood et al. 2021). However, to the best of our knowledge, BoFLC1a and BoFLC1b have not been separately analyzed. This may be due to the following reasons: (1) the two genes are difficult to distinguish due to their high sequence homology (98.6% identity of CDS between BoFLC1a-1 and BoFLC1b), (2) they are recognized as a single QTL in linkage analysis because they are only 17.5 kb apart, and (3) the expression level of BoFLC1b is higher than that of BoFLC1a as shown in this study [ Fig. 5b, BoFLC1a: BoFLC1b = 12.0:93.2 ('nfc' TPM), 2.7:32.4 ('T15' TPM)]. The results of this study indicated that both BoFLC1a and BoFLC1b were expressed in 'nfc', 'T15', and 'W1' (Figs. 5b, d, 6). The expression patterns of BoFLC1a and BoFLC1b were similar within cultivars (Fig. 6). In Arabidopsis, the expression of a small gene cluster including FLC and its flanking genes is coordinately regulated in response to vernalization or changes in genetic background (Finnegan et al. 2004(Finnegan et al. , 2005. On the other hand, a comprehensive survey of the nine B. napus FLC homologs revealed that BnaFLC. C9a and BnaFLC.C9b, corresponding to BoFLC1a and BoFLC1b, had the highest coding sequence similarity among the nine BnaFLC homologs, but their expression dynamics during vernalization were less similar (Calderwood et al. 2021). Further studies are needed to examine if the expression dynamics of BoFLC1a and BoFLC1b are coordinated or independent using several other cultivars. Furthermore, among the three isoforms of BoFLC1a, we speculated that BoFLC1a-1 functions as a floral repressor, but the function of BoFLC1a-2 and BoFLC1a-3 was not clear (Figs. S9, S10). In the future, the function of each of the three BoFLC1a isoforms should be investigated by the overexpression or knockdown experiment.
BoFLC2 has been considered to be the primary FLC homolog for determining flowering time in B. oleracea (Okazaki et al. 2007;Ridge et al. 2015;Irwin et al. 2016;Li et al. 2022). However, BoFLC1, BoFLC2, and BoFLC3 have been demonstrated to function as floral repressors (Itabashi et al. 2019;Lin et al. 2005;Irwin et al. 2016). BoFLC1 and BoFLC3 have been detected as a QTL that determine the flowering time in comparisons of lines without polymorphism in BoFLC2 (Abuyusuf et al. 2019;Lin et al. 2018). The current study suggested that BoFLC1 and BoFLC2 have comparable negative effects on the flowering date of 'kairan' × 'nfc' F 2 population and that they act additively ( Fig. 4e-g). Additionally, the contribution ratio of the QTL to days to flowering was 27.6% and 15.0% in 'nfc' × 'W1' F 2 and 'kairan' × 'nfc' F 2 , respectively (Table 2), and all non-flowering plants (7%, 8/109) in 'nfc' × 'W1' F 2 had the 'nfc' genotype for the closest CAPS marker to the peak LOD interval (Fig. 3d). Therefore, the QTL including BoFLC1a and BoFLC1b detected in this study was considered to have a strong floral repressive effect. In B. napus, total FLC expression dynamics explain the differences in Page 19 of 23 41 vernalization requirement among cultivars rather than the expression of specific FLC homologs (Calderwood et al. 2021). Therefore, it is likely that flowering time in B. oleracea is also determined by a combination of the expression dynamics of each BoFLC homolog. The floral response of B. oleracea has markedly diversified, probably resulting from the differences in the harvested organs in different varieties [leaves (cabbage), stems (kohlrabi), and flower buds (broccoli, cauliflower)] and the adaptation to the local climate. For a comprehensive understanding of the floral response of B. oleracea, we believe that the genomic polymorphisms and expression patterns of BoFLC1a/1b/2/3/5 related to the floral response need to be investigated and organized using a wide range of varieties and cultivars.

Conclusion and future prospects
The results of this study suggested that the upregulation of tandem duplicated BoFLC1a and BoFLC1b contributed to the non-flowering trait in cabbage (Fig. 7). This provides important insights into the floral induction mechanism in B. oleracea. Since flower bud development and bolting significantly reduce the yield and commercial value of cabbage, breeding cultivars with a late bolting trait is important for stable production (Jung and Müller 2009). By introducing a genomic region containing BoFLC1a and BoFLC1b from 'nfc' into commercial cultivars by using grafting and repeating selection using DNA markers, it is possible to produce unprecedentedly extremely late bolting cultivars. On the other hand, the non-flowering trait is in a trade-off relationship with seed production efficiency. For the efficient seed production of extremely late bolting cultivars with a genomic region derived from 'nfc', it is necessary to further develop a stronger floral induction method such as the grafting floral induction method and to search for a condition that reduces the expression of BoFLC1a and BoFLC1b of 'nfc'. In the future, the identification of the causal mutation in 'nfc' and the elucidation of the gene network will provide further insights into the regulatory mechanisms of flowering in B. oleracea and its application in late-bolting breeding.

Fig. 7
Schematic diagram of the non-flowering mechanism in 'nfc' revealed in this study. Green, blue, and orange bars indicate the expression levels of BoFLC1a/1b, BoFLC2/3, and BoFT, respectively (color figure online)