In the current study, F₁ and F₂ populations were developed by crossing the clubroot-resistant turnip inbred line BrT18-6-4-3 (P₁) with the susceptible Chinese cabbage DH line Y510 (P₂) followed by the selfing of F_1. Parental lines, F₁, and F₂ populations, were inoculated with P. brassicae race four isolate ‘XY-2’ (Yuan et al., 2017). The F₂ population was cultivated between September 2021 and November 2021, and the inheritance pattern of the resistant gene was determined using a Chi-square test (χ ²). Furthermore, 26 disease-resistant and 29 susceptible lines were used to analyze mutations in the candidate gene ( Supplementary Table 1 ). All materials used in this study were provided by the Institute of Vegetables, Henan Academy of Agricultural Sciences.

The disease index was assessed by inoculating the plants on the 20th day after sowing using a root irrigation method with a concentration of 10⁷ spores/mL of P. brassicae suspension (Luo et al., 2014). AKIMEKI and ECD05 were utilized for disease index estimation as controls for disease resistance and susceptibility, respectively (Matsumoto et al., 2012). The inoculated plants were maintained at 22°C ± 2°C and 16-h/8-h (light/dark) photoperiod (Yuan et al., 2021). Samples were collected and immediately frozen in liquid nitrogen and then stored at −80°C for subsequent use. Root samples from resistant and susceptible materials were collected at different time points: 0 days after inoculation (DAI; indicating no infection), 3 DAI (representing cortical infection), 9 DAI (indicating early onset of disease), and 20 DAI (indicating late onset of disease). One-centimeter-long root segments were cut from the junction of the rhizome and fixed with Formalin-Aceto-Alcohol (FAA) solution (containing 5 mL of acetic acid + 5 mL of formalin + 50 mL of 95% alcohol + 35 mL of sterile water) for 24 h. After fixation with FAA, the roots were rinsed with sterile water, dried with absorbent paper, and stained with 0.5% fluorescent pink for 3 h. The infected roots were observed under a light microscope (Li et al., 2022).

DNA was extracted from resistant and susceptible parents, along with 30 highly resistant (DI = 0; R-pool) and susceptible F₂ plants (DI = 5–7; S-pool) ( Supplementary Table 1 ) following the cetyl trimethyl ammonium bromide method. Bulked segregant analysis (BSA) sequencing (BSA-seq) was performed on four DNA pools (Huang et al., 2017), and sequencing was conducted using the Illumina HiSeq platform, with analysis carried out using NovoGene (http://www.novogene.com). The raw reads underwent quality control trimming to remove low-quality paired reads with the adapter and retained clean reads. The clean reads from the four DNA bulks were combined using BWA software (V0.7.17) and compared with B. rapa (Chiifu-401) genome version 3.0 (http://brassicadb.org/brad/index.phpM) (Kathiresan et al., 2014). SAMtools software (V1.3.1) was employed to identify SNPs and insertion/deletion (InDel) variations between the R and S pools. The ΔSNP index for all genomic positions in the R and S libraries was calculated ( Supplementary Table 2 ) using a sliding window analysis with a window width of 1 Mb and a sliding window step size of 10 kb.

Sequences of the genes were downloaded from the Brassica database (http://brassicadb.org/brad/downloadOverview.php), and the SNPs of candidate genes were validated. Thirty-five pairs of primers for KASP markers were designed using DNAMAN software (V6.0) to verify the parents and F₂ population ( Supplementary Table 3 ), and, ultimately, 21 pairs of polymorphic primers were selected (Xiang et al., 2007).

The KASP markers were utilized for genotyping the F₂ population consisting of 98 individuals. Linked loci were identified using JoinMap4.0 to form linkage groups (Voorrips, 2002), and a genetic linkage map was constructed according to the Kosambi mapping function (Kosambi, 1944).

An additional large F₂ population of 1,489 individuals was planted for recombination selection, from which 438 susceptible individuals were selected and used to map candidate genes.

DNA and cDNA sequences of CRA8.1.6 candidate genes were cloned from the parents using Phanta Super-Fidelity DNA polymerase Mix (Vazyme, Nanjing, China) ( Supplementary Table 4 ). Standard PCR was conducted with 30-s denaturation at 98°C followed by 35 cycles of 98°C for 10 s, 58°C for 5 s, and 72°C for 60 s.

PCR products were sequenced by Sangon Biotech Co., Ltd. (Zhengzhou, China) and analyzed using DNAMAN software (V6.0). The complete coding sequence (CDS) of the CRA8.1.6 gene from BrT18-6-4-3 and Y510 has been submitted to GenBank under the accession number AB605024.1. Protein’s physical and chemical properties were determined with Expasy (http://web.expasy.org/protparam/).

Total RNA was extracted from root samples collected at 0 DAI, 3 DAI, 9 DAI, and 20 DAI of BrT18-6-4-3 and Y510 lines using the RNA Prep Pure Plant Kit (Beijing Day, China) (Wei et al., 2021). Single-stranded cDNA was synthesized using Trans-script after a two-step removal of genomic DNA (gDNA) (Trans, Beijing, China).

The quantitative real-time PCR (qRT-PCR) was performed in a Roche Light Cycler 480-II system (Roche Applied Sciences, Beijing, China). GAPDH was used as an internal reference gene to calculate relative expression levels using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). B. rapa sequence information was utilized to design RT-PCR primers. Standard PCR was conducted by denaturing for 5 min at 94°C, followed by 30 s at 94°C, 30 s at 55°C, and 60 s at 72°C for 25 cycles. PCR products were visualized after electrophoresis on a 1% agarose gel ( Supplementary Table 4 ).

Specific full-length primers were designed to incorporate restriction sites and protective bases. This was achieved using the plasmid DNA of the target gene as a template. The PCR amplification procedure was conducted as follows: 95°C for 5 min; 30 cycles of 95°C for 30 s, 50°C for 45 s, and 72°C for 258 s; followed by 72°C for 10 min; and 16°C for 30 min.

The purified PCR product was digested with pBWA (V) HS expression vector (Bio Run, Hu Bei, China), and the digested product was further purified and ligated with the vector using T4 ligase to obtain the recombinant plasmid pBWA (V) HS-BrT18-6. The ligation product was then transformed into competent DH5α cells and PCR-screened positive clones. After sequencing and verification, Arabidopsis thaliana was transformed using dip method (Hamann et al., 2002). The resulting T₀ transgenic A. thaliana was identified (Pang et al., 2022), and positive plants were selected for future planting to obtain T₂ generation seeds.

The T₂ generation seeds of A. thaliana were sown on a medium containing hygromycin to select homozygous strains based on positive reactions on specific media (Tang, 2007). T₂ generation A. thaliana seeds were planted (Suwabe et al., 2006). The P. brassicae strain used in this study was obtained from a clubroot-infected Chinese cabbage field (B. rapa) in Xinye County, Henan Province, China (113.97°E, 35.05°N), and the Williams system was used to confirm the P. brassicae strain as race 4 (Yuan et al., 2017).

Twenty-five-day-old transgenic and wild-type A. thaliana were inoculated with P. brassicae. Similarly, another set was inoculated with fresh water without any P. brassicae spores, which served as the control group. After 30 days of inoculation, the wild-type and transgenic A. thaliana were examined in control and experimental groups. The frequency of clubroot disease and expression of candidate genes were observed (Keita et al., 2012).

After 30 days of infection with P. brassicae, all F₁ individuals resisted clubroot disease. We inoculated two fractions of the F₂ population with P. brassicae to account for the nature of disease-resistance gene action. In a small fraction of the F₂ population, 341 plants were resistant and 132 were susceptible, whereas, in a large fraction, 1,134 plants were resistant and 335 were susceptible. Both F₂ populations exhibited a 3:1 segregation, indicating Mendelian genetics for a single dominant gene ( Table 1 ).

Microscopic examination revealed abundant spores attached to root hairs in all test plants at 1 DAI, indicating sufficient spores to initiate root infection. However, no cortical infection was detected in the resistant line (BrT18-6-4-3) at 1 DAI, 2 DAI, and 3 DAI or in the control group. In contrast, root hair infection was observed in the susceptible line (Y510) at 1 DAI, with cortical infection and free spores in epidermal cells observed at 3 DAI. No such symptoms were observed in the susceptible lines of the control group.

At 9 DAI, the roots of the resistant line appeared normal, whereas the roots of the susceptible line were slightly swollen. Numerous spores were found in the swollen root cortical cells of the susceptible line. In contrast, only a few spores were observed in the cortical cells of the resistant line ( Figure 1 ), indicating cortical infection occurred between 3 DAI and 9 DAI.

BSA-seq analysis generated 36,416,118 and 38,958,860 clean reads from the R and S pools, respectively ( Supplementary Table 2 ). These reads were aligned to the B. rapa genome, identifying 218,150 SNPs and 53,040 InDels between the R and S pools. △ (SNP index) for each locus was calculated using sliding window analysis. Subsequently, a new quantitative trait locus (QTL) called CRA8.1.6 was identified on chromosome A08, likely between 11.04 Mb and 16.10 Mb ( Supplementary Figure 1 ).

Thirty-five KASP markers were created on the basis of SNP variants within the candidate interval and polymorphism between the resistant (BrT18-6-4-3) and susceptible (Y510) lines, with 21 markers selected for genotyping of 98 F₂ plants to develop markers for chromosome A08. Within the reported CRA8.1.6 QTL, one recombinant individual was found each for markers L-CR10 and L-CA01 with a genetic distance of 0.9 cM and 1.2 cM, respectively ( Figure 2A ). The order of markers on the genetic map was consistent with the physical map.

To refine the CRA8.1.6 QTL, 438 susceptible homozygous F₂ individuals were screened for polymorphisms between markers L-CR03 and L-CA22, identifying 36 recombinants. All 51 recombinants were further genotyped using markers L-CR10 and L-CA01, identifying nine recombinants. Five KASP markers were developed, two of which showed polymorphism. Ultimately, a candidate gene for the CRA8.1.6 QTL was determined at a position (12,255–12,365 Mb) on chromosome A08 marked by markers L-CR11 and L-CR12 ( Figure 2B ).

Detailed mapping of the CRA8.1.6 QTL region revealed 15 genes, with BraA08g015220.3.5C being the sole gene encoding the TIR-NBS-LRR domains ( Table 2 ). We utilized qRT-PCR to examine the expression of the hypothesized candidate gene. Among the 15 genes, only BraA08g015220.3.5C exhibited differential expression between the parental lines ( Figure 3A ). The expression of the BraA08g015220.3.5C gene was analyzed across various periods. The BraA08g015220.3.5C gene exhibited a four-fold increase in expression in the inoculated resistant line (BrT18-6-4-3) compared to the non-inoculated resistant line and both the inoculated and non-inoculated susceptible (Y510) lines at 3 DAI. Similarly, at 9 DAI, there was a two-fold increase in expression in the inoculated resistant line compared to the other lines ( Figure 3B ). The trend of gene expression analyzed by RT-PCR was consistent with that of qPCR ( Figure 3C ). Hence, the BraA08g015220.3.5C gene was identified as a potential candidate gene responsible for imparting resistance against clubroot disease.

PCR products of the BraA08g015220.3.5C gene were utilized to determine gDNA and CDS. Specific primers amplify disease-resistant (BrT18-6-4-3) and susceptible (Y510) lines. In the resistant line, the genome and CDS sequences of BraA08g015220.3.5C were 4295 bp and 3675 bp, respectively. Conversely, in the susceptible parent, the genome and CDS sequences were 9586 bp and 3644 bp, respectively ( Supplementary Figure 2 ). Conserved domain analysis unveiled several domains within the BraA08g015220.3.5C gene, including a Toll and IL-1 domain (TIR, amino acids 70–235), a nucleotide binding site (NBS-ARC, amino acids 260–484), and an LRR (amino acids 701–908). Compared to disease-resistant and susceptible materials, we identified 249 SNPs, seven insertions, and six deletions in the susceptible line. Among these mutations, the insertion of the long terminal repeat (LTR) retrotransposon does not influence the CRA8.1.6 function, as it does not alter the CDS sequence ( Figure 4A , Supplementary Figure 2 ). The sequencing results revealed six non-synonymous mutations that did not influence the resistance of the wild, because these mutations were contained in both resistant and susceptible materials ( Supplementary Figures 3 , 4 ). Furthermore, sequence alignment analysis revealed that the fourth exon of the susceptible line harbored two deletions and an insertion, leading to premature termination of the C-terminus and destruction of the LRR domain ( Figure 4B ). This premature termination consequently led to the loss of CRA8.1.6 function, rendering Y510 susceptible to clubroot disease.

Using the sequence variation observed in the fourth exon between the parental lines, we designed upstream and downstream primers for 2,855 bp and 3,395 bp, respectively. Functional markers, including CRA08-InDeL ( Supplementary Table 4 ), were also developed to distinguish individuals of BrT18-6-4-3, Y510, and F₂ populations. The amplified fragment size was 541 bp in the resistant materials and 447 bp in the susceptible materials. The results indicated complete concordance between the phenotypes of resistant, susceptible, and heterozygous F₂ plants with their respective genotypes. This suggests that CRA08-InDel is a functional marker capable of accurately distinguishing clubroot-resistance phenotypes within segregating populations ( Figure 5A ). Within a natural population, the genotypes of disease-resistant materials consistently matched their phenotypes. However, the identification of susceptible materials was not sufficiently consistent. Hence, we sequenced materials that did not match this marker and identified a trait-associated SNP at 3,209 bp. Subsequently, we developed the KASP functional marker CRA08-KASP1 ( Supplementary Table 4 ). Upon screening F₂ populations with this marker, we observed that the expected amplified products for resistant, susceptible, and heterozygous F₂ individuals perfectly matched the anticipated results. This suggests that the marker CRA08-KASP1 segregated in the F₂ population with a clubroot-resistant phenotype ( Figure 5B ). The CRA08-KASP1 markers were also screened for 24 disease-resistant and 26 disease-susceptible materials. Among these, the phenotype and genotype of all 23 resistant materials were matched entirely, whereas the phenotype and genotype of 22 susceptible materials were matched. However, only five materials (DH40, BrT127, BrT47-1-1-1, BrT133, and BrT71-1-1-4-3) did not correspond with the genotype ( Figure 5C ). In our future research, we will clone the genes of these five non-matching materials identified by the CRA08-KASP1 marker to explore other variations that might cause clubroot disease.

PCR was conducted using specific full-length primers to amplify the target gene from the subclone vector, which has a length of 4,000 bp. The correct cloning vector was digested using XhoI + EcoRI. Subsequently, the exogenous fragment was ligated into the pBWA (V) HS vector using T4 DNA ligase to obtain the recombinant plasmid ( Supplementary Figure 4 ). The recombinant plasmid was confirmed to be correct, and the overexpression vector was successfully constructed. Overexpression Arabidopsis lines were identified in hygromycin-resistant medium, and homozygous lines were obtained through consecutive selfing over two generations. DNA extracted from T₁ transgenic A. thaliana was identified ( Supplementary Figure 5 ). The recombinant plasmid was the positive control, whereas water was the negative control. The PCR bands observed in the eight transgenic plants were consistent with the bands of the positive control. Consequently, eight T₂ lines (CR-1, CR-2, CR-3, CR-6, CR-7, CR-8, CR-11, and CR-12) were identified as overexpressing the target gene ( Figures 6A, B ).

Transgenic lines CR-6, CR-7, and CR-12, along with wild-type Arabidopsis, were chosen for clubroot phenotypic identification. One set of 25-day-old Arabidopsis plants was inoculated with P. brassicae spores, whereas the other set was inoculated with water only as an experimental control. The experiments were performed with three biological replicates.

After 30 days of inoculation, the control wild-type and transgenic Arabidopsis plants displayed no disease symptoms. In the experimental group, within the overexpression line, CR-6, 22, 5, and 7 out of the 34 plants displayed disease grades 0, 1, and 3, respectively, resulting in a disease index of 10.92%. For the CR-7 line, 19, 5, and 6 out of the 30 plants exhibited disease grades 0, 1, and 3, respectively, with a disease index of 10.95%. Regarding the CR-12 overexpression line, 20, 3, and 9 out of the 32 plants showed disease grades 0, 1, and 3, respectively, yielding a disease index of 13.39% ( Supplementary Table 6 ). In wild-type Arabidopsis, the frequency of disease stages in CR-6, CR-7, and CR-12 was 5 or 7, with disease rates of 85.71%, 89.01%, and 79.59%, respectively ( Figure 6C ).

The disease index in the overexpression transgenic Arabidopsis was significantly lower than that of the wild-type Arabidopsis in the experimental group ( Figure 6D ). RNA was extracted from transgenic Arabidopsis plants and reverse transcribed into cDNA to analyze the expression of candidate genes. Criterion 1 was used to evaluate the expression level of the candidate genes in transgenic Arabidopsis based on their expression levels in the control group. The relative expression of the candidate gene in the wild A. thaliana from the experimental and control groups was only 0.04. In contrast, gene expression levels were over 40 in transgenic lines, significantly higher than those in the control ( Figure 6E ).