For the clubroot resistance study, we used B. rapa (ECD04) derived from the European Clubroot Differential Set, which genome has been completed assembled in our previous study (Yang et al., 2022a). ECD04 roots were inoculated with P. brassicae strain collected from Zhijiang (Hubei, China) (Chen et al., 2016), root sample treated with water was used as mock. In detail, the collected galls were thawed at room temperature, mixed with buffer, crushed and filtered, and the spore concentration was diluted to 10⁷ resting spores per milliliter in sterile distilled water used for the infection experiments. The germinated seeds were sown into the soil and after 7 days 1 ml of spore suspension was applied to their roots. Samples were collected at 18, 24, or 30 days in triplicates, resulting in 18 cDNA libraries including I18, I24 and I30 (inoculated), and M18, M24 and M30 (mock). Total RNA was extracted using TRIzol (Invitrogen Co Ltd, USA) according to the manufacturer’s instructions. The same set of samples were used for RNA-seq and sRNA-seq.

At about 35 days after inoculation, disease symptoms can be classified as 0-3 grade based on their root morphology, roots with no gall were classified as grade 0, lateral roots with small galls were classified as grade 1, lateral roots with large galls were classified as grade 2, and main roots with large galls were classified as grade 3. Based on the number of individual plants in different disease classes, the disease index (DI) of the material was calculated to assess the disease resistance of the material (Strelkov et al., 2006).

Where n is the number of plants in each grade; N is the total number of plants; and 0, 1, 2, and 3 are the disease symptom grade.

cDNA library preparation and sequencing were performed using Hifair II First Strand cDNA Synthesis Kit (Beijing Genome Institute, China). Messager RNA was isolated using oligo (dT) primers, fragmented and cDNA was synthesized with random primers. Second-strand cDNA was synthesized using buffer, deoxyribonucleotide triphosphates (dNTPs), RNase H, and DNA polymerase I. The synthetic cDNA was subjected to end repair by polymerase with Seloxa adapters, and suitable fragments were retrieved by agarose gel electrophoresis. Polymerase chain reaction (PCR) amplification was carried out to enrich the purified fragments. RNA-Seq was carried out using SEQ-500 platform (BGI, China).

Gene expression levels were analyzed by R v3.6.2 package “DESeq2” (1.10.1). DEGs were identified based on 1) adjusted P<0.05, 2) Fragments Per Kilobase of transcript per Million mapped reads (FPKM) fold change ≥2.0. Function annotation of the DEGs were retrieved from known databases (Yang et al., 2022a). The Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis of all DEGs were performed using the “TBtools” software package (Chen et al., 2020). The cut-off value for GO and KEGG terms were p-value< 0.05. Weighted gene co-expression network analysis was implemented with the WGCNA package in R (Langfelder and Horvath, 2008). Weighted networks with scale-free topology was constructed using a soft threshold power of 9 with R²-value > 0.84. Heatmap and module clustering were used to visualize the co-expression network of all DEGs. Correlations between different modules and conditions (Inoculation vs. Mock) were calculated, selected modules were visualized by cytoscape software (Michael et al., 2010).

Small RNAs in the size range of 18-30 nt were gel purified and ligated to adapters and reverse transcribed to cDNA, then PCR amplification was carried out to enrich the purified fragments. sRNA-seq was carried using SEQ-500 platform (BGI, China). Raw reads obtained from sequencing were processed for quality control and low-quality reads were removed, reads were trimmed and filtered for primer contamination and poly A tail using Trimmomatic v0.35. Reads shorter than 18 nt and longer than 30 nt were discarded. The clean reads were further screened against tRNA, rRNA, snRNA sequences, redundant reads were combined into unique sequences for miRNA prediction. Unique reads were mapped to plant miRNAs from miRbase (https://www.mirbase.org/) to identify known miRNAs. The alignment was done using blastn and these unique reads were also mapped to the ECD04 genome using bowtie2 (Langmead and Salzberg, 2012). For the identified precursor, novel miRNAs were identified using mirDeep2 (https://github.com/rajewsky-lab/mirdeep2). The expression of miRNAs was determined by the R package “DESeq2”. DEMs were identified by 1) adjusted P<0.05, 2) transcript per million mapped (TPM) fold change ≥2.0.

We used only samples at 24 dpi (I24 and M24) for degradome sequencing. The purified cDNA library was sequenced on an Illumina HiSeq 2000 instrument. After removing low-quality data, further mapped to the B. rapa genome to obtain cDNA sense and antisense tags. The tags that mapped to cDNA or mRNA sequences were used for cleavage sites prediction by CleaveLand (Addo-Quaye et al., 2009), categories definition (containing 0–4 categories, 0 for best match) and T-plot figures. psRNATarget was used to predict miRNA targets (Dai et al., 2018).

To validate the RNA-seq and sRNA-seq results, 6 DEGs and 3 DEMs were subjected to qRT-PCR ( Supplementary Table 9 ). Bra-UBC10 and Bra-U6 were used as reference genes, respectively. Total RNA was extracted from the same samples used for RNA-seq or sRNA-seq. mRNA quantification was performed with HiScript II Q RT SuperMix for qPCR Kit (Vazyme Co Ltd, China) and ChamQ Universal SYBR qPCR Master Mix Kit (Vazyme, China). MiRNA quantification was carried out by miRcute Plus miRNA First-Strand cDNA Synthesis Kit (TIANGEN Co Ltd, China) and miRcute Plus miRNA qPCR Detection Kit (TIANGEN, China). CFX384 Touch Real-Time PCR Detection System (Bio-Rad Co Ltd, USA) was used for qRT-PCR.

To generate the Bra-APS4 CRISPR/Cas9 constructs, one 20bp sgRNA targeting two positions, one targeting the first exon of BnaA06G0440400ZS and the other targeting the first exon of BnaC07G0222000ZS, was cloned into the pYLCRISPR/Cas9 expression vector (Ma et al., 2015). These constructs were transferred into Agrobacterium tumefaciens strain GV3101 (Tsingke Co Ltd, China) and transformed to B. napus 409S by hypocotyl transformation. Disease phenotypes were investigated by pot inoculation with PbZj (Zhijiang, Hubei, China) pathogen, using 409R and 409S as resistant and susceptible controls, respectively. 409S and 409R is a pair of NILs, 409R contained clubroot resistance locus from Chinese cabbage material (Li et al., 2021).

ECD04 is a resistant material that showed broad resistance to different races of clubroot pathogen P. brassicae (Yang et al., 2022a). We performed RNA-seq to analyze global transcriptome changes in order to explore its defense response at different stages of P. brassicae infection ( Figure 1 ). We chose three time points (18 dpi, 24 dpi and 30 dpi, dpi=days post infection) and un-infected plants as mock control. RNA-seq samples created 40-43 million raw reads and 39-42 million clean reads with good quality ( Supplementary Table 1 ). On average, 98% of clean data were mapped to the ECD04 reference genome, among which 93% were uniquely matched ( Supplementary Table 1 ). These transcriptome data were used for further analysis.

Differential gene expression was made first pairwise at each time points between the mock group (M) and the inoculation (I) group ( Figures 1B, C ), i.e. M18 vs. I18, M24 vs. I24, and M30 vs. I30. As shown in Figure 1B , 461 (237 up and 224 down), 1186 (857 up and 329 down), and 441 (285 up and 156 down) DEGs (differential expressed genes) were detected at 18 dpi, 24 dpi and 30 dpi, respectively, after overlapping there are 1751 DEGs with threshold of |log2FoldChange| ≥ 1, and Padj< 0.05 ( Figure 1B ; Supplementary Table 2 ). Compared to 18 dpi and 30 dpi, 24 dpi had the most DEG, and more genes were found to be up-regulated then down-regulated ( Figure 1B ). As shown on Venn diagram, a small fraction of DEGs was shared at different time points, 20 genes were differentially expressed at three time points ( Figure 1C ). Moreover, among these 20 common DEGs, 13 were sustainable up-regulated from 18 dpi to 30 dpi, including those involved in SAR (NIMI1, KIC, PME29, WRK, CHI4) and one is down-regulated (PR1) ( Figure 1D ). In particular, there are 3 DEGs that showed different trends, including Beta-galactosidase 1 (BGAL1) with “down-up-down” for 18-24-30 dpi in between mock and inoculation samples, NRX1 with “down-up-up” for 18-24-30 dpi, and NRT26 with “down-down-up”, respectively ( Figure 1D ). Most of these common 20 DEGs are previously shown to be associated with pathogen infection. Overall, these results indicate stage- specific response upon P. brassicae infection in ECD04.

GO enrichment was performed on 1751 DEGs identified for all three time points between mock and inoculation samples ( Supplementary Table 2 ). Since the Biological process is enriched to more terms, we set the value of p-value to 0.01. A total of 68 GO terms were identified in three categories: biological process (BP, 48 terms), molecular function (MF, 15 terms), cellular component (CC, 5 terms) ( Figure 2 ). In particular, the 461 DEGs at 18 dpi were assigned to 16 terms (BP, 1 term; MF, 11 terms; and CC, 4 terms. Likewise, 1186 DEGs at 24 dpi were assigned to 46 terms (BP 35 terms, MF 5 terms, and CC 5 terms). 441 DEGs at 30 dpi were assigned to 52 terms (BP 44 terms, MF 5 terms, and CC 2 terms) ( Figure 2 ).

DEGs from 18 dpi specially enriched in oxidation-reduction process, oxidoreductase activity, and antioxidant activity ( Figure 2 ). 24 dpi always show higher p-value than other times point, like DNA-binding transcription factor activity and membrane, and it specially enriched in calcium ion binding. And 30 dpi was specially enriched in response to biotic stimulus and chitin binding ( Figure 2 ). The common pathways identified for DEGs for all three time points includes the following: DNA-binding transcription factor activity (MF term), membrane (CC term) and extracellular region (CC term) ( Figure 2 ). These results suggest a dynamic change in resistant B. rapa material ECD04 upon P. brassicae infection that involves genes that play crucial roles in clubroot resistance in different pathways.

Defense response involves substantial transcriptional reprogramming in plants. Weighted gene coexpression network analysis (WGCNA) was performed to understand the correlation between the sample expression profile and the treatment. A scale-free WGCNA network was constructed with the soft threshold power of 20 with threshold R² > 0.84 ( Figure 3A ). Out of the total 1751 DEGs, 1283 were assigned into 10 distinct modules ( Figure 3A ; Supplementary Figure 1 ). Each module represents a cluster of highly interconnected genes with similar expression changes.

KEGG pathway analysis were performed on the genes in these modules, we found only “blue”, “brown”, “black” and “turquoise” showed significant enrichment with p-value< 0.05 ( Supplementary Table 3 ). The genes in “blue” module was mostly enriched in glutathione metabolism, the “brown” module genes was mostly enriched in glycolysis/gluconeogenesis ( Supplementary Table 3 ). Based on the correlation values in inoculation samples between 18 dpi, 24 dpi and 30 dpi time points ( Supplementary Figure 1 ), we chose the “turquoise” module for KEGG, GO and gene-interaction network construction ( Figures 3B-E ). The “turquoise” module contained a total of 781 genes, The gene-interaction network constructed by Cytoscape of turquoise module showed hub genes including DTX50 (detoxification, transporter), XTH23 (xyloglucan endoglucosylase/hydrolase), BH019 (Transcription factor), TIFY10A (repressor of jasmonate responses) that might be important for P. brassicae defense response. In turquoise module, the most significantly enriched KEGG pathway was plant-pathogen interaction ( Figure 3B ), and the most significantly enriched GO pathway was DNA-binding transcription factor activity and transcription regulator activity ( Figure 3C ). Furthermore, the interaction network of genes in KEGG term associated with “plant-pathogen interaction” showed additionally two of CDPKS (Calcium-dependent protein kinase 28) and PTI13 (PTI1-like tyrosine-protein kinase 3) during this process. In total, there are 7 gene were brought out after WGCNA and interaction network analysis, they might be the hub genes for clubroot resistance in ECD04 ( Figures 3D, E ).

To investigate the regulation between genes and miRNAs in disease-resistant turnip ECD04, the same set of libraries mentioned above were used for sRNA-seq to capture micro RNA expression changes ( Figure 4 ). Above 21.5 million reads were obtained on average, after comparing with rRNA, snRNA, tRNA, 15 million filtered reads were used for miRNA annotation ( Supplementary Table 1 ). A filtering pipeline was used to identify first the known miRNAs. 157 known miRNA were obtained from the 18 sRNA libraries deposited in miRBase. Besides, 51 potential novel miRNA candidates were identified by mirDeep ( Supplementary Table 4 ). The length of known miRNAs and novel miRNAs ranged from 18-25 nt ( Figure 4A ), the majority of 23 nt long miRNAs were characterized by a 5’-uridine ( Figure 4B ). 5’-uridine has been observed in many plant species and is a characteristic attribute of DCL1 cleavage and AGO1 association (Reinhart et al., 2002; Rajagopalan et al., 2006). On the other hand, the 157 known miRNAs were classified into 60 miRNA families, with miR156 being the largest family with 14 members ( Supplementary Figure 2 ).

According to the TPM values of these miRNAs, we identified DEMs (different expressed miRNAs) at 18, 24 and 30 dpi in between mock and inoculation conditions ( Figure 4C ; Supplementary Table 5 ). DEMs were selected based on p-value< 0.05 and |fold change| > 2. As result, 2, 4 and 8 DEMs were found at 18dpi, 24dpi and 30 dpi, respectively ( Figure 4C ). Further, to interrogate the miRNAs involve in P. brassicae resistance, we looked into DEMs that exhibited similar patterns at all three times points. Indeed, 1, 1 and 4 miRNAs were constantly up-regulated at 18, 24 and 30 dpi, respectively; meanwhile 6 and 11 miRNAs were down-regulated at 24 and 30 dpi, respectively ( Figure 4D ). Notably, bra-miR395a/b/c/d-3p and two novel miRNAs (nov_miR34/35, also members of MIR395) were down-regulated both at 24 dpi and 30 dpi, with more significant downregulation at 24dpi ( Figure 4D ). Overall, these results indicate stage-specific regulation on miRNA expression as well in ECD04 upon P. brassicae infection, possibly leading to differential expression of downstream target genes.

Plant miRNAs usually have perfect or near-perfect complementarity with their targets. To investigate the targets of miRNAs identified above, we used degradome sequencing and generated 22 and 31 million reads for Mock and Inoculation samples, respectively ( Supplementary Tables 6, 7 ). After filtering, 21.6 and 30.5 million reads were obtained, and 108 and 93 miRNA target pairs were identified for mock and inoculation sample, respectively ( Supplementary Tables 6, 7 ). The maximum cleavage sites belonged to category 0 (average 62%) and minimum to category 4 (average 6%). After filtering with category< 2, 90 and 81 miRNA target pairs were left for mock and inoculation conditions. In total, 128 miRNA-target pairs were identified from degradome sequencing ( Supplementary Tables 6, 7 ). Annotation of the target genes suggests most to be disease resistance genes like NBS-LRR, and TFs like auxin response factor (ARF) and scarecrow-like protein (SCL).

As for target genes, 61 targets were detected both in mock and inoculation conditions, whereas 44 and 31 were only for mock or inoculation ( Supplementary Tables 6, 7 ). The maximum number of targets were obtained for miR1885 family, where 20 targets were identified and most of them were NBS-LRR genes ( Supplementary Tables 6, 7 ). The second group of targets was related to plant hormone signal transduction, like bra-miR160/bra-miR167/nov_miR39 and their target ARF, bra-miR171 and its target encoding scarecrow-like protein, bra-miR172 and its target RAP/TOE2/AP1 etc ( Supplementary Tables 6, 7 ). Other miRNA-target pairs involve different metabolic pathways, such as bra-miR164 and target NAC, bra-miR168 with target AGO1 and bra-miR162 with target DCL1 for RNA-mediated gene silencing, and miR395 with its target APS4 (ATP sulfurylase 4) for sulfur metabolism ( Supplementary Tables 6, 7 ).

Combing with DEMs to transcript levels of target genes in response to P. brassicae, only bra-miR395a, bra-miR160a and nov_miR7 were detected in degradome sequencing. Bra-miR395a was down-regulated in 24/30 dpi and its target, BraA06g023661E (APS4) although not DEG, was up-regulated at 24 dpi. The cleavage of Bra-APS4 by bra-miR395a was confirmed by degradome sequencing ( Figures 5B, C ). Bra-miR160a was down-regulated in 30 dpi but its target are not DGEs, among its targets, BraA07g035249E (ARF) and BraA02g030057 (BAG1) were slight up-regulated in 30 dpi. As for nov_miR7, although this miRNA was down-regulated in DEM, but its predicted targets BraA07g034603E and BraA09g004346E (FBX6) could not be detected in RNA sequencing, probably due to too low expression level. Finally, we used psRNATarget server and predicted a total of 2319 miRNA-target pairs involving 208 miRNAs (157 known and 51 novel). KEGG enrichment on the target genes identified terms in Monobactam biosynthesis, Sulfur metabolism and Selenocompound metabolism ( Supplementary Table 8 ). The expression of both P. brassicae responsive miRNAs (from sRNA-seq) and their targets (from RNA-seq) were integrated and inferred 15 critical miRNA-target pairs during P. brassicae infection ( Table 1 , cutoff value p< 0.05 and |log₂FC| > 1). Among these, 1, 5 and 6 pair were identified from 18 dpi, 24 dpi and 30 dpi, respectively ( Table 1 ). Three pairs were identified both for 24 and 30 dpi, which involved potential disease resistance protein (DRL23) and Transcription factor (WRKY51) ( Table 1 , Figure 5A ).

To validation the results from sRNA-seq and RNA-seq for the selected miRNA-target pairs, we selected 5 DEGs and 3 DEMs from the miRNA-target interaction network, and 1 gene (BraA07g033342E) that was shared as DEG for all three time points ( Figure 1C , Table 1 ). The qRT-PCR results were highly consistent with the results from RNA-Seq/sRNA-Seq data for levels of these mRNAs/miRNAs, indicating the inferred miRNA-target interaction to be reliable and might be important during the defense response upon clubroot infection ( Supplementary Figure 3 ).

Based on the correlation analysis of miRNAs expression and degradome sequencing, we identified bra-miR395a to be associated with P. brassicae infection, which indeed the members of the same microRNA family has been reported in previous studies (Li et al., 2021). Bra-miR395a targets BraA06g023661E (APS4), which encodes for a protein for sulfur transfer (Hatzfeld et al., 2000). Bra-APS4 transcript levels were slightly up-regulated at 24 dpi and 30 dpi during P. brassicae infection ( Figure 6A ). Since ECD04 is difficult to transform, we checked its homologs in Brassica napus ZS11, and found the Bna-APS4 has two copies (BnaA06G0440400ZS and BnaC07G0222000ZS). In order to investigate the role of APS4 for clubroot resistance, we used CRISPR/Cas9 to knock out the two copies of Bna-APS4 in a susceptible B. napus material 409S, using one sgRNAs specific for BnaA06G0440400ZS and BnaC07G0222000ZS, respectively ( Figure 6B ). Three individual edited mutants were obtained, all leading to Bna-APS4 loss-of-function ( Supplementary Figure 4 ). When disease phenotypes were examined upon PbZj infection, the DI of resistant and susceptible control and three independent mutant lines were 0, 86, 45, 68, 67, respectively ( Figures 6C, D ). The lowest DI was observed in bna-aps4-1, shown by less gall formation on lateral and main roots at 35 dpi compared to susceptible control 409S ( Figure 6C ). Therefore, knocking out Bna-APS4 in 409S background changed disease phenotype from very susceptible to a moderate susceptible/more resistant phenotype, suggesting a functional Bna-APS4 gene might lead to a less resistance phenotype, i.e. Bna-APS4 may serve as a negative regulator of clubroot resistance in Brassica napus.

In this study, we used RNA-seq, sRNA-seq and degradome to follow the infection procedure of P. brassicae in resistant European Turnip ECD04, in order to pinpoint the resistance mechanism to facilitate molecular breeding using key genes and regulators. We have chosen three time points, 18 dpi, 24 dpi and 30 dpi, with uninfected plants as control samples, and identified 318, 224 and 892 DEGs among 20 were shared. Out of these 20 genes, 14 were consistently up-regulated by clubroot infection, therefore they might be positive regulators that are needed for the development of disease resistance. These genes include transcription factors such as WRKY and NAC, as well as cell wall associated genes, and genes involved in Ca²⁺ influx and signaling ( Figure 7 ). GO enrichment analysis found similar results for the DEG sets at 24 dpi and 30 dpi, which is also supported by the WGCNA analysis using all 1751 DEGs. Further, the turquoise module identified from WGCNA revealed good correlation between gene expression and time after infection ( Supplementary Figure 1 ), KEGG enrichment of this module also suggested plant-pathogen interaction as the most enriched pathway. From WGCNA we identified four hub genes from the turquoise module (BraA03.DTX50, BraA01.XTH23, BraA09.BH019 and BraA08.TI10A), by co-expression analysis associated with pathogen-interaction we obtained an additional three hub genes (BraA07.CDPKS, BraA09.CDPKS and BraA09.PTI13). Expression analysis on these 7 genes showed that they were consistently up-regulated throughout 24-30 dpi ( Supplementary Table 2 ). Functional annotation revealed that BraA01.XTH23 participates on cell wall biosynthesis, BraA09.BH019 is a transcription regulator, BraA08.TI10A functions as a repressor of jasmonate responses (Chung and Howe, 2009), while both BraA07.CDPKS and BraA09.CDPKS are encoding Calcium-dependent protein kinase 28, they are negative immune regulator subjected to ubiquitination and proteasome degradation (Liu et al., 2021). BraA09.PTI13 is protein kinase that interacts with CLV1 and functions in CLE peptide signaling pathway in root development.

Most of the hub genes from WGCNA analysis are associated with PTI response, yet the defense response in ECD04 most likely occur on both PTI and ETI levels, which all starts from the initial recognition of clubroot pathogen by R proteins. Therefore, we expanded the dataset to include DEGs potentially related to immune response, we identified numerous genes involved in Calcium signaling pathway, formation of reactive oxygen species (ROS), plant-pathogen interaction, MAP kinases, Cell wall related genes, plant hormone and Transcription factors as DEGs upon Plasmodiophora infection in ECD04 ( Figure 7 ). Some of these genes have been reported to be associated with plant immune response, e.g. CPK (Liu et al., 2021), RBOH (Wang et al., 2020), EDS1 (Bhattacharjee et al., 2011), MPK3 (Meng and Zhang, 2013). Besides, phytohormones including JA and SA, and Calcium influx as well as between cellular comportments are also reported to influence clubroot resistance (Lemarié et al., 2015; Wang et al., 2023). DEGs associated with cell walls might regulate pathogen penetration by increasing the physical barrier capacity by strengthened cell walls. One of the PR genes, PR1, was consistently induced at all three time points (18 dpi, 24 dpi and 30 dpi) after clubroot infection, suggesting this gene to be of particular importance. Notably, we also identified NLR genes that were induced in later stages of infection, although not as significantly high levels, indeed these genes might be endogenously regulated to present super induction that might be harmful to the host cell itself. These genes serve as valuable source for future functional validation.

During recent years, due to the advance in epigenetics and sequencing technologies, people start to realize that small RNAs play important roles as post-transcriptional regulators and could be used for crop genetic improvement (Navarro, 2016; Tang and Chu, 2017; Raza et al., 2023). In the current study, we used high-through-put sequencing and degradome to investigate changes on both transcript levels and microRNAs. We identified 157 known and 52 novel miRNAs, among which 7 were at least differentially expressed at one time point. Indeed, from 18 dpi to 24 dpi and 30 dpi, most DEGs were found at 24 dpi ( Figure 1B ), while the numbers of DEMs rose continuously from 18-30 dpi ( Figure 4C ). Most of the plant miRNAs regulate target gene expression by transcript cleavage, and the cleavage mostly occur at the 10^th nucleotide within complementary region. Through target gene annotation, we could deduce the function of microRNAs. And the target genes are obtained normally by two approaches: 1) by software prediction, which predict target genes by sequence complementarity with the mature miRNAs; 2) by degradome sequencing, which requires library construction mediated by miRNA-cleaved fragments. Obviously, parameter settings are also important for target gene filtering. In this study we filtered miRNA-target pairs based on category number (0, 1 and 2) in between mock and inoculation conditions. Consequently, we identified 128 miRNA-target pairs. Within these, the DEM known miRNAs showed miR160-ARF and miR395-APS4 ( Figure 7 ), which fits previous study for miR160-ARF involved in plant growth and immune response (Yu et al., 2019), and miR395a-APS4 with sulfite synthesis (Li et al., 2021). For the non-DEMs, we also identified miR164-NAC and miR1885a-NLR, these gene pairs have also been reported with HR response (Lee et al., 2017) and R gene expression (Cui et al., 2020). We also used the psRNATarget online prediction to find those miRNA-target pairs with opposite changes upon clubroot infection at three time points. We obtained 15 pairs of miRNA-targets using this approach, including nov_miR7/25/33-OFUT20, bra-miR395a-LAC7/DRL23, nov_miR34/35-WRKY and nov_miR15-APR3/MD37E ( Figure 5A ). Among those, LAC7 is associated with lignin in cell wall, APR3 is involved in sulfate assimilation, MD37E might play roles in disease response, in addition we also have nov_miR26-HMA8 which participate in copper-transporting. Therefore, sRNA-seq combined with degradome sequencing could help us to find more reliable miRNA-target pairs that could explain trait phenotypes and provide new viewpoint for genetic breeding.

In comparison with previous studies on microRNA and target pairs, we found our result to be similar with the ones found for miR395a-APS4 in B. napus (Hatzfeld et al., 2000), which also showed down-regulation on miR395 upon clubroot infection in both resistant B. rapa or B. napus materials and concurrent up-regulation of APS4. Combined with our data, we hypothesis that miR395a-APS4 might play regulatory role in clubroot resistance in B. rapa ECD04. Our genetic data on CRISPR-edited Bna-APS4 showed that loss-of-function mutants in 409S background showed more resistant disease phenotype to its control (409S). This APS4-CRISPR result is indeed surprising since it is contradictory to our anticipation for APS4 as a positive factor for clubroot resistance when miR395a was down-regulated after inoculation ( Table 1 ) and APS4 is the target of miR395a ( Figure 5 ). We hypothesize that there might be additional mechanisms of APS4 that may influence disease outcome, such as the accumulation of sulfate and other metabolites due to Bna-APS4-loss-of-function, as reported for Xoo infection in rice (Yang et al., 2022b). These hypotheses need to be validated further experimentally. APS4 is a sulfate adenylyltransferase that catalyze the first step of sulfate assimilation, leading to APS (adenosine 5’-phosphosulfate) from sulfate formation (Mugford et al., 2009). On one side, APS could be reduced back to sulfite to form primary metabolites such as glutathione; on the other side, APS could be an active sulfate form to form secondary metabolites such as glucosinolates (GSH) (Mugford et al., 2009, 2010). These studies showed that some sulfur-containing compounds are also components in the plant immune system, they could function as antioxidants, induce defense gene expression, or limit the growth of pathogenic bacteria (Wang et al., 2022a). For example, GSH is highly relevant for plant immune response, in rapeseed glucosinolate level influence clubroot resistance (Rodrigues and Percival, 2019). In rice, APS4 mutant results in accumulation of sulfate and elevated resistance towards Xanthomonas oryzae pv. oryzae (Xoo) and X. oryzae pv. oryzicola (Xoc) (Rodrigues and Percival, 2019; Zamani-Noor et al., 2021; Yang et al., 2022b). In tomato, cysteine as the key precursor of vascular sulfur accumulation has also been shown to be involved in Verticillium dahlia defense (Klug et al., 2014).

We used transcriptome, sRNA-seq and degradome in addition with genetic manipulation of target gene by CRISPR-Cas9 technology, to study the resistant mechanism in ECD04 with superior clubroot resistance. We identified 7 hub genes associated with infection, 15 miRNA-target pairs were predicted and two pairs of particular importance, miR164-NAC and miR395a-APS4. We validated the function of APS4 by mutant disease phenotyping. The clubroot-responding miRNA-target network constructed in this study would facilitate further functional studies for some of the hub genes as well as microRNAs, and provide gene candidate to be utilized for clubroot resistance breeding in cruciferous crops in the future.