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

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. Supplementary Information The online version contains supplementary material available at 10.1186/s13072-021-00396-6.

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 hemimethylated 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 (BSseq), 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].

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 methyltransferasecomplex 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).
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 wildtype 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].
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 wildtype 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 wildtype 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 (See figure on next page.) Fig. 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 mm 2 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 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  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. 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 chromoshadow 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 (See figure on next page.) Fig. 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 Kramer et al. Epigenetics & Chromatin (2021) 14:21 repressed in ∆Hp1, and 1617 genes were induced and 781 were repressed in ∆Dim5 when compared with wildtype (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. 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 H3K9me3deficient 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 decondensation. 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 cellcycle 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.

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), 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 10 6 conidiospores per mL. Subsequently, 10 µL of conidiospore suspension, containing 10 4 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 × 10 3 conidiospores on PDA without supplement, or on PDA supplemented with 1 M NaCl, 1 M Sorbitol, 2 mM H 2 O 2 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 10 6 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: ). TEs were categorized as methylated if any cytosine in the sequence showed evidence for methylation in wild-type JR2.

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 10 6 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 resuspension, 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].
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.