Leonardo Galindo-González
Victor Manolii
Sheau-Fang Hwang
Stephen E. Strelkov- Department of Agricultural, Food & Nutritional Science, University of Alberta, Edmonton, AB, Canada
Clubroot, caused by the obligate parasite Plasmodiophora brassicae, is an important disease of the Brassicaceae and poses a significant threat to the $26.7 billion canola/oilseed rape (Brassica napus) industry in western Canada. While clubroot is managed most effectively by planting resistant host varieties, new pathotypes of P. brassicae have emerged recently that can overcome this resistance. Whole genome analyses provide both a toolbox and a systemic view of molecular mechanisms in host-pathogen interactions, which can be used to design new breeding strategies to increase P. brassicae resistance. We used RNA-seq to evaluate differential gene expression at 7, 14 and 21 days after inoculation (dai) of two B. napus genotypes with differential responses to P. brassicae pathotype 5X. Gall development was evident at 14 dai in the susceptible genotype (the oilseed rape ‘Brutor’), while gall development in the resistant genotype (the rutabaga (B. napus) ‘Laurentian’) was limited and not visible until 21 dai. Immune responses were better sustained through the time-course in ‘Laurentian’, and numerous genes from immune-related functional categories were associated with salicylic acid (SA)-mediated responses. Jasmonic acid (JA)-mediated responses seemed to be mostly inhibited, especially in the resistant genotype. The upregulation of standard defense-related proteins, like chitinases and thaumatins, was evident in ‘Laurentian’. The enrichment, in both host genotypes, of functional categories for syncytium formation and response to nematodes indicated that cell enlargement during P. brassicae infection, and the metabolic processes therein, share similarities with the response to infection by nematodes that produce similar anatomical symptoms. An analysis of shared genes between the two genotypes at different time-points, confirmed that the nematode-like responses occurred earlier for ‘Brutor’, along with cell metabolism and growth changes. Additionally, the susceptible cultivar turned off defense mechanisms earlier than ‘Laurentian’. Collectively, this study showed the importance of SA in triggering immune responses and suggested some key resistance and susceptibility factors that can be used in future studies for resistance breeding through gene-editing approaches.
Introduction
Clubroot, caused by the obligate parasite Plasmodiophora brassicae Wor., is a soilborne disease of the Brassicaceae and represents a major threat to canola/oilseed rape (Brassica napus L.) production. In Canada, where the canola crop contributes $26.7 million annually to the national economy (Canola Council of Canada – Industry overview, 2017), clubroot has been spreading rapidly through the Prairies and has become a significant challenge for many farmers (Strelkov et al., 2015; Strelkov et al., 2019). Disease development is associated with the production of large numbers of pathogen resting spores, which can persist in the soil for many years and serve as inoculum for the infection of subsequent crops (Howard et al., 2010; Hwang et al., 2012). Strategies for the management of clubroot include crop rotation, sanitization of field equipment, and the application of fungicides and soil amendments (Donald and Porter, 2009). Unfortunately, many of these strategies are not practical or cost-effective (Hwang et al., 2014), and as such farmers have relied mostly on the planting of clubroot resistant (CR) cultivars. In recent years, however, new pathotypes of P. brassicae, capable of overcoming the resistance in most CR canola cultivars, have been identified with increasing frequency (Strelkov et al., 2016b; Strelkov et al., 2018). The emergence of these new pathotypes highlights the need for additional sources of resistance and a better understanding of the host response to infection.
The clubroot pathogen invades plants through the roots, causing a primary infection in the root hairs and then a secondary infection of the cortical tissues (Kageyama and Asano, 2009; Hwang et al., 2012). The secondary plasmodia that develop in the cortex trigger host regulation of growth related hormones (auxins, brassinosteroids, cytokinins), resulting in hypertrophy and hyperplasia that lead to gall formation (Kageyama and Asano, 2009; Ludwig-Müller et al., 2009; Schuller et al., 2014; Ludwig-Müller et al., 2017). Cell enlargement has a two-fold benefit for P. brassicae, housing the enlarged plasmodia which divide concomitantly with the host cells, and providing increased nutrients for pathogen growth and reproduction (Ludwig-Müller et al., 2009; Jahn et al., 2013; Schuller et al., 2014). Therefore, early studies on the molecular basis of clubroot were focused on hormonal control, but now have expanded further into analyzing whole genome defense responses.
Genetic and molecular studies to increase resistance to P. brassicae have traditionally focused on finding resistance genes or resistance loci that can be introgressed into cultivars by classical breeding approaches using diverse marker techniques. Two such genes (Crr1 and Crr2) were first identified in Brassica rapa L. in an F2 segregating population, showing that they acted together to confer resistance (Suwabe et al., 2003). Crr1 was subsequently characterized as consisting of two genes, Crr1a corresponding to a TIR-NBS-LRR protein and a second locus (Crr1b) that was not characterized, but which seems to have a minor effect on resistance (Hatakeyama et al., 2013). Another resistance gene (CRa) was originally introgressed from turnip (B. rapa ssp. rapifera) and mapped by linkage analysis in Chinese cabbage (B. rapa ssp. pekinensis) (Matsumoto et al., 1998), and then fully characterized as a TIR-NBS-LRR protein (Ueno et al., 2012). The Rcr1 locus, which was linked to clubroot resistance in B. rapa ssp. chinensis, was associated with differential expression of defense mechanisms in genotypes with or without the allele (Chu et al., 2014). The region encompassing this marker includes at least 35 genes, of which two (R-genes) have the larger number of polymorphisms associated with differences in resistance (Yu et al., 2016)
Mapping in B. napus has shown different sources of resistance located on diverse chromosome locations depending on the specific genotype/cultivar-isolate interaction. Such resistance can come from a major gene (Manzanares-Dauleux et al., 2000; Hasan and Rahman, 2016; Zhang et al., 2016a), or from several QTLs (Jung, 2008; Li et al., 2016) that do not necessarily match resistance genes found in progenitor genomes or in other genomes from the family. This fact is probably related to genetic rearrangements which are common during polyploidization events (Szadkowski et al., 2011; Xu et al., 2012). For example, in Brassica oleracea L., some QTLs were found to co-localize with previously identified CR loci for that species, but are not analogous to CR loci from B. rapa (Lee et al., 2016). In Arabidopsis thaliana (L.) Heynh, several QTLs related to resistance do not map to previously characterized monogenic resistance in a different genotype (Jubault et al., 2008), and QTL regions encompass R-genes, as well as genes related to the cell wall, antimicrobial compounds, the oxidative burst, auxin metabolism, cell expansion and root architecture. These data support a mix of quantitative and qualitative resistance that is highly dependent on the host-pathogen interaction, and highlights the need for exploring key molecular regulatory mechanisms beyond R-genes.
Transcriptome studies in the clubroot pathosystem have allowed the generation of a molecular landscape in susceptible and resistant interactions, which can serve as a scaffold for understanding the genomic basis of resistance. Studies analyzing host transcriptomic changes first used microarrays, but later switched to RNA-seq analysis, which allows discovery of novel transcripts and has a larger dynamic range for transcript level detection. Microarray analysis of the clubroot susceptible A. thaliana ecotype Col-0 challenged with P. brassicae showed transcriptional changes in cell growth, carbohydrate and defense mechanisms, and modulation of genes related to auxin, cytokinin and brassinosteroid metabolism (Siemens et al., 2006; Agarwal et al., 2011; Schuller et al., 2014). An alternate approach using microarrays was undertaken by using A. thaliana Bur-0, which showed partially resistant and susceptible reactions, respectively, to two different P. brassicae isolates (Jubault et al., 2013). This study indicated that the partially resistant interaction delays pathogen control of host metabolism and induces stronger or earlier defense responses.
The incremental shift towards use of RNA-seq has validated some of the results from microarrays, but also increased the level of information in genome-scale studies. One of the early studies with RNA-seq used B. rapa ssp. chinensis lines segregating for the clubroot resistance gene Rcr1 (Chu et al., 2014). Resistant lines carrying the R-gene upregulated jasmonate and ethylene metabolism, reinforced cell walls through callose deposition, and used indole-containing metabolites as defense mechanisms. More recent studies with B. rapa ssp. pekinensis highlight how resistant interactions involve receptors of pathogen recognition, pathogenesis-related proteins, salicylic acid signalling and secondary metabolism as key elements of defense (Chen et al., 2016a; Jia et al., 2017). In Arabidopsis, responses at 24 and 48 h after inoculation (period of primary infection by P. brassicae) of the susceptible ecotype Col-0 included secondary metabolite synthesis (flavonoids), lignification, hormone signalling (auxins and cytokinins) and receptor-like kinase activation (Zhao et al., 2017). Surveys of the same interaction at 17, 20 and 24 days post-inoculation (dpi) in shoots and roots showed downregulation of primary metabolism (photosynthesis, carbohydrate metabolism) and cell wall modification, while an upregulation of secondary metabolism-related genes (glucosinolates, camalexin, phytoalexins) was detected mainly in shoots (Irani et al., 2018). Roots showed features of disease development through upregulation of cell wall modification enzymes and auxin response genes, resulting in cell growth, division and expansion.
As RNA-seq becomes an important tool in understanding the molecular mechanisms of defense against P. brassicae, specific interactions between particular pathotypes and host genotypes can now be studied using this technology. Addressing specific pathotype-host genotype interactions is necessary to dissect constitutive and induced defenses and how these mechanisms may be altered when new interactions arise. Furthermore, a focus on understanding the responses of B. napus to specific pathotypes would enable a direct assessment of important genes in this host, rather than looking for homologous genes in close relatives. Selecting candidate genes in B. napus for functional validation, or alteration of resistance through gene editing, can shorten the time required for the development of new genotypes/cultivars with increased host resistance.
In the current study, we followed disease progression through a time-course in two B. napus genotypes having a contrasting interaction with pathotype 5X of P. brassicae, as classified on the Canadian Clubroot Differential set (Strelkov et al., 2018). This pathotype can overcome the resistance in most CR canola cultivars and is one of numerous novel pathotypes of the clubroot pathogen identified in western Canada in recent years (Strelkov et al., 2016a). Our study shows that the host (especially the moderately resistant genotype) relies heavily on salicylic acid-mediated immunity and that several defense mechanisms mimic defense mechanisms of plants when faced with infection by nematodes. Our research presents a complete description of the molecular responses in this specific interaction, outlining previously reported key defense genes, and shedding light on candidate genes that behave as defense and susceptibility factors.
Materials and Methods
Plant Material, Inoculations and Harvesting
Seeds of the rutabaga (B. napus subsp. napobrassica ‘Laurentian’) and the oilseed rape (B. napus var. napus ‘Brutor’) were germinated on moistened filter paper in Petri dishes for 1 week prior to inoculation. ‘Laurentian’ is characterized as resistant (index of disease, ID <50%) to pathotype 5X, while ‘Brutor’ is susceptible (ID >50%) to this pathotype (Strelkov et al., 2018).
Inoculum of pathotype 5X was prepared from galls collected initially from the canola (B. napus var. napus) ‘L135C' near Westlock, Alberta, Canada, in 2013 (Strelkov et al., 2016a). One hundred grams of galled root tissue were ground in 1 L of sterile distilled (sd)-water in a blender and filtered through eight layers of cheesecloth. The resting spore concentration in the filtrate was estimated using a haemocytometer and adjusted to 1 × 107 spores/ml with sd-water.
One-week old seedlings used as controls were transferred directly to pots filled with water-saturated Sunshine LA4 potting mix (SunGro Horticulture, Vancouver, BC, Canada), while plants used for inoculation were placed in a Petri dish with the inoculum suspension for 15 s before transferring them to the pots. An additional 1 ml of the spore inoculum was added directly to the potting mix for the inoculated plants to ensure strong disease pressure.
The plants were placed in insect cages (47.5 cm × 47.5 cm × 93.0 cm) to avoid potential arthropod infestations that might interfere with the results, and maintained in a greenhouse under long day conditions (16 h) and an average temperature of 22°C. Six biological replicates were collected at 7, 14 and 21 days after inoculation (dai) for each of the control vs. inoculated treatments and for both genotypes studied. Biological replicates comprised 12 plants, which were pooled at the time of harvest for each time-point and treatment. Five replicates were harvested for downstream expression analysis and individual plants from the last replicate were used to collect samples for microscopy. The collected roots were washed thoroughly under running water, dried on a paper towel, placed in 50 ml Falcon tubes (Fisher Scientific, Hampton, NH, U.S.), and transferred immediately to liquid nitrogen. Root sections (0.5-cm long, up to 4 cm from the root base) were collected for microscopy analysis as described below.
Microscopy
For the microscopy analysis, root tip sections (0.5 cm) were placed in fixative solution (10% neutral buffered formalin) for a minimum of 2 days. Samples were then placed in plastic cassettes and dehydrated in a 50–70–90–100% ethanol series in a Leica TP1020 tissue processor (Leica, Nussloch, Germany). Tissues were embedded in paraffin blocks using a TISSUE TEK II embedding center (Sakura, Torrance, CA, USA). Eight to 12 µm sections were obtained with a RM2125 microtome (Leica, Nussloch, Germany), and dried on microscope slides overnight at 37°C.
Staining was performed with Hematoxylin–Eosin according to the following procedure: slides were dewaxed in toluene for 5 min twice and then washed with 100% ethanol twice for 2-min, followed by a 2-min wash in each of three ethanol solutions (90, 70, 50%). The slides were washed in distilled water for 1–2 min, Hematoxylin Gill III solution (Leica, Nussloch, Germany) was added for 2 min, and then rinsed with distilled water. Next, the slides were placed under running tap water for 15 min, washed for 2 min with 70% ethanol, and Eosin (Leica) was added for exactly 30 s. The Eosin was washed with 100% ethanol twice for 2 min, followed by two cleaning steps using toluene for 2 min. DPX mounting medium (EMS, Hatfield, PA. U.S.A.) was added to the slides to fix the preparations and inhibit stain fading, and the sections were covered with a coverslip and incubated at 37°C overnight. The stained sections were examined under an Optika B-290TB light microscope (Ponteranica, Italy), and images were captured with the Optika Vision Lite v2.13 digital USB camera system.
RNA Extraction
Frozen roots were placed in 15 ml Falcon tubes and ground using 5–7 mm-diameter steel beads at 1,100 to 1,200 rpm (3–4 rounds of 1 min each) in a Genogrinder (SPEX SamplePrep, Metuchem, NJ, U.S.A), followed by additional grinding in a mortar with a pestle in the presence of liquid nitrogen. Total RNA was extracted with a combination of TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, U.S.A.) and the RNeasy kit (QIAGEN, Venlo, Netherlands) to increase the quantity and purity of the RNA. DNase treatment was performed on-column with the RNase-Free DNase set (QIAGEN, Venlo, Netherlands). The quantity of RNA was assessed with a Nanodrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, U.S.A.) and quality was checked in a 2200 TapeStation (Agilent, Santa Clara, CA, U.S.A.).
RNA-Seq
Three micrograms of RNA per pooled sample corresponding to each biological replicate (see above), with RNA Integrity Numbers (RIN) > 9, were sent for sequencing at Oklahoma State Genomics (Stillwater, OK, U.S.A.). Libraries were prepared with a TruSeq Stranded mRNA library preparation kit for 48 samples (Illumina, San Diego, CA, U.S.A) and samples were sequenced using TG Nextseq 500/550 High Output Kit v2 (75 cycles), in a NextSeq500 equipment (Illumina, San Diego, CA, U.S.A). The raw data reads were filtered using the bcl2fastq conversion software provided by Illumina. This software uses read indexes to assign each read to the corresponding sample and then removes any adapters, barcodes, or primers using a trimmomatic-like algorithm. Any reads with an average Q-score <30 were removed and trimming of read ends was performed for bins of 15 bp that had a Q-score <30. Sequencing reads were deposited in the NCBI Sequence Read Archive (SRA) under accession number: PRJNA597078.
The reads from all treatments were mapped to the B. napus reference genome (AST_PRJEB5043_v1) using the genome and gene annotation deposited in Ensembl (Chalhoub et al., 2014; Ensembl, 2018; Kersey et al., 2018). Genome indexes were generated with Bowtie2 v2.3.3.1 (Langmead and Salzberg, 2012) and transcriptome indexes with TopHat v2.1.1 (Trapnell et al., 2009). Reads mapped with TopHat were used along gene models (obtained from the Ensembl website) as input for Cufflinks v.2.2.1 to generate all potential transcripts in each sample. Cufflinks was run using the GTF-guide option, multi-read correction and fragment bias correction (multi-read correction and fragment bias correction were also used in the subprograms downstream). GTF files generated from all treatments and replicates were used to generate a consensus merged transcript gtf file with Cuffmerge. This consensus file was employed as a reference for re-mapping reads for quantification. Cuffquant was run with the merged gtf file to create ready-to-use files with quantified transcripts, to perform the necessary comparisons of differential expression among treatments. Cuffdiff was used to compare cuffquant files, contrasting controls with inoculated plants at each time-point and for each genotype separately. Expression levels were quantified as Fragments Per Kilobase of transcript per Million mapped reads (FPKM) and significant differential expression for each evaluated transcript was assessed after multiple comparison testing with the Benjamini-Hochberg correction (Benjamini and Hochberg, 1995).
The number of mapped reads from all treatments was compared using a one way ANOVA implemented in R. FPKM values of all genes from treated and inoculated replicates were log-transformed to build PCA plots using the R-package ggfortify (https://github.com/sinhrks/ggfortify) (Son et al., 2018).
cDNA Synthesis and Quantitative Reverse Transcription PCR (qRT-PCR)
Five hundred nanograms of RNA per sample were subjected to DNAse treatment using RNase-free-DNase I (Thermo Fisher Scientific, Waltham, MA, U.S.A.). cDNA was synthesized using the RevertAid H Minus Reverse transcriptase with oligo dT (18), following manufacturers' specifications (Thermo Fisher Scientific). To test for residual DNA contamination, a 1:100 dilution of cDNA was used to perform PCR with primers corresponding to a constitutively expressed Clathrin Adaptor Complex (CAC) gene (Supplementary Table S1), producing 125 bp and 288 bp bands on cDNA and DNA, respectively. PCR analysis was conducted with 1× Buffer + KCl, 2.5 mM MgCl2, 0.2 mM of each dNTP, 0.2 µM or each primer, and 1 unit of high fidelity Taq polymerase. Cycling conditions were 94°C for 3 min, followed by 35 cycles of 94°C for 30 s, 60°C for 30 s, and 72°C for 1 min, finishing with an extension step of 72°C for 10 min, and a 4°C step until samples were recovered from the thermocycler.
To validate the differential expression detected by RNA-seq analysis, 15 primer pairs from genes showing significant differential expression in at least one time-point or host genotype were evaluated for efficiency (range of 90 to 110%) using a 5-point dilution series (1:4 serial dilutions) of a mix of all cDNA samples. Ten target primers within the expected efficiency range were used for final RNA-seq validation (Supplementary Table S1). Five reference gene primers with valid efficiencies were tested for stability across all samples using BestKeeper (Pfaffl et al., 2004), and the three genes showing the best correlation and lower standard deviation were chosen for relative quantification (Supplementary Table S1). Primers were designed using Primer3 (NIH, 2012; Untergasser et al., 2012) with the following parameters: Na+ concentration 50 mM, Mg++ concentration 3 mM, dNTPs concentration 0.8 mM, oligo length: min 18, opt 22, max 30 bp, melting temperature: min 60, opt 62, max 64°C, GC content: min 35, opt 50, max 65%, amplicon length: min 70, opt 100, max 250 bp.
Relative quantification of expression was performed using the 2−ΔΔCT method (Livak and Schmittgen, 2001). All qRT-PCR experiments were performed with four biological replicates and three technical replicates per biological replicate. The experiments were run on a ViiA7 Real-Time PCR system (Applied Biosystems-Life Technologies, Carslbad, CA, USA). Reactions were performed in 10 µl with 5 µl of an in-house made SYBR-green reagent mix, 2.5 µl of the pair of mixed primers (3.2 µM), and 2.5 µl of a 1:100 dilution of cDNA. Cycling conditions were 95°C for 2 min followed by 40 cycles of 95°C for 30 s and 60°C for 1 min. A melting curve stage was added to verify product specificity: 95°C for 15 s, 60°C for 1 min, and 95°C for 15 s.
Log2-fold change values from the qRT-PCR analysis were compared with the RNA-seq values of the respective genes to validate the RNA-seq data. Correlation coefficients were calculated using the CORREL function from Microsoft Excel (Microsoft, Redmon, WA, U.S.A.) for each day and host genotype using the 10 target genes.
Bioinformatic Analyses
Sequences from all transcripts were obtained using a Perl script from transdecoder (https://github.com/TransDecoder/TransDecoder/releases), which uses the coordinates from the merged GTF transcript file of all treatments to extract the sequences from the reference genome (https://plants.ensembl.org/Brassica_napus/Info/Index). Annotation of transcripts was performed with BLASTx against the A. thaliana protein database release 10 (TAIR10) using a maximum of 20 hits and an e-value of 1e−10. The XML BLAST output was loaded to BLAST2GO (Conesa et al., 2005; Conesa and Stefan, 2008) to obtain a consensus annotation for each transcript derived from the 20 top hits. The top A. thaliana ID (AtID) was also stored in the database for further analyses.
To find general patterns of gene expression response, we performed functional categorization of genes that were significantly differentially expressed with log2-fold changes >1 (upregulated) and ≤1 (downregulated) for each time-point and genotype. The respective AtIDs were used as input for functional enrichment analysis using DAVID (Database for Annotation, Visualization and Integrated Discovery) (Huang et al., 2007; Huang et al., 2009). DAVID was set to obtain gene ontology functional enrichment analysis (biological process level 3), protein domain enrichment analysis, and KEGG pathway enrichment analysis. All transcript AtIDs from the RNA-seq study were used as the background population and a minimum count of five hits was established to find enriched categories bearing an EASE score <0.05. Enriched categories used to create heatmaps were selected after a Benjamini correction for multiple comparisons (p-value <0.05). Heatmaps were built with the matrix visualization software Morpheus (https://software.broadinstitute.org/morpheus/) using the corrected p-values for enriched gene ontology (GO) biological process ‘level 3' categories, Interpro protein enriched domains, and KEGG enriched pathways.
Log2-fold expression values of genes from selected enriched categories from DAVID were used to observe changes through the time-course for both genotypes using Multi-Experiment Viewer (MeV4.9) (Saeed et al., 2003). Hierarchical clustering for genes was applied using Pearson correlation and average linkage. Sub-clusters were generated when necessary to display co-regulated genes.
To find genes with similar or divergent regulation between treatments we used the online tool jvenn (Bardou et al., 2014).
Results and Discussion
Disease Development
We compared the anatomical changes that occurred in roots of ‘Laurentian’ (resistant) and ‘Brutor’ (susceptible) following inoculation with P. brassicae. The differential response of these hosts was first noted when they were inoculated with resistance-breaking field isolates of the clubroot pathogen recovered from canola crops in Alberta, Canada, in 2013 (Strelkov et al., 2016a). These isolates are classified as pathotype 5 on the differentials of Williams (1966) or as pathotype 5X on the CCD Set (Strelkov et al., 2018). This pathotype is now one of several ‘new’ resistance-breaking pathotypes of concern for canola farmers in western Canada (Strelkov et al., 2016a).
On the susceptible ‘Brutor’, root galls were visible as early as 14 dai, while no galls were observed on the resistant host ‘Laurentian’ at this time-point. At 21 dai, all ‘Brutor’ plants had developed large galls, while small galls were observed on half of the ‘Laurentian’ seedlings (Figure 1). These results indicate that clubroot develops more rapidly and severely on the susceptible genotype. Nonetheless, while there was a differential response between the host genotypes to pathotype 5X, clubroot severity on ‘Laurentian’ was greater than previously reported (Strelkov et al., 2018). Hence, while ‘Laurentian’ showed some resistance to pathotype 5X, it was not immune.
Figure 1 Root gall formation following inoculation with pathotype 5× of Plasmodiophora brassicae. No or only incipient galls were observed on the moderately resistant host ‘Laurentian’ (A) 21 days after inoculation (dai), while larger and more developed galls had formed on the susceptible host ‘Brutor’ (B) at the same time-point.
Roots of both hosts were stained with Haematoxylin-Eosin. At 14 dai, ‘Laurentian’ showed few or no signs of infection in the examined cells, while ‘Brutor’ presented what appeared to be young secondary plasmodia. However, at 21 dai, different stages of plasmodial development could be observed inside the cortical cells of both genotypes (Figure 2).
Figure 2 Plasmodiophora brassicae infection on roots 21 dai. A cross-section through a large root gall of Brassica napus ‘Brutor’ (A), showing the presence of secondary plasmodia, which are visible in more detail under higher magnification (B, C). Young and advanced secondary plasmodia were also visible in other sections of the susceptible ‘Brutor’ (D–F), as well as in the resistant B. napus ‘Laurentian’ (G–I). Mature secondary plasmodia (MSP), young secondary plasmodia (YSP), resting spores (RS).
RNA-Seq Data Analysis and Validation
RNA-seq was performed to detect contrasting gene expression patterns between the resistant and susceptible genotypes over a time-course that covered secondary infection by P. brassicae. Sequencing of three biological replicates per each inoculated or control treatment for each of ‘Laurentian’ and ‘Brutor’ over three time-points (7, 14 and 21 dai) resulted in over 30 million reads on average per treatment. While ANOVA showed no significant differences in the amount of mapped reads among all treatments (p = 0.645), there was a trend showing that the percentage of mapped reads was always lower for the inoculated treatments than for the control treatments (Supplementary Table S2). Also, the percentage of mapped reads for the inoculated condition was always lower for ‘Brutor’ than for ‘Laurentian’, and decreased drastically at 21 dai for both hosts, although to a larger extent for ‘Brutor'. These trends are consistent with an increasing number of P. brassicae cells colonizing the root tissues of both genotypes (part of the sequenced reads come from P. brassicae and therefore are not mapped to B. napus).
The reads were mapped to the B. napus reference genome and corresponding gene models, to find expression patterns in existing and novel genes. Fragments Per Kilobase per Million mapped reads (FPKM) values of 105,538 predicted transcripts were used for PCA to test the consistency of the biological replicates. Principal components analysis in each genotype showed the expected clustering of the three replicates from each treatment in the two host genotypes; only one control replicate of ‘Brutor’ at 21 days seemed closer to the replicates from the 14 dai control samples from the same genotype. A larger biological variability among samples can be due to uncontrolled environmental conditions or inherent variations of development among sampled plants. Nevertheless, there was always clear separation between inoculated and control samples in each host and time-point, demonstrating that the pattern of expression is influenced by the treatment and comparisons we established (Figures 3A, B). A second PCA analysis was performed to compare control samples from the two genotypes, separated from inoculated samples from the two genotypes. The analysis revealed that control samples were well separated in the first component between genotypes and that the samples from the same day in the two genotypes corresponded in the second component (Supplementary Figure S1). The inoculated samples were also separated similarly to the control samples in the first component, but there was separation between the samples of the two genotypes for the corresponding day in the second component. This could indicate some matching regulation between the two genotypes in different time-points, as discussed below.
Figure 3 Principal Component Analysis (PCA) plots for RNA-seq data in each genotype. The three biological replicates corresponding to each treatment (control vs. inoculation) in the three time points tested (7, 14 and 21 dai) show consistent clustering for both host genotypes, Brassica napus ‘Laurentian’ and ‘Brutor’. (A) Laurentian PCA plot. (B) Brutor PCA plot. ID codes: LA, Laurentian; BR, Brutor; I, Inoculated; C, Control; number after I or C stands for number of days after inoculation (dai).
Of 105,538 predicted transcripts, 85,956 (81.44%) were annotated using the TAIR 10 database. After correcting for multiple comparisons, 5,538, 2,747 and 4,349 genes were regulated at 7, 14 and 21 dai, respectively, for ‘Laurentian’, and 3,834, 3,415 and 6,570 for ‘Brutor’. The number of regulated genes for ‘Laurentian’ was higher than for ‘Brutor' at 7 dai, and there were more upregulated than downregulated genes for both hosts at this time-point. At 14 and 21 dai, the number of downregulated genes was greater than the number of upregulated genes in both hosts; ‘Brutor’ showed the largest number of regulated genes at 21 dai, with most being downregulated (Table 1). In the clubroot susceptible A. thaliana ecotype ‘Columbia' (Col), an increased number of genes in shoots and roots were regulated at later time-points of infection (Irani et al., 2018), and in B. rapa ssp. pekinensis, a larger number of genes was regulated in the susceptible interaction when compared with the resistant interaction at 30 dai (Jia et al., 2017). To increase the stringency of downstream analyses, we generated an arbitrary cut-off for significant differentially expressed genes at log2-fold changes >1 and ≤1 for up and downregulated genes, respectively (Table 1).
Table 1 Significant differentially expressed genes in both genotypes throughout the time course, without cut-off (up or down) or with log2-fold change >1 or ≤1.
To validate the trends in DEGs detected via the RNA-seq analysis, we evaluated 10 target genes across both genotypes by qRT-PCR analysis. The genes selected were differentially expressed and significant in at least one host, in at least one time-point. Comparison of log2-fold changes between the RNA-seq and qRT-PCR data showed correlations of >0.8 for the three time-points in the two genotypes (Supplementary Table S3).
Major Transcriptional Changes in the Resistant and Susceptible Interactions
Functional enrichment analysis with DAVID using significant differentially expressed genes (log2-fold changes >1 or ≤1) across all treatments identified 119 ‘level 3’ GO biological process enriched functional categories (Figure 4), 73 enriched Interpro domains and 28 enriched KEGG pathways (Supplementary Figure S2).
Figure 4 DAVID GO functional category enrichment analysis. Gene ontology (GO) enriched biological process ‘level 3' functional categories are shown across three time-points (7, 14 and 21 days) for upregulated (up) and downregulated (down) genes in the resistant and susceptible Brassica napus genotypes, ‘Laurentian' (LA) and ‘Brutor' (BR), respectively. Significantly enriched categories are color-coded (red more significant) according to the EASE score, which constitutes a modified Fisher's exact test (Huang et al., 2009). Most significant EASE value = 2.94E−27, least significant EASE value = 4.95E−02. Non-enriched bins are color coded gray.
General Stress Response
A wide scale stress response is expected upon pathogen challenge. Functional categories where all or most bins were regulated through the time-course for the two host genotypes included response to oxygen-containing compounds, response to inorganic/organic substances, response to external biotic stimulus, response to hormones, response to acid chemicals, response to osmotic stress and water deprivation, response to oxidative stress/oxidation-reduction, and response to bacterium (Figure 4). The enrichment of the response to hormone category supports the central role that hormones play in this pathosystem, where growth hormones like auxins and cytokinins have been directly involved in the hypertrophy and hyperplasia of host tissues invaded by P. brassicae (Ludwig-Müller, 2009; Ludwig-Müller et al., 2009). Drought stress enrichment likely reflects the inhibition of water uptake and transport during secondary infection, due to disruption of the roots and vascular system (Ludwig-Müller, 2009). Finally, the oxidative response is part of a defense mechanism that includes the activation of oxidative burst-related enzymes that act as signalling molecules, but can also have direct antimicrobial effects or result in localized cell death. This behavior is supported by several studies of the P. brassicae–Brassicaceae interaction (Agarwal et al., 2011; Jubault et al., 2013; Zhang et al., 2016b; Zhao et al., 2017).
Immunity in Response to P. brassicae Is Regulated Extensively by Salicylic Acid
Innate immune responses in plants can be activated in roots by pathogen-associated molecular patterns (PAMPs) leading to PAMP-triggered immunity (PTI), or by effectors leading to effector-triggered immunity (ETI) (Miller et al., 2017). Molecule recognition usually leads to an oxidative burst (via the production of reactive oxygen species (ROS)), where calcium takes part in a signal transduction cascade to modulate transcription factors that in turn activate defense responses (Stael et al., 2015). In our study, a clear distinction was observed in the innate immune response functional category. This category was enriched for upregulated genes in ‘Laurentian’ at 7, 14 and 21 dai (Figure 4). In contrast, ‘Brutor’ showed enrichment for this category only for upregulated genes at 7 dai and for downregulated genes at 21 dai. While the expression level of many genes in this functional category followed the pattern of enrichment (Figure 5A), other immune response genes where mostly upregulated (Figure 5B) or downregulated (Figure 5C) for both genotypes throughout the time course.
Figure 5 Immunity-related gene regulation. Genes from this functional category clustered according to expression patterns. Heatmaps were built with MeV with log2-fold expression limits of 4 (red) to −4 (green). Hierarchical clustering of gene expression patterns was performed using Pearson correlation. Clusters (A–C) represent different trends in gene expression co-regulation.
Many of the genes that displayed increased transcriptional levels in ‘Laurentian’ at 7, 14 and 21 dai, and in ‘Brutor’ at 7 dai were key receptors/signalling genes or were related to salicylic acid (SA)-mediated responses (Figure 5A). Other immunity-related genes displayed more heterogeneous patterns of regulation (Supplementary Figure S3).
Kinase regulators are key in signal transduction of pathogen signals. Among these, RBK1 (ROP-binding kinase 1) showed sustained expression in the resistant host ‘Laurentian’ (Figure 5A). RBK1 is phosphorylated by the mitogen-activated protein kinase 1 (MPK1), which negatively influences auxin-dependent cell expansion (Enders et al., 2017). Therefore, RBK1 may be a key component in downstream regulation of auxin-dependent responses.
Downstream of signal transduction, the functionally redundant transcription factors involved in SA signalling, Calmodulin-Binding Protein 60g (CBP60g) and Systemic-Acquired Resistance Deficient 1 (SARD1) (Wang et al., 2011a), were also regulated in the current study (Figure 5A). Two transcripts corresponding to CBP60g had their highest upregulation at 14 dai in ‘Laurentian’, while a transcript matching SARD1 was upregulated throughout the time course (Figure 5A). For the susceptible genotype both genes were upregulated at 7 dai and downregulated at 14 and 21 dai. The susceptibility of A. thaliana to Pseudomonas syringae was enhanced when these two genes were mutated, with a concomitant reduction of both, salicylic acid (SA), and SA-mediated marker genes like PR1 (Wang et al., 2011a). Two genes involved in regulation and transport of SA, isochorismate synthase (ICS) and enhanced disease susceptibility 5 (EDS5), act downstream of CPB60g/SARD1 (Zhang et al., 2010; Wang et al., 2011a). Isochorismate synthase is translocated to the chloroplasts and is involved in the transformation of chorismite into isochorismate, a precursor of SA, while EDS5 allows the translocation of the SA produced in the chloroplast into the cytoplasm. One transcript corresponding to gene ICS2 (BnaC08g18420D), which is redundant in function with ICS1, was highly expressed across the three time-points in ‘Laurentian’ (Supplementary Table S4), while expression of EDS5 peaked at 21 dai in this genotype (Supplementary Figure S3B). Once SA accumulates in the cytoplasm, it is regulated via different mechanisms. For example, the pool of SA can be increased when cytoplasmic methyl esterases transform MeSA into SA (Dempsey et al., 2011). Transcripts matching methyl esterase 9 (MES9) showed high expression in both genotypes, with a stronger regulation in ‘Laurentian’ (Figures 5A, B). MES9 from B. oleracea was also upregulated in both resistant and susceptible genotypes upon P. brassicae infection, but with a stronger regulation in the susceptible genotype (Manoharan et al., 2016). Another gene, PBS3 (AvrPphH susceptible 3), promotes accumulation of SA by suppressing proteins that conjugate SA to amino acids (Lee et al., 2007; Nobuta et al., 2007; Dempsey et al., 2011). Three transcripts of PBS3 were predominately upregulated throughout the time course in ‘Laurentian’ (Figure 5A). This gene confers resistance to the biotrophic pathogen P. syringae and acts synergistically with NPR1 upstream of SA to mediate defense responses (Lee et al., 2007; Nobuta et al., 2007).
Regulation of responses to different levels of SA are mediated by an array of different genes. For example the EDS1 (enhanced disease susceptibility 1)-PAD4 (phytoalexin deficient 4)/SAG101 (senescence-associated gene 101) complex (Wagner et al., 2013) is involved in a positive feedback loop for SA accumulation that includes ICS expression (Seyfferth and Tsuda, 2014). While our results indicated regulation of transcripts corresponding to these genes, their transcriptional levels did not produce a clear trend for either host genotype or time-point (Supplementary Figure S3C and Supplementary Table S4). Another family of genes, the non-expressor of PR (NPR) genes, are involved in regulating SA levels. NPR1 is the main transducer of SA signals, interacting as a co-factor of transcription factors to alter defense (Cao et al., 1997; Dong, 2004). NPR1interacts specifically with TGA factors to bind the promoter of the important defense marker PR1 (Chen et al., 2016b), to modulate its activation (Weigel et al., 2005). While our study did not find regulation of the transcript corresponding to NPR1, PR1 was highly upregulated throughout the time course in both host genotypes (Figure 5B). A second transcript of the NPR family, NPR3, was significantly upregulated at 14 and 21 dai in ‘Laurentian’ and at 7 dai in ‘Brutor’ (Figure 5B). NPR3 and its paralog NPR4 serve as adaptors for a CUL3 ligase, targeting NPR1 for degradation depending on SA concentration; with low concentrations promoting an NPR1-NPR4 association and with high SA concentrations promoting an NPR1-NPR3 association (Fu et al., 2012; Moreau et al., 2012). Another mechanism to regulate SA through degradation is provided by the nitrogen limitation adaptation (NLA) gene. This gene acts as a negative regulator of SA accumulation upon P. syringae infection, and seems to be independent of the ICS1-mediated SA accumulation (Yaeno and Iba, 2008). It was argued that NLA codes for a ring-type ubiquitin E3 ligase, and mutations in this gene result in excessive levels of SA (Yaeno and Iba, 2008), suggesting that NLA is involved in SA-mediated immune response regulation via ubiquitination. Transcripts of NLA were downregulated in ‘Brutor’ at 14 and 21 dai and upregulated in ‘Laurentian’ throughout the time-course (Figure 5A). Keeping SA levels consistent (not under- or over-regulated) may be necessary for the basal defense response of the resistant ‘Laurentian’, hence the observed upregulation of NLA throughout the time course by this host. Supporting this view, silencing of NLA via microRNA activation increased susceptibility of A. thaliana to a cyst nematode (Hewezi et al., 2016). One more gene that could have a similar function to NLA was the carbon/nitrogen insensitive 1 gene (CNI1). This gene was upregulated in ‘Laurentian’ throughout the time course and downregulated in ‘Brutor’ at 14 and 21 dai (Figure 5A). CNI1 is also known as Arabidopsis toxicos en levadura 31 (ATL31), a RING-type ubiquitin ligase that is upregulated by bacterial elicitors and increases plant resistance when overexpressed (Maekawa et al., 2012). We speculate that the resulting SA-immune response inferred from our differential expression analysis (especially in the resistant genotype) is potentially accompanied by increasing SA levels, which can cause the activation of genes that fine tune the SA-mediated plant defense.
Evidence of the regulation of immunity through SA in our study is consistent with a typical biotrophic interaction, with longer lasting effects in the resistant but not in the susceptible host. In A. thaliana, the partially resistant ecotype Bur-0 activated SA signalling through secondary pathogen infection, while the susceptible Col-0 did not show SA-mediated responses (Lemarie et al., 2015). The exogenous application of SA reduced clubroot symptom development on B. oleracea and A. thaliana (Lovelock et al., 2013; Lemarie et al., 2015), highlighting the importance of this hormone for defense against P. brassicae. A model of SA-mediated immunity responses is presented in Figure 6. While many of the genes found in our study are regulated in the same manner in both hosts, key genes are mainly upregulated in the resistant genotype, which supports a stronger SA-mediated immune response in ‘Laurentian’.
Figure 6 Model of SA-mediated immunity. SA-related genes (bold-italics-underlined) were regulated differentially in the interaction between Plasmodiophora brassicae and the Brassica napus hosts ‘Laurentian' and ‘Brutor', as indicated by the left and right arrows, respectively, next to each gene (red: upregulated, green: downregulated, yellow: mixed/undetermined regulation; expression trend of the arrows is based on all time-points and associated transcripts of the specific gene with log2-fold changes >1 or ≤1). Detection of PAMPs and effectors by RLKs and R-genes triggers a change in calcium status and an oxidative burst leading to signal transduction. Calcium sensor proteins (CaM) can bind transcription factors such as CBP60g/SARD1 or CBP60a, which positively or negatively regulate expression of genes like ICS. ICS is translocated to the chloroplast to transform chorismate into isochorismate, a SA precursor. SA is then translocated into the cytoplasm by the EDS5 transporter, where additional SA can be synthesized from methyl salicylate thanks to the action of MES9. Under low SA levels, NPR1 binds NPR4, but as SA increases NPR4 binds SA and dissociates from NPR1, allowing NPR1 to fine-tune the response of defense genes (like PR1) via its interaction with transcription factors (TGA). If SA levels keep increasing, NPR3 binds SA and associates with NPR1, promoting NPR1 degradation via the proteasome. At the same time, feedback loops between associated homomeric and heteromeric complexes between PAD4/SAG101 and EDS1, and the gene PBS3, regulate SA levels. Finally, NLA/CNI1 help to control the levels of SA tagging unknown genes for degradation via the proteasome. Abbreviations: PAMPs (pathogen-associated molecular patterns), ROS (reactive oxygen species), RLKs (receptor like kinases), R-genes (resistance genes), CaM (calcium-modulated protein/calmodulin), CBP60g (calmodulin-binding protein 60g), SARD1 (systemic-acquired resistance deficient 1), CBP60a (calmodulin-binding protein 60a), ICS1 (isochorimate synthase 1), ICS2 (isochorismate synthase 2), IPL (isochorismate pyruvate lyase), EDS5 (enhanced disease susceptibility 5), MES9 (methyl esterase 9), NPR1 (non-expressor of pathogenesis related genes 1), NPR3 (non-expressor of pathogenesis related genes 3), NPR4 (non-expressor of pathogenesis related genes 4), NLA (nitrogen limitation adaptation), CNI1 (carbon/nitrogen insensitive 1), TGA (TGA-binding bZIP transcription factors), PR1 (pathogenesis-related protein 1), EDS1 (enhanced disease susceptibility 1), SAG101 (senescence-associated gene 101), PAD4 (phytoalexin deficient 4), PBS3 (AvrPphB susceptible 3).
SA and SAR Are Antagonistic to JA/ET in Both Genotypes
After establishing the mechanisms of SA-mediated immunity, we analyzed the role of other DEGs that were also part of the immune response functional category. Transcription of some of these genes is influenced by SA regulation, while others regulate defenses through crosstalk with other hormones.
We found a strong induction of transcripts for plant natriuretic peptides (PNPs) in both hosts throughout the time course (Figure 5B). PNPs are signalling molecules that respond to biotic and abiotic stresses (Wang et al., 2011b). In response to a P. syringae effector, transcriptional regulation of PNP-A was shown to be dependent on EDS1 and therefore related to SAR (Breitenbach et al., 2014). Furthermore, co-expression analysis showed that A. thaliana PNP-A expresses together with genes that are involved in SAR, and also induced by SA or analogue molecules (Meier et al., 2008). Among these genes, the transcripts corresponding to PR1 PR2, NIMIN1 and WRKY70 were mostly upregulated throughout the time course in both of our hosts (Figures 5A, B, Supplementary Figure S3B). These genes are important components of plant defense responses and WRKY70 was shown to specifically control AtPNP-A (Meier et al., 2008). The transcription factor WRKY70 was also shown to be the key regulator of antagonistic responses between SA and JA (Li et al., 2006). Modulated by NPR1, WRKY70 potentially binds promoters of negative regulators of JA signalling or induces factors that block positive regulators (Li et al., 2006). WRKY70 was activated in the A. thaliana paritally resistant genotype Bur-0 upon P. brassicae infection, in parallel with strong activation of SA-mediated defenses and repression of JA-mediated responses (Jubault et al., 2013). The upregulation of WRKY70 in our study, plus the downregulation of genes like jasmonate resistant 1 (JAR1) (Figure 5C and Supplementary Figure S3C), a gene induced by jasmonate and catalyzer in the generation of active JA (Wasternack and Hause, 2013), support a SA-mediated response which is antagonistic to JA responses.
It is also known that JA acts synergistically with ethylene (ET), usually in response to necrotrophic pathogens (Wasternack and Hause, 2013). In this case, we would also expect an inhibition of ethylene mediated responses in our study. The key marker of JA/ET-mediated defense responses, PDF1.2 (Plant defensin) (Penninckx et al., 1996), was not regulated in either of the genotypes in our study, while a peptide elicitor of PDF1.2 (precursor of peptide 1) that is activated by methyl jasmonate and ethylene (Huffaker et al., 2006) was downregulated (Figure 5C). Additionally, the ethylene response factor 4 (ERF4), which is regulated via JA (McGrath et al., 2005; Yang et al., 2005), and ERF2, which acts as a transcriptional activator of ethylene responsive genes (Fujimoto et al., 2000), were mainly downregulated (Figure 5C and Supplementary Figures S3A, C).
Pathogenesis-Related (PR) Proteins and Susceptibility Factors Contribute to Defense Against Clubroot
PR proteins are generally activated due to induction of SA and SAR. Among these, chitinases and thaumatins have been shown to have a direct effect on pathogen cell walls and membranes, arresting pathogen spread, or producing oligomers that act as signals to initiate reprogramming of defenses (Grover et al., 2012). Most chitinases were mainly upregulated in both genotypes throughout the time course (Figures 5A, B), and chitin-binding domains were also enriched for upregulated ‘Laurentian’ genes 21 dai (Supplementary Figure S2). Upregulation of these genes was detected in CR B. rapa, B. oleracea and B. macrocarpa (Chen et al., 2016a; Zhang et al., 2016b; Wang et al., 2019). The thaumatin-like protein 1 (TLP1) gene, which was upregulated at 7 and 14 dai in both genotypes and downregulated at 21 dai in ‘Brutor’ (Figure 5A), has been shown to be involved in defense against nematodes (Hamamouch et al., 2012), organisms which usually cause symptoms in their host resembling clubroot symptoms (e.g. enlarged cells).
Additional PR genes (PR1, PR2, PR4) were consistently upregulated in both host genotypes (Figure 5B). PR1, PR2 and PR5 are known to be SA-responsive, and PR2 and PR5 have been shown to be induced during secondary infection by P. brassicae of partially resistant A. thaliana (Lemarie et al., 2015).
While upregulation of defense genes is commonly seen as a strategy to deter a pathogen, downregulation of certain genes (susceptibility factors) may also contribute to immunity. In A. thaliana triple mutants, mlo2, mlo6 and mlo12 confer resistance to powdery mildew via the activation of multiple basal defense mechanisms (Kuhn et al., 2017). In our study, the mildew resistance locus 12 gene (MLO12) had one associated transcript with complete downregulation in all treatments (Figure 5C), and a second transcript slightly upregulated at 7 dai in both genotypes (Supplementary Figure S3A), while a transcript matching MLO2 showed a trend similar to the second MLO12 transcript. Another gene, CAD1 (constitutively activated cell death 1) negatively regulates SA-mediated responses, with knockouts of the gene resulting in activation of PR1 (Asada et al., 2011). Downregulation of this gene in our study, mainly in the case of ‘Laurentian' (Figure 5C and Supplementary Figure S3C), supports steady SA-mediated immunity.
Similar Symptoms, Similar Responses: P. brassicae and Nematodes Activate Comparable Pathways
The formation of a syncytium—an enlarged multinucleated cell mass—is usually triggered by parasitic nematode infection of plants (Rodiuc et al., 2014). The formation of syncytia shares some parallels with the root swelling and enlarged cells that are associated with clubroot development. Cell hypertrophy and cell wall modification constitute common processes for both the clubroot interaction and nematode infections (Rodiuc et al., 2014; Chen et al., 2016a). Root-knot nematodes induce galls with giant cells, where mitosis progresses without cytokinesis to generate the syncytium, whereas cell enlargement and accelerated cell division occur during clubroot development. These processes result in disruption of water and nutrient transport through the plant, generating sinks that favor parasite development. In the present study, the category of syncytium formation was enriched for upregulated genes at 21 dai for ‘Laurentian’ and at 14 and 21 dai for ‘Brutor’ (Figure 4). Many of the transcripts within this functional category matched expansins that were upregulated in both hosts, especially at 14 and 21 dai (Figure 7A).
Figure 7 Nematode-like gene regulation responses. Two functional categories in the interaction between Plasmodiophora brassicae and Brassica napus indicated responses similar to nematode infections: (A) Syncytium formation. (B) Response to nematode. Heatmaps were built with MeV with log2-fold expression limits of 4 (red) to −4 (green). Hierarchical clustering of gene expression patterns was performed using Pearson correlation.
‘Laurentian’ showed an enrichment of upregulated genes only at 21 dai for the category response to nematode (Figure 4). Some of the B. napus transcripts that matched the A. thaliana IDs from this functional category suggest a potential effector-related manipulation of the host (Irani et al., 2018). In a susceptible interaction between A. thaliana and P. brassicae, the highest gene expression corresponded to a glucose-6-phosphate/phosphate translocator (GPT2) (Irani et al., 2018). This is the exact same glucose transporter found in our study, which matched three transcripts showing different levels of upregulation in both hosts throughout the time course (Figure 7B). Amino acids are also important for pathogen development, and an amino acid permease (AAP1) displayed higher expression at 21 dai in ‘Laurentian’ and at 14 dai in ‘Brutor’ (Figure 7B). In A. thaliana plants infected by the cyst nematode Heterodera schachtii, AAP1 along with other five AAP genes were upregulated to provide amino acids in the sink syncytial tissue, where the nematodes were developing (Elashry et al., 2013); mutants of these genes decreased the incidence of infection. Therefore, the expression of GPT2 and AAP1 in P. brassicae infected plants could provide amino acids and carbohydrates for pathogen nutrition and metabolic activity.
Another interesting gene in the ‘response to nematode' functional category was WRKY23. Transcripts from this gene were upregulated at 14 and 21 dai in both host genotypes, and to a lesser extent in ‘Brutor’ at 7 dai (Figure 7B). This transcription factor has been linked to polar transportation of auxin (Prat et al., 2018), a hormone that is at the core of increased host cell growth due to P. brassicae infection (Ludwig-Müller, 2009; Ludwig-Müller et al., 2009). Knocking down of this gene in A. thaliana resulted in reduced levels of infection by a pathogenic nematode (Grunewald et al., 2008). Since auxin is key for the establishment and colonization of both nematodes and P. brassicae, we hypothesize that both interactions trigger similar host responses.
Cell Wall Metabolism: Transcriptional Repression of Specific Cell Wall Genes can Contribute to Plant Defense
Cell wall modification is a common mechanism of defense against pathogen progression, and alterations to the cell wall structure may contribute to stopping disease. Two functional categories related to cell wall organization or biogenesis were enriched for downregulated genes at 7 and 14 dai in ‘Laurentian’ and at 7, 14 and 21 dai in ‘Brutor’ (Figure 4). Genes related to cell wall organization or biogenesis were upregulated at 7 dai and downregulated at 14 dai in resistant wild cabbage (Zhang et al., 2016b), and in A. thaliana roots challenged with the clubroot pathogen, the cell wall organization functional category was enriched with downregulated genes (Irani et al., 2018).
While most genes related to cell wall metabolism were downregulated in both hosts throughout the time course (Supplementary Figure S4A), other genes including expansins and xyloglucan endotransglucosylase/hydrolase (XTHs) showed a mixed response (Supplementary Figures S4B, C). Expansins and XTHs are important genes promoting cell growth through cell wall loosening, which favors P. brassicae infection (Devos et al., 2005; Siemens et al., 2006). These genes have been found to be upregulated in susceptible interactions (Agarwal et al., 2011; Irani et al., 2018), while symptomless B. oleracea plants and resistant A. thaliana interactions show downregulation of these genes (Jubault et al., 2013; Ciaghi et al., 2019). Except for a few transcripts annotated as cellulose synthase-like genes (Supplementary Figure S4C), cellulose synthases (CESAs) were largely downregulated in both hosts in the current study (Supplementary Figure S4A). Other genes involved in cell wall deposition, including cobra-like proteins (COBL) and irregular xylem genes (corresponding to CESAs and glycosyltransferases), and genes involved in cell wall lignification like peroxidase 25 (PRX25), were also mostly downregulated (Supplementary Figure S4A). Cellulose synthases (CESAs) control cellulose deposition in primary and secondary cell walls. Counterintuitively perhaps, A. thaliana CESAs mutants have been shown to be resistant to several pathogens (Ellis et al., 2002; Hernandez-Blanco et al., 2007; Malinovsky et al., 2014). The reaction of these mutants seems to be concomitant with a constitutive plant immune response, potentially through reallocation of resources to defense mechanisms. Overall, the regulation of cell wall genes indicates common mechanisms between the susceptible and moderately resistant host genotypes, with little marked differences between the two.
Contrasting and Shared Gene Expression Leads to Clues Regarding Susceptibility vs. Resistance
Contrasting Patterns of Gene Expression Between Cultivars at the Same Time Point
We created Venn diagrams for each time-point to identify genes with contrasting expression patterns. These genes can be key regulators and mark differences between susceptible and resistant interactions. Such genes are good candidates for gene editing-based mutagenesis to either validate gene function or increase host resistance. Genes downregulated in ‘Laurentian’ and upregulated in ‘Brutor’ could be susceptibility factors, while genes upregulated in ‘Laurentian' and downregulated in ‘Brutor’ could be involved in resistance.
At 7 dai, three genes were significantly downregulated in ‘Laurentian’ and upregulated in ‘Brutor’ (Figure 8). One of these transcripts corresponded to PRX37, a peroxidase potentially involved in cell wall cross-linking via ROS generation. This gene was coregulated with PDF1.2 in an A. thaliana overexpressor of the transcription factor DEWAX, which showed increased resistance to Botrytis cinerea (Ju et al., 2017).
Figure 8 Venn diagram of DEGs 7 dai. The diagram depicts the number of genes with common, opposite, and distinct expression patterns in Brassica napus ‘Laurentian' and ‘Brutor' in response to Plasmodiophora brassicae pathotype 5X. The annotation of genes that demonstrated opposite expression between the genotypes is shown. LA, Laurentian; BR, Brutor.
Most genes upregulated in ‘Laurentian’ and downregulated in ‘Brutor’ at 7 dai matched heat shock proteins (HSPs) (Figure 8). HSPs (e.g., HSP90) maintain the tridimensional integrity of proteins including membrane and cytosolic receptors and transductors of pathogenic signals (e.g., RPM1), (Hubert et al., 2003; Park et al., 2015). Complexes including cochaperones (HOP3) and HPS70 and HSP90 may be key to allowing correct protein folding in the ER during periods of high abiotic or biotic stress (Fernández-Bautista et al., 2017). The activation of these genes in the moderately resistant genotype (‘Laurentian’) is notable since, to our knowledge, there are only two reports of regulation of HSPs in this pathosystem (Zhao et al., 2017; Luo et al., 2018).
At 14 dai, the pectin methylesterase inhibitor (PMEI11) was upregulated in the resistant host and downregulated in the susceptible one (Figure 9). Previously, PMEI11 was found to be upregulated in response to two necrotrophic pathogens, A. thaliana challenged with B. cinerea (Lionetti et al., 2017), and B. napus infected with Sclerotinia sclerotiorum (Zhao et al., 2007). A second gene upregulated in ‘Laurentian’ and downregulated in ‘Brutor’ at 14 dai, anthranilate synthase alpha subunit (ASA1), was also regulated in response to S. sclerotiorum in the same study (Zhao et al., 2007). Along with PMEI11, ASA1 was one of 12 candidate genes with high expression and linkage to resistance QTLs in the B. napus–S. sclerotiorum interaction (Zhao et al., 2007).
Figure 9 Venn diagram of DEGs 14 dai. The diagram depicts the number of genes with common, opposite, and distinct expression patterns in Brassica napus ‘Laurentian’ and ‘Brutor’ in response to Plasmodiophora brassicae pathotype 5X. The annotation of genes that demonstrated opposite expression between the genotypes is shown. LA, Laurentian; BR, Brutor.
At 21 dai, nine genes were upregulated in the susceptible host and downregulated in the resistant one (Figure 10). One of these genes, protein phosphatase 2A (PP2A), is a key regulator of root polar auxin transport (Rashotte et al., 2007) and was implicated in brassinosteroid regulation upon P. brassicae infection (Schuller et al., 2014). This gene is an important candidate susceptibility factor, given the relationship between auxin regulation and clubroot development. Another gene upregulated in ‘Brutor’ and downregulated in ‘Laurentian’ at 21 dai was a tonoplast intrinsic protein 2 (TIP2). This gene was shown to interact specifically with an effector from a root knot nematode (Xue et al., 2013), suggesting an additional similarity between nematode-host and P. brassicae-host interactions.
Figure 10 Venn diagram of DEGs 21 dai. The diagram depicts the number of genes with common, opposite, and distinct expression patterns in Brassica napus ‘Laurentian’ and ‘Brutor’ in response to Plasmodiophora brassicae pathotype 5X. The annotation of genes that demonstrated opposite expression between the genotypes is shown. LA, Laurentian; BR, Brutor.
Some of the transcripts upregulated in ‘Laurentian’ and downregulated in ‘Brutor’ 21 dai (Figure 10) were discussed previously (e.g., SARD1), while the remaining genes indicated alterations in development, resource allocation and potential hijacking of the host molecular machinery. These included: i) a NAC-containing protein 83 (NAC083), also known as VND-interacting 2, which interacts with the VASCULAR-RELATED NAC-DOMAIN7 (VND7) to negatively regulate xylem formation (Yamaguchi et al., 2010); ii) a transport protein particle (TRAPPC10) that acts as a tethering factor to join vesicles that are placed on the cell plate during cytokinesis (Jaber et al., 2010); iii) a nitrate transporter (NPF7.2 also known as NRT1.8) that coordinates nitrogen allocation to the xylem (Zhang et al., 2014); and iv) a myo-inositol oxygenase (MIOX4), which reduces myo-inositol levels and thus galactinol levels that can activate defenses against the nematode H. schachtii (Siddique et al., 2014). The upregulation of these genes at 21 dai suggests that some of the defense processes in the moderately resistant host may have started to abate at the later stages of infection.
Collectively, our results indicate candidate genes that can be tested for validation of their potential function in resistance, and putative susceptibility factors that may increase resistance if mutated. Future candidate genes in studies comparing other genotypes with other pathotypes may be intersected with the genes found in the current study, to identify core elements of defense against clubroot. The common regulation of these genes will shed light on the basic or constitutive defense mechanisms used by most hosts against P. brassicae.
Shared Genes Between Cultivars at Different Time Points
Since we observed that the moderately resistant genotype presented delayed symptoms based on our phenotypical data, we performed comparisons of the number of genes shared between consecutive time-points between the two genotypes, to contrast these numbers with the common genes at the same time-point. For comparisons where the number of genes in consecutive time-points between cultivars was higher than the number of genes shared in the same time-point, functional enrichment analysis of Biological Process ‘level 3’ gene ontology categories was performed for shared genes in the same time-point and for shared genes in consecutive time-points. After eliminating common functional categories between these two analyses, the unique functional categories corresponding to shared genes in consecutive time-points were obtained.
The number of shared genes, when comparing upregulated or downregulated transcripts between the two genotypes at 7 dai, was always greater than the number of shared genes from either cultivar at 7 dai when compared with the other cultivar at 14 dai (Supplementary Table S5). The number of shared upregulated genes at 14 dai was 205, which was higher than the number of shared upregulated genes between ‘Laurentian’ at 14 dai and ‘Brutor’ at 21dai, but not higher than the number of shared upregulated genes between ‘Brutor’ at 14 dai and ‘Laurentian’ at 21 dai (247 genes). The functional enrichment analyses showed that most of the unique categories from genes in common between ‘Brutor’ at 14 dai and ‘Laurentian’ at 21 dai corresponded to nematode-related responses, cell growth and development and cell wall metabolism (Supplementary Table S5). This further supports similarities between the processes involved in nematode infection and clubroot development. In addition, the categories corresponding to cell growth and cell wall metabolism confirm that these processes are taking place earlier in ‘Brutor’, consistent with our phenotypic and histological observations. Finally, the number of shared downregulated genes between the two genotypes at 14 dai was smaller than the number of shared downregulated genes between ‘Laurentian’ at 14 dai and ‘Brutor’ at 21 dai, or ‘Brutor’ at 14 dai and ‘Laurentian’ at 21 dai (Supplementary Table S5). Unique functional categories being downregulated earlier in ‘Laurentian’ than in ‘Brutor’ indicate potential downregulation of carbohydrate metabolism, which usually allows better allocation of resources to defense mechanisms. In contrast, ‘Brutor’ seemed to turn off genes related to defense mechanisms (e.g., ROS metabolic process, defense response to other organism), earlier than ‘Laurentian’.
These analyses show how the resistant genotype is still affected by P. brassicae infection, but its molecular defense mechanisms may play a role in delaying symptom development.
Data Availability Statement
The datasets generated for this study can be found in the NCBI Sequence Read Archive (SRA) accession number PRJNA597078.
Author Contributions
LG-G designed and conducted all the experiments, performed all analyses and wrote the manuscript. VM assisted with experiments and provided expertise with respect to the inoculations and greenhouse material. S-FH helped secure project funding and contributed to the original experimental context. SS provided project guidance and edited several versions of the manuscript.
Funding
The authors thank the Canola Council of Canada and Agriculture and Agri-Food Canada for support provided via the Growing Forward 2 and Canadian Agricultural Partnership programs. Funding from Alberta Canola and in-kind support by the University of Alberta is also gratefully acknowledged.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We would like to thank Dr. T. Cao (University of Alberta) for assistance with material collection. Thanks to T. Locke from the Molecular Biology Service Unit, University of Alberta, for his advice on qRT-PCR validation experiments.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020.01025/full#supplementary-material
References
Agarwal, A., Kaul, V., Faggian, R., Rookes, J. E., Ludwig-Müller, J., Cahill, D. M. (2011). Analysis of global host gene expression during the primary phase of the Arabidopsis thaliana – Plasmodiophora brassicae interaction. Funct. Plant Biol. 38, 462–478. doi: 10.1071/FP11026
Asada, Y., Yamamoto, M., Tsutsui, T., Yamaguchi, J. (2011). The Arabidopsis NSL2 negatively controls systemic acquired resistance via hypersensitive response. Plant Biotechnol. 28, 9–15. doi: 10.5511/plantbiotechnology.10.0913a
Bardou, P., Mariette, J., Escudié, F., Djemiel, C., Klopp, C. (2014). jvenn : an interactive Venn diagram viewer. BMC Bioinf. 15, 293. doi: 10.1186/1471-2105-15-293
Benjamini, Y., Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
Breitenbach, H. H., Wenig, M., Wittek, F., Jordá, L., Maldonado-Alconada, A. M., Sarioglu, H., et al. (2014). Contrasting roles of the apoplastic aspartyl protease APOPLASTIC, ENHANCED DISEASE SUSCEPTIBILITY1-DEPENDENT1 and LEGUME LECTIN-LIKE PROTEIN1 in arabidopsis systemic acquired resistance. Plant Physiol. 165, 791–809. doi: 10.1104/pp.114.239665
Canola Council of Canada - industry overview (2017). Canola Counc. Canada. Available at: https://www.canolacouncil.org/markets-stats/industry-overview/.
Cao, H., Glazebrook, J., Clarke, J. D., Volko, S., Dong, X. (1997). The Arabidopsis NPR1 gene that controls systemic acquired resistance encodes a novel protein containing ankyrin repeats. Cell 88, 57–63. doi: 10.1016/S0092-8674(00)81858-9
Chalhoub, B., Tang, H., Wang, X., Chiquet, J., Belcram, H., Tong, C., et al. (2014). Early allopolyploid evolution in the post-Neolithic Brassica napus oilseed genome. Science 345, 950–953. doi: 10.1126/science.1253435
Chen, J., Pang, W., Chen, B., Zhang, C., Piao, Z. (2016a). Transcriptome analysis of Brassica rapa near-isogenic lines carrying clubroot-resistant and – susceptible alleles in response to Plasmodiophora brassicae during early infection. Front. Plant Sci. 6, 1–14. doi: 10.3389/fpls.2015.01183
Chen, T., Bi, K., He, Z., Gao, Z., Zhao, Y., Fu, Y. (2016b). Arabidopsis mutant bik1 exhibits strong resistance to Plasmodiophora brassicae. Front. Physiol. 7, 1–13. doi: 10.3389/fphys.2016.00402
Chu, M., Song, T., Falk, K. C., Zhang, X., Liu, X., Chang, A., et al. (2014). Fine mapping of Rcr1 and analyses of its effect on transcriptome patterns during infection by Plasmodiophora brassicae. BMC Genomics 15, 1166. doi: 10.1186/1471-2164-15-1166
Ciaghi, S., Schwelm, A., Neuhauser, S. (2019). Transcriptomic response in symptomless roots of clubroot infected kohlrabi (Brassica oleracea var. gongylodes) mirrors resistant plants. BMC Plant Biology 19, 288. doi: 10.1186/s12870-019-1902-z
Conesa, A., Stefan, G. (2008). Blast2GO : A Comprehensive Suite for Functional Analysis in Plant Genomics. Int. J. Plant Genom. 2008, 619832. doi: 10.1155/2008/619832
Conesa, A., Gotz, S., Garcia-Gomez, J. M., Terol, J., Talon, M., Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Dempsey, D. A., Vlot, A. C., Wildermuth, M. C., Klessig, D. F. (2011). Salicylic acid biosynthesis and metabolism. Arab. B. 9, e0156. doi: 10.1199/tab.0156
Devos, S., Vissenberg, K., Verbelen, J., Prinsen, E. (2005). Infection of Chinese cabbage by Plasmodiophora brassicae leads to a stimulation of plant growth : impacts on cell wall metabolism and hormone balance. New Phytol. 166, 241–250. doi: 10.1111/j.1469-8137.2004.01304.x
Donald, C., Porter, I. (2009). Integrated control of clubroot. J. Plant Growth Regul. 28, 289–303. doi: 10.1007/s00344-009-9094-7
Dong, X. (2004). NPR1 , all things considered. Curr. Opin. Plant Biol. 7, 547–552. doi: 10.1016/j.pbi.2004.07.005
Elashry, A., Okumoto, S., Siddique, S., Koch, W., Kreil, D. P., Bohlmann, H. (2013). The AAP gene family for amino acid permeases contributes todevelopment of the cyst nematode Heterodera schachtii in roots of Arabidopsis. Plant Physiol. Biochem. 70, 379–386. doi: 10.1016/j.plaphy.2013.05.016
Ellis, C., Karafyllidis, I., Wasternack, C., Turner, J. G. (2002). The Arabidopsis mutant cev1 links cell wall signaling to jasmonate and ethylene responses. Plant Cell 14, 1557–1566. doi: 10.1105/tpc.002022.1558
Enders, T. A., Frick, E. M., Strader, L. C., Drive, B., Louis, S. (2017). An Arabidopsis kinase cascade influences auxin-responsive cell expansion. Plant J. 92, 68–81. doi: 10.1111/tpj.13635
Ensembl (2018). EnsemblPlants. Brassica napus (AST_PRJEB5043_v1). Available at: http://plants.ensembl.org/Brassica_napus/Info/Index.
Fernández-Bautista, N., Fernández-Calvino, L., Muñoz, A., Castellano, M. M. (2017). HOP3 a new regulator of the ER stress response in Arabidopsis with possible implications in plant development and response to biotic and abiotic stresses. Plant Signal. Behav. 2324, 1–3. doi: 10.1080/15592324.2017.1317421
Fu, Z. Q., Yan, S., Saleh, A., Wang, W., Ruble, J., Oka, N., et al. (2012). NPR3 and NPR4 are receptors for the immune signal salicylic acid in plants. Nature 486, 228–232. doi: 10.1038/nature11162
Fujimoto, S. Y., Ohta, M., Usui, A., Shinshi, H., Ohme-Takagi, M. (2000). Arabidopsis ethylene-responsive element binding factors act as transcriptional activators or repressors of GCC box – mediated gene expression. Plant Cell 12, 393–404. doi: 10.1105/tpc.12.3.393
Grover, A., Taylor, P., Grover, A. (2012). Plant chitinases: genetic diversity and physiological roles. Crit. Rev. Plant Sci. 31, 57–73. doi: 10.1080/07352689.2011.616043
Grunewald, W., Karimi, M., Wieczorek, K., Van De Cappelle, E., Wischnitzki, E., Beeckman, T., et al. (2008). A role for AtWRKY23 in feeding site establishment of plant-parasitic nematodes. Plant Physiol. 148, 358–368. doi: 10.1104/pp.108.119131
Hamamouch, N., Li, C., Hewezi, T., Baum, T. J., Mitchum, M. G., Hussey, R. S., et al. (2012). The interaction of the novel 30C02 cyst nematode effector protein with a plant b-1,3-endoglucanase may suppress host defence to promote parasitism. J. Exp. Bot. 63, 3683–3696. doi: 10.1093/jxb/err313
Hasan, M. J., Rahman, H. (2016). Genetics and molecular mapping of resistance to Plasmodiophora brassicae pathotypes 2, 3, 5, 6, and 8 in rutabaga (Brassica napus var. napobrassica). Genome 59, 805–815. doi: 10.1139/gen-2016-0034
Hatakeyama, K., Suwabe, K., Tomita, R. N., Kato, T., Nunome, T., Fukuoka, H., et al. (2013). Identification and characterization of Crr1a , a gene for resistance to clubroot disease ( Plasmodiophora brassicae Woronin) in Brassica rapa L. PloS One 8, e54745. doi: 10.1371/journal.pone.0054745
Hernandez-Blanco, C., Feng, D. X., Hu, J., Sanchez-Vallet, A., Deslandes, L., Llorente, F., et al. (2007). Impairment of cellulose synthases required for Arabidopsis secondary cell wall formation enhances disease resistance. Plant Cell 19, 890–903. doi: 10.1105/tpc.106.048058
Hewezi, T., Piya, S., Qi, M., Balasubramaniam, M., Rice, J. H., Baum, T. J. (2016). Arabidopsis miR827 mediates post-transcriptional gene silencing of its ubiquitin E3 ligase target gene in the syncytium of the cyst nematode Heterodera schachtii to enhance susceptibility. Plant J. 88, 179–192. doi: 10.1111/tpj.13238
Howard, R. J., Strelkov And, S. E., Harding, M. W. (2010). Clubroot of cruciferous crops - new perspectives on an old disease. Can. J. Plant Pathol. 32, 43–57. doi: 10.1080/07060661003621761
Huang, D. W., Sherman, B. T., Tan, Q., Collins, J. R., Alvord, W. G., Roayaei, J., et al. (2007). The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 8, R183. doi: 10.1186/gb-2007-8-9-r183
Huang, D. W., Sherman, B. T., Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211
Hubert, D. A., Tornero, P., Belkhadir, Y., Krishna, P., Takahashi, A., Shirasu, K., et al. (2003). Cytosolic HSP90 associates with and modulates the Arabidopsis RPM1 disease resistance protein. EMBO J. 22, 5679–5689. doi: 10.1093/emboj/cdg547
Huffaker, A., Pearce, G., Ryan, C. A. (2006). An endogenous peptide signal in Arabidopsis activates components of the innate immune response. Proc. Natl. Acad. Sci. U. S. A. 103, 10098–10103. doi: 10.1073/pnas.0603727103
Hwang, S.-F., Strelkov, S. E., Feng, J., Gossen, B. D., Howard, R. J. (2012). Pathogen profile Plasmodiophora brassicae: a review of an emerging pathogen of the Canadian canola (Brassica napus) crop. Mol. Plant Pathol. 13, 105–113. doi: 10.1111/J.1364-3703.2011.00729.X
Hwang, S. F., Howard, R. J., Strelkov, S. E., Gossen, B. D., Peng, G. (2014). Management of clubroot (Plasmodiophora brassicae) on canola (Brassica napus) in western Canada. Can. J. Plant Pathol. 36, 49–65. doi: 10.1080/07060661.2013.863806
Irani, S., Trost, B., Waldner, M., Nayidu, N., Tu, J., Kusalik, A. J., et al. (2018). Transcriptome analysis of response to Plasmodiophora brassicae infection in the Arabidopsis shoot and root. BMC Genomics 19, 23. doi: 10.1186/s12864-017-4426-7
Jaber, E., Thiele, K., Kindzierski, V., Loderer, C., Rybak, K., Jürgens, G., et al. (2010). A putative TRAPPII tethering factor is required for cell plate assembly during cytokinesis in Arabidopsis. New Phytol. 187, 751–763. doi: 10.1111/j.1469-8137.2010.03331.x
Jahn, L., Mucha, S., Bergmann, S., Horn, C., Staswick, P., Steffens, B., et al. (2013). The clubroot pathogen ( Plasmodiophora brassicae) influences auxin signaling to regulate auxin homeostasis in Arabidopsis. Plants 1, 726–749. doi: 10.3390/plants2040726
Jia, H., Wei, X., Yang, Y., Yuan, Y., Wei, F., Zhao, Y., et al. (2017). Root RNA-seq analysis reveals a distinct transcriptome landscape between clubroot-susceptible and clubroot-resistant Chinese cabbage lines after Plasmodiophora brassicae infection. Plant Soil 421, 93–105. doi: 10.1007/s11104-017-3432-5
Ju, S., Go, Y. S., Choi, H. J., Park, J. M., Suh, M. C. (2017). DEWAX transcription factor is involved in resistance to Botrytis cinerea in Arabidopsis thaliana and Camelina sativa. Front. Plant Sci. 8, 1210. doi: 10.3389/fpls.2017.01210
Jubault, M., Lariagon, C., Simon, M., Delourne, R., Manzanares-dauleux, M. J. (2008). Identification of quantitative trait loci controlling partial clubroot resistance in new mapping populations of Arabidopsis thaliana. Theor. Appl. Genet. 117, 191–202. doi: 10.1007/s00122-008-0765-8
Jubault, M., Lariagon, C., Gravot, A., Delourme, R., Manzanares-dauleux, M. J. (2013). Partial resistance to clubroot in Arabidopsis is based on changes in the host primary metabolism and targeted cell division and expansion capacity. Funct. Integr. Genomics 13, 191–205. doi: 10.1007/s10142-013-0312-9
Jung, J. S.Æ.C. (2008). Genetic mapping of clubroot resistance genes in oilseed rape. Theor. Appl. Genet. 116, 363–372. doi: 10.1007/s00122-007-0674-2
Kageyama, K., Asano, Æ.T. (2009). Life Cycle of Plasmodiophora brassicae. J. Plant Growth Regul. 28, 203–211. doi: 10.1007/s00344-009-9101-z
Kersey, P. J., Allen, J. E., Allot, A., Barba, M., Boddu, S., Bolt, B. J., et al. (2018). Ensembl Genomes 2018 : an integrated omics infrastructure for non-vertebrate species. Nucleic Acid Res. 46, D802–D808. doi: 10.1093/nar/gkx1011
Kuhn, H., Lorek, J., Kwaaitaal, M., Consonni, C., Becker, K. (2017). Key components of different plant defense pathways are dispensable for powdery mildew resistance of the Arabidopsis mlo2 mlo6 mlo12 triple mutant. Front. Plant Sci. 8, 1006. doi: 10.3389/fpls.2017.01006
Langmead, B., Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–360. doi: 10.1038/nmeth.1923
Lee, M. W., Lu, H., Jung, H. W., Greenberg, J. T. (2007). A Key Role for the Arabidopsis WIN3 Protein in Disease Resistance Triggered by Pseudomonas syringae That Secrete AvrRpt2. Mol. Plant-Microbe Interact. 20, 1192–1200. doi: 10.1094/MPMI-20-10-1192
Lee, J., Izzah, N. K., Choi, B., Joh, H. J., Lee, S., Perumal, S., et al. (2016). Genotyping-by-sequencing map permits identification of clubroot resistance QTLs and revision of the reference genome assembly in cabbage (Brassica oleracea L.). DNA Res. 23, 29–41. doi: 10.1093/dnares/dsv034
Lemarie, S., Robert-seilaniantz, A., Lariagon, C., Lemoine, J., Manzanares-Dauleux, M. J., Gravot, A., et al. (2015). Both the jasmonic acid and the salicylic acid pathways contribute to resistance to the biotrophic clubroot agent Plasmodiophora brassicae in Arabidopsis. Plant Cell Physiol. 56, 2158–2168. doi: 10.1093/pcp/pcv127
Li, J., Brader, G., Kariola, T., Palva, E. T. (2006). WRKY70 modulates the selection of signaling pathways in plant defense. Plant J. 46, 477–491. doi: 10.1111/j.1365-313X.2006.02712.x
Li, L., Luo, Y., Chen, B., Xu, K., Zhang, F., Li, H., et al. (2016). A genome-wide association study reveals new loci for resistance to clubroot disease in Brassica napus. Front. Plant Sci. 7, 1483. doi: 10.3389/fpls.2016.01483
Lionetti, V., Fabri, E., De Caroli, M., Hansen, A. R., Willats, W. G. T., Piro, G., et al. (2017). Three pectin methylesterase inhibitors protect cell wall integrity for Arabidopsis immunity to Botrytis. Plant Physiol. 173, 1844–1863. doi: 10.1104/pp.16.01185
Livak, K. J., Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Lovelock, D. A., Donald, C. E., Conlan, X. A., Cahill, D. M. (2013). Salicylic acid suppression of clubroot in broccoli (Brassicae oleracea var. italica) caused by the obligate biotroph Plasmodiophora brassicae. Asutralas. Plant Pathol. 42, 141–153. doi: 10.1007/s13313-012-0167-x
Ludwig-Müller, J., Prinsen, E., Rolfe, S. A., Scholes, J. D. (2009). Metabolism and plant hormone action during clubroot disease. J. Plant Growth Regul. 28, 229–244. doi: 10.1007/s00344-009-9089-4
Ludwig-Müller, J., Auer, S., Jülke, S., Marschollek, S. (2017). Manipulation of auxin and cytokinin balance during the Plasmodiophora brassicae-Arabidopsis thaliana interaction Methods Mol. Biol. 1569, 41–60. doi: 10.1007/978-1-4939-6831-2
Ludwig-Müller, J. (2009). Plant defence – what can we learn from clubroots ? Australas. Plant Pathol. 38, 318–324. doi: 10.1071/AP09020
Luo, Y., Dong, D., Su, Y., Wang, X., Peng, Y., Peng, J., et al. (2018). Transcriptome analysis of Brassica juncea var. tumida Tsen responses to Plasmodiophora brassicae primed by the biocontrol strain Zhihengliuella aestuarii. Funct. Integr. Genomics 18, 301–314. doi: 10.1007/s10142-018-0593-0
Maekawa, S., Sato, T., Asada, Y. (2012). The Arabidopsis ubiquitin ligases ATL31 and ATL6 control the defense response as well as the carbon / nitrogen response. Plant Mol. Biol. 79, 217–227. doi: 10.1007/s11103-012-9907-0
Malinovsky, F. G., Fangel, J. U., Willats, W. G. T. (2014). The role of the cell wall in plant immunity. Front. Plant Sci. 5, 178. doi: 10.3389/fpls.2014.00178
Manoharan, R. K., Shanmugam, A., Hwang, I., Park, J.-I., Nou, I.-S. (2016). Expression of salicylic acid-related genes in Brassica oleracea var. capitata during Plasmodiophora brassicae infection. Genome 59, 379–391. doi: 10.1139/gen-2016-0018
Manzanares-Dauleux, M. J., Delourne, R., Baron, F., Thomas, G. (2000). Mapping of one major gene and of QTLs involved in resistance to clubroot in Brassica napus. Theor. Appl. Genet. 101, 885–891. doi: 10.1007/s001220051557
Matsumoto, E., Yasui, C., Ohi, M., Tsukada, M. (1998). Linkage analysis of RFLP markers for clubroot resistance and pigmentation in Chinese cabbage (Brassica rapa ssp. pekinensis). Euphytica 104, 79–86. doi: 10.1023/A:1018370418201
McGrath, K. C., Dombrecht, B., Manners, J. M., Schenk, P. M., Edgar, C. I., Udvardi, M. K., et al. (2005). Repressor- and activator-type ethylene response factors functioning in jasmonate signaling and disease resistance identified via a genome-wide screen of Arabidopsis transcription factor gene expression. Plant Physiol. 139, 949–959. doi: 10.1104/pp.105.068544.2001
Meier, S., Bastian, R., Donaldson, L., Murray, S., Bajic, V., Gehring, C. (2008). Co-expression and promoter content analyses assign a role in biotic and abiotic stress responses to plant natriuretic peptides. BMC PlantBiol. 8, 24. doi: 10.1186/1471-2229-8-24
Miller, R. N. G., Costa-Alves, G. S. C., Van-Sluys, M. (2017). Plant immunity: unravelling the complexity of plant responses to biotic stresses. Ann. Bot. 119, 681–687. doi: 10.1093/aob/mcw284
Moreau, M., Tian, M., Klessig, D. F. (2012). Salicylic acid binds NPR3 and NPR4 to regulate NPR1-dependent defense responses. Cell Res. 22, 1631–1633. doi: 10.1038/cr.2012.100
NIH (2012). Primer3. Available at: http://bioinfo.ut.ee/primer3-0.4.0/ Accessed 25 September 2018.
Nobuta, K., Okrent, R. A., Stoutemyer, M., Rodibaugh, N., Kempema, L., Wildermuth, M. C., et al. (2007). The GH3 acyl adenylase family member PBS3 regulates salicylic acid-dependent defense responses in Arabidopsis. Plant Physiol. 144, 1144–1156. doi: 10.1104/pp.107.097691
Park, C., Seo, Y., Park, C. (2015). Heat shock proteins: a review of the molecular chaperones for plant immunity. Plant Pathol. J. 31, 323–333. doi: 10.5423/PPJ.RW.08.2015.0150
Penninckx, I. A. M. A., Eggermont, K., Terras, F. R. G., Bart, P. H. J., De Samblanx, G. W., Buchala, A., et al. (1996). Pathogen-induced systemic activation of a plant defensin gene in Arabidopsis follows a salicylic acid-independent pathway. Plant Cell 8, 2309–2323. doi: 10.1105/tpc.8.12.2309
Pfaffl, M. W., Tichopad, A., Prgomet, C., Neuvians, T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper–Excel-based tool using pair-wise correlations. Biotechnol. Lett. 26, 509–515. doi: 10.1023/B:BILE.0000019559.84305.47
Prat, T., Hanjny, J., Grunewald, W., Vasileva, M., Molnar, G., Tejos, R., et al. (2018). WRKY23 is a component of the transcriptional network mediating auxin feedback on PIN polarity. PloS Genet. 14, e1007177. doi: 10.1371/journal.pgen.1007177
Rashotte, A. M., DeLong, A., Muday, G. K. (2007). Genetic and chemical reductions in protein phosphatase activity alter auxin transport, gravity response, and lateral root growth. Plant Cell 13, 1683. doi: 10.2307/3871394
Rodiuc, N., Vieira, P., Banora, M. Y., Engler, J. D. A. (2014). On the track of transfer cell formation by specialized plant-parasitic nematodes. Front. Plant Sci. 5, 1–14. doi: 10.3389/fpls.2014.00160
Saeed, A. I., Sharov, V., White, J., Li, J., Liang, W., Bhagabati, N., et al. (2003). TM4: a free, open-source system for microarray data management and analysis. Biotehcniques 34, 374–378. doi: 10.2144/03342mt01
Schuller, A., Kehr, J., Ludwig-Müller, J. (2014). Laser microdissection coupled to transcriptional profiling of Arabidopsis roots inoculated by Plasmodiophora brassicae indicates a role for brassinosteroids in clubroot formation. Plant Cell Physiol. 55, 392–411. doi: 10.1093/pcp/pct174
Seyfferth, C., Tsuda, K. (2014). Salicylic acid signal transduction: the initiation of biosynthesis , perception and transcriptional reprogramming. Front. Plant Sci. 5, 1–10. doi: 10.3389/fpls.2014.00697
Siddique, S., Bohlmann, H., Endres, S., Tenhaken, R., Sobczak, M., Radakovic, Z. S., et al. (2014). Myo-inositol oxygenase is important for the removal of excess myo-inositol from syncytia induced by Heterodera schachtii in Arabidopsis roots. New Phytol. 201, 476–485. doi: 10.1111/nph.12535
Siemens, J., Keller, I., Sarx, J., Kunz, S., Schuller, A., Nagel, W., et al. (2006). Transcriptome analysis of Arabidopsis clubroots indicate a key role for cytokinins in disease development. Mol. Plant-Microbe Interact. 19, 480–494. doi: 10.1094/MPMI-19-0480
Son, K., Yu, S., Shin, W., Han, K., Kang, K. (2018). A simple guideline to assess the characteristics of RNA-Seq Data. BioMed. Res. Int. 2018, 1–9. doi: 10.1155/2018/2906292
Stael, S., Kmiecik, P., Willems, P., Van Der Kelen, K., Coll, N. S., Teige, M., et al. (2015). Plant innate immunity – sunny side up ? Trends Plant Sci. 20, 3–11. doi: 10.1016/j.tplants.2014.10.002
Strelkov, S., Manolii, V. P., Harding, M. W. (2015). The spread of clubroot on canola in Alberta in 2014. Can. Plant Dis. Surv. 95, 155–158.
Strelkov, S. E., Hwang, S., Manolii, V. P., Cao, T., Feindel, D. (2016a). Emergence of new virulence phenotypes of Plasmodiophora brassicae on canola (Brassica napus) in Alberta , Canada. Eur. J. Plant Pathol. 145, 517–529. doi: 10.1007/s10658-016-0888-8
Strelkov, S. E., Manolii, V. P., Harding, M. W., Hwang, S., Manolii, E., Zuzak, K., et al. (2016b). Occurrence and spread of clubroot on canola in Alberta in 2016. Can. Plant Dis. Surv. 97, 164–167.
Strelkov, S. E., Hwang, S. F., Manolii, V. P., Cao, T., Fredua-Agyeman, R., Harding, M. W., et al. (2018). Virulence and pathotype classification of Plasmodiophora brassicae populations collected from clubroot resistant canola (Brassica napus) in Canada. Can. J. Plant Pathol. 40, 284–298. doi: 10.1080/07060661.2018.1459851
Strelkov, S. E., Manolli, V. P., Lageyre, J., Hwang, S. F., Harding, M. W., Daniels, G. C. (2019). Occurrence and spread of clubroot on canola in Alberta in 2018 - Canadian Plant Disease Survey.
Suwabe, K., Tsukazaki, H., Iketani, H., Hatakeyama, K., Fujimura, M., Nunome, T., et al. (2003). Identification of two loci for resistance to clubroot (Plasmodiophora brassicae Woronin) in Brassica rapa L. Theor. Appl. Genet. 107, 997–1002. doi: 10.1007/s00122-003-1309-x
Szadkowski, E., Eber, F., Huteau, V., Lodé, M., Coriton, O., Jenczewski, E., et al. (2011). Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus. New Phytol. 191, 884–894. doi: 10.1111/j.1469-8137.2011.03729.x
Trapnell, C., Pachter, L., Salzberg, S. L. (2009). TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105–1111. doi: 10.1093/bioinformatics/btp120
Ueno, H., Matsumoto, E., Aruga, D. (2012). Molecular characterization of the CRa gene conferring clubroot resistance in Brassica rapa. Plant Mol. Biol. 80, 621–629. doi: 10.1007/s11103-012-9971-5
Untergasser, A., Cutcutache, I., Koressaar, T., Ye, J., Faircloth, B. C., Remm, M., et al. (2012). Primer3–new capabilities and interfaces. Nucleic Acids Res. 40, e115. doi: 10.1093/nar/gks596
Wagner, S., Stuttmann, J., Rietz, S., Guerois, R., Brunstein, E., Bautor, J., et al. (2013). Structural basis for signaling by exclusive EDS1 heteromeric complexes with SAG101 or PAD4 in plant innate immunity. Cell Host. Microbe 14, 619–630. doi: 10.1016/j.chom.2013.11.006
Wang, L., Tsuda, K., Truman, W., Sato, M., Nguyen, L. V., Katagiri, F., et al. (2011a). CBP60g and SARD1 play partially redundant critical roles in salicylic acid signaling. Plant J. 3000, 1029–1041. doi: 10.1111/j.1365-313X.2011.04655.x
Wang, Y. H., Gehring, C., Irving, H. R. (2011b). Plant natriuretic peptides are apoplastic and paracrine stress response molecules. Plant Cell Physiol. 52, 837–850. doi: 10.1093/pcp/pcr036
Wang, S., Yu, F., Zhang, W., Tang, J., Li, J., Yu, L., et al. (2019). Comparative transcriptomic analysis reveals gene expression changes during early stages of Plasmodiophora brassicae cab infection in cabbage (Brassica oleracea var . capitata L .). Can. J. Plant Pathol. 41, 188–199. doi: 10.1080/07060661.2019.1567592
Wasternack, C., Hause, B. (2013). Jasmonates: biosynthesis, perception, signal transduction and action in plant stress response, growth and development. An update to the 2007 review in Annals of Botany. Ann. Bot. 111, 1021–1058. doi: 10.1093/aob/mct067
Weigel, R. R., Pfitzner, U. M., Gatz, C. (2005). Interaction of NIMIN1 with NPR1 modulates PR gene expression in Arabidopsis. Plant Cell 17, 1279–1291. doi: 10.1105/tpc.104.02744
Williams, P. H. (1966). A system for the determination of races of Plasmodiophora brassicae that infect Cabbage and Rutabaga. Phytopathology 56, 624–626.
Xu, Y., Xu, H., Wu, X., Fang, X., Wang, J. (2012). Genetic changes following hybridization and genome doubling in synthetic Brassica napus. Biochem. Genet. 50, 616–624. doi: 10.1007/s10528-012-9505-5
Xue, B., Hamamouch, N., Li, C., Huang, G., Hussey, R. S., Baum, T. J., et al. (2013). The 8D05 Parasitism Gene of Meloidogyne incognita Is Required for Successful Infection of Host Roots. Phytopathology 103, 175–181. doi: 10.1094/phyto-07-12-0173-r
Yaeno, T., Iba, K. (2008). BAH1/NLA, a RING-type ubiquitin E3 ligase, regulates the accumulation of salicylic acid and immune responses to Pseudomonas syringae DC3000. Plant Physiol. 148, 1032–1041. doi: 10.1104/pp.108.124529
Yamaguchi, M., Ohtani, M., Mitsuda, N., Kubo, M., Ohme-Takagi, M., Fukuda, H., et al. (2010). VND-INTERACTING2, a NAC domain transcription factor, negatively regulates xylem vessel formation in Arabidopsis. Plant Cell 22, 1249–1263. doi: 10.1105/tpc.108.064048
Yang, Z., Tian, L., Latoszek-green, M., Brown, D., Wu, K. (2005). Arabidopsis ERF4 is a transcriptional repressor capable of modulating ethylene and abscisic acid responses. Plant Mol. Biol. 58, 585–596. doi: 10.1007/s11103-005-7294-5
Yu, F., Zhang, X., Huang, Z., Chu, M., Song, T. (2016). Identification of genome-wide variants and discovery of variants associated with Brassica rapa clubroot resistance gene Rcr1 through bulked segregant RNA sequencing. PloS One 11, e0153218. doi: 10.1371/journal.pone.0153218
Zhang, Y., Xu, S., Ding, P., Wang, D., Ti, Y., He, J., et al. (2010). Control of salicylic acid synthesis and systemic acquired resistance by two members of a plant-specific family of transcription factors. Proc. Natl. Acad. Sci. 107, 18220–18225. doi: 10.1073/pnas.1005225107
Zhang, G.-B., Yi, H.-Y., Gong, J.-M. (2014). The Arabidopsis ethylene/jasmonic acid-NRT signaling module coordinates nitrate reallocation and the trade-off between growth and environmental adaptation. Plant Cell 26, 3984–3998. doi: 10.1105/tpc.114.129296
Zhang, H., Feng, J., Hwang, S., Strelkov, S. E., Falak, I., Huang, X., et al. (2016a). Mapping of clubroot (Plasmodiophora brassicae) resistance in canola (Brassica napus). Plant Pathol. 65, 435–440. doi: 10.1111/ppa.12422
Zhang, X., Liu, X., Fang, Z., Zhansheng, L., Yang, L., Zhuang, M., et al. (2016b). Comparative transcriptome analysis between broccoli (Brassica oleracea var. italica) and wild cabbage (Brassica macrocarpa Guss.) in response to Plasmodiophora brassicae during different infection stages. Front. Physiol. 7, 1929. doi: 10.3389/fpls.2016.01929
Zhao, J., Wang, J., An, L., Doerge, R. W., Chen, Z. J., Grau, C. R., et al. (2007). Analysis of gene expression profiles in response to Sclerotinia sclerotiorum in Brassica napus. Planta 227, 13–24. doi: 10.1007/s00425-007-0586-z
Keywords: Brassica napus, Plasmodiophora brassicae, clubroot, RNA-seq, salicylic acid, nematodes, immunity
Citation: Galindo-González L, Manolii V, Hwang S-F and Strelkov SE (2020) Response of Brassica napus to Plasmodiophora brassicae Involves Salicylic Acid-Mediated Immunity: An RNA-Seq-Based Study. Front. Plant Sci. 11:1025. doi: 10.3389/fpls.2020.01025
Received: 18 February 2020; Accepted: 23 June 2020;
Published: 09 July 2020.
Edited by:
Dilantha Fernando, University of Manitoba, CanadaReviewed by:
Jacqueline Batley, University of Western Australia, AustraliaRobert Malinowski, Institute of Plant Genetics (PAN), Poland
Copyright © 2020 Galindo-González, Manolii, Hwang and Strelkov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Stephen E. Strelkov, stephen.strelkov@ualberta.ca