Identification of two recessive etiolation genes (py1, py2) in pakchoi (Brassica rapa L. ssp. chinensis)

Background Leaf color is a major agronomic trait, which has a strong influence on crop yields. Isolating leaf color mutants can represent valuable materials for research in chlorophyll (Chl) biosynthesis and metabolism regulation. Results In this study, we identified a stably inherited yellow leaf mutant derived from ‘Huaguan’ pakchoi variety via isolated microspore culture and designated as pylm. This mutant displayed yellow leaves after germination. Its etiolated phenotype was nonlethal and stable during the whole growth period. Its growth was weak and its hypocotyls were markedly elongated. Genetic analysis revealed that two recessive nuclear genes, named py1 and py2, are responsible for the etiolation phenotype. Bulked segregant RNA sequencing (BSR-Seq) showed that py1 and py2 were mapped on chromosomes A09 and A07, respectively. The genes were single Mendelian factors in F3:4 populations based on a 3:1 phenotypic segregation ratio. The py1 was localized to a 258.3-kb interval on a 34-gene genome. The differentially expressed gene BraA09004189 was detected in the py1 mapping region and regulated heme catabolism. One single-nucleotide polymorphism (SNP) of BraA09004189 occurred in pylm. A candidate gene-specific SNP marker in 1520 F3:4 yellow-colored individuals co-segregated with py1. For py2, 1860 recessive homozygous F3:4 individuals were investigated and localized py2 to a 4.4-kb interval. Of the five genes in this region, BraA07001774 was predicted as a candidate for py2. It encoded an embryo defective 1187 and a phosphotransferase related to chlorophyll deficiency and hypocotyl elongation. One SNP of BraA07001774 occurred in pylm. It caused a single amino acid mutation from Asp to Asn. According to quantitative real-time polymerase chain reaction (qRT-PCR), BraA07001774 was downregulated in pylm. Conclusions Our study identified a Chl deficiency mutant pylm in pakchoi. Two recessive nuclear genes named py1 and py2 had a significant effect on etiolation. Candidate genes regulating etiolation were identified as BraA09004189 and BraA07001774, respectively. These findings will elucidate chlorophyll metabolism and the molecular mechanisms of the gene interactions controlling pakchoi etiolation.


Background
The photosynthetic pigment chlorophyll (Chl) is ubiquitous in cyanobacteria and the chloroplasts of higher plants. Chl converts the energy of sunlight into bioavailable chemical energy which drives carbohydrate biosynthesis [1]. Chl is an essential component of leaf color which influences dry matter accumulation and crop yield. In general, leaves appear green because Chl predominates and has a vital role in them. When the Chl content changes in plants, various leaf color mutant phenotypes result including chlorina, virescent, albino, yellow-green, and stay-green [2]. Leaf color mutants develop from the inhibition of genes regulating Chl biosynthesis and chloroplast development. Downregulation of these genes directly or indirectly influences Chl synthesis and degradation and produces the leaf color mutations [3][4][5]. Thus, leaf color mutants may be ideally suited for the elucidation of the mechanisms of photosynthesis, Chl biosynthesis, chloroplast development, and the expression and regulation of the genes associated with these processes [6][7][8][9]. Leaf color mutants have been characterized in Arabidopsis [10], rice [11,12], wheat [13], Brassica napus [14], Brassica oleracea [15], barley [16], kale [17], tobacco [18], soybean [19], cotton [20], and cucumber [5,21]. Much research has been invested in the analysis of the genetics, physiology, and molecular mechanisms of Chl biosynthesis and chloroplast development via leaf color mutants.
Several studies in genetic analysis have categorized leaf color mutation inheritance as nuclear and cytoplasmic. Most leaf color mutations are recessively inherited and conferred by a single nuclear gene [21][22][23][24]. Leaf color mutations involving two recessive genes are rare. Moreover, their inheritance is complex and it is difficult to apply genetic analysis and gene mapping on them. Wu et al. identified the light color mutant ws1 in Nicotiana tabacum and determined that this phenotype was controlled by the recessive nuclear genes ws1a and ws1b localized by different BC 1 F 2 groups to linkage groups 5 and 24, respectively [18]. BnChd1-1 and BnChd1-2 are responsible for the light green leaf mutant phenotype in Brassica napus. Fine mapping of BnChd1-1 was achieved using the BC 3 F 1 population. Candidate gene prediction suggested that BnChd1-1 encodes a subunit of the nicotinamide adenine dinucleotide phosphate (NADPH) complex in the thylakoid lumen [25]. Chl-deficient mutant phenotypes in durum wheat [26] and Brassica juncea [27] are also controlled by two recessive genes. Cytoplasmic mutants are uncommon compared to nuclear mutants. However, they have been reported for tobacco [28], barley [29], and Brassica campestris [30].
In plants, Chl biosynthesis comprises 15 enzymatic steps regulated by at least 27 genes [31]. Inactivation mutations of the Chl biosynthetic genes usually results in Chldeficient mutants [11,22,32,33]. Mutations in the genes governing Chl degradation metabolism generally produce stay-green mutants which retain their green leaf phenotype even during senescence [24,34,35]. Chl and heme biosynthesis are two types of tetrapyrrole formation and share a common metabolic pathway from 5aminolevulinic acid (ALA) to protoporphyrin IX (Proto IX) [36]. Heme is essential for both respiration and photosynthesis. In contrast, excessive heme accumulation inhibits glutamyl-tRNA reductase activity and ALA synthesis, reduces the rate of tetrapyrrole biosynthesis, and affects Chl biosynthesis [37]. Leaf color mutants arising from abnormal heme metabolism were identified for Arabidopsis [38], rice [39][40][41], pea [42], and maize [43].
In a previous study, we developed a pakchoi (Brassica rapa L. ssp. chinensis) yellow leaf mutant (pylm) derived from the 'Huaguan' pakchoi variety by isolated microspore culture. This strain is a double haploid (DH) with a stable yellow leaf phenotype [44]. In the present study, we conducted genetic analysis on pylm using bulked segregant RNA sequencing (BSR-Seq) with linkage analysis to map the corresponding genes. Then, the candidate genes associated with the mutant phenotype were predicted. The information derived from this work may help facilitate the cloning of etiolation genes and elucidate the molecular mechanisms of gene interactions.

Phenotypic characterization of mutant pylm
The wild type 'CK-51' and pylm were both obtained from isolated microspore culture of the 'Huaguan' pakchoi variety. However, the latter displayed yellow leaves at germination and this phenotype was stable throughout its lifetime (Fig. 1). The mutant had a slender phenotype and weak growth. However, its yellow leaf color was nonlethal. Moreover, relative to 'CK-51', pylm displayed an elongated hypocotyl at the seedling stage ( Fig. 1c) and early flowering at the bolting stage (Fig. 1b).
Genetic analysis of mutant pylm F 1 and F 2 populations were constructed from crosses between pylm and the Chinese cabbage DH line 'FT' (Additional file 1: Figure S1). The F 1 individuals from the reciprocal crosses had the same green leaf phenotype as 'FT'. Therefore, inheritance of the etiolation phenotype in pylm is nuclear rather than cytoplasmic. Segregation statistics data for the green and yellow leaf phenotypes of the F 2 population accorded with the expected Mendelian ratio of 15:1 (χ 2 < χ 2 0.05 = 3.84). Thus, the Chl deficiency trait is controlled by two recessive nuclear genes. The BC 1 progeny was obtained from F 1 separately backcrossed with pylm and 'FT'. Segregation statistics data for the green and yellow leaf phenotypes of the BC 1 F 1 population from the cross of F 1 with pylm fitted the expected Mendelian ratio of 3:1 (χ 2 < χ 2 0.05 = 3.84). This finding confirms that the mutant trait is conferred by two recessive nuclear genes. They were designated as py1 and py2. Neither gene alone can induce the yellow leaf phenotype. Phenotypic data for all generations are listed in Table 1.

Segregation of py1 and py2
According to the genetic analysis of pylm, the reduced-Chl phenotype is controlled by the recessive nuclear genes py1 and py2. In consequence, the F 2 and BC 1 F 1 populations could not be used to map these genes separately. To separate py1 from py2, F 2 individuals with green leaves may be randomly selected and selfpollinated to produce F 2:3 progeny. Green-colored F 2:3 plants with a statistical segregation ratio of 3:1 (greencolored:yellow-colored) may be self-pollinated to generate F 3:4 progeny. In theory,~2/3 of the F 3:4 families should display the expected Mendelian segregation ratio of 3:1. Some of these could map py1 while the others could map py2 (Fig. 2).

BSR-Seq analysis
A total of 47,526,126 and 49,119,466 raw reads (150-bp) were generated from the G-pool and Y-pool, respectively. After quality evaluation and data filtering, 97% of the read pairs (46,456,174 for the G-pool and 47,581,728 for the Ypool) remained. Clean reads were mapped against the Brassica reference genome with Hisat v. 2.0.14. Of these, > 66% were uniquely mapped in both pools.  Relative to the reference genome, 154,863 and 157,022 SNPs were detected in the G-pool and Y-pool, respectively. Differential SNP loci were screened for ED 5 calculation and 412 target differential SNP loci were obtained between the pools according to the top 1% ED 5 threshold. Two distinct peaks were observed on chromosomes A07 and A09 (Fig. 3). This finding was consistent with the hypothesis that the mutant trait is controlled by two recessive nuclear genes. Thus, it was predicted that the etiolation genes were located on chromosomes A07 and A09 within five chromosome regions ( Table 2).

Identification of differentially expressed genes
RPKM was used to measure gene expression level. By setting RPKM ≥0.1, 55,250 genes were detected. These were divided into six RPKM distribution intervals (Additional file 3: Table S3). There were 181 DEGs between the G-pool and Y-pool according to the constraint (|log 2 fold change| ≥ 1 and FDR ≤ 0.05). Ninety genes were upregulated and the others were downregulated when the G-pool was compared with the Y-pool (Additional file 2: Figure S2). The DEGs are shown in Additional file 3: Table S4.

Fine mapping of py1
Ninety-six SSR markers were developed around the three predicted chromosome regions on chromosome A09. They were used to detect polymorphisms between pylm and 'FT'. After screening, thirty-seven SSR markers displayed polymorphisms between parents. They were used to test twelve green-colored and yellow-colored individuals each from the No. 1 F 3:4 family. SSRzk5 and SSRzk12 were located near the 23,811,435-27,563,122 region on chromosome A09 and showed linkage to py1 on the opposite side.
A total of 1520 yellow-colored individuals of the No. 1 F 3:4 family were selected as the py1 mapping population. A linkage analysis disclosed that py1 was located between SSRzk5 and SSRzk12 at estimated genetic distances of 3.2 cM and 1.8 cM, respectively (Fig. 4a). To identify the molecular markers tightly linked to py1 and narrow the py1 mapping interval, new SSR and Indel markers were developed between SSRzk5 and SSRzk12. The polymorphic markers SSRzk17, SSRzk28, SSRzk29, SSRzk36, Indelzk72, and Indelzk125 were linked to py1 (Additional file 3: Table S5). SSRzk17, Indelzk72, and Indelzk125 were located on one side of py1 as SSRzk5 while SSRzk28, SSRzk29, and SSRzk36 were located on the other side of py1 as SSRzk12. The py1 was mapped between Indelzk125 and SSRzk36 at 0.13 cM and 0.2 cM, respectively (Fig. 4b). Therefore, py1 was mapped in a 258.3-kb region between the most tightly linked markers (Fig. 4c). The target DNA sequences of the 258.3-kb region between Indelzk125 and SSRzk36 were obtained from the Brassica database. A genomic sequence analysis revealed that the candidate region contained 34 genes (Fig. 4c, Additional file 3: Table S6). Differential gene expression analysis disclosed only BraA09004189 in the py1 mapping region. BraA09004189 is a heme oxygenase (HO1) which participates in heme catabolism. Mutants with yellow leaf phenotype induced by defective HOs were reported in earlier studies [40,41]. BraA09004189 was predicted to be the most probable candidate py1 gene.
To confirm this hypothesis, two pairs of primers were designed to sequence BraA09004189 in pylm and 'CK-51' (Additional file 3: Table S7). The BraA sequence did not differ between parents whereas the BraB sequence in pylm presented with one SNP (Fig. 5). Based on the position of BraA09004189, an SNP marker was designed to screen 1520 yellow-colored individuals from the No. 1 F 3:4 family. The bands of whole mapping individuals cosegregated with py1.
A qRT-PCR was performed to determine BraA09004189 expression in pylm and 'CK-51'. In accordance with the differential gene expression analysis, the results indicated that BraA09004189 expression level was much higher in 'CK-51' than that in pylm (Fig. 6). This finding further supports the likelihood that BraA09004189 is the candidate for py1.

Fine mapping of py2
Considering the constructed populations size, we screened the Nos. 2-5 F 3:4 families using the same research strategy applied for SSRzk5 and SSRzk12 linked to py1. The etiolation gene py1 was identified in the Nos. 2, 4, and 5 F 3:4 families. In theory, then, the No. 3 F 3:4 family may be used to establish the py2 locus.
Forty-eight SSR markers were developed around the two predicted regions on chromosome A07 to detect polymorphisms between pylm and 'FT'. After screening, eleven SSR markers displayed polymorphisms between the parents. They were used to test twelve green-colored and twelve yellow-colored individuals of the No. 3 F 3:4 family. SSR84 and SSR103 were located around the   Table S8). There were 1860 yellow-colored individuals from the No. 3 F 3:4 family selected as the py2 mapping population. The py2 was located between SSR11 and SSR15 at estimated genetic distances of 0.24 cM and 0.02 cM, respectively (Fig. 7b). To narrow the py2 mapping interval and identify the molecular markers tightly linked to py2, new SNP markers were developed between SSR11 and SSR15. Only the polymorphic marker SNP11 was linked to py2. Based on the recombinant individuals, the py2 interval was narrowed to 14,851,951-14,896,902 and contained five genes (Fig. 7c).

Candidate py2 analysis
Annotation data for the five candidate genes in the py2 target region were obtained from the Brassica database (Additional file 3: Table S9). Primers were designed to cover the cDNA for each gene and predict the candidate genes (Additional file 3: Table S10). There were no differences between pylm and 'CK-51' in terms of BraA07001775, BraA07001776, or BraA07001777. After PCR amplification, the BraA07001773 sequence was disordered and the sequence comparisons were inconsistent over serial repetitions. There was SNP variation between parents for the first exon in BraA07001774 (Fig. 8). It caused a single amino acid mutation from Asp (GAT) in the wild type to Asn (AAT) in pylm (Fig. 9). Therefore, BraA07001774 was taken as the most probable candidate gene for py2.
BraA07001774 is an embryo defective 1187 (emb 1187) and a phosphotransferase. The albino mutants (pds1, pds2) phenotypes in Arabidopsis thaliana may be caused by emb 71 [45]. For Arabidopsis seeds with silique defects, hypocotyl elongation was characterized during the development of F 2 generation mutant seedlings [46]. We proposed that the mutant phenotype is determined by mutations in BraA07001774. To validate our prediction, BraA07001774 expression in pylm and 'CK-51' was analyzed by qRT-PCR. BraA07001774 was dramatically  Table S6 downregulated in pylm (Fig. 10). Thus, it probably is the candidate gene for py2.

Discussion
Mutations in leaf color are widespread in nature. The main type of leaf color mutation is Chl deficiency. Dwarfism, retarded growth, attenuated photosynthetic capacity, low yield, and death are associated with this defect [25,47,48]. Here, we identified the pakchoi yellow leaf mutant pylm from isolated microspore culture. Unlike previously reported Chl deficient mutants, pylm presented with substantially elongated hypocotyls at the seedling stage and early flowering at the bolting stage. The etiolation phenotype in pylm was nonlethal and stable throughout the growth period. The Chl deficiency in pylm was controlled by two recessive genes. These characteristics suggested that pylm was of high value for  Map-based cloning is an effective gene isolation strategy. It has been extensively used for plant gene function analysis [49][50][51]. However, it is contingent upon fine mapping of the target gene. For most leaf color mutants, the traits are recessively inherited and controlled by a single nuclear gene. The F 2 populations are instinctively applied to map the target gene [5,41,48]. With regard to the character conferred by two recessive nuclear genes, F 2 populations may also be used in preliminary mapping. An efficient way to isolate allele pairs from each other and separately map them is to construct advanced backcrosses and other populations. The recessive white stem (ws) loci in Nicotiana tabacum and the Chl deficiency (Bnchd1) loci in Brassica napus were successfully mapped using constructed BC 1 F 2 and BC 3 F 1 populations, respectively [18,25]. In the present study, genetic analysis revealed that the recessive nuclear genes designated as py1 and py2 were responsible for the etiolation trait. We successfully segregated py1 from py2 and constructed an inheritance model for the Chl deficiency trait in pakchoi. Twenty F 3:4 families with a phenotypic segregation ratio of 3:1 were constructed. Various F 3:4 families were successfully used to map the py1 and py2 loci separately. Compared to using advanced backcross populations to map pairs of recessive nuclear genes, creating and using F 2:3 or F 3:4 families avoid the selection errors and interference in genetic analysis caused by the incomplete emasculation of Brassica rapa.
BSR-Seq efficiently combines the respective superiorities of bulk segregation analysis (BSA) and RNA sequencing (RNA-Seq) for rapid gene mapping [52][53][54][55]. BSR-Seq is targeted at the mRNA level. It selects phenotypically opposite individuals from segregated populations and constructs two RNA mixing pools to find SNPs at the transcript level. The transcriptome data localize the target genes and detect potentially associated DEGs [56]. BSR-Seq has been extensively applied to map the causal genes linked to a single target trait [24,57,58]. Two independently inherited traits may also be localized by BSR-Seq. Tan et al. applied BSR-Seq to locate the genes controlling male sterility on chromosome A05 and white petal on chromosome A02 [59]. Here, the mutant pylm and Chinese cabbage DH 'FT' lines were chosen as parents to construct the F 2 separation population. The phenotypes of the wild type and mutant individuals significantly differed. Thus, extreme mixed pools could be accurately and conveniently created for BSR-Seq. Release of Brassica rapa genomic data enhanced the reliability of these populations in BSR-Seq Fig. 7 Genetic and physical maps of py2 and candidate gene analysis. a: Linkage map of chromosome A07 was constructed with 1860 individuals bearing the pylm phenotype in the No. 3 F 3:4 family. The py2 was preliminarily mapped between SSR133 and SSR103; b: Fine mapping of py2. The py2 was restricted to the region between SSR11 and SSR15. The number of recombinants between the markers and py2 is shown under the genetic map. The distance above the linkage map is in centimorgan (cM) units; c: Candidate py2 region and annotated genes in the Brassica database. The py2 locus was narrowed to a 4.4-kb region containing the five predicted genes BraA07001773-BraA07001777. Arrows indicate the direction of gene expression. Detailed information on these genes is presented in Table S9 applications. Five candidate regions related to the yellow leaf phenotype in pylm were identified on chromosomes A07 and A09. Molecular markers were developed according to the locations of the candidate regions. The etiolation genes py1 and py2 were separately mapped with different F 3:4 families. The findings confirmed the feasibility of BSR-Seq for mapping two recessive nuclear genes. They also showed that BSR-Seq simplifies molecular marker development and screening in traditional mapping methods and greatly improves their efficiency.
New molecular markers were developed near the target regions based on preliminary py1 mapping by BSR-Seq. The py1 was mapped between the markers Indelzk125 and SSRzk36 on A09 chromosome over a 258.3-kb localization interval containing 34 predicted genes. No new polymorphic SSR or Indel markers were available to limit the localization interval. The gene expression patterns determined by BSR-Seq disclosed only one differentially expressed gene (BraA09004189) within the py1 mapping region. A gene annotation referenced from the Brassica database indicated that BraA09004189  HO1). This enzyme plays a vital role in phytochrome chromophore metabolism, the photoresponse mechanism, adventitious root formation, and oxidative damage mitigation [60][61][62][63]. HO1 stabilizes and maintains the heme content by transforming heme into BV-IXα [64]. As embranchments of tetrapyrrole biosynthesis, Chl and heme biosynthesis share a common metabolic pathway from ALA to Proto IX. Excessive heme accumulation caused by abnormal heme metabolism leads to feedback inhibition of Chl biosynthesis [36]. Therefore, a decrease in HO1 activity may influence Chl biosynthesis. The hy1 mutant of Arabidopsis and the yellow-green leaf 2 mutant of rice presented with the reduced-Chl phenotype because of free heme inhibition resulting from HO1 mutations [40,60]. HO1 defects strongly affected thylakoid development in rice [41]. Davis et al. found that the abnormally elongated hypocotyl phenotype of mutant Arabidopsis seedlings may be associated with HO1 defects [60]. Thus, py1 may encode HO and a mutation thereof may influence Chl biosynthesis and leaf color. In this study, qRT-PCR demonstrated that BraA09004189 was downregulated in pylm. This finding was consistent with those obtained by BSR-Seq. The SNP BraA09004189 was detected between pylm and 'CK-51'. A candidate gene-specific SNP marker in 1520 F 3:4 yellow-colored individuals co-segregated with py1. Thus, BraA09004189 corresponds to the yellow leaf locus py1 in pylm.
It was already known that certain Chl deficiency traits are controlled by two recessive nuclear genes. However, there was a lack of appropriate mapping populations or reliable molecule markers. Therefore, they were either approximately mapped without definite locations [18,26] or only one of the pair could be localized [25]. In previous studies, little progress was made in the simultaneous fine mapping or accurate prediction of the candidate genes. Here, we used the same mapping strategy as that for py1 to accomplish fine mapping for py2. The linkage analysis disclosed that py2 was mapped between SSR11 and SSR15 on chromosome A07. The mapping interval was narrowed to 4.4 kb by the SNP11 marker  linked to py2. Sequence analysis of the five genes in the py2 localization interval showed that only BraA07001774 expression significantly differed between pylm and 'CK-51'. For pylm, BraA07001774 had a SNP missense mutation on the first exon such that the wild type had an Asp residue whereas pylm had an Asn. The qRT-PCR revealed that BraA07001774 was downregulated in the mutant relative to the wild type. Gene annotation in the Brassica database indicated that BraA09004189 encodes emb 1187 and a phosphotransferase. In Arabidopsis, the emb genes are essential for seed development [65]. EMB genes encode various proteins. Thirty percent of them are active in the plastids [66]. Most emb mutations result in albinism or etiolated seeds and embryos which are secondary effects of mutations in chloroplast biogenesis and function [67]. The albino mutants (pds1, pds2) and hypocotyl elongation phenotypes in Arabidopsis may be related to mutations in EMB genes [45,46]. Thus, BraA07001774 is a candidate gene for py2.

Conclusions
We reported the identification of a Chl deficiency mutant pylm in pakchoi in this study. The etiolation trait was controlled by two recessive nuclear genes py1 and py2. We successfully segregated py1 from py2 by constructing F 3:4 families and achieved fine mapping and predictions for the two etiolation genes. Candidate genes regulating etiolation were identified as BraA09004189 and BraA07001774, respectively. These discoveries may help elucidate the molecular mechanisms underlying the trait controlled by two recessive nuclear genes. In future studies, functional validation will be conducted to clarify the functions of these candidate genes. In this manner, the molecular mechanism of gene interactions may be better understood.

Plant materials and mapping population development
The DH line pylm was obtained by isolated microspore culture of the 'Huaguan' pakchoi variety introduced from Japan musashino seed company. This strain was characterized by yellow leaves and elongated hypocotyls [44]. The parent used for the segregating population development with pylm in this study was 'FT', a DH line derived from Chinese cabbage variety 'Fukuda 50' screened by Shenyang greenstar Chinese cabbage research institute (Shenyang, China), which exhibits folded green leaves [68]. The pylm was reciprocally crossed with 'FT' to produce the F 1 , F 2 , and BC 1 generations. Twenty F 2 individuals with green leaves were selfpollinated to produce F 2:3 progenies. Eight green F 2:3 individuals per group were randomly selected from the corresponding populations and self-pollinated to produce the F 3:4 populations. Those with character segregation were used in linkage analysis and gene mapping. All plants were grown in the greenhouse at Shenyang Agricultural University, China.

Detection of variations by BSR-Seq
One hundred individuals with extreme leaf color phenotype at three-leaf stage were separately collected from the F 2 progeny and pooled for RNA extraction. Total RNA of each sample was extracted using TRIzol Reagent (Invitrogen, USA). The RNA concentration and integrity were analyzed with an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). The extreme mixed pools green-leaf (G-pool) and yellow-leaf (Ypool) were constructed by mixing equal amounts of each RNA sample. RNA-Seq library preparations were constructed according to the manufacturer's protocol (NEBNext® Ultra™ RNA Library Prep Kit for Illu-mina®). Sequencing was run on an Illumina HiSeq 2500 by GENEWIZ Suzhou Biological Technology Co., Ltd., China.
The quality of the raw RNA-Seq reads was assessed with FastQC (v. 0.10.1). Adapter sequences and lowquality reads containing N and < 70 were deleted from the raw reads with Cutadapt (v. 1.9.1). Low-quality bases at the 5′ or 3′ end were filtered out. Those with mean quality < 20 were trimmed by the 4-bp sliding window method. Clean data were aligned to the Brassica database (http://brassicadb.org/brad/) with Hisat v. 2.0.14 [69]. Candidate single nucleotide polymorphisms (SNPs) between the pools were obtained with the Mpileup module in SAM Tools (v. 0.1.18) SNP loci with depth coverage > 3× were screened for differential SNP analysis in the mutant and wild type pools. Euclidean distances (ED) for the differential SNP loci were calculated. The ED for each differential SNP locus was raised to the power five, namely, ED 5 , to eliminate background noise [70]. All ED 5 were sorted and the differential SNP loci with ED 5 in the top 1% were screened and mapped to specific chromosome regions based on the SNP locus distributions. Chromosome regions associated with the target traits were predicted according to the distributions of the ED 5 for the differential SNP loci on the chromosomes.

Differential gene expression analysis
To detect differentially expressed genes (DEGs) between the pools, a gene expression level analysis was performed with Htseq (v. 0.6.1). The reads per kilobases per million mapped reads (RPKM) were calculated [71]. DEGs were screened using a preset threshold (|log 2 fold change| ≥ 1 and false discovery rate (FDR) ≤ 0.05).