The rapeseed mutant (df74) with Chl-deficient sepals was obtained from the conventional variety Ningyou 18 (NY18), whose seeds were treated with 1% EMS for 12 h. The mutant with a pale-yellow sepal phenotype was identified and was selfed to generate an inbred line. The df74 is smaller and shorter than NY18. The investigation of the F2 individuals derived from the cross of df74 and NY18 showed that the dwarfism and sepal-specific chlorophyll-deficiency are not linked to each other. So, we focused on the sepal-specific chlorophyll-deficiency phenomenon in this study. The cross between df74 and NY18 was carried out, and the F₂ population was derived from the self-pollination of F₁ plants. The BC₁ was derived by the backcrosses of F₁ to df74. The phenotype of the reciprocal hybrid F₁ and the segregation ratio of F₂ and BC₁ population were used to detect the genetic pattern of the Chl-deficient sepal mutant. Another conventional variety Zhongshuang11 (ZS11) with normal green sepals was used to cross with df74 to construct a F₂ population for fine mapping the candidate gene. All the plants were grown in fields located in Jiangsu Academy of Agricultural Sciences, Nanjing, Jiangsu Province, China, under normal cultivation.

For pigment extraction, 200 mg fresh weight of sepals were harvested from df74 and NY18, respectively, and then extracted with 5 ml 95% ethanol. Specific Chl contents were measured at wavelengths of 665 and 649 nm with UV-2450 UV/Vis Spectrophotometers (Shimadzu Corporation, Kyoto, Japan) according to the method of Lichtenthaler (1987). Each measurement involved three biological replicates.

For transmission electron microscopy (TEM) analysis, sections were manufactured as described by Wang et al. (2016). Sepals at the same stage on df74 and NY18 were fixed in 2.5% glutaraldehyde and further fixed in 1% osmic acid, dehydrated in gradient acetone, embedded in Epon812 and sectioned, and then double stained with 2% uranyl acetate and lead citrate. The sections were observed and photographed under a Hitachi H-7650 TEM (Hitachi, Tokyo, Japan).

Two contrasting DNA pools were constructed with equal amounts of DNA from 25 green-sepal plants (G-pool) and 25 yellow-sepal plants (Y-pool) from the F₂ populations derived from df74 × NY18. The two bulks and two parents were sequenced and the data were generated using the Illumina HiSeq™ PE150 (Illumina, Inc.; San Diego, CA, United States) platform. The sequence depth of the parents were ×30 genome coverage and the two bulks × 45. Sequencing and data analysis were performed by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). The clean reads of each sample were aligned to the B. napus “Darmor-bzh” reference genome (Chalhoub et al., 2014) using BWA software (Li and Durbin, 2009), and multiple read pairs were removed by the SAMtools command “rmdup” (Li et al., 2009). Sequence variants including single-nucleotide polymorphisms (SNPs) and small insertions and deletions (INDELs) were called using the Unified Genotyper function in GATK software (McKenna et al., 2010). SNPs with read depth <8 or quality score in scale <20 were filtered out. Only the SNPs that were polymorphic between the parents and homozygous in either parent were selected for the further analysis. Using df74 as the reference parent, SNP-index of the two bulks were calculated and the corresponding SNP-index graphs were plotted to find the candidate regions responsible for phenotypic variation. The Δ(SNP-index) can be estimated by the formula: Δ(SNP-index) = SNP-index (G-pool) − SNP-index (Y-pool). The SNPs with significantly different Δ(SNP-index; p < 0.01) were considered closely linked to the causal gene (Yaobin et al., 2018).

A F₂ population with 3,911 individuals was derived from the cross between df74 and ZS11. The yellow-sepal individuals were selected for genotyping. According to the results of the sequencing of df74 and NY18, the SNPs between the two parents flanking the candidate interval were screened. The penta-primer amplification refractory mutation system (PARMS) was used to screen recombinant plants for fine-mapping. The primers were designed as described by Lu et al. (2020). The master mix for PARMS markers was purchased from Gentides Biotech Co. Ltd. (Wuhan, China). The parents, F₁, 10 yellow-sepal individuals, and 10 green-sepal individuals were used to detect the reliability of PARMS markers. Then, two polymorphic markers (M1 and M2) were used to screen the yellow-sepal individuals of the F₂ population. The recombinant individuals were selected to conduct the next round of genotyping. New PARMS SNP markers (M3–M10) were designed step by step based on the physical map of ZS11 (Song et al., 2020) to narrow the mapping interval gradually.

The sepals of NY18 and df74 were harvested and immediately frozen in liquid nitrogen and then stored at −80°C for RNA extraction and transcriptome analysis. Three biological replicates were performed. Total RNA was extracted using a MiniBEST Plant RNA Extraction Kit (TaKaRa, Dalian, China). RNA-sequencing (RNA-seq) was carried out by Novogene Bioinformatics Technology Co. Ltd. The RNA quality was determined by a Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States), RNA Nano 6000 Assay Kit, and Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). Then, six sequencing libraries were constructed and sequenced on Illumina Hiseq 2000 platform, and 100 bp paired-end reads were generated. These methods were described by Wang et al. (2020).

Raw RNA-seq reads were filtered by removing reads containing the adapter, and reads containing ploy-Ns and low-quality reads. Then, the clean reads were aligned to the B. napus “Darmor-bzh” reference genome using the hisat2 program (Wen, 2017). DESeq2 was used for filtering the DEGs with |logFC(log₂fold change)| ≥ 1 and fdr < 0.05. The Gene Ontology (GO) annotation of DEGs was performed by the GOseq R package (Young et al., 2010), and GO terms with corrected p < 0.05 were considered as significantly enriched terms. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation was carried out with a BLAST search against the KEGG database.¹ Pathway enrichment analysis was performed by KOBAS 2.0.²

In the seedling stage, the leaves of df74 were as green as NY18, but smaller than NY18 (Figure 1A). The df74 entered the bolting stage and flowering stage 3–5 days later than NY18, and had smaller floral organs (Figures 1B–F). The sepals of df74 were pale yellow, but the other green tissues including leaves, stems, and carpels, had the same green color as NY18. The Chl content of sepals in df74 was significantly different from NY18. Sepals of df74 had 82, 71, and 80% reduction of Chl a, Chl b, and total Chl, respectively, compared with NY18 (Figure 2).

To investigate how the mutation affects chloroplast development, we observed the ultrastructure of the chloroplasts in the sepals of df74 and NY18 with TEM. In the cells of NY18, the chloroplasts were neatly arranged beside the cell walls and contained well-developed thylakoid membrane systems (Figures 3A,C). However, the cells of df74 showed plasmolysis, the number of chloroplasts was significantly less than that of NY18, and the chloroplasts were occupied by disordered thylakoids and 1–2 giant starch grains (Figures 3B,D).

The phenotypes of the generations obtained by the crosses of df74 and NY18 were thoroughly investigated. The reciprocal F₁ plants exhibited green sepals as wild type, which indicated that the Chl-deficient sepal trait was controlled by nuclear genes. An F₂ population with 322 individuals contained 239 green-sepal and 83 yellow-sepal plants. A chi-squared test indicated that the segregation pattern agreed with the Mendelian segregation ratio of 3:1 (χ² = 0.067, p > 0.05). In addition, the ratio of green-sepal plants to yellow-sepal plants in the BC₁ progenies was approximately 1:1 (χ² = 0.019, p > 0.05). These results indicated that the phenotype of df74 was controlled by a single recessive nuclear gene.

Solexa sequencing of the two pools and two parental lines generated about 172 Gb clean data, with Q20 ≥ 97.97% and Q30 ≥ 93.08%. The data were aligned to the B. napus reference genome “Darmor-bzh.” The average read depth was 39.61-fold in NY18, 26.77-fold in df74, 35.94-fold in the G-pool, and 36.89-fold in the Y-pool. A total of 26,900 differential and homozygous SNP loci were detected between two parental lines. At a 99% significance level, a single peak in the 26.22–37.54 Mb region on chromosome C08 had an average Δ(SNP-index) of 0.55. The largest Δ(SNP-index) among this genomic region was 0.81 (Figure 4). These data confirmed the presence of a major gene located in this genomic region controlling Chl-deficient sepal phenotype, which was designated as BnC08.cds.

To fine-map the BnC08.cds gene, the candidate region derived from BSA-seq was aligned to the reference genome sequence of B. napus “ZS11.” According to the result, two flanking PARMS SNP markers, M1 (physical position of the B. napus “ZS11” genome: 44420751) and M10 (50875750) were first designed to investigate the recombinants among the 900 yellow-sepal individuals in the (ZS11 × df74) F₂ population. According to the result, 369 individuals containing recombinants between the two markers were selected. Furthermore, eight polymorphic PARMS SNP markers (M2–M9) were designed step by step according to the genome sequencing data (Table 1), and the recombinant individuals were analyzed. The BnC08.cds was finally mapped in an interval between M5 and M6 (Figure 5), both resulting in one recombinant. The region had a 228.72 kb physical distance and 18 candidate genes (Table 2) corresponding to the B. napus “ZS11” genome (Song et al., 2020).

To further understand the molecular mechanisms of the difference between NY18 and df74, high-throughput RNA-seq was performed using Illumina technology. After the raw data were trimmed, 19.90–23.54 million clean reads for each sample were aligned to the B. napus “Darmor-bzh” genome. A total of 4,108 DEGs between NY18 and df74 were identified, including 1,550 genes upregulated and 2,558 downregulated in the sepals of df74. To investigate the function of the DEGs, GO enrichment analysis was performed. DEGs were mainly categorized and highly enriched as follows, “cellular process,” “metabolic process,” and “single-organism process” in the “Biological Process” category, and “binding” and “catalytic activity” in the “Molecular Function” category. The most abundant terms in the “Cellular Component” category were “cell,” “cell part,” “membrane,” “organelle,” and “membrane part” (Supplementary Additional File 1). Functional annotation of the DEGs was performed against the KEGG database to further investigate the functions of DEGs. The 1,522 up- and 2,527 downregulated DEGs were classified into 116 and 117 pathways, respectively. The result showed the most significantly enriched KEGG pathways were pentose and glucuronate interconversions, starch and sucrose metabolism, and metabolic pathways (Supplementary Additional File 2).

As color change is related to Chl content, we focused on the porphyrin and Chl metabolism pathways. Three DEGs that encode proteins participating in Chl breakdown were upregulated in df74, and one DEG participating in Chl biosynthesis, was downregulated. The three upregulated genes included BnaA08g30360D, BnaC01g20010D, and BnaC01g01800D. The first two genes were both homologous genes of Methylesterase family member 16 (MES16) in Arabidopsis, which catalyzes the demethylation of fluorescent Chl catabolite (Li and Pu, 2016). The third, BnaC01g01800D, was homologous to accelerated cell death gene ACD2 in Arabidopsis, which encodes a Chl breakdown enzyme to protect cells from programmed cell death caused by porphyrin-related molecules (Pattanayak et al., 2012). The downregulated gene was BnaC09g47630D, which is the homologous gene of HEMC that encodes the porphobilinogen deaminase function in Chl biosynthesis (Lim et al., 1994). These pathway annotations of the DEGs revealed the functions of BnC08.cds.

According to the fine-mapping results, the BnC08.cds was located in an interval with a 228.72 kb physical distance and 18 candidate genes (Table 2). The RNA-seq analysis indicated that none of the 18 genes had significantly different expression levels between the sepals of df74 and NY18. Six genes (BnaC08G0441400ZS, BnaC08G0441600ZS, BnaC08G0441700ZS, BnaC08G0441800ZS, BnaC08G0442000ZS, and BnaC08G0442200ZS) were not expressed in either df74 or NY18. These genes were excluded from the candidate genes.

In this interval, there was an enrichment region of basic helix–loop–helix (bHLH) DNA-binding superfamily genes, including BnaC08G0441900ZS, BnaC08G0442000ZS, and BnaC08G0442200ZS. The bHLH superfamily is one of the largest transcription factor gene families in Arabidopsis (Heim, 2003). Members in the phytochrome-interacting factor family that belongs to the bHLH superfamily VII can negatively regulate Chl biosynthesis (Huq et al., 2004; Shin et al., 2009). In view of this previous research, we cloned and aligned the sequence of the three genes with bHLH domain in df74 and NY18. There was no difference in the three genes between df74 and NY18; therefore, these genes were excluded as candidate genes for BnC08.cds.

BnaC08G0441400ZS is homologous to CYP708A3 in Arabidopsis, which belongs to the cytochrome P450 gene family and functions in heme binding. We sequenced the BnaC08G0441400ZS in df74 and NY18. There is only one base mutation in the intron region, so it was also excluded as a candidate gene for BnC08.cds.

According to the result of resequencing, there was only one missense mutation in the 18 genes. In BnaC08G0442100ZS, there was a single nucleotide substitution (G–A) in the seventh exon, which changed serine into asparagine at the 232nd amino acid site. The mutation was confirmed by sequencing the gene in df74 and NY18. The orthologous gene of BnaC08G0442100ZS in A. thaliana is AT5G43745, which encodes an ion channel POLLUX-like protein and is located in the chloroplast and chloroplast envelope. Therefore, we speculated that BnaC08G0442100ZS may be the most likely candidate gene of BnC08.cds.