Two B. rapa subsp. perviridis cultivars (komatsuna) were selected, one susceptible “Misugi” (Sakata Seed Corporation, Japan) and one resistant “Nanane” (Takii & Co., Ltd., Japan) to A. candida Mibuna isolate WMB01 (Figure 1). Seeds of “Misugi” and “Nanane” cultivars were sown on soil and kept under 16 h light and 8 h dark at 21°C. Seven-day-old plants were inoculated by spraying the A. candida (WMB01) with a concentration of 1 × 10⁵ zoosporangium/ml. Mock inoculation with water was also performed for control samples. To ensure successful inoculation, plants were incubated in a dark growth chamber for 24 h at 22°C with 100% humidity. Then, the plants were moved to a growth chamber and kept under growth conditions of 16 h light and 8 h dark at 21°C, with regular irrigation. For DNA methylation and gene expression studies, one cotyledon of each plant was harvested after 72 h of inoculation, snap-frozen in liquid nitrogen, and stored at −80°C until further use.

Genomic DNA of infected and control samples of susceptible and resistant cultivars was extracted using the CTAB method (Murray and Thompson, 1980). The construction of the bisulfite sequencing library was performed by the BGI Genomics company. In brief, genomic DNA (1,000 ng) was sheared to 200–300 bp using sonication (Covaris^® LE220). DNA end-repair, 3′-dA overhang, ligation of methylated sequence adaptors, and bisulfite treatment were performed using the ZYMO EZ DNA Methylation-Gold kit (ZYMO RESEARCH, USA) according to the manufacturer's instructions. The qualified library was sequenced paired-end (150 bp) on an Illumina HiSeqXTEN System. Two biological replicates were used in this study.

The reads from whole-genome bisulfite sequencing (WGBS) were trimmed using Trimmomatic-0.39 and quality control was performed using FastQC. Trimmed and high-quality reads were mapped to the reference genome of B. rapa V3 (Zhang et al., 2018) using Bowtie2 version 2.2.5 and Bismark-v0-14.3 (Krueger and Andrews, 2011). PCR duplicates were removed as one of the reads, which align to the same position, and are randomly retained by using duplicate_bismark script in the bismark package. To calculate the methylation level of CG, CHG, and CHH contexts, the number of methylated and unmethylated reads were extracted at each cytosine position using a bismark methylation extractor script with a paired-end parameter. The methylation level at each cytosine was calculated by dividing the number of methylated cytosine reads by the total number of reads. A binomial test (Lister et al., 2008) was performed for the classification of methylated and unmethylated cytosine sites. The sequence context-specific error rates of bisulfite conversion were estimated from the mapping results of the unmethylated phage genome and used for the binomial test. The test was applied only to the cytosine sites with greater than or equal to 3 WGBS coverage. A significant cut-off of two-tailed P < 0.01 was used to detect methylated cytosine.

To detect the differentially methylated regions (DMRs), the genome was divided into 500 bp windows with no overlap. The total number of cytosines in each context (CG, CHG, and CHH) in a given window were calculated and based on the distribution frequency, and windows containing ≥ 10 cytosines in all CG, CHG, and CHH contexts were kept for downstream analysis. The methylation level of CG, CHG, and CHH was calculated in each window by dividing the number of methylated cytosines in each given context by the total number of cytosines in the given context. To identify the DMRs in each cultivar between control and infected samples, an F_ST-like approach was performed. F_ST statistics are widely used in population genetics to measure the level of population differentiation (Holsinger and Weir, 2009). Here, first, highly differentiated windows between the biological replicates were removed and the remaining windows were used for F_ST calculation, where the methylation level of a certain window was denoted as X1: methylation level of a given cultivar under control condition for replicate 1, X2: methylation level of the given cultivar under control condition for replicate 2, X3: methylation level of the given cultivar after infection with A. candida replicate 1 and X4: methylation level of the given cultivar after infection with A. candida for replicate 2; the variance between replicates (V) and absolute methylation difference between treatments (δ) was then calculated as

where the max (X, Y) function is to take X, if X > Y, otherwise to take Y. Using V and δ value, the relative methylation difference in the given window was calculated as

F_ST was calculated for each methylation context (i.e., CG, CHG, and CHH) separately in both resistant and susceptible cultivars. To detect the regions with maximum methylation difference between control and infected samples while ensuring the minimum differences between biological replicates only, windows with top 1% F_ST and top 1% δ values were assigned as putative DMRs and used for downstream analysis. DMRs were associated with their proximal genes where DMRs were located at the gene body (i.e., between the start and stop codon including all exons and introns of a gene), and/or located up to 2 kb upstream and downstream of genes. Then, DMR-associated genes were used for gene ontology (GO) enrichment analysis. GO enrichment analysis was performed using the R package topGO (v 2.40.0). All the heatmaps for comparative analysis of DNA methylation levels were constructed using the R package ComplexHeatmap (V 2.2.0).

To verify the WGBS results, we picked four regions, including hyper- and hypo-methylated DMRs in both cultivars. A total of 500 ng of DNA was fragmented by sonication and the fragments were ~300–800 bp in length. MethylCode Bisulfite Conversion Kit (Thermo Fisher Scientific, Inc, USA) was used for chemical bisulfite reaction and PCR was performed using bisulfite-treated DNAs as templates. PCR conditions were 95°C for 2 min followed by 40 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s. Amplified PCR fragments were gel-purified using GENECLEAN III Kit (MP Biomedicals, USA) and cloned into pGEM-T Easy vector (Promega Co., USA). Ten independent clones were sequenced. Primers used for bisulfite sequencing are listed in Supplementary Table 1.

Total RNAs of infected (72 h after A. candida inoculation) and control samples of resistant and susceptible cultivars were extracted using SV Total RNA Isolation System (Promega). RNA sequencing was performed using paired-end Nextseq500 (75 bp read length). The number of clean reads and the percentage of mapped reads are shown in Supplementary Table 2. Low-quality reads were filtered using FASTX-Toolkit v. 1.4.5 and HISAT2 was used to align the filtered reads to the reference genome of B. rapa V3 (Zhang et al., 2018). The expression levels (fragments per kilo-base per million—FPKM) were scored using Cuffdiff. Differentially expressed genes were also identified based on two criteria of two-fold difference (|log 2 ratio| ≥ 1.0) and 95% confidence. Differentially expressed genes located close to DMRs were also detected and their protein sequences were searched against Pfam, SMART, and PRINTS databases using InterProScan (https://www.ebi.ac.uk/interpro/) for domain identification.

Seven genes were used for the validation of RNA-seq results. cDNA was synthesized from 500 ng total RNA using the ReverTra Ace^® qPCR RT Master Mix with gDNA Remover (TOYOBO Co., Ltd., Japan). The specificity of the primer set of each gene was first tested by electrophoresis of RT-PCR amplified products using QuickTaq^®HS DyeMix (TOYOBO) on 1.5% agarose gel in which single products were observed. RT-PCR conditions were 94°C for 2 min followed by 35 cycles of 94°C for 30 s, 55°C for 30 s, and 68°C for 30 s. The absence of genomic DNA contamination was confirmed by the PCR of no RT control. Real-time RT-PCR (qPCR) was performed using a LightCycler 96 (Roche Molecular Systems, Inc., USA). cDNA was amplified using FastStart Essential DNA Green Master (Roche). qPCR conditions were 95°C for 10 min followed by 40 cycles of 95°C for 10 s, 60°C for 10 s, and 72°C for 10 s, and a melting program of 65–97°C at 0.1°C/s. After amplification cycles, each reaction was subjected to melt temperature analysis to confirm the presence of single amplified products. The relative expression level of each gene relative to ACTIN (Bractin) was automatically calculated using automatic CQ calling according to the manufacturer's instructions (Roche) (Fujimoto et al., 2006a). The data presented are the average and standard error of three biological and experimental replicates. The primer sets are listed in Supplementary Table 1.

WGBS was performed to detect the single-based resolution and relative amount of 5-methylcytosines (5-mCs) changes across the genome of susceptible (“Misugi”) and resistant (“Nanane”) cultivars of B. rapa in response to A. candida. The high-quality reads were mapped to the B. rapa reference genome (Zhang et al., 2018) (Table 1). The mapping efficiency of all samples was between 35.9 and 43.7% (Table 1). The bisulfite conversion error of CG, CHG, and CHH of all samples was between 0.002135 and 0.004833 (Table 1). In both “Misugi” and “Nanane,” DNA methylation occurs predominantly at the CG sites, followed by CHG and CHH, ranging between 40 and 51%, 19 and 27%, and 5 and 10% of sites being methylated, respectively (Figure 2).

DMRs of each sequence context (CG, CHG, and CHH) were detected between control and infected samples of “Misugi” and “Nanane.” Regions were assigned as putative DMRs where maximum methylation difference occurs between control and infected samples while having the minimum differences between biological replicates. In total, 233 and 275 DMRs were detected in the susceptible cultivar “Misugi” and the resistant cultivar “Nanane,” respectively (Figure 3A and Supplementary Table 3). No overlapping DMRs were found in all three cytosine contexts between “Misugi” and “Nanane.” We defined the heterochromatic region as having more than 0.4 of the density of TEs. Approximately 30% of DMRs were found in heterochromatic regions in both “Misugi” and “Nanane,” and more than 40% of CHG DMRs were in heterochromatic regions (Supplementary Table 3).

In both cultivars, despite the higher percentage of methylated cytosine at CG sites through the genome, the number of CG DMRs (“Misugi:” 20, “Nanane:” 28) was noticeably lower than DMRs for CHG (“Misugi:” 120, “Nanane:” 129) and CHH (“Misugi:” 93, “Nanane:” 118) (Figure 3A and Supplementary Table 3).

The comparative analysis of the methylation level and status of DMRs (i.e., hypo- and hyper-methylation) showed in both cultivars, more than half of DMRs were hyper-methylated, 53.64 and 62.90% in “Misugi” and “Nanane,” respectively. In the susceptible cultivar “Misugi,” most of the CG (80.00%) and CHH (69.89%) DMRs were hyper-methylated, while the majority of CHG DMRs (63.33%) were hypo-methylated. In the resistant cultivar “Nanane,” half of the CG DMRs (50.00%) and CHH (52.54%) were hypo- and hyper-methylated and CHG (75.19%) was predominantly hyper-methylated after infection with the pathogen (Figure 3B, Supplementary Figure 1, and Supplementary Table 3). Bisulfite sequencing confirmed the reliability of the performed WGBS (Supplementary Figure 2).

Genes associated with DMRs were screened where DMRs were located in gene bodies (i.e., from the start codon to stop codon) and/or up to 2 kb up/downstream of genes. Out of the 233 DMRs in the susceptible cultivar “Misugi,” 126 (54.07%) were proximate to 136 genes (CG: 21, CHG: 56, CHH: 59) and in the resistant cultivar “Nanane,” out of 275 DMRs, 155 (56.36%) were proximate to 178 genes (CG: 29, CHG: 72, CHH: 77) (Figure 4, Table 2, and Supplementary Table 3). In both cultivars, the majority of genes were linked to the CHH DMRs (Figure 4).

We then identified DMR locations relative to genes. In “Misugi” (susceptible cultivar), 61 out of 136 DMRs (44.85%) were located in gene bodies, which also included the highest number of all three DMRs' contexts (CG: 17, CHG: 20, CHH: 24). In “Nanane” (resistant cultivar), 65 out of 178 (36.51%) DMRs were located upstream of genes with the highest number of CHH (32) and CHG (25) DMRs and 62 out of 178 (34.83%) DMRs were located downstream of genes with the highest number of CHG [29 (40.27%)] and CHH [29 (37.66%)] DMRs (Figure 5A and Table 2). Despite the difference in CHG and CHH DMRs distribution between the two cultivars, in both cultivars, most of the CG DMRs are located in gene bodies (“Misugi:” 80.95%, “Nanane:” 58.62%) (Figure 5A and Table 2).

In both “Misugi” and “Nanane,” 59.52 and 58.06% of DMRs associated with genes were hyper-methylated, respectively. In both cultivars, the hyper-methylation was predominant in DMRs within the gene bodies and downstream of genes; however, in the upstream of genes, the majority of DMRs were hyper-methylated in “Misugi” (16.17%) and hypo-methylated (19.10%) in “Nanane” (Figure 5B and Table 2).

In “Misugi,” CG DMR-genes tended to be hyper-methylated, dominantly located within the gene body (66.66%). No hypo-methylated CG DMRs were identified in the downstream of the genes (Figure 5C and Table 2). CHH DMRs tended to be hyper-methylated at the upstream and downstream of the genes (Figure 5C and Table 2). In “Nanane,” CG DMRs tended to be hypo-methylated at the upstream and downstream of the genes, while CG DMRs tended to be hyper-methylated within the gene body (Figure 5C and Table 2). In addition, CHG and CHH DMRs were dominantly hyper-methylated in the gene body (Figure 5C and Table 2). In the upstream region, CHG DMRs tended to be hyper-methylated, while CHH DMRs tended to be hypo-methylated (Figure 5C and Table 2).

GO enrichment analysis of genes associated with DMRs indicates the involvement of genes in all three GO categories, that is, biological process, molecular function, and cellular component (Figure 6 and Supplementary Table 4). DMR-associated genes in “Misugi” were highly enriched for biological process and and molecular function. The top two highly enriched (P < 0.05) classes for the biological process were DNA integration (GO: 0015074) and RNA-dependent DNA biosynthetic process (GO: 0006278) and top two of highly enriched (P < 0.05) terms for molecular function were aspartic-type endopeptidase activity (GO: 0004190) and RNA-directed DNA polymerase activity (GO: 0003964) (Figure 6 and Supplementary Table 4). In “Nanane,” DMRs-associated genes were highly enriched for the plastid membrane (GO: 0042170) from the cellular component category and for the sieve element enucleation (GO: 0090602) from the biological process category (Figure 6 and Supplementary Table 4). DMR genes were also compared against B. rapa resistance gene analogous (RGAs), which were reported by Tirnaz et al. (2020a), and the results showed that DMR genes in “Misugi” and “Nanane” included two receptor-like kinases (RLKs) genes (Supplementary Table 3).

DNA methylation is one of the main regulatory mechanisms controlling TEs (Frost et al., 2005; Fujimoto et al., 2008b; Sasaki et al., 2011), movements, and activity throughout the genome (Tirnaz and Batley, 2019a). In this regard, here we investigated the occurrence of DMRs within TEs. In the susceptible cultivar (“Misugi”), 117 out of 233 DMRs (50.21%) were associated with TEs (DMR-TEs), and in “Nanane,” 129 out of 275 DMRs (46.90%) were associated with TEs (Supplementary Table 3). About 40% of DMRs associated with TEs were located in the heterochromatic regions in “Misugi” and “Nanane” (Supplementary Table 3). Consistent with other studies (Cokus et al., 2008; Li et al., 2012; Regulski et al., 2013; Song et al., 2013; Wang et al., 2015), we found TE methylation in all three contexts (Supplementary Figure 3). From total DMRs in each context in “Misugi,” Ten percentage of CG DMRs, 52.50% of CHG DMRs, and 55.91% of CHH DMRs were associated with TEs. Similarly, in “Nanane,” a low percentage of CG DMRs (21.42%) were associated with TEs. However, as opposed to “Misugi,” CHG DMRs (55.81%) showed a higher percentage of TE association than CHH DMRs (43.22%). In addition, the results revealed that in “Misugi,” 100% of CG-DMR associated with TEs were hyper-methylated, while the majority of CG-DMR-associated TEs in “Nanane” (83.33%) were hypo-methylated (Figure 7). DMRs in CHG context also showed the opposite trend of methylation between the two cultivars, where in “Misugi,” the large number of CHG-DMR associated with TEs were hypo-methylated (69.84%) and in “Nanane,” they were hyper-methylated (76.38%) (Figure 7). In both cultivars, the higher percentage of CHH-DMR associated with TEs were hyper-methylated (Figure 7). We also found 17 DMRs (CG: 1, CHG: 8, CHH: 8) and 30 DMRs (CG: 3, CHG: 13, CHH: 14) in association with TEs located at up to 2 kb upstream of genes in “Misugi” and “Nanane,” respectively (Supplementary Table 3). In addition, 20 DMRs (CG: 1, CHG: 7, CHH: 12) and 16 DMRs (CG: 2, CHG: 7, CHH: 7) in association with TEs located within the genes in “Misugi” and “Nanane,” respectively (Supplementary Table 3). About 35 and 18% of DMR-TEs were Gypsy-type in “Misugi” and “Nanane,” respectively, and half of CHH-DMR associated with TEs were Gypsy-type in “Misugi” (Supplementary Table 5).

Here, RNA-seq analysis on control and infected cotyledons (72 h after A. candida inoculation) of “Misugi” and “Nanane” cultivars were performed. The results showed that, in total, 4,104 and 2,157 genes were differentially expressed in “Misugi” (down: 1,979, up: 2,125) and “Nanane” (down: 1,112, up: 1,045), respectively (Supplementary Table 6). In total, seven genes were selected and the expression levels of these genes in control and infected samples of “Misugi” and “Nanane” cultivars were examined by qPCR. There was a high correlation (r = 0.97 in “Misugi,” r = 0.96 in “Nanane”) of the ratio of expression levels between the control and infected samples of these seven genes observed between qPCR analysis and RNA-seq data (Supplementary Table 7). The association study of differentially expressed genes with DMRs showed that, in total, 13 genes were associated with DMRs in “Misugi” (8 genes) and “Nanane” (5 genes). The methylation context of DMRs between the two cultivars appeared to be opposite, while CG and CHH contexts were predominated in “Misugi” and CHG context was predominated in “Nanane” (Table 3). Four (1 in “Nanane” and 3 in “Misugi”) differentially expressed genes overlapped with DMR-TEs (Supplementary Table 3).

In total, 22 domains were identified across 10 genes and no domains have been identified for 3 genes (Supplementary Table 8). Some of the domains are reported to be involved in defense mechanisms such as flavin-dependent monooxygenases (FMOs) (Mishina and Zeier, 2006; Thodberg and Jakobsen Neilson, 2020) and type III polyketide synthase-like protein (PKS) (Tanjung et al., 2020). Similarly, genes containing haloacid dehalogenase-like hydrolase and Proline-rich domains were also expressed in sweet orange in response to citrus blight (Fu et al., 2019). In addition, the Arabidopsis thaliana orthologs of these 13 genes have been identified, and 8 genes out of 13 had orthologs in A. thaliana (Supplementary Table 8). A. thaliana orthologs, including AT1G55690 (Sec14p-like phosphatidylinositol transfer family protein), AT1G03870 (FASCICLIN-LIKE ARABINOOGALACTAN 9, FLA9), and AT1G24510 (CCT5, CHAPERONIN CONTAINING T-COMPLEX POLYPEPTIDE-1 SUBUNIT 5), have multiple roles in plant growth, development, and signaling pathways (Kamal et al., 2020; Wu et al., 2020), but their involvement in response to pathogens is not well-understood.