Comparative transcriptome analysis of canola carrying a single vs stacked resistance genes against clubroot

Pyramiding resistance genes may expand the efficacy and scope of a canola variety against clubroot (Plasmodiophora brassicae), a serious threat to canola production in western Canada. However, the mechanism(s) of multigenic resistance, especially the potential interaction among clubroot resistance (CR) genes, are not well understood. In this study, transcriptome was compared over three canola (Brassica napus L.) inbred/hybrid lines carrying a single CR gene in chromosome A03 (CRaM , Line 16) or A08 (Crr1rutb , Line 20), and both genes (CRaM +Crr1rutb , Line 15) inoculated with a field population (L-G2) of P. brassicae pathotype X, a new variant found in western Canada recently. The line16 was susceptible, while lines 15 and 20 were partially resistant. Functional annotation identified differential expression of genes (DEGs) involved in biosynthetic processes responsive to stress and regulation of cellular process; The Venn diagram showed that the partially resistant lines 15 and 20 shared 1,896 differentially expressed genes relative to the susceptible line 16, and many of these DEGs are involved in defense responses, activation of innate immunity, hormone biosynthesis and programmed cell death. The transcription of genes involved in Pathogen-Associated Molecular Pattern (PAMP)-Triggered and Effector-Triggered Immunity (PTI and ETI) was particularly up-regulated, and the transcription level was higher in line 15 (CRaM + Crr1rutb ) than in line 20 (Crr1rutb only) for most of the DEGs. These results indicated that the partial resistance to the pathotype X was likely conferred by the CR gene Crr1rutb for both lines 15 and 20 that functioned via the activation of both PTI and ETI signaling pathways. Additionally, these two CR genes might have synergistic effects against the pathotype X, based on the higher transcription levels of defense-related DEGs expressed by inoculated line 15, highlighting the benefit of gene stacking for improved canola resistance as opposed to a single CR gene alone.


Introduction
Clubroot, caused by the soil-borne protist Plasmodiophora brassicae Woronin, is an important disease of brassica crops worldwide (Dixon, 2009), and a serious threat to canola production in western Canada where the crop contributes billions to the annual economy (Rempel et al., 2014;Strelkov and Hwang, 2014;Gossen et al., 2015).In the realm of managing the disease on canola, genetic resistance has proven to be both efficacious and practical, particularly when used with extended crop rotations exceeding a 2-year break from a canola crop that drastically diminishes P. brassicae inoculum in the soil (Peng et al., 2014;Peng et al., 2015).Despite the effectiveness, repeated uses of a single clubroot resistance (CR) gene has led to rapid resistance erosion in Canada, a consequence of the heightened selection pressure from a diverse pathogen population (Strelkov et al., 2018;Mcdonald et al., 2020).Novel P. brassicae pathotypes have been identified since the introduction of resistant canola cultivars in 2009, including pathotype X.These 'new' pathotypes defeated the single CR genes used in earlier resistant canola cultivars released in Canada (first-generation) within 3-4 years (Strelkov et al., 2018;Sedaghatkish et al., 2019).
It has been recognized that CR resources are limited (Hirai, 2006).To date 18 major CR loci have been reported, with most of them being originated from European turnip (B.rapa).Many of these CR loci belong to the toll-interleukin receptor (TIR)nucleotide-binding site (NBS) family (Karim et al., 2020;Mehraj et al., 2020).Two cloned CR genes, CRa and Crr1, also encode TIR-NBS-LRR (leucine-rich repeat) proteins (Ueno et al., 2012;Hatakeyama et al., 2013).In Canada, one of the CR genes derived from the winter rapeseed cultivar 'Mendel' (A03) appeared to be present in most first-generation resistant canola cultivars post-2009 (Fredua-Agyeman andRahman, 2016).While effective against initial five pathotypes found in Canada, this gene was observed to have lost the effectiveness in some field by 2013 (Strelkov et al., 2018;Karim et al., 2020).Additional CR genes (A08) derived from rutabaga (Hasan and Rahman, 2016) were introduced later, which showed resistance or moderate resistance to most of the novel pathotypes, including pathotype X (Tonu et al., 2023).The rutabaga-derived CR gene(s) have also been pyramided with those on A03 (second-generation CR) for a broader range of efficacy.
While stacking resistance genes has proven effective against multiple pathogen races in other crops, including rice against bacterial blight (Li et al., 2001), potentially prolonging the useful life of individual genes, the situation in canola/rapeseed, especially concerning resistance to blackleg (Leptosphaeria maculans), appears complex (Balesdent et al., 2022).In comparison, combining major R genes with quantitative resistance consistently enhanced the resistance durability (Brun et al., 2010).Stacking three CR genes in Chinese cabbage (B.rapa) expanded the resistance against multiple P. brassicae pathotypes (Matsumoto et al., 2012), although the resistance mechanisms associated with CR-gene stacking remain unexplored.In a recent study, Tonu et al. (2023) reported that canola varieties carrying stacked CR genes (A03 and A08) had greater resistance durability than those carrying either CR gene alone, when exposed repeatedly to a field P. brassicae pathotype X population.It was unclear, however, whether these stacked CR genes would provide more sophisticated resistance mechanisms than a single CR gene alone against novel pathotypes virulent towards the first-generation resistant canola cultivars.
Studying molecular mechanisms, particularly those conferred by both single and stacked CR genes, may help understand the basis of resistance for developing genetic resources efficiently.As single and stacked genes provide resistance against particular pathotypes/ strains of the pathogen, understanding the modes of action may allow for a better design in CR breeding, enabling the development of canola cultivars with tailored resistance profiles.This type of study may also help identify gene combinations that confer more durable resistance, aiding in judicious deployment of genetic resources.For instance, breeders can use complementary CR genes which offer the best protection against prevailing pathogen populations based on the understanding of synergistic or additive effects of gene stacking.
Transcriptome analysis, aimed at identifying differentially expressed genes (DEGs) between susceptible and resistant plants, serves as a powerful tool for unraveling the biological pathways in the host against pathogen attacks (Orshinsky et al., 2012;Soneson and Delorenzi, 2013).The comparative analysis can also shed light on defense pathways triggered by specific CR genes or their combinations.The primary objective of this study was to the novel pathotype X as a model to unravel the molecular mechanisms conferred by double (A03 and A08) vs. single (A03 or A08) CR genes, probing for potential interactions in mediating the resistance.This was achieved through: 1) Identifying DEGs between canola lines carrying single and double CR genes inoculated with the pathotype X; 2) conducting functional analysis of DEGs to identify biological pathways crucial for the resistance; and 3) determining distinct defense mechanisms and transcription levels activated by specific CR genes or gene combinations.The information holds practical value for designing CR-gene stacking in canola breeding and for aiding deployment of cultivars equipped with optimized CR genes for enhanced resistance performance and durability.

Plant materials, inoculation and disease assessment
Three commercial canola inbreed/hybrid lines/varieties, designated as line 15, 16 and 20 in this study, were provided by Nutrien Ag Solutions, Saskatoon, Saskatchewan, Canada.Line 16 carries a single CR gene in Chromosome A03 that was derived from the winter oilseed rape cultivar Mendel.This CR gene was originally from a fodder turnip (B.rapa) used in the European Clubroot Differentials (ECD) designated as ECD4.This ECD differential was believed to carry three CR genes but two of them might have been lost during backcrossing with B. napus in production of the oilseed rape cultivar Mendel (Diederichsen et al., 2006).Fredua-Agyeman and Rahman (2016) showed that this CR locus was physically close to that of CRa/CRb Kato located in Chromosome A03 of B. rapa (Matsumoto et al., 1998;Kato et al., 2012;Ueno et al., 2012;Kato et al., 2013).Recent evidence showed that this CR gene from 'Mendel' is identical to CRa (Hu et al., 2024).
Line 20 carries CR gene(s) in Chromosome A08 derived from a variety of rutabaga (B.napus var.napobrassica); Composite Interval Mapping analysis showed that this CR locus can be flanked by the SSR markers sS1702 and A08_5024 (Hasan and Rahman, 2016), which would indicate a region in A08 where the CR gene Crr1 is also located (Suwabe et al., 2012;Hatakeyama et al., 2013).Based on fine mapping of a 1.6 cM region on A08, Suwabe et al. (2012) provided the evidence that the Crr1 region would carry two CR genes, i.e., Crr1a and Crr1b in a very close range.However, it is unclear whether line 20 carries only a single or both alleles in A08.To avoid nomenclature confusion, the CR gene(s) carried in lines 16 and 20 were respectively referred to as CRa M and Crr1 rutb in this paper to reflect their origins and relatedness to CRa and Crr1.
Line 15, on the other hand, was produced by the hybridization of inbred lines 16 and 20, which would carry at least two dominant CR genes derived from A03 and A08, respectively, based on marker analysis (data not shown).However, it remains undetermined whether a single or two CR genes were present in A08 of the hybrid.These inbred/hybrid lines (15, 16 and 20) were used to assess and compare the transcriptomic responses of CRa M and Crr1 rutb , both individually and in combination.
A field collection (L-G2) of P. brassicae pathotype X was used throughout the study to inoculate all canola lines.Pathotype X was the first strain characterized to be virulent on the first generation of resistant canola cultivars in Canada (Strelkov et al., 2016).The clubroot reaction of the three lines to L-G2 was compared against a mock treatment as a negative control and the susceptible canola cultivar Westar as a positive control.
A resting spore suspension of L-G2 was prepared and mixed into soil-less Sunshine #3 potting mix (pH = 6.2.SunGro Horticulture, Vancouver, BC) to reach 1 x 10 6 spores/g soil concentration.Subsequently, each line was planted in infested growth medium (20 cm in diameter, 15 cm deep) with 20 seeds per pot.The pots were placed in a growth room set at 22°C/16°C (day/night) with a 16 h photoperiod until root sampling.Following seeding, the growth medium in each pot was initially saturated for one week and then maintained in a moist condition through regular watering.
Root tissues were collected from each plant at 14 days post inoculation (dpi), with 15 plants per cultivar per pot (replicate) and three replicates in total for each line.The soil was rinsed off from the roots with tap water to preserve the whole root system, and then the entire root system was cut, flash frozen in liquid nitrogen and stored at -80°C until use.Five remaining plants from each pot were maintained in the growth room and assessed for clubroot severity on a standard 0-3 scale at 35 dpi to confirm the success of inoculation.A disease severity index (DSI) was calculated for the five remaining plants assessed for each replicate (Tonu et al., 2023).

RNA-seq analysis
RNA extraction from root samples was performed using the RNeasy Plant Mini Kit (Qiagen, Toronto, ON) following the manufacturer's instructions, with DNase digestion using the RNase-Free DNase Set from Qiagen.The quality and concentration of RNA were evaluated using the Experion RNA StdSens Analysis Kit on the Experion automated electrophoresis system (Bio-Rad, Montreal, QC) and Nanodrop 2000c (Fisher Scientific, Toronto, ON) to ensure sufficient RNA quality and quantity for cDNA library preparation.
cDNA libraries were prepared using the TruSeq RNA Sample Preparation Kit v2 (Illumina; San Diego, CA).The quality and concentration of cDNA were assessed also with the Experion DNA Analysis Kit and Nanodrop 2000c.These libraries were sent to the Geńome Quebec Innovation Centre, McGill University for RNA sequencing on the Illumina Hiseq 2500 platform.The raw data of RNA-seq has been deposited at the National Center for Biotechnology Information (NCBI) under this ID number: SUB14143825.
RNA-seq data processing and DEG analysis followed the protocol outlined by Chu et al. (2014).Briefly, CLC Genomics Workbench v10.1.1 (Qiagen) was used to process and analyze raw sequencing data (FASTQ files); the reads underwent quality control checks, followed by trimming to remove Illumina adapters and low-quality reads.Clean reads were aligned to the reference genome of Brassica napus (v4.1; http:// brassicadb.org/brad).The level of gene expression was quantified in Reads Per Kilobase of exon model per Million reads (RPKM).There are concerns about incorrect uses of RPKM for some samples (Evans et al., 2018;Zhao et al., 2020), and we have been keeping all plants grown under the same environment and using the same RNA-seq protocol on all samples.To minimize unnecessary data analyses due to excessively intricate interactions, DEGs were identified at RPKM |>4| to better focus on most highly differentiallyexpressed genes relevant to CR functions.Such strategy has been used in other RNA-seq data analyses (Rosli et al., 2013;Chhiba et al., 2017;Iqbal & Kumar, 2022;Jayavelu et al., 2023).A false discovery rate (FDR) threshold was set at P ≤ 0.01.DEG comparisons were made between the inoculated treatment and the negative control (mock) for each canola variety/line.

Annotation of DEGs
The list of DEGs underwent gene ontology (GO) annotation using Blast2GO Pro (Conesa et al., 2005) with BLASTx algorithms against the non-redundant protein database provided by the National Center for Biotechnology Information (http:// www.ncbi.nlm.nih.gov).DEGs were classified into three GO classes: biological process, molecular function and cellular component.The enrichment analysis was conducted by using a built-in tool, Fisher's Exact Test, in Blast2GO Pro.Heatmaps and gene clustering were performed using a R package (https:// CRAN.Rproject.org-/package=pheatmap),and Log-transformation was applied prior to plotting.The annotations were further analyzed using MapMan software (Thimm et al., 2004) to categorize abiotic and biotic pathways associated with selected DEGs.

Verification of RNA-seq data quality
To ensure the reliability of RNA-seq data, droplet digital PCR (ddPCR) was conducted using a QX200 ™ System (Bio-Rad).Reagents and consumables were also from the same supplier, including droplet generator oil, DG8TM cartridges and gaskets, droplet reader oil and ddPCR EvaGreen supermix.For each sample, 20-µl reaction mix containing 10-µl EvaGreen supermix, 1 µl each of forward and reverse primers (22.5 µM), 2-µl sample complementary DNA and 6-µl double distilled H 2 O was partitioned into aqueous droplets in oil via a droplet generator.The droplet suspension was transferred to a 96-well PCR plate, and a thermocycling process carried out in a conventional thermal cycler (Thermo Fisher Scientific Canada, Toronto, ON) at 95°C for 5 min, 95°C for 30 s (40 cycles), 56-62°C (depending on specific genes) at ramp rate: 2°C/s for 60 s (40 cycles), 4°C for 5 min, 90°C for 5min and then, 4°C infinite.Then the PCR plate was determined for the fraction of PCR-positive droplets in each sample on a droplet reader, and the data analyzed using the QuantaSoft ™ Software.The expression patterns of 10 selected DEGs, involved in different biological process, were analyzed using droplet digital PCR (ddPCR).The primers for amplifying these genes were listed in Supplementary Table 1.

Statistical analysis
DSI data were transformed (arcsine square root) and normal distribution confirmed based on Shapiro-Wilk Test (PROC UNIVARIATE) using SAS v9.3 (SAS Institute, Cary, NC, USA) prior to statistical analyses.The homogeneity of variance was assessed using Levene's Test.Fisher's Protected LSD was used to separate treatment means when ANOVA was significant (P ≤ 0.05).

Clubroot severity on inbred and hybrid canola lines
In this investigation, we evaluated the impact of CR genes CRa M and Crr1 rutb , individually and in combination, against the L-G2 collection of pathotype X. Line 16, carrying the single CR gene CRa M , was susceptible with severe clubroot symptoms visible at 35 dpi.In contrast, lines 20 and 15, carrying Crr1 rutb and both CRa M and Crr1 rutb , respectively, displayed partial resistance with only slight to moderate symptoms (Figures 1A, B).The disease severity indices (DSIs) were 68% for line 16, 26% for line 20, and 24% for line 15 (Figure 1C).

RNA-seq analysis
RNA-seq analysis was conducted on root samples from the three B. napus lines carrying CRa M , Crr1 rutb , or both CR genes collected at 14 dpi with three biological replicates.A total of 1.08 billion reads were generated through 100 bp paired-end sequencing from 18 cDNA libraries (Supplementary Table 2).Approximately one-third of the total reads were obtained from each line: 368,201,810 from line 15 (166,577,030 for control and 201,624,780 for inoculated, respectively), 337,590,744 from line 16 (175,260,710 for control and 162,330,034 for inoculated), and 377,931,694 from line 20 (189,199,162 for control and 188,732,532 for inoculated).
The GC content ranged from 46-48%, and all libraries contained >80% of sequences with an average PHRED score >30 (Q30%), ensuring the quality and accuracy of the sequencing data for further analyses and gene expression comparisons.Sequences were trimmed to remove impurities (i.e., low-quality reads, adapters etc.) before the analysis of RNA-seq data.Around 70% of pairedsequence reads for each cDNA library were mapped to the B. napus reference genome version 4.1 (http://brassicadb.org).A summary of RNA-seq data quality checks and analyses is presented in Supplementary Table 2.
Principal Component Analysis (PCA) revealed that biological replicates were grouped together for each treatment in the same dimension, indicating no major variations in the RNA-seq data within each treatment (Figure 2; Supplementary Figure 1).Additionally, the three biological replicates of control and inoculated treatments were closely clustered for line 16, suggesting no major differences between the control and inoculated samples.Conversely, the replicates of control and inoculated treatments for lines 15 and 20 were in different dimensions, indicating significant differences in the RNA-seq data.

Annotation of DEGs
The threshold for the identification of DEGs was set at an absolute fold change of ≥ 4 and a False Discovery Rate (FDR) Pvalue ≤ 0.01.In line 15, a total of 6,000 DEGs were identified, comprising 3,163 up-regulated genes and 2,837 down-regulated genes in the inoculated plants compared to the control.For line 16, 953 DEGs were identified, with 546 up-regulated and 407 downregulated genes.Line 20 exhibited 3,613 DEGs, including 2,286 up-regulated and 1,327 down-regulated genes (Figure 3; Supplementary Table 3).
The Blast2GO analysis revealed that up-and down-regulated DEGs were mainly involved in cellular processes, metabolic processes, response to stimuli, and biological regulation.The total number of up-and the DEGs across the three canola lines are shown in Supplementary Figure 2. GO analysis revealed DEGs associated with Biological Processes, Molecular Functions and Cellular Components; most DEGs from inoculated plants were involved in the biological processes of biosynthesis, response to stress and regulation of cellular processes (Figure 4A).Line 15 had the most DEGs assigned to each biological process compared to line 20 or line16.For molecular functions, a substantial proportion of DEGs were associated with organic cyclic compound binding, heterocyclic compound binding, ion binding and transcription factor activity (Figure 4B).Regarding cellular components, DEGs were relatively abundant in intracellular parts, followed by intracellular organelles and membrane-bounded organelles (Figure 4C).Furthermore, enrichment analysis was performed on identified DEGs using Fisher's Exact Test in Blast2GO, with The top 20 enriched GO terms identified in Supplementary Figure 3.The highest number of DEGs were linked to cellular anatomical entities.Interestingly, DEGs associated with cellular processes ranked as the second most common in lines 15 and 20, but were not present in significant numbers in line 16.Lines 15 and 20 (partially resistant) appeared to share many common DEGs and functions in similar biological processes, while the susceptible line 16 lacked most of these DEGs against the infection by pathotype X; Venn diagram analysis showed that 1,896 DEGs were shared between lines 15 and 20 (Figure 5A; Supplementary Table 4).The enrichment analysis of these DEGs revealed their involvement in the biological process of cell communication and responses to hormones and defense responses, with some of the DEGs functioning in programmed cell death and the regulation of systemic acquired resistance (SAR) (Figure 5B).
Further analysis of DEGs involved in defense response categories, shared by lines 15 and 20 but not line 16, showed the activation of genes involved in signaling pathways PTI and ETI, respectively (Tables 1, 2; Supplementary Tables 5, 6).Additionally, PTI signaling pathways involved pattern recognition receptors (PRRs), mitogen-activated protein kinase (MAPK) cascade, and transcription factors (TF).Key PRRs included wall-associated receptor kinase-like 10, cysteine-rich receptor-like protein kinase Principal component analysis of RNA-Seq data over all replicates of root samples from three canola lines inoculated with L-G2 of Plasmodiophora brassicae pathotype X. Overlapping regions correspond to the number of DEGs identified in different biological replicates.Red marks represent line 15, purple marks for line 16, and green marks for line 20.M represents mock inoculation, and L represents L-G2 of pathotype X inoculation.Samples used for RNAseq were collected at 14 dpi.
In the MAPK cascade, 95% of 61 associated genes were intensely up-regulated, indicating a strong response caused by the inoculation.The main TF involved in PTI, including probable WRKY TF70, 38 and 51 isoform ×1, as well as NAC TF29, showed the largest fold changes in line 15, followed by line 20, and the smallest or occasionally insignificant changes in line 16 (Table 1).Collectively, these results indicated that PTI signaling pathways were activated in lines 15 and 20, but not in line 16.
ETI pathways implicated disease resistance proteins, such as Enhanced Disease Susceptibility 1 (EDS1) and isochorismate synthase, with 13 associated DEGs identified.Their transcription levels were notably up-regulated more in lines 15 and 20 compared to line 16.(Table 2; Supplementary Figure 4B), including two isochorismate synthases and two isochorismate synthase chloroplastic-like proteins.This would suggest the differential activation of ETI in lines 15 and 20 by pathotype X, but not in line 16.When DEGs involved in abiotic and biotic signaling responses were analyzed for inoculated plants using the MapMan, more MAPK-associated genes, TF and PR genes were up-regulated in lines 15 and 20 than in line 16 (Supplementary Figures 5-7).
Furthermore, a total of 224 DEGs were identified as common between mock (M)-15/inoculated (I)-15 and M-16/I-16, while only 80 DEGs were shared between M-20/I-20 and M-16/I-16.Additionally, 241 DEGs were shared among M-15/I-15, M-16/I-16 and M-20/I-20.Subsequently, GO analysis was conducted on these shared DEGs (Supplementary Figures 8-10).DEGs associated with molecular functions were categorized into small molecule/ organic cyclic compound binding, transferase, catalytic and hydrolase activities.Biological processes were linked to metabolic processes, cellular responses to stimuli and biological regulation.

Validation of RNA-seq data quality
The expression patterns of 10 selected DEGs, involved in pathogen recognition, signal transduction, transcriptional regulation and disease responses, were analyzed using droplet digital PCR (ddPCR).To elucidate the alterations in gene expression following inoculation, relative fold changes in expression levels compared to the mock control were computed based on data obtained from both RNA-Seq and ddPCR (Figure 6).
The results showed that the expression patterns of the chosen DEGs were similar between ddPCR and RNA-Seq, although the relative fold changes may vary slightly for certain genes between the two methods.

Discussion
This study used canola lines carrying different CR genes (CRa M and/or Crr1 rutb ) against P. brassicae pathotype X (L-G2) as a model system to gain insights into molecular clubroot resistance responses.The evaluation of disease severity illustrated that line 16, with the single CR gene CRa M alone, exhibited susceptibility, while lines 20 and 15, carrying Crr1 rutb alone and CRa M + Crr1 rutb , respectively, displayed moderate levels of resistance.This observation aligns with previous reports on the effectiveness of certain CR genes and gene combinations in conferring resistance to specific P. brassicae pathotypes (Strelkov et al., 2018;Sedaghatkish et al., 2019;Tonu et al., 2023).
Stacking the CR genes CRa M and Crr1 rutb (line 15) provided only partial resistance to the filed population of pathotype X (L-G2) in the current study.It is possible that this partial resistance is derived from Crr1 rutb in line 20, as the inbred line carrying CRa M alone was susceptible.In a prior study, Tonu et al. (2023) reported that line 20 produced in different batches may vary in resistance to pathotype X, with some batches showing almost immunity.Mapping of Crr1 rutb (Hasan and Rahman, 2016) showed that this CR gene is located in the region in A08 where Crr1 (Suwabe et al., 2012;Hatakeyama et al., 2013) and Rcr3 (Karim et al., 2020) also reside.Rcr3 can be flanked by two SNP markers that are 231.6Kb apart; this interval contains 32 genes, three of which (Bra020951, Bra020974, Bra020979) are associated with disease resistance.Previous studies also showed that Crr1 possibly consists of two loci; a major locus, Crr1a, which encodes a TIRNB-LRR protein, and a minor locus, Crr1b (Suwabe et al., 2012;Hatakeyama et al., 2013).Although it cannot be determined unequivocally whether Crr1 rutb is homologous to Crr1/Rcr3 due to lacking sequence information on Crr1 rutb and Rcr3, it is possible that more than one allele are involved in Crr1 rutb for its strong resistance efficacy against pathotype X (Tonu et al., 2023).It is also possible that one or more of the Crr1 rutb alleles had been lost during the production of lines 15 and 20, as some batches of line 20 showed almost immunity to pathotype X (Tonu et al., 2023).This loss of allele during hybridization has also happened with Crr1 (Fredua- Agyeman et al., 2018).Despite the partial resistance, line 15, carrying the stacked CR genes CRa M and Crr1 rutb , showed stronger resistance resilience than line 16 or line 20 in a previous study (Tonu et al., 2023), with little evidence for resistance erosion or soil inoculum buildup over five cycles of continuous exposure to pathotype X.
The RNA-seq analysis provided a comprehensive understanding of the transcriptomic changes in response to the pathotype X (L-G2) inoculation across these three B. napus inbred/ hybrid lines.The reason for extraction of root-tissue RNA at 14 dpi was that the cortical infection by P. brassicae would have been completed (McDonald et al., 2014;Irani et al., 2018).Despite the absence of visible symptoms, the expression of genes associated with pathogen recognition, signal transduction and defense responses would have already peaked (Chu et al., 2014).This early time frame is particularly suitable for the RNA-seq analysis, providing valuable insights into plant defense mechanisms against the pathogen attack.
The high-quality sequencing data and subsequent analysis including PCA and ddPCR, ensured the reliability and reproducibility of results.The PCA results indicated tight clustering of biological replicates within each treatment, affirming the consistency of the RNA-seq data.The ddPCR results further supported the transcriptomic data; the generally consistent gene expression patterns observed in both platforms underscored the robustness of DEGs identified.This careful experimental design enhanced the confidence in the transcriptomic changes observed, allowing for meaningful interpretation of the results.The slight differences in gene expression levels (fold changes) between RNA- The comparisons were between mock (M) and inoculated (I) samples.
seq and ddPCR can be due to more sensitive detection of lowabundance transcripts by ddPCR or lower susceptibility of RNA-seq to issues related to PCR efficiency.It is also noteworthy that few distinctions were observed among the DEGs between the control and inoculated line 16 carrying only CRa M .In contrast, notable disparities were evident between the control and inoculated lines 15 (CRa M +Crr1 rutb ) or line 20 (Crr1 rutb ), culminating in a moderately resistant response to L-G2.The annotation of DEGs revealed substantial transcriptional changes in response to infection by P. brassicae.The stringent criteria for DEG identification (fold change ≥ 4, FDR ≤ 0.01) ensured the selection of robust candidates.The GO analysis provided insights into the Biological Processes, Molecular Functions, and Cellular Components associated with the DEGs.Notably, the abundance of DEGs involved in biosynthetic processes, response to stress and regulation of cellular processes highlights the complex interplay of molecular events during the plant-pathogen interaction.The comparison of DEG profiles among these three lines shed light on the common and distinct transcriptional responses; lines 15 and 20, both partially resistant to L-G2, shared a substantial number of DEGs, suggesting some common defense mechanisms.Conversely, line 16 exhibited a distinct transcriptional profile, which failed to mount effective resistance.
Plant defense responses involve pathogen recognition, signal transduction and defense activation (Andersen et al., 2018).In the first stage, receptor-like proteins (RLPs) and receptor kinases (RKs) recognize immunogenic molecular patterns of pathogens (Zhou and Zhang, 2020).Up-regulation of several RLPs and RKs, including wall-associated receptor kinase-like 10, cysteine-rich receptor-like protein kinase 36, receptor-like protein 12, and Gtype lectin S-receptor-like serine threonine-protein kinase rlk1, was observed in lines 15 and 20, but not in line 16, upon infection (Table 1; Supplementary Figure 4A).Additionally, the downstream signaling event of pattern recognition receptor (PRR)-based pathogen recognition often involves the mitogen-activated protein kinase (MAPK).Activation of MAPK was strongly evident for lines 15 and 20, with 95% of 61 associated genes being intensely upregulated and 5% of them down-regulated.In comparison, only 79% of these genes were up-regulated and 21% down-regulated in line 16, with much smaller fold changes (Supplementary Table 5).
This upregulation of PRRs, MAPK cascade and TF was generally more pronounced in lines 15 and 20 than in line 16.This strongly indicates the activation of PTI, suggesting robust pathogen recognition and potent PTI-triggered defense responses (Bigeard et al., 2015) against L-G2.In contrast, apparently weaker responses of line 16 in MAPK may suggest a compromised defense mechanism.The high-fold changes of probable WRKY transcription factor would contribute to the activation of PTI responses.In Arabidopsis, activated MAPKs phosphorylate substrate proteins, including ethylene-responsive TF AtERF104 and AtWEKY70 (Bigeard and Hirt, 2018).In the current study, the fold changes for three genes encoding possible WRKY transcription factor 70 were the highest in line 15, lower in line 20 and below the threshold level of significance in line 16 (Table 1; Supplementary Figure 4A).This emphasizes the strong activation of PTI associated with the double CR genes.There may be The comparisons were between mock (M) and inoculated (I) samples.complementary interaction between two CR genes (Yu et al., 2022) in response to specific strains of P. brassicae, although these enhanced molecular responses did not translate into substantially reduced disease severity (Figure 1B).The analysis of defense responses also revealed the activation of ETI signaling pathways in lines 15 and 20, but not in line 16.The identification of DEGs encoding plant intracellular immune receptors, known as resistance proteins, as well as isochorismate synthases (ICS), supported the notion of ETI.The resistance proteins are believed to trigger ETI upon the detection of pathogen effectors; it appears that EDS1 and ICS are the key components in ETI downstream responses.Up-regulation of EDS1 and two ICS genes was the highest in line 15, followed by line 20.Additionally, DEGs associated with hypersensitive response (HR), characteristic of ETI, were also identified in lines 15 and 20.EDS1 is a general signaling component in the down-stream of NBS-LRR, and is activated by the presence of cleaved products released by resistance proteins (Zhou and Zhang, 2020).EDS1 can combine with phytoalexin deficient 4 (PAD4) or senescence-associated gene 101 to form a complex that induces salicylic acid (SA) biosynthesis or results in cell death (Joglekar et al., 2018).Up-regulation of EDS1 was the highest in line 15 and lowest in line 16 (Table 2; Supplementary Figure 4B).The high transcription level of EDS1 indicated that it was activated following infection of lines 15 and 20, but not line 16.In Arabidopsis, SA is synthesized by two ICS genes, ICS1 and ICS2 (Chen et al., 2009), which were up-regulated in lines 15 and 20 (BnaA08g22340D and BnaC08g18420D), but not in line 16 (Table 2; Supplementary Figure 4B).In total, thirteen disease resistance-related DEGs involved in ETI were identified, with more pronounced up-regulation in the resistant canola, especially line 15.Validation of the RNA-Seq data using ddPCR.The golden line in the chart is intended solely to differentiate the expression of genes in the ddPCR testing from those represented in the RNA-seq data (bar chart).Each ddPCR value is independent, similar to that in the RNA-seq data.
A distinguishing characteristic of ETI that sets it apart from PTI is HR (Hatsugai et al., 2017).In our investigation across the three canola lines, a total of 105 genes associated with HR were identified.Notably, 90% of them exhibited up-regulation, while only 10% showed down-regulation in lines 15 and 20.In contrast, line 16 displayed a slightly different pattern, with 80% of the genes upregulated and 20% down-regulated.However, the fold changes in these genes were markedly higher in lines 15 and 20 compared to line 16 (Supplementary Table 6).Collectively, these findings signify the activation of both PTI and ETI for lines 15 and 20 in response to L-G2 inoculation, whereas line 16 failed to exhibit robust activation of these immune responses.Upon the activation of PTI and ETI, the transcriptional reprogramming of plant defense response genes ensued, involving a range of biological processes, including hormone biosynthesis, regulation of SAR and cell wall organization.These processes possibly contributed to the defense mechanism against the L-G2.The recently published study by Wang et al. (2023) identifies WeiTsing (WTS) as a broadspectrum CR gene inhibiting the colonization of P. brassicae in the root stele.Activation of WTS induces the expression of many defense genes associated with plant immunity, notably increasing the levels of RLPs and wall-associated kinases (WAKs).These findings suggest that the activation of WTS by P. brassicae likely simultaneously triggers both PTI and ETI signaling pathways, thereby inducing plant immune responses.

Conclusion
Understanding the molecular mechanisms of CR genes is crucial for efficiently harnessing genetic resources.Single and stacked CR genes offer resistance against specific P. brassicae pathotypes, enabling the development of breeding strategies with tailored CR profiles.The comparative transcriptome analysis of inbred/hybrid canola lines carrying CRa M , Crr1 rutb , and CRa M + Crr1 rutb provides a comprehensive understanding of the hostpathogen interaction involving pathotype X (L-G2).While few distinctions were observed for DEGs in control and inoculated line 16 (CRa M ), noticeable differences were found for lines 20 (Crr1 rutb ) and 15 (CRa M + Crr1 rutb ).It appears that the observed resistance may predominantly stem from Crr1 rutb for both lines 15 and 20, given their substantial overlap in DEGs, even though transcription levels were often higher in the double CR-gene line 15.The up-regulation of several RLPs and RKs suggests the activation of PTI during pathogen recognition.Furthermore, activation of ETI is also suggested by the pronounced upregulation of EDS1, ICS and HR-related genes, especially in line 15.These findings indicated both PTI and ETI contribute to the resistance CRa M and Crr1 rutb against pathotype X. Upon pathogen recognition, transcriptional reprogramming seems to occur, with several biological processes being activated, including hormone biosynthesis, SAR regulation, and cell-wall organization.This study provides information on intricate molecular interactions involved in CR-mediated resistance and contributes insights for future breeding strategies aimed at sustainable clubroot resistance and management.
reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

FIGURE 3
FIGURE 3Identification of DEGs filtered by an absolute fold change > 4 and Pvalue ≤ 0.01 between inoculated and control root samples of three canola lines.The red portions at the top of each bar represent down-regulated genes in the inoculated plants relative to the control, and the blue portions at the bottom represent upregulated genes.
FIGURE 5 Distribution of DEGs and their biological process analysis for three canola lines: (A) Venn diagram of the distribution of DEGs; (B) Enrichment analysis of the DEGs in disease resistance-related biological processes shared by lines 15 and 20, but not 16.

TABLE 1
List of top 10 Pattern Recognition Receptors (PRRs) and transcription factors involved in Pattern-Triggered Immunity (PTI) based on fold changes in gene regulation associated with inoculation of three canola lines.

TABLE 2
The primary genes involved in Effector-Triggered Immunity (ETI) based on fold changes in gene regulation associated with the inoculation of three canola lines.