QTL Mapping and Inheritance of Clubroot Resistance Genes Derived From Brassica rapa subsp. rapifera (ECD 02) Reveals Resistance Loci and Distorted Segregation Ratios in Two F2 Populations of Different Crosses

In this study, Brassica rapa subsp. rapifera (ECD 02) which exhibits broad-spectrum resistance to many Canadian Plasmodiophora brassicae isolates was crossed with two clubroot-susceptible B. rapa accessions to produce two F2 populations. The F2 plants were screened against P. brassicae pathotypes 3H, 5X, and 5G. The Chi-square goodness of fit test showed that the vast majority (≈75%) of the crosses that produced the F2 populations showed segregation ratios of 9R:7S, 7R:9S, 13R:3S, 3R:13S, 5R:11S, 11R:5S, and 1R:15S. These were modifications of the 15R:1S ratio expected for the inheritance of two dominant major clubroot resistance (CR) genes from ECD 02. The distorted segregation ratios suggest that the two resistance genes are on different chromosomes and that two genes interact in an epistatic manner to confer resistance. Genotyping was conducted with 144 PCR-based markers in the two F2 populations. Linkage and QTL analysis with the polymorphic markers identified two QTLs on chromosome A03 to be associated with resistance to P. brassicae pathotypes 5X and 5G in Popl#1 while only the second QTL on chromosome A03 was associated with resistance to pathotypes 5X and 5G in Popl#2. The QTLs clustered in genomic regions on the A03 chromosome of B. rapa where the CRa/CRbKato gene(s) are mapped. In addition, the Crr1 gene on the A08 chromosome of B. rapa was detected in the two F2 populations. Therefore, the phenotypic and molecular data confirm the existence of two CR genes in ECD 02. This is the first study that shows that major dominant genes in Brassica interact in a non-additive manner to confer resistance to different P. brassicae pathotypes. Key Message: This study provides knowledge on the inheritance and type of gene action for clubroot resistance derived from Brassica rapa subsp. rapifera (ECD 02). The results indicated that duplicate recessive and recessive suppression epistatic interactions, digenic additivity and complementary gene action between the CRa/CRbKato gene(s) on the A03 and the Crr1 gene on the A08 chromosome of B. rapa controlled clubroot resistance to P. brassicae pathotypes 3H, 5X and 5G.

In this study, Brassica rapa subsp. rapifera (ECD 02) which exhibits broad-spectrum resistance to many Canadian Plasmodiophora brassicae isolates was crossed with two clubroot-susceptible B. rapa accessions to produce two F 2 populations. The F 2 plants were screened against P. brassicae pathotypes 3H, 5X, and 5G. The Chisquare goodness of fit test showed that the vast majority (≈75%) of the crosses that produced the F 2 populations showed segregation ratios of 9R:7S, 7R:9S, 13R:3S, 3R:13S, 5R:11S, 11R:5S, and 1R:15S. These were modifications of the 15R:1S ratio expected for the inheritance of two dominant major clubroot resistance (CR) genes from ECD 02. The distorted segregation ratios suggest that the two resistance genes are on different chromosomes and that two genes interact in an epistatic manner to confer resistance. Genotyping was conducted with 144 PCR-based markers in the two F 2 populations. Linkage and QTL analysis with the polymorphic markers identified two QTLs on chromosome A03 to be associated with resistance to P. brassicae pathotypes 5X and 5G in Popl#1 while only the second QTL on chromosome A03 was associated with resistance to pathotypes 5X and 5G in Popl#2. The QTLs clustered in genomic regions on the A03 chromosome of B. rapa where the CRa/CRb Kato gene(s) are mapped. In addition, the Crr1 gene on the A08 chromosome of B. rapa was detected in the two F 2 populations. Therefore, the phenotypic and molecular data confirm the existence of two CR genes in ECD 02. This is the first study that shows that major dominant genes in Brassica interact in a non-additive manner to confer resistance to different P. brassicae pathotypes.

INTRODUCTION
Clubroot is a soil-borne disease of the Brassicaceae caused by the obligate parasite Plasmodiophora brassicae. Disease development is associated with the formation of large galls on the roots of susceptible plants, which interfere with water and nutrient uptake and lead to significant yield losses in Brassica crops (Hwang et al., 2012;Dixon, 2014). Yield losses of 20-100% have been reported worldwide including in Canada (Tewari et al., 2005;Rahman et al., 2014), China (Chai et al., 2014), and India (Bhattacharya et al., 2014). The clubroot pathogen survives as resting spores that can persist in the soil for many years (Dixon, 2009). Given the longevity of P. brassicae in infested soils and the significant economic value of Brassica crops, the management of clubroot has been a focus of agricultural researchers for decades. In recent years, clubroot has emerged as an important constraint to canola (Brassica napus; oilseed rape) production in western Canada, further increasing interest in this disease . While many strategies have been employed for clubroot control , the planting of clubroot resistance (CR) canola cultivars is the most effective and environmentally friendly approach to manage the disease (Rahman et al., 2014). The identification of effective resistance is the first step in breeding for this trait, with Brassica rapa (2n = 20, AA) considered a superior source of dominant major clubroot resistance genes than Brassica oleracea (2n = 18, CC) (Toxopeus et al., 1986;Hirai, 2006;Piao et al., 2009). Over the past 20 years, at least 15 CR genes have been identified in B. rapa, including CRa , CRb (Piao et al., 2004), CRb Kato  , CRk , Crr3 (Hirai et al., 2004;Saito et al., 2006), Rcr1 (Chu et al., 2014;Yu et al., 2016), CRc (Sakamoto et al., 2008;Matsumoto et al., 2012), Crr1, Crr2, Crr4 (Suwabe et al., 2003, CRd (Pang et al., 2018), BraA.Cr.a, BraA.Cr.b, BraA.Cr.c (Hirani et al., 2018), and CrrA5 (Nguyen et al., 2018). Clubroot resistance genes from B. rapa have been introgressed into several European B. napus oilseed cultivars, including 'Mendel' and 'Tosca' (Frauen, 1999).
In Canada, different Brasscia genotypes have been used as resistance donors in the breeding of CR canola/oilseed rape (B. napus). Fredua-Agyeman and  found that the B. napus cultivar 'Mendel' contained one dominant CR gene effective against pathotype 3H of P. brassicae, as defined on the Canadian Clubroot Differential Set . This was the dominant pathotype in Alberta, Canada (Strelkov et al., 2006(Strelkov et al., , 2007, at least prior to the introduction of CR canola cultivars beginning in 2009 . However, the planting of CR canola in short rotations over large acreages led to the rapid development of new pathotypes of P. brassicae , and 'Mendel' resistance was eroded and is no longer a good choice for breeding cultivars resistant to the new pathotypes. Rutabaga (Brassica napus ssp. napobrassica) is another potential donor of clubroot resistance genes. The rutabaga cultivars 'Brookfield' and 'Polycross' possessed excellent resistance to Canadian P. brassicae pathotypes 2, 3, 5, 6, and 8 Rahman, 2011, 2016;Hasan and Rahman, 2013). The ratio of resistant (R) to susceptible (S) in the F 2 generation derived from crossing both Brookfield and Polycross with susceptible B. napus lines was 3R:1S while segregation in the test cross family of the latter deviated from a 1R:1S ratio. This suggested that CR in Brookfield was controlled by a single dominant gene while resistance in Polycross was more complex (Hasan and Rahman, 2011).
With the recent emergence of new pathotypes of P. brassicae and their ability to overcome clubroot resistance in most CR canola cultivars, additional sources of resistance are needed. The CRa resistance gene was first detected in the European Clubroot Differential (ECD) 02 (B. rapa ssp. rapifera line AAbbCC) (Buczacki et al., 1975;Matsumoto et al., 1998). Zhang et al. (2016) reported one marker (i3e4) that was tightly linked to CRa. While ECD 02 appears to be resistant to all P. brassicae pathotypes identified in Canada (Strelkov et al., 2006Leboldus et al., 2012), and was used as a resistance source in studies from Japan Sakamoto et al., 2008), its application in clubroot resistance breeding in Canada has not yet been reported. The objectives of the current study were to introgress clubroot resistance from ECD 02 (male parent) into two susceptible female B. rapa genotypes CR 2599 and CR 1505, to evaluate the genetic basis of resistance to three isolates of P. brassicae representing different pathotypes and identify molecular markers associated with resistance to the new pathotypes.

Development of Mapping Populations
The parent B. rapa subsp. rapifera line AAbbCC (ECD 02) was resistant to all 17 P. brassicae pathotypes identified in Canada . In contrast, B. rapa accessions CR 2599 and CR 1505 ('Emma') were susceptible to these same pathotypes (Fredua-Agyeman et al., 2019) and served as the susceptible parents. ECD 02 is a winter-type while the two susceptible parents are spring-type.
To produce the F 1 plants, genetic crosses were carried out between June 2016 and January 2017 by emasculation followed by hand-pollination, with the plants kept in a growth chamber maintained under an 18 h photoperiod and temperatures of 21/18 • C (day/night). Vernalization of ECD 02 was carried out as described by Fredua-Agyeman et al. (2019). The susceptible parents were seeded much later to ensure that they flowered around the same time as ECD 02. Seeding, vernalization and the self-pollination of single F 1 plants to produce F 2 seeds were carried out in a growth chamber, cold room and greenhouses, respectively, at the Crop Diversification Centre North, Alberta Agriculture and Forestry, Edmonton, Alberta from March 2017 to June 2018. The ECD 02 × CR 2599 and ECD 02 × CR 1505 derived F 2 populations will be designated here as Popl#1 and Popl#2, respectively, for convenience.

Inoculum Preparation
Isolates of P. brassicae were stored as galled roots at −20 • C until needed. To prepare resting spore inoculum, the galls were homogenized in sterile distilled water in a Waring LB10G blender (Cole-Parmer) following . The resulting slurry was passed through eight layers of cheesecloth into a beaker to remove plant debris and other detritus, and the filtrate was collected. The resting spore concentration was estimated using a hemocytometer and adjusted to 1 × 10 8 spores/mL with sterile distilled water. Inoculum was kept at 4 • C and used within 24 h of preparation.

Clubroot Assays
Seedlings were geminated in Petri dishes (100 mm × 15 mm) on moistened Whatman no. 1 filter paper for 7 days at room temperature and a 12 h photoperiod. Inoculations were carried using the root-dip method (Nieuwhof and Wiering, 1961;Strelkov et al., 2006), with additional inoculum added by the pipette method (Lamers and Toxopeus, 1977 as cited by Voorrips and Visser, 1993). Briefly, the rootlets were dipped into the pathogen resting spore suspension for about 10-20 s and then planted in 8 × 4 flat trays filled with Sunshine Mix #4 Aggregate Plus (Sun Gro Horticulture, Seba Beach, AB, Canada) potting medium. This was followed by pipetting 1 ml of inoculum and dispensing into the soil surrounding each seedling.
The seedlings were transferred to a greenhouse maintained at 20-25/15-18 • C day/night with a 16 h photoperiod, and watered daily with slightly acidified water (pH ≈ 5.5-6.5, adjusted with HCL). Beginning at 3 weeks after inoculation, the plants were fertilized once a week with 20 N: 20 P: 20 K Classic Fertilizer with micronutrients (Plant Products Brampton, ON, Canada).
Eight weeks after inoculation, the plants were gently removed from the potting medium, washed in water, and assessed for clubroot severity on a 0-3 scale, where: 0 = no galls; 1 = a few small or bead-sized galls on <1/3 of the roots; 2 = medium galls on 1/3-2/3 of the roots, and 3 = large galls on >2/3 of the roots (Kuginuki et al., 1999;Xue et al., 2008). The susceptible B. napus 'Westar' was included as a positive control in all of the assays.

DNA Extraction
Three hundred sixty-eight leaves were collected from 46 F 2 individuals resistant (disease score = 0) and 46 F 2 individuals susceptible (disease score = 3) to each of pathotypes 5X and 5G of Popl#1 and Popl#2. Genomic DNA was extracted from the leaves of the parents and the 368 F 2 individuals using the cetyltrimethyl ammonium bromide (CTAB) method (Sambrook and Russell, 2001). The DNA concentration was quantified with a ND-2000c spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, United States) and the template DNA diluted to 20-25 ng/µL for PCR.
No molecular work was carried out on F 2 plants inoculated with pathotype 3H (for the two populations) because almost all previous genetic mapping studies in Canada have utilized this pathotype. Therefore, molecular analyses were completed on four of the six plant population/P. brassicae pathotype combinations.

PCR and SSR Genotyping
PCR amplification was carried out in a 12 µL reaction volume containing 2.5 µL of 5 × Taq buffer, 1.0 µL of 25 mM MgCl 2 , 0.25 µL of 10 mM dNTPs mix, 0.25 µL of 25 nM of each primer, 0.25 µL of 25 nM of fluorescently labeled M13 primer, 1.0 µL of 20-25 ng DNA template, and 1.25 U of Taq polymerase (Promega, Madison, WI, United States). The PCR cycling conditions consisted of an initial denaturation at 95 • C for 5 min, followed by 35 cycles at 95 • C for 1 min, 58 • C for 1 min and 72 • C for 1.5 min, and a final elongation step at 72 • C for 10 min. Aliquots of the PCR products were analyzed on an ABI PRISM 3730 × l DNA analyzer (Applied Biosystems). Amplified products from resistant and susceptible plants which differed by >200 bp were separated on 3% agarose gels. The gels were stained with SYBR Safe DNA Gel Stain (Thermo Fisher Scientific, Carlsbad, CA, United States) and visualized with a Typhoon FLA 9500 Variable Mode Laser Scanner Image Analyzer (GE Healthcare Life Sciences, Mississauga, ON, Canada).
Bulk segregant analysis as described by Michelmore et al. (1991) was carried using 144 PCR-based markers with DNA of the parents and two resistant and two susceptible DNA bulks from each of the two F 2 populations (Popl#1 and Pop#2). The markers included 13 markers spanning the A03 chromosome of B. rapa, 65 markers linked to six previously reported clubroot resistant genes: CRk, Crr3, CRb, CRb Kato , CRa, and Crr1, as well as 66 markers designed in this study. Of the 66 markers, 50 were SSR markers designed from the identified QTL regions while the remaining 16 were designed based on the CRa gene sequence of Ueno et al. (2012). The polymorphic markers were used to genotype the parents as well as the 384 F 2 individuals of Popl#1 and Pop#2, and a genotype matrix was constructed for pathotypes 5X and 5G for Popl#1 and Pop#2.

Map Construction and QTLs Analysis
Genetic linkage analysis was performed for each of the four combinations using MAPMAKER/EXP 3.0 (Lincoln et al., 1992). The logarithm of odds (LOD) score was set for a minimum of 3.0 and a recombination fraction (O) of 0.40. The Kosambi map function was used to convert the recombination frequencies into genetic distances (in centimorgans, cM) (Kosambi, 1944). Genetic linkage maps were constructed using MapChart v. 2.32 (Voorrips, 2002). Inclusive Composite Interval Mapping-Additive (ICIM-ADD) and Single Marker Analysis (SMA) methods were performed based on stepwise regression of simultaneous consideration of all markers using the IciMapping software V4.1 with 1 cM walking speed, 1000 permutations and P < 0.05 (Meng et al., 2015). To confirm a QTL, Composite Interval Mapping (CIM) and SMA were performed at the same mapping parameters as the above (Churchill and Doerge, 1994) using WinQTL Cartographer software v2.5 (Wang et al., 2011). A definitive QTL was declared if peaks with LOD score ≥ 3.0 were identified at the same position by both the IciMapping software V4.1 and WinQTL Cartographer software V2.5. Smaller peaks in close proximity to major peaks were regarded as artifacts. The QTLs identified in this study were named using a modified gene nomenclature system for the Brassica genus proposed by Østergaard and King (2008). The QTL name was in the order genus (1 letter), species (2 letters), genome (1 letter), chromosome number (1 letter), pathotype name (3 letters), closest published gene(s) (3-6 letters), and QTL number (2 letters). The additive effects of individual QTLs and epistatic interaction (QTL × QTL) were estimated as the percentage of phenotypic variation explained (PVE) by the IciMapping software V.4.1 using the ICIM-ADD and ICIM-EPI methods, respectively.

Statistical Analysis
To test for the inheritance of clubroot resistance in the two F 2 populations, the phenotypic data from different crosses of the same parents were subjected to Chi-square (χ 2 ) tests of homogeneity. In addition, individual data from the different crosses found to be uniform were pooled. The resulting data for each population were subjected to χ 2 goodness-of-fit for different segregation ratios using SAS v. 9.4 (SAS Institute, United States).
Eleven plants of ECD 02 showing absolute resistance (disease score = 0) to pathotype 5G were used as donor parents in crosses with the two susceptible parents CR 2599 and CR 1509. The F 1 plants produced from the two crosses were winter-types and required vernalization to induce flowering. Approximately 76% (25 of 33) and 25% (2 of 8) of the crosses carried out between ECD 02 and CR 2599 and CR 1509, respectively, resulted in siliques. Of the successful crosses, 24% (6 of 25) and 50% (1 of 2) produced 2-30 good quality rounded seeds per silique. All 27 of the F 1 plants screened with P. brassicae (2-9 pathotypes) were resistant to the pathogen (23 disease score = 0, 4 disease score = 1) (Supplementary Table S1). The 23 F 1 plants with a disease score = 0 were used as donors of clubroot resistance genes for the various crosses.
The development of two F 2 mapping populations from the crosses between ECD 02 and the two susceptible parents is summarized in Figure 1. One thousand five hundred and fortyeight and 710 F 2 individuals derived from 6 and 1 clubroot resistant F 1 families from Popl#1 and Popl#2, respectively, were evaluated for resistance to the P. brassicae pathotypes 3H (SACAN-ss1), 5X (LG-1), and 5G (CDCS). The frequency distribution of disease scores to the three pathotypes in the two F 2 populations is presented in Figure 2. All six distributions were bimodal, suggesting oligogenic or polygenic resistance to clubroot in the two F 2 populations.

Phenotypic Variation in F 2 Mapping Subpopulations
The χ 2 tests of homogeneity indicated that the phenotypic data from the 534 F 2 plants of Popl#1 inoculated with pathotype 3H and the 524 F 2 plants inoculated with pathotype 5X were significantly different (P < 0.00001) and hence could not FIGURE 1 | Development of F 2 mapping population for the study of the inheritance pattern of clubroot resistance introgressed from ECD 02.
be pooled (Supplementary Tables S2, S3). In contrast, no significant differences were found in the phenotypic data from the 490 F 2 plants of Popl#1 inoculated with pathotype 5G (Supplementary Table S4).
In the case of Popl#2, the 227, 244, and 239 F 2 plants inoculated with pathotypes 3H, 5X, and 5G, respectively, were produced by self-pollination of a single F 1 family from the cross ECD 02 × CR 1505 and hence, no χ 2 tests of homogeneity were conducted. Instead, the F 2 data for Popl#2 were tested directly for χ 2 goodness-of-fit for different segregation ratios.

Inheritance Pattern of Clubroot Resistance Derived From ECD 02
The segregation analysis was carried out in two ways. The first considered only those plants with a disease score = 0 as resistant (R), and all others as susceptible (S) (Tables 1, 2). The second considered plants with disease scores = 0 or 1 as R, and those with disease scores = 2 or 3 as S (Tables 1, 2). The Chi-square goodness of fit test showed that the segregation of clubroot resistance in the F 2 populations largely deviated from the expected Mendelian segregation ratios of 3:1, 15:1, and 63:1 for resistance controlled by a single, two or three dominant genes, respectively. Instead, deviations from the normal ratios, such as those observed in the event of non-allelic or linked interactions, were obtained.

Linkage Analyses and QTL Mapping
Of 144 PCR-based markers screened by bulk segregant analysis, 49 markers detected polymorphism between the parents of Popl#1 (i.e., ECD 02 and CR 2599), while 45 markers detected polymorphism between the parents of Popl#2 (i.e., ECD 02 and CR 1505). Twenty-seven of the 49 and 23 of the 45 markers detected polymorphisms in Popl#1 and Popl#2, respectively. Sixteen (14 from A03 + 2 from A08) of these markers detected polymorphism in the two F 2 populations, while 11 and 7 markers detected polymorphism only in Pop#1 and Popl#2, respectively. Table 3 provides the sequence information and origin of the 34 (16 + 11 + 7) polymorphic markers used to genotype the two F 2 populations.
At a LOD threshold of 3.0 and a recombination fraction (O) of 0.40, 14 of the 25 A03 chromosome markers used to screen the F 2 plants of Popl#1 inoculated with pathotype 5X were linked, while 11 of the A03 chromosome markers remained unlinked. Based on the ICIM-ADD method, two QTLs from the A03 chromosome of B. napus were detected for resistance to pathotype 5X. The SSR markers KB69N05 and B4732 flanked the first QTL BraA3P5X.CRa/b Kato 1.1 [LOD score = 3.6, located between 13.4 and 21.3 cM on Figure 3A and ≈24274312-24348056 base pairs (bp) on the physical map of B. napus]. The SSR markers BGA06 and KB29N19 bordered the second QTL BraA3P5X.CRa/b Kato 1.2 (LOD score = 15.9, located between 33.8-41.2 cM on Figure 3A and ≈24426905-24637310 bp on the physical map of B. napus). These two QTLs explained 4.5 and 51.0% of the phenotypic variance, respectively. Single Marker Analysis (SMA) by the IciMapping software showed that the SCAR marker GC2360-1 had the highest LOD and PVE (phenotypic variation explanation) scores of 14.2 and 48.9%, respectively, followed by marker KB59N06 with a LOD and PVE scores of 5.7 and 23.8%, respectively. Overall, the linked markers explained 55.4 and 72.7% of the phenotypic variance with the ICIM-ADD and SMA methods, respectively.
Similarly, in the F 2 plants from Popl#1 inoculated with pathotype 5G, 18 of the 25 markers on the A03 chromosome were linked, while three of the markers remained unlinked. The ICIM-ADD method identified two QTLs from the A03 chromosome of B. napus for resistance to pathotype 5G. The SSR markers KB59N06 and B4732 bordered the first QTL BraA3P5G.CRa/b Kato 1.1 (LOD score = 17.0, located between 32.7-46.6 cM on Figure 3B and ≈24262454-24348056 bp on the physical map of B. napus). The SSR markers CRaJY and BGB41 flanked the second QTL BraA3P5G.CRa/b Kato 1.2 (LOD score = 14.0, located between 76.8-92.0 cM on Figure 3B and ≈ 24557499-24579679 bp on the physical map of B. napus). These two QTLs explained 23.3 and 21.8% of the phenotypic variance, respectively. The SMA by the IciMapping software showed that SSR marker KB59N03 had the most significant association with LOD and PVE scores of 6.8 and 27.4%, respectively. This was followed by marker GC2360-1 with LOD and PVE scores of 6.3 and 25.7%, respectively. Overall, based on the IciMapping software and using the ICIM-ADD and SMA methods, the PVE by the linked markers was 45.2 and 53.1%, respectively.
For Popl#2 F 2 plants inoculated with pathotype 5X and 5G, the linkage analysis revealed that almost all (at least 19 of the 21) of the polymorphic markers on the A03 chromosome were linked. However, only the QTLs BraA3P5X.CRa/b Kato 1.2 and BraA3P5G.CRa/b Kato 1.2 were detected for resistance by the ICIM-ADD method in the case of the F 2 plants of Popl#2 inoculated with pathotype 5X (LOD score = 16.7, Figure 4A) and 5G (LOD score = 17.4, Figure 4B), respectively. The QTL was located in the interval between SSR markers CRaJY and KB29N19 in the genomic region located approximately 24557499-24637310 bp on the physical map of B. napus. The percent of phenotypic variance explained by the ICIM-ADD method was 35.6 and 32.5% for the F 2 plants inoculated with pathotype 5X and 5G, respectively. The SMA by the IciMapping software showed that SSR marker GC2360-1 had the most significant association in response to infection by P. brassicae pathotypes 5X and 5G. The LOD scores for the SMA was 6.6 and 4.7 for the F 2 plants inoculated with pathotype 5X and 5G while the PVE scores was 27.0 and 29.3%, for pathotype 5X and 5G, respectively.
The LOD profiles and PVE determined by the CIM and SMA methods implemented in the WinQTL Cartographer software for the F 2 plants from Popl#1 and Popl#2 inoculated with pathotypes 5X and 5G showed the same pattern as those determined with the IciMapping software (Supplementary Figures S1, S2). For example, the LOD score values for BraA3P5X.CRa/b Kato 1.1 and BraA3P5X.CRa/b Kato 1.2 for pathotype 5X were 4.2 and 8.2, respectively for the CIM method compared with 3.6 and 15.9 by the ICIM-ADD method. The LOD score values for BraA3P5G.CRa/b Kato 1.1 and BraA3P5G.CRa/b Kato 1.2 for pathotype 5G were 21.3 and 6.5, respectively for the CIM method compared with 17.0 and 14.0 by the ICIM-ADD method. In the case of Popl#2, the LOD scores for the QTLs BraA3P5X.CRa/b Kato 1.2 and BraA3P5G.CRa/b Kato 1.2 were 6.2 and 6.9 for pathotypes 5X and 5G, respectively, for the CIM method compared with 16.7 and 17.4 by the ICIM-ADD method. Thus, the LOD scores determined with the ICIM software were in general about two times higher than the values obtained with the WinQTL cartographer software. , which showed the most significant association with clubroot resistance in the F 2 plants of Popl#1 and Popl#2 inoculated with pathotypes 5X and 5G ranged from 11.5 to 18.8% and 12.5 to 22.9%, respectively. In the case of the two co-segregating markers on the A08 chromosome, recombination between the BRMS-088 allele and clubroot resistance in the F 2 plants of Popl#1 inoculated with pathotype 5X and 5G was 10.4 and 15.6%, respectively. Similarly, recombination between marker allele and phenotype was 7.3 and 25.0% in the F 2 plants of Popl#2 inoculated with pathotype 5X and 5G, respectively. Recombination between the SSR marker A08-5021 allele and clubroot resistance in the F 2 plants of Popl#1 inoculated with pathotypes 5X and 5G was 9.4 and 28.1%, respectively, while the recombination in Popl#2 was 7.3 and 11.5%, respectively. Thus, the recombination between the two A08 markers FIGURE 3 | QTL likelihood profile and partial linkage map of the A03 chromosome of Brassica napus and Brassica rapa showing resistance to P. brassicae pathotypes 5X (A) and 5G (B) in Popl#1. The LOD scores are indicated on the y-axis while the marker names and the genetic distances (in cM) are indicated on the x-axis. Two QTLs were detected for resistance to pathotypes 5X (A) and 5G (B). The first QTL (BraA3P5X.CRa/b Kato 1.1 and BraA3P5G.CRa/b Kato 1.1) was detected in the genomic region where the CRb Kato gene was mapped  while the second QTL (BraA3P5X.CRa/b Kato 1.2 and BraA3P5G.CRa/b Kato 1.2) was detected in the genomic region where the CRa gene was mapped Hayashida et al., 2008;Ueno et al., 2012). and clubroot resistance to pathotype 5X was much smaller (7.3-10.4%) compared with that for resistance to pathotype 5G (11.5-28.1%).
The additive effects detected by the IciMapping and the WinQTL cartographer software were positive for the QTLs, which suggested that the alleles conferring clubroot resistance were derived from the resistant parent ECD 02. Epistatic interaction (Q × Q) analysis showed PVE at levels ranging from 0.0 to 17.6% for the A03 × A03 QTLs, 35.6-53.4% for the A03 × A08 QTLs and 25.8-37.3% for the A08 × A08 QTLs. The results of the epistatic QTL analysis (Figure 5) suggested that genes from both the A03 and the A08 chromosomes are needed for resistance to be effective.

DISCUSSION
Clubroot is widespread in many of the Brassica growing areas of the world. Both quantitative and qualitative disease resistance have been reported in Brassica crops (Piao et al., 2009). The erosion of single dominant CR genes in commercial canola cultivars and the emergence of new virulent isolates of P. brassicae have been reported in Canada  and Europe (Oxley, 2007;Diederichsen et al., 2014;Wallenhamar et al., 2014;Zamani-Noor, 2017). The erosion of the effectiveness of CR genes has also occurred in cruciferous vegetables in Asia (Bhattacharya et al., 2014;Chai et al., 2014). The elevated infection in clubroot-resistant cultivars and volunteers would lead to increased spore load of the pathogen in the soil (Zamani-Noor and Rodemann, 2017). Hence, clubroot remains a huge problem and poses by far the most significant threat to cruciferous crop production worldwide.
One of the strategies to combat clubroot caused by the many pathotypes of P. brassicae is to deploy cultivars with multiple CR genes. The inheritance of different encoding genes provides a buffer when one resistance mechanism is overcome (Lagudah, 2011). However, the combined effects of inheriting gene combinations seems to be very complex. In durum wheat (Triticum turgidum L. ssp. durum), strong additive effects of the Lr34/Yr18 gene combination to leaf rust, stripe rust and powdery mildew were observed compared to the resistance effects of the Lr46/Yr29 gene combination to the three diseases (Lillemo et al., 2008). Therefore, in canola different gene combinations are expected to confer different levels of resistance to clubroot.
In this study, the inheritance of multiple CR genes was examined in F 2 plants derived from the crossing between B. rapa L. spp. rapifera line ECD 02 (turnip) (the resistant parent) and two B. rapa accessions (the susceptible parents). The fact that all the F 1 plants from the above crosses were highly resistant to pathotypes 3H, 5X, and 5G suggested that the resistance gene in ECD 02 is dominant and could be due to the combined effect of two or more genes. The F 1 means of all crosses were not significantly different from the resistant parent ECD 02 but were significantly different from the mean measurements of the two susceptible parents. These results indicate the complete dominance genes controlling clubroot resistance genes came from ECD 02. The segregation of F 2 plants for clubroot resistance (Tables 1, 2) and the fact that the F 2 means of all crosses were significantly different from the mean measurements of both parents suggested that the genetic variance consisted of both additive and dominance effects. The 3R:1S and 1R:3S ratios are consistent with the inheritance of a trait controlled by a single dominant major gene. Digenic dominant epistasis, represented as 12R:4S and 4R:12S, gives the same 3R:1S and 1R:3S ratios, respectively. Therefore, based on 3R:1S/1R:3S ratios, we could not conclude whether there was one or two CR genes in ECD 02. On the other hand, the 9R:7S, 7R:9S, 5R:11S, 11R:5S, 13R:3S, 3R:13S, and 1R:15S ratios are modifications of the 15R:1S segregation ratio expected for a trait controlled by two dominant major genes.
The distorted segregation ratios suggested that the resistance genes were on different chromosomes (Hayman and Mather, 1955). The 9R:7S ratio confirms the existence of two genes with duplicate recessive epistasis, which suggested that the dominant allele at the two loci were necessary to control the clubroot caused by pathotypes 5X and 5G. In other words, individuals with double recessive at either locus or both loci were susceptible. The 13R:3S ratio shows recessive suppression in which the presence of two recessive alleles of one gene are needed to suppress the effects of the dominant gene. The 5R:11S, 11R:5S, 7R:9S, 3R:13S, and 1R:15S ratios indicated two-gene control of clubroot resistance with digenic additive epistasis or quantitative control of the resistance. In the current study, most of the segregation ratios fit the two-gene model although segregation patterns of 3R:1S, 1R:3S, and 1R:63S were also possible. However, the mapping population was of the F 2 generation and hence it is not clear if the preponderance of non-allelic gene effects will alter at advanced generations. Pioneering studies based on segregation ratios suggested that at least two dominant major genes (originally designated A and C) controlled clubroot resistance in ECD 02 (Wit and Van De, 1964;Buczacki et al., 1975). Therefore, the results of the greenhouse screening work was in agreement with the published literature.
In spite of the substantial contribution of non-additive effects to the variation of complex traits, gene effects controlling clubroot resistance have not been studied. Yu et al. (2017) reported that the Rcr8 gene on the A02 chromosome and Rcr9 gene on the A08 chromosome conferred resistance to pathotype 5X. However, the induced resistance to pathotype 5X was not correlated with the resistance to pathotypes 2F, 3H, 5I, 6M, and 8N conferred by the Rcr4 gene on the A02 chromosome (Yu et al., 2017). On the other hand, our study shows that resistance conferred by the CRa/CRb Kato gene(s) on the A03 and the Crr1 gene on the A08 chromosome also conferred resistance to pathotype 3H. In addition, Yu et al. (2017) did not show the different interactions involving the multiple (Rcr4, Rcr8, and Rcr9) resistance genes they identified. To the best of our knowledge, our study is the first report that demonstrates that the two CR genes in ECD 02 interact in a non-additive manner to confer resistance to clubroot. Such non-allelic interactions of multiple genes have been reported in the response of many plants to fungi, bacteria, virus and insect attack. For example, barrel clover (Medicago truncatula) to aphid attack (Kamphius et al., 2020); mungbean (Vigna radiata) to mungbean yellow mosaic virus attack (Akbar et al., 2017); soybean (Glycine max) to rust infection (Pierozzi et al., 2008) and cotton (Gossypium hirsutum) to insect pest and virus attacks (Ahuja et al., 2007).
In summary, clubroot tests, linkage analysis and QTL mapping carried out in this study demonstrated that the CRa/CRb Kato and the Crr1 genes on the A03 and the A08 chromosomes of ECD 02 interact in a non-allelic manner to confer resistance to pathotypes 5X and 5G. Based on the QTL analysis, the genetic control against virulent P. brassicae pathotypes may also involve additional genes modulating the action of the two major genes. The presence of at least the two dominant genes complementing each other might explain why ECD 02 confers strong and highly stable qualitative resistance to many P. brassicae pathotypes from around the world. Knowledge of gene effects controlling clubroot resistance offers the possibility of exploiting ECD 02 resistance for the breeding of clubroot resistant canola cultivars and cruciferous vegetables. In addition, the genomic regions identified in this study will provide additional resources for marker-assisted selection in Brassica breeding programs.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.