Introduction

Tubulins are major components of the microtubules that are involved in many cellular processes, such as cell division, ciliar or flagellar motility and intracellular transport in eukaryotic organisms. In general, tubulins comprise of the α-, β-, γ-, δ-, ε- and η-tubulin families1,2. The α-, β- and γ-tubulins are ubiquitous and present in all the eukaryotic organisms. The α- and β-tubulins assemble in a head-to-tail heterodimers to form the basic building block of the microtubule1. The γ-tubulins are mainly found in the microtubule organizing center and play essential roles in the initiation of microtubule assembly1,3. Interestingly, although they are well conserved in eukaryotes, fungal β-tubulins are the molecular targets of benomyl or MBC fungicides that are effective in controlling many plant diseases caused by ascomycetous fungi4,5,6. Treatments with these fungicides targeted at β-tubulins inhibit microtubule assembly and hyphal growth. Unlike the α-, β- and γ-tubulin genes, the δ-, ε- and η-tubulin genes have only been found in animals and some protists1,2. To date, their homologs have not been reported in fungi and are assumed to be lost during fungal evolution7,8. However, previous tubulin studies were mainly focused on several model organisms. With fungal species across the Kingdom Fungi being sequenced, it is important to thoroughly examine the distribution and expansion of various tubulin families in different fungal lineages.

Most animals have multiple α- and β-tubulin genes but no more than three γ-tubulin genes1,2. For example, the human genome contains at least 15 α-tubulin genes and 21 β-tubulin genes but only 3 γ-tubulin genes3. Fungi have much fewer tubulins genes than animals. Many of them contain only a single α-, β-, or γ-tubulin gene. However, some fungi such as Aspergillus nidulans9, Saccharomyces cerevisiae10, Schizosaccharomyces pombe11 and Fusarium graminearium12 are known to have two α- or two β-tubulin genes. In the budding yeast S. cerevisiae, either Tub1 or Tub3, two divergent alpha-tubulins, could perform all the functions. However, only Tub1, not Tub3, is essential for normal growth in the haploid strains13. Similarly, only one (nda2) of the two α-tubulin genes (nda2 and atb2) is essential for growth in the fission yeast S. pombe11. In the model filamentous fungus A. nidulans, the two α-tubulin genes, tubA and tubB, are functionally interchangeable but deletions of these two genes results in different phenotypes14. Two β-tubulin genes with different roles in hyphal growth and fungicide resistance also have been characterized in the wheat scab fungus F. graminearum15. These observations suggest that tubulin genes seem to undergo functional divergence after gene duplication events in different fungi. However, the underlying molecular mechanism driving the functional divergence is unknown. The evolutionary relationship among the tubulin paralogs in fungi is also not clear. Furthermore, tubulin genes are often used as the hallmark to construct the phylogenetic tree of fungi16,17,18 and identify fungal species19,20. A systematic analysis of the tubulin genes and their evolution in fungi is important for understanding the limits of these applications with various tubulin genes.

The aims of this study are to systematically analyze different families of tubulin genes and their evolutionary relationships and to understand molecular mechanism driving functional divergence of fungal tubulins. We performed a comprehensive survey of the tubulin genes across the fungal kingdom and studied the evolution of each tubulin family. The α- and β-tubulin genes were found to undergo multiple independent duplications and losses in different fungal lineages and formed distinct paralogous and orthologous clades. Molecular evolutionary analysis indicated that α1, α2 and β2-tubulins have been under strong divergent selection pressure and adaptive positive selection. Amino acid residues that are likely involved in functional diversification were identified by functional divergence analysis. Many of the positively selected sites are at or adjacent to functional important sites and likely contribute to functional diversification. Moreover, we experimentally confirmed the divergent functions of two β-tubulin genes in F. graminearum that may interact with different microtubule associate proteins and identified type II variations in FgTub2 that are likely responsible for its functional shifts. Our results are important for better understanding of the evolution of fungal tubulin genes and their mutational dynamics under natural selection in fungi.

Results and Discussion

Distribution of fungal α-, β- and γ-tubulin genes

To systematically analyze fungal tubulin genes, we selected 59 fungi that have genome sequences available to the public and are representative of main fungal phyla of Ascomycota, Basidiomycota, Zygomycota and Chytridiomycota (Fig. 1; Supplementary Table S1). All these fungal species examined in this study have α-, β- and γ-tubulin genes although they differ in the copy numbers of each tubulin gene families (Fig. 1). While γ-tubulin is at the minus-end of microtubules, α- and β-tubulins form chains of heterodimers in the microtubule1,3. Majority (42 out of 53) of ascomycetes and basidiomycetes contain two α-tubulin genes. In contrast, the three Chytridiomycota species examined has only a single α-tubulin gene but the Zygomycota Rhizopus oryzae and Rhizomucor miehei each has 5 α-tubulin genes and another Zygomycota Mucor circinelloides has 3 α-tubulin genes (Fig. 1). For the β-tubulin genes, majority of ascomycetes (28 out of 38) contain one but all basidiomycetes examined have at least two (Fig. 1). All the three Zogomycota and three Chytridiomycota species each contains at least 3 β-tubulin genes. Unlike the α- and β-tubulin genes, all these fungi have a single γ-tubulin gene except Candida albicans and Allomyces macrogynus that has two (Fig. 1).

Figure 1
figure 1

Identification of tubulin genes in 59 representative fungi.

The species tree was drawn based on the phylogeny of a-tubulins. The copy number of each tubulin type was indicated. For the definition of α1, α2, α3, β1 and β2 see the text.

Full size image

Zoospore-producing chytridiomycetes retain tubulins functionally related to flagella

Members of the δ-, ε- and η-tubulin families were thought to be lost in fungi7,8. However, we found that all three chytridiomycetes, A.macrogynus, Spizellomycespunctatus and Batrachochytriumdendrobatidis, have ε-tubulin genes (Fig. 1; Supplementary Fig. S1). A. macrogynus also has δ- and η-tubulin genes but B.dendrobatidis has only one η-tubulin gene. The zygomycetes, ascomycetes and basidiomycetes lack orthologs of δ, ε, or η tubulins that are associated with the centrioles and basal bodies in other eukaryotic organisms21,22. Among the true fungi, only chytridiomycetes produce zoospores with flagella. Therefore, the δ, ε, or η-tubulins in chytridiomycetes may be involved in the assembly of the basal bodies (kinetosomes), which is required for the anchorage of flagella. Other fungi do not produce zoospores23 and lack developmental stages with a flagellum24. It is likely that the δ-, ε- and η-tubulins were present in the ancestral fungi but subsequently lost in other fungi after their divergence from zoospore-producing chytridiomycetes due to the lack of selection pressure on flagellum formation.

Multiple independent origins of α-tubulin paralogs in fungi

Phylogenetic analysis showed that α-tubulin genes from ascomycetes formed two distinct clades (Fig. 2). Members of one clade (named α1-tubulins) were found in all ascomycetes examined and clustered with α1-tubulins from basidiomycetes. The other clade (named α2-tubulins) formed a sister group to the α1-tubulins and basidiomycetes α-tubulins, suggesting that it present in the last common ancestor of ascomycetes and basidiomycetes. Member of this clade were lost in fungi belonging to Basidiomycota, Saccharomycotina and Taphrinomycotina during subsequent evolution.

Figure 2
figure 2

Phylogeny of fungal α-tubulins.

The phylogenetic tree was constructed with the Maximum-likelihood (ML) approach. It had similar topology with the tree generated with the Bayesian inference (BI) method. Numbers on major branches indicate SH-like approximate likelihood ratio test (SH-aLRT) probabilities/Bayesian posterior probabilities. Wathet and blue branches indicate α1- and α2-tubulins of ascomycetes, respectively. Green and red branches indicate α1- and α3-tubulins of agaricomycetes. The tubulin genes highlighted in yellow have been functionally characterized in previous studies11,13,14. For abbreviations, see Additional file Table S1 and supplementary Fig. S1 legend.

Full size image

The ascomycetes α1- and α2-tubulins were further duplicated in several species during evolution. For examples, Tub1 and Tub3 of S. cerevisiae13 were clustered together in the α1-tubulin clade (Fig. 2), suggesting that they were originated from species-specific gene duplication. The essential (nad2) and non-essential (atb2) alpha-tubulin genes from the fission yeast S. pombe11 also belong to the α1-tubulin clade and were clustered with their counterparts from S. japonicus (Fig. 2), suggesting that duplications occurred in their last common ancestor. The two previously reported α-tubulins of A. nidulans, TubA and TubB14, belong to the α1- and α2-tubulin clades, respectively. Moreover, we found A. nidulans has two additional copies of α2-tubulin genes (Fig. 2). Among the other ascomycetes examined, Botryotinia fuckeliana and Penicillium chrysogenum also have one extra copy of α2-tubulin gene. These results suggest that frequent duplications of the α1- or α2-tubulin gene occurred in different lineages of ascomycetes during evolution.

In basidiomycetes, most of the rust and smut fungi have a single α-tubulin gene except Melampsora laricis-populina that has two (Fig. 2). In agaricomycetes, a gene duplication event occurred in their last common ancestor, which led to a new clade of α-tubulins (named α3-tubulins) in fungi belonging to Agaricomycotina (Fig. 1). Fomitiporia mediterranea and Cryptococcus neoformans are the only two agaricomycetes that may have lost the α3-tubulin gene during evolution (Fig. 2).

Interestingly, R. oryzae has five α1-tubulin genes that clustered together (Fig. 2), suggesting that they are originated from recent gene family expansions. Indeed, it has been reported that the R. oryzae genome underwent a high level of duplication25. However, all three fungi belonging to Chytridiomycota have a single α-tubulin gene (Fig. 1; Fig. 2). Therefore, it is likely that the common fungal ancestor may contain only one α-tubulin gene. Multiple independent gene duplication events have occurred in different fungal lineages on varying evolutionary timescales, resulting in the formation of distinct α-tubulin paralogs.

Most ascomycetes have only β1-tubulins but basidiomycetes have β1 and β2-tubulins

Phylogenetic analysis showed that β-tubulins underwent complex gene duplication and gene loss events in different fungal lineages (Fig. 3). In lower fungi, each species encode multiple copies of β-tubulins and most of them were originated by species-specific gene duplications. The two β-tubulin genes of basidiomycetes formed two distinct clades named β1- and β2-tubulins. In contrast, most of the ascomycetes have a single β-tubulin gene that forms the third clade (Fig. 3). Interestingly, two additional β-tubulins from ascomycetous yeast Yarrowia lipolytica and Galactomyces geotrichum, were clustered with the β2-tubulins of basidiomycetes. The branch length of the β2-tubulins is remarkably long, suggesting that accelerated evolution occurred in these genes. It is likely that the duplication event lead to the evolution of β2-tubulins occurred before the split of ascomycetes and basidiomycetes although β2-tubulin was lost in most ascomycetes. However, it remains possible that Y. lipolytica and G. geotrichum gained the β2-tubulin gene from basidiomycetous species by horizontal gene transfer (HGT).

Figure 3
figure 3

Phylogeny of fungal β-tubulins.

The phylogenetic tree was constructed with the Maximum-likelihood approach. Numbers on major branches indicate SH-like approximate likelihood ratio test (SH-aLRT) probabilities/Bayesian posterior probabilities. Green and red branches indicate β2- and β1-tubulins of basidiomycetes, respectively. The tubulin genes highlighted in yellow have been functionally characterized in previous studies12,15. For abbreviations, see supplementary Fig. S1 legend and Additional file Table S1.

Full size image

In filamentous ascomycetes, the phylogenetic pattern of additional β1-tubulin genes is complicated (Fig. 3). The two β1-tubulin genes were clearly originated from gene duplication events in some fungi such as Tuber melanosporum and Cordyceps militaris. However, additional β1-tubulin genes were more likely derived from HGT events in other fungi. For example, one of the two β1-tubulin genes in Metarhizium acridum (NCBI gene ID: 322700420) clustered with its homolog from C. militaris (Fig. 3), which is not consistent with their species phylogeny. One reasonable explanation is that this β1-tubulin gene was originated from a HGT event between these two entomopathogenic species. Similarly, the additional β1-tubulin genes in Glomerella graminicola (310790380 and 310801747) and Fusarium solani (302882768 and 302905990) may be generated by HGTs (Fig. 3).

Gamma-tubulins are well-conserved in other fungi but highly divergent in saccharomycetes

Although most fungi examined encode a single γ-tubulin gene, A. macrogynus and C. albicans have two γ-tubulin genes. Phylogenetic analysis suggested that these two organisms gained the additional γ-tubulin gene through species-specific gene duplication events (Supplementary Fig. S2). In addition, we found that γ-tubulins from saccharomycetes are highly divergent in comparison with their orthologs from other fungi (Supplementary Fig. S2), which may be due to their adaption to the specialized microtubule attachment structures in yeasts26.

Tubulin-like proteins misatos are conserved in all the fungi except pucciniomycetes

Because misatos share sequence similarity with tubulins and contain some structural features similar to those of tubulins such as N-terminal GTP-binding sites27,28, we also investigated their evolution in fungi. In general, each fungus contains a single misato gene in their genomes and the phylogenetic tree of misatos is consistent with the fungal species phylogeny (Supplementary Fig. S3), suggesting that they are highly conserved in fungi. Surprisingly, misatos were absent in the genomes of sequenced pucciniomycetes, including M. laricis-populina, Puccinia graminis, P. triticina, and P. striiformis. In S. cerevisiae, misato is functionally related to mtDNA inheritance and partition of mitochondria27. In Drosophila, misato is required for the kinetochore-driven microtubule formation29. We also searched the genome of powdery mildew fungus Blumeria graminis and found that it contains a misato gene, thus the lack of misato gene is unlikely related to the obligate biotrophic lifestyle. The lack of misatos in Pucciniomycotina may be linked to the specific life cycle of rust fungi.

Both the α1 and α2-tubulins in ascomycetes are under divergent selective pressure

Gene duplications could shift the functional constrains and give rise to functional divergence30. Many tubulins showed functional divergence after gene duplication and account for different microtubular structures in animals1,31,32,33. To understand evolutionary mechanisms of duplicated α-tubulin paralogs in fungi, we used codon-based models implemented in PAML34 to analyze the selective pressure acting on different tubulin paralogs by estimating the ratio (ω) of nonsynonymous (dN) to synonymous (dS) substitutions rates (ω = dN/dS). We found that the mean ω estimated using one-ratio model (M0) was 0.03 for α-tubulins (Table 1), indicating that they are under strong purifying selection during evolution, which is in accordance with their important biological roles in cellular processes. However, results from the two-ratio branch model analysis suggested that the ω values of branches leading to α1- and α2-tubulins of ascomycetes and α1- and α3-tubulins of agaricomycetes differ significantly from those of the background branches (Table 1), indicating asymmetrical evolution after gene duplication.

Table 1 Parameter Estimates of Codon-Substitution Evolutionary Models for α-Tubulins
Full size table

To examine if the α-tubulins are under divergent selective pressure, we detected the selective pressure among paralogous clades of the α-tubulin phylogeny using Clade model D30 implemented in PAML. Moreover, functional divergence of gene could be related to the site-specific change in protein sequence during evolution35. Thus, we also examined type I36 and type II functional divergence (FD)36 after duplications of α-tubulins using DIVERGE36. Results of clade model D analysis showed that the selective pressure acting on α3-tubulins was similar to that acting on α1-tubulins in agaricomycetes (Table 1). Moreover, both α1- and α3-tubulins showed no difference in comparison with single copy α-tubulins of other basidiomycetes (Table 1). Consistent with this analysis, type I or type II divergence was not detected between α1- and α3-tubulins of agaricomycetes or α-tubulins of other basidiomycetes (Supplementary Table S2). These results together with data from clade model D analysis suggest that α1- and α3-tubulins in agaricomycetes lack significant functional divergence.

In contrast, results from clade model D suggested that α1- and α2-tubulin paralogs are under divergent selective pressure in ascomycetes (Table 1). In comparison with the α-tubulins from chytridiomycetes and basidiomycetes, the divergent selective pressure is stronger on α2-tubulins than on α1-tubulins (Table 1). Furthermore, significant type I and type II functional divergences was identified between α1- and α2-tubulins of ascomycetes (Supplementary Table S2). In total, we identified 19 sites that are likely responsible for type I functional divergence. These sites are well conserved in one duplicate cluster but highly variable in the other duplicate cluster, suggesting that they may have experienced the shifted functional constrains after gene duplication. We also identified 35 sites that are likely responsible for type II functional divergence (Supplementary Table S2). These sites are highly conserved in both duplicate clusters but vary between the two clusters. Some of them may be involved in functional specifications. In comparison with α-tubulins from chytridiomycetes and basidiomycetes, α2-tubulins have strong type I and type II FD in ascomycetes but α1-tubulins have only weak type I FD and lack type II FD. Results from both clade D model and DIVERGE analyses suggest that the two α-tubulin clades in ascomycetes may both have undergone functional divergence and divergent selective pressure but the driving force is more intense on α2-tubulins than on α1-tubulins.

The β1- and β2-tubulins are under strong purifying selection and divergent selective pressure, respectively, in basidiomycetes

The ω value estimated using one-ratio model (M0) was 0.028 for β-tubulins (Table 2), indicating that they are also under strong purifying selection during evolution. Although the two ratio model showed no divergent selective pressure acting on the branch leading to the β2-tubulins, results from the clade model analysis indicated that there is significant divergence in selective pressure between β1- and β2-tubulins (Table 2). Furthermore, the β1- and β2-tubulin genes were found to be under type I functional divergence in basidiomycetes. In comparison with β-tubulin genes of ascomycetes, both β1- and β2-tubulin genes of basidiomycetes are under divergent selective pressure (Table 2). However, unlike β2-tubulins, β1-tubulins have a lower ω value than the background in basidiomycetes. In addition, the β2- but not β1-tubulins have significant type I variations. These results suggest that whereas the β1-tubulin genes are under purifying selection, the β2-tubulin genes are under divergent selective pressure in basidiomycetes.

Table 2 Parameter Estimates of Codon-Substitution Evolutionary Models for β-Tubulins
Full size table

Positive selections act on the functionally diverged tubulin paralogs

To determine whether positive selection (PS) drive functional divergence, the branch-site model (model A) was used to detect the selective pressure acting on pre-specified branches leading to α1-, α2-, or β2-tubulins, respectively. We found that branches leading to ascomycetous α1-, α2-, or β2-tubulins contain positively selected amino acid sites (Table 3). In total, 87 sites evolving under positive selection were identified in the branch leading to α2-tubulins, but only 21 were detected in the branch leading to α1-tubulins, which is consistent with their different extents of functional divergence (Supplementary Fig. S4). We also identified 18 sites that are under positive selection in the branch leading to β2-tubulin. These results suggested that positive selection contributes to the functional divergence of α1-, α2 and β2-tubulins.

Table 3 Results of Model A analyses with α1 and α2-tubulins of ascomycetes, β2-tubulins of basidiomycetes and β2-tubulins of Fusarium
Full size table

Many PS or FD sites are involved in or adjacent to functionally important sites

We then examined the distribution of amino acid sites that were identified as PS or under FD in tubulins. The α1- and α2-tubulins from ascomycetes shares 10 PS sites. Three of them, including amino acid residues at sites 28, 235 and 330 (Supplementary Fig. S4) are conserved in α1- or α2-tubulins but vary between the α1- and α2-tubulin clades. These common PS sites may contribute significantly to the functional divergence. Thirty-four sites contributing to the type I or type II FD also were identified to be under positive selection. For examples, in α2-tubulins from ascomycetes, type I FD site 12 (adjacent to the GTP binding site 13) and type II FD sites 86 (GTP binding site) and 88 (α- and β-interface) are under positive selection.

Only seven PS or FD sites, such as site 88 in ascomycetous α2-tubulins, are directly involved in tubulin functions by affecting GTP binding and the interface between α- and β-tubulins (Fig. 4; supplementary Fig. S4 and S5) and these mutations may affect the assembly of α/β-tubulin heterodimer. Fifteen PS or FD sites are located next to functional important residues (Fig. 4; supplementary Fig. S4 and S5). For example, PS sites of amino acids 12, 19, 185 and 224 in α2-tubulins are adjacent to residues that are involved in GTP binding. Moreover, a number of PS or FD sites are located in the C-terminal domain, which probably constitutes the binding surface for microtubule associated proteins (MAPs) and motor proteins1,37. Therefore, functional divergence may not affect core functional or structural features but allows duplicated tubulin genes to acquire new peptide binding abilities.

Figure 4
figure 4

The 3D structures of α- and β-tubulins showing positively selected sites.

(A) Positively selected sites in α2-tubulins. (B) Positively selected sites in β2-tubulins. Sites are mapped onto PDB ID: 4FFB, with α-tubulin in green and β-tubulin in blue. Only positively selected sites that are at or adjacent to sites involved in GTP binding and α/β connection are shown. Sites related to the tubulin function and/or structure were obtained according to the conserved domains of α-tubulin (cd02186) and β-tubulin (cd02187) in CDD of NCBI. The Arabic numbers refers to the amino acid positions in the sequence alignments of α-tubulins and β-tubulins presented in supplementary Fig. S4 and S5, respectively.

Full size image

Different evolutionary scenarios for functional diversification of α-tubulin paralogs in ascomycetes and β-tubulin paralogs in basidiomycetes

Both α1- and α2-tubulins from ascomycetes were found to contain substantial PS sites and have functional divergence in comparison with α-tubulins from basidiomycetes (Supplementary Table S2, Supplementary Fig. S4). Therefore, the ancestors of either α1- and α2-tubulins may have undergone the first phase of ancestral function sub-functionalization under divergent selective pressure and subsequently gained new function through the fixation of adaptive amino acid changes under positive selection. Because both the divergent selective pressure and positive selection acting on α2-tubulins are more intense than those acting on α1-tubulins, α1-tubulins may have retained more ancestral functions of α-tubulins. This observation is supported by the fact that the α1-tubulins are constitutively expressed in different developmental stages and play more critical roles in fungi11,13,14.

Of the two clades of basidiomycetous β-tubulins, the β1-tubulin genes are much more conserved than the β2-tubulin genes that have been evolving rapidly (Fig 3), suggesting that the duplication of β-tubulin genes relaxed constrains only on one of them. Moreover, we found that the β1-tubulins are under strong purifying selection but the β2-tubulins are under the divergent selective pressure and positive selection in basidiomycetes (Table 2). Thus, the β1-tubulin genes likely have retained and/or improved the ancestral β-tubulin functions. In contrast, the β2-tubulin genes may explore or improve the sub-functions of ancestral β-tubulins and possibly gain novel functions.

Two β-tubulins in Fusarium are under purifying selection and divergent selective pressure, respectively

The ascomycetes F. graminearum has two β-tubulin genes, FgTub1 and FgTub2 (Fig. 3). Previous studies showed that, of the two β-tubulin genes, mutations in FgTUB2 but not in FgTUB1 have been shown to confer resistance to benomyl or MBC fungicides12. Furthermore, the two β-tubulin genes differ in their functions in hyphal growth15. To investigate the functional divergence of these two β-tubulin genes in F. graminearum, we then identified their orthologs in other sequenced Fusarium species. Phylogenetic analysis showed that these β-tubulins formed two distinct clades (Supplementary Fig. S6). Whereas the Fusarium Tub1 or FgTub1 clade is in accordance with the species tree, the Tub2 clade falls to the base of the Sordariomycete cluster, suggesting that the Tub1 clade consists of the canonical tubulin gene while Tub2 is the additional copy.

Results of clade D model showed that both Tub1 and Tub2 proteins are under divergent selective pressure in comparison with β-tubulins from other ascomycetes. However, members of the Tub1 clade have a lower ω value than the background. In contrast, members of the Tub2 clade have a higher ω value than the background (Table 2). These results suggest that the Fusarium Tub1 tubulins are under strong purifying selection but the Tub2 tubulins are under divergent selective pressure. The branch-site model failed to identify positively selected sites in the branch leading to Tub2 tubulins, which may be due to the limitation of ML methods in analyzing high similar sequences38.

FgTub1 and FgTub2 play distinct roles in vegetative growth and sexual reproduction in F. graminearum

To experimentally confirm the functional divergence of two β-tubulin genes in F. graminearum, we deleted the FgTUB1 and FgTUB2 genes respectively in F. graminearum PH-1 strain and examined the phenotypes of the Fgtub1 and Fgtub2 deletion mutants. In comparison with the wild type PH-1, growth rate was slightly reduced in the Fgtub1 mutant but significantly reduced in the Fgtub2 mutant (Fig. 5A). However, the Fgtub1 mutant was blocked in ascosporogenesis although it still produced small perithecia that were sterile. The Fgtub2 mutant was slightly reduced in perithecium formation but produced normal perithecia. It had no obvious defects in the production of asci and ascospores (Fig. 5A). These results showed that FgTUB1 and FgTUB2 differ significantly in their functions during sexual reproduction, which was not reported in the previous study15. Overall, it appears that FgTUB1 plays an essential role in sexual development but FgTUB2 is more important for vegetative growth, confirming that functional divergence occurred between these two β-tubulin genes in F. graminearum.

Figure 5
figure 5

Phenotypes of the Fgtub1 and Fgtub2 deletion mutants of Fusarium graminearum.

(A) Three day-old PDA cultures, perithecia and asci in cracked perithecia of the wild type strain (PH-1) and the Fgtub1 and Fgtub2 deletion mutants. White bar = 1 mm; Black bar = 50 µm. (B) The subcellular localization of FgTub1-GFP and FgTub2-GFP in the wild type strain PH-1 and Fgkin1 deletion mutant. In PH-1, both FgTub1 and FgTub2 localized to the microtubules. In the Fgkin1 mutant, FgTub1 became aggregated in the nucleolus but FgTub2 still localized to the microtubules. Bar = 10 μm.

Full size image

Tub2 and Tub1 in Fusarium may have different interactions with MAPs

The DIVERGE examination showed that, unlike Tub1 and other canonical β-tubulins of ascomycetes, Tub2 has significant type II variations (Fig. 6 and Supplementary Table S2). Interestingly, four substitutions occurred at amino acid residues that are at or adjacent to the GTP binding sites or the interface between α- and β-tubulins (Fig. 7A). These changes may affect the interactions between FgTub2 and α-tubulin monomers. Notably, four residues under type II variations, A153, A291, M331 and N424, are located on the surface of FgTub1. Their counterparts in FgTub2, S153, H291, T331 and H424 (Fig. 7B) are still on the protein surface but undergo dramatic changes in amino acid properties. For example, residue 424 is the polar, uncharged N in FgTub1 but the positively charged H in FgTub2. These changes may affect the interaction of FgTub2 with other proteins such as FgTub1- or FgTub2-specific MAPs.

Figure 6
figure 6

Amino acid sites under Type II functional divergence in Fusarium Tub2 tubulins.

The names of Tub1 tubulins are in blue and those of Tub2 are in red. Residues identified in multiple sequence alignment to vary among Tub1 and Tub2 tubulins are highlighted in different colors. Red, blue and green arrows mark residues involved in GTP binding, α/β interface and β/α interface, respectively. “+” symbols mark sites that contribute to Type II functional divergence. Sites related to the tubulin function and/or structure were obtained according to the conserved domains of α-tubulin (cd02186) and β-tubulin (cd02187) in CDD of NCBI, respectively. Those boxed in red are amino acid sites reported to be associated with the fungal resistance to fungicide12,41,42. Abbreviations: Fg, Fusarium graminearum; Fv, F. verticillioides; Fo, F. oxysporum; Ff, F. fujikuroi; Fs, F. solani; Fp, F. pseudograminearum.

Full size image
Figure 7
figure 7

The 3D-structures of FgTub1 and FgTub2 tubulins showing amino acid residues contributed to the type II functional divergence.

(A) Type II sites directly related to or adjacent to GTP binding or α/β connection. (B) Distribution of the four type II sites (highlighted in red) on the surface of FgTub1 and FgTub2. The 3D structures of FgTub1 and FgTub2 were predicted with HHpred58 using the HHblits multiple sequence alignment method.

Full size image

In the mutant deleted of the FgKIN1 kinase gene that likely regulates microtubule stability via phosphorylation of MAPs39, the subcellular localization of FgTub1 but not FgTub2 was affected40. Instead of localization to the microtubules, FgTub1 was over-aggregated in the nucleolus in the Fgkin1 mutant (Fig. 5B). These results suggested that the interactions of FgTub1 and FgTub2 with MAPs have been diverged and the FgKin1 kinase may specifically regulate FgTub1 via its specific MAPs. Interestingly, of the five amino acid sites reported to be associated with resistance to benomyl12,41,42, three are conserved among Tub1 and Tub2 proteins (Fig. 6). Site 81 differ between Tub1 and Tub2 in different Fusarium species.

Overall, results from our evolutionary analysis and experimental characterization suggested that the two β-tubulins in Fusarium are likely functional divergent and residues under type II variations most probably contribute to their functional divergence. These, together with functional divergence occurred between α-tubulins that have been reported in S. pombe, S. cerevisiae and A. nidulans11,13,14 support that different paralogs of α- or β-tubulin genes are functionally divergent and divergent selective pressure and positive selection likely account for their functional divergence.

Implications for tubulin genes in the reconstruction of fungal species tree

Tubulin sequences are frequently used to construct the fungal tree of life16,17,18. Reconciliation of the different tubulin gene trees with their species trees revealed multiple gene duplication and loss events in α-/β-tubulin evolution (Supplementary Fig. S7), raising cautions for their applications in the reconstruction of fungal species tree. α-tubulins can be used for the reconstruction of fungal species tree but sequences from different paralogous clades should be avoided. β-tubulin genes are not suitable for the species tree reconstruction of filamentous ascomycetes where they underwent a complicated recent gene duplication/loss events. Intriguingly, γ-tubulin gene tree fit well with the fungal tree of life and most fungi contain only a single γ-tubulin gene, thus, it may be more suitable for studying the fungal species evolution.

Conclusions

In conclusion, we systematically identified fungal tubulins and analyzed their distribution and evolution in 59 representative fungi across the fungal kingdom in this study. We reported the identification of tubulin genes belonging to the δ, ε and η-tubulin sub-families in three Chytriodiomycetes A.macrogynus, S.punctatus and B. dendrobatidis, which have not been reported in fungi. Phylogenetic analysis showed that α- and β-tubulin genes underwent multiple independent duplications and losses in ascomycetes and basidiomycetes and formed distinct paralogous and orthologous clades. Our results indicate that the last common ancestor of basidiomycetes and ascomycetes likely possessed at least two paralogs each of the α-tubulin and β-tubulin genes (α1 and α2; β1 and β2). Whereas the α2-tubulin genes were lost in all the ascomycetes, the β2-tubulin genes were retained in a few ascomycetous yeasts during evolution. Molecular evolutionary analysis indicated that the α1, α2 and β2-tubulin genes have been under strong divergent selection pressure and adaptive positive selection. Many of these positively selected sites are directly involved or adjacent to functionally important sites such as GTP binding and the α/β or β/α interface and likely contribute to their functional diversification. Our results suggest that α-tubulins in ascomycetes and β-tubulins in basidiomycetes have different evolutionary mechanisms. Amino acid residues identified in this study to be under PS or FD may be involved in function shifts of fungal tubulin genes. Further characterization of the effects of genetic variations at these sites will be important to determine functional differences between various tubulin paralogs in fungi. Furthermore, we characterized the phenotypes of mutants deleted of two β-tubulin genes in F. graminearum in sexual development and hyphal growth. Our experimental data confirmed the stage-specific functional differences between these two β-tubulin paralogs and we identified a number of sites that are under type II divergence and may be responsible for functional specifications. Overall, results from this study indicated that different molecular mechanisms are responsible for driving functional diversity in fungal tubulin genes. It will be important to functionally characterize the amino acid residues that are identified in this study to be related to functional divergence in fungal tubulins.

Methods

Data collection and identification of fungal tubulins

For the predicted proteomes of 59 fungi used in this study, 47 were obtained from the GenBank of NCBI, 9 were downloaded from the Fungal Genome Initiative (FGI) site at the Broad Institute (http://www.broadinstitute.org/science/projects/projects) and 3 were downloaded from the DOE Joint Genome Institute (http://genome.jgi.doe.gov/programs/fungi/index.jsf). Tubulins of G. geotrichum and T. viride reported in previous studies43,44 were obtained from the GenBank.

We used the Hmmsearch program in HMMER 3.045 to search for the family-specific HMM profiles of tubulins (PF00091.20) and Misatos (PF10644_misato) downloaded from Pfam database46 with each of fungal predicted proteomes as queries. Fungal proteins identified as tubulins were then grouped into different families by phylogenetic analysis.

Phylogenetic analyses

Multiple sequence alignments were performed using the M-Coffee program, which combines the output of popular aligners47. Sequence alignments used for phylogenetic analysis were trimmed by trimAL48 with the gappyout model.

We used the Maximum likelihood (ML) and Bayesian inference (BI) methodologies to construct phylogenetic trees. The ML trees were constructed with PhyML 3.149 using the best-fit model selected by ProtTest350, with SPRs algorithms and 16 categories of gamma-distributed substitution rates. The reliability of internal branches was evaluated with SH-aLRT supports. The BI tree was constructed with MrBayes-3.251 using mixed models of amino acid substitution with 16 categories of gamma-distributed substitution rates, performing two runs for each of four Monte Carlo Markov Chains (MCMCs). Trees were sampled every 1000th iteration over 1.1 million generations after a burn-in of 101 samples.

The reconciled trees were built using most parsimonious reconciliations methods and viewed by Primetv52. The species tree was built based on the alpha tubulin gene tree and manually edited according to the previous studies17,53,54.

Detection of positive selection in tubulins

The CODEML program from the PAML 4.7 package34 was used to estimate the ω values of tubulins and test for positive selections in tubulin foreground lineages versus the background lineages. To detect different selective pressures that affect specific lineages, we used one-ratio model (model 0) and two-ratio model (model 2). Model 0 assumes that the ω ratio is constrained along all branches in the phylogeny, while model 2 allows a different ω ratio for the foreground lineages. To detect the different selective pressure affecting specific clades in which functional divergence may have occurred after gene duplication, we compared the site-specific discrete model 3 that assumes two classes of sites with different ω ratios and the clade model D that allows selective pressure at one class of sites (foreground clade) to be different from the rest of the phylogeny30. We then detected positive selection that affects some sites on the specific lineages by comparing the branch-site model A, assuming one class of sites of the foreground lineage ω >1 and null model A with fixed ω = 155. Trees used for PAML analysis were presented in additional files as supplemental figures S8–S11.

Protein modeling and 3D visualization

The 3D-structural models of α- and β-tubulins were presented according to that of αβ-tubulins of S. cerevisiae (PDB ID: 4FFB)56 and displayed with Chimera 1.8.157. The 3D structures of F. graminearum FgTub1 and FgTub2 were predicted with HHpred58 using the HHblits multiple sequence alignment method.

Culture conditions and phenotype characterization of F. graminearum mutants

All the stains used in this study were cultured at 25°C on Potato Dextrose Agar (PDA) as described59. Colony morphology and growth rate were assayed with PDA cultures as described60. The mutant of FgTub1 and FgTub2 were generated through the split-marker approach as described60. For sexual reproduction assays, aerial hyphae of 6-day-old carrot agar cultures were pressed down with 0.5 ml of sterile 0.1% Tween 20. Perithecium formation and ascospore production were assayed after incubation at 25°C for another 2 weeks.