RNA-Seq Bulked Segregant Analysis of an Exotic B. napus ssp. napobrassica (Rutabaga) F2 Population Reveals Novel QTLs for Breeding Clubroot-Resistant Canola

In this study, a rutabaga (Brassica napus ssp. napobrassica) donor parent FGRA106, which exhibited broad-spectrum resistance to 17 isolates representing 16 pathotypes of Plasmodiophora brassicae, was used in genetic crosses with the susceptible spring-type canola (B. napus ssp. napus) accession FG769. The F2 plants derived from a clubroot-resistant F1 plant were screened against three P. brassicae isolates representing pathotypes 3A, 3D, and 3H. Chi-square (χ2) goodness-of-fit tests indicated that the F2 plants inherited two major clubroot resistance genes from the CR donor FGRA106. The total RNA from plants resistant (R) and susceptible (S) to each pathotype were pooled and subjected to bulked segregant RNA-sequencing (BSR-Seq). The analysis of gene expression profiles identified 431, 67, and 98 differentially expressed genes (DEGs) between the R and S bulks. The variant calling method indicated a total of 12 (7 major + 5 minor) QTLs across seven chromosomes. The seven major QTLs included: BnaA5P3A.CRX1.1, BnaC1P3H.CRX1.2, and BnaC7P3A.CRX1.1 on chromosomes A05, C01, and C07, respectively; and BnaA8P3D.CRX1.1, BnaA8P3D.RCr91.2/BnaA8P3H.RCr91.2, BnaA8P3H.Crr11.3/BnaA8P3D.Crr11.3, and BnaA8P3D.qBrCR381.4 on chromosome A08. A total of 16 of the DEGs were located in the major QTL regions, 13 of which were on chromosome C07. The molecular data suggested that clubroot resistance in FGRA106 may be controlled by major and minor genes on both the A and C genomes, which are deployed in different combinations to confer resistance to the different isolates. This study provides valuable germplasm for the breeding of clubroot-resistant B. napus cultivars in Western Canada.


Introduction
Clubroot, caused by Plasmodiophora brassicae Woronin, is an important soilborne disease of cruciferous crops worldwide [1,2].Infection by P. brassicae results in excessive growth and division of the host root cells, resulting in the formation of root galls and an eventual reduction in the plant's capacity for water and nutrient uptake [3,4].The cruciferous genus Brassica is known for its economically important agricultural and horticultural crops [5].These include Chinese cabbage, turnip, Polish canola, and other crops belonging to the species Brassica rapa (A genome); cabbage, cauliflower, broccoli, kale, Brussels sprouts, and others classified as B. oleracea (C genome), and rutabaga and canola/oilseed rape, which are B. napus (AC genome) [1,5,6].Globally, average yield losses caused by clubroot are estimated at 10% to 15% but may be as high as to 30% to 100% under favourable conditions [1,7,8].The clubroot pathogen survives as resting spores that can persist in the soil for many years, making the management of this disease difficult [9,10].
In Alberta and other Canadian provinces, clubroot has emerged as a constraint to canola (B.napus var.napus L.) production [11][12][13].The number of P. brassicae-infested fields in Alberta has increased from 12 in 2003 [14] to 3894 individual fields by 2022 [15].Although clubroot-resistant canola varieties represent the most effective and environmentally friendly strategy for clubroot management [13,16], P. brassicae populations show high diversity in terms of virulence and can quickly adapt to overcome host resistance [11,12,17].Over the past decade, 'resistance-breaking' pathotypes have been documented in hundreds of fields across Alberta [11,12,18].Forty-three pathotypes of P. brassicae, as classified on the Canadian Clubroot Differential (CCD) set [11], have been reported to date from Canadian collections of the pathogen [12].A majority of these pathotypes are highly virulent on canola cultivars carrying 'first-generation'-type resistance [12], which appears to be derived from the European oilseed rape cv.'Mendel' [19].
Genetic mapping is important for the identification of clubroot resistance (CR) gene loci and for the development of molecular markers for marker-assisted selection (MAS).Conventional PCR-based markers, such as amplified fragment length polymorphisms (AFLPs), cleaved amplified polymorphic sequences (CAPSs), random amplification of polymorphic DNA (RAPD), restriction fragment length polymorphisms (RFLPs), sequence characterized amplified regions (SCARs), sequence tagged sites (STSs), and simple sequence repeats (SSRs), were used widely for linkage-based identification and mapping clubrootresistance gene loci before the era of sequencing technologies [20][21][22][23][24][25][26].Next-generation sequencing (NGS) for genetic mapping has facilitated the development and application of genomics tools in plant breeding, such as genotyping-by-sequencing, single-nucleotide polymorphism (SNP) arrays, and bulked segregant analysis (BSA) [27].Bulked segregant RNA sequencing (BSR-seq) is one of the most cost-effective methods for mapping genes of interest via BSA.This process evaluates two bulk DNAs or RNAs of plants with different phenotypes for differentially expressed genes (DEGs) and maps QTLs by variant calling [28].Multiple CR gene loci against Canadian pathotypes of P. brassicae, such as Rcr1, Rcr2, Rcr3, Rcr6 and Rcr9 wa , have been identified via BSA or BSR-seq [29][30][31][32].
Rutabaga (B.napus ssp.napobrassica) could be a good source of clubroot resistance genes for emerging virulent P. brassicae pathotypes.Old rutabaga varieties like 'Wilhelmsburger' are known to carry resistance effective against most Canadian P. brassicae pathotypes [11,33].Hasan and Rahman [34] found a Canadian rutabaga cv.'Brookfield', which was resistant to all five 'old' pathotypes (2F, 3H, 5I, 6M, and 8N) found in Canada prior to the introduction of clubroot-resistant canola.Wang et al. [35] used a rutabaga cv.'Polycross' as a resistance donor to breed for canola populations resistant to three Canadian pathotypes.Fredua-Agyeman et al. [21] observed that 87.9% of 124 rutabaga accessions from Nordic countries showed resistance to at least one of 16 Canadian P. brassicae pathotypes.
In this study, an F 2 population derived from rutabaga accession FGRA106 that was reported to be resistant to 17 isolates representing 16 pathotypes of P. brassicae [21] was evaluated for its reaction to pathotypes 3A, 3D, and 3H.The inheritance of the resistance was determined based on segregation ratios.The genomic regions that co-segregated with resistance were determined based on the number of differentially expressed genes (DEGs), and the quantitative trait loci (QTLs) were mapped by BSR-seq.

Clubroot Tests
The chi-square tests of homogeneity indicated that the phenotypic data of F 2 plants inoculated with all three pathotypes were not significantly different in the three replicates (Table S1).Therefore, data for the same pathotype were pooled for analysis.The frequency distribution of disease ratings to the three pathotypes is presented in Figure 1.The inoculation conducted on the F 2 plants with P. brassicae pathotype 3A showed that 12.3%, 7.4%, 32.1%, and 48.3% (n = 408) of the plants were rated as 0, 1, 2, or 3, respectively (Table S1).The screening results with pathotype 3D indicated that 29.4%, 11.4%, 17.8%, and 41.4% of the F 2 plants (n = 411) exhibited disease ratings of 0, 1, 2, or 3, respectively (Table S1).In response to P. brassicae pathotype 3H, 23.0%, 9.8%, 13.0%, and 54.2% of the F 2 population showed disease ratings of 0, 1, 2, or 3 (n = 439) (Table S1).Based on Fisher's LSD test on the disease reaction data, the virulence of the pathotypes on the F 2 population was in the order of 3A > 3H = 3D.

Figure 1.
Frequency distribution of clubroot disease severity ratings in an F2 population derived from FGRA106 (♀) × FG769 (♂) to isolates representing pathotypes 3A, 3D, and 3H of Plasmodiophora brassicae.Plants were grown under greenhouse conditions and evaluated for clubroot severity on a 0-3 disease severity scale as described by Kuginuki et al. [36] and Strelkov et al. [3] at 7 weeks following inoculation with each pathotype.

Inheritance of Clubroot Resistance in F2 Populations
The chi-square goodness-of-fit test was carried out in two ways following the protocol of Fredua-Agyeman et al. [37].The first method grouped plants with a disease rating = 0 or 1 as resistant (R) and those with scores 2 or 3 as susceptible (S); meanwhile, the second method grouped plants with a disease rating = 0 as R and all others as S. The results obtained with the first method showed that the segregation of clubroot resistance in the F2 population was not significantly (p < 0.05) different from the expected Mendelian segregation ratios (R:S) of 3:13, 7:9, and 5:11 for pathotypes 3A, 3D, and 3H, respectively, all of which fit the two-gene models (Table 1).The second method indicated that the R:S ratios for pathotypes 3D and 3H were not significantly different from 5:11 (two-gene model) or 1:3 (one-gene model), respectively, while the ratio for pathotype 3A significantly deviated from all assumptions (Table 1).

Inheritance of Clubroot Resistance in F 2 Populations
The chi-square goodness-of-fit test was carried out in two ways following the protocol of Fredua-Agyeman et al. [37].The first method grouped plants with a disease rating = 0 or 1 as resistant (R) and those with scores 2 or 3 as susceptible (S); meanwhile, the second method grouped plants with a disease rating = 0 as R and all others as S. The results obtained with the first method showed that the segregation of clubroot resistance in the F 2 population was not significantly (p < 0.05) different from the expected Mendelian segregation ratios (R:S) of 3:13, 7:9, and 5:11 for pathotypes 3A, 3D, and 3H, respectively, all of which fit the two-gene models (Table 1).The second method indicated that the R:S ratios for pathotypes 3D and 3H were not significantly different from 5:11 (two-gene model) or 1:3 (one-gene model), respectively, while the ratio for pathotype 3A significantly deviated from all assumptions (Table 1). 1 Plants with clubroot disease severity ratings of 0 and 1 were regarded as resistant (R), and those with ratings of 2 and 3 as susceptible (S). 2 Plants with a clubroot disease severing rating of 0 were regarded as R, and those with ratings of 1, 2, and 3 as S.

RNA Sequencing, Filtering, and Sequence Alignment
The raw RNA sequences of the resistant and susceptible bulks were filtered, and the adapters were removed.Subsequently, the number of clean reads retained in the resistant bulks ranged from 24.2 to 28.7 Gb, 23.5 to 26.5 Gb, and 22.3 to 25.0 Gb, while in the susceptible bulks, this number ranged from 21.5 to 42.9 Gb, 20.7 to 23.2 Gb, and 21.8 to 23.8 Gb for pathotypes 3A, 3D, and 3H, respectively (Table S2).Therefore, the RNA sequencing data yielded 20× to 30× the genome size of B. napus.The GC content ranged from 47% to 48%.Approximately 89.7% to 93.1% of these reads were mapped to the B. napus cv.'ZS11' reference genome v2.0, 68.6% to 71.9% of which mapped only to exonic gene regions (Table S2).The mismatch rate per base ranged from 0.9% to 1.1% (Table S2).Therefore, the sequencing data were of adequate quality for the subsequent analysis.

SNP Calling and Marker Distribution
A total of 338,177, 331,344, and 325,623 SNP markers were obtained for the comparisons between the reference (B.napus cv.'ZS11') genome and the resistant and susceptible bulks from the inoculation experiments with pathotypes 3A, 3D, and 3H, respectively (Table 2 and Figure 2).The SNP marker densities on all 19 B. napus chromosomes are presented in Table 2 and Figure 2. The SNP densities for the pathotype 3A, 3D, and 3H bulks on the different chromosomes ranged from 190.78 to 732.31, 187.07 to 719.08, and 189.14 to 692.25 SNPs/Mb, respectively (Table 2).The highest SNP density was found on chromosome A10 for the three pathotypes, while the lowest occurred on chromosome C02 (Table 2).Consistently high SNP densities were observed at the beginning of chromosomes A01, A02, and A03, as well as at the end of chromosome A10, across all three pathotypes.Conversely, a region of low coverage density was noted on chromosome C09 (Figure 2).

Identification of QTLs Associated with Clubroot Resistance
Based on ∆(SNP-index) statistics at a 99% confidence interval (CI) from the variant calling between the R and S bulks of the three pathotypes, a total of 12 QTLs associated with resistance to 3 three pathotypes were detected on 7 of the 19 chromosomes of B. napus (Figure 6, Table S6).These QTLs were located on chromosomes A01 (1), A05 (1), A08 (4), C01 (3), C07 (2), C08 (1), and C09 (1) (Table S6).The peak |∆(SNP-index)| values ranged from ~0.23 to 0.53 (Table S6).The higher the |∆(SNP-index)|, the stronger the correlation between marker SNPs and traits.In this study, the QTLs were classified as major if the |∆(SNP-index)| > 0.32 and the peak was clearly above the 99% confidence interval; conversely, QTLs were classified as minor if the |∆(SNP-index)| < 0.30 and the peak fell between the 95% and 99% confidence intervals (Table S6).from ~0.23 to 0.53 (Table S6).The higher the |Δ(SNP-index)|, the stronger the correlat between marker SNPs and traits.In this study, the QTLs were classified as major if |Δ(SNP-index)| > 0.32 and the peak was clearly above the 99% confidence interv conversely, QTLs were classified as minor if the |Δ(SNP-index)| < 0.30 and the peak between the 95% and 99% confidence intervals (Table S6).The QTLs were named following the Brassica gene nomenclature system propos by Østergaard and King [39], as modified by Fredua-Agyeman et al. [37].For examp the QTL on chromosome A01 was designated BnaA1P3D.CRX1.1,where the first lett denotes the genus (Brassica), the second and third letters the species (napus), the fou The QTLs were named following the Brassica gene nomenclature system proposed by Østergaard and King [39], as modified by Fredua-Agyeman et al. [37].For example, the QTL on chromosome A01 was designated BnaA1P3D.CRX1.1,where the first letter denotes the genus (Brassica), the second and third letters the species (napus), the fourth letter the genome (A), the fifth letter the chromosome (1), and the sixth, seventh, and eighth letters (P3D) the pathotype of P. brassicae used for inoculation.These are followed by the name(s) of the closest published CR gene(s) (3-8 letters) or the letter X if no previous markers have been reported (CRX), and finally, the number of the QTL number (two digits, 1.1).

Discussion
The German rutabaga cv.'Wilhelmsburger' (ECD 10) was originally proposed as a differential host by Williams [33] and subsequently included in both the European Clubroot Differential (ECD; [44]) and CCD [11] sets.However, Yu et al. [22] and Fredua-Agyeman et al. [21] observed different resistance phenotypes in seven 'Wilhemsburger' accessions from Denmark, FGRA106, FGRA107, FGRA108, FGRA109, FGRA110, FGRA111, and FGRA112, when these were challenged with the same set of isolates representing 16 different P. brassicae pathotypes from Canada.Only the 'Wilhemsburger' accession FGRA106 from Denmark showed broad-spectrum clubroot resistance comparable with that of 'Wilhelmsburger' (ECD 10) from Germany; both were resistant to the 16 pathotypes tested by Fredua-Agyeman et al. [21].Given the proximity of Denmark to Germany as neighboring countries and the potential for germplasm movement, it is plausible that the 'Wilhelmsburger' accession FGRA106, based on its reactions, could be equivalent to 'ECD 10' and thus might harbor the same CR gene(s).
'Wilhelmsburger' has been used as a clubroot resistance donor in breeding programs worldwide for many decades [45].Lammerink [46] evaluated F2 progenies derived from 'Wilhelmsburger' with a P. brassicae isolate designated 'Race B' and suggested that the

Discussion
The German rutabaga cv.'Wilhelmsburger' (ECD 10) was originally proposed as a differential host by Williams [33] and subsequently included in both the European Clubroot Differential (ECD; [44]) and CCD [11] sets.However, Yu et al. [22] and Fredua-Agyeman et al. [21] observed different resistance phenotypes in seven 'Wilhemsburger' accessions from Denmark, FGRA106, FGRA107, FGRA108, FGRA109, FGRA110, FGRA111, and FGRA112, when these were challenged with the same set of isolates representing 16 different P. brassicae pathotypes from Canada.Only the 'Wilhemsburger' accession FGRA106 from Denmark showed broad-spectrum clubroot resistance comparable with that of 'Wilhelmsburger' (ECD 10) from Germany; both were resistant to the 16 pathotypes tested by Fredua-Agyeman et al. [21].Given the proximity of Denmark to Germany as neighboring countries and the potential for germplasm movement, it is plausible that the 'Wilhelmsburger' accession FGRA106, based on its reactions, could be equivalent to 'ECD 10' and thus might harbor the same CR gene(s).
'Wilhelmsburger' has been used as a clubroot resistance donor in breeding programs worldwide for many decades [45].Lammerink [46] evaluated F 2 progenies derived from 'Wilhelmsburger' with a P. brassicae isolate designated 'Race B' and suggested that the resis-tance was controlled by one dominant gene based on a 3R:1S segregation ratio.Similarly, Ayers and Lelacheur [47] reported that, based on segregation ratios of an F 2 population derived from 'Wilhelmsburger', the resistance to P. brassicae race 2 (sensu Williams, 1966 [33]) was controlled by two dominant genes, whereas resistance to race 3 was controlled by one dominant gene.In a separate study involving an F 2 population, Gustafsson and Falt [48] reported that the resistance of 'Wilhelmsburger' to a less virulent isolate 'Pb3' may have been conferred by two genes, whereas resistance to a highly virulent isolate 'Pb7' appeared to involve only one gene.In contrast, Crute et al. [49] suggested that 'Wilhelmsburger' possesses three clubroot resistance genes.These observations indicate that 'Wilhelmsburger' may carry multiple CR genes that could be differentially effective depending on the virulence of specific P. brassicae isolates.In the current study, the segregation ratios of F 2 plants inoculated with Canadian P. brassicae isolates representing various pathotypes were analyzed.The results suggested that the resistance inherited from FGRA106 to pathotypes 3A and 3D was likely determined by two genes, whereas resistance to pathotype 3H was conferred by either one or two genes.
Two major effect QTLs on the C genome were mapped to chromosomes C01 and C07.The first major QTL, BnaC1P3H.CRX1.2, which conferred resistance to pathotype 3H, was distant from the QTL Rcr_C01-1 reported on chromosome C01 of B. oleracea [58] (Figure 10).Two minor effect QTLs, BnaC1P3H.CRX1.1 and BnaC1P3H.CRX1.3 on chromosome C01, which also conferred resistance to pathotype 3H, were identified as proximal and distal, respectively, to a major effect QTL BnaC1P3H.CRX1.2.The second major effect QTL, BnaC7P3A.CRX1.1, which conferred resistance to pathotype 3A, was located on the bottom half of chromosome C07.A minor QTL, BnaC7P3D.CRX1.1, was also located on the bottom half of chromosome C07.These results confirm the bottom half of chromosome C07 as a genomic hotspot for several clubroot resistance genes, including Rcr7, qCRc7-1, qCRc7-2, qCRc7-3, and qCRc7-4 [52,53].Therefore, the clubroot resistance derived from the donor FGRA106 in this study appears to be conferred not only by major QTLs on the A genome, but also by major and minor QTLs located on the C genome.
J. Mol.Sci.2024, 25, x FOR PEER REVIEW 20 o C07 as a genomic hotspot for several clubroot resistance genes, including Rcr7, qCRc7 qCRc7-2, qCRc7-3, and qCRc7-4 [52,53].Therefore, the clubroot resistance derived from donor FGRA106 in this study appears to be conferred not only by major QTLs on th genome, but also by major and minor QTLs located on the C genome.The plant immune system is generally considered to comprise two layers: pathog associated molecular patterns (PAMPs)-triggered immunity (PTI), which offers ba protection against many pathogens, and effector-triggered immunity (ETI), which resu The plant immune system is generally considered to comprise two layers: pathogenassociated molecular patterns (PAMPs)-triggered immunity (PTI), which offers basal protection against many pathogens, and effector-triggered immunity (ETI), which results in robust and localized responses against specific pathogens [59].The GO and KEGG enrichment analyses conducted in this study suggested that the DEGs between the R and S bulks were associated with both of these layers of defense.The GO analysis revealed the differential expression of genes involved in bioprocesses related to sucrose, actin filaments, and sulfur amino acids, while enriched KEGG pathways were associated with metabolic pathways, ribosomes, and tRNA.These processes and pathways have been implicated as common initial defense signaling processes in eukaryotes [60][61][62][63][64]. Additionally, pathways involving inositol phosphate, glycine, serine, and threonine were also identified and have been implicated in the recognition of PAMPs by pattern recognition receptors (PRRs) in the host and in subsequent defense responses [60,[65][66][67][68]. Several GO processes and KEGG pathways were also identified that have been implicated in ETI, including GO bioprocesses related to methionine, phospholipid and chitin and KEGG pathways associated with glutathione, phenylpropanoid, and flavonoids [42,62,[69][70][71][72].
Some of the genes identified in the QTL regions through the SNPs and differential gene expression analyses have also been associated with PTI and ETI.For instance, heat shock proteins (HSPs) function as chaperones, playing roles in protein folding, assembly, translocation, and degradation during both abiotic and biotic stress.These processes are vital for the formation of PRRs and intracellular responsive proteins essential for resistance [73].Polyubiquitin modulates cellular protein turnover and homeostasis in basal host defense to abiotic and biotic stresses, and it is involved in the responsive modification of proteins in both PTI and ETI [74].The ethylene-responsive transcription factor ERF109 plays a key role in ethylene-mediated defense pathways during ETI [75].GLABROUS1 enhancer-binding protein-like, UDP-glucosyltransferase, and 60S ribosome proteins are implicated in ETI, and silencing of these genes can activate plant defense pathways [76][77][78][79].Additionally, downy mildew resistance 6 serves as a resistance gene in ETI [80].In a recent study of differential gene expression in 'Wilhelmsburger' in response to inoculation with pathotype 3A, Zhou et al. [81] found that salicylic acid and ethylene-mediated defense were involved in the host reaction.

Phenotyping Assays
The parents and F 2 population were screened against one single-spore isolate each of P. brassicae pathotypes 3A and 3H and a field isolate of pathotype 3D, as classified based on the CCD set [11].A total of 1620 F 2 individuals were tested in three rounds of bioassays, in which 180 plants were inoculated with each pathotype in each experiment.The inoculations were conducted as described previously by Strelkov et al. [3,83] with slight modifications.Briefly, clubbed roots were blended in sterile water and filtered through two layers of cheesecloth to generate a resting spore suspension.The spore concentration was then adjusted to 1 × 10 7 resting spores/mL with sterile water.For inoculation, the roots of 7-day-old seedlings were dipped into the spore suspension for about 10 s and then planted in plastic pots (6 cm × 6 cm × 6 cm) filled with Sunshine Mix #4 potting mixture (Sun Gro Horticulture, Seba Beach, AB, Canada) at a density of one seedling per pot.One millilitre of the resting spore suspension was then pipetted in the potting mix around each seedling to ensure infection and minimize disease escape.In addition to the resistant and susceptible parents FGRA106 and FG769, respectively, the susceptible B. napus cv.'Westar' was also included as a positive control in all experiments.
The inoculated plants were kept in a greenhouse maintained at 25 • C/18 • C day/night with a 16 h light period (natural light supplemented with artificial lighting) and assessed for clubroot symptoms at 7 weeks after inoculation (wai) on a 0-3 disease severity scale as described by Kuginuki et al. [36] and Strelkov et al. [3], where: 0 = no galling, 1 = slight galling on side roots, 2 = moderate galling on main and side roots, and 3 = severe galling with almost no observable side roots.

Bulk Construction and RNA Extraction
RNA extraction from the R (resistant) and S (susceptible) plant pools was based on the phenotypic reactions of individual plants to the three pathotypes.For each pathotype, 15 plants with a disease rating of 0 were pooled into an R bulk, while each S bulk consisted of 15 plants with a rating of 2 or 3. Three biological replicates of both the R and S bulks were assigned for each pathotype.
Leaf samples of each bulk collected at 7 weeks after inoculation were mixed and ground into powder in liquid nitrogen.The total RNA from each bulk replicate was extracted from 0.1 mL (~100 mg) of powdered root tissue of each sample using an RNeasy Plant Mini Kit (Qiagen; Toronto, ON, Canada) and purified using an RNase-Free Dnase kit (Qiagen; Toronto, ON, Canada).The RNA concentration was measured with a NanoDrop 2000c spectrophotometer (Thermo Scientific; Waltham, MA, USA), and its quality (RNA integrity numbers (RIN) ≥ 6.5 and a 28S/18S ratio ≥ 1.0) was confirmed using an Agilent 2200 TapeStation system (Agilent, Santa Clara, CA, USA).

RNA Sequencing
The cDNA library preparation and RNA sequencing were performed by the Oklahoma Medical Research Foundation NGS Core (Oklahoma City, OK, USA) with an IDT xGen RNA Library kit (Integrated DNA Technologies, Inc., San Diego, CA, USA) and an Illumina NovaSeq 6000 S4 platform (Illumina; San Diego, CA, USA).Pair-end read sequences (2 × 150 bp) were generated in 'fastq' format for further analysis.
The variant calling was performed using the GATK v4.2.2.0 function 'HaplotypeCaller', with the SNPs detected being subsequently filtered using the GATK 'VariantFiltration' function under proper standards ('QD < 2.0||MQ < 40.0||FS > 60.0||SOR > 3.0').The final .vcffiles were converted to .tableformat using the GATK 'VariantsToTable' tool for analysis in R 4.2.1 [84,91].SNP frequency calling was carried out on resistant and susceptible bulks as described by Liu et al. [28], Yu et al. [29], and Wu et al. [92] with modifications.Comparisons of SNP variants between the bulks were performed with the R package 'QTLseqr' [93] with the following filter settings: refAlleleFreq = 0.20, minTotalDepth = 50, maxTotalDepth = 500, minSampleDepth = 80, minGQ = 99.The detection of QTLs was based on the SNP-index and the ∆(SNP-index) in a 1 Mb sliding window [38].The SNPindex statistic calculates marker association differences in the genotype frequencies of mixed pools, where a value of 0 indicates that the short reads contain genomic fragments from the reference parent, while a value of 1 indicates that all of the short reads represent the genome from the other parent [38].A ∆(SNP-index) graph was used to detect the differences between the 'highest' and 'lowest' pools of extreme phenotypes, where ∆(SNP-index) = 1 and −1 indicates bulk DNA from one parent and the other parent, respectively, while ∆(SNP-index) = 0 if both parents have the same SNP-indices at the genomic regions [38].The physical positions of the detected QTLs were visualized with MapChart v2.3.2 [94] and compared with previously reported QTLs for clubroot resistance.

Statistical Analyses
The phenotypic data from different replicates of the same inoculum (pathotype) were subjected to chi-square tests of homogeneity.A chi-square goodness-of-fit test was performed to determine the segregation of the phenotypic data for all three pathotypes using R 4.2.1.Other built-in R packages including ggplot2, reshape2, and ggrepel were also used for data analysis or visualization.

Conclusions
The CR donor FGRA106 and the resistant F 2 progeny evaluated in this study were confirmed to carry resistance to P. brassicae pathotypes 3A, 3D, and 3H, which are predominant in canola in Western Canada [12].The resistance donor FGRA106 exhibited reactions similar to ECD10 and was previously reported to be resistant or moderately resistant to 17 isolates representing 16 pathotypes of P. brassicae [21].Based on the DEGs, QTLs, and associated GO terms and KEGG pathways, gene loci conferring resistance to pathotype 3A were mapped to chromosomes A05 and C07, while major QTLs for resistance to pathotypes 3D and 3H co-segregated to at least three genomic regions on chromosome A08.Another major QTL on chromosome C01 was required for resistance to pathotype 3H.The CR donor and the SNP markers identified in this study may serve as valuable resources for clubroot resistance breeding in Canadian canola.

Figure 4 .
Figure 4. Volcano plots from the log2 fold change (x−axis) and −log10 P adj (y−axis) values for the selection of 428, 67, and 98 differentially expressed genes (DEGs) in resistant (R) or susceptible (S)

Figure 7 .
Figure 7.The single-nucleotide polymorphisms (SNPs) and QTLs identified from clubroot-resistant and susceptible bulks of F2 plants inoculated with Plasmodiophora brassicae pathotypes 3A, 3D, and 3H on chromosomes A01, A05, A08, C01, C07, C08, and C09.The positions and names of SNPs are denoted on each chromosome (left and right, respectively), while the flanking regions of identified QTLs are represented as red bars, both on and to the right of each chromosome.

Figure 7 .
Figure 7.The single-nucleotide polymorphisms (SNPs) and QTLs identified from clubroot-resistant and susceptible bulks of F 2 plants inoculated with Plasmodiophora brassicae pathotypes 3A, 3D, and 3H on chromosomes A01, A05, A08, C01, C07, C08, and C09.The positions and names of SNPs are denoted on each chromosome (left and right, respectively), while the flanking regions of identified QTLs are represented as red bars, both on and to the right of each chromosome.

Figure 9 .
Figure 9. KEGG pathways associated with differentially expressed genes between clubroot-resistant and susceptible bulks of F 2 plants derived from FGRA106 (♀) × FG769 (♂) and tested with Plasmodiophora brassicae pathotypes 3A (a), 3D (b), and 3H (c).Visualization of pathways was sorted by fold enrichment.The colors indicate a −log10 false discovery rate (FDR) as per the scale on the right-hand side.

Figure 10 .
Figure 10.Comparison between previously reported clubroot-resistant (CR) gene loci and Q identified in this study on chromosomes A05, A08, C01, and C07.The position and name of SN identified in this study are indicated on each chromosome (left and right, respectively), and flanking regions of previously reported CR gene loci and QTLs identified in this study are deno by black and red bars, respectively, both on and to the right of each chromosome.

Figure 10 .
Figure 10.Comparison between previously reported clubroot-resistant (CR) gene loci and QTLs identified in this study on chromosomes A05, A08, C01, and C07.The position and name of SNPs identified in this study are indicated on each chromosome (left and right, respectively), and the flanking regions of previously reported CR gene loci and QTLs identified in this study are denoted by black and red bars, respectively, both on and to the right of each chromosome.

Table 2 .
Distribution 1 # SNP, SNP count located on A and C chromosomes or scaffolds.2SNP/Mb,SNP density per million base pairs.NA, not applicable.
Note: The QTLs are denoted as major or minor in parentheses below each QTL name based on peak |∆(SNPindex)| values.Gene IDs and names were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/,accessed on 1 December 2023), where gene names denote descriptions of gene functions.Overlapping QTLs denote the same QTL identified in bulks tested with different pathotypes.
Note: Gene IDs, symbols, and names were obtained from the NCBI database (https://www.ncbi.nlm.nih.gov/,accessed on 1 December 2023), where gene names denote descriptions of gene functions.