The association population of rapeseed including 472 accessions, which were collected from 23 countries worldwide of four continents, has been previously reported (Li et al., 2014). All the seeds of these inbred lines were conserved in the National Mid-term Gene Bank for Oil Crops in Wuhan, China. To carry out the artificial inoculation in greenhouse (GH), the infected roots of rapeseed were collected from the infected field (IF) of Dangyang, China, where the pathogen was reported as pathotype 4 based on Williams classification (Ren et al., 2012).

All of the 472 accessions were planted in field infected seriously by P. brassicae at Dangyang, China on October 2013. The field experiments followed a randomized complete block design with two rows and three replications. Before sowing, we partitioned compartments with width of 2.0 m in consideration of field drainage and soil fertility factors. The plants were sown in strip with a 0.3 m spacing between rows, singling 12 plants each row after sowing for 1 month. Infected rapeseed roots collected from Dangyang were stored at -20°C, which were used for artificial inoculation trials in GH. The artificial inoculation of the population was performed at Shenyang on October of 2015. Firstly, nutrient soil was dried and powdered, while the disease roots were unfrozen and decomposed for one night. Then we squeezed the disease roots with distilled water to make P. brassicae suspension. The suspension was mixed with the dried nutrient soil and adjusted to 10⁸ resting spores per gram (g) dry soil. The soil humidity which reached to the extent of kneading with clumps and touching with scatter was advisable, and then was sealed to ferment for 2 days in dark. The hole tray was filled with nutritional soil and then irrigated enough water. About 3 g fermentative soil per hole was placed, and sowed two seeds in the fermentative soil. The artificial inoculation followed a randomized complete block design with three replications, and inoculated 24 seeds for each replication. The plants grew under a 16 h L/8 h D photoperiod at an average temperature of 20–25°C.

At flowering time of IF trials and 50 days after inoculation in greenhouse, plants were all taken out with roots and the disease severity was assessed using a 0–3 scale as follows, 0 = no clubs, 1 = a few small clubs, 2 = moderate clubbing, and 3 = severe clubbing (Kuginuki et al., 1999). A disease index (DI) was calculated according to the formula DI = [(n₁+2n₂+3n₃)/(N_T × 3)] × 100, where n₁ to n₃ were the number of plants with different disease severity of 1–3 scale and N_T represented the total number of identified plants, respectively. An incidence rate (IR) was calculated with a formula IR (%) = [(n₁+n₂+n₃)/N_T] × 100. The broad-sense heritability was calculated as h²= σ²_g/(σ²_g+σ²_ge/n+σ²_e/nr). Where σ²_g is the genetic variance, σ²_ge is the interaction variance of the genotype with environment, σ²_e is the error variance, n is the number of environments and r is the number of replications. The estimates of σ²_g, σ²_ge, and σ²_e were obtained from the analysis of variance (ANOVA) procedure in SPASS 20.0 and the frequency distribution of phenotype was carried out using Excel 2010. The multiple comparisons of phenotypic average value from different origins were tested with package of SAS 9.2 at a significance level of P < 0.05. Also, the difference analysis of genotypic effects with containing different No. of favorable SNPs was performed online¹.

The young leaves from the representative plants of all accessions were collected, and extracted their genomic DNAs using a modified cetyltrimethyl ammonium bromide (CTAB, Murray and Thompson, 1980) method. SNP genotyping of 472 accessions was performed using Brassica 60K Illumina Infinium^® HD Assay SNP arrays by Emei Tongde, Co. (Beijing) according to the manufacturer’s protocol². The quality control process of SNP array genotyping data was the same to our previous study (Li et al., 2014). The LD was calculated using the software TASSEL3.0 (Bradbury et al., 2007). The principal component analysis (PCA), population structure and kinship were analyzed by Li et al. (2014). To confirm the physical position of each SNP, the probe sequences of 26,841 high-quality SNPs previously selected (Li et al., 2014) were used to perform a BlastN search against the B. napus (Darmor-bzh) reference genome by Li et al. (2016). Only the top blast hits were considered to be mapped in the genome with an E-value threshold of 1E–15, and blast matches to multiple loci with the same top E-value were not considered to be mapped successfully. Totally, 19,945 SNPs were assigned to B. napus chromosomes.

To evaluate the effects of population structure (Q, PC), the trait-SNP association analysis was performed using TASSEL 4.0 with six models: (i) GLM, controlling for nothing, (ii) Q model, controlling for Q, (iii) PCA model, controlling for PC, (iv) K model, controlling for K, (v) Q+K model, controlling for both Q and K, (vi) P+K model, controlling for both PC and K. The GLM, Q, K, and PCA models were performed using a general linear model with the following equation: y = Xα+ e. The Q+K and P+K models were performed using a mixed linear model with the following equation: y = Xα+Kμ+e. In above equations, y represented phenotype, X represented genotype, α was a vector containing fixed effects, K was the relative kinship matrix, μ was a vector of random additive genetic effects, and e was the unobserved vector of random residual. The distribution of observed -log₁₀(p) for each SNP from marker-trait associations was compared with the expected distribution in a quantile-quantile plot (Q-Q plot). The significance of associations between SNPs and traits was based on a threshold P < 5.01 × 10^-5 [P = 1/N, where N = the number of SNPs used, -log₁₀(1/19,945) = 4.3]. Meanwhile, to avoid ignoring the effects of minor loci, the threshold P was lowered to 2/N = 1.00 × 10^-4. Q-Q plots and Manhattan plots were drawn using the R package qqman³.

To compare the QTLs detected in this study with previous studies on CR, the sequences of QTLs/genes which have been reported in previous studies were collected and performed BlastN or e-PCR on the B. napus reference genome (Darmor-bzh) with a threshold of 1E–10. By that, we could ensure the homoeologous fragments’ physical positions of QTLs/genes reported in previous studies on B. napus genome.

It was regarded as a LD block that the region containing all SNPs which r² > 0.4 with the most significant SNP. To ascertain regions that potential candidate genes of each QTL identified in this study maybe located, we considered the flanking markers that outside and adjacent to the LD block were considered as the candidate region’s boundary. To further identify the candidate genes and the participated pathways of the identified loci, the gene function annotation and enrichment analysis were carried out. All B. napus genes which located in the candidate regions were searched against the NCBI non-redundant (Nr) protein database using BlastP with an E-value ≤ 1E–05. The gene ontology (GO) terms of the B. napus genes were annotated by merging the Blast2GO and InterPro annotation results (Wu et al., 2016). GO enrichment analysis provided all of the GO terms that were significantly enriched in candidate genes compared with the genome background using Blast2GO with a false discovery rate (FDR) ≤ 0.01. The published differential expressed genes (DEGs) information in B. rapa were downloaded (Chen et al., 2015), and searched the homoeologous genes in B. napus with a threshold of 1E–50.

The 472 rapeseed accessions were planted in two environments of IF in Dangyang at the year of 2013–2014, and GH in Shenyang at the year of 2015, respectively. DI and IR were used to evaluate extent of CR. The frequency distribution diagrams of these two indicators in the two environments all showed continuous distributions (Supplementary Figure S1), which indicated that DI and IR were both controlled by multiple loci. Extensive variations of these two indicators were observed among the population in both two environments (Table 1). In IF environment, DI ranged from 9.26 to 75.00, with an average value of 31.39 ± 5.92 (± SE); IR ranged from 18.18 to 100% with an average value of 48.11 ± 8.10% (± SE). The broad-sense heritability (h²) of DI and IR in IF were 64.8 and 55.0%, respectively (Table 1). In GH environment, DI was from 36.11 to 100.00, with an average value of 84.77 ± 9.67 (±SE); IR was from 58.26 to 100%, with an average value of 96.49 ± 7.64% (± SE). Moreover, the h² of DI and IR in GH were 78.2 and 69.4%, respectively (Table 1). The results thus indicated that the stability of CR is high, especially in a control environment. The correlation analysis of these two indicators between the two environments showed that the correlation reached to highly significant level, but the value of which (0.24 and 0.19) is too low (Supplementary Table S1). The results indicated that genetic factor and environment factor both played an important role in CR.

To assess the effects of origin including different continent (Europe, Asia, America, and Oceanica), different country (China and abroad), ecotype (spring, semi-winter, and winter type), and breeding eras in China (the years of 1950–1970, 1970–1980, 1980–1990, 1990–2000, and 2000–2011) on CR in the population, the accessions were classified by the above categories, respectively (Supplementary Figure S2). The DI and IR were both compared among different categories, and the results showed that the accessions from China were more resistant than that of from abroad both in DI and IR of two environments. Nevertheless, the absolute difference of these two indictors between China and abroad were very small (Supplementary Figure S2A). The similar tendency was emerged in the comparison among Europe, Asia, America, and Oceanica. The accessions from Asia and Europe were more resistant (Supplementary Figure S2B). Moreover, the semi-winter ecotype accession was the most resistant among the different ecotypes (Supplementary Figure S2C). The resistance of accessions from different breeding eras had significant difference (Supplementary Figure S2D), which indicated that artificial selection could impact the CR significantly. Therefore, we can further reinforce the resistance breeding in B. napus.

In order to select suitable model for association mapping, GWAS was performed using four general liner models (GLM, Q, K, and P) and two mixed liner models (Q+K, and P+K) for IF-DI (DI in IF), IF-IR (IR in IF), GH-DI (DI in GH), and GH-IR (IR in GH), respectively. As the Q-Q plot showed, the distribution of observed -log₁₀(p) from P+K model was the closest to the expected distribution for all four traits, which led to a low level of false-positive signals (Supplementary Figure S3). Therefore, the association signals were identified with P+K model subsequently.

To avoid ignoring the micro effective locus, significant associated SNP (sSNP, P < 5.01 × 10^-5) and potential association SNP (pSNP, 5.01 × 10^-5 < P < 1.00 × 10^-4) were introduced in the study. If a confidence interval of a locus included one or more sSNPs, we called it significant association locus. Similarly, if a confidence interval of a locus did not have sSNP, and only included two or more pSNPs, we called it potential association locus. In total, seven significant association loci and two potential association loci were identified with these two indicators (Figure 1).

For IF-DI, one significant association locus on BnC03 (named MCR-C3) was identified, with a peak SNP (highest significant) Bn-scaff_17521_1-p419499 which contributed to 4.81% of phenotypic variance. One potential association locus on BnC09 (named MCR-C9) was detected, with a peak SNP Bn-scaff_15576_1-p660538 which explained 4.21% of phenotypic variance (Table 2). For IF-IR, one significant association locus on BnA04 (named MCR-A04) was mapped, with a peak SNP Bn-A04-p16156157 which contributed 4.72% to phenotypic variance. It was noteworthy that MCR-C3 in IF-DI was also detected in IF-IR. For GH-DI, one significant association locus SCR-C6 was detected on BnC06 with a peak SNP Bn-scaff_16064_1-p26703 which explained 5.26% of phenotypic variance. Another significant association locus SCR-A10a was identified on BnA10 with a peak SNP Bn-A10-p16087066 which accounted for 4.71% of phenotypic variance. For GH-IR, three significant association loci were identified in total. Among that two loci were both on BnC04 (named SCR-C4a and SCR-C4b), which explained 5.72 and 4.54% of phenotypic variance, respectively. The rest one named SCR-A10b was on BnA10, with a peak SNP Bn-A10-p3966740 which explained 4.87% of phenotypic variance. In addition, one potential association locus SCR-C3 was mapped on BnC03 with a peak SNP Bn-scaff_18559_1-p166394 which accounted for 4.12% of phenotypic variance (Table 2).

To understand the effects of allelic variations on CR in rapeseed, the combined effects of DI of different QTLs were analyzed in IF and GH, respectively (Figure 2). The genotype and phenotypic variance (R²) of each QTL was substituted by that of each peak marker of the corresponding QTL. On this basis, 472 accessions were grouped into three classes which contained 1, 2, and 3 favorable alleles in IF, respectively (Figure 2A). Because of too less accession contained in one group if all accessions were grouped by containing 1–6 favorable alleles in GH, all these six loci were merged and divided into three groups which contained 1–2, 3–4, and 5–6 favorable alleles to statistical analyze the DI in GH (Figure 2B). The results indicated that the more favorable alleles (resistant alleles) had, the lower DI (more resistance to P. brassicae) was demonstrated (Figure 2). In IF, the DI of the accessions held three resistant alleles can be reduced by 20.62% compared with that of holding one resistant allele. Similarly, the DI of the accessions held 5–6 resistant alleles can be reduced by 10.56% compared with that of only holding 1–2 resistant alleles in GH. The above results revealed that the CR was mainly controlled by additive effect. Therefore, the CR for B. napus can be improved by polymerization of all identified resistant alleles.

To compare the GWAS results in this study with the previous reported CR QTLs/genes, the sequences of CR genes and markers of QTLs related to CR were collected in B. napus, B. rapa, and B. oleracea. And then the blastn or e-PCR with a threshold of 1E–10 was performed to B. napus reference genome (Darmor-bzh) to search their homoeologous regions in B. napus (Figure 3). In B. rapa, two homoeologous regions of CRa were identified on BnA03 (22,864,716–22,877,171 bp) and BnC07 (38,863,628–38,883,622 bp), respectively (Figure 3). Three homoeologous regions of Crr1 were found on BnA08 (9,456,084–9,467,947 bp), BnC03 (54,084,701–54,109,095), and BnC07 (38,863,628–38,869,983 bp), respectively (Figure 3). One homoeologous region of CRb was identified on BnA03 (21,682,961–22,238,961 bp). Two homoeologous regions of Rcr1 were found on BnA03 (22,797,049–23,009,129 bp) and BnC07 (36,229,328–39,036,510 bp) in B. napus. Homoeologous regions of other QTLs in B. rapa were not identified, because of their limited marker information or the physical positions of forward and reverse markers linked with the QTLs were not on the same chromosome in B. napus genome. Because of the limited reports on CR studies in B. oleracea, only three QTLs CRQTL-YC, CRQTL-GN_1, and CRQTL-GN_2 were searched. Two homoeologous regions of CRQTL-YC were identified on BnA03 and BnC03 respectively. Also, two homoeologous regions of CRQTL-GN_1 were detected on BnA02 (22,576,577–22,970,418 bp) and BnC02 (42,737,365–43,546,906 bp). And one homoeologous region of CRQTL-GN_2 was mapped on BnC03 (1,185,066–2,835,468 bp; Figure 3).

Moreover, it was attempted to ensure the QTL physical position in B. napus according to the providing markers information in previous reports. Markers sequences of five QTLs reported on previous study in B. napus were searched and the physical position of PbBn-01:07-1 and PbBn-korp-4 were confirmed successfully. The physical positions of these two QTLs, BnA03 (5,386,730–5,386,845 bp) and BnC06 (27,067,129–27,067,246 bp), were proved to be accordance with the previous reports. However, the physical positions of other three QTLs were not accordance with the previous studies. For example, the QTL PbBn-a-1 located on BnA08 was mapped on BnC03, and the PbBn-01:07-3 located on BnC03 was mapped on BnA03 in this study (Figure 3).

It was found that the three QTLs identified in IF in our study were all novel, in which regions there were no homoeologous regions of CR QTLs reported in previous studies (Figure 3). The QTL SCR-C3 identified in GH was not located on the confidence intervals of one homoeolog of Crr1, but was close to it. We considered it as an old CR locus. The same situation was on SCR-C6, which closed to the physical position of PbBn-korp-4. The other four QTLs identified in this study were novel. So the GWAS had stronger ability to detect new loci on a given trait.

To predict the candidate genes of CR QTLs identified in this study, the candidate region (confidence interval) for each QTL was confirmed. There were 542 predicted genes in all nine candidate regions, 471 genes of which had functional annotations in total. To gain insights into the functionality of above genes, we performed GO enrichment analysis using Blast2GO (Conesa et al., 2005). The result showed that eight pathways of the molecular function related to ADP binding and methylthiopropyl-desulfoglucosinolate sulfotransferase activity were enriched (Supplementary Figure S4), the pathways of which were highly significant with TIR-NBS-LRR genes. It was interesting that two TIR-NBS gene clusters held 28 TIR-NBS genes which located in SCR-C6 (10 TIR-NBS genes) and MCR-C9 (18 TIR-NBS genes; Table 3) were participated in the function of ADP binding. The results indicated that the TIR-NBS gene family may associated with CR.

In order to get more evidence to predict candidate genes, the homoeologous genes of DEGs from the transcriptome data after inoculation P. brassicae in B. rapa were identified in B. napus genome (Figure 4). For example, there were five SNPs involved in the MCR-C3 LD block, the candidate gene region was 21.72–21.92 Mb (194.6 Kb) in BnC03, where 44 genes were included, and four of which were the homoeologous genes of four DEGs identified in B. rapa (Figure 4A). Also, seven SNPs were involved in the LD block for SCR-C3, which the candidate region was 57.86–58.10 Mb (339.3 Kb) in BnC03, and 16 genes were predicted in this region. Only one homoeologous gene BnaC03g68270D of DEGs corresponding in B. rapa was searched (Figure 4A). Similarly, 9, 1, 2, 1, 6, 2, and 4 candidate genes obtained from the B. rapa homoeologous genes were identified in the candidate regions of MCR-A4, SCR-A10a, SCR-A10b, SCR-C4a, SCR-C4b, SCR-C6, and MCR-C9, respectively (Figure 4; Supplementary Table S2; Supplementary Figure S5). Overall, the 30 homoeologous genes of DEGs from B. rapa were likely candidate genes as these nine QTLs (Supplementary Table S2). However, the more evidence needs to be obtained by functional analysis of these genes.

The virulence of pathogen and resistance of host plant also changed along with the continuous game between them. There has been an abundant of research on characterizing the virulence of P. brassicae (Xue et al., 2008; Cao et al., 2009), and many races (pathotypes) have been identified by two identification system (Williams classification and European Clubroot Differential). At present, most CR resources identified in Brassica were pathotype specific. However, the resistance to single race controlled by single gene was unstable, which was continuously eroded by the pathogen’s adaptation. The above dynamics and instability made the genetic regulation of host plant became more complexity.

Until now, there is no certain conclusion about genetic regulation mode of CR. Genetic analysis showed that both qualitative and quantitative manners existed in the important Brassica crops due to the different resistant sources. In B. rapa, CR was controlled in a qualitative and quantitative manners both, whereas it was quantitative traits in B. oleracea (Piao et al., 2009; Lee et al., 2016). Genetic studies on CR of rapeseed were most based on the sources, of which the resistance derived from B. rapa (Zhang et al., 2015). Also, the inheritance of CR in rapeseed was governed by one or two single independent dominant genes (Piao et al., 2009). The difference of genetic regulation mode in CR among different species may be related to the various races. So enhancing the study to physiological race, including developing of new identification system, is the basis to perform intensive study.

There were three and six QTLs were identified in IF and GH, respectively (Table 2). But the same QTL between two environments were not detected, which may because that we performed the evaluation at the seedling stage in GH while at the harvest stage in IF. There were maybe different genetic mechanisms between vegetable stage and adult stage in rapeseed. It also gave us an ideal that it can improve the rapeseed CR both in vegetable stage and adult stage by polymerization of all identified CR loci/genes by MAS. The identified CR QTLs were showed additive effect (Figure 2), which provided theoretical basis for putting forward the resistant breeding strategy of polymerization of multiple loci, similarly. The SNPs detected in this study could be transformed to PCR markers (Yang et al., 2012), which could prove and consummate the method of screening the resistant germplasms resources and also provide strong and powerful measures for MAS. Meanwhile, we can pyramid different loci that resistant to different pathotypes for breeding to increase the resistant durability (Matsumoto et al., 2012; Tomita et al., 2013).

The results of integration of all the CR QTLs/genes on B. napus also gave us revelations that the closely relate species of B. napus (e.g., B. rapa, B. oleracea, and B. juncea) can be used for CR improvement through distant hybridization between these species (Figure 3). Therefore, we will have two approaches to improve the rapeseed CR. Firstly, all QTLs/genes which come from B. napus and closely relate species can be polymerized by MAS or interspecific hybridization (Yu et al., 2016); Secondly, we can get more resistant B. rapa and B. oleracea by MAS and phenotypic identification, and then the highly resistance rapeseed can be obtained by artificially synthesized (Diederichsen et al., 2006).

In this study, 6 of 9 identified QTLs were on C-genome, and most of them were novel, which illustrated that the C-genome (B. oleracea, B. carinata) had more potential on CR improvement (Figure 3). Researchers also found that there were more resistant loci in C-genome than A-genome of B. napus in the study of QTL mapping of resistance to Sclerotinia sclerotiorum (Wu et al., 2013; Wei et al., 2016). Wang et al. (2014) thought that breeding selection loci of C-genome were more than A-genome variation analysis of rapeseed core germplasms, which further verified that our research results had crucial contribution to resistant breeding in B. napus. Therefore, it is necessary for us to expand the resistance sources to discover more resistance genes/loci.

The CR QTL mapping results in B. rapa showed that 6 of 11 CR genes/QTLs were on BrA03, which indicated that there were CR hotspots existed in B. rapa. Until now, the genome sequences of B. rapa, B. oleracea, and B. napus have been sequenced which will accelerate their genetic improvements. All reported QTLs/loci information (linkage markers, candidate genes, and positions information) related to CR were collected in these three species (Figure 3). Most of the reported CR QTLs/genes were mapped on B. napus genome through bioinformatic analyses. The integration results illustrated the CR hotspots also existed in B. napus, which were the regions of BnA03 (21–23 Mb), BnC03 (1–2 Mb), and BnC07 (36–39 Mb; Figure 3). It was interesting that the top of BnA03 was homoeologous with the top of BnC03, and the bottom of BnA03 was homoeologous with the bottom of BnC07 (Chalhoub et al., 2014). Therefore, the evolution and origin of BnA03, and the genes related to disease-resistance on BnA03 are worth to deeply study. It was a pity that we did not detect any loci in the above hotspot regions of CR. The reason may be that the above hotspot regions were summarized from the population based on bi-parents in B. rapa. These hotspots may belong to rare variation sites which were difficult to detect in our natural population, or they were peculiar in B. rapa. Another reason maybe was that the previous mapping results were most aimed at the race/pathotype 3, the results of this study was based on the pathotype 4.

So far, there were more than 70 R genes cloned from plants (Liu et al., 2007). Most of the cloned R genes were NBS-LRR family (Belkhadir et al., 2004; McHale et al., 2006). It was interesting that two R genes controlled CR in B. rapa were all TIR-NBS-LRR genes, which encoded large and abundant proteins involved in the detection of diverse pathogens. Similarly, the enrichment analysis also found eight pathways which were highly significant with TIR-NBS-LRR genes (Supplementary Figure S4), and there were two TIR-NBS gene clusters in our candidate regions (Table 3). NBS-LRR proteins could recognize the specialized pathogen effectors of ETI (also called avirulence proteins). Therefore, TIR-NBS-LRR gene family is most likely to play an important role in the process of the clubroot disease, which needs to deeply study in Brassica. Besides that, plant hormones, mainly salicylic acid (SA), jasmonic acid (JA), auxins and cytokinins, also played a role in compatible interactions between Arabidopsis and P. brassicae (Siemens et al., 2006; Lemarié et al., 2015). It also enlightens us on prediction candidate genes of above pathways and understanding the resistant molecular mechanism. Up to now, only one paper was reported on transcriptome analysis of CR in Brassica (Chen et al., 2015). It was a pity that two NILs of CRb were used for RNA-Seq in B. rapa, not two accessions or varieties. Therefore, these 30 candidate genes identified through the data from RNA-Seq in B. rapa just only provided clues for candidate genes prediction and understanding the molecular regulation of CR in B. napus. However, it needs to carry out more experiments to obtained more evidence to identify or evidence these candidate genes. For example, the RNA-Seq analysis can be carried out by using of some CR and susceptible lines from the natural population in this study.