Two Chinese cabbage DH lines, white-flowered Y640-288 and yellow-flowered Y641-87, were used as parental lines to generate F₁ and F₂ populations for inheritance analysis and gene mapping. Additionally, three F₂ populations, (Y640-288 × SD369)-F₂, (Y640-288 × Chiifu)-F₂, (Y66-83 × R16-11)-F₂, were generated for marker validation by crossing the white-flowered DH lines Y640-288 and Y66-83 with the yellow-flowered DH lines SD369, Chiifu and R16-11. Furthermore, ten white-flowered and ten yellow-flowered DH lines (Supplementary Table 1) were also used to analyze mutations in the candidate gene. All the materials used in this study were provided by Institute of Horticulture, Henan Academy of Agricultural Sciences.

Petals from Y640-288 and Y641-87 flowers at the flowering stage were used for transmission electron microscopy (TEM) analysis, which was performed according to Ariizumi et al. (2004).

Carotenoid composition was measured by MetWare (Wuhan, China). Petals from 10 white-flowered F₂ plants were combined to form one replicate W-bulk and petals from 10 yellow-flowered F₂ plants were included in the Y-bulk. In total, three replicates were assessed. Fresh petals were frozen in liquid nitrogen and stored at −80°C until needed. The direct extraction steps were performed according to Zhou et al. (2020). The saponified extraction steps were performed according to Inbaraj et al. (2008) with some modification. The direct and saponified extracts were then analyzed using an LC-APCI-MS/MS system (UHPLC, ExionLCs AD; MS, Applied Biosystems 6500 Triple Quadrupole). The chromatographic conditions and parameters for API 6500 Q TRAP LC-MS/MS System were performed as previously reported (Liu et al., 2020; Wang et al., 2020; Zhou et al., 2020). Specific procedures for extraction, identification and quantification of carotenoids were supported in Supplementary Material 1.

Carotenoids were identified by comparing their retention times and ion pair information (Supplementary Table 2). In saponified extracts, the integrated peak area was substituted into the linear equations of standard (Sigma, St. Louis, MO, United States) curves for content calculation (Supplementary Table 3); finally, the absolute content data for the carotenoids in the actual samples were obtained. Carotenoid content (μg/g) = B^∗C/1000/D, where B is the concentration (μg/mL) obtained by substituting the integrated peak area of a carotenoid in the sample into the corresponding standard curve, C is the resuspension volume (μL), and D is the mass of the weighed sample (g) (Liu et al., 2020; Wang et al., 2020; Zhou et al., 2020).

Using the modified cetyltrimethylammonium bromide (CTAB) method, total genomic DNA was isolated from young leaves of the parents and F₂ (Y640-288 × Y641-87) plants (Liu et al., 2003). For BSA-seq, two DNA pools were constructed by mixing equal amounts of DNA from 50 white-flowered F₂ individuals (W-pool) and 50 yellow-flowered F₂ individuals (Y-pool). The Illumina HiSeq X Ten platform was used to generate 150-base paired-end reads for Y640-288, Y641-87, the W-pool and the Y-pool by Anoroad Biotech Co., Ltd. (Beijing, China). The raw data were deposited in the Sequence Read Archive (SRA) in NCBI as PRJNA682710.

The clean reads of Y640-288, Y641-87, the W-pool and the Y-pool were aligned to the B.rapa Chiifu reference genome version 1.5 using BWA software (Li and Durbin, 2010). SAMtools software package (version 1.3.1) (Li et al., 2009) was then used to call single-nucleotide polymorphism (SNP) and insertion/deletion (InDel) variants based on alignment files. To identify candidate regions associated with the white flower trait, the SNP-index and △(SNP-index) were calculated for all genomic positions in the W-pool and Y-pool. The SNP-index was estimated from the proportion of reads harboring SNPs among the entire number of reads compared to the reference genome sequence (Abe et al., 2012). Then △(SNP-index) was calculated by subtracting the SNP-index of the Y-pool from that of the W-pool (Takagi et al., 2013). We filtered out SNPs with SNP-index <0.3 or >0.8 simultaneously in the two bulked pools. Furthermore, SNPs with heterozygous genotypes in the parental lines were also excluded. A 1-Mb sliding window with a 50-kb increment was applied to slide across the genome, and Δ(SNP-index) graphs were plotted using the average Δ(SNP-index) against the positions of each sliding window. We calculated the statistical confidence intervals of Δ(SNP-index) among all SNP positions with given read depths under the null hypothesis of no major genes, and the 95% and 99% confidence intervals of the Δ(SNP-index) were then generated for each read depth according to Takagi et al. (2013).

To map the BrWF3 gene, we extracted the 70-bp upstream and downstream sequences of the selected SNP for KASP marker development. For each selected SNP, two allele-specific forward primers and one common reverse primer were designed using the Primer Premier 5.0 program (Singh et al., 1998) according to the standard KASP guidelines. The two allele-specific primers were added with the standard FAM (5′-GAAGGTGACCAAGTTCATGCT-3′) and HEX (5′-GAAGGTCGGAGTCAACGGATT-3′) tails at the 5′ end. The developed KASP markers were first validated in the two parental lines and F₁ plants for polymorphism screening. Then, the polymorphic KASP markers (Supplementary Table 4) were employed to genotype the F₂ population using 135 individuals. KASP assays were performed as described by Yang et al. (2020). The genetic linkage map was constructed using JoinMap 4.0¹ software. Recombination values were converted into genetic map distances (cM) following the Kosambi mapping function (Kosambi, 1944).

To clone the DNA and cDNA sequences of the putative candidate genes, primers (Supplementary Table 5) were designed according to the B.rapa reference genome. PCR amplification was performed in a total volume of 25 μL according to the manual supplied with Phanta Max Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China). The PCR products were sequenced by Sunya Biotech Co., Ltd. (Zhengzhou, China), and sequence alignments were performed using DNAMAN software. The complete coding sequences of candidate gene from Y641-87 and Y640-288 were submitted to GenBank under the accession numbers: MW362118 (Y641-87) and MW362119 (Y640-288), respectively.

Total RNA was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) according to the manufacturer’s instructions. First-strand cDNAs were synthesized in a 20 μL volume containing approximately 7 μg RNA and oligo (dT) primers using the TransScript One-Step gDNA Removal and cDNA Synthesis Kit (Trans, Beijing, China). qRT-PCR was performed with 2 × TB Green Premix Ex Taq II (TaKaRa, Japan) on a Roche LightCycler 480-II System (Roche Applied Sciences, Beijing, China). The Bractin gene was used as an internal control (Fujimoto et al., 2006; Takada et al., 2019). Each qRT-PCR experiment was performed in triplicate and the resultant mean value was used for qRT-PCR analysis (Zhang et al., 2013). Relative expression levels were calculated using the 2^{–Δ Δ Ct} method (Livak and Schmittgen, 2001).

The W-bulk and Y-bulk each with three replicates used for carotenoid analysis were also used for transcriptome analysis. Six cDNA libraries were constructed and sequenced on the Illumina HiSeq X Ten platform at Metware Biotech Co., Ltd. (Wuhan, China). Raw reads were filtered by removing low-quality reads and reads containing the adapter or ploy-N using in-house Perl scripts available from Metware Biotech Co., Ltd. (Wuhan, China). The clean reads were aligned to the B.rapa V1.5 reference genome using HISAT2 software (Kim et al., 2015). Differentially expressed genes (DEGs) were identified using the DESeq2 package (v1.30.0) (Love et al., 2014). The P-value of the DEGs between samples was adjusted using the Benjamini & Hochberg method (Benjamini and Hochberg, 2000). Genes with an adjusted P-value ≤0.05 and | log2 (fold change)| ≥ 1 were recognized as DEGs. To determine the biological significance of the DEGs, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was implemented using KOBAS software (Wu et al., 2006).

The phenotypic analyses showed significant differences in flower color between the two parental lines. In Y641-87, the flower organs, particularly the petal tissue, showed stable yellow coloration at the flowering stage, whereas those of Y640-288 showed white coloration (Figures 1A,B). The flower color in Y640-288 was traditionally called white color. It was not like white papers as in B.napus and B.oleracea (Zhang B. et al., 2015; Han et al., 2019), but similar to that in B.juncea (Zhang et al., 2018a,b), the flowers of which still had pale yellow pigments. TEM analysis showed that the petals of Y641-87 had normal chromoplasts with numerous fully developed PGs, whereas the petals of Y640-288 showed abnormal chromoplasts with only a few irregular and small PGs (Figures 1C,D).

To investigate the genetic inheritance of white flowers in Chinese cabbage, we performed reciprocal crosses between Y641-87 and Y640-288. The resulting F₁ plants all displayed a yellow flower phenotype. The flowers of F₂ plants exhibited two types of colorations, corresponding to the coloration of either Y641-87 or Y640-288. Among 200 F₂ individuals, 142 exhibited yellow flowers, and 58 showed white flowers, corresponding to a segregation ratio of 3:1 by the Chi-square test (Table 1). In a larger F₂ population, the segregation ratio was also 3:1 (1775 yellow:596 white, χ² = 0.02). These results demonstrated that the inheritance of the white flower trait in Y640-288 follows a monogenic recessive pattern. We named this white flower gene BrWF3.

To investigate whether the lower pigmentation in white flowers was due to decreased carotenoid accumulation, we analyzed the carotenoid profiles of a white petal pool (W-bulk) and a yellow petal pool (Y-bulk). We detected 20 carotenoid peaks in the Y-bulk in the direct extracts (Figure 2 and Supplementary Table 6). Among these peaks, 9 peaks represented esterified carotenoids with retention times ranging from 5.5-7.5 min, as these peaks were not detected in the saponified sample (Figure 2). The esterified carotenoids are mostly derived from lutein and violaxanthin (Supplementary Table 6). In contrast, the composition and content of carotenoid esters in the W-bulk were much less than those in the Y-bulk (Figure 2). Analysis with saponification identified 10 carotenoids in both the W-bulk and Y-bulk (Figure 2 and Table 2). Two major xanthophylls, violaxanthin and lutein, accounted for approximately 83 and 91% of the total carotenoids in the Y-bulk and the W-bulk, respectively. The total average content of violaxanthin and all carotenoids in Y-bulk was about 2.76 and 1.70 times higher than that in the W-bulk, whereas lutein content did not significantly differ between Y-bulk and W-bulk (Table 2). Taken together, these results indicated that white petals accumulate less xanthophylls esters (likely violaxanthin esters) than yellow petals, resulting in lower carotenoid accumulation and color pigmentation.

To map the BrWF3 gene, we conducted BSA-seq using two pooled samples, which comprised 50 white-flowered (W-pool) and 50 yellow-flowered (Y-pool) F₂ plants, and two parental lines, Y640-288 and Y641-87. In total, 204, 209, 120, and 86 million raw data were generated for the W-pool, Y-pool, Y640-288 and Y641-87 (Supplementary Table 7), representing approximately 63-, 64-, 37- and 26-fold genome coverage, respectively, based on the estimated genome size of 485 Mb (Wang et al., 2011). The clean reads of each sample were mapped to the reference genome of the Chiifu cultivar. After filtering, a total of 358,141 SNPs and 54,500 InDels, which were distributed merely on ten chromosomes, were identified between the W-pool and the Y-pool. The Δ(SNP-index) of each position was calculated for sliding window analysis. According to the null hypothesis, a total of five adjacent regions on chromosome A02 (Supplementary Figure 1 and Supplementary Table 8) exhibiting significant linkage disequilibrium were identified as the candidate region for white flower trait at a 95% significant level. These results were not consistent with the assumption that the white flower trait is controlled by a single recessive nuclear genetic locus. However, most of the genomic regions on other chromosomes exhibited a Δ(SNP-index) = 0. In theory, the SNP-index of the W- and Y-pools should be the same in the genomic regions that are not related to the phenotypic diference (SNP-index = 0.5), and Δ(SNP-index) should equal to 0 (Islam et al., 2016; Wang et al., 2017). Therefore, the most likely chromosome where BrWF3 located was A02.

Based on BSA-seq analysis, 46 KASP markers previously available in our group (Yang et al., 2020) and 65 newly developed KASP markers, which were uniformly distributed across chromosome A02, were used to identify polymorphisms between the two parental lines (Y640-288 and Y641-87). The results showed that 29 KASP markers (Supplementary Table 4) exhibited polymorphism. These polymorphic markers were further genotyped in 135 F₂ plants for linkage analysis (Supplementary Table 9). The 135 F₂ plants is a sub-set from the ‘whole’ population containing 200 plants in Table 1. The results revealed no recombinant individuals between BrWF3 and markers TXBH57, TXBH58, TXBH62 and TXBH30 and 3 recombinant individuals between TXBH64 and BrWF3. The genetic distances between the BrWF3 locus and TXBH30 and TXBH64 were 0.7 and 2.0 cM, respectively (Figure 3A). The order of the markers in the genetic map is consistent with that in the physical map (Figure 3A).

To fine-map the BrWF3 locus, we screened 596 white-flowered F₂ individuals using flanking markers (TXBH46 and TXBH31) and identified 52 recombinants. All the 52 recombinants were further genotyped using TXBH57, TXBH58, TXBH62, TXBH 30, TXBH64 and TXBH65, based on which 10 recombinants (type 2 and type 7) were identified (Figure 3B). Then, 21 new markers were developed, and seven of which exhibited polymorphism in the two parents (Supplementary Table 4). These seven new polymorphic markers were further used to screen all the 10 recombinants using KASP assay. The results delimited the BrWF3 gene to a 105.6 kb interval between markers TXBH85 and TXBH84, each with one recombinant (type 3 and type 4) (Figure 3B). Three markers, SEQ14, TXBH76 and TXBH83, co-segregated with the BrWF3 gene (Figure 3B).

Based on the fine mapping results for BrWF3, DNA sequences in the 105.6 kb interval were analyzed in the Brassica database² and comparative gene annotation in Arabidopsis thaliana was performed. As a result, 13 annotated or predicted genes were identified in the mapping region (Table 3). Four of these genes, Bra032956, Bra032957, Bra032958, and Bra032959, are homologs of AT3G26840 (PES2) in Arabidopsis thaliana, which encodes a protein with phytyl ester synthesis and diacylglycerol acyltransferase activities and was previously reported to regulate carotenoid esterification (Zhang et al., 2018a,b; Kishimoto et al., 2020). Next, we examined the expression of these four candidate genes via RNA-seq and qRT-PCR analysis. RNA-seq analysis revealed that only Bra032957 was differentially expressed among the four genes with the expression level decreasing approximately three fold in white petals compared with yellow petals (Supplementary Table 10). qRT-PCR analysis showed that the expression of Bra032957 in petals was much higher than that of Bra032956 and Bra032959 and was significantly upregulated in yellow petals compared to white petals (Figure 4A). The results of qRT-PCR analysis and RNA-seq analysis were consistent. Subsequently, expression analysis in different tissues (root, stem, leaf, petal, sepal, stamen, and pistil) showed that Bra032957 was predominantly expressed in petals (Figure 4B). Taken together, the results indicated that the Bra032957 gene was the most likely candidate gene for BrWF3.

To characterize the sequence of the candidate genes in the white-flowered parental line Y640-288 and the yellow-flowered parental line Y641-87, the genomic sequence (gDNA) and coding sequence (CDS) of Bra032957 were amplified and sequenced using specific primers (Supplementary Figure 2 and Supplementary Table 5). The results showed that the Bra032957 gene of Y641-87 was 5162 bp in length and contained 14 exons and 13 introns. The CDS of the Bra032957 gene in Y641-87 was 2106 bp in length. Sequence alignment revealed a base deletion (G) at 477 bp of the CDS in the third exon in Y640-288 (Supplementary Figure 3 and Figure 5A). The SNP deletion caused a frameshift mutation in the Bra032957 gene and a premature stop codon in 168 a.a. residues (Figure 5A). Conserved domain analysis in NCBI showed that the Bra032957 gene contained an α/β hydrolase-fold (amino acids 123-380) and a lysophospholipid acyltransferases (LPLAT) domain (amino acids 430-666), which can transfer acyl groups to acceptors, such as glycerol 3-phosphate (Figure 5B). The SNP deletion mutation caused a loss of the two conserved domains and ultimately caused the loss of function of the BrWF3 protein.

Based on the identified SNP deletion, we designed a KASP marker TXBH83 to screen the other three F₂ populations (Y640-288 × SD369-F₂, Y640-288 × Chiifu-F₂, Y66-83 × R16-11-F₂), including a total of 282 individuals. The results showed that TXBH83 co-segregated with the flower color phenotype (Figures 6A-C and Supplementary Table 11). Furthermore, 10 white-flowered and 10 yellow-flowered DH lines were genotyped for TXBH83, which also showed a 100% consistency between the flower color phenotype and genotype (Figure 6D and Supplementary Table 11). Overall, these findings suggest that the Bra032957 gene is the most promising candidate gene for the white flower gene BrWF3 in Chinese cabbage.

To identify the gene regulatory networks involved in petal coloration, we performed comparative transcriptome analysis between W-bulk and Y-bulk. Approximately 300.8 million raw reads were generated for the six cDNA libraries, ranging from 43.1 to 57.4 Gb reads per library (Supplementary Table 12). All the raw reads were deposited in the NCBI Short Read Archive (SRA) database under accession number PRJNA682761. Among the clean reads, 85.2-87.46% were uniquely mapped to the reference genome (Wang et al., 2011). The Pearson correlation coefficients among the three replicates of each petal pool ranged from 0.98 to 0.99, indicating a high consistency among the three replicates (Supplementary Figure 4). In total, we identified 6,009 differentially expressed genes (DEGs) between the W-bulk and Y-bulk, among which 2,913 genes were up-regulated and 3,816 were down-regulated in the Y-bulk compared with the W-bulk. Pathway enrichment analysis based on the KEGG database revealed that carotenoid biosynthesis was the most significantly enriched pathway (Figure 7A). In the Y-bulk, genes involved in carotenoid biosynthesis (Figure 7B), such as PSY (Bra006391 and Bra008569), PDS (Bra010751), BCH1 (Bra013912) and ZEP (Bra037130), were significantly up-regulated. Genes involved in carotenoid degradation (Figure 7B), such as NCED3 (Bra027336 and Bra001552) and NCED4 (Bra013378), were also upregulated. Pathways of linoleic acid metabolism, alpha-linolenic acid metabolism, glycerophospholipid metabolism and arachidonic acid metabolism were also enriched, and most of the genes in these pathways were downregulated (Figure 7A and Supplementary Table 12). For example, genes encoding glycerophosphodiester phosphodiesterase (Bra040704, Bra027481, Bra035967, Bra020395, Bra002676, and Bra006785), phospholipase A (Bra015531 and Bra010327), phosphatidate phosphatase (Bra029774) were downregulated (Supplementary Figure 5 and Supplementary Table 13). However, genes participating in fatty acid elongation, such as 3-ketoacyl-CoA synthase (Bra011936, Bra004513, Bra024749, and Bra004034) and very-long-chain enoyl-CoA reductase (Bra008657), were significantly upregulated (Supplementary Figure 5 and Supplementary Table 13). These results suggested that saturated but not unsaturated fatty acids might be the main acyl group donors for esterification of xanthophylls. FIBRILLIN (FBN) and ACTIVITY OF BC1 COMPLEX KINASE (ABC1K) are the most abundant proteins in PGs (van Wijk and Kessler, 2017). In this study, FBN1b (Bra013602) and ABC1K8 (Bra024339) were significantly upregulated in Y-bulk (Figure 7C). Furthermore, FBN1b is highly expressed. The average fragments per kilobase million (FPKM) value of FBN1b was as high as 8383 in Y-bulk, 3148 in the W-bulk. VITAMIN E DEFICIENT 1 (VTE1), which encodes a key enzyme in tocopherol biosynthesis, was upregulated in Y-bulk but with no significant difference (Figure 7C and Supplementary Table 13). Lipoxygenase (LOX) are proteins recruited to chloroplast PGs and participated in jasmonate biosynthesis (van Wijk and Kessler, 2017). Genes encoding LOXs were all downregulated in yellow petals (Figure 7C and Supplementary Table 13), indicating that LOX proteins were not indispensible for PG development and formation in chromoplasts.