Fine-mapping and transcriptome analysis of a candidate gene controlling plant height in Brassica napus L.

Background Brassica napus provides approximately 13–16% of global vegetable oil for human consumption and biodiesel production. Plant height (PH) is a key trait that affects plant architecture, seed yield and harvest index. However, the genetic mechanism of PH in B. napus is poorly understood. Results A dwarf mutant df59 was isolated from a large-scale screening of an ethyl methanesulphonate-mutagenized rapeseed variety Ningyou 18. A genetic analysis showed that the dwarfism phenotype was controlled by one semi-dominant gene, which was mapped on C9 chromosome by quantitative trait loci sequencing analysis and designated as BnaDwf.C9. To fine-map BnaDwf.C9, two F2 populations were constructed from crosses between conventional rapeseed cultivars (Zhongshuang 11 and Holly) and df59. BnaDwf.C9 was fine-mapped to the region between single-nucleotide polymorphism (SNP) markers M14 and M4, corresponding to a 120.87-kb interval of the B. napus ‘Darmor-bzh’ genome. Within this interval, seven, eight and nine annotated or predicted genes were identified in “Darmor-bzh”, “Ningyou 7” and “Zhongshuang 11” reference genomes, respectively. In addition, a comparative transcriptome analysis was performed using stem tips from Ningyou 18 and df59 at the stem elongation stage. In total, 3995 differentially expressed genes (DEGs) were identified. Among them, 118 DEGs were clustered in plant hormone-related signal transduction pathways, including 81 DEGs were enriched in auxin signal transduction. Combining the results of fine-mapping and transcriptome analyses, BnaC09g20450D was considered a candidate gene for BnaDwf.C9, which contains a SNP that co-segregated in 4746 individuals. Finally, a PCR-based marker was developed based on the SNP in BnaC09g20450D. Conclusions The combination of quantitative trait loci sequencing, fine-mapping and genome-wide transcriptomic analysis revealed one candidate gene located within the confidence interval of 120.87-kb region. This study provides a new genetic resource for semi-dwarf breeding and new insights into understanding the genetic architecture of PH in B. napus.

of stearic acid [1]. Vegetable oils are triglyceride sources for biodiesel production. In Europe, biodiesel has been produced mainly from rapeseed oil [2]. In addition, rape straw is an abundant lignocellulosic material for the production of liquid biofuel, particularly ethanol. In China, the cultivated area of rapeseed is ~ 67 million hectares, with an annual seed yield of ~ 4.5 million tons every year [3], and the collectable production of rape straw was 38.17 million tons in 2013 [4]. Plant height (PH) is a key trait that affects the plant architecture, seed yield, dry weight and harvest index [5]. Moreover, PH is the major contributor to lodging tolerance, a serious abiotic stress during rapeseed production. Lodging makes B. napus unsuitable for mechanical harvesting and causes dramatic decreases in yield and seed quality [6]. Therefore, it is important to understand the genetic bases of PH to breed new cultivars with an ideal plant architecture and to maximize B. napus' economic benefits as an oil and bioenergy crop.
Plant growth and development consist of very precise and complicated procedures. The molecular mechanisms regulating PH in the model plant Arabidopsis thaliana and rice (Oryza sativa) have been well recorded, and most of the genes are involved in phytohormone-related pathways [7][8][9][10][11]. Gibberellic acid (GA) is an essential endogenous regulator of PH. For example, the rice 'green revolution' gene sd1 encodes GA20ox2, which is an oxidase enzyme involved in the biosynthesis of GA [12]. The wheat 'green revolution' gene Rht-B1/Rht-D1 encodes a DELLA protein that acts as a negative regulator in the GA-signaling pathway [13]. Auxins play pivotal functions in developmental processes, because they are involved in controlling virtually every aspect of plant biology [14]. Indole-3-acetic acid (IAA) is the key auxin in most plants, and it is mainly biosynthesized from tryptophan (Trp) through a Trp-dependent pathway [14,15]. Lossof-function mutations in IAA-related Trp-dependent biosynthetic genes (e.g., TAA1 or YUCCA ) can seriously affect PH [16,17]. In the auxin-signaling pathway, transport inhibitor resistant1/auxin-signaling F-box (TIR1/ AFB), auxin/indole acetic acid proteins (Aux/IAAs) and auxin response factors (ARFs) are three key factors regulating auxin-induced gene expression [18,19]. At low auxin levels, the transcriptional repressor Aux/IAA proteins and the corepressor TOPLESS interact with ARF proteins, resulting in auxin-induced gene repression. At high auxin levels, auxin interacts with TIR1, leading to the formation of TIR1/AFB-Aux/IAA-ARF complexes, then ubiquitin-ligase (E3) ubiquitinates AUX/IAAs, forming the ubiquitin 26S proteasome responsible for the degradation of AUX/IAAs [16,18,20]. Free ARF proteins can activate the transcription of auxin-response genes. Mutations that occur in the genes involved in auxin signaling may result in typical dwarf phenotypes, such as IAA7 in B. napus [21].
Brassica napus is a tetraploid species, which originated from a hybridization between B. rapa (AA, 2n = 20) and B. oleracea (CC, 2n = 18) around 7500 years ago [22]. In B. napus, many studies have focused on quantitative trait loci (QTLs) detection using linkage mapping or genomewide association studies (GWASs), and hundreds of QTLs have been identified across all 19 chromosomes [23][24][25][26][27][28][29][30][31]. However, most QTLs explain a small percentage of the total phenotypic variance (PV) [23][24][25][26][27][28][29][30][31]. Only a few QTLs affecting PH have been fine-mapped and cloned in B. napus. Wang et al. obtained two dwarf mutants "Bndwf1" and "Bndwf1/dcl1", and fine-mapped the associated QTLs on the A9 chromosome to a 152-kb interval and on the C5 chromosome to a 175-kb interval, respectively [32,33]. Liu et al. mapped the semi-dominant gene ds-1 (BnaA06.RGA ) on the A6 chromosome. It encodes a DELLA protein that functions as a repressor in GA signaling [34]. A single proline (P)-to-leucine (L) change was identified in the VHYNP motif of a DELLA protein in ds-1 that leads to a gain-of-function and caused a dwarf phenotype [34]. Subsequently, the DS-3 gene, which is syntenic to ds-1 with a similar gene function but a weaker effect on PH, was mapped on the C7 chromosome [35]. In the auxin-signaling pathway, only two genes that encode Aux/IAA proteins were identified and functionally validated in B. napus, BnaA3.IAA7 on the A3 chromosome and its homolog BnaC05.IAA7 on the C5 chromosome [21,36]. The molecular mechanisms regulating the PH of B. napus remain elusive, and elucidating the mechanism of a new dwarf gene has important scientific significance and applicable value.
The rapeseed ideotype is a semi-dwarf stature with a plant height of ~ 120-140 cm, narrow branch angles (< 30°), and a middle-long silique length [21,37]. A dwarf mutant of df59 with a PH of ~ 65 cm was obtained from ethyl methanesulphonate (EMS)-mutagenized Ningyou 18 (NY18). The df59 is an excellent germplasm resource for semi-dwarf breeding, with the average PH of F 1 hybrid lines from crosses between NY18 and df59 being 126.75 ± 4.3 cm, which is in accordance with the ideotype criteria of B. napus [21,37]. The aims of the present study were to: (1) fine-map the gene responsible for dwarf architecture in df59 using QTL sequencing (QTL-seq) and map-based cloning strategies; (2) elucidate the patterns of gene expression between NY18 and df59 using comparative transcriptomic analyses; and (3) develop a stable single-nucleotide polymorphism (SNP) marker tightly linked to the dwarf gene that could be used for marker-assisted selection. The present study provides a new gene source for the semi-dwarf breeding of new varieties, and the findings contribute to a better understanding of the molecular mechanisms underlying dwarfism.

Phenotypic variation and genetic analysis of plant height
At the seedling stage, df59 was already significantly shorter than NY18, with smaller leaves and shorter petioles (Fig. 1a, b). The internode length and PH were different between df59 and NY18 at the mature stage (Fig. 1c), while the petals and siliques of df59 were smaller than those of NY18 (Fig. 1d, e). At maturity, agronomic traits and seed yield-related traits were investigated in NY18, df59 and the F 1 of their cross (Table 1). Among the 15 traits, 12 were significantly higher in NY18 than in df59 (Table 1), with the exception of the first effective branch number, pod number of main inflorescence and silique number per plant. The PH of the F 1 was 126.75 ± 4.3 cm, which meets the ideotype requirement of ~ 120-140 cm and will be beneficial to high-density planting (Fig. 1f, Table 1). In addition, the seed yield per F 1 plant was 28.4 ± 2.20 g, which was 84.02% of the NY18 yield (33.8 ± 2.51 g). The phenotypic values of seed oil content and seed fatty acid concentrations for the NY18, df59 and the F 1 are provided in Additional file 1: Table S1. F 1 and  NY18 had similar C16:0, C18:0, C18:2, C18:3, C20:1 and  oil content, but had different C18:1 and C22:1 content. Theoretically, if the planting density of the F 1 increased to 1.5 times, the seed yield would increase ~ 26% per unit of area compared with the normal planting density of NY18.
Root-related traits, including total root length, root surface area, root volume and number of root tips, were measured for NY18, df59 and their F 1 at 10 days after germination (Additional file 2: Figure S1). The values of the four traits in df59 were significantly less than in NY18, and the F 1 had approximate mean values of the two parents, which were in accordance with the PH value ( Table 1).
The 165 individuals of the F 2 population from a cross between NY18 and df59 (named NY-DF) observed contained 43 tall plants, 88 medium plants and 34 dwarf plants. A Chi-squared test indicated that the segregation pattern agreed with the Mendelian segregation ratio of 1:2:1 (χ 2 = 1.719, P > 0.05). The mixed major-gene plus polygenes inheritance model in the software package SEA-G4F 2 [38], was used to perform the genetic analysis of PH for the NY-DF F 2 population along with the two parents and the F 1 . The results showed that PH is controlled by one major gene with additive-dominant effects. The major gene heritability was 80.04%, with an additive effect of 28.73 and a dominant effect of − 3.12 (Additional file 3: Table S2). The high heritability of the major gene indicated that PH in df59 was relatively stable and not greatly influenced by environment; therefore, it can be selected in the early generations of a breeding process.

QTL-seq of the NY-DF population
Two contrasting DNA pools were constructed from 16 extremely tall lines (T-pool) and 24 extremely dwarf lines (D-pool) in NY-DF F 2 population. Illumina highthroughput sequencing generated a total of 92.986 Gb clean data for the two pools and two parental lines, with average Q20 ≥ 95.66% and Q30 ≥ 93.14%. When the obtained reads were matched to the B. napus "Darmor-bzh" reference genome, the results showed that the depths of sequencing coverage for T-pool, D-pool, NY18 and df59 were 35.85-fold, 31.27-fold, 16.79-fold and 17.84-fold, respectively.
Based on the genotyping, 4520 polymorphic SNPs with haplotype differences between the two parents were identified. The region (17.00-27.30 Mb) on chromosome C9 had an average SNP-index of 0.726 in the D-pool, with the highest being 1 (Fig. 2a), while the SNP-index in T-pool was 0.128, with the lowest being 0 (Fig. 2b). Then, the genome sequence of the parent NY18 was used as a reference to calculate the Δ(SNP-index) of the 4520 SNPs by combining the SNP-indices of the T-and D-pools. At a 95% significance level, the genomic region from 17.00 to 27.30 Mb had an average Δ(SNP-index) of 0.60 ( Fig. 2c), suggesting that this region harbored a major QTL for PH, which was designated BnaDwf.C9.

Fine-mapping the BnaDwf.C9 locus
To fine-map the BnaDwf.C9 locus, two flanking pentaprimer amplification refractory mutation system (PARMS) SNP markers, M1 (physical position of the B. napus 'Darmor-bzh' genome: 17,001,732) and M11 (27,299,112), were first designed to investigate recombinants in the F 2 population (HO-DF) from a cross between Holly and df59 (Additional file 4: Table S3). The results showed that 81 recombinants were obtained among the 2356 F 2 individuals, including 56 tall recombinants and 25 dwarf recombinants (Additional file 4: Table S3). Based on the re-sequencing information of NY18 and df59 at the target region, nine polymorphic SNP markers (M2-M10) between M1 and M11 were designed and used for genotyping the 81 recombinant individuals. Recombinant genotypes showed that M2 was located on the left side of BnaDwf.C9, while six SNP markers (M5-M10) were mapped on the right side ( Fig. 3, Additional file 4: Table S3). In addition, M3 and M4 co-segregated with BnaDwf.C9 and were consistent with the PH phenotype. Eventually, the location of BnaDwf.C9 was narrowed down to the interval between M2 and M5 in a region of 770.72 kb in the 'Darmor-bzh' genome. No polymorphic SNP markers were able to narrow the interval in the HO-DF population.
For the further fine-mapping of BnaDwf.C9, a population of 2210 F 2 individuals was constructed from a cross between Zhongshuang 11 (ZS11) and df59 (ZS-DF population). Using the methods described above, 34 recombinant individuals were detected between M2 and M5 in the ZS-DF population (Table 2). Subsequently, three new polymorphic SNP markers in this region (M12-M14) were designed, together with co-segregating SNP markers (M3 and M4) in the HO-DF population and used to screen the 34 recombinants (Additional file 5: Table S4), resulting in 6, 5, 4, 0 and 2 recombinants, respectively ( Fig. 3, Table 2). These results suggested that the locus BnaDwf.C9 was fine-mapped to the region between M14 and M4, corresponding a 120.87 kb interval in the B. napus 'Darmor-bzh' genome ( Fig. 3).

Genome-wide transcriptomic analyses of NY18 and df59
High-throughput RNA sequencing (RNA-seq) generated 20.32-26.83 million raw reads for each sample. After trimming the low-quality sequences, 16.29-22.37 million clean reads were aligned to the B. napus 'Darmor-bzh' genome (Additional file 6: Table S5). A total of 3995 differentially expressed genes (DEGs) were identified between NY18 and df59 stem tip transcriptomes using the threshold false discovery rate < 0.005 and at least a 2.0-fold expression changed. In total, 1266 genes were significantly upregulated and 2729 genes were significantly down-regulated in df59 compared with NY18 (Additional file 7: Table S6).
To understand gene functions associated with the dwarf phenotype in df59, a Gene ontology (GO) enrichment analysis of the DEGs was performed. For both upand down-regulated genes, the most three significantly enriched GO terms in the 'biological process' , 'cellular component' and 'molecular function' groups were the same (Additional file 8: Table S7). For example, "metabolic process", "cellular process" and "single-organism process" were the three most enriched GO terms of the "biological process" group.
To discern DEG functions, we also conducted a Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The 1266 up-and 2729 down-regulated genes were enriched in 92 and 114 pathways, respectively (Additional file 9: Table S8). The top 20 enriched pathways for all of the up-and down-regulated DEGs are shown in Fig. 4. The five most significantly enriched pathways were carbohydrate metabolism, signal transduction, biosynthesis of other secondary metabolites, amino acid metabolism, and global and overview (Fig. 4). As the PHs in A. thaliana and rice were mainly regulated by plant hormones [7][8][9][10][11], we focused on the DEGs involved in signal transduction pathway (138 genes). Among them, 118 DEGs were clustered in plant hormone-related signal transduction pathways. In total, 81 DEGs were enriched in auxin signal transduction, accounting for 68.6% of the 118 DEGs (Additional file 10: Table S9). The remaining 37 DEGs were enriched in brassinosteroid (10 genes), abscisic acid (seven genes), GA (five genes), ethylene (five genes), jasmonic acid (four genes), salicylic acid (four genes) and cytokinin (two genes) signaling pathways (Additional file 10: Table S9). This result strongly indicated that the mutation leading to the dwarf phenotype of df59 may occur in an important gene in the auxin signal transduction. The KEGG pathway enrichment analysis provided an important clue for identifying candidate genes.

Identification of a candidate gene
Candidate genes underlying the 17.42-17.54 Mb of BnaDwf.C9 were analyzed based on the B. napus "Darmor-bzh" reference genome [22]. The segment harboring the BnaDwf.C9 locus contained seven annotated or  (Table 3). An RNA-seq analysis demonstrated that none of the expression levels among the seven genes significantly differed between NY18 and df59 at the stem elongation stage (Additional file 11: Table S10). Interestingly, BnaC09g20480D, BnaC09g20490D and BnaC09g20500D were not expressed at all in either NY18 or df59. Three other genes, BnaC09g20440D, BnaC09g20460D and BnaC09g20470D, showed no sequence differences between NY18 and df59 in the open reading frames (ORFs) based on the resequencing results. To confirm this result, the ORFs of the three genes in NY18 and df59 were amplified. DNA sequencing also revealed that the ORFs in NY18 and df59 were absolutely accordant. Therefore, BnaC09g20440D, BnaC09g20460D and BnaC09g20470D were not candidate genes as there were neither sequence differences in the ORFs nor expression differences between NY18 and df59.
In addition, "Ningyou 7" [39] and "ZS11" [40] reference genomes were used for the candidate gene analysis within the 120.87-kb interval harboring BnaDwf.C9. A new gene (chrC09g002459) in the "Ningyou 7" genome that was not in "Darmor-bzh" was identified (Additional file 11: Table S10). However, the gene was not expressed in either NY18 or df59. In the "ZS11" genome, BnaC09G0251200ZS was a new gene with no expression, while BnaC09G0251600ZS was a new predicted gene with a sequence length of only 240 bp, which showed neither a sequence nor expression difference between NY18 and df59 (Additional file 11: Table S10). Therefore, these genes were excluded as candidate genes for BnaDwf.C9.
BnaC09g20450D, named chrC09g002455 in "Ningyou 7" and divided into BnaC09G0251300ZS and BnaC09G0251400ZS in "ZS11", is homologous to AT2G01190 in Arabidopsis. AT2G01190 belongs to the octicosapeptide/Phox/Bem1p family of proteins and encodes a protein of unknown function ( Table 3). The PB1 domain (aa 69-167) is an important functional domain in BnaC09g20450D (Additional file 12: Figure   Fig. 4 KEGG pathway categories of differentially expressed genes between NY18 and df59 at the stem elongation stage  [41,42]. The KEGG pathway enrichment analysis showed that most of the DEGs related to signal transduction were enriched in the auxin signal transduction pathway. The RNA-seq analysis revealed that BnaC09g20450D had similar transcript levels between NY18 and df59 (Additional file 11: Table S10), while the re-sequencing analysis revealed a SNP in NY18 and df59 (Fig. 3). Subsequently, we amplified and sequenced the promoter region (2.0-kb upstream of the start codon) and the ORF of BnaC09g20450D from NY18 and df59. No polymorphism was identified in the promoter sequences. The ORF is 1887-bp in length, encoding a protein of 628 amino acids, and a single nucleotide substitution (C to T) was identified in the second exon, which changed the proline at the 585th position to serine acid (P585S) (Fig. 3). However, the SNP did not occur in PB1 domain or any other domain (Additional file 12: Figure S2). The SNP in BnaC09g20450D was also developed as PARMS SNP marker M3, which cosegregated with BnaDwf.C9 in 2536 individuals of the HO-DF population and in 2210 individuals of the ZS-DF population (Fig. 3). Therefore, we speculated that BnaC09g20450D was the most likely candidate gene of BnaDwf.C9.

Development of a PCR-based SNP marker specific for BnaDwf.C9
The SNP in the candidate gene BnaC09g20450D was targeted to develop a molecular marker. A SNP marker named BnaPHC9-SNP containing four primers (BnaM3pcr-F/BnaM3pcr-R/BnaM3pcr-Fc/ BnaM3pcr-Rt) was designed based on the 400-bp flanking sequence of the SNP. BnaPHC9-SNP was first used to amplify NY18, df59 and their F 1 , which produced a 351-bp fragment in NY18, a 179-bp fragment in df59, and both 351-bp and 179-bp fragments in the F 1 . Subsequently, 21 (tall PH), 21 (dwarfism PH) and 22 (medium PH) individuals were randomly selected from the ZS-DF F 2 population to test the SNP marker amplification. The results of the agarose gel electrophoresis showed that only a 179-bp PCR fragment was present in individuals with a dwarf PH (Fig. 5a), only a 351-bp fragment was amplified in individuals with a tall PH (Fig. 5b), and both 179-bp and 351-bp fragments were produced in individuals with a medium PH (Fig. 5c). These results suggested that BnaPHC9-SNP could produce specific PCR amplicons from different alleles of the SNP and that it can be used for the rapid identification of the BnaDwf.C9 locus, which confers the dwarfing trait, in breeding programs.

Discussion df59 is an elite genetic resource for semi-dwarf breeding
In the 1960s, the introduction of dwarfing traits into wheat and rice, combined with the application of improved cultivation methods, led to spectacular increases in grain yields, the so-called "green revolution" [12,13]. B. napus is the third most important oilseed crop, providing 13-16% of vegetable oil globally [43], and rapeseed oil has a strong potential for use in biodiesel production [44]. In China, up to 70% of the total rapeseed cultivation areas were planted with hybrid rapeseed, because they normally produce a 25% greater seed yield and have greater yield stability [37]. However, the utilization of heterosis in rapeseed also led to the PH increasing significantly, resulting in an increased risk of lodging and a decrease in the harvest index. Owing to the lack of excellent dwarf B. napus germplasm resources, there is no high-yielding semi-dwarf variety widely planted, as there is for wheat and rice. df59 is a dwarf mutant obtained from a large-scale screening of EMS-mutagenized NY18. The PH of F 1 individuals of a cross between df59 and NY18 was 126.75 ± 4.3 cm (Table 1), and the seed yield was 84% that of NY18, suggesting that the F 1 are suitable for high-density planting and mechanized harvest [37]. The average harvest index in B. napus of 0.28 [45,46], was much lower than those of wheat and rice (0.4-0.6) [45]. The harvest index of F 1 individuals was 0.41, which was significantly greater than that of NY18 (0.34), indicating that the F 1 more efficiently utilized water and soil resources ( Table 1). In addition, the PH of ZS11 was 162.6 ± 5.5 cm, while the F 1 PH of a ZS11 × df59 cross was 114.4 ± 3.9 cm (data not shown), suggesting that allele in df59 showed a dwarfing effect in multiple genetic backgrounds. In summary, df59 is an outstanding male parent for breeding new semi-dwarf hybrid varieties that can be densely planted and machine harvested.

BnaDwf.C9 is a new locus associated with plant height in B. napus
In B. napus, hundreds of QTLs associated with PH have been identified and located on all 19 chromosomes [23][24][25][26][27][28][29][30][31]. Several QTLs on the A3, A6, A9, C5 and C7 chromosomes have been fine-mapped or cloned [21,[32][33][34][35][36] [26]. Using a population containing 520 diverse rapeseed accessions, Sun et al. identified 68 loci, which were distributed over the chromosomes, except for the C8 and C9 chromosomes, significantly associated with PH using a GWAS [29]. Luo et al. and Li et al. also performed GWASs for PH, and no locus significantly associated with PH was identified on the C9 chromosome [27,28]. Thus, it appeared that there was no major QTL controlling PH on the C9 chromosome in B. napus. In the present study, BnaDwf.C9 was fine-mapped to a 120.87kb region on the C9 chromosome (Fig. 3). Although it difficult to directly compare BnaDwf.C9 with those reported QTLs owing to the lack of common markers, we still believe that BnaDwf.C9 is a potentially new locus, because it showed an obvious effect on PH.

The candidate gene may regulate plant height of df59 through the auxin-signaling pathway
The plant hormone auxin plays an essential role in most aspects of growth and developmental processes [18,19,47]. In most plants, auxin regulates the transcription of auxin-responsive genes through the well-established TIR1/AFB-Aux/IAA-ARF pathway [18,19]. The Aux/ IAAs contain four functional domains [18,19]: Domain I is a repressor domain that contains a leucine repeat EAR motif [48]; Domain II is an internal domain that contains the degron motif GWPPV; Domains III and IV share high homology and contain C-terminal regions that form a PB1 protein-protein interaction domain [49]. In B. napus, a G-to-A mutation in the GWPPV motif of domain II causes multiple phenotypic alterations, including a reduction in PH and branch angles [21]. Among the seven genes underlying the mapped interval of BnaDwf.C9 (Table 3, Additional file 11: Table S10), only BnaC09g20450D has a SNP (C-to-T mutation) between NY18 and df59. BnaC09g20450D also contains a PB1 domain (aa 69-167) (Additional file 12: Figure S2); however, the P585S substitution did not occurred in the PB1 domain or any other domain. The PB1 domain can facilitate the formation of ARF-ARF, ARF-Aux/IAA, and Aux/IAA-Aux/IAA homo-and hetero-oligomers, because the ARF and Aux/IAA proteins contain similar PB1 domains [49]. We speculate that the P585S substitution in BnaC09g20450D disrupts the normal function of the PB1 domain in an unknown way and then affects . PCR products were analyzed by electrophoresis in 2.5% agarose gel auxin signaling resulting in the dwarf stature of df59. The RNA-seq analysis also showed 81 DEGs in the auxinsignaling pathway between NY18 and df59 (Additional file 10: Table S9). Combining the results of fine-mapping and RNA-seq, we hypothesized that the candidate gene of BnaDwf.C9 regulated PH through the auxin-signaling pathway, but the molecular mechanism remains unknown and requires further study.

The developed SNP marker can facilitate the application of BnaDwf.C9 in breeding
Traditional breeding, which largely depends on the breeder's experience and subjective judgment, is a relatively resource-and time-consuming process. Markerassisted breeding is an effective and accurate approach to perform target-trait selection in breeding, because traits can be examined at any developmental stage, providing results without incurring environmental impacts [50,51]. In the present study, we determined that the dwarf B. napus mutant df59 has important potential application value in breeding. However, it is important to determine how to efficiently and extensively use the dwarfing gene in semi-dwarf breeding. To address this question, a molecular marker BnaPHC9-SNP was developed for BnaDwf.C9 based on a SNP that co-segregated in two populations, containing 4746 individuals in total ( Table 2, Additional file 4: Table S3). BnaPHC9-SNP was a dominant, allelespecific functional marker, which could amplify specific PCR products in df59 (179-bp fragment), NY18 (351-bp fragment) and their F 1 (both 351-bp and 179-bp fragments). The marker was subsequently confirmed using individuals with different PHs in the ZS-DF F 2 population (Fig. 5), suggesting that BnaPHC9-SNP can be used to select the target PH without morphological characterization, which would speed up the breeding process. To date, a variety of high-throughput SNP genotyping methods has been developed, including the TaqMan system [52] and Kompetitive allele-specific PCR. However, for most breeders, the equipment needed in these methods is prohibitory and more expensive than using normal primers. The SNP marker developed in this study requires no fluorescence-tagged probes or real-time PCR instruments, and PCR products can be correctly detected using agarose gel electrophoresis. Therefore, BnaPHC9-SNP can be widely used in marker-assisted selection of BnaDwf.C9 and speed up the breeding process.

Conclusions
Brassica napus provides not only edible vegetable oil for human consumption, but also a triglyceride source for biofuel and lubricant production. In the present study, we isolated the dwarf mutant df59 from an EMS-mutagenized NY18. The main agronomic traits in NY18, df59 and their F 1 showed that df59 is an elite genetic resource for semi-dwarf breeding. Subsequently, the combination of QTL-seq and fine-mapping revealed a candidate gene located within an interval of 120.87 kb on the C9 chromosome. The transcriptome analysis suggested that the PH of df59 was most likely influenced by a gene involved in auxin signal transduction. In addition, a comprehensive analysis revealed that BnaC09g20450D was the most likely candidate gene. Then, a molecular marker was developed based on the SNP in BnaC09g20450D. These results enrich our knowledge of the genetic architecture underlying PH in B. napus and also provide valuable resources for semi-dwarf breeding.

Plant materials
The dwarf mutant df59 was isolated from EMS-mutagenized lines of NY18. NY18 is a variety cultivated by Jiangsu Academy of Agricultural Sciences, China. Mature seeds of NY18 were mutagenized with 1.0% EMS solution (W/V, Sigma-Aldrich) at pH 7.0 phosphate buffer for 12 h at 25 °C, according to the descriptions of Li et al. with modifications [53]. The mutagenized seeds (M 1 generation) were rinsed with water for 4 h and sown in the field. Individual plants were bagged at the flowering stage. Self-pollination seeds of approximately 10,000 individual plants were harvested. Each M 2 seeds were then sown into independent line. Dwarf lines were bagged to harvest seeds for evaluation in the M 3 generation.
The NY-DF F 2 population containing 165 individual lines derived from a cross between NY18 and df59 was used for genetic inheritance and QTL-seq. Two conventional rapeseed cultivars, Holly and ZS11, were used as parental lines to develop segregating populations for finemapping of the QTLs associated with dwarf architecture. Holly is a Canadian spring variety, while ZS11 is an elite Chinese semi-winter rapeseed cultivar. The HO-DF F 2 population derived from a cross between Holly and df59 contained 2536 lines, while the ZS-DF F 2 population derived from a cross between ZS11 and df59 contained 2210 lines.

Trait measurement
NY-DF, HO-DF and ZS-DF populations, as well as their parents and the F 1 were all planted in the field of Jiangsu Academy of Agricultural Sciences, Nanjing, Jiangsu Province, China. No specific permissions were required for the field trials. The field experiments were conducted in accordance with Wang et al. with 20 plants per row and 40 cm between the rows [54]. The PH values for all the materials were measured at the mature stage.
At maturity, five open-pollinated plants each of NY18, df59 and their F 1 growing in the middle of the plot were selected. Agronomic traits and seed yield-related traits were measured in accordance with the description of Zhao et al. [55] and Wang et al. [56], including PH, biomass yield, seed yield, thousand seed weight, first effective branch height, first effective branch number, internode length, harvest index, length of main inflorescence, pod number of main inflorescence and silique number per plant. The harvest index was calculated as the ratio of seed yield to biomass yield [45]. In total, 10 well-developed siliques, which were randomly selected from the first branch adjacent to the main inflorescence, were used to determine the silique-related traits of NY18, df59 and their F 1 , including silique length, silique breadth, seed number per silique and silique volume [56]. The seed oil content was measured by nuclear magnetic resonance using standard methods and fatty acid compositions were determined by near infrared reflectance spectroscopy in accordance with Chen et al. [57].
To investigate root traits at the seeding stage, four plump seeds each of NY18, df59 and the F 1 were planted in a seed germination pouch (CYG-19LB; Pheno Trait Technology, Co., Ltd, Beijing, China) at 25 °C, and three repetitions were conducted. After 10 days, total root length (cm), root surface area (cm 2 ), root volume (cm 3 ) and number of root tips for each seeding were measured using LA-S Root Analysis software (Wanshen Ltd, Hangzhou, China).

Statistical and QTL-seq analyses
The mixed major-gene plus polygenes inheritance model of the software package SEA-G4F 2 was used to identify the inheritance of PH in the NY-DF F 2 population [38].
From the 165 lines of the NY-DF F 2 population, 16 extremely tall lines and 24 extremely dwarf lines were selected. Genomic DNA was extracted from young leaves using a Plant Genomic DNA Kit (Tiangen, Beijing, China), and T-pool and D-pool were constructed by mixing equal ratios of appropriate individual DNAs. Sequencing libraries of the two bulks and two parents were constructed, and sequence data were generated using the Illumina HiSeq ™ PE150 (Illumina, Inc; San Diego, CA, USA) platform. Both data sequencing and data analyses were performed by Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). The raw data were filtered through a series of quality controls which resulted in the removal of reads with ≥ 10% unidentified nucleotides, > 50% bases having phred quality < 5 or > 10 nt aligned to the adapter. The clean reads of each sample were aligned against the B. napus "Darmor-bzh" reference genome [22] using BWA software [58], and the SAMtools command "rmdup" was used to remove multiple read pairs [59]. Variant calling was performed for all the samples using the Unified Genotyper function in GATK software [60]. SNPs were determined using the Variant Filtration parameter in GATK. Using NY18 as the reference parent, the above two bulks SNP-indices were calculated as the proportion of reads containing SNPs not found in NY18. The Δ(SNP-index) was calculated using the formula: Δ(SNP-index) = SNP-index (D-pool) − SNP-index (T-pool). A QTL was considered a candidate associated with PH if the Δ(SNP-index) was significantly different (P < 0.05).

Fine-mapping of the QTL for plant height
The major QTL for PH was identified based on QTL-seq and named BnaDwf.C9, and SNPs underlying the confidence interval of BnaDwf.C9 were obtained. Total DNAs of HO-DF and ZS-DF F 2 individuals were extracted from fresh leaves using a modified cetyl-trimethylammonium bromide method [61]. PARMS was used to screen recombinant plants for fine-mapping among HO-DF and ZS-DF F 2 populations. The principle of PARMS SNP genotyping is similar to that of Kompetitive allelespecific PCR assays [62]. The master mix for PARMS markers was purchased from Gentides Biotech Co., Ltd. (Wuhan, China). Detailed information for conducting the qPCR-based PARMS assay was previously published by Zhang et al. [62] and Liu et al. [63].
PARMS SNP markers flanking the confidence interval of BnaDwf.C9 were screened from the two parents, the F 1 , and 10 individuals each of the extremely tall and dwarf lines of the HO-DF population. Then two polymorphic markers (M1 and M11) were used to screen the 2536 individuals of the HO-DF F 2 population (Fig. 3). Combined with phenotypes, plants containing recombinants between the two markers were selected. Furthermore, the recombinant plants were analyzed with newly developed polymorphic PARMS SNP markers (M2-M10), and BnaDwf.C9 was finally mapped to an interval between M2 and M5 (Fig. 3).
For the further fine-mapping of BnaDwf.C9, M2 and M5 were used to screen the 2210 individuals of the ZS-DF F 2 population. Using the same method as described above, plants containing recombinants between M2 and M5 were selected. Three new polymorphic PARMS SNP markers (M13-M15) together with co-segregating SNP markers (M3 and M4), were used to analyze the selected recombinant plants (Fig. 3, Additional file 5: Table S4). Finally, we narrowed down BnaDwf.C9 to a genomic region between the markers M14 and M4.

RNA library construction and sequencing
Equivalent amounts of stem tips from NY18 and df59 at the stem elongation stage were collected for RNA extraction, and two biological replicates were performed. Total RNA was extracted using a MiniBEST Plant RNA Extraction Kit (TaKaRa, Dalian, China). RNA-seq and data analysis were carried out by Novogene Bioinformatics Technology Co. Ltd. In brief, the RNA concentration was measured using a Qubit ® RNA Assay Kit and a Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). The RNA integrity was assessed using an RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Each of the two sequencing libraries for NY18 and df59 were constructed using an NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following the manufacturer's instructions, and the library quality was assessed on the Agilent Bioanalyzer 2100 system. The four libraries were sequenced on an Illumina HiSeq 2000 platform, and 100 bp paired-end reads were generated. These methods were described by Yu et al. with modifications [64].

RNA-seq data analysis
Raw RNA-seq reads were processed to remove reads containing the adapter, reads containing ploy-Ns and low-quality reads. Paired-end clean reads were aligned to the B. napus "Darmor-bzh" reference genome [22] using TopHat version 2.0.6 [65]. The read numbers mapped to each gene were counted using HTSeq version 0.6.1 [66], and reads per kilobase of exon per million reads of each gene were calculated. Differential expression analyses of NY18 and df59 (two biological replicates per sample) were performed using the DESeq R package (1.10.1). Genes with a false discovery rate < 0.005 and log 2 (fold change) ≥ 1 were declared DEGs.
The GO annotation of DEGs was performed using the GOseq R package [67], and GO terms with corrected P < 0.05 were considered to be significantly enriched terms. In addition, DEGs were submitted to the KEGG (https ://www.genom e.jp/kegg/) website, and KEGG enrichment pathways of DEGs were determined using KOBAS online analysis database [68].

Development of an allele-specific marker for traditional PCR amplification
Based on the results of the fine-mapped BnaDwf.C9, the PARMS SNP marker M3 co-segregated with the BnaDwf. C9 gene associated with PH. To confirm this SNP, the coding region of the gene that contained the SNP was amplified and sequenced from NY18 and df59. A region, including 200 bp both upstream and downstream of the SNP, was considered the target region. A SNP marker containing four primers and named BnaPHC9-SNP was designed for the allele-specific amplification of the SNP. The primers were a forward locus primer (BnaM3pcr-F: 5′-GAG AAA TAC TCC GCA ACC TACG-3′), reverse locus primer (BnaM3pcr-R: 5′-ATG TTC CGA AAC CAA CCA GAG-3′), allele primer 1 (BnaM3pcr-Fc: 5′-TAT GAA TAT GTG GAA AAT GAGC-3′) and allele primer 2 (BnaM3pcr-Rt: 5′-GCG TGT AGT ATA CCT GCT TGGA-3′). BnaM3pcr-F began amplifying from 157-bp upstream of the SNP, while BnaM3pcr-R began amplifying from 193-bp downstream of the SNP. BnaM3pcr-Fc and BnaM3pcr-Rt were upstream and downstream, respectively, with 3′-terminal bases of C and A (according to the allele of NY18 or df59), respectively. A mismatch base was introduced at the third-to-last base of the primer BnaM3pcr-Fc. All the primer oligonucleotides were synthesized by Tsingke Biological Technology Co., Ltd. (Wuhan, China).