Loss of CG methylation enables an augmented defense response in weakly primed plants

We first confirmed the function of DNA methylation in the defense response by testing the growth of a bacterial pathogen Pseudomonas syringae pathovar (pv.) tomato strain DC3000 (DC) in wild-type plants (Col-0 and Ler-0) and DNA methylation mutants (Fig. 1A). The DNA demethylase mutants had increased DNA methylation and compromised resistance (Fig. 1A, white bar; ros1 and ros1 dml2 dml3 (rdd) [26, 27, 44]), whereas the met1 mutation, which decreases CG, significantly enhanced resistance (Fig. 1A, white bar; met1). Consistent with a previous report [23], CHH demethylation (drm1 drm2 plants) and nonCG (CHG and CHH) demethylation (drm1 drm2 cmt3 plants (ddc) and drm1 drm2 cmt2 cmt3 plants (ddcc)) suppressed pathogen growth (Fig. 1A, white bar; drm1 drm2, ddc). Specific loss of CHG methylation did not affect resistance (Fig. 1A, white bar; cmt3). These results suggest that global CG hypomethylation (met1) or CHH hypomethylation (met1, drm1 drm2, ddc, and ddcc) enhances resistance to DC, whereas loss of CHG methylation (cmt3) has a limited effect.

Fig. 1

Augmented defense response against the biotrophic pathogen in CG-hypomethylated mutants. A Bacterial growth of Pseudomonas syringae pv. tomato DC3000 (DC) in DNA methylation mutants with or without sub-optimal priming. For sub-optimal priming, water (mock control) or β-amino butyric acid (BAsub; 10 ml of 30 μg/ml per plant) was applied near roots 2 days prior pathogen inoculation. Bacterial growth was measured 3 days post-inoculation (dpi) from at least five biological replicates (n ≧ 5). ddc, ddcc, and rdd indicate drm1 drm2 cmt3, drm1 drm2 cmt2 cmt3, and ros1 dml2 dml3 plants, respectively. p-values by Student’s t test. cfu, colony-forming units. B Expression of WRKY genes in Col-0, ddm1, and met1 plants before (0 dpi; mock and BAsub) and after (1 dpi; mockDC and BAsubDC) DC inoculation (n ≧ 3). p-values by Student’s t test. Expression levels for each gene were normalized by that of an internal control UBQ1, which were further normalized with respect to those in Col-0 under the mock condition. Log2(normalized expression levels) are presented. C Expression of WRKY genes in Col-0 and ddm1 plants at early time points after DC inoculation (0, 6, and 12 h post-inoculation, hpi). UBQ1 expression was used as an internal control (n = 3). A, C Different letters indicate significant differences at p≦0.05, from two-way (A) and one-way (C) analysis of variance (ANOVA) with Tukey’s correction (α = 0.05). A–C Error bars indicate the standard error of the mean (SEM)

DDM1 mediates both CG and nonCG methylation, and ddm1 knockout thus partially compromises DNA methylation in all cytosine contexts [10, 15]. A previous study showed that ddm1 plants exhibit enhanced resistance to DC [23]. However, we did not observe altered immunity in ddm1 plants (Fig. 1A, white bar, ddm1), suggesting that the partial DNA methylation loss was insufficient to induce disease resistance under our experimental conditions. Thus, we assisted the defense response by applying BABA. High-dose BABA treatment (10 ml of 35 μg/ml; BAopt) for 3 days can prime an immune response in Col-0 plants [41,42,43] (Additional file 1: Fig. S1A and S1B). However, lowering the dose (10 ml of 30 μg/ml BABA) and shortening the treatment time (2 days), named “sub-optimal” (BAsub), did not enhance resistance in Col-0 plants (Fig. 1A and Additional file 1: Fig. S1C); therefore, whether the DNA methylation mutants suppress bacterial growth can be effectively evaluated in BAsub. BAsub treatment enhanced resistance only in met1 and ddm1 plants (Fig. 1A, red bars) and did not affect resistance in drm1 drm2, cmt3, ddc, and ddcc plants, suggesting that CG hypomethylation and BABA-induced chemical priming can synergistically induce resistance independently of CHG and CHH methylation.

We next examined the expression levels of three WRKY defense marker genes (WRKY29, WRKY70, and WRKY75) during defense response [45,46,47] (Fig. 1B and Additional file 1: Fig. S1D). In a water-treated mock control (mock), WRKY29 and WRKY70 expression was already upregulated in met1 plants compared to Col-0 (Fig. 1B). When we inoculated mock- and BAsub-treated plants with DC (mockDC and BAsubDC), WRKY expression was further hyperactivated at 1 day post-inoculation (1 dpi) in met1 plants (Fig. 1B and Additional file 1: Fig. S1D). By contrast, the expression of WRKY genes was not changed in ddm1 compared to Col-0 after mock, BAsub, and mockDC treatments, but hyperactivated after BAsubDC treatment (Fig. 1B and Additional file 1: Fig. S1D). These data are consistent with the resistance phenotype (Fig. 1A) because ddm1 plants were more resistant to DC only after being weakly primed with BABA, whereas met1 plants showed enhanced resistance under the mockDC condition, and BAsub treatment further fortified the immunity of plants. The WRKY expression level increased from 6 h post-DC inoculation (6 hpi) in weakly primed ddm1 plants, whereas Col-0 plants and non-primed ddm1 plants showed no or subtle upregulation of WRKY expression at 6 hpi (Fig. 1C). Therefore, the loss of ddm1 induced faster and stronger transcription of WRKY genes during a weakly primed defense response.

Effector-triggered immune response genes are hyperactivated in weakly primed ddm1 plants

To gain insight into how DDM1 modulates the defense transcriptome, we measured gene expression under four treatment conditions (mock, BAsub, mockDC, and BAsubDC) in Col-0 and ddm1 plants using microarray analyses (Fig. 2A and Additional file 1: Fig. S2A and S2B). Among 26,327 genes on the microarray, 11,138 were identified as differentially expressed genes (DEGs) between ddm1 and Col-0 under at least one treatment condition (Additional file 1: Fig. S2A). We categorized these genes into four groups based on their differential expression patterns under mock and BAsubDC conditions: (1) basally upregulated (ddm1 basal up; 1,718 genes in Additional file 1: Fig. S2A and S2B) or (2) downregulated genes (ddm1 basal down; 1148 genes) in ddm1 plants (vs. Col-0) under the mock condition; and (3) hyperactivated (ddm1 hyper; 3125 genes in Fig. 2A) or (4) suppressed (ddm1 sup; 3061 genes) genes in ddm1 plants (vs. Col-0) under the BAsubDC condition, but not under the mock condition (Figs. 2A and Additional file 1: Fig. S2A). Intriguingly, ddm1 hyper and sup genes showed no significant changes in ddm1 (vs. Col-0) under BAsub, but apparent up- or downregulation under mockDC, which was further boosted under BAsubDC, indicating potent priming effects (Fig. 2A). To investigate which group of genes contributed to enhanced resistance in ddm1, we analyzed the enriched Gene Ontology (GO) terms and found that ddm1 hyper genes were enriched with “defense response,” “response to external biotic stimulus,” and “defense response to fungus” while the other groups had no enrichment of defense-related GO terms (Additional file 1: Fig. S2D).

Fig. 2

Effector-triggered immune responsive genes were hyperactivated in weakly primed ddm1 plants. A Heatmaps of gene expression in Col-0 and ddm1 plants before (0 dpi; mock and BAsub) and after (1 dpi; mockDC and BAsubDC) DC inoculation. For Col-0 and ddm1 columns, log2 fold change (FC) against Col-0 mock samples were plotted. For ddm1/Col-0 columns, log2 fold change in ddm1 plants against Col-0 under four different conditions were plotted. B Heatmaps of gene expression changes after inoculation of avirulent (avrRpt2 and avrRpm1) or virulent (DC) pathogens. Log2 fold change (FC) against mock control at each time point was plotted. C Venn diagram of hyperactivated genes in ddm1 (vs. Col-0; ddm1 hyper), hyperactivated genes after avirulent pathogen inoculation (avrRpt2 or avrRpm1 vs. DC; ETI hyper), and upregulated genes in met1 (vs. Col-0; met1 up). D Observed to expected ratio of the number of overlaps between ETI hyper genes and other gene groups in C. ddm1 basal up indicates basally upregulated genes in ddm1 (vs. Col-0) under mock condition. E Venn diagram of suppressed genes in ddm1 (vs. Col-0; ddm1 sup), suppressed genes after avirulent pathogen inoculation (avrRpt2 or avrRpm1 vs. DC; ETI sup), and downregulated genes in met1 (vs. Col-0; met1 down). F Observed to expected ratio of the number of overlaps between ETI sup genes and other gene groups in E. ddm1 basal down indicates basally downregulated genes in ddm1 (vs. Col-0) under mock condition. D, F p-values by Fisher’s exact test

Avirulent pathogens possessing effectors trigger faster and stronger activation of defense-related genes than virulent pathogens, eventually fortifying plant immunity [48,49,50]. This boosted immune response is called effector-triggered immunity (ETI). To examine whether ddm1 hyper genes are also activated during ETI, we reanalyzed previously reported datasets [40] generated during ETI. Col-0 plants were inoculated with avirulent (Pseudomonas syringae pv. tomato avrRpt2 and avrRpm1) or virulent (DC) pathogens, and gene expression was measured from 1 to 24 hpi (Figs. 2B and Additional file 1: Fig. S2C). We then examined temporal gene expression patterns of ddm1 hyper genes in these datasets. ddm1 hyper genes were activated as early as 4 hpi with avrRpt2 and avrRpm1, whereas their expression peaked at 24 hpi with DC (Fig. 2B, top). ddm1 sup genes were quickly suppressed compared with DC (Fig. 2B, bottom). Both basal up and down genes in ddm1 were suppressed by inoculation with avirulent pathogens (Additional file 1: Fig. S2C). We next identified hyperactivated (ETI hyper) or suppressed (ETI sup) genes after treatment with avrRpt2 and avrRpm1 (vs. DC, differentially expressed at 6 or 9 hpi; p ≦ 0.01) and compared them with the above four gene groups in ddm1. We found that 46% of ETI hyper genes overlapped with ddm1 hyper genes, which was 3.86-fold more than the expected portion of random overlap (p < 10–5, Fig. 2C, D, Additional file 2: Supplementary data 1), and 34% of ETI sup genes overlapped with ddm1 sup genes (2.94-fold and p < 10–5, Fig. 2E, F, Additional file 2: Supplementary data 1). In contrast, only 9% of ddm1 basal up genes overlapped with ETI hyper genes, and 11% of ddm1 basal down genes overlapped with ETI sup genes (Fig. 2D, F, Additional file 2: Supplementary data 1). The early activation of ddm1 hyper genes and their significant overlap with ETI hyper genes suggest potential associations of ddm1 hyper genes with induced immunity in ETI.

Moreover, met1, drm1 drm2, ddc, and ddcc plants showed enhanced resistance to DC without BAsub (Fig. 1A). Comparison of the resistance between mutants and their corresponding wild-type plants under mockDC condition revealed that both CG demethylation (met1) or CHH demethylation (met1, drm1 drm2, ddc, and ddcc) enhanced resistance to DC, whereas CHG demethylation (cmt3) has little effect on resistance to DC (Fig. 1A). Consistent with resistance in met1 plants, up- and downregulated genes in met1 plants under the mock condition [51] significantly overlapped with ETI hyper and ETI sup genes (Fig. 2C–F). However, comparison of the resistance between mockDC and BAsubDC showed that BAsub treatment boosted resistance to DC only in mutants showing CG demethylation (met1 and ddm1), consistent with significant overlaps of ETI-responsive genes with ddm1 hyper or sup genes in weakly primed conditions. These results suggest that global loss of CG methylation in met1 and ddm1 plants contributed to the enhanced disease resistance with weak BAsub priming, whereas loss of CHH methylation enhanced resistance to DC independent of priming.

DDM1 mediates both TE-like and gene body-like methylation

To investigate how BAsub treatment boosted immunity in CG-hypomethylated plants, we focused on ddm1 because met1 plants also showed enhanced resistance without BAsub. We first produced genome-wide DNA methylation profiles in Col-0 and ddm1 plants with no treatment of BAsub or DC. We then identified CG or nonCG methylated regions from the profiles in Col-0 by selecting 50-bp windows with significant DNA methylation in CG or nonCG (CHG or CHH) contexts and then merging nearby windows (within 500 bp). The distributions of CG, CHG, and CHH methylation levels in the merged windows are shown in Additional file 1: Fig. S3A. Based on these distributions, the merged windows having more than 10% CG, 5% CHG, and 2% CHH methylation were defined as TE-like methylation (TEM-like) regions while those having more than 10% CG but less than 5% CHG and 2% CHH methylation were as gene body-like methylation (GBM-like) regions. In Col-0 plants, 67% of methylated regions were TEM-like, and 33% were GBM-like ones (Additional file 1: Fig. S3B).

As expected, DNA methylation was severely disrupted at TEM-like sites in ddm1 plants, indicating the crucial role of DDM1 in maintaining DNA methylation at TEs (Fig. 3A) [10, 15]. For integrated analysis of TEM-like in ddm1, met1, and drm2, we reanalyzed previous bisulfite datasets [10, 52] generated from met1 and drm2 plants to identify TEM-like in these mutants. When we isolated DDM1- and MET1-mediated TEM-like that simultaneously lost CG, CHG, and CHH methylation in ddm1 and met1 plants (Additional file 1: Fig. S3C), demethylated TEM-like in ddm1 (TEM-likeddm1) had high levels of H3K9me2 and histone H1 in Col-0. In contrast, demethylated TEM-like in met1 (TEM-likemet1) had low levels of H3K9me2 and H1 (Additional file 1: Fig. S3D) [10, 53], significantly overlapping with CHH demethylated TEM-like in drm2 (TEM-likedrm2) (Additional file 1: Fig. S3C and S3E).

Fig. 3

Transposon-like methylation (TEM-like) and gene body-like methylation (GBM-like) patterns in ddm1 plants. A Violin plots of DNA methylation levels of TEM-like regions in Col-0 and ddm1 plants. B Violin plots of CG methylation levels of GBM-like regions in Col-0 and ddm1 plants. C Scatter plot of CG methylation level of GBM-like regions in Col-0 (x-axis) and ddm1 (y-axis) plants. Each dot represents a gene body methylation locus. Orange dots indicates GBM-likeddm1, and blue dots indicate GBM-likeweak. GBM-likeddm1; strongly demethylated GBM-like regions in ddm1. GBM-likeweak; GBM-like regions that their methylation levels are weakly affected or not affected in ddm1. D Histone H1.1 and H1.2 level at GBM-likeddm1 and GBM-likeweak. p-values from Student’s t test. E, F Gene body CG DNA methylation level at GBM-likeddm1 (E) and GBM-like.weak (F) in Col-0, h1, ddm1, and h1ddm1 plants. G Scatter plots of CG methylation change at TEM-like regions in relation to histone H1 and H2A.Z change in ddm1. Magenta box indicates highly demethylated TEM-like regions with increased DNA accessibility in ddm1. H Scatter plots of CG methylation change at GBM-like regions in relation to histone H1 and H2A.Z change in ddm1. G, H Dot color shows log2 ratio of DNA accessibility in Col-0 to ddm1

On the other hand, consistent with a previous finding [15], DDM1 had a limited effect on GBM-like regions (Fig. 3B). The distribution of CG methylation levels in Col-0 and ddm1 showed two groups of GBM-like regions (Figs. 3C and Additional file 1: Fig. S3F, Additional file 2: Supplementary data 1): (1) GBM-likeddm1 group showing significantly lost CG methylation in ddm1 (Fig. 3C, orange dots in bottom box) and (2) GBM-likeweak group showing relatively weak or no CG methylation changes (Fig. 3C, blue dots in top box). As ddm1 plants can exhibit random fluctuations of DNA methylation, we tested whether GBM-like demethylation in ddm1 occurred randomly. CG methylation changes in these two groups were consistent to those observed in other independent data sets (Additional file 1: Fig. S3G for GBM-likeddm1 and Additional file 1: Fig. S3H for GBM-likeweak). Previous studies suggested that DDM1 enables the access of DNA methylases to the H1-rich DNA. Hence, ddm1 plants lose DNA methylation in H1-rich regions, and the lost methylation is rescued by knockdown of h1 [10, 15]. Consistent with this, GBM-likeddm1 regions had higher H1 levels than GBM-likeweak regions (Fig. 3D), and regained DNA methylation in h1ddm1 (Fig. 3E) [51]. In contrast, DNA methylation at GBM-likeweak regions showed subtle changes in ddm1 or h1ddm1 plants (Fig. 3F) [10, 15]. These results imply that DDM1 mediates GBM-like at specific loci, e.g., H1-rich regions.

DDM1, a chromatin remodeler, affects DNA accessibility, DNA methylation, H3K9me2 levels, and H2A.W deposition at heterochromatin [10, 51, 54], and perturbation of DNA methylation has been suggested to indirectly modulate H1 and H2A.Z distribution in Arabidopsis [33, 55]. Therefore, we next investigated whether changes of DNA methylation in ddm1 correlated with those of H1, H2A.Z, H2A.W, and H3K9me2. To this end, we performed chromatin immunoprecipitation (ChIP) sequencing of H1 and H2A.Z in Col-0 and ddm1 plants and used previously reported ChIP sequencing data for H2A.W and H3K9me2 [54]. For the TEM-like, changes of CG methylation levels between ddm1 and Col-0 showed positive correlations for highly demethylated TEM-like in ddm1 (more than 0.5 loss of CG methylation in ddm1) with those of H1, H2A.W, and H3K9me2 (Figs. 3G and Additional file 1: Fig. S3I, magenta box). On the other hand, there was no significant correlation for H2A.Z (Fig. 3G). However, when the same analysis was done using a previously reported H2A.Z profile generated from Col-0 and met1 [33, 52], CG methylation changes in met1 (vs. Col-0) showed a stronger correlation for highly demethylated TEM-like (Additional file 1: Fig. S3K, magenta box) compared to in ddm1, potentially due to suppression of H2A.Z accumulation by the remaining CG methylation after the limited demethylation in ddm1, unlike extensive demethylation in met1 (compare CG methylation levels of ddm1 (Fig. 3A) and met1 (Additional file 1: Fig. S3L)). These data suggest that demethylation at TEM-like in ddm1 is independent to H2A.Z changes and partially correlated with changes of H1, H2A.W, and H3K9me2 for highly demethylated TEM-like. Moreover, we also examined the association of CG methylation changes with those of H1, H2A.Z, H2A.W, and H3K9me2 in all GBM-like regions. Unlike TEM-like, CG methylation changes in GBM-like regions showed no correlation with those of and H1, H2A.Z, H2A.W, and H3K9me2 (Figs. 3H and Additional file 1: Fig. S3J). Differently from the negative correlation between CG methylation and H2A.Z changes for highly demethylated TEM-like, in case of GBM-like, there was no significant correlation in met1 regardless of demethylation extents (Additional file 1: Fig. S3K). Therefore, at GBM-like regions, ddm1 knockout specifically regulated DNA methylation, but independently with histone composition changes, unlike the diverse correlation patterns of ddm1 knockout with histone composition changes at TEM-like.

DDM1-mediated TE-like methylation at gene bodies and transcription start sites potentiates transcription

TEs can exist in the promoter and gene body to affect expression of nearby genes. Thus, we next investigated whether TEM-likeddm1 is associated with gene expression in a weakly primed defense response. To precisely evaluate the unique effect of DDM1-mediated TEM-like (TEM-likeddm1) at promoter, transcription start site (TSS), or gene body to gene expression, we selected the following three mutually exclusive groups of genes (Additional file 1: Fig. S4A): (1) 4702 genes with TEM-like only at the promoter (pTEM-like), (2) 1250 genes with TEM-like only at TSS (TSS TEM-like), (3) 1187 genes with TEM-like only at gene body (gbTEM-like). We then examined the DNA methylation patterns of these three gene groups in Col-0 and ddm1 plants over the gene structure under four treatment conditions (mock, BAsub, mockDC, and BAsubDC) (Fig. 4A). The CG and nonCG methylation patterns around genes were similar across treatment conditions (Fig. 4A and Additional file 1: Fig. S4B–S4D). Correspondingly, the methylation difference between ddm1 and Col-0 over the gene structure were similar across different conditions with high correlation coefficients (r = 0.73 ~ 0.94; Additional file 1: Fig. S4E–S4G).

Fig. 4

DDM1-mediated transposon-like methylation (TEM-like) at transcription start sites (TSS) and gene bodies is associated with augmented gene expression during defense response under sub-optimal priming conditions (BAsubDC). A Averaged CG DNA methylation level around genes with TEM-like at TSS, gene bodies, or promoters. B–D Boxplots of gene expression of DEGs with strongly demethylated TEM-like in ddm1 (TEM-likeddm1) at TSS, gene body, and promoter (TSS TEM-likeddm1, gbTEM-likeddm1, pTEM-likeddm1). Expression fold change (log2 (ddm1/Col-0)) of all DEGs at each condition was plotted as control. DEG, differentially expressed gene. E CG, CHG, CHH methylation, histone H1, H2A.Z, and gene expression profiles of genes with TEM-likeddm1 at TSS, gene body, and promoter (TSS TEM-likeddm1, gbTEM-likeddm1, pTEM-likeddm1) and TEM-likeweak at TSS, gene body, and promoter (TSS TEM-likeweak, gbTEM-likeweak, pTEM-likeweak) in Col-0 and ddm1 under the four treatment conditions. TEM-likeddm1; strongly demethylated TEM-like in ddm1, TEM-likeweak; TEM-like that their methylation levels are weakly affected or not affected in ddm1. B–E p-values by Student’s t test

Next, we analyzed the effects of the three gene groups on gene expression. To this end, we first selected all DEGs between ddm1 and Col-0 plants in each treatment condition (e.g., BAsubDC). From the three gene groups, we then selected the DEGs with demethylated pTEM-like (pTEM-likeddm1), TSS TEM-like (TSS TEM-likeddm1), and gbTEM-like (gbTEM-likeddm1) in ddm1 (vs. Col-0) under the treatment condition (e.g., BAsubDC). We next compared the expression changes of all DEGs and DEGs with each of the three TEM-likeddm1 (Fig. 4B–D and Additional file 2: Supplementary data 1) in ddm1 (vs. Col-0) under the treatment condition (e.g., BAsubDC). When one of the three TEM-likeddm1 had strong effects on gene expression, DEGs with the TEM-likeddm1 would show higher expression changes than all DEGs. This procedure was performed for all four treatment conditions (mock, BAsub, mockDC, and BAsubDC). DEGs with TSS TEM-likeddm1 were hyperactivated under BAsub, mockDC, and BAsubDC conditions, although they showed no changes under the mock condition, compared to all DEGs (Fig. 4B). Also, DEGs with gbTEM-likeddm1 were hyperactivated predominantly in the BAsubDC condition (Fig. 4C). In contrast, pTEM-likeddm1 did not affect gene expression, except under the mockDC condition (Fig. 4D).

We then integrated DNA methylation, H1, H2A.Z, and gene expression levels at TEM-like under the four conditions to examine correlations of DNA methylation, H1, and H2A.Z changes with expression changes of genes with the above three groups of TEM-likeddm1 or TEM-likeweak. Among the three methylation contexts, genes with TSS TEM-likeddm1 showed hyperactivation, which was much stronger than genes with TSS TEM1-likeweak (Fig. 4E, expression and the boxplot on top), and CG demethylation was stronger in these genes than CHG demethylation while CHH methylation showed no significant changes (Fig. 4E, CG, CHG, and CHH), consistent with the finding in Fig. 3A. Moreover, H2A.Z levels showed no correlation with hyperactivation of genes with TSS TEM-likeddm1 while H1 showed a partial correlation (Fig. 4E and Additional file 1: Fig. S4H, H1 and H2A.Z), consistent with the finding in Fig. 3G. Taken together, TSS TEM-likeddm1 and gbTEM-likeddm1, to a lesser extent, augmented transcription especially under the BAsubDC condition.

DDM1-mediated GBM potentiates transcription

We next investigated the effect of GBM-likeddm1 on gene expression. GBM-like can also occur in the promoter and TSS to affect expression of target genes. Similar to genes with TEM-likeddm1, we thus selected three mutually exclusive groups of genes (Additional file 1: Fig. S5A): (1) 2157 genes only with pGBM-like, (2) 83 genes only with TSS GBM-like, and (3) 6816 genes only with gbGBM-like (i.e., GBM). The TSS GBM-like gene group, including only 83 genes, was excluded for the downstream analyses because they may provide insufficient statistical power. Like TEM-like, CG methylation levels of the GBM-like were reduced in ddm1 plants and showed similar patterns over gene structure across the treatment conditions (Fig. 5A).

Fig. 5