Methylome Response to Proteasome Inhibition by Pseudomonas syringae Virulence Factor Syringolin A

Plants are sessile organisms and must constantly adapt their gene expression to react to constantly changing environments. Among other environmental stresses, the infection of plants by pathogens, such as bacteria or fungi, is a constant threat to plant fitness and, thus, leads to significant losses in crop production. Plants have evolved several ways to defend themselves against pathogens, starting with physical barriers, such as the cuticula, to specific intra-cellular pathways that recognize pathogens. These intra-cellular pathways can be divided into PTI (pattern-triggered immunity) or ETI (effector-triggered immunity), depending on the cellular localization of the plant receptor. For PTI, the receptor is located on the cell-surface, whereas, for ETI, the receptor is found intra-cellularly. The ETI involves the specific recognition of effector molecules that are secreted by the pathogen and counteract the plant primary immune response (Jones and Dangl 2006).

Interestingly, several virulence factors produced by bacteria, viruses, and fungi were found to inhibit the host proteasome machinery (Chiu et al. 2010; Dudler 2013; Groll et al. 2008; Jin et al. 2007; Marino et al. 2012; Sahana et al. 2012; Sorel et al. 2019; Üstün and Börnke 2014; Üstün et al. 2016; Verchot 2016) and, more generally, to regulate protein degradation (Langin et al. 2020). Thus, proteasome inhibition seems to represent a common strategy used by plant pathogens. Additionally, host proteasome inhibition has been shown to suppress systemic acquired resistance (SAR), which is responsible to increase pathogen resistance in the whole plant following a localized infection. The inhibition of SAR by proteasome inhibition leads to increased plant susceptibility to bacterial pathogens (Üstün et al. 2016). Beyond classical bacterial type III effectors (a class of virulence factors directly injected into host cells) (Shames and Finlay 2012) such as XopJ and HopM1 (Üstün et al. 2013; Üstün et al. 2016), which inhibit proteasome activity by interacting with proteasome subunits, a virulence factor called syringolin A (SylA) inhibits the proteasome via an irreversible covalent binding to the proteasome catalytic subunits (Groll et al. 2008). SylA is a small tri-peptide derivative, which is secreted by the bacteria Pseudomonas syringae pv. syringae. Bacteria deficient for SylA secretion are less virulent than their wild-type counterpart. SylA secretion is important for the bacteria to overcome stomata closure and consequently gaining better leaf penetration (Schellenberg et al. 2010). Similarly, SylA also plays a role for wound entry and colonization of the host via the vascular tissue (Misas-Villamil et al. 2013). Interestingly, SylA was found to accumulate in the nucleus, which led the authors to suggest that it could act preferentially on the nuclear proteasome (Kolodziejek et al. 2011).

The proteasome is one of the main protein degradation machineries in eukaryotic cells. Protein targeting to the proteasome relies in part on polyubiquitination by E3 ubiquitin ligases (Sadanandom et al. 2012). Due to its central role in plant physiology, the proteasome machinery is involved in several biological processes, such as hormonal signaling, circadian clock, as well as plant stress response (Sadanandom et al. 2012; Vierstra 2009). Additionally, the proteasome is also involved in the regulation of chromatin and transcription (Geng et al. 2012). The role of the proteasome in the regulation of chromatin has not only been linked to the regulation of chromatin components, such as histones and their assembly into nucleosomes, but also to chromatin regulators such as histone chaperones and histone modifiers. (Bach and Hegde 2016; Jeong et al. 2011; Kinyamu et al. 2008; Lee et al. 2011). In addition, proteins involved in DNA methylation were found to be regulated by the proteasome.

DNA methylation in Arabidopsis is characterized by the apposition of a methyl group to cytosine residues. In plants, DNA methylation occurs in three distinct sequence contexts: CG, CHG, and CHH, where H stands for any nucleotide except cytosine, and each methylation context is regulated by its own cognate pathway (Law and Jacobsen 2010). CG methylation is maintained by the maintenance DNA methyltransferase MET1 (METHYLTRANSFERASE 1). Despite the lack of evidence for a proteasome regulation of the MET1 protein, other important proteins regulating CG methylation, such as the VIM proteins, are E3 ubiquitin ligases (Johnson et al. 2007; Kraft et al. 2008; Woo et al. 2007), suggesting a link between ubiquitination and CG methylation. The main enzyme linked to CHG methylation is the DNA methyltransferase CMT3 (CHROMOMETHYLASE 3) (Lindroth et al. 2001). CMT3 itself was found to be targeted for proteasomal degradation by an atypical JmjC domain protein, JMJ24, which has E3 ubiquitin ligase activity (Deng et al. 2016). Recently, DRM2 (DOMAINS REARRANGED METHYLASE 2), the main DNA methyltransferase involved in de novo cytosine methylation at CHH sites was also found to be regulated by an E3 ubiquitin ligase named CFK1, leading to its targeting for degradation by the proteasome (Chen et al. 2021). DNA methylation in all sequence contexts can, therefore, be directly or indirectly linked to ubiquitin-dependent proteasome degradation.

DNA methylation has previously been shown to be involved in the response against bacterial and fungal pathogens (Deleris et al. 2016; Zhu et al. 2016). Indeed, mutants with lower DNA methylation levels, such as met1, were shown to display an increased resistance to Pseudomonas bacteria (Dowen et al. 2012), while mutants exhibiting global DNA hypermethylation, such as ros1, were shown to be more susceptible to Pseudomonas, Fusarium, and Hyaloperonospora infection (Le et al. 2014; López Sánchez et al. 2016; Schumann et al. 2019; Yu et al. 2013). DNA de novo methylation is known to be targeted by small RNA (sRNA) molecules of principally 24 nucleotides (nt) in size. This pathway is referred to as RNA-directed DNA methylation (RdDM) (Matzke and Mosher 2014). Several members of the RdDM pathway were also found to be involved in Botrytis, Plectosphaerella, and Pseudomonas resistance (Agorio and Vera 2007; López et al. 2011), and Pseudomonas bacterial infection was shown to affect the genome-wide DNA methylation pattern (Dowen et al. 2012). Additionally, DNA methylation is involved in transgenerational memory of bacterial infection, a phenomenon referred to as transgenerational acquired resistance (TAR) (Luna and Ton 2012; Luna et al. 2012).

Noting the converging indications for a role of DNA methylation during bacterial infection as well as proteasome inhibition during this process, we investigated the effect of the exposure by the virulence factor and proteasome inhibitor SylA on global DNA methylation levels. Here, we show that treatment with SylA affects genome-wide DNA methylation and results in a moderate genome hypermethylation in all sequence contexts. This hypermethylation is mainly affecting the centromeric region of Arabidopsis chromosomes. We identified 10,101 100-bp differentially methylated regions (DMRs), of which 6,304 are CHH-hypermethylated DMRs. We could not identify major changes in sRNA composition upon SylA treatment, suggesting that sRNAs might not have a direct effect on the observed DNA methylation changes. Furthermore, the Arabidopsis transcriptome is strongly affected by SylA-mediated proteasome inhibition. We propose that the induction of atypical members of the RdDM pathway, such as DRM1, AGO9, and AGO3, might explain the changes observed upon SylA treatment.

Results

SylA induces moderate centromeric hypermethylation

To investigate the effect of proteasome inhibition on the Arabidopsis methylome, we performed genome-wide bisulfite sequencing (GWBS) of Arabidopsis plantlets treated with the bacterial virulence factor SylA, a known proteasome inhibitor (Groll et al. 2008). Considering all three sequence contexts, we obtained DNA methylation levels for 14,269,258 and 10,867,987 cytosines for mock- and SylA-treated GWBS libraries, respectively, that reached more than 10-fold coverage. The obtained data were plotted as a circular heat map using 50-kb genomic bins (Fig. 1A). We observed that the distribution of the genome-wide DNA methylation did not change between treated and the control samples in all sequence context CG, CHG, and CHH (Fig. 1A). As previously observed (Cokus et al. 2008; Lister et al. 2008; Stroud et al. 2012), DNA methylation is higher in centromeric and pericentromeric regions of the Arabidopsis genome, coinciding with increased transposable element (TE) densities. As expected, the percentage of DNA methylation is higher in CG context (28.36 and 29.25%) than in CHG context (12.05 and 13.62%) and, finally, is at its lowest in the CHH context (3.16 and 3.54%) in both control and treated samples. To refine our analyses, we focused on methylated windows (i.e., >40% methylation for CG, >20% for CHG, and >10% for CHH) (Fig. 1B). We could observe a slight but consistent increase in DNA methylation levels upon SylA treatment in all sequence contexts. To investigate in which chromosomal environment this increase of DNA methylation occurs, we plotted the percentage of changes in DNA methylation across the genome. We observed that, except for a few discreet loci on chromosome arms, most of the increase in DNA methylation occurred in the centromeric and pericentromeric regions of all Arabidopsis chromosomes (Fig. 1C; Supplementary Fig. S1). We conclude that proteasome inhibition by SylA treatment does not result in a global reprograming of genome-wide DNA methylation but, rather, in a moderate hypermethylation of the (peri-)centromeric region.

Download as PowerPoint

SylA treatment results mainly in CHH DMRs

To investigate the changes in DNA methylation on the local scale, we identified DMRs of 100 bp in size for each DNA methylation context. To satisfy the criteria as a DMR, the respective 100-bp genomic bins had to differ by at least 10% between treatments and exhibit a q value < 0.05. We identified a total of 10,101 individual genomic bins satisfying DMR criteria (Fig. 2A). The vast majority of DMRs were found in the CHH context (9,173), whereas only a minority was detected in CG (409) and CHG (672) contexts. Consistent with our global observations, we observed a bias towards hyper DMRs (gain of DNA methylation) rather than hypo DMRs (loss of DNA methylation). CHH DMRs mostly corresponded to TEs (51%), untranslated regions (UTRs) (14%), promoters (14%), and intergenic sequences (15%) and were rather depleted in genic sequences (2%) (Fig. 2B). Similarly, CG and CHG DMRs were also preferentially associated to TEs (Supplementary Fig. S2). The overlap between the DMRs of different contexts is limited yet significantly greater than expected by chance (permutation-based P value using 10⁴ repetitions for all dual overlaps and triple overlaps <10⁻⁴). Considering the frequent association of DMRs of all contexts with TEs, we observed an especially pronounced overlap when regarding DMR-associated TEs (Supplementary Fig. S3), suggesting that most significant changes in DNA methylation occurs within or in the proximity of TEs. TEs associated with DMRs were not enriched for a particular TE length or family (Supplementary Fig. S4A and B). To visualize in which genomic region DMRs are located, we plotted the DMRs along the five Arabidopsis chromosomes (Fig. 2C). Both hypo and hyper CHH DMRs are enriched in centromeric and pericentromeric regions, mirroring the TE density distribution in the Arabidopsis genome. To better characterize the location of DMRs and their potential effect on gene expression, we computed the distribution of their mean distance to transcriptional start sites (TSSs) and compared it with that of methylated regions not qualifying as DMRs. For this purpose, we separated our analysis between the DMRs occurring in the pericentromeric regions and the DMRs occurring on chromosomal arms. Interestingly, SylA-induced DMRs in chromosomal arms are closer to a TSS when compared with randomly sampled methylated regions for both DMRs associated with TEs (Fig. 2D) and genes (Fig. 2E). This increased proximity to a TSS is more pronounced for CHH DMRs (Fig. 2D and E) than for CG and CHG DMRs (Supplementary Fig. S5A and B). No significant differences were observed for DMRs situated in pericentromeric regions (Supplementary Fig. S5C and D). In summary, our results show that SylA treatment results principally in CHH DMRs situated in the centromeric and pericentromeric regions and that DMRs situated in chromosomal arms are near TSSs.

Download as PowerPoint

As previously mentioned, DNA methylation changes were observed upon Arabidopsis infection by the bacteria Pseudomonas syringae pv. tomato (Dowen et al. 2012). To evaluate how similar SylA DMRs are from P. syringae pv. tomato–induced DMRs, we analyzed the overlap between P. syringae pv. tomato and SylA DMR-associated features (Supplementary Fig. S6A, B, and C). Although no significant overlap could be observed for CG DMRs or CHG DMRs (Supplementary Fig. S6A and B), a significant overlap could be observed concerning CHH DMRs, with 413 overlapping loci (P value < 10⁻⁴) (Supplementary Fig. S6C). The overlap between P. syringae pv. tomato and SylA DMRs suggests that these DMRs could be linked to proteasome inhibition taking place during P. syringae pv. tomato infection. Further studies will be required to confirm this hypothesis.

SylA treatment does not influence sRNA population

DNA methylation in CHH context is closely linked to the RNA-directed DNA methylation pathway (RdDM). This pathway depends on sRNAs to induce DNA methylation at complementary DNA loci. In its canonical form, RdDM relies on the loading of 24-nt sRNAs by AGO4 or AGO6 and the recruitment of the de novo methyltransferase DRM2 to the DNA (Law and Jacobsen 2010; Matzke and Mosher 2014). To investigate potential changes in sRNA populations between mock-treated and SylA-treated plantlets, we profiled the sRNA populations using next-generation sequencing. As expected in both conditions, 24-nt sRNAs are the most abundant size fraction, mostly starting with a 5′A, and map predominantly to TEs (Fig. 3A, B, and C). In both samples, the second most abundant are 21-nt sRNAs. They mainly start with a 5′U and likely correspond to micro RNAs (miRNAs). The 22-nt sRNAs are the least abundant of the three sizes, predominantly starting with a 5′G and mapping to TEs. We conclude that the overall sRNA composition is very similar between the control and the treated samples and comparable to previously published Arabidopsis sRNA compositions (Mi et al. 2008). To investigate the sRNA population in more detail, we selected 100-bp genomic bins showing differential expression (exhibiting at least 10 reads across the two treatments and showing at least a fourfold change between the treatments). These differentially expressed bins are subsequently referred as SylA-enriched (Supplementary Fig. S7). We found that, despite being very similar, several bins (n = 1,539) mostly comprising 24-nt sRNA (n = 1,016) were enriched in the SylA-treated sample (Fig. 3D). Bins enriched in 21- and 22-nt sRNA were not as abundant (21-nt = 326, 22-nt = 197) (Fig. 3D). For all sRNA sizes, we identified only a low number of sRNA-depleted bins (n = 192). As expected, the SylA-enriched 24-nt sRNA were predominantly starting with a 5′A and mapped to TEs (Supplementary Fig. S7A, B, and C). The SylA-enriched 21- and 22-nt sRNA had an unusual 5′C and 5′U enrichment, respectively, and mapped predominantly to TEs (Supplementary Fig. S7A, B, and C). In comparison to the total library populations, 21-nt sRNA, mostly representing miRNAs, were underrepresented in the SylA-enriched sRNA population, indicating that miRNAs are not activated upon SylA exposition. Bins enriched in 24-nt sRNA are located principally in pericentromeric regions (Supplementary Fig. S8). However, only 49 of the 1,055 24-nt SylA-enriched bins overlapped with CHH DMR bins (Fig. 3E, P = 0.9437). The same observation was made for 21- and 22-nt–enriched sRNA bins (Supplementary Fig. S9A, B, and C). Additionally, we could not find a significant correlation between the increase in 24-nt sRNAs and the CHH methylation levels (Fig. 3F). The absence of significant correlation was also observed for 21- and 22-nt–enriched bins (Supplementary Fig. S9D, E, and F). We conclude that the sRNA composition is mostly unaffected after 24-h SylA treatment and that the moderate changes in sRNA composition or abundance do not correlate with the changes in CHH DNA methylation and is likely not the cause of these changes.

Download as PowerPoint

SylA treatment predominantly affects transcription in chromosomal arms

DNA methylation is a well-known regulator of gene and TE expression. To investigate a potential link between the observed changes in DNA methylation upon SylA treatment and genome-wide transcriptional activity, we subjected our samples to RNA sequencing. Upon SylA treatment, we could detect a major change in transcriptional activity (Fig. 4A and B), identifying a total of 3,241 differentially expressed transcripts (adjusted P value < 0.01 and absolute log fold change value > 2), among them 2,945 differentially expressed genes (DEGs). We noted a significant bias toward upregulated DEGs (1,675) compared with downregulated DEGs (1,270), showing that SylA treatment leads to transcriptional upregulation rather than downregulation (χ² = 55.696, P = 8.4586E-14). Interestingly, gene ontology (GO) term enrichment analysis (Fig. 4C and D; Supplementary Table S1) for the upregulated DEGs revealed an enrichment in “positive regulation of RNA polymerase II initiation” (GO:0045899, P = 9.4 × 10⁻⁶) (Fig. 4D; Supplementary Table S1). Our results suggest that an inhibition of the proteasome using SylA results in increased transcriptional activity. We can observe a similar GO term enrichment in “positive regulation of RNA polymerase II initiation” (GO:0045899, P = 2.7 × 10⁻⁹), using previously published transcriptomic data using MG132, a non-covalent proteasome inhibitor (Supplementary Table S1) (Gladman et al. 2016). This suggests that this feature is not only linked to SylA itself but more generally to proteasome inhibition. Increased transcriptional activity in response to proteasome inhibition seems to be widely conserved, as it was also observed in human cells treated with MG132 (Kinyamu et al. 2008).

Download as PowerPoint

As previously reported upon SylA treatment (Michel et al. 2006), we observed an enrichment in “proteasome assembly” (GO:0043248, P value = 10⁻⁹). We hypothesize that, to compensate for the inhibition of the proteasome, proteasome-related genes are transcribed at higher levels. This phenomenon was also observed using MG132 (Supplementary Table S1) (Gladman et al. 2016). Similarly, an increase of proteasome-related gene expression also occurs during P. syringae pv. tomato infection (Üstün et al. 2018), further supporting that proteasome inhibition genuinely occurs during P. syringae pv. tomato infection. In addition, upregulated DEGs are enriched in GO terms describing several abiotic and biotic stresses (Fig. 4D). Downregulated DEGs are clearly enriched in GO terms related to photosynthesis as well as response to light (Fig. 4C; Supplementary Table S1). A link between ubiquitin-dependent proteasome degradation and light response in Arabidopsis has been long-established with the discovery of the COP9 signalosome in plants (Chamovitz et al. 1996; Cope and Deshaies 2003; Wei and Deng 2003). In addition to modifying gene expression, SylA treatment also affects the expression of TEs, 217 of them exhibiting significantly differential expression (Fig. 4A and B), with a strong tendency towards being upregulated rather than downregulated (145 up, 72 down). Surprisingly and similarly to DEGs, we predominantly detected differentially expressed TEs located on chromosomal arms rather than pericentromeric or centromeric region (Fig. 4B) despite the majority of TEs being located there. Upregulated TEs are enriched in long terminal repeat copia TE families (Supplementary Fig. S4C and D), whereas downregulated TEs tend to be enriched in DNA-MuDR transposons (Supplementary Fig. S7C and E). To investigate whether SylA induced DMRs could be responsible for DEG expression levels, we performed a correlation analysis between SylA DMRs methylation level (DMR situated within the gene annotation extended by ±1 kb) and the log fold change of SylA-induced DEGs (Fig. 4E, F, and G). We could not detect any correlation in all three DNA methylation contexts. Additionally, we did not observe any significant change in the mean DNA methylation level across up- and downregulated genes or TEs (Supplementary Fig. S10A and B). Despite these analyses, loci fitting the assumption “loss of DNA methylation leads to increased expression” can still be found. A snapshot of such a locus is illustrated in Supplementary Figure S10C.

Our result shows that SylA treatment has a significant effect on the Arabidopsis transcriptome, which is likely linked to the role of SylA as a proteasome inhibitor rather than its direct effect on DNA methylation.

SylA treatment induces the expression of atypical DNA methylation factors

To understand how SylA could influence DNA methylation beyond its role in inhibiting the proteasome, we investigated if it could act by modulating the expression of key actors of DNA methylation. We have focused our analyses on the key members of the different DNA methylation pathways and homologous genes encoded in the Arabidopsis genome (Fig. 5A), additional genes of interest can be found in Supplementary Figure S11A to D. Methylation on CG sites requires the principal CG DNA methyltransferase MET1 (Finnegan et al. 1996; Kankel et al. 2003; Saze et al. 2003). Additionally, on centromeric and pericentromeric sequences, the chromatin remodeler DDM1 is required for the maintenance of CG methylation (Jeddeloh et al. 1998). In addition, the Arabidopsis genome encodes three other MET genes, namely, MET2a, MET2b, and MET3 (Jullien et al. 2012). None of these genes affecting CG methylation showed a sufficient expression fold change to meet DEG criteria. Similarly, the CHG DNA methyltransferases (CMT1, CMT2, and CMT3) (Henikoff and Comai 1998; Lindroth et al. 2001) and the DNA demethylase (ROS1, DME, DML2, and DML3) (Choi et al. 2002; Gong et al. 2002; Penterman et al. 2007) were not among the DEGs. However, we observed a significant upregulation for some members of the de novo DNA methylation RdDM pathway, which elicits CHH methylation. We analyzed the de novo DNA methyltransferases DRM1, DRM2, and DRM3 and observed a significant induction of DRM1 [log₂(FC) = 5.9]. In the RdDM pathway, sRNAs loaded in proteins called Argonautes are targeting de novo methylation on the complementary DNA sequence. The Argonautes associated with the RdDM pathway are AGO4, AGO6, AGO9, AGO8, and AGO3 (Duan et al. 2015; Havecker et al. 2010; Huang et al. 2016; Stroud et al. 2012; Zilberman et al. 2003). We found that AGO3 [log₂(FC) = 3.4] and AGO9 [log₂(FC) = 2.2] were significantly upregulated in response to SylA treatment. Performing quantitative PCR (qPCR) on Arabidospis complementary DNA (cDNA), we could confirm that SylA treatment leads to an upregulation of AGO3, AGO9, and DRM1 (Fig. 5B). To investigate whether AGO3, AGO9, and DRM1 could be regulated by a common transcription factor, we analyzed their promoter sequences for putative transcription factor binding sites, using the Gene Analysis Tool from AthaMap (Supplementary Table S2) (Steffens et al. 2005). We could find one potential Arabidopsis transcription factor targeting all three genes, the AT-Hook AHL12(1) (AT1G63480); however, AHL12 was not found among the SylA-DEGs. Further studies would be required to investigate whether AHL12 might be regulated at the protein level by the proteasome and is thereby responsible for AGO3, AGO9, and DRM1 upregulation upon SylA treatment.

Download as PowerPoint

Discussion

Our study shows that SylA, a virulence factor secreted by P. syringae pv. syringae strains, induces a moderate change of the Arabidopsis methylome. Differentially methylated regions are mostly located in the centromeric and pericentromeric regions of the Arabidopsis chromosomes. Upon SylA treatment, substantial changes in the Arabidopsis transcriptome are observed, which are most likely linked to SylA action on the proteasome rather than on direct effects on DNA methylation. Indeed, changes in transcription are not globally correlated with changes in DNA methylation. However, several of the enriched GO-terms we observed are also observed upon MG132 treatment, suggesting, rather, a link through proteasome inhibition. Surprisingly, we observed an enrichment for bacterial defense-related genes in our set of upregulated DEGs. This is rather counterintuitive, knowing that SylA acts as a known virulence factor (Schellenberg et al. 2010; Wäspi et al. 2001). With our current analysis, it is difficult to assess whether the upregulated disease-related genes corresponding to the GO term GO:0042742 have a positive or negative effect on plant defense. One possibility could be that, overall, the strong induction of these genes might rather be detrimental for the plant defense.

Beyond SylA treatment, dynamic changes in centromeric DNA methylation seem to also occur during Arabidopsis infection by P. syringae pv. tomato, which in contrast to P. syringae pv. syringae does not secrete SylA. Changes in DNA methylation levels during P. syringae pv. tomato infection were previously reported (Dowen et al. 2012; Pavet et al. 2006). Pavet et al. (2006) showed that DNA methylation decreases in centromeric regions during early infection (1 day postinfection [dpi]) with P. syringae pv. tomato, whereas Dowen et al. (2012) report an increase of centromeric DNA methylation at 5 dpi. At 5 dpi, similarly to SylA-treated plantlets, the increase of DNA methylation was clear in the CHH and less-pronounced in CG and CHG contexts. Despite not secreting SylA, P. syringae pv. tomato also uses type III effectors, such as HopM1, to inhibit the proteasome during infection (Üstün et al. 2016). It is, therefore, likely that at least parts of the DNA methylation changes observed during P. syringae pv. tomato infection at 5 days after infection could be linked to proteasome inhibition by type III effectors, as we have observed using SylA. Considering the temporal effect of the infection on centromeric DNA methylation, one may speculate that an early response of plants to infection triggers DNA demethylation, as suggested by centromeric demethylation at 1 dpi and the preferential reactivation and, thus, upregulation of TEs observed upon flagellin treatment (Yu et al. 2013). This DNA demethylation has been associated with the expression of nucleotide binding site-leucine rich repeat (NBS-LRR) genes, which are often in close proximity to TEs. At later stages of infection, however, the action of bacterial effectors such as proteasome inhibitors might counteract this DNA demethylation and induce centromeric hypermethylation.

Interestingly, in contrast to the changes in DNA methylation mainly occurring in centromeric and pericentromeric regions, we observed the majority of differentially expressed TEs and DEGs to be located on chromosomal arms. If linked, such phenomena would imply a trans effect of the centromere DNA methylation or organization on the expression of genes situated on the chromosomal arms or both. Such a potential trans effect has been noted previously in response to P. syringae pv. tomato infection (Cambiagno et al. 2018). Indeed, loss of centromeric and pericentromeric DNA methylation in early infection has been hypothesized to be linked to the activation of NBS-LRR disease resistance genes situated on chromosomal arms (Cambiagno et al. 2018). It was hypothesized that centromeric hypomethylation would lead to a recruitment of the RdDM RNA machinery to the centromere to remethylate centromeric TE-dense regions. As a consequence, the RdDM machinery would be depleted from the chromosomal arm, facilitating pattern recognition receptor and NBS-LRR receptor gene expression. Another hint for a trans effect during infection came from epigenetic quantitative trait loci (QTL) linked with increased resistance to the fungi Hyaloperonospora spp. These epigenetic QTL are located in pericentromeric regions but are associated to increased priming of resistance genes situated on chromosomal arms (Furci et al. 2019). It was hypothesized that this trans effect might be due to changes in the three-dimensional organization of the genome that would affect transcriptional activity.

However, to date, it is still not clear how either centromeric and pericentromeric DNA methylation, organization, or both would affect the expression of resistance genes in trans. In the present study, we observe an accumulation of CHH DMRs at the pericentromeric region upon SylA treatment and a change in defense gene expression on the chromosomal arm. We do not observe a significant overlap between the differential sRNA loci and the DEGs at our timepoint, despite allowing for two mismatches. Indeed, it was shown that DNA methylation could be induced even in the presence of mismatches (Fei et al. 2021; Long et al. 2021). As hinted previously, timing of observation could be an important factor explaining the lack of overlap. Indeed, the increase of centromeric 24-nt sRNAs could be a transient but essential process to change pericentromeric CHH methylation patterns and defense gene expression. The trans effect observed would then be linked to the temporal action of sRNAs. An alternative hypothesis to explain the trans effect observed would be, as previously mentioned, a change in nuclear chromatin architecture or of epigenetic marks, such as histone modifications. Indeed, this change in pericentromeric DNA methylation could only be the tip of the iceberg, and infection could have a more profound effect on nuclear chromatin organization and modification. Further studies will be required to evaluate if this trans effect is causative or just correlative and to investigate the potential molecular mechanisms of this trans effect.

Here, we have shown that SylA treatment induces the expression of atypical RdDM components, namely, DRM1, AGO9, and AGO3. DRM1 was previously shown to be mainly active during sexual reproduction (Jullien et al. 2012). Being preferentially involved in CHH methylation, the significant upregulation of DRM1 may potentially explain the high abundance of CHH DMRs compared with the two other sequence contexts. Interestingly, AGO3 and AGO9 were also characterized for their specificity to the reproduction phase of the Arabidopsis life cycle (Jullien et al. 2020; Olmedo-Monfil et al. 2010). So far, their role during bacterial infection remains to be investigated. The variation in genome-wide DNA methylation pattern upon proteasome inhibition by SylA might be linked to the upregulation of these atypical RdDM components. We could neither confirm nor reject this hypothesis using McrBC-qPCR analysis. Looking at other factors involved in RdDM and histone modification (Supplementary Fig. S11A to D), we could also observe an upregulation of JMJ14, SUVH8, and MEA. It is interesting that all those factors, apart from JMJ14, seem specific of the reproductive stages of the Arabidopsis life cycle (Supplementary Fig. S10E). It is therefore tempting to speculate that an epigenetic reprogramming akin to the one happening during sexual reproduction might also happen during bacterial infection and be part of the plant defense mechanism.

Materials and Methods

Plant materials and growth condition

Arabidopsis thaliana wild-type Colombia-0 (Col-0) 1-week-old seedlings were used throughout this study. Seeds were obtained from the Nottingham Arabidopsis Stock Centre (number N22681), were stratified at 4°C for 2 to 3 days, and were germinated in a sterile Petri dish containing half Murashige and Skoog (1/2 MS) basal medium (0.8% micro agar and 0.215% MS in pH 5.7 Milli-Q H₂O). Seedlings were grown in the growth chamber under long-day conditions of 15 h of light at 25°C and 60% humidity and 10 h of night at 21°C and 75% humidity.

Treatment with SylA

Three to five 1-week-old plantlets were transferred in a 24-well plate with 500 μl of liquid 1/2 MS. The plate was then returned to the growth chamber and was placed under agitation. Seedlings were allowed to recover for 24 h. After recovery, the seedlings were treated with SylA at a final concentration of 20 μM for 24 h. SylA stock solution (10 mM in water) was stored at −20°C. After 24 h, seedlings were collected, were immediately frozen in liquid nitrogen, and were used for subsequent extraction. Purified SylA was obtained from the R. Dudler laboratory of the University of Zürich (Waspi et al. 1999, 2001).

Sample preparation, qPCR, and sequencing

Total RNA was extracted from seedlings using QIAzol lysis reagent (Qiagen). All samples were treated with DNase I (ThermoScientific) at 37°C for 30 min. DNAse I was subsequently inactivated by the addition of EDTA and heat treatment (65°C for 10 min). First-strand cDNA was synthesized using 1 μg of DNase-treated total RNA and the Maxima first strand cDNA synthesis kit (ThermoScientific), containing both oligo-dT and random hexamer primers. qPCR tests were performed with a QuantStudio 5 (ThermoScientific), using SYBR green (KAPA SYBR FAST qPCR master mix). qPCR mixes were prepared according to manufacturer protocol (KAPA Biosystems). An RNA equivalent of 25 ng of the cDNA templates was distributed for each reaction. The qPCR program was as follows: 95°C for 3 min, followed by 45 cycles of 95°C for 5 s and 60°C for 30 s. ACTIN2 (AT3G18780) expression was used to normalize the transcript level in each sample. For each condition, RNA abundance of target genes was calculated from the average of three independent biological replicates with three qPCR technical replicates. Real-time PCR primers used in this study are listed in Supplementary Table S3. Expression pattern snapshots from Supplementary Figure S11 were obtained on eBar using data from Klepikova et al. (2016) (Fucile et al. 2011).

For genome-wide bisulfite sequencing, genomic seedling DNA was extracted using a DNeasy plant mini kit and was subsequently processed and sequenced by Novogene (https://en.novogene.com/). Total sRNA was trizol-extracted and was then processed into sequencing libraries and sequenced by Fasteris (http://www.fasteris.com). Total RNAs were extracted and DNAse I–treated as previously mentioned. mRNA libraries were prepared and sequenced by Novogene.

Bioinformatics

mRNA profiling

Paired-end raw mRNA sequencing reads from two biological replicates per treatment regime were aligned to the Arabidopsis thaliana TAIR10 genome assembly using HISAT2 (Kim et al. 2019), using default settings but setting maximal intron length at 10 kb. Aligned reads were sorted using the samtools sort command (Li et al. 2009). Aligned and sorted mRNA sequencing reads were mapped to genomic features using featureCounts (Liao et al. 2014). Multi-mapping reads were counted by assigning fractional counts to the respective features (“–fraction” option). Feature positional information was provided by a custom-made SAF file containing gene (and corresponding exon), transposable_element_gene, transposable_element, pseudogene, miRNA, ncRNA, snoRNA, tRNA, rRNA features retrieved from TAIR10 gff annotation (TAIR10_GFF3_genes_transposons.gff obtained from www.arabidopsis.org). Differential expression has been analyzed using the R package edgeR (Robinson et al. 2010), using the exactTest function. Features with less than 15 reads across all samples were omitted. Differentially expressed features were defined by an absolute log fold change >2 and an adjusted P value < 0.01. All plots have been generated using R-base functions. GO enrichment analysis was performed using the R package topGO (Alexa et al. 2006), using the GO term dataset “athaliana_eg_gene” retrieved from EnsemblPlants (www.plants.ensembl.org). GO enrichment has been performed separately for up- and downregulated DEGs.

sRNA profiling

Raw sRNA sequencing reads were trimmed and were subsequently filtered to remove reads mapping two ribosomal DNA loci (Chr2:1-10000 and Chr3:14194000-14204000), which have previously been repeatedly observed as a source of high numbers of sRNA in Arabidopsis thaliana seedlings and, thus, compromise later analysis. Filtered sRNA reads were aligned to the Arabidopsis TAIR10 genome assembly using bowtie (Langmead et al. 2009) with the following options: -a –best –strata -m 10000. This setting includes two mismatches over a 28-nt seed sequence, chooses only alignments with the minimal number of mismatches, and allows up to 10,000 multiple alignments. To minimize the occurrence of false negatives, no normalization was performed on multi-mapping reads; hence, one read could map to multiple genomic locations and each of these alignments was counted as one (for alignment statistics are provided in Supplementary Table S4). The aligned and sorted sRNA reads were split by size using a custom AWK script and were subsequently mapped to 100-bp genomic bins using HiCdat (Schmid et al. 2015). Only genomic bins, in which at least one of the samples had ≥10 reads were kept for further analysis. Both the control and SylA had near identical effective library sizes (44,810,122 and 44,827,416, respectively, for sRNA between 17 and 30 nt), thus no further normalization was performed. Genomic bins, which differed by more than fourfold counts were considered differential sRNA bins for further analysis. To assign specific genomic features to differential sRNA bins, all non-ambiguous gff annotations were extended by 1 kb up- and downstream of the start and end site and sRNA bins falling within these intervals were associated with the respective annotation unit. The 5′ nucleotide identity was assessed by a custom AWK script, taking into account all aligned sRNA reads of a given length within the respective 100-bp sRNA bins, whereas the sRNA read start had to lie within the bin. To associate differential sRNA bins with specific feature types, we employed bedtools intersect, using custom-made bed files for different feature types, including promoter (defined as 1 kb upstream of gene annotation TSSs), intergenic (excluding promoters), transfer RNA, miRNAs, pseudogenes, TEs (combination of transposable_element and transposable_element_gene, UTR (excluding promoter), intron, coding sequences (gene coordinates – introns). The initial coordinates for the features were retrieved from the publicly available TAIR10 gff file (described above).

Cytosine methylation analysis

Illumina sequencing from bisulfite-treated genomic DNA was aligned using Bismark (Krueger and Andrews 2011) with the following parameters: -q –score-min L,0,-0.2 –ignore-quals –no-mixed –no-discordant –dovetail –maxins 500. Sorted alignment files were analyzed for their cytosine methylation levels, employing the R package methylKit (Akalin et al. 2012). Thereby, aligned reads were processed with the processBismarkAln() function setting the minimal coverage at 10 and minimal quality score at 20. Subsequently, methylation information per 100-bp genomic bin was generated using the tileMethyCounts() function. DMRs were extracted using a q value threshold of 0.05 and methylation difference cutoff of 10%. For further statistical analyses, methylated regions were defined with a cutoff of 40% cytosine methylation for CG context, 20% for CHG context, and 10% for CHH context. To associate a specific genomic feature (e.g., gene identifier) with a DMR, the DMR had to be positioned within a range of 1 kb 5′ of the annotation start to 1 kb 3′ of the annotation end. To associate DMRs with specific feature types, we employed bedtools intersect, using custom-made bed files for different feature types previously used to analyze feature type distribution in differential sRNA bins.

Acknowledgments

We thank the following people for their help: J. Sekulovski for support concerning plant growth, G. Schott and O. Voinnet for technical and financial support in the early stage of this study, R. Dudler and Z. Hasenkamp for providing syringolin A, and U. Grossniklaus for hosting S. Grob's research group.