Clubroot resistance derived from the European Brassica napus cv. ‘Tosca’ is not effective against virulent Plasmodiophora brassicae isolates from Alberta, Canada

In this study, clubroot resistance in the resynthesized European winter Brassica napus cv. ‘Tosca’ was introgressed into a Canadian spring canola line ‘11SR0099’, which was then crossed with the clubroot susceptible spring line ‘12DH0001’ to produce F1 seeds. The F1 plants were used to develop a doubled haploid (DH) mapping population. The parents and the DH lines were screened against ‘old’ pathotypes 2F, 3H, 5I, 6M and 8N of the clubroot pathogen, Plasmodiophora brassicae, as well as against the ‘new’ pathotypes 5X, 5L, 2B, 3A, 3D, 5G, 8E, 5C, 8J, 5K, 3O and 8P. Genotyping was conducted using a Brassica 15K SNP array. The clubroot screening showed that ‘Tosca, ‘11SR0099’ and the resistant DH lines were resistant to three (2F, 3H and 5I) of the five ‘old’ pathotypes and four (2B, 3O, 8E and 8P) of the 12 ‘new’ pathotypes, while being moderately resistant to the ‘old’ pathotype 8N and the ‘new’ pathotypes 3D and 5G. ‘Tosca’ was susceptible to isolates representing pathotype 3A (the most common among the ‘new’ pathotypes) as well as pathotypes 6M, 5X, 5L, 5K and 8J. Linkage analysis and QTL mapping identified a ca. 0.88–0.95 Mb genomic region on the A03 chromosome of ‘Tosca’ as conferring resistance to pathotypes 2F, 3H, 5I, 2B, 3D, 5G, 8E, 3O and 8P. The identified QTL genomic region housed the CRk, Crr3 and CRd gene(s). However, the susceptibility of ‘Tosca’ to most of the common virulent pathotypes makes it unattractive as a sole CR donor in the breeding of commercial canola varieties in western Canada.

www.nature.com/scientificreports/ its spread across Alberta 11 and into Saskatchewan and Manitoba 12,13 . The disease is managed primarily by the planting of clubroot resistant cultivars, which have allowed the continued cultivation of canola even in fields that are heavily infested by P. brassicae 14 . In recent years, however, 'new' virulent pathotypes of P. brassicae capable of overcoming this resistance have emerged in Alberta 8,15,16 ; in Saskatchewan and Manitoba, most isolates are still avirulent on clubroot resistant canola cultivars, although a virulent pathotype was recently confirmed in Manitoba 16 . Clubroot resistance in Canadian canola varieties was derived from the European winter B. napus cultivar 'Mendel' 17,18 . Since most current commercial canola varieties do not possess resistance to isolates representing the 'new' P. brassicae pathotypes, there is a need to identify and utilize additional resistance sources for development of the next generation of clubroot resistant cultivars. This task is especially daunting in Alberta, where various novel pathotypes have become widespread 8,15,16 . In this study, clubroot resistance (CR) derived from the resynthesized Swedish winter Brassica napus cv. 'Tosca' 19 was evaluated against 18 isolates representing 'old' and 'new' pathotypes of P. brassicae from Alberta. 'Tosca' was developed through many breeding cycles and hence is a stable cultivar. To understand the genetic basis of the resistance, 'Tosca' was used in genetic crosses with a Canadian spring canola line to develop a clubroot resistant spring-type canola. A doubled haploid mapping population developed from F 1 plants of the clubroot resistant spring line and a clubroot susceptible Canadian spring canola line was genotyped with a Brassica 15K SNP array, and linkage analysis and QTL mapping were conducted to identify genomic regions associated with clubroot resistance from 'Tosca' .
Genetic linkage mapping. The initial filtering steps removed 10,437 (76.1%) of the 13,714 SNP markers.
These comprised 445 (3.2%) SNP markers which failed to amplify genomic DNA in the parents, 492 (3.6%) SNP markers that were monomorphic in the parents, 2149 (15.7%) SNP markers that were monomorphic in the DH population, and 7351 (53.6%) markers that had minor-allele frequency ≤ 5% and were missing data points for > 5% in the DH population. Chi-square tests on the remaining 3277 (23.9%) SNP markers showed that 2365 (17.3%) SNP markers fit the 1:1 segregation ratio expected for a DH population (p < 0.05), 785 (5.7%) of the markers showed 'minor' segregation distortion (p < 1.67 × 10 -5 ), and 127 (0.9%) were highly distorted and hence could be discarded. Therefore, only 23.0% of the initial markers used for screening the DH population were used for linkage map construction.
Using only the 'Mendelian' markers, the linkage group lengths ranged from 13.1 cM (linkage group 4) to 189.4 cM (linkage group 20), while the total length was 2211.5 cM (Table S3). The number of 'Mendelian' markers per chromosome ranged from 5 to 204 and averaged 123.5 markers, while the marker density per cM ranged from 0.1 to 2.7 and averaged 1.1 markers per cM (Table S3). In the case of the use of the 'Mendelian' and 'distorted' markers, the linkage group lengths ranged from 14.2 cM (linkage group 5) to 188.2 cM (linkage group 6), while the total length was 2114.8 cM ( Table 3). The number of 'Mendelian' and 'distorted' markers per chromosome ranged from 20 to 289 and averaged 164.9 markers, whereas the marker density per cM ranged from 0.7 to 2.6 and averaged 1.4 markers per cM (Table 3).
Additive-effect QTL analysis. QTL analysis conducted by the CIM method with 829 Bin 'Mendelian' markers (Table S4) detected 15 coincident QTL on chromosome A03, which were significantly associated with resistance to P. brassicae pathotypes 2F, 3H, 5I, 2B, 3D, 5G, 8E, 3O and 8P. Based on the R 2 values, five of the QTL were major-effect QTL, nine were moderate-effect QTL and one was a minor-effect QTL. The peak LOD values for the QTL ranged from 6.8 to 48.1 (Table S4). The SNP markers Bn_A03_p14784764 and Bn_A03_p15704830, which were within two-LOD confidence intervals and spanned 20.6 cM (at position 56.6 to 77.2 cM), bordered the genomic region conferring resistance to the nine P. brassicae pathotypes according to the use of the 'Mendelian' markers ( Fig. 3).

Discussion
In this study, the European winter B. napus cv. 'Tosca' exhibited high clubroot resistance to the P. brassicae pathotypes 2F, 3H, 5I, 2B and 8E, which was comparable to levels previously reported in 'Mendel' 8,15,25 . In addition, 'Tosca' exhibited higher levels of resistance to pathotypes 3O and 8P, which were reported to cause increased disease (MR and S, respectively) on 'Mendel' 8,15 . In contrast, 'Tosca' was susceptible to pathotypes 6M, 5X (L-G1 and L-G2), 5L, 3A, 5K and 8J, which caused only minor or moderate disease on 'Mendel' 8,15,25 . 'Tosca' exhibited moderate resistance to pathotypes 8N, 3D and 5G, as opposed to the complete resistance shown by 'Mendel' to these pathotypes. Both 'Tosca' and 'Mendel' were susceptible to pathotype 5C. Collectively, the results of this Table 2. Chi-square tests for 1:1 segregation ratio for doubled haploid (DH) lines produced from F 1 plants obtained from the cross '11-99' (Brassica napus cv. 'Tosca' (clubroot resistant) × '12-1' (clubroot susceptible)) screened for resistance to 18 Plasmodiophora brassicae isolates under greenhouse conditions. Plasmodiophora brassicae pathotype designations are based on the Canadian Clubroot Differential (CCD) set 15 . Pathotypes 2F, 3H, 5I, 6M and 8N are single-spore isolates collected prior to the introduction of clubroot resistant canola (Strelkov et al. 20 ; Xue et al. 21 ). Pathotypes 5X (L-G1 and L-G2) and 5L are field isolates collected in 2013 (Strelkov et al. 18,20 ). Pathotypes 2B, 3A, 3D, 5C, 5G and 8E are field isolates collected in 2014 (Strelkov et al. 20,22 ). Pathotypes 5K and 8J are field isolates collected in 2015 (Strelkov et al. 20,23  www.nature.com/scientificreports/ study showed that 'Tosca' was resistant to 7 isolates, moderately resistant to 4, and susceptible to another 7 isolates. In contrast, 'Mendel' was resistant to 10 isolates, moderately resistant to 5, and susceptible to 3 isolates 8,15,25 . The differences in the resistance phenotypes of 'Tosca' and 'Mendel' to the same pathotypes in this study and in Fredua-Agyeman et al. 25 , respectively, suggest that loci controlling clubroot resistance in the two cultivars might be different. However, different loci can confer resistance to the same P. brassicae pathotypes. QTL mapping is usually carried out with markers that follow expected 'Mendelian' segregation ratios, which in a DH population should be 1:1 for resistance and susceptibility. Xu et al. (2008) 26 reported that QTL mapping could benefit from using all available ('Mendelian' + 'distorted') marker resources. Recently, Coulton et al. 27 reported that markers that showed extreme segregation distortion affected the estimation of recombination between marker pairs and hence should be discarded. In this study, Chi-square goodness of test was used to measure deviation from a 1:1 ratio. By adjusting for p-value using the Bonferroni correction, we retained an additional 785 SNP markers for the QTL analysis. The use of the additional 'distorted' markers did not result in much improvement in the QTL profiles compared with the use of only the 'Mendelian' markers at a minimum significance threshold of p < 0.05 (Figs. 2 and Fig. S1). In addition, the use of the markers with lower levels of segregation distortion did not affect the order of the genetic map. Based on the 'Mendelian' and 'distorted' markers, however, the genomic region conferring resistance to the nine P. brassicae pathotypes mapped to a narrower (15.5 cM) region compared with the use of only the 'Mendelian' markers (20.6 cM). Therefore, as previously reported, it was beneficial for QTL mapping to include both the 'Mendelian' and low 'distorted' markers.
The genomic region identified to confer resistance to P. brassicae pathotypes 2F, 3H, 5I, 2B, 3D, 5G, 8E, 3O and 8P mapped to the top half of the A03 chromosome in B. rapa and B. napus. Fredua-Agyeman et al. 18 A closer inspection of the QTL region indicated that the different-sized fragments of the A03 chromosome were responsible for the resistance to the different P. brassicae pathotypes (Fig. 3). By use of the 'Mendelian' and 'distorted' markers, the genomic region of the QTL could be partitioned into at least two CR 'hotspots' . The first CR hotspot comprised the region between the SNP markers Bn_A03_p14758285 (57.9 cM) and Bn_A03_ p15004059 (64.8 cM), while the second comprised the region between the SNP markers Bn_A03_p14968153 (67.0 cM) to Bn_A03_p15351982 (73.4 cM). The first hotspot conferred resistance to all nine (2F, 3H, 5I, 2B, 3D, 5G, 8E, 3O and 8P) aforementioned pathotypes while the second conferred resistance to pathotypes 2B, Disease resistance is a complex trait and may involve the interaction of subunits of the same gene or different genes. In this study, the QTL region contained several genes including LRR kinases. Mutation studies in Arabidopsis showed that the interaction between the LRR and the kinase domains of the ERECTA (ER) gene were required for resistance to rot caused by Plectosphaerella cucumerin 33 . In addition, the QTL region identified in this study contained transcriptional factor family proteins. In rice, the interaction of transcription activatorlike effector (TALE) proteins and transcription factor IIA small subunit was reported to determine resistance or susceptibility to bacterial leaf blight and bacterial leaf streak 34 . Epistatic interaction between the CR genes  www.nature.com/scientificreports/ CRa/CRb Kato on the A03 chromosome and the Crr1 genes on the A08 chromosome of B. rapa and B. napus was reported to confer clubroot resistance in Brassica 35 . The fact that the QTL region contained several genes suggests the possibility of various interactions within subunits of the same genes and also amongst different genes. This is further complicated by the identification of the CR genes CRk 29 , Crr3 30,31 and CRd 32 in the QTL region introgressed from 'Tosca' . Therefore, mutation studies are needed to confirm whether the CR gene(s) introgressed from 'Tosca' are three different genes or alleles of the same gene. The genome region conferring clubroot resistance derived from 'Mendel' was reported by Fredua-Agyeman and Rahman 18 to be located on the A03 chromosome at positions 24,376, 817 to 24,684,311b in B. rapa and 40,936,414 to 41,929,968 b in B. napus. The CR loci from 'Mendel' conferred resistance against the old P. brassicae pathotypes 2F, 3H, 5I, 6M and 8N from Alberta, Canada. In contrast, the CR loci derived from 'Tosca' in this study was located upstream (14,396,161,430 nt in B. rapa and 24,338,876-26,070,712 nt in B. napus) of the genomic region confering clubroot resistance in 'Mendel' . Thus, the genomic hotspot regions reported in 'Tosca' (this study) and that reported in 'Mendel' 18 are different.
In conclusion, the Swedish B. napus cv. 'Tosca' is resistant to multiple P. brassicae pathotypes, including isolates representing the 'old' pathotypes 2F, 3H and 5I as well as the 'new' pathotypes 2B, 3O, 8E and 8P. This host also exhibited moderate resistance to isolates representing the 'old' pathotype 8N and the 'new' pathotypes 3D and 5G. Unfortunately, 'Tosca' was susceptible to isolates representing the 'old' pathotype 6M and the 'new' pathotypes 5X (L-G1 and L-G2), 5L, 3A, 5K and 8J. This is the first report on the genomic loci controlling clubroot resistance in the Brassica napus cv. 'Tosca". The resistance was shown to be different from the clubroot resistance derived from 'Mendel' . The increased clubroot severity on 'Tosca' , especially in response to pathotypes 3A and 3D, which constitute the bulk of the virulent pathotypes (note: pathotype 3H is still most prevalent of all pathotypes) in Alberta, makes the cultivar unattractive as the sole CR donor in the breeding of commercial canola varieties in Canada. However, the CRk, Crr3 and CRd gene(s) present in 'Tosca' could be stacked with other CR genes present in additional resistance resources such as 'Mendel' , ECD 02 and ECD 04.

Materials and methods
Plant materials. One hundred sixteen doubled haploid (DH) lines obtained from F 1 plants of the cross '11SR0099' (clubroot resistant) × '12DH0001' (clubroot susceptible) were used as the mapping population. The CR parent '11SR0099' is a spring-type canola line derived from a spring canola × winter canola cv. 'Tosca' cross, while the CS parent '12DH0001' is a spring-type canola line with good agronomy and quality characteristics. Seeds of the DH parents and lines were provided by Corteva AgriScience (Caledon, ON, Canada), while seeds of a 'Tosca' , used as a resistant (negative) control in the inoculation experiments, were obtained from Prof. Ann-Charlotte Wallenhammar (Swedish University of Agricultural Sciences, Skara, Sweden). The European Clubroot Differential (ECD) 04 (B. rapa subsp. rapifera), which exhibits broad-spectrum resistance to many Canadian isolates of P. brassicae 25 , was included as resistant (negative) control in the inoculation experiments, while ECD 05 (B. rapa var. pekinensis 'Granaat') 36 and B. napus cv. 'Westar' 25 were included as susceptible (positive) controls.  www.nature.com/scientificreports/ 17 pathotypes) were maintained as frozen (− 20 °C) root galls of the universally susceptible host ECD 05 until needed. Clubroot inoculum was prepared following Fredua-Agyeman et al. 25 by macerating root galls using a variable-speed blender, filtering the spore suspension through three layers of cheesecloth, and adjusting the final resting spore concentration to 1 × 10 7 spores/mL with sterile distilled water. Each batch of inoculum was stored at 4 °C and used within 24 h after preparation. Statistical analyses of phenotypic data. Statistical analyses of the disease severity data for the individual experiments and the combined data were conducted with SAS v. 9.4 (SAS Institute, United States) as described by Fredua-Agyeman et al. 25,38 . In brief, the PROC CORR function was used to determine the correlation among the mean ID values for each DH line, parents and checks for each pathotype in the three experiments. Broad sense heritability (H 2 ), which is the ratio of total genetic variance to phenotypic variance, was estimated from variance components from Analysis of Variance (ANOVA) 41,42 . The PROC MEANS function was used to calculate the mean ID, standard error of the mean (SEM), minimum and maximum ID for each genotype and isolate investigated. The PROC FREQ function was used to count the number of accessions that were resistant (Mean ID ± SEM ≤ 30%), moderately resistant (30 < ID ± SEM ≤ 50%) or susceptible (ID ± SEM > 50%) to each P. brassicae isolate based on the combined data, while SigmaPlot (SYSTAT Software Inc., San Jose, California, USA) was used to create histograms. The Shapiro-Wilk test was used to test for normality in the phenotypic data 43 . The 1R: 1S ratio suggested for segregation in a DH population was determined through Chi-Square goodness of fit tests (χ 2 ) at p ≤ 0.05 for each of the 18 P. brassicae isolates. Differences in the mean ID values of all DH lines to pairs of the 18 P. brassicae isolates were compared with Tukey's test at p ≤ 0.05.

Evaluation of DH lines for clubroot resistance.
Genotyping with SNP markers. The parental lines and the DH population were genotyped using the Brassica 15 K array, which contained 13,714 SNPs, at SGS TraitGenetics GmbH, Gatersleben, Germany, according to the manufacturer's protocols 44 . The software TASSEL v5.2.2.5 45 was used to perform SNP filtering by deleting failed SNP reactions, setting minor allele frequency (MAF) to ≤ 0.05 and removing markers missing data for > 5% of the accessions. Segregation distortion was determined through a χ 2 test for goodness-of-fit for the 1:1 ratio expected for a DH population. A minimum significance threshold of p < 0.05 was used for markers that followed expected 'Mendelian' ratio, while adjusted Bonferroni correction p-values (α/n, where α = level of significance, n = number of markers) were used to select markers with 'minor' segregation distortion 46 . Markers that showed extreme segregation distortion were discarded.
Construction of genetic linkage maps. Two 'draft' linkage maps were constructed using the minimum spanning tree map (MSTMap) program with the following parameters: logarithm of odds (LOD) value of 10.0, a maximum distance between markers of 15 cM and Kosambi mapping function 47,48 . The first linkage map was constructed only with markers that fit the expected 1: 1 Mendelian ratio, while the second linkage map was constructed using markers that fit the 1:1 Mendelian ratio expected for a DH population and markers that showed 'minor' segregation distortion. In both cases, markers that mapped to the same position were placed in the same bin, with only one of the markers retained for linkage analysis. MAPMAKER/EXP v. 3.0b 49 was run at a logarithm of odds (LOD) score ≥ 3.0 and recombination fraction (ϴ) value ≤ 0.40 to 'refine' the marker order obtained by the MSTMap software. Recombination fractions were converted to centiMorgans (cM) using the ID (%) = (n × 0 + n × 1 + n × 2 + n × 3) N × 3 × 100  47 . Linkage groups were assigned to chromosomes based on the SNP sequence information provided by SGS TraitGenetics GmbH.
Additive effect QTL mapping. Quantitative trait loci (QTL) mapping was conducted separately with only the 'Mendelian' markers and with 'Mendelian' and 'distorted' markers together. The QTL analyses were carried out by composite interval mapping (CIM) 50 using WinQTL Cartographer v. 2.5 51 . The program was run at a walking speed of 1 cM and with the following settings: forward-backward regression method, a window size 2.0 cM, five background markers as cofactors, 1000 permutations and p < 0.05. The significance level required to declare a QTL was set at LOD ≥ 5.0. Locations of putative QTL were estimated based on two-LOD support intervals for a 99% confidence interval (CI) 52 .
The QTL designations were of the order genus (1 letter), species (1 letter), genome (1 letter), chromosome number (1 letter), pathotype name (3 letters), closest published gene(s) (3-6 letters) and QTL number (2 letters) 35,53 . The percentage of phenotypic variation (R 2 ) explained by each QTL was calculated. QTL were arbitrary assigned as major-, moderate-or minor-effect QTL when the R 2 explained > 50%, 25-50% or < 25% of the phenotypic variation, respectively. The additive effects of each QTL were calculated by deducting the phenotypic average of all individuals with the susceptible DH parent allele from all individuals with the resistant DH parent allele.
Identification of candidate genes. The physical positions of the SNP markers in the QTL CI region were mapped to the B. napus, B. rapa and B. oleracea reference genomes deposited in the National Center for Biotechnology Information (NCBI) GenBank database (www. ncbi. nlm. nih. gov). Candidate genes present in the QTL two-LOD confidence interval were identified by BlastN searches (E-value ≤ E−20, minimum identity of sequence ≥ 95%) of the three Brassica genome sequences.
Ethical approval. On behalf of all authors, the corresponding author states that the collected seeds and all conducted experiments in this study complied with relevant guidelines/regulations of the University of Alberta, the Canada Food Inspection Agency and International Treaty for Plant Genetic Resources guidelines and legislation.

Data availability
The datasets generated during the current study are available in the manuscript or the supplementary materials.