Phylogenetic and functional analyses of N⁶-methyladenosine RNA methylation factors in the wheat scab fungus Fusarium graminearum

To resolve phylogenetic relationships of genes encoding an MT-A70 protein domain (Pfam domain: PF05063), we collected 1,568 protein sequences from the UniProt database (accessed on 21 May 2023; https://www.uniprot.org/) and curated the list of MT-A70-containing proteins, excluding redundant entries or sequences with an E-value of smaller than 10⁻⁵ for PF05063 (see Section Materials and Methods). Following the manual curation, a maximum likelihood tree was reconstructed for 1,116 protein sequences containing an MT-A70 domain in 829 diverse fungal species. This tree revealed the presence of three distinct clades, each including the three previously characterized m⁶A writers in fungi: IME4, KAR4, and MTA1 (Fig. 1A). Within the IME4 and KAR4 clades, there was a branch exhibiting an early divergence, which exclusively consists of sequences derived from early-diverging fungi, such as species belonging to the phyla Chytridiomycota, Mucoromycota, and Zoopagomycota (see clades labeled with red asterisks in Fig. 1A). Notably, distinct patterns were observed in the phyletic distribution of IME4, KAR4, and MTA1 among different taxonomic groups in fungi (Fig. 1B). Some early-diverging fungi and species within the Pucciniomycotina subphylum were found to possess all the three MT-A70-containing proteins, IME4, KAR4, and MTA1, although many are lacking one of these proteins. The Basidiomycota species (mushroom-forming fungi) and Saccharomycotina species (budding yeasts) typically had IME4 and KAR4, whereas nearly all of the Pezizomycotina species that are comprised of entirely filamentous fungi were found to possess MTA1 as a sole MT-A70-containing protein.

Independently reconstructed maximum likelihood trees were generated for the IME4, KAR4, and MTA1 clades including their respective homologs of METTL3, METTL14, and METTL4, which were found in humans (Homo sapiens), mouse (Mus musculus), zebrafish (Danio rerio), fruit fly (Drosophila melanogaster), and a model plant (Arabidopsis thaliana) serving as an outgroup (Fig. 1C through E; see Fig. S1 for bootstrap values and protein IDs). The phylogenetic patterns of the IME4 and KAR4 clades exhibited resemblance, encompassing early-diverging fungi and species from Basidiomycota and Saccharomycotina. Among species possessing either IME4 (294 species) or KAR4 (286 species), three-quarters of the species (221 species) had both m⁶A writers (Table S1). As previously mentioned, some early-diverging fungi exhibited the presence of two copies of IME4 and KAR4. One of these copies appears to have undergone divergence before the emergence of Dikarya (commonly known as “higher fungi”), forming distinct clades here dubbed IME4-like and KAR4-like with robust 100% bootstrap support (Fig. S1). The phylogeny of MTA1 exhibited monophyly and was largely congruent to the phylogeny of Pezizomycotina (Fig. 1E).

Species belonging to Pezizomycotina possessed MTA1 as a sole MT-A70-containing protein. On the other hand, many early-diverging fungi and species from Pucciniomycotina exhibited the presence of three m⁶A writers, IME4, KAR4, and MTA1 (Fig. 2A). We took a closer look at 26 species that featured the three m⁶A writers. The increased abundance of MT-A70 domain-containing proteins in early-diverging fungi can be attributable to the existence of IME4-like and KAR4-like genes that diverged from MTA1 at earlier points in time. Although the distribution patterns of MT-A70 domain-containing proteins in early-diverging fungi that belong to the class Chytridiomycetes and Zoopagomycetes are highly variable, they all contained IME4 (or IME4-like), KAR4 (or KAR4-like), and MTA1 (Fig. 2B). Species that belong to the Mucoromycota (except for species from the Endogonales) all lacked IME4 and KAR4, instead having IME4-like and KAR4-like, suggesting that some lineages in the Mucorales have lost IME4 and KAR4 that can be found in species in the Glomeromycetes (the fungal Class sister to the Mucoromycetes). Interestingly, some species in the Mucorales and Mortierellales have two copies of KAR4-like, suggesting gene duplication events. Species that belong to the Pucciniomycotina, the earliest-diverging subphylum within Basidiomycota, have IME4, KAR4, and MTA1 but are devoid of IME4-like and KAR4-like genes (Fig. 2B), suggesting that these early-diverged MT-A70 domain-containing genes were completely lost in Dikarya.

Recent studies revealed the components of yeast m⁶A methyltransferase complexes were highly conserved with those found in mammals, insects, and plants, including IME4, KAR4, MUM2 (a homolog of WTAP), and VIR1 (32, 33). Thus, we conducted a search in the UniProt database for relevant Pfam domains, specifically PF17098 for the WTAP/Mum2p family and PF15912 for the N-terminal domain of virilizer (hereafter coined VIRN). We identified 172 and 53 fungal proteins that harbor PF17098 for WTAP and PF15912 for VIRN, respectively, in the current UniProt database (accessed on 4 June 2023; Table S1). After manual curation, we reconstructed phylogenetic trees for WTAP and VIRN, together with homologs found in humans, mice, zebrafish, fruit flies, and Arabidopsis thaliana as an outgroup. WTAP orthologs including Mum2 in yeasts were found only in early-diverging fungi and species from Basidiomycota, and Saccharomycotina, lacking in Pezizomycotina species (Fig. S1). Only 37 fungal species that belong to early-diverging fungi and Basidiomycota were shown to have genes encoding a VIRN domain (Fig. S1). The paucity of homologous sequences in fungi may be attributable to the fact that WTAP and VIRN sequences have diverged significantly from the known sequences, making their identification challenging, especially in Ascomycetous fungi, such as Saccharomycotina and Pezizomycotina. Indeed, a homolog of the virilizer protein, VIR1, was recently found in the budding yeast, which lacks the apparent Pfam domain for VIRN but it exhibits a 3D structure similar to the human homolog VIRMA (32).

Among RNA-binding proteins specifically recognizing and binding to m⁶A, the best studied readers are YTH domain-containing proteins. PHO92 was found as the only protein containing the YTH domain in the budding yeast, while animals can have up to five such proteins and plants can have more than 10. A recent study revealed target transcripts of PHO92 and its important roles during meiosis in budding yeast (37). However, little is known about m⁶A readers in other fungal taxa. Therefore, we searched for proteins containing a YTH domain (PF04146) in the UniProt database (accessed on 21 May 2023) and found 2,029 fungal sequences. To shorten the list of potential m⁶A readers, entries with irrelevant protein names (e.g., DNA repair protein RAD51) were excluded, and protein sequences below 300 or above 1,000 amino acids in length were filtered out. After manual curation, we reconstructed a phylogenetic tree for the final list of 868 YTH domain-containing proteins in 695 diverse fungal species (Fig. 3A; see Fig. S1 for bootstrap values and protein IDs). The tree was rooted to YTH domain-containing proteins outside of fungi, including human YTHDF2. The tree was composed of three main clades that were strongly supported by high bootstrap values. One clade contained PHO92 homologs from species belonging to Saccharomycotina and Basidiomycota, as well as YTH domain-containing proteins found in mammals, an insect, and a plant. Another clade exclusively comprised species from Pezizomycotina and was placed sister to the PHO92 clade. The third clade, here dubbed YTH1, consisted of early-diverging fungi and species from Basidiomycota and Ascomycota, displaying a topology that largely matched the species phylogeny. Approximately half of the Agaricomycotina species (67 out of 128) and one-fifth of the Pezizomycotina species (98 out of 498) examined in this study were found to possess both PHO92 and YTH1 while early-diverging fungi and species belonging to Pucciniomycotina, Ustilaginomycotina, and Wallemiomycotina, were found to have only PHO92 (Fig. 3B and C). While the majority of species belonging to Saccharomycotina (30 out of 37) were found to possess PHO92, there were several species including Yarrowia lipolytica (a lipid-producing yeast) that lacked PHO92 and instead exhibited the presence of YTH1 (Fig. 3C; Table S1).

2OG-Fe[II] dioxygenase AlkB-like domain-containing proteins were known to play roles as m⁶A erasers by removing the methyl group from m⁶A in mammals and plants (26, 48–50). However, there has been no information on AlkB homologs (ALKBHs) in fungi. Thus, we searched for protein sequences containing a 2OG-Fe(II) dioxygenase superfamily domain in the genome of F. graminearum strain PH-1. We identified five genes harboring a 2OG-Fe(II) dioxygenase domain (PF13532) with an E-value of smaller than 10⁻¹⁰ (gene ID: FGRRES_01255, FGRRES_09872, FGRRES_16456, FGRRES_16652, and FGRRES_20373) (Table 1). 2OG-Fe(II) dioxygenases are involved in diverse cellular processes, such as demethylation of DNA/RNA and modification of histone proteins, as well as secondary metabolite production. To get insights into evolutionary relationships of different types of ALKBH, we reconstructed a phylogenetic tree of ALKBH proteins found in 27 selected fungi, including several functionally characterized ALKBH proteins in humans and Arabidopsis (Fig. 4). Early-diverging fungi, such as Basidiobolus meristosporus and Spizellomyces punctatus, possessed 10 and 9 2OG-Fe(II) dioxygenases, whereas yeasts that belong to Saccharomycotina tend to have only one or two (Table S2). ALKBH2 (FGRRES_09872), ALKBH3 (FGRRES_16456), and ALKBH4 (FGRRES_01255 were closely related to human ALKBH3 and Arabidopsis ALKBH2, which are responsible for repairing DNA damage caused by alkylation (51, 52). ALKBH5 (FGRRES_20373) was placed sister to the clade containing human ALKBH6, which plays a role in DNA repair (53, 54). ALKBH1 (FGRRES_16652) formed a well-supported clade with human and Arabidopsis ALKBH1s, known for their versatile functions in removing methyl groups from DNA/RNA, as well as modifying histone H2A (55–57). It was shown that human ALKBH5 and Arabidopsis ALKBH9B and ALKBH10B demethylate m⁶A (48, 58–60). These m⁶A erasers in mammals and plants formed a clade ancestral to clades including ALKBH1 and ALKBH5 in F. graminearum (Fig. 4).

To examine the functions of m⁶A factors in F. graminearum during sexual development, we conducted a search for Pfam domains associated with m⁶A writer, reader, and eraser, and discovered several potential m⁶A factors (Table 1). Among them, F. graminearum possessed genes homologous to MTA1 and YTH1, along with the previously mentioned five ALKBHs. In addition, we identified two genes that encoded hypothetical proteins harboring PF17098 for WTAP and PF15912 for VIRN, respectively, with E-values of greater than 10⁻⁵ (Table 1). Despite the low sequence similarity to WTAP and VIRN found in fungi, we included these potential m⁶A writers in the knockout study because their expression levels were observed to increase during sexual development (see Fig. 7 in the section “Roles of MTA1 in conidial germination and hyphal growth”). We individually deleted m⁶A factors in F. graminearum strain PH-1 [wildtype (WT)]. The authenticity of the resulting knockout mutants was verified by diagnostic PCRs (Fig. 5A). The knockouts displayed normal growth and produced perithecia, which are characterized by their dark purple coloration and flask-shaped or spherical structures (Fig. 5B). Morphologies of ascospores (the products of meiosis) enveloped in sac-like structures called asci were also normal in the knockouts, indicating that m⁶A factors did not affect sexual development in F. graminearum (Fig. S2).

Since m⁶A writer, IME4, in the budding yeast, plays a crucial role in initiating sexual development, we carefully examined the genotype of Δmta1. Although genotyping confirmed the deletion of MTA1 in Δmta1, we were intrigued to find that a PCR band specific to MTA1 could still be amplified from Δmta1. This suggested that Δmta1 is heterokaryotic, indicating that it still contains one or more “wild-type” nuclei. To obtain homokaryons, we performed additional single-spore isolations, using Δmta1, and examined 10 single-spored isolates derived from Δmta1. Diagnostic PCR analysis confirmed the deletion of MTA1 in all 10 isolates; yet, a fragment of MTA1 was still amplified from these isolates (Fig. S3). The expression levels of MTA1 were markedly reduced in Δmta1, suggesting a lower abundance of wild-type nuclei compared to the WT (Fig. 5C). This observation was further supported by RNA-seq analysis. When examining the mapped reads on the MTA1 locus in both WT and Δmta1, we found a decrease in the number of reads mapped to the MTA1 locus in Δmta1 compared to the WT (Fig. S3). The presence of RNA-Seq reads mapping to the MTA1 locus, which should have been absent in homokaryotic Δmta1, served as additional proof for the heterokaryotic nature of Δmta1 and indicated that Δmta1 is, in fact, a knockdown strain. Consequently, we will refer to Δmta1 as the MTA1-KD1 strain hereafter.

Extensive efforts were unsuccessful in obtaining knockouts completely lacking MTA1 (Fig. S3), presumably due to its indispensable role in F. graminearum. To get better insight into the roles of m⁶A RNA methylation in F. graminearum, we generated strains overexpressing MTA1. RT-PCR analysis confirmed that two strains, MTA1-OE4 and MTA1-OE5, exhibited overexpression of MTA1 (Fig. 6A). It was notable that MTA1 was overexpressed to a greater extent in the MTA1-OE5 strain, compared to MTA1-OE4. Both semi-quantitative RT-PCR and RNA-seq analyses indicated that the expression level of the MTA1-OE5 strain was approximately 100 times greater than that of the WT (Fig. S4). The WT and MTA1-OE4 strains produced normal asci containing mature ascospores 6 days after the induction of sexual development, while it was not until 11 days after the induction that the MTA1-OE5 strain produced fully developed asci (Fig. 6B). The difference in the observed phenotypes in the two MTA1-overexpressing strains could likely be attributable to the variation in the degrees of MTA1 expression, and thus the MTA1-OE5 strain which displayed substantial overexpression of MTA1 was used for further analysis and investigation. To determine whether MTA1 is required for m⁶A RNA methylation in F. graminearum, we compared total m⁶A RNA methylation levels between the WT and MTA1-OE5 strains. The amount of m⁶A RNA in WT was 0.075% ± 0.015% of the total RNA, which was approximately half of that in the MTA1-OE5 strain (0.152% ± 0.033% of the total RNA), which indicated that the m⁶A RNA methylation level was significantly increased in the MTA1-OE5 strain compared to the WT (Fig. 6C).

In previous studies, we obtained transcriptome data for conidial germination on the Bird agar medium (61), as well as for perithecial development on the Carrot agar medium (62). By integrating these data sets with the RNA-seq data obtained in the present study (freshly harvested conidia), we conducted a comprehensive analysis of expression level changes in m⁶A factors throughout the life cycle of F. graminearum, encompassing various stages of vegetative growth and perithecial development (Fig. 7). The m⁶A writer MTA1 and a potential m⁶A reader YTH1 demonstrated a synchronized expression pattern, as indicated by a high Pearson correlation coefficient of 0.94. Their expression levels exhibited a significant increase during stage 1 (S1, 15 minutes after incubation of conidia on Bird agar medium), followed by a decline during hyphal growth (S2 and S3), and remained at a basal level during perithecial development (S4–S9). Interestingly, the expression levels of ALKBH1 and ALKBH5, the two most probable m⁶A erasers in F. graminearum, reached their peak at stage 2 (S2, 3 h after incubation), at which conidia germinated and hyphae started to extend (Fig. S5). During active hyphal growth at stage 3 (S3, 11 h after incubation), the major components of the fungal cytoskeleton, actin (FGRRES_07735) α-tubulin (FGRRES_00639), and β-tubulin (FGRRES_09530) that play a crucial role in polarity establishment, maintenance, and polar growth, displayed the greatest induction during the life cycle of F. graminearum.

Taking into account the expression profile, we proceeded to investigate the phenotypic characteristics of the MTA1-OE5 strain during conidial germination and hyphal growth. In contrast to the WT, the hyphal tips did not exhibit elongation in an overnight culture of the MTA1-OE5 strain cultivated in quarter-strength potato dextrose broth (q-PDB) medium (Fig. 7B). Furthermore, in the MTA1-OE5 strain growing on potato dextrose agar (PDA) medium, the hyphae at the leading edge displayed a tendency to have fewer branches compared to the WT (Fig. S6). Next, we conducted measurements of the total biomass of the WT and MTA1-OE5 strains cultured in q-PDB medium for a duration of 8 days. The MTA1-OE5 strain exhibited a slower growth rate compared to the WT, particularly during the initial 2 days of cultivation (Fig. 7C). To investigate genes responsible for the delayed conidial germination observed in the MTA1-OE5 strain, we performed differential gene expression analyses between the WT, MTA1-KD1, and MTA1-OE5 strains, using RNA-seq data obtained from freshly harvested conidia (stage 0, S0). In this analysis, we identified a total of 20 differentially expressed (DE) genes [fold change > 4 at a false discovery rate (FDR) of 5%], when comparing the WT and MTA1-OE5 strains, as well as 21 DE genes when comparing the MTA1-KD1 and MTA1-OE5 strains (Fig. 7D). Notably, 16 genes were found to be commonly differentially expressed in both comparisons, including MTA1. There were no DE genes observed between the WT and MTA1-KD1 strains (Fig. S7). Among the DE genes, we found that 17 genes were transcriptionally upregulated, and eight genes were downregulated in the MTA1-OE5 strain (Table 2). The DE genes included nine genes encoding hypothetical proteins without any predicted protein domain and 16 functionally annotated genes encoding two transcription factors, two transporters, and diverse enzymes. FGRRES_00725 and FGRRES_05926 encoding GAL4-type transcription factors were significantly downregulated in the MTA1-OE5 strain, while FGRRES_02139 and FGRRES_04188 encoding an ABC multidrug transporter and major facilitator superfamily, respectively, were highly upregulated in the MTA1-OE5 strain (Table 2). In consistent with the expression pattern of MTA1, the expression levels of two transporters were elevated at S1 and reduced thereafter. For two transcription factors, the expression levels of FGRRES_05926 and FGRRES_00725 began to decline at S2 and S1, respectively, and remained basal, suggesting that their roles in early stages of conidial germination and hyphal growth.

In the phylum Ascomycota, the roles of m⁶A writers in budding yeasts (Saccharomycotina) and filamentous fungi (Pezizomycotina) may have evolved independently. Species belonging to the Pezizomycotina possess only a single m⁶A writer MTA1, a homolog of METTL4, while species belonging to the Saccharomycotina have two m⁶A writers, IME4 and KAR4 that are homologs of METTL3 and METTL14, respectively. In budding yeast, IME4 and KAR4 are crucial for the initiation of meiosis for ascospore formation (29, 63). However, the m⁶A writer MTA1 is dispensable for perithecial development and ascospore production in filamentous fungi, F. graminearum and M. oryzae (46). The relative abundance of m⁶A in M. oryzae (0.069% ± 0.003%) is comparable to that in F. graminearum (0.075 ± 0.015%) (47), suggesting that the degree of m⁶A levels of total mRNA is similar in these two filamentous fungi. In the rice blast fungus M. oryzae, a knockout strain lacking MTA1 showed a defect in appressorium formation during the infection process on the epidermal cells of rice (47). Expression levels of autophagy-related genes in M. oryzae were changed according to the degree of m⁶A level, which is important for normal appressorium formation. Although the homozygous Δmta1 strain in M. oryzae was viable, our extensive efforts to obtain homozygous Δmta1 were unsuccessful in F. graminearum, suggesting its essential role. Given the significant impact on conidial germination and hyphal growth observed in the MTA1-overexpressing strain in F. graminearum, the role of MTA1 appears to have diverged in the two plant pathogenic fungi.

Most importantly, F. graminearum does not form appressoria and utilizes distinct strategies to penetrate the epidermal cells of the host plants. Both F. graminearum and M. oryzae belong to the Sordariomycetes but fall into different orders, Hypocreales and Magnaporthales, respectively (64). RNA-seq analyses indicated significant divergence in gene expression including autophagy-related genes during conidial germination and infection processes between F. graminearum and M. oryzae (61). The m⁶A writers IME4 and KAR4 were crucial for the initiation of the sexual cycle and the formation of ascospores in the budding yeast. However, the m⁶A writer MTA1 does not appear to be responsible for perithecial development and ascospore production in F. graminearum and M. oryzae. Although the MTA1-OE5 strain exhibited delayed perithecial development, it produced normal ascospores. In the Pezizomycotina, A-to-I RNA editing, one of the eukaryotic RNA modifications, is crucial for sexual development and ascospore formation (43–45, 65–67).

Several ancestral traits observed in early-diverging fungi are shared with metazoans or unicellular opisthokonts, which have been subjected to extensive parallel loss across the Dikarya lineages (68, 69). As with metazoans, many species of early-diverging fungi and species that belong to the Pucciniomycotina (the earliest-diverging subphylum within Basidiomycota) possess all three m⁶A writers, IME4, MTA1, and KAR4. These suggested that budding yeasts (Saccharomycotina) may have lost the m⁶A writer MTA1. Conversely, species belonging to Pezizomycotina, which entirely consists of filamentous fungi, may have lost the m⁶A writers IME4 and KAR4. The phyletic distribution of m⁶A writers in fungi suggested that Pezizomycotina species including F. graminearum and M. oryzae have likely embraced m⁶A machineries, such as MTA1, as a means to regulate filamentous growth and facilitate host penetration. This adoption of m⁶A machinery may have played a crucial role in the development and maintenance of the filamentous morphology that is characteristic of these fungi, allowing them to effectively colonize and interact with their respective hosts.

In contrast to well-characterized yeast m⁶A methyltransferase complexes including IME4 and KAR4 m⁶A writers, our understanding of m⁶A machinery in filamentous fungi is still limited. Despite our functional characterization of potential WTAP and VIR1 homologs in F. graminearum, it is unlikely that these represent genuine homologs. This is because the homologs of WTAP/MUM2 and VIR1 are missing in species within all the Pezizomycotina species. The failure to identify homologs of WTAP and VIR1 in our study is likely due to significant sequence divergence of these m⁶A factors within the Pezizomycotina. The sequence differences might have hindered their detection using conventional homology search methods, indicating that these m⁶A factors in Pezizomycotina may have undergone substantial changes in the primary amino acid sequences. Until recently, the existence and function of VIR1 in the budding yeast had remained unknown, but through the application of protein folding prediction tools, its identity was revealed (32). Recent advancements in protein 3D structure prediction tools to uncover hidden homologs (70, 71) hold promise for identifying homologous members of the m⁶A methyltransferase complexes in Pezizomycotina species. Alternatively, it is possible that components of m⁶A machinery associated with MTA1 are completely distinct from those found in the yeast m⁶A methyltransferase complex associated with IME4 and KAR4. In this case, it will be necessary to identify components that interact with MTA1 via co-immunoprecipitation.

In addition, it is worth mentioning, considering the number of putative m⁶A writers and m⁶A erasers found in early-diverging fungi and species from the Pucciniomycotina (up to five m⁶A writers in a species), not many m⁶A readers containing an YTH domain were identified, suggesting that there may be different types of RNA-binding proteins that recognize m⁶A modification. It would be interesting to study the possibly divergent roles of m⁶A factors in these relatively understudied fungal taxa.

The levels of MTA1 and YTH1 expression were significantly elevated at 15 minutes after placing conidia on the Bird medium (S1) and then decreased sharply to a basal level during polar growth and hyphal branching (S2–S3). These findings suggest that the potential m⁶A writer and m⁶A reader likely have a role during the initial stage of conidial germination. In line with this hypothesis, we observed slower conidial germination and less hyphal branching in the MTA1-OE5 strain. Although m⁶A demethylase activities of ALKBH genes have never been confirmed in fungi, the three putative m⁶A erasers, ALKBH1, ALKBH4, and ALKBH5, were induced during polar growth (S2) or hyphal branching (S3). These possible m⁶A erasers might be involved in maintaining low m⁶A levels during active vegetative growth. The significant increase in expression of actin and tubulin genes at S2–S3 reflected extensive hyphal growth and branching. Specifically, the expression level of ALKBH4 reached its peak at S0 (fresh conidia) but experienced a sharp decline at S1 (the initial stage of conidial germination). However, it was gradually induced during polar growth and hyphal branching stages (S2-S3), exhibiting an opposite expression pattern compared to MTA1 and YTH1. Although we did not observe any phenotypic changes in the knockout mutant lacking ALKBH1, it would be possible that other ALKBH genes may be involved in the regulation of conidial germination and hyphal growth by maintaining low m⁶A level in F. graminearum. Among the five putative m⁶A erasers, ALKBH2 showed a dramatic increase in expression during the ascospore formation. These expression dynamics of m⁶A factors highlight the complex regulation and involvement of m⁶A factors throughout different stages of the F. graminearum life cycle.

The largest subphylum Pezizomycotina in fungi underscores their importance and their complex interactions with humans, encompassing a wide range of fungi that have significant influences on humans, both negative and positive. Examples include the opportunistic human pathogen Aspergillus fumigatus and the dermatophyte Trichophyton rubrum, as well as many plant pathogenic fungi causing serious diseases on our staple crops, such as F. graminearum and M. oryzae. The Pezizomycotina also includes several ecologically significant species, encompassing those that play crucial roles in wood and litter decay processes, as well as those that form symbiotic associations with other organisms, including lichens. In filamentous fungi, the m⁶A writer MTA1 appears to have evolved to have important roles in adapting to diverse ecological niches, particularly in relation to filamentous growth. To identify target genes of MTA1 that caused delayed conidial germination and slower hyphal growth in the MTA1-OE5 strain, we are currently investigating potential m⁶A sites in transcripts, using an Oxford Nanopore direct RNA sequencing technology.

We downloaded protein sequences from the UniProt database (https://www.uniprot.org/) that possessed pfam domains associated with m⁶A factors. These domains include PF05063 (MT-A70), PF04146 (YTH), PF13532 (ALKB), PF17098 (WTAP/MUM2), and PF15912 (Virilizer, N-terminal). However, the initial lists of potential m⁶A factors contained duplications and likely misannotated sequences. To address this issue, we conducted a manual inspection of the lists, excluding sequences that were either too short or too long, as well as sequences with small E -values, applying specific thresholds (e.g., E -value > 10^–5). Whenever a sequence was removed from the lists, we also eliminated all entries from the corresponding species to ensure an accurate estimation of the number of m⁶A factors per species. The sequences for each pfam domain were aligned using MAFFT (v7.310) with the “auto” setting (72). Poorly aligned regions of the resulting multiple sequence alignment were trimmed, using the Trimal program, with the parameter setting “–gappyout” (73). To determine the best protein substitution model for each pfam domain, we used a perl scrip that can be found on the following website (https://github.com/stamatak/standard-RAxML/blob/master/usefulScripts/ProteinModelSelection.pl). We selected protein substitution models GAMMALG for m⁶A writers, GAMMAJTTF for m⁶A readers, GAMMALGF for ALKBHs and VIRN, and GAMMAJTT for WTAP. Maximum likelihood trees were constructed using the RAxML program (v8.2) (74). The outgroup was set to the homologue of m⁶A factors in Arabidopsis thaliana and nodal supports were evaluated by 1,000 bootstrap replications.

To generate gene deletion mutants in the F. graminearum PH-1 strain, we employed a split marker strategy (75). This involved amplifying the left and right flanking regions of the target genes and combining them with a minimal gene cassette that carried the hygromycin phosphotransferase (HPH) gene under the control of the trpC promoter from Aspergillus nidulans. To achieve this, we conducted fusion PCR, following the previously described method (76), and the specific primers used for targeted gene deletion can be found in Table S3. In brief, we amplified the left and right flanking regions of the target genes separately, using L5 and L3 primer pairs and R5 and R3 primer pairs, respectively. The L3 and R5 primers contained 27-nucleotide (nt)-long overhang sequences that were complementary to the 5′ and 3′ ends of the minimal HPH cassette (1,376 bp in length). The HPH cassette was obtained from the pCB1004 plasmid (77), and amplified using HYG-F and HYG-R primers. Subsequently, we merged the PCR amplicons through overlap extension, assembling the left flanking region and the HPH cassette, or the right flanking region and the HPH cassette. The split marker constructs were obtained by amplifying the fused amplicons using nested primer pairs (N5 and HY-R primers for the left-half construct and YG-F and N3 primers for the right-half construct). Finally, we introduced the two split marker constructs into protoplasts through polyethylene glycol-mediated transformation (78). Following transformation, transformants resistant to 200 mg/mL of hygromycin were examined for replacement of the target gene with the HPH cassette by diagnostic PCR checks, in which L5 and R3 primers that aneal to flanking sequences of the homologous recombination event were used to confirm homologous integration of the HPH cassette into the target loci.

For the generation of MTA1-overexpressing strains, the coding sequence of MTA1 was cloned to pDS23 plasmid digested with BglII and HindIII restriction enzymes (79), using the In-Fusion HD Cloning kit (Takara Bio, Otsu, Japan). Five micrograms of the plasmid harboring MTA1 was transformed into protoplasts of the WT strain. Following transformation, transformants resistant to 200 mg/mL of nourseothricin (Jena Bioscience, Jena, Germany) were examined for introduction of an additional copy of MTA1 by diagnostic PCR check using primer pair, MTA1_OE_RT_fwd and MTA1_RT_rev. Primers used in this study are listed in Table S3.

Carrot agar plates (60 mm in diameter) (80) were inoculated by placing an agar block containing hyphae of F. graminearum at the center. The plates were incubated at room temperature under constant light. Six days after incubation, mycelia were removed by gently scraping the surface with a spatula, and then 0.9 mL of 2.5% Tween 60 (Sigma-Aldrich, St. Louis, MO, USA) was applied to the surface to assist the formation of perithecia. Sexual development in knockout mutants was observed by a stereomicroscopy for size and number formed. Squash mounts of young developing perithecia in water were examined using a compound microscope to check the morphology and maturity of ascospores.

To obtain macroconidia of the WT and MTA1-overexpressing strains, small agar blocks containing each strain were inoculated and cultured in carboxymethylcellulose medium for 3 days at 20°C (81). Macroconidia were harvested by filtration through two layers of Miracloth. Concentrations of macroconidia were adjusted to 1 × 10⁸ spores per mL, and 40 µL of spore suspensions was inoculated into 100 mL of q-PDB. Cultures were shaken at 150 rpm at room temperature. Mycelia were collected by filtration through one layer of Miracloth and oven-dried at 55°C for 2 days.

The m⁶A RNA methylation level was assessed using the EpiQuikTM m6A RNA methylation quantitative kit (Epigentek, Farmingdale, NY, USA). Briefly, 200 ng of total RNA was added and bound with the antibody in the strip wells, followed by the process of washing, capturing, and detecting antibody. Finally, the signals were detected colorimetrically by reading the absorbance at 450 nm. m⁶A levels were eventually calculated based on a constructed standard curve.

Total RNA was extracted from fresh conidia, hyphae, and perithecial tissues ground in liquid nitrogen using TRIzol reagent (Thermo Fisher Scientific) according to the manufacturer’s instructions together with the following extraction steps: two phenol (pH 4.6)-chloroform-isoamyl alcohol (25:24:1) extraction steps followed by two chloroform extraction steps after the initial TRIzol-chloroform phase separation. RNA pellets were dissolved in 88 µL of nuclease-free water and subjected to genomic DNA digestion with DNase (Qiagen Inc.). RNA samples were then concentrated using RNA Clean & Concentrator (Zymo Research). For reverse transcription (RT)-PCR analysis of MTA1, 200 ng of RNAs was reverse transcribed and amplified using the OneStep RT-PCR kit (Qiagen Inc.). MTA1 and EF1α (FGRRES_08811) were amplified for 34 cycles in the experiment for Fig. 4C, or 30 cycles in the experiment for Fig. 5A, respectively. The annealing temperature was 62°C. The primer set for amplifying EF1α was designed to amplify flanking exons, including two introns to detect possible genomic DNA contamination. For semi-quantitative RT-PCR assays, first-strand cDNA synthesis was prepared from 200 ng of total RNA using the iScript cDNA Synthesis Kit (Bio-Rad), according to the manufacturer’s instructions. Real-time RT-PCR analyses were performed using the CFX96 Touch Real-Time PCR Detection System (Bio-Rad). RT-PCR mixtures were composed of 1.5 µL of each primer (10 µM), 5 µL of SYBR Green Supermix (Bio-Rad), and 2 µL of cDNA (100 ng/mL). The PCR conditions consisted of an initial denaturing step at 95°C for 3 minutes, a denaturation step at 95°C for 10 s, both annealing and extension steps at 65°C for 30 s for 40 cycles, and 65 to 95°C with a 0.5°C increment and each temperature for 5 s to obtain the melting curve. The quantification of the relative expression levels was performed with the comparative cycle threshold method normalization (82), in that the expression level of MTA1 was normalized against the reference gene EF1α. The averages of the three biological replicates and standard deviations of the relative expression values were presented (Fig. S4). Primers used in the expression analysis are listed in Table S3.

Total RNA samples were extracted from fresh conidia of the WT, MTA1-KD1, and MTA1-OE5 strains that have been harvested from carboxymethylcellulose medium 4 days after incubation at 20°C (81). Three separate experiments were performed to produce conidia, and the conidia samples from each experiment were used as replicates for RNA-Seq analysis. Two micrograms of total RNA was sent to Macrogen Inc. (Seoul, Korea) for cDNA library construction, using the TruSeq Stranded Total RNA Library Prep Gold Kit (Illumina, San Diego, CA, USA), and for sequencing on the HiSeq4000 platform (Illumina). Raw reads (paired-end, 100 bp) were further processed and filtered, using the TrimGalore (v0.6.6) (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Filtered reads were mapped to the genome sequence of F. graminearum (NCBI accession: GCA_000240135.3), using the HISAT2 program (v2.1.0). The gene annotation file used in this study was Ensembl annotation v.32 (83). Mapped reads on genomic features, such as exon and intron, were calculated, using the htseq-count program. Gene expression levels in reads per kilobase per million mapped reads (RPKM) values were computed and normalized by effective library size estimated by trimmed mean of M values, using the edgeR R package (v3.26.8). For differential expression analysis, only genes with CPM values greater than 1 in at least three samples were kept for further analyses (9,792 out of 16,001 gene loci). Then, DE genes showing greater than a fourfold difference at an FDR of 5% were identified between the WT, MTA1-KD1, and MTA1-OE5 strains, using the limma R package (v3.28.21).

The RNA-seq data generated in the present work have been deposited in NCBI’s Sequence Read Archive and are accessible through SRA accessions from SRX21157089 to SRX21157097, which belong to the BioProject (accession, PRJNA998561). RNA-seq data for investigating the expression level changes in m⁶A factors throughout the life cycle of F. graminearum can be found in NCBI’s Gene Expression Omnibus GSE109088 for the conidial germination stages (S1–S3) and GSE109094 for the sexual development (S4–S9).