Introduction

Pseudomonas syringae inflicts leaf spots and cankers on plants globally, and its adequate control poses significant challenges. More than 60 pathovars have been identified, which infect almost all economically important crops in the world (1, 2). P. syringae pv. phaseolicola 1448A (Psph), one of the classical strains of P. syringae, can cause severe halo blight of common beans (Phaseolus vulgaris), which raises it from a common pathogen to a molecular plant-pathogen bacterium(3). P. syringae pv. tomato DC3000 (Pst) and P. syringae pv. syringae B728a (Pss) are two other model strains whose natural host plants are tomatoes and beans, respectively (4, 5). The most important weapon of P. syringae, and the first to be characterised, is the Type III secretion system (T3SS) encoded by hrp and hrc gene clusters, which are flanked by the conserved effector locus and deliver T3 effectors into host cells (6).

The expression of T3SS genes is inhibited in nutrient media such as King’s B (KB), whereas it is induced in minimal medium or plant cells (79). The HrpRSL pathway is the primary regulator of the T3SS in P. syringae, with HrpRS activating hrpL expression. HrpL then binds to the hrp box to induce the translation of downstream T3 effectors (1012). HrpRSL also impacts other virulence-related mechanisms, including motility, biofilm formation, siderophore production, and oxidative stress resistance, which are under the control of complicated regulatory networks, allowing P. syringae to infect plants (1319). However, the methylome and function of DNA methylation in P. syringae pathogenesis and metabolism remain largely unknown.

In bacteria, DNA methylation is the primary level of epigenetic regulation because these prokaryotes lack the histones and nucleosomes of eukaryotic cells. Bacterial DNA has three primary forms of methylation: N6-methyladenine (6mA), N4-methylcytosine (4mC), and N5-methylcytosine (5mC). The latter is the most common type in eukaryotes, while 6mA is the dominant form in prokaryotes. DNA methylation in bacteria is mainly the product of methyltransferase (MTase) enzymes, which transfer a methyl group from S-adenosyl-l-methionine to various positions on target bases, depending on the modification (20). MTases originate from the restriction-modification (R-M) system, which protects bacterial cells from invading DNA by recognising and cleaving specific unmethylated motifs (21, 22). There are three main types of R-M systems, Types I, II, and III, categorised according to the related subunits and the precise site of DNA restriction (23, 24). Additionally, orphan MTases are an emerging group of MTases without cognate restriction that are involved in the regulation of DNA replication and gene expression (21). The Type I R-M system contains three host-specificity determinant (hsd) subunits, restriction (R), modification (M), and specificity (S), which are encoded by hsdR, hsdM, and hsdS, respectively. Genome sequencing and bioinformatics analyses have revealed that the HsdMSR system exists in around half of all bacteria and archaea species (24). Nonetheless, the biological roles and specific motifs alongside their targets of most MTases remain unknown, especially those related to 6mA (25).

A new sequencing technology, single-molecule real-time (SMRT) sequencing, can detect all three forms of DNA methylation, but particularly 6mA (26). SMRT-seq has been applied for complete genome sequencing and insertion sequence profiling in different P. syringae pv. actinidiae (Psa) strains (2729). In P. aeruginosa, Type I R-M systems, along with their specific motifs, have been identified in the model strain PAO1 and two clinical strains via SMRT-seq. These MTases are considered to play important roles in P. aeruginosa virulence or drug resistance (3032).

To study the targets and functions of DNA modifications in P. syringae, SMRT-seq was used to identify global methylation sites and reveal their conserved and divergent functions in the three model strains in this study. We found that HsdMSR in Psph was involved in modulating three pathways, namely T3SS, biofilm production, and the metabolism-related gene expression of ribosomal proteins. Moreover, we found that HsdMSR-mediated transcriptional regulation depends on the full methylation of both DNA strands. Taken together, our results provide insights into the involvement of DNA methylation in the phenotypic traits of P. syringae and transcription regulatory mechanisms that may apply to other pathogenic bacteria.

Results

Genome-wide methylome profiling of three model P. syringae strains

According to REBASE (33) database prediction, 0 to 2 Type I R-M systems and 3 to 4 Type II R-M systems exist in the three pathovars of P. syringae, indicating that Type I is more conserved than Type II in this pathogen (Table S1). To characterise their methylomes, we performed SMRT-seq to obtain the methylation patterns on a genome-wide scale using DNA extracted from the stationary phase of the three wild-type (WT) strains. Our results revealed a total of 10,302 modified bases, of which 3,849, 4,646, and 1,807 nucleotides were significantly identified as being 6mA, 5mC, and 4mC modified, with an interpulse duration ratio (IPD) of > 1.5 throughout the genome of Psph (Fig. 1A). After identifying the significantly modified sites, we orientated the gene locus tags within the methylated bases. Most genes harboured only 6mA modifications, although 339 genes harboured all three DNA methylation types (Fig. 1B).

Genome-wide identification of P. syringae DNA methylation.

(A) The circle map displays the distribution of 6mA, 5mC, and 5mC in Psph WT. (B) The Venn plot reveals overlapped genes within three types of DNA methylation of Psph WT. (C) The circle map displays the distribution of 6mA, 5mC, and 5mC in Pst WT. (D) The Venn plot reveals overlapped genes within three types of DNA methylation of Pst WT. (E) The circle map displays the distribution of 6mA, 5mC, and 5mC in Pss WT. (F) The Venn plot reveals overlapped genes within three types of DNA methylation of Pss WT.

In Pst, a total of 6,002, 2,242, and 3,514 methylation sites were significantly identified as 6mA, 4mC, and 5mC, respectively, using the same cut-off (IPD >1.5) (Fig. 1C). Among these modified bases, a higher degree of overlap was observed across all three types of methylation than in Psph, with 591 (17.7%) genes exhibiting modifications in their sequences (Fig. 1D). Additionally, the methylome atlas of Pss revealed a lower incidence of methylation than those of Psph and Pst, particularly in terms of 6mA modifications, which were only seen in 457 significant occurrences (Fig. 1E). Comparative analysis of the three methylations in Pss showed that the highest numbers of genes were modified by 4mC and 5mC simultaneously (Fig. 1F).

Characteristics of modification loci and GC contents in P. syringae

To further uncover the distribution of the modifications throughout the genome, we calculated how many of the three kinds of modification sites were in these three strains. The majority of modifications occurred in coding sequence (CDS) regions, accounting for at least 80% (Fig. S1A), and less than 20% were in intergenic regions and non-coding RNA (tRNA and rRNA). However, compared with 4mC and 5mC, 6mA was more likely to be in intergenic regions, which suggests that it potentially functions in transcription regulation in P. syringae.

It is known that 5mC occurs within CpG sites in CpG islands, which contain higher GC percentages than is typical in the human genome (34). Additionally, the GC architecture can influence DNA methylation in eukaryotic cells (35). To determine the GC content characteristics of modification sites in P. syringae, we extracted the sequences 50-bp upstream and downstream from the modified bases and calculated their GC content. The distribution of modification sites and their GC content are shown in density plots. Compared with 4mC and 5mC, the GC contents of 6mA sites were the lowest and the closest to the average GC percentage of P. syringae (58% for Psph and Pst, 59.2% for Pss) (Fig. S1B). In contrast, 5mC modification bases had the highest GC content, especially in Pss (Fig. S1C). Furthermore, 4mC sites in Psph had a higher GC content than those in Psa and Pss (Fig. S1D). Similar phenomena have been observed in various other bacterial species. For instance, the Escherichia coli MTase Dcm can catalyse the 5′CCWGG3′ motif (36). Furthermore, in Spiroplasma sp. strain MQ-1, more than 95% of 5mC modifications are found in high-CG sequences, which is similar to the pattern in eukaryotes (37).

Methylated genes revealed highly functional conservation among three P. syringae strains

When we compared the three methylation types in the three tested strains, we noticed different levels of similarity among them (Fig. 2A-B and Table S2). Notably, about 25% to 45% of methylated genes were conserved in more than two of the strains. Interestingly, 5mC had the highest conservation level (39.1%), which might be explained by the rare occurrence of 6mA in Pss. Conversely, an obviously higher degree of conservation was observed in virulence-related genes, including T3SS and alginate biosynthesis-related genes, between Psph and Pst (Table S2). Additionally, some methylated genes (n = 739) harboured sites for different modification types in the three strains (Fig. 2A). For example, the modification sites of a Cro/CI family transcriptional factor (TF) in PSPPH_1319 were for different modification types in the three P. syringae strains. Psph carried 6mA and 5mC, while Pst and Pss had 5mC and 4mC. This difference might result from differences in the MTases and their specific motifs (Table S2).

Functional enrichment analysis of methylation sites in three P. syringae strains.

(A) Repartition of the total pool of modified genes among strains. (B) Proportion of methylated genes detected in one, two, or three genomes for all P. syringae strains and conserved DNA methylation sites with detected genes. (C) The dot plot reveals the functional annotation enrichment of genes within different modified bases among three strains in GO and KEGG terms. The size of the dots indicates the number of related genes.

To elucidate the functions of genes methylated with these three modification types, functional enrichment analyses were performed based on gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (38, 39). Analysis revealed the shared functional characteristics of genes with 5mC and 4mC modifications in Pss and Pst. In contrast, genes with 6mA modifications exhibited greater conservation between Psph and Pst (Fig. 2C). For instance, Pst and Pss contained 4mC-methylated genes enriched in signalling transduction and TF binding, whereas Psph exhibited enrichment in oxidation–reduction process and nucleic acid binding genes. Remarkably, significant conservation of functional terms was observed for genes with 6mA modifications, with at least 21 terms conserved between Psph and Pst. These terms included ATP binding, catalytic activity, integral component of membrane, transferase activity, and phosphorylation. We also conducted Cluster of Orthologous Group (COG) protein function analysis (40, 41), which assigned the methylated genes from the three strains to 20 categories with diverse functions (Fig. S2A-C). In Psph, the top three COG classifications encompassed inorganic ion transport, amino acid transport, and transcription (Fig. S2A). Furthermore, the abundance of modified genes related to cell wall/membrane/envelope biogenesis (M category) in Pss was lower than that of Psph and Pst. Despite the slight variations in COG distribution, the overall pattern exhibited considerable similarity and conservation across the three strains.

Five conserved sequence motifs of 6mA and 4mC were identified in P. syringae

Obtaining the whole methylomes of the three strains led us to investigate the associated MTases and their target motifs. The PacBio motif finder and MEME suite (42) were used to determine the specific sequence motifs within 50 bp upstream and downstream of the significantly modified bases. Using the results of both software packages with corresponding cut-offs, we identified two 6mA motifs in Psph: CAGCN(6)CTC and RAGTACTY (Fig. 3A). The bolded bases indicate the methylation sites in these motifs. The two motifs occurred a total of 2,998 and 300 times in double strands, and more than 98% and 77% of those occurrences were methylated, respectively, revealing the high methylation rates of these two motifs and implying their crucial roles in Psph (Fig. 3B). However, the analysis did not identify any credible motifs for cytosine modifications.

DNA methylation motifs in P. syringae.

(A) 6mA methylation motifs found in Psph using SMRT-seq. The black arrows indicate the site of adenine methylation. (B) Bar plot shows the abundance of methylated numbers throughout all motif sites in Psph. (C) 6mA and 4mC methylation motifs found in Pst using SMRT-seq. The black arrows indicate the site of base methylation. (D) Bar plot shows the abundance of methylated numbers throughout all motif sites in Pst.

Three motifs for 6mA and one motif for 4mC were identified in Pst (Fig. 3C). The first 6mA motif, CAARGAA, was found 2,022 times in the genome, with 2,019 of the occurrences methylated (99.8%) (Fig. 3D). The second 6mA motif, GAAN(4)RTRCC, was fully methylated in both strands. Similarly, the last 6mA motif, GYAGN(5)CTRC, exhibited a high methylation level (96%) in the genome of Pst. Conversely, the only motif identified for 4mC in Pst displayed a considerably lower methylation level than the 6mA sites, with 2,207 occurrences of methylation out of a total of 8,048 sites (approximately 27%). Notably, no credible motifs were found in Pss because of its noticeably lower modification levels compared with the other two strains. Overall, the application of SMRT-seq to the three model P. syringae strains identified specific sequence motifs, including 6mA motifs, exhibiting extensive methylation statuses.

In vivo validation revealed the activity and specificity of Type I R-M system MTase in Psph

The first motif, CAGCN(6)CTC, identified in Psph is similar to the counterpart modified by a Type I R-M MTase in P. aeruginosa (26) (Fig. S3), suggesting that the putative HsdMSR is responsible for the observed motif. As no other motifs exhibited a discernible association with the Type II R-M MTases, we opted to focus our investigation on the putative HsdMSR in Psph. To determine the potential functions of the putative MTases, we constructed MTase knock-out mutants and complementary strains of HsdMSR. Growth curve experiments showed the adverse effects on ΔhsdMSR compared to the WT strain (Fig. S4).

In accordance with previous studies showing that growth conditions affect the bacterial methylation status, we applied dot blot experiments to detect the 6mA levels of Psph in the logarithmic and stationary phases. The results revealed that 6mA levels in the stationary phase were much higher than those in the logarithmic phase (Fig. S5A). To confirm the activity of HsdMSR, a dot blot assay was used to detect the 6mA intensity in WT and ΔhsdMSR. The results shown in Fig. S5B illustrate that 6mA levels were significantly decreased in ΔhsdMSR compared to the WT, but were restored in the complemented strain. In addition, SMRT-seq was performed on ΔhsdMSR, which further confirmed that HsdMSR and its specificity were crucial for maintaining methylation levels in P. syringae. ΔhsdMSR lacked all the 6mA sites within the identified CAGCN(6)CTC motif found in the WT (Fig. S5C). Taken together, the results demonstrated that 6mA levels change with the bacterial growth phase, and HsdMSR is responsible for maintaining 6mA sites within the sequence motif of CAGCN(6)CTC.

Transcriptomic analysis profiling of differentially regulated genes associated with virulence and metabolism in the HsdMSR mutant

To explore the regulatory influence of HsdMSR on gene expression, RNA-seq was applied to analyse WT and ΔhsdMSR in the stationary growth phase, which ensured that the strains had sufficient DNA methylation levels. Differentially expressed genes (DEGs) were obtained with the cut-off of |log2FC| > 1 and an adjusted p-value < 0.05. We identified 395 DEGs between WT and ΔhsdMSR under these experimental conditions. Among these, 218 genes were upregulated, while 177 genes were downregulated compared with the WT (Fig. 4A and Table S3). To investigate the functional differences between WT and ΔhsdMSR, GO and KEGG databases were used to perform functional enrichment analyses based on the DEGs. When using adjusted p-value < 0.05 as a significant cut-off, the upregulated functional terms included alginate biosynthesis, fructose and mannose metabolism, and the oxidation– reduction process (Fig. 4B). The downregulated pathways included the citrate cycle, oxidative phosphorylation, ribosome structure, and translation (Fig. 4B).

Transcriptional changes profiling of hsdMSR mutant in Psph.

(A) The volcano plot reveals DEGs between Psph WT and ΔhsdMSR. The DEGs were (|log2 FC| > 1 and adjusted P value < 0.05) in blue (Ribosomal protein), pink (Alginate biosynthesis), and yellow (T3SS). Each dot represents one gene. (B) Bubble plot shows enriched GO (Red) and KEGG (Blue) terms based on the DEGs between Psph WT and ΔhsdMSR. The x-axis shows the significance of functional annotation terms (−log10 adjusted P-value), and the y-axis indicates the Z-score of terms. The bubble size represents the gene number of terms. (C) COG classification and distribution of up-regulated DEGs. (D) COG classification and distribution of down-regulated DEGs. COG terms are highlighted in different colors. (E) Venn plot reveals the overlapped genes between DEGs and genes within the HsdMSR motif.

We also analysed the function alteration via COG category analysis of DEGs, which revealed that the upregulated genes were more likely to be involved in lanes E (amino acid transport and metabolism), G (carbohydrate transport and metabolism), and I (lipid transport and metabolism) (Fig. 4C). Downregulated DEGs were significantly involved in lanes F (nucleotide transport and metabolism) and J (translation, ribosomal structure, and biogenesis) (Fig. 4D). To investigate the correlation between DEGs and HsdM-recognising motif sites, we performed comparative analysis to reveal 76 genes that overlapped (Fig. 4E). Eight DEGs harboured the HsdMSR methylation motif at the putative promoter region (within 100-bp upstream of the gene) (Table S4).

HsdMSR was required for T3SS and biofilm formation in P. syringae

To confirm the RNA-seq results, we performed quantitative real-time PCR (RT-qPCR) experiments on T3SS genes. Several T3SS effector genes were identified as upregulated, including hopAE1, hrpF, hrpA2, hrpK1, and hrpW1. RT-qPCR was applied under the same conditions as RNA-seq, and the transcriptional levels of all genes were confirmed as being significantly higher in the ΔhsdMSR than the WT (Fig. 5A). These gene expression levels were restored in the complemented strain. To determine whether the loss of HsdMSR affected the Psph virulence phenotypes because of alterations to T3 effectors, we allowed the WT, ΔhsdMSR, and complemented strains to infiltrate the primary leaves of bean plants. At 6 days post-inoculation, ΔhsdMSR was observed to induce more severe symptoms than the WT and complemented strains (Fig. 5B). Overall, HsdMSR inhibited the expression of T3 effectors under KB conditions, and the loss of this modification system enhanced the pathogenicity of the strain during plant infection.

Influence of HsdMSR on virulence and metabolism in Psph.

(A) HsdMSR negatively regulated T3SS-related genes. Data are shown as means ± SD (n = 3). (B) Disease symptoms caused by Psph WT, ΔhsdMSR, and the complemented strains, photographing 6 days after inoculation of 105 CFU/mL bacteria. (C) HsdMSR negatively regulated alginate biosynthesis-related genes. Data are shown as means ± SD (n = 3). (D) The quantification of biofilm production in the Psph WT, ΔhsdMSR, and the complemented strains using a crystal violate staining assay. (E) HsdMSR positively regulated ribosomal protein-related genes. Data are shown as means ± SD (n = 3). (F) The scatterplot shows the TE and mRNA change between Psph WT and ΔhsdMSR. The x-axis presents the log2FC of the mRNA level, and the y-axis shows the log2FC of the TE. The greater TE is represented in red, whereas the lesser TE is displayed in blue.

Besides T3SS, alginate biosynthesis-related genes were also observed among the DEGs between WT and ΔhsdMSR. Alginate biosynthesis proteins have been demonstrated to be essential for Extracellular polymeric substances (EPS) production, which enhances biofilm formation during the bacterial infection process (43, 44). RT-qPCR experiments were performed on the relevant genes, including alg8, algE, alg44, algF, and algD. The expression levels of these genes in ΔhsdMSR were significantly increased compared with the WT and complemented strain (Fig. 5C). To further investigate the influence of HsdMSR on biofilm formation, we performed a crystal violet staining assay to detect biofilm production by the three strains. The intensity of crystal violet staining of biofilms formed by the WT and complemented strains was significantly less strong than that of ΔhsdMSR cells (Fig. 5D), supporting the role of HsdMSR in controlling biofilm formation. We therefore concluded that HsdMSR regulates the virulence of P. syringae by tuning its expression of T3SS and alginate biosynthesis genes.

HsdMSR regulated ribosomal protein synthesis and translational efficiency

We noticed that HsdMSR was important for bacterial growth (Fig. S4), while the expression of many ribosomal protein genes was reduced with the deletion of hsdMSR (Fig. 5E). Regulation of ribosomal gene expression is essential to the integrity of the ribosome structure, which affects translation (45). Therefore, we hypothesised that the low expression level of ribosome proteins in ΔhsdMSR would result in delayed growth. We selected four genes encoding different ribosome subunits (rplS, rpmD, rpsS, and rpsJ) to verify the RNA-seq results through RT-qPCR. While rplS and rpmD encode 50S ribosomal proteins (L19 and L30), rpsS and rpsJ are responsible for 30S ribosomal proteins (S19 and S10). The RT-qPCR results demonstrated that the mRNA expression levels of these genes were significantly lower in ΔhsdMSR than in WT, while levels in the complemented strain were restored to the WT level (Fig. 5E). Given the important role of ribosomal proteins in translation, we performed ribosome profiling, also termed Ribo-seq, of the Psph WT and ΔhsdMSR strains to detect the influence of HsdMSR on translational efficiency (TE). When we compared the WT and ΔhsdMSR, the TE of 162 genes was significantly different (|log2FC| >1.5) (Fig. 5F). Of these genes, 73 genes had enhanced TE, while 89 had supressed TE (Fig. 5F). Most of these genes were related to transmembrane transport and transcription regulation (Fig. S6A-B). Remarkably, the TE of genes linked to the oxidation-reduction process tended to be supressed (Fig. S6B). These results showed that HsdMSR plays an important role in the regulation of metabolism by promoting ribosomal protein synthesis and modulating TE in Psph.

HsdMSR regulation of hrpF was dependent on the full methylation of both strands

To investigate whether and how HsdMSR directly affects gene expression, we focused on those DEGs whose upstream regions harbour its motif. Interestingly, we noticed a methylation motif in the upstream region of hrpF (encoding the pathogenicity factor HrpF), which increased its expression level in the Psph WT strain (Fig. 5A and 6A). To further explore the effects of HsdMSR on the transcription of hrpF, we constructed a lux-reporter plasmid carrying the motif and extended it by 50 bp, which covered the upstream region of hrpF. After transferring the plasmid into the WT and ΔhsdMSR strains, we detected a significant difference in the expression level of hrpF-lux between these two strains (Fig. 6B). The hrpF-lux signal in ΔhsdMSR increased along with culture time and reached a peak at 24 hours, when the methylation level was also elevated in the late growth stage (Fig. 6C). This result of the lux-reporter assay confirmed the RNA-seq and RT-qPCR results showing that hrpF was remarkably upregulated in ΔhsdMSR compared with the WT strain.

6mA methylation regulates gene transcription based on fully methylated.

(A) Upstream region sequence of hrpF carrying the HsdMSR motif. Adenine methylation is highlighted in red. (B) Constantly lux activity detection of the hrpF between Psph WT and ΔhsdMSR. (C) Lux activity of hrpF between Psph WT and ΔhsdMSR at the time point 24 hours. (D) Lux activity of hrpF between Psph WT and ΔhsdMSR with point mutation. In the HsdMSR motif, an “A” was replaced by a “C”, highlighted in bold and underlined. (E) Lux activity of hrpF between Psph WT and ΔhsdMSR with double point mutations. In the HsdMSR motif, two “A” was replaced by a “C”, highlighted in bold and underlined.

To detect the influence of the 6mA site within the HsdMSR motif on the expression of hrpF, we induced a point mutation on the sense strand in the reporter plasmid to change the 6mA to C (CCGCN6CTC). This resulted in a low hrpF expression level similar to that of the WT sequence carrying the A base (Fig. 6D), indicating that a hemi-methylated pattern is insufficient for transcriptional alteration. To verify our hypothesis, we constructed a reporter with two A mutations (CCGCN6CGC), transferred it to the WT and ΔhsdMSR strains, and detected signal changes. As expected, the reporter carrying the HsdMSR motif without two A bases resulted in a significantly higher signal than the hemi-mutated and non-mutated reporters in the WT (Fig. 6E). Taken together, the results showed that HsdMSR regulation of hrpF transcription was based on both strands being fully methylated.

Discussion

Recent advances in sequencing technology have expanded our understanding of DNA methylation in bacteria. Numerous DNA MTases have been extensively characterised in different pathogenic bacteria, and data on their distribution patterns, recognition motifs, and biological functions have been collated (31, 32, 4649). However, most research has concentrated on animal pathogenic bacteria, with limited insights into phytopathogenic bacteria except for Xanthomonas(5052). Consequently, in this study, we employed SMRT-seq to profile the methylome patterns and specific conserved sequence motifs of three P. syringae strains. Notably, we discovered that the virulence and metabolism of Psph are regulated by Type I R-M MTase HsdM. The nonfunctional nature of the subunit restriction endonuclease HsdR was acknowledged, although its potential functionality relies on the presence of MTase (30).

Comprehensive analysis of the methylome atlas revealed that 6mA is the most prevalent modification type in Psph and Pst, mirroring the trends observed in other bacteria. Conversely, Pss displays the lowest occurrences of 6mA throughout its genome (Fig. 1). This discrepancy might be attributed to the absence of the Type I R-M MTase responsible for 6mA in Pss. Similarly, the presence of two Type I R-M MTases in Pst possibly contribute to it having a higher number of 6mA sites than Psph. Additionally, all three pathovars had a similar pattern, with 4mC being the least frequent modification type. Despite belonging to the same P. syringae species, the three strains tested displayed remarkably divergent methylation patterns, which is reminiscent of the phenomenon observed in Xanthomonas spp. (50).

We further identified two (both 6mA) and four (three 6mA, one 4mC) motifs in Psph and Pst, respectively, but none in Pss (Fig. 3). In fact, more motifs were predicted by the MEME and PacBio motif finders, but we chose only those motifs identified by both algorithms for higher accuracy. We found that the 6mA motifs had extremely high methylation levels, ranging from 77% to 100%, in the stationary phase of P. syringae. We subsequently found that the 6mA methylation levels increased during Psph growth, similar to the observations in P. aeruginosa (30). However, DNA methylation is much more stable under different conditions and growth phases in Helicobacter pylori and Salmonella typhimurium (53, 54).

MTases have been reported to play important roles in bacterial growth; for example, they are involved in cell cycle processes such as chromosomal replication in E. coli and Caulobacter crescentus (5557). DNA methylation can activate DnaA and signal the initiation of chromosome replication to affect the growth of gamma-proteobacteria (58, 59). Previous studies revealed that DNA methylation influences bacterial growth through genes related to diverse carbohydrate transport mechanisms (51). We report that the 6mA modification catalysed by HsdM in Psph positively affects its growth, which can be explained by the downregulation of ribosomal proteins and the alteration of metabolism-related gene TE. However, many MTases are thought to be uncorrelated to growth in bacterial species such as P. aeruginosa and Klebsiella pneumoniae (30, 31, 48).

Many MTases and DNA methylations are involved in bacterial virulence through their modulation of gene expression (30, 31, 48, 49, 51, 52, 6062). In this study, RNA-seq analysis of the HsdMSR mutation strain revealed 395 DEGs. Remarkably, not all DEGs harboured the HsdMSR methylation motif in their promoter region. There were no significant correlations between the HsdM-recognised motif sites in the promoter region and DEGs, even though DNA methylation is generally known to affect gene expression by altering the interactions between DNA and proteins such as TFs, which compete with MTases at specific motif sites and thus influence downstream transcription (63). However, recent studies have revealed that DNA methylation within CDS can also alter gene expression in bacteria; for example, 5mC in CDS can enhance transcription while blocking the transcription of CpG islands in promoters in eukaryotic cells (64, 65). We hypothesise that 6mA in Psph can regulate gene expression directly and indirectly. We confirmed that the 6mA motif in the promoter region of hrpF can directly and negatively regulate gene transcription in a manner dependent on full methylation of the strands (Fig. 6). In addition, DNA methylation can change DNA curvature to decrease the thermodynamic stability of the double helix (66). Therefore, those DEGs carrying modified sites, including alginate biosynthesis-related genes and T3 effectors, are believed to be caused by the alteration of nucleoid topology, as seen in Salmonella and H. pylori (49, 67). Furthermore, the 6mA sites, including those in virulence-related hrc/hrp and alg genes, are conserved in Psph and Pst, implying that similar phenotypes occur in the latter strain. Altogether, 6mA modifications in Psph were observed to act as important epigenetic regulators of gene expression.

P. syringae uses its diverse abilities to infect plants. TFs play important roles in virulence regulation and participate in sensory crosstalk by sensing environmental signals and responding. Previous studies have illustrated these complex networks of TFs in P. syringae (16, 18, 68) ; however, the degree to which genes are regulated at the epigenetic level remains unknown. Here, we have provided novel insights into the roles of 6mA in virulence-related gene regulation in Psph.

As more studies apply SMRT-seq to investigate bacterial methylomes, it has become evident that the epigenetic regulation of gene expression is highly prevalent among bacteria. Understanding the mechanisms through which methylation functions can provide novel insights into how strain-specific epigenetic modifications shape the adaptive responses of bacteria to distinct environmental challenges. As the repertoire of MTases varies among P. syringae strains, this approach will help us to better understand the diversity of DNA methylation and epigenetic patterns among P. syringae species.

Materials and methods

Bacterial strains and culture conditions

The bacterial strains, plasmids, and primers used in this study can be found in Table S5. All three model strains and mutations were cultured at 28 ℃ in KB medium with shaking at 220 rpm or LB agar plates. The E. coli strains were grown at 37 ℃ in LB broth (Luria-Bertani) with shaking at 220 rpm or on LB agar plates. Antibiotics were used at the following concentrations: kanamycin at 50 μg/mL, spectinomycin at 50 μg/mL, and rifampin 25 μg/mL.

DNA extraction and SMRT sequencing

Three WT model strains and HsdMSR knockout Psph strains were cultured to the stationary phase. The genomic DNA of the strains was extracted using a TIANamp Bacteria DNA kit (Tiangen Biotech), using the manufacturer’s standard protocols. The SMRT sequence was performed at Abace Biotechnology company using Pacific Biosciences sequel II and Ile sequencer (PacBio, Menlo Park, CA, USA). SMRT-seq reads were aligned to the genome reference of Psph 1448A ( GCF_000012205.1), Pst DC3000 (GCF_000007805.1), and Pss B728a (GCF_000012245.1), respectively. SMRTLink software v13.0 was used to perform DNA methylation analysis. A modification quality value (QV) score of 50 and 100 was used to call the modified bases A and C, respectively.

Deletion mutant and complemented strains construction

The restriction enzymes digested the pK18mobsacB suicide plasmid(69) listed in Table S5. The upstream (∼1500-bp) and downstream (∼1000-bp) fragments of HsdMSR open reading frame (ORF) were amplified from the Psph genome and digested. Then, the digested upstream and downstream fragments were ligated with T4 DNA ligase (NEB). The ligated fragments were inserted into the digested plasmid using ClonExpress MultiS One Step Cloning Kit (Vazyme). The recombinant plasmids were transformed into the Psph WT strain in the KB plate with rifampin and kanamycin. The single colonies were picked to a sucrose plate and then cultured in both KB with kanamycin/rifampin and KB with rifampin alone. Loss of kanamycin resistance indicated a double crossover. Finally, the deletion strains were further confirmed by qRT-PCR to detect the mRNA level of HsdMSR. For the complemented plasmids, the ORF of HsdMS was amplified by PCR from the Psph genome and cloned into the pHM1 plasmid.

Growth curve measure assay

Overnight cultures of Psph and ΔhsdMSR were diluted to an OD600 of 0.1 in fresh KB liquid medium for use as the inoculum. One hundred μL of the inoculum was aliquoted into a 96-well microtiter plate in triplicate and incubated at 28°C for growth. OD600 values were recorded per 2 h for 8 times for plotting of growth curves.

Dot blotting assay

Dot blotting was performed as previously described(70). Briefly, DNA samples were denatured at 95°C for 10 minutes and cooled down on ice for 3 minutes. Samples were spotted on the nylon membrane and air dry for 5 minutes, followed by heat-crosslink at 80 ℃ for 2 hours. Membranes were blocked in 5% nonfat dry milk in TBST for 1 hour at room temperature, incubated with N6-mA antibodies (1:1000) overnight at 4°C. After 5 washes with TBST, membranes were incubated with HRP-linked secondary anti-rabbit IgG antibody (1:5,000) for 1 hour at room temperature. Signals were detected with ChemiDoc Imaging Systems (Bio-Rad). After imaging, incubate the membrane with methylene blue staining buffer for 15 mins with gentle shaking.

RNA-seq analyses

Psph and ΔhsdMSR strains were first cultured to the stationary phase in KB at 28°C. Cultures were collected, and total RNA was extracted using a RNeasy mini kit (Qiagen) for subsequent purification with DNaseI (NEB) treatment. After removing rRNA by using the MICROBExpress kit (Ambion), mRNA was used to generate the cDNA library according to the NEBNext® UltraTM II RNA Library Prep Kit protocol (Illumina), which was then sequenced using the Illumina HiSeq 2000 system, generating 150-bp paired-end reads. Two biological replications have been performed. For each RNA-seq sample, raw reads were trimmed using trim_galore (version 0.6.7) and mapped to the Psph genome using hisat2 (version 7.5.0)(71). DEGs were identified using DESeq2(72). Function enrichment analysis of DEGs was conducted using the R package clusterProfiler(73).

Ribo-seq library construction and analysis

The construction of the Ribo-seq library followed the previous protocol(74). Overnight cultures of strains were transferred into fresh KB medium and culture for 6 hours, and chloramphenicol was added before centrifugation. The pellet was resuspended in lysis buffer and fast-frozen in liquid nitrogen. MNase then digested the supernatant. Sephacryl S400 MicroSpin columns were used to purify the MNase-digested products. The sRNA was separated using a Zymo RNA kit, and the rRNA was removed. The final library was constructed using the NEBNext Small RNA Library Prep Set. The translational efficiency was calculated by dividing the normalized Ribo-seq counts by the normalized RNA counts.

Real-time quantitative PCR (RT-qPCR) verification

For RT-qPCR, RNA was purified using the RNeasy minikit (Qiagen). The cDNA synthesis was performed using a FastKing RT Kit (Tiangen Biotech). The assay was performed by SuperReal Premix Plus (SYBR Green) Kit (Tiangen Biotech) according to the manufacturer’s instruction. Each sample was repeated thrice with 600 ng cDNA and 16S rRNA as the internal control. The fold change represents the relative expression level of mRNA, which can be estimated by the threshold cycle (Ct) values of 2-(ΔΔCt).

Plant infection assay

The bean (Phaseolus vulgaris cv. Red Kidney) plants were used for the pathogenicity assay. The plant was grown in a greenhouse, as described previously(75). Overnight bacterial cultures were diluted to OD600 = 0.2 × 10−3 and were hand-inoculated into the primary leaves of week-old bean plants for 6 days.

Biofilm formation assay

Biofilm production was detected, as previously reported(76). In brief, overnight bacterial cultures of Psph WT and HsdMSR mutants were diluted into OD600 = 0.1 and transferred to a 10 ml borosilicate tube containing 1 ml KB medium. Then, the cultures grow statically at 28°C for 3-5 days. Then, 0.1% crystal violet was used to stain the biofilm adhered to the tube tightly for 30 min, and other components bound to the tube loosely were washed off with distilled deionized water. The remaining crystal violet was fully dissolved in 1 ml 95%∼100% ethanol with constant shaking. Then, it was transferred to a transparent 96-well plate to measure its optical density at 590 nm using a Biotech microplate reader. The experiment was repeated using three independent bacterial cultures.

Data analysis and statistics

The graphs in this paper were plotted using ggplot2 in R 4.2.0 software and GraphPad Prism 10.0.3 (GraphPad Inc.) Differences between groups were analyzed using students’ two-tailed t-tests. The results of all statistical analyses are shown as mean ± standard deviation (SD). All experiments were repeated independently at least three times with similar results.

Data availability

The Ribo-seq and RNA-seq data were uploaded to the National Center for Biotechnology Information SRA database as part of BioProject PRJNA1055550. Codes were available upon reasonable request.

Acknowledgements

This study was supported by grants from Guangdong Major Project of Basic and Applied Basic Research (2020B0301030005), Shenzhen Science and Technology Fund (JCYJ20210324134000002), the National Natural Science Foundation of China (32172358), General Research Funds of Hong Kong (11103221, 11101722, 11102223), the Sichuan Province Science and Technology Planning Project (2021YFSY0005). The funders had no role in study design, data collection, interpretation, or the decision to submit the work for publication.

Author contribution

Conceptualization, J.H., and X.D; Methodology, J.H.; Validation, J.H., F.C., B.L., Y.Y., and Y.S.; Investigation, J.H., F.C. and B.L.; Writing & Editing, J.H., and X.D; Supervision X.D

Declare of interest

The authors declare no competing interests.

Figure. S1 Distribution patterns of methylation sites in three model strains. (A) Bar plot shows the distribution of modified sites located in different regions of CDS, intergenic region, and non-coding RNA. (B) GC contents distribution of 6mA between three P. syringae strains. (C) GC contents distribution of 5mC between three P. syringae strains. (D) GC contents distribution of 4mC between three P. syringae strains.

Figure. S2 COG analysis of methylation sites in P. syringae. (A) COG classification of three DNA methylations in Psph. (B) COG classification of three DNA methylations in Pst. (C) COG classification of three DNA methylations in Pss.

Figure. S3 Comparative genomic analysis of Type I RM system among bacterial species. Different species are highlighted in different colors. The color key indicates the %identity compared to the MTases in Psph.

Figure. S4 Growth curve of Psph WT, ΔhsdMSR, and complemented strains at 28°C in KB medium.

Figure. S5 Identification of Type I DNA methyltransferase in Psph. (A) Dot blot reveals stationary phase of Psph has higher 6mA methylation level. (B) HsdMSR can catalyze 6mA formation in Psph. (C) SMRT-seq reveals HsdMSR catalyzed the 6mA motif CAGCN(6)CTC.

Figure. S6 Number of genes with significant TE changes. (A) The GO terms for the genes in which TE changed more than 1-FC. (B) The GO terms for the genes in which TE changed less than 1-FC.

Table S1 Restriction-modification systems predicted in P. syringae.

Table S2 Modified genes conserved in three P. syringae strains.

Table S3 DEGs genes between Psph WT and ΔhsdMSR.

Table S4 DEGs carrying HsdMSR motif in their putative promoter regions.

Table S5 Genes with changed TE between Psph WT and ΔhsdMSR.

Table S6 Strains, plasmids, and primers.