The two B. rapa accessions, BrT24 (resistant, R-line) and Y510-9 (susceptible, S-line) having contrasting resistance to P. brassicae were obtained from the Institute of Horticulture, Henan Academy of Agricultural Sciences, China. The P. brassicae strain used in this study was collected from clubroot of B. rapa in Xinye County, Henan Province, China and identified as race 4 by the Williams system (Yuan et al., 2017). The P. brassicae spore suspension, diluted to 1 × 10⁷/ml was used as inoculum (Luo et al., 2014).

The seeds of R- and S-lines were germinated on wet filter paper followed by transplanting into 50-well plastic trays and incubating at a temperature of 25°C/20°C and photoperiod of 16 h/8 h (light/dark). After 20 days, each well was inoculated with 20 ml of P. brassicae solution with the control group receiving 20 ml of sterile water. Root samples were collected on four different time points representing four stages of disease development namely, before inoculation (0 days after inoculation, DAI), at the time of cortex infection (3 DAI), early onset (9 DAI), and later stage (20 DAI) of the disease.

For microscopic observation, the fresh tissue was fixed with fixed liquid for more than 36 h. The tissue was removed from the fixed liquid and trimmed with a scalpel in the ventilation cupboard; the trimmed tissue and the label were placed in the dehydration box; and the dehydration box was put into the dehydrator in order to dehydrate with gradient alcohol. The wax-soaked tissue was embedded in the embedding machine. Finally, placed the trimmed tissue wax block on the freezing table to cool at −20°C and sliced the modified tissue chip wax block on the paraffin slicer (with slice thickness: 4 μm). The tissue was flattened when the slice floated on the 40°C warm water of the spreading machine, and the tissue was picked up by the glass slides and baked in the oven at 60°C. After the water-baked dried wax was melted, it was taken out and stored at room temperature.

Toluidine blue staining was performed according to the following procedure: the slices were kept in xylene I for 20 min and in xylene II for 20 min. Then, they were kept in 100% ethanol I—100% ethanol II–75% ethanol, each step took around 5 min, and rinsed with tap water. Then the plant tissue slices were treated with toluidine blue for 2–5 min and rinsed with tap water. Then they were microscopically inspected, and according to the color of the tissues, checked for their differentiation. After tap water wash, they were dried in the oven. Finally, they were kept in xylene for 10 min, sealed with neutral gum, and followed microscope inspection, image acquisition, and analysis.

The RNA extraction and sequencing was performed by Beijing Biomarker Biotechnology Co., Ltd., Beijing, China. There, the RNA samples were assessed by 1% agarose gel electrophoresis for RNA degradation and impurities, Nanodrop 2000 for RNA purity and concentration, and Agilent 2100 RNA 6000 Nano Kit for RNA integrity (Agilent, Santa Clara, CA, United States). The Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, United States) was used to remove rRNA. Sequencing libraries were generated using NEBNext Ultra^TM RNA Library Prep Kit for Illumina^® [New England Biolabs (NEB), Ipswich, MA, United States] following the recommendations of the manufacturer, and index codes were added to attribute sequences to each sample. mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First-Strand Synthesis Reaction Buffer (5×. First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments, preferentially of 240 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, United States). About 3 μl USER Enzyme (NEB) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Sequencing was done using the Illumina HiSeq^TM 2500 platform for paired-end sequencing. The raw data were cleaned by removing the adapter sequence, N-sequence, and low-quality reads. Then, Q20, Q30, GC contents, and sequence repetition levels were calculated. The obtained clean data were then compared with the reference genome of Chinese cabbage (version 1.5).¹ Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped (FPKM) (Florea et al., 2013). Differential expression analysis of two samples was performed using the edgeR (Wang et al., 2009). The false discovery rate (FDR) of <0.01 and fold change of ≥2 was set as the threshold for determining significant differential expression. The raw reads were further processed with a bioinformatic pipeline tool, BMKCloud.² Co-expression network analysis using WGCNA version 1.61 software package in R software (Langfelder and Horvath, 2008). The interaction network of hub-genes in module was visualized using Cytoscape 3.7.1.

The samples were ground with liquid nitrogen. About 80 ± 5 mg of finely ground samples were then taken in a 2-ml centrifuge tube prior to adding 50 μl of internal standard solution and 1 ml of acetonitrile aqueous solution (1% FA). The mixtures were mixed well by shaking for 2 min and left for 12 h at 4°C in the dark for extraction. After that the samples were centrifuged at 14,000 × g for 10 min, and 800 μl of the supernatant were blow dried with nitrogen prior to reconstituting with 200 μl of acetonitrile water (1:1, v/v). This was further centrifuged at 14,000 × g for 10 min, and the supernatant was taken for analysis. The samples were separated by Agilent 1290 Infinity LC ultra-high-performance liquid chromatography system and analyzed by mass spectrometry. Selective reaction/multiple reaction monitoring (MRM) technology mode was used to detect the ion pair to be tested. Data processing and analysis were conducted using the Thermo Scientific Xcalibur 2.1 software package (Thermo Fisher Scientific, San Jose, CA, United States). The contents of plant hormones in the samples were estimated based on the standard curve.

Upon inoculation, no obvious clubroot symptom was observed in both the R-line “BrT24” and S- line “Y510-9” till 9 DAI. At 20 DAI, conspicuous irregular swelling was observed in the roots of S-line, whereas no signs of swelling were found in the roots of R-line (Figure 1). These results indicate that clubroot develops more rapidly and severely on the susceptible genotype.

Cross sections of paraffin-embedded roots indicate no obvious changes in the S- and R-lines at 3 DAI (Figures 2B,b) and there is no significant change compared with the uninoculated period (Figures 2A,a). At 9 DAI, irregularly enlarged and squeezed cells are observed in the roots of S-line (Figure 2C), whereas in the R-line, the root cells were not compressed or deformed (Figure 2C). At 20 DAI, S-line was severely infected by P. brassicae, R-line was slightly infested, and there is no obvious cell swelling and squeezing in R-line (Figure 2d). However, in S-line, the enlarged cells contain a large number of dormant spores, the root cortex cells are compressed by the enlarged cells, and the deformed cells are arranged disorderly (Figure 2D).

A total of 24 libraries from the three replicates of the root RNA samples of R- and S-line at 0, 3, 9, and 20 DAI were constructed and analyzed. A total of 291.23 Gb clean data (10.42 Gb per sample) was obtained. The Q30 values were all above 94.64%, indicating that the quality and accuracy of sequencing data were sufficient for subsequent analyses (Supplementary Table 1). Principal component analysis (PCA) plot of gene expression levels of the samples of R-line indicated that samples of 3, 9, and 20 DAI clustered far from the samples of 0 DAI (separated by PC1). In S-line, the samples of 0, 3, and 9 DAI were clustered far from the samples of 20 DAI (separated by PC1). This indicated that the S-line did not respond promptly to the root disease on the 3 and 9 DAI (Figure 3A).

A total of 41,087 transcripts were detected across all samples (Supplementary Table 2). Compared to 0 DAI, at 3 DAI, 3507, and 3,362 genes were upregulated and downregulated, respectively, in R-line, whereas 2,204 and 2,356 genes were upregulated and downregulated, respectively, in S-line (Figure 3B). At 9 DAI, 3,912 genes were upregulated, and 4,151 genes were downregulated in R-line, whereas only 3,045 genes were upregulated and 2,843 genes were downregulated in S-line. At 20 DAI, 3,880 and 3,660 genes were upregulated and downregulated in R-line, whereas 3,558 and 4,866 genes were upregulated and downregulated in S-line (Figure 3B). The FPKM values of the genes are shown in Supplementary Table 3.

Overall, the change of gene expression patterns of 12 selected differentially expressed genes (DEGs) (Supplementary Table 4) were examined by quantitative real-time PCR (qRT-PCR) to validate the RNA-seq results. The high correlation coefficients of log2 fold changes obtained from RNA-seq and qRT-PCR results suggested that the RNA-seq data in this study is reliable (Supplementary Figure 1).

The genes that are contrastingly expressed in R- and S-lines may play key roles in plants responses toward P. Brassicae infection. Genes downregulated in R-line and upregulated in S-line could be susceptibility factors, whereas genes upregulated in R-line and downregulated in S-line could be involved in resistance. Such genes are good candidates for gene editing-based mutagenesis to either validate gene function or increase host resistance.

At 3 DAI, 29 genes were significantly downregulated in R-line but upregulated in S-line. On the contrary, 60 genes were upregulated in R-line and downregulated in S-line at 3 DAI (Supplementary Figure 2A). At 9 DAI, 30 genes were upregulated in the S-line and downregulated in R-line, 35 genes were upregulated in R-line and downregulated in S-line at 9 DAI (Supplementary Figure 2B). At 20 DAI, 209 genes were upregulated in S-line and downregulated in R-line, 868 genes were upregulated in R-line and downregulated in S-line at 20 DAI (Supplementary Figure 2C).

In order to exclude the background difference between resistant and susceptible materials, we analyzed the relative transcriptional changes over time. In R- or S-lines, DEGs at 3, 9, and 20 DAI were identified by comparing the expression levels at 3 DAI with those at 0 DAI, the level at 9 DAI with those at 3 DAI, and the level at 20 DAI with those at 9 DAI (Supplementary Figure 3A). In R-line, 3,507, 1,119, and 136 DEGs were upregulated, and 3,362, 1,391, and 67 DEG were downregulated at 3, 9, and 20 DAI, respectively. In S-line, 2,204, 1,285, and 2,648 DEGs were upregulated, and 2,356, 589, and 4,549 DEGs were downregulated at 3, 9, and 20 DAI, respectively (Supplementary Figure 3A). Compared with S-line, 2,475, 2,182, and 4,292 relative DEGs were identified in R-line at 3, 9, and 20 DAI, respectively (Supplementary Figure 3B). The above results indicate that the number of DEGs in 3 DAI of R-line is relatively large compared with other time points, indicating that 3 DAI is the key time point for R-line to respond to clubroot. It also indicates that the R-line activates the defensive response earlier compared to that of S-line.

The weighted gene co-expression network analysis (WGCNA) of 6,092 valid genes (Supplementary Table 3) was carried out. The genes were assigned into 10 distinct modules (Figure 4A). A module is a cluster of highly interconnected genes with similar expression changes in a physiological process. The modules were then associated with each of the samples, which identified seven modules namely, blue, darkgray, darkorange, greenyellow, darkgreen, tan, and brown to be highly correlated with R- or S-line (Figures 4B,C). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these seven modules revealed that the genes of three modules (blue, darkgray, and greenyellow) highly significantly enriched six pathways namely, Plant–pathogen interaction, Plant hormone signal transduction, Phenylalanine metabolism, Phenylpropanoid biosynthesis, Taurine and hypotaurine metabolism and Nitrogen metabolism (Table 1).

Gene interaction network analysis of the genes of these three modules was done using Cytoscape 3.7.1, which revealed a total of 15 hub genes (Figures 5A–C and Table 2). Among these, the notable genes include cleavage site for pathogenic type III effector avirulence factor avr 4 (RIN4), auxin-responsive protein gene (IAA16) of plant hormone signal transduction pathway, universal stress protein family (Bra019995), peroxidase 44 (PER44) under phenylpropanoid biosynthesis pathway, protein detoxification 9 (DTX9) gene, HSF-type DNA-binding gene (HSFA6b), probable E3 ubiquitin-protein ligase (ARI16), b-cell receptor-associated 31-like protein (Bra001681), glutamate decarboxylase 1 (GAD1), glutamate decarboxylase 2 (GAD2), and s-adenosylmethionine synthetase (SAMS3) in cysteine and methionine metabolism, carbohydrate-binding protein of the ER (SIRK), and Brassica_rapa_new Gene_17466.

RIN4 is a negative regulator of PTI, since it interacts with Pseudomonas syringae type III effector molecules, and its overexpression has been shown to increase PAMP-triggered cell wall thickening in Arabidopsis (Mackey et al., 2002; Hou et al., 2011). This RIN4 may have a similar role in CR in B. rapa. IAA16 promotes root responsiveness to ABA and is involved with growth of roots, cotyledons, and hypocotyls of Arabidopsis seedlings (Schmid et al., 2005; Rinaldi et al., 2012). Peroxidase-generated apoplastic reactive oxygen species (ROS) were found to compromise cuticle integrity and contribute to damage-associated molecular pattern (DAMP)-elicited defenses in Arabidopsis. It will be interesting to investigate in future if PER44 and another detoxification-related gene, DTX9, identified as hub genes in our study, have any roles in resistance against CR (Mantas et al., 2016). Overexpression of glutamate decarboxylase conferred resistance to the northern root-knot nematode in tobacco (Mantas et al., 2016). The two hub glutamate decarboxylase genes, GAD1 and GAD2, identified in our study, may have similar roles against clubroot. Overall, these hub genes may play pivotal roles in conferring resistance against clubroot disease in B. rapa.

The KEGG enrichment analyses of the upregulated and downregulated DEGs in R- and S-lines (Supplementary Figure 4), contrasting DEGs between R- and S-lines (Supplementary Figure 2) and the seven modular genes (Table 1), were conducted to understand the biological mechanisms of CR. The key significantly enriched pathways of all these KEGG enrichment analyses are summarized in Figure 6 and the list of genes under each significant KEGG pathways are shown in Supplementary Table 5.

Overview of all the significantly enriched KEGG pathways revealed that the plant hormone signal transduction pathway is only enriched by the upregulated genes of R-line at 3 DAI (but not at 9 and 20 DAI). Contrastingly, this pathway was enriched by the downregulated genes of 9 and 20 DAI in R-line and all time points in S-line. This contrasting signal transduction in R- vs. S-lines clearly indicates that hormone signal transduction was activated and thereby, plants immune responses were initiated at an earlier stage of infection in R-line. Besides, enrichment of the pathway by upregulated DEGs of R-line and downregulated DEGs of S-line at 9 and 20 DAI indicate that hormone signaling based immune response is still in operation in R-line which is absent in S-line. This is also evident from the enrichment of the plant–pathogen interaction pathway by the upregulated DEGs at all time-points in R-line whereas only by the upregulated DEGs at 3 DAI in S-line. This indicates that plants interaction with the pathogen is occurring at all times in R-line whereas this immune response is inhibited at later time points in S-line. To exclude the background difference, we performed the KEGG enrichment analysis of the relatively DEGs. For relative DEGs, the genes involved in CR regulation are indeed enriched in plant–pathogen interaction, plant hormone signal transduction, and phenylpropanoid biosynthesis (data not shown).

Brassinosteroid biosynthesis and Zeatin biosynthesis were enriched by the downregulated DEGs of R-line at different time points. These pathways were not enriched by the up- or downregulated DEGs of S-line. Isoquinoline alkaloid biosynthesis is only significantly enriched in genes that are upregulated in R-line and downregulated in S-line by 3 DAI.

Contrastingly, glyoxylate and dicarboxylate metabolism pathway is enriched by the downregulated DEGs of S-line at initial stages of infection, and GSL biosynthesis is only enriched by the downregulated DEGs of S-line at a later stage of infection. Glyoxylate and dicarboxylate metabolism is known to play important functions in the detoxification during stress and contribute to redox balance in plants (Allan et al., 2009). GSLs, a group of sulfur-containing plant secondary metabolites found in the cruciferous family, play an important role in the resistance of Brassica root disease (Ludwig-Müller et al., 2009; Zhang et al., 2016; Li et al., 2019).

Upregulated DEGs (at all time points) of both R- and S-lines are enriched in glutathione metabolism pathway, whereas phenylpropanoid biosynthesis, phenylalanine metabolism, and phenylalanine, tyrosine, and tryptophan biosynthesis pathways were enriched by the upregulated DEGs of all time points in R-line and by the upregulated DEGs of initial time points and downregulated DEGs of later time points in S-line (Figure 6). This indicates that these pathways can increase the resistance during the secondary infection process of clubroot. A similar observation was also observed by Irani et al. (2019).

The pathway enrichment analysis of DEGs showed that “Plant hormone signal transduction,” “Alanine, aspartate and glutamate metabolism,” and “Phenylpropanoid biosynthesis” pathways are involved in the immune response to P. brassicae. Phytohormones play an important role in the different growth and development processes of plants and various biotic and abiotic stress responses. Early studies have shown that different pathogens infect plants can cause changes in the levels of different phytohormones (Adie et al., 2007; Robert-Seilaniantz et al., 2007). In order to further explore the role of plant hormones in defense against P. brassicae, we quantified the contents of phytohormone and examined the changes in expression of genes of related pathways in the contrastingly resistant lines at different time points after inoculation, which is discussed thoroughly in the following section.