Skip to main content
Download PDF

Three putative DNA methyltransferases of Verticillium dahliae differentially contribute to DNA methylation that is dispensable for growth, development and virulence

Epigenetics & Chromatin volume 14, Article number: 21 (2021) Cite this article

Abstract

Background

DNA methylation is an important epigenetic control mechanism that in many fungi is restricted to genomic regions containing transposable elements (TEs). Two DNA methyltransferases, Dim2 and Dnmt5, are known to perform methylation at cytosines in fungi. While most ascomycete fungi encode both Dim2 and Dnmt5, only few functional studies have been performed in species containing both.

Methods

In this study, we report functional analysis of both Dim2 and Dnmt5 in the plant pathogenic fungus Verticillium dahliae.

Results

Our results show that Dim2, but not Dnmt5 or the putative sexual-cycle-related DNA methyltransferase Rid, is responsible for the majority of DNA methylation under the tested conditions. Single or double DNA methyltransferase mutants did not show altered development, virulence, or transcription of genes or TEs. In contrast, Hp1 and Dim5 mutants that are impacted in chromatin-associated processes upstream of DNA methylation are severely affected in development and virulence and display transcriptional reprogramming in specific hypervariable genomic regions (so-called adaptive genomic regions) that contain genes associated with host colonization. As these adaptive genomic regions are largely devoid of DNA methylation and of Hp1- and Dim5-associated heterochromatin, the differential transcription is likely caused by pleiotropic effects rather than by differential DNA methylation.

Conclusion

Overall, our study suggests that Dim2 is the main DNA methyltransferase in V. dahliae and, in conjunction with work on other fungi, is likely the main active DNMT in ascomycetes, irrespective of Dnmt5 presence. We speculate that Dnmt5 and Rid act under specific, presently enigmatic, conditions or, alternatively, act in DNA-associated processes other than DNA methylation.

Background

Transcriptional control is important for regulating developmental processes and environmental responses. In eukaryotes, transcriptional control is achieved through transcription factor-mediated and epigenetic mechanisms, the latter affecting DNA accessibility and altering interactions between DNA and various proteins [1,2,3]. Eukaryotic DNA associates with histone–protein complexes to form nucleosomes that are the main constituents of chromatin, a highly ordered DNA-structure [4]. DNA accessibility for the transcriptional machinery is regulated in part by chemical modifications to histones that can alter chromatin structure or nucleosome positioning, and by direct DNA modifications that can alter transcription factor-binding sites [5]. One such DNA modification is mediated by DNA methyltransferases (DNMT) that covalently add a methyl group to the 5th carbon of a cytosine residue (5-methylcytosine, 5mC) [6]. Cytosine methylation can occur in symmetric CG or CHG genomic contexts, or in the asymmetric CHH genomic context, where H stands for either A, C or T. In general, 5mC occurs more commonly at symmetric sites because maintenance methylation can cause methylation of daughter strands during DNA-replication, whereas asymmetric sites require de novo methylation [7]. In mammals, DNA methylation is largely restricted to CG sites, while plants and fungi show methylation in each of the genomic contexts [8].

Compared to animal and plant genomes, fungi typically have smaller and less complex genomes, and they serve as important eukaryote models for various cellular processes including DNA methylation [9]. Much of the initial research on DNA methylation in fungi was performed in the saprophytic ascomycete fungus Neurospora crassa. In N. crassa, DNA methylation is restricted to transposable elements (TEs) and is dependent on a single DNMT, Deficient In Methylation-2 (Dim2), an ortholog of Human Dnmt1 that performs de novo as well as maintenance methylation [10]. Dim2 operates in a complex with Heterochromatin Protein-1 (Hp1) that recognizes and directs DNA methylation to genomic regions marked by tri-methylation of histone 3 lysine 9 (H3K9me3) that is deposited by the histone methyltransferase Deficient In Methylation-5 (Dim5) [11, 12]. Besides Dim2, N. crassa encodes another DNMT domain-containing protein of the fungal-specific class Dnmt4, named Repeat-Induced Point Mutation (RIP)-Defective (Rid), which is only active during sexual reproduction [13, 14]. However, Rid has not been shown to methylate DNA, but is required for the RIP mechanism that can induce C to T mutations in duplicated genomic regions, including TEs [13, 14]. Similar to N. crassa, the ascomycete plant pathogenic rice blast fungus Magnaporthe oryzae encodes orthologues of Dim2 and Rid. However, in contrast to N. crassa Rid, M. oryzae Rid displays DNA methylation activity, albeit with lower activity than Dim2 [15, 16]. The opportunistic human pathogenic basidiomycete Cryptococcus neoformans encodes neither Dim2 nor Rid, but relies on an ortholog of Human Dnmt5 for DNA methylation [17]. C. neoformans Dnmt5 can methylate DNA through direct binding to H3K9me3 or through association with the Hp1 homolog Swi6 [18]. Additionally, C. neoformans Dnmt5 performs maintenance methylation through association with the 5mC-reader Uhrf1 that recognizes hemi-methylated CG sites [18]. Recent phylogenetic analyses of DNMTs across the fungal kingdom revealed extensive diversity in the DNMT repertoires, with only few (less than 10%) species containing either both Dim2 and Rid, or only Dnmt5, whereas many contain the combination of Dim2, Rid and Dnmt5 [19]. Thus, our knowledge on DNA methylation in fungi has been primarily based on species that are not representative for the typical DNMT repertoire of most fungi.

Verticillium dahliae is a xylem-invading, soil-borne ascomycete fungus that causes Verticillium wilt disease on hundreds of plant species [20, 21]. Sexual reproduction has not been reported for V. dahliae that is presumed to mainly reproduce asexually [22]. Recently, we demonstrated that DNA methylation in V. dahliae requires Hp1 and is restricted to H3K9me3-enriched TEs that localize mainly in evolutionary stable core genomic regions that are typically shared across different V. dahliae strains, including centromere regions [23, 24]. In contrast to stable core regions, genomic regions that are important for adaptation show extensive presence–absence polymorphisms between V. dahliae strains, and are therefore designated as adaptive genomic regions [22, 23, 25,26,27]. Many genes that play critical roles in host colonization reside in adaptive genomic regions [22, 26,27,28]. Adaptive genomic regions are enriched in TEs that typically lack DNA methylation, which corresponds with increased transcriptional activity when compared with TEs in the core genome [23]. Interestingly, the transcriptional activity of TEs seems instrumental for the evolution of adaptive genomic regions [26], indicating that TEs in the core genome may carry DNA methylation to supress their transcriptional activity and to prevent genomic alterations that might reduce fitness. In this study, we investigated the contribution of various putative components of the methylation machinery on the physiology and biology of V. dahliae by performing bisulfite sequencing (BS-seq), transcriptomic analysis (RNA-seq), and functional studies on DNA methylation-associated genes.

Results

The genome of V. dahliae encodes three putative DNA methyltransferases

Putative DNMTs in V. dahliae were identified using homology searches to known fungal DNMTs. We selected representative basidiomycete, ascomycete and phycomycete fungi that were previously shown to have DNA methylation, as well as ascomycete Fusarium species that are related to Verticillium. The predicted proteomes of the selected species were searched with a Hidden Markov Model (HMM) pfam model (PF00145) that is characteristic for Dnmt1, Dim2, Rid and Dnmt5. Whereas N. crassa, M. oryzae and C. neoformans possess either a combination of Dim2 and Rid, or only Dnmt5, our analyses showed that several ascomycete species, including all ten species of the Verticillium genus, encode all three DNMTs (Fig. 1, Additional file 2: Figure S1). Thus, Verticillium spp. encode the most commonly shared DNMT complement as observed in ascomycete fungi [19].

Fig. 1
figure 1

Presence of putative 5mC DNA methyltransferases in various fungi. Phylogenetic tree showing a phycomycete (black line), basidiomycetes (green lines) and ascomycetes (blue lines). Filled circles indicate presence of the corresponding DNA methyltransferase as indentified in Additional file 2: Figure S1

Full size image

DNA methyltransferase mutants are not affected in growth and virulence

As V. dahliae encodes three potential DNA methyltransferases, we sought to determine their activity and impact on development and virulence. To this end, we constructed deletion mutants for each DNMT gene, ∆Dim2, ∆Dnmt5 and ∆Rid, as well as the ∆Dim2∆Dnmt5 double mutant in V. dahliae strain JR2. We furthermore generated the H3K9 histone methyltransferase deletion mutant ∆Dim5 that is H3K9me3 deficient (Additional file 2: Figure S2), and used ∆Hp1 [23], the DNA methyltransferase-complex member that recognizes H3K9me3 [11, 12].

Growth impacts were assessed for each strain under axenic growth as determined by colony size, spore production and morphology. Whereas all DNMT mutants displayed similar growth rates, spore production and colony morphology when compared with the wild-type strain, both ∆Hp1 and ∆Dim5 displayed decreased radial growth when compared with the wild-type and complementation strains (Fig. 2a, c, Additional file 2: Figure S3). However, whereas ∆Hp1 produced statistically significant fewer spores, ∆Dim5 produced similar amounts of spores as wild-type V. dahliae when also considering their respective colony sizes (Fig. 2b). This is likely due to ∆Hp1 growing relatively flat, similar to wild-type V. dahliae, while ∆Dim5 colonies display a severely crinkled surface, leading to an increased surface area on the same area of cultivation medium (Fig. 2c). Both ∆Hp1 and ∆Dim5 displayed reduced pigmentation when compared with wild-type V. dahliae (Fig. 2c).

Fig. 2
figure 2

DNMT mutants of Verticillium dahliae do not show altered growth under axenic conditions, stress, or host colinization. a Radial growth of wild-type and mutants over 12 days, with b number of spores produced per mm2 of colony and c pictures showing representative colony morphology after 12 days of growth. d Colony area of wild-type and mutants subjected to various stress agents, relative to average colony area of wild-type (see Additional file 2: Figure S4). e Representative pictures of infected tomato plants at 21 days after inoculation, with f biomass of wild-type and mutants, relative to wild-type infection. Statistically significant differences from wild-type (Wilcoxon signed rank, p < 0.01) are indicated with asterisks. For a, statistical tests were only performed on colony diameter at 12 dpi

Full size image

The deletion mutants were also assessed for growth under abiotic stress conditions by axenically culturing all the strains at elevated temperature, or in the presence of osmotic, oxidative and genotoxic stress agents. Under these conditions, the DNMT mutants grew similar as the wild-type strain, while both ∆Hp1 and ∆Dim5 displayed reduced growth (Fig. 2d, Additional file 2: Figure S4). Interestingly, however, ∆Hp1 grew similar as the wild-type strain when exposed to the genotoxic compound phleomycin despite its growth retardation under all other conditions tested (Fig. 2d, Additional file 2: Figure S4).

The ability to infect tomato plants was also assessed for all mutants. Tomato plants inoculated with any of the DNMT mutants displayed severe stunting at a level similar to plants inoculated with the wild-type strain (Fig. 2f). Fungal biomass measurements on the infected plants confirmed that fungal colonization by the DNMT mutants was similar to that of the wild-type strain (Fig. 2e). In contrast, ∆Hp1 and ∆Dim5 displayed significantly reduced tomato infection, evidenced by a similar canopy area of plants inoculated with these mutants when compared with mock-inoculated plants, as well as by the finding that inoculated plants contained only low amounts of fungal biomass (Fig. 2e). Arguably, the observation of significantly reduced plant infection for both ∆Hp1 and ∆Dim5 should be attributed to their compromised growth characteristics (Fig. 2a).

Dim2 is the main DNA methyltransferase in V. dahliae

To determine the role of the putative DNMTs in cytosine methylation in V. dahliae, whole-genome bisulfite sequencing was conducted on the wild-type strain, along with the DNMT and Hp1 mutants. We recently reported that wild-type V. dahliae displays relatively low levels of DNA methylation, with an average of ~ 0.4% methylation in CG and CHG context and essentially no DNA methylation in CHH context [23]. DNA methylation in V. dahliae is restricted to particular inactive TEs that locate in condensed, H3K9me3-enriched, chromatin regions in the core genome, including those localized in centromeres (Fig. 3a) [23, 24]. We furthermore showed that the ∆Hp1 mutant lost all DNA methylation, indicating that Hp1 is required for cytosine methylation and V. dahliae DNMTs cannot methylate DNA independently [23].

Fig. 3
figure 3

Dim2 is the main DNA methyltransferase in V. dahliae. a Whole-chromosome plot displaying the fraction of methylated cytosines for non-overlapping 10-kb windows for wild-type, and DNMT and Hp1 deletion mutants with chromosome 5 as an example. Grey boxes, displayed below the DNA methylation tracks, indicate the hypomethylated windows compared to the wild-type strain in CG and CHG context from Table 1. Previously defined adaptive genomic regions [23] are highlighted in yellow. b Overlap of hypomethylated windows in mutant strains showing severe loss of methylation. c Expression (TPM values) of DNA methyltransferase genes Dim2, Dnmt5 and Rid, as well as Hp1 and Dim5 of V. dahliae strain JR2 cultured in Czapec–Dox medium (CZA), half-strength Murashige–Skoog medium (MS) and potato dextrose broth (PDB), and during Arabidopsis infection at 21 days post inoculation (Arab), in triplicates

Full size image

To study the extent to which the different V. dahliae mutants lost DNA methylation, we compared the bisulfite sequencing patterns over the genome in 10-kb windows and assessed the amount and location of hypomethylated windows when compared with the wild-type methylation pattern. Of the DNMT deletion mutants, ∆Dim2 showed considerable loss of cytosine methylation, having 100 and 61 hypomethylated windows in the CG and CHG context, respectively (Table 1). As there is little methylation in CHH context, we combined the methylation data for CG and CHG context and also determined hypomethylation for the contexts simultaneously. This combination optimizes the number of potential methylated cytosines per window and therefore better captures differential methylation. In the combined contexts we observed 97 hypomethylated windows that locate at regions that have relatively high methylation percentages in wild-type V. dahliae (Table 1, Fig. 3, Additional file 2: Figure S8). Notably, additional regions showed reduced DNA methylation in ∆Dim2, yet these were not classified as hypomethylated because the methylation level was already low in the wild-type and therefore did not meet our criteria for calling hypomethylated regions (see “Materials and methods” for details). In contrast to the results for ∆Dim2, ∆Dnmt5 and ∆Rid largely retained DNA methylation with only 3 and 12 windows being hypomethylated in CG context and 15 and 17 windows in CHG context, respectively (Table 1). When assessing CG and CHG methylation combined, ∆Dnmt5 and ∆Rid have one and eight hypomethylated windows, respectively. Additionally, their genome-wide DNA methylation patterns are similar to the wild-type with no obvious loss of DNA methylation peaks (Fig. 3a). Thus, both Dnmt5 and Rid may contribute to DNA methylation in V. dahliae, albeit to a much lower degree as Dim2. The ∆Dim2∆Dnmt5 double mutant as well as ∆Hp1 showed similar cytosine methylation levels over the genome as ∆Dim2 and had hypomethylation of 93 and 99 windows in CG context and 59 and 65 windows in CHG context, respectively, and 89 and 92 hypomethylated windows when combining CG and CHG methylation data (Table 1, Fig. 3a). Loss of methylation in the ∆Dim2, ∆Dim2∆Dnmt5 and ∆Hp1 mutants largely occurred in the same genomic regions, as 84 of the hypomethylated windows were shared between the mutants (Fig. 3b). Even though the chromosome plots of the ∆Dim2∆Dnmt5 double mutant as well as ∆Hp1 are similar to those of ∆Dim2, the few bins with slightly elevated methylation levels in ∆Dim2 have decreased further (Fig. 3a). This finding suggests that Dnmt5 has DNA methylation activity on particular genomic regions, albeit at a lower level. However, ∆Dnmt5 does not display reduced methylation at the regions that remain slightly methylated in ∆Dim2 (Fig. 3a). Thus, if Dnmt5 has DNA methylation activity, it is redundant and secondary to the DNA methylation activity of Dim2. No windows were hypomethylated for CHH in any of the mutants (Table 1, Additional file 2: Figure S7). The few bins with low levels of DNA methylation that locate in adaptive genomic regions behave similar as those in core regions, in that they are hypomethylated in ∆Dim2, ∆Hp1 and ∆Dim2∆Dnmt5 (Fig. 3a). These results show that the methyltransferase Dim2 is responsible for the vast majority of detectable DNA methylation in V. dahliae.

Table 1 Number of 10-kb windows that are hypomethylated in Verticillium dahliae DNMT mutants relative to those in the wild-type strain
Full size table

Previous research shows that methylated TEs have accumulated more C–T mutations than non-methylated TEs [23, 26]. As we observe that methylation mainly occurs in CG and CHG contexts, we investigated whether cytosine methylation is directly involved in C–T mutation. If methylated cytosines are more likely to mutate, we would expect that methylated TEs have specifically lost cytosines in CG and CHG context, while CHH sites remain intact. Unexpectedly, however, we observe that methylated TEs only display significant depletion of CHG sites, whereas CG sites occur as frequently as would be expected based on the sequence composition of the TE (Additional file 2: Figure S9). In contrast, non-methylated TEs have slightly fewer CG sites than expected based on their sequence composition, yet they do not have reduced CHG sites (Additional file 2: Figure S9). For both methylated, as well as non-methylated TEs, we observe that CHH sites occur as frequent as would be expected based on sequence composition. Thus, although methylated TEs show increased C–T mutations [23], these mutations do not affect all methylated cytosines, as they are largely restricted to cytosines in CHG context.

We compared V. dahliae Dnmt5 to the homolog in C. neoformans, where it is the sole active DNA methyltransferase [18]. The two proteins share only 18% sequence similarity, but do share similar domain structures, except that the V. dahliae Dnmt5 lacks the N-terminal chromo-shadow domain found in C. neoformans Dnmt5 (Additional file 2: Figure S10). This domain is responsible for the direct binding to H3K9me3, and this histone mark is required for DNA methylation, which could explain why we observed little Dnmt5 contribution to DNA methylation in V. dahliae. However, C. neoformans Dnmt5 can also bind H3K9me3 indirectly through Hp1 [18], and it is not clear if this is also the case for V. dahliae. The lack of DNA methylation by Dnmt5 cannot be explained by transcriptional activity, as Dnmt5 is expressed higher than Dim2 during cultivation in PDB (Fig. 3b).

Loss of DNA methylation does not affect transcriptional regulation

While the Dim2 mutant loses nearly all DNA methylation (Table 1, Fig. 3a), it displays wild-type-like growth in vitro, under stress conditions as well as during infection (Fig. 2), suggesting that DNA methylation is not essential under these conditions. However, given that DNA methylation is mainly restricted to TEs [23], we anticipated that loss of DNA methylation could result in activated transcription at TEs. To address this, we performed RNA sequencing on axenically grown cultures of all DNMT mutants, as well as the Hp1 and Dim5 mutants. Consistent with the lack of DNA methylation at coding regions, all three single DNMT mutants and the double mutant showed differential expression of only few genes (< 10) (Fig. 4a, Additional file 1: Table S1). Unanticipatedly, the four DNMT mutant strains similarly showed differential expression of only a few TEs (< 10) (Fig. 4b). In contrast, the ∆Hp1 and ∆Dim5 mutant strains showed considerable differential expression of genes and TEs (Fig. 4a). In total, 1661 genes were induced and 663 repressed in ∆Hp1, and 1617 genes were induced and 781 were repressed in ∆Dim5 when compared with wild-type (Fig. 4a). Furthermore, 261 TEs were induced and 23 were repressed in ∆Hp1, whereas 241 TEs were induced and 47 were repressed in ∆Dim5 when compared with wild-type (Fig. 4b). Analysis of the induced genes and TEs revealed a large overlap between the mutants, with 1207 out of 1617 (~ 75%) of the induced genes and 166 out of 241 (~ 69%) of the induced TEs shared between the two mutants (Fig. 4a, b). This overlap is likely related to the functional link between these heterochromatin components, as Hp1 directs DNA methylation to H3K9me3 deposited by Dim5 [11]. To study whether the genes activated in ∆Hp1 and ∆Dim5 represent specific biological functions, we performed GO enrichment analysis on the 1207 induced genes. Interestingly, genes encoding secreted proteins and proteins involved in transport and metabolic processes were overrepresented among the induced genes (Fig. 4c), suggesting that the mutants impact expression of genes with roles in responses to the environment.

Fig. 4
figure 4

Genes and TEs that are induced in the Verticillium dahliae Hp1 and Dim5 mutants do not associate with H3K9me3-marked chromatin. Differentially expressed genes (a) and TEs (b) in the mutants relative to wild-type. Induced genes and TEs are indicated in blue, repressed genes and TEs in red. The number of genes and TEs that are induced and repressed in both the Hp1 and Dim5 mutants are indicated by opaque coloring in black rectangles. c Gene Ontology (GO) terms that are enriched (fold change > 1, p.adj < 0.01) in the set of genes that are induced in both the Hp1 and Dim5 mutant. (d). Whole-chromosome plot displaying the location of induced genes (in blue) and TEs (in red) on chromosome 5 as an example. Clusters of induced genes and TEs are indicated as blue and red rectangles, respectively. H3K9me3-ChIP signal along the chromosome is indicated in green in the upper track. adaptive genomic regions [23] are highlighted in yellow. The minimal distance of genes and TEs to H3K9me3-enriched genomic regions (e) and to adaptive genomic regions (f). Asterisks indicate statistical differences (Wilcoxon signed rank test, p < 0.01) between genes and TEs induced in both the ∆Hp1 and ∆Dim5 mutants and those that are not induced in the mutants

Full size image

As Hp1 functions downstream of Dim5 during DNA methylation, we expect that the induction of genes as well as TEs in both mutants may be due to reduced recruitment of repressive complexes to previously silenced chromatin regions because either H3K9me3 or Hp1 is lacking. Based on H3K9me3 ChIP-seq, we found that approximately 2.1 Mb (~ 6%) of the genome is associated with H3K9me3 that occurs in 621 enriched genomic regions, of which 38 are larger than 10 kb (Additional file 2: Figure S12). To study whether the induced genes and TEs localize in these H3K9me3 domains, we investigated the occurrence of physical clustering of the 1207 induced genes and the 166 induced TEs. Our hypothesis was that Dim5 deposited H3K9me3 and associated Hp1 mediate transcriptional silencing in physical proximity to H3K9me3 domains. As such, we expected that genes and TEs induced in ∆Hp1 and ∆Dim5 occurred in clusters and that these clusters would be in proximity to H3K9me3-enriched genomic regions. We identified 58 clusters containing 526 of the 1207 (~ 44%) induced genes and four clusters containing 37 of the 166 (~ 22%) induced TEs, which is more than expected by chance, as measured from 1000 random sets of 1207 genes (p < 0.001) (Additional file 2: Figure S12) and 166 TEs (p = 0.024) (Additional file 2: Figure S13). H3K9me3 domains contain numerous TE copies (1034 out of 2574, ~ 40%) and only few genes (76 out of 11,426, ~ 0.6%) (Fig. 4d, Additional file 2: Figure S14) [23]. Next, we calculated the distance of each gene and TE to the closest H3K9me3 domain and associated this to induced and non-induced genes and TEs. When considering the smallest distance to H3K9me3 domains, we observed that induced genes are slightly closer to H3K9me3 domains than non-induced genes (Fig. 4e). In contrast, induced TEs are slightly further from H3K9me3 domains than non-induced TEs (Fig. 4e). Additionally, when considering presence of TEs in H3K9me3 domains, we observed that significantly fewer induced TEs, 48 out 165 (29.1%), than non-induced TEs, 987 out of 2409 (41.0%), locate in H3K9me3 domains (Fisher’s exact test, p = 0.0126). The relatively minor enrichment of induced genes near H3K9me3 domains and the enrichment of induced TEs away from H3K9me3 domains suggests that the transcriptional changes in ∆Hp1 and ∆Dim5 are not due to reduced H3K9me3 and Hp1 association. As we demonstrated that the induced genes and TEs occur clustered in the genome (Additional file 2: Figures S12, S13), we asked whether the clusters localize in specific genomic regions. As we found that induced genes are enriched for genes encoding secreted proteins and proteins involved in metabolic processes (Fig. 4c), we speculated that the induced genes may be involved in processes related to plant infection. Such genes are typically located in adaptive genomic regions of V. dahliae, which are enriched in TEs but are not associated with H3K9me3 [22, 23]. Therefore, we also tested whether the induced genes and TEs locate in proximity to adaptive genomic regions. Intriguingly, genes and TEs induced in ∆Hp1 and ∆Dim5 are significantly closer to adaptive genomic regions than non-induced genes and TEs (Fig. 4d, f). Additionally, 180 out of 1207 (14.9%) and 59 out of 165 (35.8%) of the induced genes and TEs locate in adaptive genomic regions, which is significantly more than the 818 out of 10,219 (8.0%) and 474 out of 2409 (19.7%) of non-induced genes (Fisher’s exact test, p < 0.00001) and TEs (Fisher’s exact test, p < 0.00001). Consequently, both genes and TEs that are induced in ∆Hp1 and ∆Dim5 reside significantly closer to adaptive genomic regions than non-induced genes and TEs (Fig. 4f). Since adaptive genomic regions are not associated with H3K9me3, these findings suggest that the observed transcriptional changes are not directly related to loss of Hp1 binding at H3K9me3 domains, but rather through pleiotropic effects affecting transcription throughout the genome, and especially at adaptive genomic regions.

Discussion

DNA methylation is essential for proper functioning of nuclear processes in many organisms [6], but various fungal species have lost or degraded their machinery for DNA methylation [19]. The most commonly found combination of DNMTs in ascomycete genomes is the presence of Dim2, Rid and Dnmt5 [19]. As Dim2 is the main DNA methyltransferase gene in fungal species that lack Dnmt5, and vice versa, it is relevant to study the importance of these DNA methyltransferase genes in fungal species that carry both. The fungal pathogen Z. tritici carries all three DNMTs and, similar to V. dahliae, loss of Dim2 almost completely abolishes DNA methylation [29], indicating that Dnmt5 and Rid have little to no DNA methylation activity. However, a low residual DNA methylation signal remains in Dim2 mutants of V. dahliae and Z. tritici, which may be due to a low degree of Dnmt5 activity (Fig. 3, [29]). Our results indicate that Dnmt5 is more highly expressed during growth in nutrient-limited media, a type of environmental stress. It is possible Dnmt5 may be more active and cause differential DNA methylation during specific growth conditions not tested here, an occurrence that has been observed for DNMTs in several plant and animal species [30, 31]. Differential expression of DNMTs has previously been observed in the entomopathogenic ascomycete fungus Cordyceps militaris that contains orthologues of Dim2 and Rid [32]. Whether V. dahliae Dnmt5 plays such a role requires further study.

The importance of DNA methylation in fungi that are able to perform DNA methylation remains unclear. Deletion of functional components of DNA methylation did not result in clear phenotypic alterations in N. crassa or the necrotrophic plant pathogenic fungus Botrytis cinerea, while deletion of Dim2 in M. oryzae leads to aberrant colony morphology and compromised conidiospore formation [10, 16, 33]. In Z. tritici, strains collected in the center of origin of its wheat host carry DNA methylation, while strains collected in Europe contain mutated Dim2 copies that lack DNA methylation [29]. The Z. tritici strains that lack a functional copy of Dim2 are at least as virulent as strains that perform DNA methylation [34], suggesting that the recent loss of DNA methylation in these Z. tritici strains does not negatively affect their infection biology. Our study reveals that DNA methylation in V. dahliae is not essential for growth and infection. Moreover, we show that loss of DNA methylation does not result in altered expression of genes or TEs, an observation that could be explained by DNA methylation co-localizing with H3K9me3, which is likely sufficient for heterochromatin formation and transcriptional silencing in the absence of DNA methylation. This is further supported as the H3K9me3-deficient V. dahliae Dim5 mutant showed significant differential expression of genes and TEs. Interestingly, some genes and TEs induced in this mutant were located in adaptive genomic regions that are not labeled with H3K9me3 in wild-type V. dahliae, suggesting that the removal of H3K9me3 leads to pleiotropic effects in unrelated genomic regions. Similar effects on gene and TE expression occurs in the V. dahliae Hp1 mutant, indicating that the differential expression observed in the Hp1 and Dim5 mutants are related to disrupted Hp1 functioning. In N. crassa, and also fission yeast Schizosaccharomyces pombe that lacks DNA methylation, Hp1 was found to be involved in the formation of H3K9me3-associated heterochromatin [11, 35, 36]. As such, it is possible that the transcriptional changes in the V. dahliae Hp1 and Dim5 mutants are due to pleiotropic effects of chromatin de-condensation. In line with such pleiotropic effects on chromatin architecture, deletion of Dim5 in Z. tritici, and Dim5 and Hp1 in N. crassa leads to re-localization of H3K27me3 to previous H3K9me3 domains [37, 38]. Previously, we showed that the adaptive genomic regions in V. dahliae are enriched for H3K27me3 [23], suggesting that the observed transcriptional induction of adaptive genomic region-localized genes and transposons in the V. dahliae Hp1 and Dim5 mutants may be due to altered localization patterns of H3K27me3.

Considering that experimental and natural loss of DNA methylation in various fungi does not seem to affect their proliferation, it is remarkable that the vast majority of fungal species have retained DNA methylation. One explanation for the role of DNA methylation in fungi, which accounts for the lack of reported phenotypes, is that it serves in maintaining genome integrity during evolution. In this way, DNA methylation does not functionally regulate transcription per se, but works in conjunction with H3K9me3 to minimize the impact of TEs in the genome. One possible mechanism is that DNA methylation may have persistent effects on TE activity through spontaneous deamination of methylated cytosines, resulting in C to T mutations [39]. The deamination process is considered an important driver of mutations in Z. tritici TEs as recently shown in an experimental evolution experiment in which a DNA methylation competent strain had increased in C to T mutations compared to the strain lacking DNA methylation [29]. Interestingly, in V. dahliae we previously observed that TEs that carry DNA methylation contain more C to T mutations than unmethylated TEs [23]. Typically, such C to T mutations are also caused by the RIP mechanism, which relies on the Rid DNA methyltransferase that is active during sexual cycles in N. crassa [13, 14]. However, since V. dahliae is presumed to reproduce asexually, it may be more likely that C to T mutations in TEs are caused by spontaneous deamination. These results support that the main role for DNA methylation in fungi might be to aid in TE sequence degradation over time, not to directly supress transcriptional activity. Alternatively, it is possible DNA methylation is important for inhibiting transcriptional activity of TEs during specific developmental or cell-cycle stages which have not been reported or observed to date.

Conclusion

Our results show that although V. dahliae encodes multiple DNMTs, only Dim2 seems to be essential for DNA methylation. As Dim2, Dnmt5 and Rid are wide-spread among ascomycetes, it is likely that their combined presence is an ancestral state [19]. Even though only four ascomycete species have been studied with respect to the contribution of their DNMTs to DNA methylation so far, these studies suggest that species, irrespective of the presence or absence of Dnmt5, utilize Dim2 as the main DNMT (this study; [10, 16], Möller et al., 2020). Additional research is needed to determine if Dnmt5 and Rid play a role in DNA methylation, or possibly in other DNA-associated pathways, or if their presence is the remnant of an ancestral state that is not strongly selected against.

Materials and methods

Assessment of DNMT occurrence

To assess the presence of DNA methyltransferases in a selection of fungal species with confirmed DNA methylation performance, we downloaded predicted proteomes of Aspergillus flavus strain NRRL_3357 (AFL2T), Coprinopsis cinerea strain Okayama-7#130 (CC1G), Cryptococcus neoformans strain H99 (CNAG), Fusarium graminearum strain PH-1 (FGSG), F. oxysporum strain 4287 (FOXG), Laccaria bicolor strain S238N-H82 (lacbi2), Magnaporthe oryzae strain MG8 (MGG), Neurospora crassa strain OR74a (NCU), Phycomyces blakesleeanus strain NRRL 1555(-) (Phybl2), Postia placenta strain MAD698 (pospl1), Tuber melanosporum strain Mel28 (Tubme1), Uncinocarpus reesii strain 1704 (URET) and the Verticillium albo-atrum PD747, V. alfalfae PD683, V. dahliae JR2, V. isaacii PD618, V. klebahnii PD401, V. longisporum PD589, V. nonalfalfae T2, V. nubilum 397, V. tricorpus PD593 and V. zaregansianum PD739. The predicted proteomes were scanned for the presence of a DNA methyltransferase domain with hmmsearch with -cut_ga option. Identified proteins were visually inspected. To check for presence of additional not-annotated DNA methyltransferase homologs of Dim-2, Dnmt5 and Rid that were initially not annotated in the predicted proteomes, we manually assessed the genomes using TBLASTN and Augustus. Phylogenetic trees were constructed using IQ-tree.

Fungal growth and mutant generation

V. dahliae strain JR2 (CBS 143773; [40] was maintained on potato dextrose agar (PDA) (Oxoid, Thermo Scientific, CM0139) and grown at 22 °C in the dark. The ∆Dim2, ∆Dnmt5, ∆Rid and ∆Dim5 single deletion mutants and the ∆Dim2-∆Dnmt5 double mutant were constructed as previously described [41]. Briefly, for all genes except Dnmt5, genomic DNA regions flanking the 5′ and 3′ ends of the coding sequences were amplified with PCR using primers listed in Additional file 1: Table S2 and cloned in to the pRF-HU2 vector [42], using USER enzyme following the manufacturer’s protocol (New England Biolabs, MA, USA). For Dnmt5, the 5′ and 3′ amplicons were cloned into vector pRF-NU2, a custom-made pRF-HU2 variant, containing the NAT-cassette for selection on nourseothricin. Sequence-verified vectors were transformed into Agrobacterium tumefaciens strain AGL1 used for V. dahliae conidiospore transformation as described previously [41]. V. dahliae transformants that appeared on hygromycin B or nourseothricin (for Dnmt5) were transferred to fresh PDA supplemented with hygromycin B or nourseothricin after 5 days. Putative transformants were screened using PCR to verify deletion of the target gene sequence (Additional file 1: Table S3) when compared with positive amplification from the wild-type strain. To further confirm integration of the selectable marker at the locus of interest, another round of PCR was conducted in which one primer was position adjacent to the deleted genomic region, and the other primer was designed to bind a portion of the inserted vector DNA (Additional file 1: Table S4). In this manner, deletion mutants were confirmed to lack the gene of interest and contain the selectable marker at the locus of interest. Generation of the Hp1 deletion mutant was conducted in the same way and described previously [23]. Complementation vectors were generated by amplifying the coding region of Dim2, Hp1 and Dim5 from genomic DNA using primers listed in Additional file 1: Table S2, and ligating the amplicons into PacI-digested pFBT-005 vector using the NEBuilder HiFi DNA Assembly Cloning Kit (New England Biolabs, MA, USA). Fungal transformations were performed as described above and obtained colonies were screened by PCR to verify presence of target gene (Additional file 1: Table S3).

Growth and inoculation assays

To check for aberrant growth phenotypes of the generated mutants, all strains were cultured as described above. To this end, conidiospores were harvested in sterile water and brought to a final concentration of 106 conidiospores per mL. Subsequently, 10 µL of conidiospore suspension, containing 104 conidiospores, was deposited in the middle of a 90 mm Petri dish containing 20 mL of PDA. Plates were stored at 22 °C in the dark and colony diameter was measured in perpendicular directions after 4, 8 and 12 days of growth. After 12 days of growth all newly formed conidiospores were harvested in 1 mL of water and counted using a hemocytometer.

Stress assays were performed by spotting 5 µL conidiospore suspension containing 5 × 103 conidiospores on PDA without supplement, or on PDA supplemented with 1 M NaCl, 1 M Sorbitol, 2 mM H2O2 or 10 µg/µL phleomycin, and on PDA adjusted to pH 10. Plates were incubated at 22 °C in the dark, apart from one set of PDA plates without supplement that was incubated at 28 °C to assess heat stress responses. Pictures were taken after 6 or 10 days, depending on wild-type colony development, and colony size was determined using ImageJ software with custom settings for each stress condition.

Infection assays were performed using root dip inoculation in a conidiospore suspension of 106 spores per mL on 10-day-old seedlings of tomato cultivar Moneymaker. Stems of infected plants were harvested at 21 days after inoculation, cut in small pieces, frozen in liquid nitrogen and ground by reciprocal shaking in a MixerMill MM 400 (Retsch, Haan, Germany). DNA was isolated incubating the ground powder with 800 µL of CTAB lysis buffer at 65 °C for 1 h, followed by addition of 400 µL chloroform/IAA (24:1), vigorous shaking and centrifuging for 5 min at ~ 13,000 RCF. DNA was precipitated from the aqueous layer with isopropanol and the precipitate was washed with 70% ethanol. The fungal biomass in the stem tissue was determined with real-time PCR using V. dahliae ITS-specific and tomato GAPDH-specific primer sets (Additional file 1: Table S5).

Bisulfite sequencing and analysis

The V. dahliae wild-type strain, ∆Dim2, ∆Dnmt5, ∆Dim2/∆Dnmt5, ∆Rid and ∆Hp1 were grown in potato dextrose broth (PDB) for 3 days, strained through miracloth (22 μm) (EMD Millipore, Darmstadt, Germany), pressed to remove excess liquid, flash frozen in liquid nitrogen and ground to powder with a mortar and pestle. Genomic DNA was isolated as described above and sent to the Beijing Genome Institute (BGI, Hong Kong, China) for bisulfite conversion, library construction and Illumina sequencing. Briefly, the DNA was sonicated to a fragment range of 100–300 bp, end-repaired and methylated sequencing adapters were ligated to 3′ ends. The EZ DNA Methylation-Gold kit (Zymo Research, CA, USA) was employed according to manufacturer’s guidelines for bisulfite conversion of non-methylated DNA. Lambda DNA was used as spike-in to determine conversion efficiency, which was > 99% for all samples. Libraries were paired-end 100 bp sequenced on an Illumina HiSeq 2000 machine.

Whole-genome bisulfite sequencing reads were analyzed using the BSMAP pipeline (v. 2.73) and methratio script [43]. The results were partitioned into CG, CHG and CHH cytosine sites for analysis. Only cytosine positions containing more than 4 sequencing reads were included for analysis. BSMAP datasets were further analyzed using MethylKit (v. 1.12.0) [44]. Methylation levels were summarized as the number of methylated cytosines divided by the total number of sequenced cytosines per 10-kb window. Hypo-methylated windows in the mutants were determined by comparing corresponding 10-kb windows between mutants and wild-type and selecting windows with meth.diff value < − 1 and a q value < 0.01. Genome plots displaying methylation data were generated using karyoploteR (v. 1.12.4) [45].

Analysis of CG, CHG and CHH site occurrence in methylated and non-methylated TEs

The observed occurrence of CG, CHG and CHH sites were compared to the expected occurrence of these sites for each of the 2574 TEs based on the sequence composition (CG expected: n(C) × n(G)/(N × N − 1), CHG expected: n(C) × n(C + A + T) × n(G)/(N × N − 1 × N − 2), CHH expected: n(C) × n(C + A + T) × n(C + A + T)/(N × N − 1 × N − 2)). TEs were categorized as methylated if any cytosine in the sequence showed evidence for methylation in wild-type JR2.

Protein extraction and western blot

Total proteins were extracted from V. dahliae wild-type, DDim2, DDim5, and DHp1 grown for 10 days in potato dextrose broth at 22 °C. Mycelium was collected by straining over a double layer of miracloth, frozen in liquid nitrogen and ground with a mortar and pestle. Approximately 0.5 g of ground mycelium was resuspended in 12 mL lysis buffer (75 mM Tris–HCl pH 7.4, 0.5 mM EDTA, 0.3 M Sucrose, 40 mM NaHSO3, 10 mM MgSO4, 0.5% NP-, 2 mM Phenylmethanesulfonyl fluoride (PMSF), 100 µM Leupeptin, 1 µg/mL Pepstatin), briefly vortexed and rotated at 4 °C for 15 min. Samples were spun at 4 °C at 10,500g for 20 min, and the pellet was resuspended in 2 mL CW buffer (10 mM Tris–HCl pH 8.0, 150 mM NaCl, 1:400 β-mercaptoethanol, 2 mM Phenylmethanesulfonyl fluoride (PMSF), 100 µM Leupeptin, 1 µg/mL Pepstatin). 0.5 mL of 0.4 M sulfuric acid was added, rotated at 4 °C for 2.5 h and centrifuged at 4 °C at 10,500g for 15 min. Supernatant was added to 25 mL acetone on ice and proteins were precipitated at − 20 °C overnight. Precipitates were collected by centrifugation 4 °C at 7500g for 15 min and resuspended in 300 µL 4 M urea. To assess H3K9 methylation status, approximately 15 uL of total protein was added to 2× Laemmli loading buffer (4% SDS, 20% glycerol, 0.004% bromophenol blue, 125 mM Tris HCL pH 6.8), boiled at 95 °C for 1 min, and separated using PAGE (15% polyacrylamide gel). Proteins were transferred to PVDF membranes, blocked in 5% BSA, washed twice in TBST, and incubated with 1:3500 anti-H3K9me3 antibody (#39161, Active Motif, Carlsbad, CA, USA). Blot was subsequently stripped, washed in TBST and incubated with 1:4000 anti-H3 antibody (ab1791, Abcam, Cambridge, United Kingdom).

Dnmt5 analysis and expression of DNA methyltransferase genes

The Dnmt5 amino-acid sequences of V. dahliae strain JR2 (VDAG_JR2_Chr1g14260) and C. neoformans var grubii strain H99 (CNAG_07552) were retrieved from EnsemblFungi, aligned using Uniprot.org/align and their domain structure was predicted using Interproscan. To assess expression of Dim2, Dnmt5, Rid, Hp1 and Dim5 in V. dahliae cultured in different growth media, wild-type V. dahliae was cultured in Czapec–Dox medium (CZA), half-strength Murashige–Skoog medium with vitamins and supplemented with 3% sucrose (MS), and potato dextrose broth (PDB) at 22 °C at 160 RPM in the dark for six days and mycelium was collected for three replicates per growth medium and ground as described above. To obtain RNA-seq data from V. dahliae grown in planta, 3-week-old A. thaliana (Col-0) plants were root dip inoculated in a conidiospore suspension of 106 spores per mL for 10 min. After root inoculation, plants were grown in individual pots in a greenhouse for 21 days, under a cycle of 16 h of light and 8 h of darkness, with temperatures maintained between 20 and 22 °C during the day and a minimum of 15 °C overnight. Three pooled samples (10 plants per sample) of complete flowering stems were used for total RNA extraction, respectively. Total RNA from in vitro cultured mycelium and was isolated using Trizol (Thermo Fisher Science, Waltham, MA, USA) following the manufacturer’s guidelines. Following RNA re-suspension, contaminating DNA was removed using the TURBO DNA-free kit (Ambion, Thermo Fisher Science, Waltham, MA, USA) and RNA integrity was determined by separating 2 μL of each sample on a 2% agarose gel and quantified using a Nanodrop (Thermo Fisher Science, Waltham, MA, USA) and stored at − 80 °C until further use. Library preparation was carried out at BGI (BGI, Hong Kong, China) and 50 bp fragments were sequenced using the BGISEQ-500 platform. Sequenced reads were mapped to V. dahliae strain JR2 gene annotation using Kallisto quant (settings: -single -l 50 -s 0.001 -pseudobam) to obtain normalized TPM values [46].

Transcriptional analysis of mutants in DNA methylation-associated genes

The V. dahliae wild-type strain JR2 [40], ∆Dim2, ∆Dnmt5, ∆Dim2-∆Dnmt5, ∆Rid, ∆Hp1 and ∆Dim5 were cultured in PDB at 22 °C at 160 RPM in the dark for 6 days. Mycelium collection, RNA-isolation and sequencing are performed for three replicates per growth condition as described above.

Differential gene and TE expression between mutants and wild-type was determined by mapping sequencing reads to the V. dahliae strain JR2 gene and TE annotation using TEtranscripts, which uses an iterative multimapping approach [23, 40, 47]. Genes and TEs were considered differentially expressed when they displayed log2FoldChange of < − 1 and a p value of < 0.05. Additionally, as transcript mapping to TEs using a multimapping approach may falsely identify transcription of highly similar sequences, we also used a mapping approach using unique reads only on TEs. The unique mapping approach resulted in slightly fewer differentially expressed TEs than the iterative multimapping approach (Additional file 1: Table S6).

Gene ontology (GO) terms were annotated to the V. dahliae JR2 proteome using Blast2GO (v1.4.4) [48]. GO enrichment analysis was performed using Ontologizer (v2.1), using Parent–Child-Union calculation and Benjamini–Hochberg multiple testing correction with 1000 resampling steps [49]. Fold change enrichment was calculated by dividing the fraction of genes annotated to each GO term in the study set by the fraction of genes annotated to each GO term in the population.

Clustering of induced genes and TEs was determined using CROC (settings, for genes: -w 50,000 -o 10,000 -m 5, for TEs: -w 100,000 -o 20,000 -m 5) [50]. To analyze whether induced genes or TEs display more clustering than in random gene or TE sets, 1000 random selections of 1280 genes and of 191 TEs were generated. These random gene and TE sets were similarly analyzed using CROC. Overlap with H3K9me3 domains and adaptive genomic regions [23] was assessed using bedtools closest (settings -d) [51].

Availability of data and materials

BSseq data of wild-type and mutants were submitted to the SRA database under the Accession numbers: PRJNA592220 and PRJNA659638. RNAseq data of in vitro cultivation were submitted to the SRA database under the Accession numbers: PRJNA592220 and PRJNA659638. RNAseq data of Arabidopsis infection were submitted to the SRA database under the Accession number: SRP149060. H3K9me3 ChIPseq data is submitted to the SRA database under the Accession number: PRJNA592220.

References

  1. Burdge GC, Hanson MA, Slater-Jefferies JL, Lillycrop KA. Epigenetic regulation of transcription: a mechanism for inducing variations in phenotype (fetal programming) by differences in nutrition during early life? Br J Nutr. 2007;97:1036–46.

    Article  Google Scholar 

  2. Chen ZJ. Genetic and epigenetic mechanisms for gene expression and phenotypic variation in plant polyploids. Annu Rev Plant Biol. 2007;58:377–406.

    Article  Google Scholar 

  3. Huang Q, Ma C, Chen L, Luo D, Chen R, Liang F. Mechanistic insights into the interaction between transcription factors and epigenetic modifications and the contribution to the development of obesity. Front Endocrinol. 2018. https://doi.org/10.3389/fendo.2018.00370.

    Article  Google Scholar 

  4. Luger K, Mäder AW, Richmond RK, Sargent DF, Richmond TJ. Crystal structure of the nucleosome core particle at 2.8 A resolution. Nature. 1997;389:251–60.

    Article  Google Scholar 

  5. Slotkin RK, Martienssen R. Transposable elements and the epigenetic regulation of the genome. Nat Rev Genet. 2007;8:272–85.

    Article  Google Scholar 

  6. Bird A. The essentials of DNA methylation. Cell. 1992;70:5–8.

    Article  Google Scholar 

  7. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11:204–20.

    Article  Google Scholar 

  8. Schmitz RJ, Lewis ZA, Goll MG. DNA methylation: shared and divergent features across eukaryotes. Trends Genet. 2019;35:818–27.

    Article  Google Scholar 

  9. Galagan JE, Henn MR, Ma LJ, Cuomo CA, Birren B. Genomics of the fungal kingdom: insights into eukaryotic biology. Genome Res. 2005;15:1620–31.

    Article  Google Scholar 

  10. Kouzminova E, Selker EU. Dim-2 encodes a DNA methyltransferase responsible for all known cytosine methylation in Neurospora. EMBO J. 2001;20:4309–23.

    Article  Google Scholar 

  11. Freitag M, Hickey PC, Khlafallah TK, Read ND, Selker EU. HP1 is essential for DNA methylation in Neurospora. Mol Cell. 2004;13:427–34.

    Article  Google Scholar 

  12. Tamaru H, Selker EU. A histone H3 methyltransferase controls DNA methylation in Neurospora crassa. Nature. 2001;414:277–83.

    Article  Google Scholar 

  13. Cambareri EB, Jensen BC, Schabtach E, Selker EU. Repeat-induced G-C to A-T mutations in Neurospora. Science. 1989;244:1571–5.

    Article  Google Scholar 

  14. Lewis ZA, Honda S, Khlafallah TK, Jeffress JK, Freitag M, Mohn F, Schübeler D, Selker EU. Relics of repeat-induced point mutation direct heterochromatin formation in Neurospora crassa. Genome Res. 2009;19:427–37.

    Article  Google Scholar 

  15. Ikeda KI, Van Vu B, Kadotani N, Tanaka M, Murata T, Shiina K, Chuma I, Tosa Y, Nakayashiki H. Is the fungus Magnaporthe losing DNA methylation? Genetics. 2013;195:845–55.

    Article  Google Scholar 

  16. Jeon J, Choi J, Lee G-W, Park S-Y, Huh A, Dean RA, Lee Y-H. Genome-wide profiling of DNA methylation provides insights into epigenetic regulation of fungal development in a plant pathogenic fungus Magnaportheoryzae. Sci Rep. 2015;5:8567.

    Article  Google Scholar 

  17. Huff JT, Zilberman D. Dnmt1-independent CG methylation contributes to nucleosome positioning in diverse eukaryotes. Cell. 2014;156:1286–97.

    Article  Google Scholar 

  18. Catania S, Dumesic PA, Pimentel H, Nasif A, Stoddard CI, Burke JE, Diedrich JK, Cook S, Shea T, Geinger E. Evolutionary persistence of DNA methylation for millions of years after ancient loss of a de novo methyltransferase. Cell. 2020;180:263–77.

    Article  Google Scholar 

  19. Bewick AJ, Hofmeister BT, Powers RA, Mondo SJ, Grigoriev IV, James TY, Stajich JE, Schmitz RJ. Diversity of cytosine methylation across the fungal tree of life. Nat Ecol Evol. 2019;3:479–90.

    Article  Google Scholar 

  20. Fradin EF, Thomma BPHJ. Physiology and molecular aspects of Verticillium wilt diseases caused by V. dahliae and V. albo-atrum. Mol Plant Pathol. 2006;7:71–86.

    Article  Google Scholar 

  21. Klosterman SJ, Atallah ZK, Vallad GE, Subbarao KV. Diversity, pathogenicity, and management of Verticillium species. Annu Rev Phytopathol. 2009;47:39–62.

    Article  Google Scholar 

  22. de Jonge R, Bolton MD, Kombrink A, van den Berg GCM, Yadeta KA, Thomma BPHJ. Extensive chromosomal reshuffling drives evolution of virulence in an asexual pathogen. Genome Res. 2013;23:1271–82.

    Article  Google Scholar 

  23. Cook DE, Kramer HM, Torres DE, Seidl MF, Thomma BPHJ. A unique chromatin profile defines adaptive genomic regions in a fungal plant pathogen. Elife. 2020;9:e62208.

    Article  Google Scholar 

  24. Seidl MF, Kramer HM, Cook DE, Lorencini Fiorin G, van den Berg GCM, Faino L, Thomma BPHJ. Repetitive elements contribute to the diversity and evolution of centromeres in the fungal genus Verticillium. MBio. 2020. https://doi.org/10.1128/mBio.01714-20.

    Article  Google Scholar 

  25. Depotter JRL, Shi-Kunne X, Missonnier H, Liu T, Faino L, Berg GCM, Wood TA, Zhang B, Jacques A, Seidl MF, et al. Dynamic virulence-related regions of the plant pathogenic fungus Verticillium dahliae display enhanced sequence conservation. Mol Ecol. 2019. https://doi.org/10.1111/mec.15168.

    Article  Google Scholar 

  26. Faino L, Seidl MF, Shi-Kunne X, Pauper M, Van Den Berg GCM, Wittenberg AHJ, Thomma BPHJ. Transposons passively and actively contribute to evolution of the two-speed genome of a fungal pathogen. Genome Res. 2016;26:1091–100.

    Article  Google Scholar 

  27. Klosterman SJ, Subbarao KV, Kang S, Veronese P, Gold SE, Thomma BPHJ, Chen Z, Henrissat B, Lee YH, Park J, et al. Comparative genomics yields insights into niche adaptation of plant vascular wilt pathogens. PLoS Pathog. 2011;7:e1002137.

    Article  Google Scholar 

  28. de Jonge R, Peter van Esse H, Maruthachalam K, Bolton MD, Santhanam P, Saber MK, Zhang Z, Usami T, Lievens B, Subbarao KV, et al. Tomato immune receptor Ve1 recognizes effector of multiple fungal pathogens uncovered by genome and RNA sequencing. Proc Natl Acad Sci USA. 2012;109:5110–5.

    Article  Google Scholar 

  29. Möller M, Habig M, Lorrain C, Feurtey A, Haueisen J, Fagundes WC, Alizadeh A, Freitag M, Stukenbrock EH. Recent loss of the Dim2 DNA methyltransferase decreases mutation rate in repeats and changes evolutionary trajectory in a fungal pathogen. PLoS Genet. 2021;17:e1009448.

    Article  Google Scholar 

  30. Dowen RH, Pelizzola M, Schmitz RJ, Lister R, Dowen JM, Nery JR, Dixon JE, Ecker JR. Widespread dynamic DNA methylation in response to biotic stress. Proc Natl Acad Sci USA. 2012;109:e2183–91.

    Article  Google Scholar 

  31. Meaney MJ, Szyf M. Environmental programming of stress responses through DNA methylation: life at the interface between a dynamic environment and a fixed genome. Dialogues Clin Neurosci. 2005;7:103–23.

    Article  Google Scholar 

  32. Wang Y, Wang Z, Liu C, Wang S, Huang B. Genome-wide analysis of DNA methylation in the sexual stage of the insect pathogenic fungus Cordyceps militaris. Fungal Biol. 2015;119:1246–54.

    Article  Google Scholar 

  33. Zhang X, Liu X, Zhao Y, Cheng J, Xie J, Fu Y, Jiang D, Chen T. Histone H3 lysine 9 methyltransferase DIM5 is required for the development and virulence of Botrytis cinerea. Front Microbiol. 2016. https://doi.org/10.3389/fmicb.2016.01289.

    Article  Google Scholar 

  34. Haueisen J, Möller M, Eschenbrenner CJ, Grandaubert J, Seybold H, Adamiak H, Stukenbrock EH. Highly flexible infection programs in a specialized wheat pathogen. Ecol Evol. 2019;9:275–94.

    Article  Google Scholar 

  35. Gessaman JD, Selker EU. Induction of H3K9me3 and DNA methylation by tethered heterochromatin factors in Neurospora crassa. Proc Natl Acad Sci. 2017;114:E9598–607.

    Article  Google Scholar 

  36. Hayashi A, Ishida M, Kawaguchi R, Urano T, Murakami Y, Nakayama J. Heterochromatin protein 1 homologue Swi6 acts in concert with Ers1 to regulate RNAi-directed heterochromatin assembly. Proc Natl Acad Sci USA. 2012;109:6159–64.

    Article  Google Scholar 

  37. Basenko EY, Sasaki T, Ji L, Prybol CJ, Burckhardt RM, Schmitz RJ, Lewis ZA. Genome-wide redistribution of H3K27me3 is linked to genotoxic stress and defective growth. Proc Natl Acad Sci. 2015;112:E6339–48.

    Article  Google Scholar 

  38. Möller M, Schotanus K, Soyer JL, Haueisen J, Happ K, Stralucke M, Happel P, Smith KM, Connolly LR, Freitag M, et al. Destabilization of chromosome structure by histone H3 lysine 27 methylation. PLoS Genet. 2019;15:e1008093.

    Article  Google Scholar 

  39. Holliday R, Grigg GW. DNA methylation and mutation. Mutat Res Fundam Mol Mech Mutagen. 1993;285:61–7.

    Article  Google Scholar 

  40. Faino L, Seidl M, Datema E, van den Berg GCM, Janssen A, Wittenberg AHJ, Thomma BPHJ. Single-molecule real-time sequencing combined with optical mapping yields completely finished fungal genome. MBio. 2015;6:e00936-15.

    Article  Google Scholar 

  41. Santhanam P. 2012. Random insertional mutagenesis in fungal genomes to identify virulence factors. In: Plant fungal pathogens. Springer, pp. 509–57.

  42. Frandsen RJN, Andersson JA, Kristensen MB, Giese H. Efficient four fragment cloning for the construction of vectors for targeted gene replacement in filamentous fungi. BMC Mol Biol. 2008;9:70.

    Article  Google Scholar 

  43. Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinform. 2009;10:232.

    Article  Google Scholar 

  44. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, Mason CE. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13:R87.

    Article  Google Scholar 

  45. Gel B, Serra E. karyoploteR: an R/Bioconductor package to plot customizable genomes displaying arbitrary data. Bioinformatics. 2017;33:3088–90.

    Article  Google Scholar 

  46. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7.

    Article  Google Scholar 

  47. Jin Y, Tam OH, Paniagua E, Hammell M. TEtranscripts: a package for including transposable elements in differential expression analysis of RNA-seq datasets. Bioinformatics. 2015;31:3593–9.

    Article  Google Scholar 

  48. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

    Article  Google Scholar 

  49. Bauer S, Grossmann S, Vingron M, Robinson PN. Ontologizer 2.0—a multifunctional tool for GO term enrichment analysis and data exploration. Bioinformatics. 2008;24:1650–1.

    Article  Google Scholar 

  50. Pignatelli M, Serras F, Moya A, Guigó R, Corominas M. CROC: finding chromosomal clusters in eukaryotic genomes. Bioinformatics. 2009;25:1552–3.

    Article  Google Scholar 

  51. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

    Article  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by a Ph.D. grant of the Research Council Earth and Life Sciences (Project 831.15.002) to HMK and by and EMBO and Human Frontier Science Program Postdoctoral Fellowship (HFSP, LT000627/2014-L) to DEC. Work in the laboratories of M.F.S and B.P.H.J.T. is supported by the Research Council Earth and Life Sciences (ALW) of the Netherlands Organization of Scientific Research (NWO) and B.P.H.J.T acknowledges funding by the Alexander von Humboldt Foundation in the framework of an Alexander von Humboldt Professorship endowed by the German Federal Ministry of Education and Research and is furthermore supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—EXC 2048/1—Project ID: 390686111.

Author information

Author notes

  1. H. Martin Kramer and David E. Cook contributed equally to this work

Authors and Affiliations

  1. Laboratory of Phytopathology, Wageningen University and Research, Droevendaalsesteeg 1, 6708 PB, Wageningen, The Netherlands

    H. Martin Kramer, David E. Cook, Grardy C. M. van den Berg, Michael F. Seidl & Bart P. H. J. Thomma

  2. Department of Plant Pathology, Kansas State University, 1712 Claflin Road, Manhattan, KS, 66506, USA

    David E. Cook

  3. Theoretical Biology & Bioinformatics, Department of Biology, Utrecht University, Utrecht, The Netherlands

    Michael F. Seidl

  4. Institute for Plant Sciences, Cluster of Excellence on Plant Sciences (CEPLAS), University of Cologne, 50674, Cologne, Germany

    Bart P. H. J. Thomma

Authors
  1. H. Martin Kramer
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. David E. Cook
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Grardy C. M. van den Berg
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Michael F. Seidl
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Bart P. H. J. Thomma
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

DEC and BPHJT conceived the study. HMK, DEC, MFS and BPHJT designed experiments. HMK, DEC and GCMB carried out laboratory experiments and HMK, DEC and MFS analyzed data. HMK, DEC, MFS and BPHJT wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Bart P. H. J. Thomma.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Table S1. Differentially expressed genes in the mutants relative to wild-type. Table S2. Primers used to generate homologous recombination knockout and complementation vectors in V. dahliae. Table S3. Primers used to confirm homologous recombination knockout and complementation vectors in V. dahliae. Table S4. Primers used to confirm gene deletion in V. dahliae. Table S5. Primers used for qPCR fungal biomass quantification of V. dahliae infecting tomato. Table S6. TE-mapping parameter.

Additional file 2:

Figure S1. Phylogenetic tree of DNA methyltransferases. Figure S2. Verticillium dahliae ∆Dim5 loses H3K9me3. Figure S3. Growth assay of complementation strains. Figure S4. Stress assay pictures at 10 dpi. Figure S5. DNA methylation in CG context. Figure S6. DNA methylation in CHG context. Figure S7. DNA methylation in CHH context. Figure S8. DNA methylation over the genome. Figure S9. Occurrence of CG, CHG and CHH sites in methylated and non-methylated transposable elements. Figure S10. Comparison of protein domain structure of C. neoformans and V. dahliae Dnmt5. Figure S11. Distribution of H3K9me3 domain lengths. Figure S12. Genes induced in Hp1 and Dim5 mutants cluster more often than expected based on chance. Figure S13. Transposons induced in Hp1 and Dim5 mutants cluster more often than expected based on chance. Figure S14. Clusters of genes and transposons over all chromosomes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kramer, H.M., Cook, D.E., van den Berg, G.C.M. et al. Three putative DNA methyltransferases of Verticillium dahliae differentially contribute to DNA methylation that is dispensable for growth, development and virulence. Epigenetics & Chromatin 14, 21 (2021). https://doi.org/10.1186/s13072-021-00396-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13072-021-00396-6

Keywords

Epigenetics & Chromatin

ISSN: 1756-8935