Broad Genomic Sampling Reveals a Smut Pathogenic Ancestry of the Fungal Clade Ustilaginomycotina

Abstract

Ustilaginomycotina is home to a broad array of fungi including important plant pathogens collectively called smut fungi. Smuts are biotrophs that produce characteristic perennating propagules called teliospores, one of which, Ustilago maydis, is a model genetic organism. Broad exploration of smut biology has been hampered by limited phylogenetic resolution of Ustilaginiomycotina as well as an overall lack of genomic data for members of this subphylum. In this study, we sequenced eight Ustilaginomycotina genomes from previously unrepresented lineages, deciphered ordinal-level phylogenetic relationships for the subphylum, and performed comparative analyses. Unlike other Basidiomycota subphyla, all sampled Ustilaginomycotina genomes are relatively small and compact. Ancestral state reconstruction analyses indicate that teliospore formation was present at the origin of the subphylum. Divergence time estimation dates the divergence of most extant smut fungi after that of grasses (Poaceae). However, we found limited conservation of well-characterized genes related to smut pathogenesis from U. maydis, indicating dissimilar pathogenic mechanisms exist across other smut lineages. The genomes of Malasseziomycetes are highly diverged from the other sampled Ustilaginomycotina, likely due to their unique history as mammal-associated lipophilic yeasts. Despite extensive genomic data, the phylogenetic placement of this class remains ambiguous. Although the sampled Ustilaginomycotina members lack many core enzymes for plant cell wall decomposition and starch catabolism, we identified several novel carbohydrate active enzymes potentially related to pectin breakdown. Finally, ∼50% of Ustilaginomycotina species-specific genes are present in previously undersampled and rare lineages, highlighting the importance of exploring fungal diversity as a resource for novel gene discovery.

phylogenomics, tree expansion analysis, CAZymes, undersampled lineages, Exobasidiomycetes, Malasseziomycetes

Introduction

With 5.1 million estimated species and an age of ∼1 By (Taylor and Berbee 2006; Blackwell 2011), Fungi represents one of the most ancient and diverse eukaryotic kingdoms. They vary widely in nutritional modes, inhabit diverse environments, and range from mutualists to pathogens of algae, plants, animals, and other fungi (Blackwell 2011). Currently, one of the most powerful approaches for making sense of fungal evolution and diversity is to combine comparative genomics with phylogenomics (Hibbett et al. 2013). While sequencing efforts have enabled this strategy for many major groups within Fungi, sampling for some lineages, such as Ustilaginomycotina (Basidiomycota)—historically known as the lineage containing the true “smut” fungi (Begerow et al. 2014)—remains poor. For example, there are nearly 800 publicly available fungal genomes in the JGI MycoCosm portal (Grigoriev et al. 2014), but as of March 2017, only 11 Ustilaginomycotina genomes, most of which are from the single order Ustilaginales, were available prior to this study. As such, extensive exploration within this subphylum has thus far been severely limited.

Smut fungi are one of the most important groups of plant pathogenic fungi. As a group they are able to infect numerous economic crops, especially cereals (Begerow et al. 2014). Exemplars are Ustilago maydis (corn smut), Tilletia spp. (common bunt on wheat), Sporisorium reilianum (smut on maize and sorghum), and S. scitamineum (sugarcane smut). Most smut fungi related to the model species U. maydis (Ustilaginales) have a complex dimorphic life cycle (supplementary fig. S1, Supplementary Material online). They are saprotrophic and yeast-like (termed sporidia) at first, reproducing asexually by budding in the haploid stage (n). Only following fertilization, which occurs between two compatible sporidia to form dikaryotic (n + n) filamentous hyphae, is the pathogenic phase of the life cycle enabled. The hyphae then colonize and exploit their host through biotrophic strategies. Smut fungi ultimately produce teliospores on their host for dispersal and perennation. Under favorable environmental conditions the teliospores germinate to form basidia, where two nuclei fuse to become diploid (2n), and then undergo meiosis to produce haploid basidiospores (n), which give rise to the sporidia that begin the cycle again. However, members of Tilletiales, Doassanssiales, and some members of Urocystidiales appear to lack the yeast-like state (Begerow et al. 2014). While smut fungi in Ustilaginales have received much attention as models for studying fungal dimorphism and biotrophic pathology, much less is known about the biology of other members of Ustilaginomycotina or about how pathogenic traits may be conserved across the subphylum.

At present, Ustilaginomycotina is divided into four classes and 15 orders (table 1). The smut fungi are found in seven Ustilaginomycotina orders (dispersed across two classes): Ustilaginales, Uleiellales and Urocystidales (Ustilaginomycetes) and Tilletiales, Doassansiales, Georgefischeriales and Entylomatales (Exobasidiomycetes) (Bauer et al. 1997; Riess et al. 2016). Two of these, Tilletiales and Ustilaginales, represent some of the most successful and diverse smut fungal groups, the members of which are specialized upon graminaceous plant hosts. In addition to smuts, Ustilaginomycotina hosts fungi with a wide range of other ecological and nutritional strategies. These include nonsmut plant pathogens (such as gall-forming fungi in the genus Exobasidium), animal-associated lipophilic fungi (such as dandruff fungi in the genus Malassezia) and anamorphic yeasts for which no pathogenic stage is known (such as members of the genera Pseudozyma and Tilletiopsis) (Bauer et al. 1997; Begerow et al. 2000, 2006). Interestingly, these varying nutritional strategies appear interspersed across Ustilaginomycotina, raising compelling (and controversial) questions regarding how smut characteristics and nutritional strategies have evolved within this subphylum (Begerow et al. 2014; Riess et al. 2016).

Table 1.

Ustilaginomycotina Higher Level Classification.

Taxonomic Rank	Name	Reference
Subphylum	Ustilaginomycotina	Bauer et al. 2006
True smut lineages
Class	Ustilaginomycetes	Bauer et al. 1997
Order	Ustilaginales	Clinton 1906; Bauer et al. 1997
	Urocystidales	Bauer et al. 1997
	Uleiellales	Riess et al. 2016
Class	Exobasidiomycetes	Begerow et al. 2006
Order	Entylomatales	Bauer et al. 1997
	Doassansiales	Bauer et al. 1997
	Georgefischeriales	Bauer et al. 1997
	Tilletiales	Bauer et al. 1997
Allied lineages
Class	Ustilaginomycetes	Bauer et al. 1997
Order	Violaceomycetales	Albu et al. 2015
Class	Exobasidiomycetes	Begerow et al. 2006
Order	Ceraceosorales	Begerow et al. 2006
	Exobasidiales	Hennings 1898; Bauer et al. 1997
	Microstromatales	Bauer et al. 1997
	Golubeviales	Wang et al. 2015
	Robbauerales	Wang et al. 2015
Class	Malasseziomycetes	Wang et al. 2014
Order	Malasseziales	Moore 1980
Class	Moniliellomycetes	Wang et al. 2014
Order	Moniliellales	Wang et al. 2014

The first representative of Ustilaginomycotina for which an entire genome was sequenced is U. maydis (Kämper et al. 2006). Its genome is one of the smallest known among fungal plant pathogens in part due to retention of relatively few repetitive elements. Being a pathogen, the genome has been utilized for the discovery of key effectors involved in infection and tumorigenesis (Skibbe et al. 2010; Döehlemann et al. 2014). It also serves as a resource for exploration of carbohydrate active enzymes (CAZymes) that interact with plant cell walls (Couturier et al. 2012), classified by their substrate specificity on cellulose, hemicellulose, lignin, and pectin (de Vries et al. 2017). While breakdown of plant cell walls can serve as an important food source for saprotrophs, it is also a critical component of fungal-plant signaling during plant infection (Nühse 2012; Spanu 2012). As the majority of Ustilaginomycotina are plant-associated, we suspect there is much potential for discovery of novel enzymes involved in plant biomass breakdown. Moreover, comparative analysis allows us to investigate whether molecular machineries associated with the smut pathology are conserved in other lineages. However, previous attempts to properly address these questions have been hampered for several reasons. Primarily, the lack of reference genomes for nearly all major Ustilaginomycotina lineages has severely limited inferences drawn from comparative studies at a subphylum-wide scale. In addition, poor phylogenetic resolution has been a major impediment to fully understanding this group of fungi in an evolutionary context (Begerow et al. 2014; Wang et al. 2015).

Here, we generated whole genome sequences for eight Ustilaginomycotina species—Ceraceosorus guamensis (Ceraceosorales), Jaminaea rosea and Pseudomicrostroma glucosiphilum (Microstromatales), Tilletiopsis washingtonensis (Entylomatales), Violaceomyces palustris (Violaceomycetales), Acaromyces ingoldii and Meira miltonrushii (Exobasidiales), and Testicularia cyperi (Anthracoideaceae, Ustilaginales)—representing families and orders for which no prior genomic data are available in any public database. Using these data, we resolved ordinal-level relationships among Ustilaginomycotina through phylogenomics. Additionally, we created a pipeline, called tree expansion analysis, for placing additional taxa without genome sequence data within our robustly resolved Ustilaginomycotina phylogeny. From this framework, we estimated divergence times and origins of teliospore formation (as a key characteristic of smut fungi). Finally, we compared several Ustilaginomycotina genomes to determine whether there are any shared properties that can be generalized for the smut pathogenic strategy, as well as explored functional gene space, particularly of plant biomass breakdown machinery.

Results

Whole Genome Sequencing and Genome Annotation

Eight Ustilaginomycotina genomes were sequenced, assembled, and annotated through the 1KFG project as described in the Materials and Methods section. All of the generated genomes appear to be haploid based on SNP consensus rate generated by ALLPATHS-LG (supplementary table S1, Supplementary Material online). Summaries of genome assembly and annotation statistics are shown in supplementary tables S1 and S2, Supplementary Material online.

Phylogenomic Reconstruction

For phylogenetic reconstruction, we included these eight newly sequenced Ustilaginomycotina genomes, as well as the publicly available genomes of seven Ustilaginomycotina species and six outgroups (two Pucciniomycotina and four Agaricomycotina). Single copy orthologs were identified using mcl clustering (Li et al. 2003), resulting in a total of 910 orthologous sets which were used for phylogenetic reconstruction (fig. 1). Most nodes in the tree were concordantly well-resolved with two different methods—concatenation and maximum likelihood (ML) and species tree estimation using average ranks of coalescences (STAR). The only existing conflict is in the placement of Malasseziomycetes. The result from the concatenation-ML places Malasseziomycetes as the earliest diverging lineage of the subphylum (fig. 1). In contrast, the STAR method suggests Malasseziomycetes as a sister class to Ustilaginomycetes, but the bootstrap support is lower than the concatenation-ML analysis (supplementary fig. S2, Supplementary Material online).

Fig. 1.

Phylogenomic tree reconstruction depicting ordinal-level relationships of Ustilaginomycotina. Twenty-one fungal genomes were included in the analyses: 15 species representing eight orders from Ustilaginomycotina, four species of Agaricomycotina and two species of Pucciniomycotina. The figure shows the phylogenomic tree topology resulting from concatenation of informative sites of 910 orthologous sets followed by ML reconstruction. Pucciniomycotina species (Mixia osmundae and Cystobasidium minutum) were selected as the outgroup. All nodes in the tree have a bootstrap value of 100. Thick lines indicate branches in the tree that are congruent with the topology from species tree estimation using the STAR method (supplementary fig. S2, Supplementary Material online). Branch lengths indicate the substitution differences between taxa from the ML tree reconstruction. Bar = 0.06 substitutions per amino acid residue.

Tree Expansion Analysis, Divergence Time Estimation, and Ancestral State Reconstruction of Teliospore Formation

We created a pipeline for placing additional taxa that lack genome sequence data (supplementary table S3, Supplementary Material online) within our robustly resolved Ustilaginomycotina phylogeny. For each phylogenetic reconstruction method, the expanded trees guided by two different phylogenomic trees (concatenation-ML and STAR) have very similar topologies except regarding the placement of Malasseziomycetes (supplementary figs. S3 and S5, Supplementary Material online for ML and supplementary figs. S4 and S6, Supplementary Material online for Bayesian tree reconstruction). The consensus tree after combining the expanded trees resulted from all analyses (supplementary figs. S3–S6, Supplementary Material online) is depicted in figure 2. Tilletiales could be precisely placed on the tree as sister to Microstromatales. The phylogenetic positions of Urocystidales (as sister to Ustilaginales) and Uleiellales (as sister to Violaceomycetales) are congruent with previous studies (Begerow et al. 2006; Riess et al. 2016). Moniliellales is most often recovered as the earliest-diverging lineage of Ustilaginomycotina. The exact placements of Doassansiales, Golubeviales, and Robbauerales could not be determined due to topological conflicts between ML and Bayesian methods. Despite this, the consensus tree indicates that Doassansiales and Robbauerales are orders in Exobasidiomycetes, and Golubeviales is clustered with Microstromatales and Tilletiales (fig. 2).

Fig. 2.

The consensus tree topology from tree expansion analyses. Forty additional taxa representing 15 orders and 4 classes in Ustilaginomycotina were included in the analyses (supplementary table S3, Supplementary Material online). Each rDNA sequence (SSU, ITS, LSU) was aligned and concatenated into a supermatrix. Then, maximum likelihood and Bayesian trees were reconstructed from the supermatrix using the phylogenomic tree topologies from figure 1 and supplementary figure S2, Supplementary Material online, as constraint backbone trees. Trees generated from all approaches were visually inspected to create the expanded consensus tree. Dashed lines refer to lineages that are not fully resolved due to topological conflicts between the two methods. Filled/unfilled circles in each node indicate the support values of ML/Bayesian approaches, respectively. Left and right halves of each circle indicate the support value from the expanded trees guided by figure 1 and supplementary figure S2, Supplementary Material online, topologies. The expanded tree topology from each analysis can be found in supplementary figures S3–S6, Supplementary Material online. Numbers and gray bars on the nodes represent the mean and 95% confidence interval range of the estimated divergence time in million years ago (Ma) unit. Details on estimated times can be found in supplementary table S4, Supplementary Material online. A triangle key mapped on each order indicates character states of teliospore formation. Character states of three deep nodes were estimated by the ancestral reconstruction analyses. Attached numbers are probabilities that a particular character state is present under the condition of Bayesian expanded tree topologies from supplementary figures S4 and S6, Supplementary Material online, respectively. An asterisk (*) indicates that the character state is statistically significant from the other. More information on the ancestral reconstruction analyses can be found in supplementary file 2, Supplementary Material online.

From the two Bayesian expanded trees guided by the two different phylogenomic trees, we estimated the divergence time of major classes and orders within the subphylum (fig. 2 and supplementary table S4, Supplementary Material online). Both expanded trees yield similar estimated times—the estimated times of some nodes are slightly different from Riess et al. (2016) with the deviation of circa ±50 My from the estimated means. These differences are likely due to differences in taxon sampling and gene loci used for analyses. According to our estimates, the crown node times of Ustilaginaceae and Tilletiales (excluding Erratomyces patelii which infects a leguminous host) are 44 and 43 Ma, respectively. These are after the divergence of major graminaceous (Poaceae) clades, which is circa 51–57 Ma (Bouchenak-Khelladi et al. 2010; Christin et al. 2014). The estimated time of the crown node of Malasseziales/Malasseziomycetes is close to the estimated time of the crown node of placental mammals, which is circa 90–100 Ma (Meredith et al. 2011; Springer et al. 2017). Seven of the 15 orders in Ustilaginomycotina are comprised of teliosporic fungi. According to the ancestral state reconstruction analyses, teliospore formation is present at the most recent common ancestor (MRCA) node of Exobasidiomycetes, Ustilaginomycetes, and Malasseziomycetes (fig. 2 and supplementary file 2, Supplementary Material online).

Genome Architecture Comparison

Compared with other subphyla of Basidiomycota, all Ustilaginomycotina genomes are relatively small (the median genome size of ca. 18 Mb) with repetitive elements <5%. It is notable that the genomes of Malassezia are much smaller than other Ustilaginomycotina species, with the genome sizes of 8.9 and 7.67 Mb for M. globosa and M. sympodialis (Xu et al. 2007; Gioti et al. 2013), respectively (supplementary fig. S7, Supplementary Material online). The genes of all Ustilaginomycotina genomes are densely packed together, with the median gene density of 389 genes per Mb (supplementary fig. S7, Supplementary Material online). This is in striking contrast to Pucciniomycotina and Agaricomycotina, where genome sizes can range from 10 Mb up to > 100 Mb with varied proportion of repetitive elements. The median gene densities of Pucciniomycotina and Agaricomycotina are circa 319 and 278 genes per Mb, respectively. In addition, while pathogenicity is often associated with genome expansion in other biotrophic fungi, the genomes of Ustilaginomycotina pathogens are relatively small and compact (fig. 3).

Fig. 3.

Comparison of genome architecture among fungal and fungal-like plant pathogens. The tree illustrates phylogenetic relationships between pathogenic species in five different clades: Oomycetes, Pucciniomycotina, Ustilaginomycotina, Agaricomycotina, and Ascomycota. The figure was partially modified from Raffaele and Kamoun (2012). Pathogenic strategy, genome size and gene density of each species were obtained from the literature (see Materials and Methods) or the JGI MycoCosm portal (Grigoriev et al. 2014). ND, no available data.

Comparative Studies on Genes Related to Pathogenic Strategies

The majority of genes that have been implicated in smut infection and pathogenesis in the model U. maydis (Skibbe et al. 2010; Döehlemann et al. 2014) are limited in distribution to Ustilaginaceae (table 2). For example, Testicularia cyperi, the smut fungus most closely related to U. maydis in our analysis, possesses only five out of 35 genes. Two of these—Pit1, encoding a membrane protein required for tumor formation (Döehlemann et al. 2011), and Clp1, encoding a regulator of cell cycle control (Heimel et al. 2010)—are also detected in V. palustris (supplementary file 3, Supplementary Material online), a recently described yeast species isolated from aquatic ferns (Albu et al. 2015). Only a gene encoding the high affinity sucrose transporter Srt1 (Wahl et al. 2010), which is required for full virulence of U. maydis, is present across most of Ustilaginomycotina species, except species of Microstromatales, Malasseziales, and Violaceomycetales. Details can be found in supplementary file 3, Supplementary Material online.

Table 2.

Orthologous Genes Related to Smut Infection and Pathogenicity.

Order/Family	Species	Validated Organ-Specific Genes During Ustilago maydis Infection and Tumorigenesis (Skibbe et al. 2010)	Published Smut Pathogenicity Genes (Döehlemann et al. 2014)
Ustilaginales
Ustilaginaceae	Ustilago maydis	24	11
	Pseudozyma hubeiensis	14	11
	Sporisorium reilianium	15	10
Anthracoideaceae	Testicularia cyperi	1	4
Malasseziales	Malassezia globosa	0	0
	M. sympodialis	0	0
Exobasidiales	Acaromyces ingoldii	0	1
	Exobasidium vaccinii	0	1
	Meira miltonrushii	0	1
Entylomatales	Tilletiopsis washingtonensis	0	1
Ceraceosorales	Ceraceosorus guamensis	0	1
Microstromatales	Jaminaea rosea	0	0
	Pseudomicrostroma glucosiphilum	0	1
Georgefischeriales	Tilletiaria anomala	0	1
Violaceomycetales	Violaceomyces palustris	0	2

CAZyme Analysis

To explore plant biomass breakdown potential, we assessed whether any CAZymes differ significantly (Bonferroni corrected p ≤ 0.05; Fisher’s exact test) in abundance in Ustilaginomycotina compared with other basidiomycete subphyla. Our results indicate that Ustilaginomycotina members lack many of the core enzymes for breakdown of lignin and cellulose (fig. 4 and supplementary table S5 and supplementary file 4, Supplementary Material online), for example, AA1_1 and AA2 (lignin), AA9 (cellulose) and CBM1 (cellulose binding). They also lack several enzymes for starch/glycogen metabolism such as CBM20 (starch binding), GH13_32 (secreted α-amylase), GT35 (glycogen phosphorylase). In particular, there is a complete loss of the enzymes that target glycogen (GH13, GH15, GH133, GT3, GT35, CBM20, and CBM21) in Malassezia species. However, Ustilaginomycotina harbors several unique CAZyme families that are absent from other Basidiomycota such as GH5_16 (β-1,6-galactanase), GH8 (broad activity hydrolase), GH42 (β-galactosidase), GT34 (α-galactosyltranferase), and AA10 (lytic polysaccharide monooxygenase). Based on substrate specificity, GH5_16 and GH42 are potentially involved in pectin decomposition.

Fig. 4.

Results of Fisher’s exact test for CAZyme abundance in Ustilaginomycotina compared with other subphyla of Basidiomycota. Tracks show gene counts within an individual genome, shaded by number of gene copies. Information of CAZyme functions and relative abundances are listed in supplementary table S5, Supplementary Material online. Genome order and underlying gene counts are provided in supplementary file 4, Supplementary Material online. Orange: Ustilaginomycotina, green: Agaricomycotina, blue: Pucciniomycotina. Labels for Plant Cell Wall Degrading Enzymes (PCWDEs) are colored purple. Only CAZyme families that show a significant difference (Bonferroni corrected p ≤ 0.05) in abundance between Ustilaginomycotina and other Basidiomycota are shown.

On average, Ustilaginomycotina has circa 16 CAZyme families considered as plant cell wall-decomposing enzymes (PCWDEs), except Malassezia species that lack all PCWDEs. Most Ustilaginomycotina species, which are primarily biotrophic or plant-associated fungi, have fewer PCWDEs (ca. 20 families on an average) compared with necrotrophic and saprotrophic fungi (ca. 100 families on an average), while at the same time have more PCWDEs than fungal species associated with animals (ca. four families on an average). Details on CAZymes and PCWDEs in each fungal group are tabulated in supplementary file 4, Supplementary Material online.

Gene Conservation

On average, each species harbors circa 532 subphylum-specific genes (fig. 5), with the highest number in V. palustris (782 genes) and the lowest number in M. sympodialis (144 genes). Species of Malasseziomycetes share the fewest number of genes with any other species included in the analyses. About 9, 13, and 503 shared orthologous sets are found between Malasseziomycetes/Ustilaginomycetes, Malasseziomycetes/Exobasidiomycetes and Ustilaginomycetes/Exobasidiomycetes, respectively (supplementary fig. S8, Supplementary Material online).

Fig. 5.

Gene conservation across Ustilaginomycotina genomes. Protein models of 15 Ustilaginomycotina genomes were run through MCL clustering to determine orthologous clusters. Each cluster was then evaluated for level of conservation using taxonomic assignments of each species according to NCBI taxonomy database. Recently described species (since 2010) are indicated by an asterisk (*).

For each species included in the analyses, there are about 948 species-specific genes on average. Half of the species-specific genes (6,609 of 14,217) are derived from only five species that were recently discovered and described within the past 5 years—M. miltonrushii (Rush and Aime 2013), V. palustris (Albu et al. 2015), C. guamensis (Kijpornyongpan and Aime 2016), P. glucosiphilum, and J. rosea (Kijpornyongpan and Aime 2017). Of the species-specific genes, 2,359 could be associated with previously known KOG genes with Tilletiopsis washingtonensis, C. guamensis, and M. miltonrushii being particularly enriched in these (309, 307, and 299, respectively). The remaining 83.4% of the species-specific genes (11,857/14,217) are not found to be homologous with any genes in the KOG database. Ceraceosorus guamensis and V. palustris, which both represent newly described species in rare lineages of Ustilaginomycotina (Albu et al. 2015; Kijpornyongpan and Aime 2016), possess the highest numbers of non-KOG species-specific genes at 1,415 and 1,461, respectively. Details on the species-specific genes are enumerated in supplementary file 5, Supplementary Material online.

Discussion

In this study, we sequenced eight new genomes representing orders and families of Ustilaginomycotina for which no prior genomic data were available. Together with data available in public databases, we reconstructed a robust phylogeny of Ustilaginomycotina. We then employed a tree expansion analysis by using the well-resolved phylogenomic tree as a guiding backbone to place additional taxa on the tree. Through these approaches, we are able to establish the relationships of most orders in Exobasidiomycetes (figs. 1 and 2), a class which has completely evaded previous attempts at resolution with single or multilocus approaches (Begerow et al. 2006; Wang et al. 2015). This demonstrates the utility of combining genomic data with the tree expansion approach utilized here for fungal systematics. As the availability of fungal genomes increases through efforts such as the 1KFG project (Grigoriev et al. 2011), this approach may be effective for phylogenetic placement of many other taxa for which single gene loci are available, such as those generated during the Assembling the Fungal Tree Of Life (AFTOL) project (Spatafora 2005).

A compact genome is a characteristic shared across all Ustilaginomycotina species sampled in this study (fig. 3 and supplementary fig. S7, Supplementary Material online). The median genome size of Ustilaginomycotina (ca. 18 Mb) is much smaller than that of most Fungi (which have a median genome size of ca. 36 Mb; the JGI fungal genome portal MycoCosm). While significant genome expansion events are observed across most other subphyla of Basidiomycota, this pattern is not, thus far, observed in Ustilaginomycotina. This reflects the general trend of a yeast-like life strategy—10 out of 15 orders of the subphylum contain yeast-like or dimorphic fungi—although a few filamentous fungi exist in this subphylum (Begerow et al. 2014; Wang et al. 2014; Albu et al. 2015). For example, small genomes also appear to be a characteristic of other yeasts within Taphrinomycotina and Saccharomycotina in Ascomycota, Microbotryomycetes, Mixiomycetes, and Cystobasidiomycetes in Pucciniomycotina and Tremellomycetes in Agaricomycotina (Dujon 2010; Nagy et al. 2014; Toome et al. 2014b).

More than 40% of described true smut fungi belong to Tilletiales and Ustilaginaceae, which are specialized on members of Poaceae (Vánky 2013; Begerow et al. 2014). Based on our time estimation results, we found that these two lineages diverged circa 10 My after the major divergence of Poaceae (fig. 2). Since Poaceae is one of the most successful angiosperm lineages in terms of dispersal and diversification (Bouchenak-Khelladi et al. 2010), specialization on this lineage can accelerate pathogen diversification according to Fahrenholz’s rule of host-parasite evolution (Fahrenholz 1913). This was previously illustrated for at least a few smut genera through a phylogenetic framework (Begerow et al. 2004; Escudero 2015).

The ancestral state reconstruction supports our hypothesis that teliospore formation—a key characteristic of true smut fungi—is synapomorphic in Ustilaginomycotina (fig. 2). This is in direct contrast to the hypothesis of Riess et al. (2016), which suggests that early diverging lineages of Ustilaginomycotina are saprotrophic. Subsequent parallel losses of teliospore formation in Microstromatales, Ceraceosorales, and Exobasidiales, where these lineages have transitioned to directly producing basidiospores, are notable (fig. 2). Pathogenic members of these orders, in contrast to the majority of smut fungi that infect herbaceous hosts, are specific to woody perennial plants (Begerow et al. 2014). Teliospores are also present in other fungal groups, such as in rust fungi (Pucciniomycetes, Pucciniomycotina), anther smut fungi and relatives (Microbotryomycetes, Pucciniomycotina) and root gall-forming fungi (Entorrhizomycetes, Entorrhizomycota) (Aime et al. 2014; Bauer et al. 2015)—some members of which can infect both woody and herbaceous plants. Teliospores in these non-Ustilaginomycotina fungi are functionally equivalent to those produced in the smut fungi as they give rise to basidia, which subsequently undergo meiosis to produce haploid basidiospores. However, it is still unclear whether teliospores in these groups of fungi are truly homologous or whether this trait independently arose several times across the fungal tree of life. More robustly resolved phylogenetic trees in other lineages are required to fully address this question.

We found that genes related to smut infection and pathogenesis identified in U. maydis are not widespread throughout Ustilaginomycotina. Even though T. cyperi belongs to the same order (but different family; Anthracoideaceae) as U. maydis (Vánky 2013), it shares much fewer orthologs than the species of Ustilaginaceae, U. maydis, S. reilianum, and P. hubeiensis (table 2). It should be noted that these species differ in thier specific hosts—U. maydis and S. reilianum commonly parasitize maize (Poaceae), while T. cyperi infects sedges (Cyperaceae). Although P. hubeiensis has thus far not been found as a plant pathogen, gene conservation indicates that it was ancestrally pathogenic on grass species. Meanwhile, the studied species of Exobasidiomycetes, such as the gall-forming phytopathogen E. vaccinii, lack almost all of these genes. While genomic data are yet to be available for members of Tilletiales, these would be important for comparative studies as their host plants are graminaceous, similar to the smut fungi in Ustilaginaceae. Considering the placement in Exobasidiomycetes (fig. 2), we hypothesize that Tilletiales lacks this set of orthologous genes and likely has developed an alternative mechanism for infection and pathogenesis on grasses.

Plant cell wall decomposing enzymes play an important role in fungal nutritional strategies. In general, our results are consistent with the previously observed trend that biotrophic pathogens tend to have fewer CAZymes compared with necrotrophic or saprotrophic fungi, and animal-associated fungi tend to have fewer CAZymes than plant-associated fungi (Zhao et al. 2013). Through comparison across different subphyla in Basidiomycota, we found that several core enzymes involved in breakdown of lignin and cellulose are absent from the entire Ustilaginomycotina. Moreover, having a limited number of CAZymes for starch catabolism indicates that Ustilaginomycotina members acquire small carbohydrates primarily from direct absorption or bioconversion, but not through decomposition of storage polysaccharide. Finally, some smut fungi uniquely harbor the CAZyme families GH5_16 and GH42 involved in pectin degradation. We speculate that these enzymes may assist in host tissue maceration in order to enlarge intercellular space during teliosporogenesis (Piepenbring et al. 1998).

Malasseziomycetes and Moniliellomycetes represent two recently elevated Ustilaginomycotina classes with equivocal phylogenetic placements (Wang et al. 2014). Our inability to confidently place these lineages suggests their extreme divergence from other members of Ustilaginomycotina, which is consistent with the observation that members of these classes are not plant-associated fungi. Members of Malasseziomycetes are lipophilic animal-associated fungi commonly found on mammal skin (Guého et al. 1998, 2011), while Moniliellomycetes comprises both animal-associated and saprotrophic fungi isolated from industrial settings (de Hoog et al. 2011). Genomes of Malassezia species are also markedly different from those in other classes of Ustilaginomycotina. For instance, the genome sizes of Malassezia spp. are half that of other studied species (supplementary fig. S7A, Supplementary Material online). Moreover, disregarding common genes in the subphylum, the number of shared orthologous genes between Malasseziomycetes and other classes account for only 4.4% of those shared between Ustilaginomycetes and Exobasidiomycetes (supplementary fig. S8, Supplementary Material online). The complete loss of CAZymes both for plant cell wall decomposition and starch/glycogen metabolism in Malassezia species correlates with their divergent nutritional modes. The overlapping crown date estimations for Malasseziomycetes and placental mammals (fig. 2), also reflects the adaptation of Malasseziomycetes to become mammal-associated.

It is notable that around half of the Ustilaginomycotina species-specific genes, most of which are yet to be clearly annotated, are found in five recently described species—M. miltonrushii, V. palustris, C. guamensis, P. glucosiphilum, and J. rosea. Moreover, C. guamensis and V. palustris, which represent underexplored lineages of the subphylum (Albu et al. 2015; Kijpornyongpan and Aime 2016), possess the largest number of species-specific genes (1,722 and 1,697 genes, respectively; fig. 5). This highlights the continued importance of discovering new fungal species as a huge resource of hidden genetic diversity. As most Ustilaginomycotina members are plant-associated fungi, further comparative genomics of this overlooked fungal subphylum may lead to the discovery of novel genes that have a potential for industrial and agricultural applications in the future.

Materials and Methods

Fungal Strains Culturing and Nucleic Acid Extraction

A total of eight fungal strains, representing six orders of Ustilaginomycotina, were selected for genome sequencing: Testicularia cyperi ATCC MYA-4640 for Ustilaginales, Meira miltonrushii CBS12591, and Acaromyces ingoldii CBS140884 for Exobasidiales, Tilletiopsis washingtonensis NRRL Y-63783 for Entylomatales, Ceraceosorus guamensis CBS139631 for Ceraceosorales, Pseudomicrostroma glucosiphilum CBS14053, and Jaminaea rosea CBS14051 for Microstromatales and Violaceomyces palustris CBS139708 for Violaceomycetales. An additional seven Ustilaginomycotina strains for which genomic data are publicly available in the JGI MycoCosm portal (http://genome.jgi.doe.gov/programs/fungi/index.jsf; last accessed April 20, 2018) were included: Exobasidium vaccinii (Döehlemann G, et al., unpublished data) for Exobasidiales, Malassezia globosa CBS7966 (Xu et al. 2007) and M. sympodialis ATCC 42132 (Gioti et al. 2013) for Malasseziales, Pseudozyma hubeiensis SY62 (Konishi et al. 2013), Sporisorium reilianum SRZ2 (Schirawski et al. 2010), and Ustilago maydis 521 (Kämper et al. 2006) for Ustilaginales and Tilletiaria anomala CBS436.72 (Toome et al. 2014a) for Georgefischeriales. Two strains of Pucciniomycotina, Mixia osmundae IAM 14324 (Toome et al. 2014b), and Cystobasidium minutum NRRL Y-63784 (Toome M and Aime MC, unpublished data) and four strains in Agaricomycotina, Tremella mesenterica ATCC24925 (Floudas et al. 2012), Cryptococcus neoformans var. grubii H99 (Janbon et al. 2014), Calocera viscosa (Nagy et al. 2016), and Wallemia sebi CBS633.66 (Padamsee et al. 2012), representing early diverging lineages of these two subphyla, were also included in the analyses for rooting purposes.

To extract nucleic acids from the eight strains selected for sequencing, each was cultured on potato dextrose agar (PDA) or in potato dextrose broth (PDB) for 4–10 days. Fungal tissues were then harvested by scraping the tissues from the agar surface for the strains on PDA, or centrifuging and rinsing with sterile water for the strains in PDB. After that, the fungal tissues were ground in liquid nitrogen, and stored at −80°C. DNA of each strain was extracted from ground and frozen tissues using the Promega Wizard Genomic Purification kit (Promega, Madison, WI) and treated with RNase A solution. Individual extraction of each strain was repeated until a total of circa 20 µg of DNA was acquired, then pooled into a single microfuge tube. DNA quantity was remeasured by Quantifluor dsDNA dye (Promega). RNA isolation and on-column DNA digestion were processed using E.Z.N.A. fungal RNA kit and RNase-free DNase set (Omega Bio-Tek, Norcross, GA). The RNA extract was quantified using a Nanodrop spectrophotometer.

Sequencing, Assembly, and Annotation Methods

Genomes in this study were sequenced using an Illumina platform. For genomes two libraries have been built: fragment (270 bp insert size) and 4 kb the long-mate-pair (LMP) libraries.

Fragment libraries were obtained from 100 ng of DNA, sheared to 270 bp using the Covaris LE220 (Covaris, Woburn, MA) and size selected using SPRI beads (Beckman Coulter, Brea, CA). The fragments were treated with end-repair, A-tailing, and ligation of Illumina compatible adapters (IDT Inc., Coralville, IA) using the KAPA-Illumina library creation kit (KAPA biosystems, Wilmington, MA).

For LMP libraries, 5 µg of DNA was sheared using the Covaris g-TUBE and gel size selected for 4 kb. The sheared DNA was treated with end repair and ligated with biotinylated adapters containing loxP. The adapter ligated DNA fragments were circularized via recombination by a Cre excision reaction (NEB, Ipswich, MA) and randomly sheared using the Covaris LE220. The sheared fragments were treated with end repair and A-tailing using the KAPA-Illumina library creation kit (KAPA biosystems) followed by immobilization of mate pair fragments on streptavidin beads (Invitrogen, Waltham, MA). Illumina compatible adapters (IDT Inc.) were ligated to the mate pair fragments and 8–10 cycles of PCR was used to enrich for the final library (KAPA Biosystems).

For transcriptomes, stranded cDNA libraries were generated using the Illumina Truseq Stranded RNA LT kit. mRNA was purified from 1 µg of total RNA using magnetic beads containing poly-T oligos, fragmented using divalent cations and high temperature, and then reverse transcribed using random hexamers and SSII (Invitrogen) followed by second strand synthesis. The fragmented cDNA was treated with end-pair, A-tailing, adapter ligation, and 10 cycles of PCR.

The prepared libraries were quantified using KAPA Biosystem’s next-generation sequencing library qPCR kit and run on a Roche LightCycler 480 real-time PCR instrument. The quantified libraries were then prepared for sequencing on the Illumina HiSeq sequencing platform utilizing a TruSeq paired-end cluster kit (v3 for genomes and v4 for transcriptomes) and Illumina’s cBot instrument to generate a clustered flowcell for sequencing. Sequencing of the flowcell was performed on the Illumina HiSeq2000 (HiSeq2500 for transcriptomes) sequencer using a TruSeq SBS sequencing kit 200 cycles following a 2×150 (for genomic fragments and transcriptome) or 2×100 (for LMP) indexed run recipes.

Each fastq file was QC filtered for artifact/process contamination and pairs of fragment and LMP Illumina data sets were subsequently assembled together with AllPathsLG (Gnerre et al. 2011). For C. guamensis MCA4658 and V. palustris SA807 with no LMP produced, the initial assembly of fragment data was performed with Velvet (Zerbino and Birney 2008) and used to create a long mate-pair library with insert 3000±300 bp in silico, which was then assembled together with the original Illumina library using AllPathsLG. Transcriptome reads were de novo assembled using Rnnotator v. 3.4.0 (Martin et al. 2010) and mapped to genome assemblies to assess genome completeness and to facilitate genome annotation. Genomes were annotated using the JGI Annotation pipeline and made available via the JGI fungal genome portal MycoCosm (Grigoriev et al. 2014; jgi.doe.gov/fungi; last accessed April 20, 2018). The genome and gene annotation data of eight Ustilaginomycotina species are deposited at DDBJ/EMBL/GenBank under the following accessions—MCHA00000000 for P. glucosiphilum, MCHB00000000 for J. rosea, MCHC00000000 for M. miltonrushii, MCHD00000000 for C. guamensis, MCHE00000000 for V. palustris, MCOH00000000 for Testicularia cyperi, MCOI00000000 for A. ingoldii, and MKCN00000000 for Tilletiopsis washingtonensis.

Single-Copy Orthology Finding, Sequence Alignment, and Phylogenomic Reconstruction

The filtered protein model of each studied taxon was downloaded from the JGI MycoCosm portal. The orthologous sets among studied taxa were evaluated using MCL clustering with inflation factor 2 (Enright et al. 2002). The total 916 single-copy orthologous sets present in all 21 species were further processed. An alignment of each orthologous set was performed in MAFFT version 6.903 (Katoh and Toh 2008) using mafft –auto mode with default parameters. Then, the informative sites for phylogenetic purposes were extracted from the alignment of each orthologous set using GBLOCKs version 0.91 b (Talavera and Castresana 2007). About 6 of 916 orthologous sets were screened out due to lack of informative sites.

Two different methods were designed for the species tree reconstruction: concatenation and maximum-likelihood (ML) and species tree estimation from gene trees using pseudolikelihood method. For concatenation-ML analyses, the informative sites of 910 orthologous sets were concatenated into a supermatrix. The ML tree was reconstructed from the supermatrix by RAxML version 7.3.0 (Stamatakis 2006) using –f a mode with 1000 bootstrap replicates. PROTGAMMAWAGF was used as a protein substitution model, and other parameters were set as default. Mixia osmundae and C. minutum were selected as the outgroup taxa for these two methods based on the recovered topologies of Toome et al. (2014b) and Padamsee et al. (2012). For species tree estimation from gene trees, the optimal protein substitution model of the informative sites of each orthologous set was evaluated using Prottest version 3.3 (Darriba et al. 2011). Each ML gene tree construction was then performed using RAxML in –f a mode with 1000 bootstrap replicates, optimized protein model and M. osmundae and C. minutum selected as the outgroup. After that, all generated gene trees were pooled for species tree estimation using average ranks of coalescences (STAR) method with bootstrap calculation (Liu et al. 2009) performed through STRAW server (Shaw et al. 2013). An overview of phylogenomic pipeline can be found in supplementary figure S9, Supplementary Material online.

Tree Expansion Analyses, Divergence Time Estimation, and Ancestral State Reconstruction

The phylogenomic tree topologies from the previous step were used as a guide for the tree expansion analysis, which included more taxa from every known order of Ustilaginomycotina to place representatives of orders not yet present on the phylogenomic tree: Moniliellales, Urocystidales, Tilletiales, Doassansiales, Golubeviales, Robbauerales, and Uleiellales. Ribosomal DNA (rDNA) sequences of the small subunit, internal transcribed spacer region, and D1/D2 region of large subunit were used in the analyses (supplementary table S3, Supplementary Material online). First, each region of rDNA sequences from 21 species present in the phylogenomic tree was aligned in MAFFT version 6.903 using the phylogenomic tree topology from figure 1 as a guide tree for the alignment. Then, the expanded taxa were included into the alignment of each rDNA sequence using mafft –add function from MAFFT online version 7 (Katoh and Standley 2013; http://mafft.cbrc.jp/alignment/server/add.html; last accessed April 20, 2018). The alignment of all rDNA sequences was manually curated and concatenated into one supermatrix.

The supermatrix was used to perform ML tree reconstruction in PAUP* version 4.0a151 (Swofford 2002). Initial analyses were performed using heuristic search with keeping the optimal trees and enforcing topological constraints using the tree topologies from figure 1 and supplementary figure S2, Supplementary Material online, as backbone constraint trees. Other parameters were set as default. Then, a set of optimal trees from the initial analyses was used as starting trees for the subsequent heuristic search with keeping 500 best trees and enforcing the topological constraint. The 50% majority-rule consensus tree was finally computed from the 500 best trees, and each node was numbered based on the frequency of occurrence.

The supermatrix was also used for Bayesian phylogenetic reconstruction of the expanded tree using BEAST package v. 1.7.5 (Drummond and Rambaut 2007; Drummond et al. 2012). The parameters for BEAST were set as following—general time reversible model with substitution-rate among sites of gamma distribution (GTR + G) was set for nucleotide substitution model, the tree topologies from figure 1 and supplementary figure S2, Supplementary Material online, were applied as a backbone topological constraint, Yule speciation process was modeled as a tree prior (Gernhard 2008), the Markov chain Monte Carlo (MCMC) was run for ten million generations, a tree was sampled every 1,000 generations, and other parameters were set as default. The 10,000 sampled trees were annotated for the maximum clade credibility tree using TreeAnnotator by discarding the first 1,000 sampled trees as burn-in, and the posterior probability of each node was annotated. Each phylogenetic reconstruction method was performed for three independent runs to ratify consistent results. Finally, the tree topologies resulted from two methods were combined to visually inspect for the most consensus species tree. The overview of tree expansion pipeline is provided in supplementary figure S10, Supplementary Material online.

To perform divergence time estimation, we deployed the estimated divergence times from Floudas et al. (2012) and used for secondary time calibration. The divergence times between Tremella mesenterica and Cryptococcus neoformans (mean = 152.8 Ma, SD = 37.5) and Ustilago maydis and Malassezia globosa (mean = 272.88 Ma, SD = 53) were selected as nodes for calibration. Since there are two possible diverging nodes between U. maydis and M. globosa, we ran time estimation analyses for both scenarios, which provide the similar results (supplementary table S4, Supplementary Material online). To set up the BEAST .xml files, we employed uncorrelated relaxed clock model with lognormal relaxed distribution. The nodes to be estimated for ages (supplementary table S4, Supplementary Material online) were constrained to be monophyletic in accordance with the consensus Bayesian tree resulted from the previous steps (supplementary figs. S4 and S6, Supplementary Material online). The priors of calibrated nodes were set as normal distribution with means and SDs as mentioned earlier. The MCMC chain was run for 50 million generations with sampling every 1,000 generations. Estimated divergence times were obtained from the .log file after discarding the first 10% of states as burn-in using Tracer v.1.6. The log files resulted from two scenarios were combined to find the final estimated time. We performed the whole analyses in three independent runs to ensure the consistency of the results.

The character states of teliospore formation of orders presented in the expanded tree were obtained from the recent literature (Begerow et al. 2014), otherwise indicated as unknown. The ancestral state reconstruction was performed in BayesTraits V2.0 (Pagel and Meade 2014). The post burn-in 9,000 trees from Bayesian phylogenetic reconstruction (mentioned earlier) were used as the set of input trees for the analysis. The initial ML analysis of a teliospore formation trait was run using default parameters to determine the distribution of transition rates between two character states (absent or present). After that, the same dataset was used for the MCMC analysis. Three nodes on the expanded tree were added in order to estimate the probability that each character state occurs: the most recent common ancestor (MRCA) node of Exobasidiomycetes, the MRCA node of Ustilaginomycetes, and the MRCA node of Exobasidiomycetes, Ustilaginomycetes, and Malasseziomycetes. The analyses were run as following—the chain was run for 1,010,000 generations with the first 10,000 generations discarded as burn-in, the estimated probability values were sampled every 1,000 generations, the priors of transition rates between two character states were set in accordance with the results from the ML analysis, and other parameters were set as default.

To test which character state is statistically significant to occur on each of the three MRCA nodes, we fixed a certain character state to a node of interest and ran the MCMC analysis using the same parameters as above. The log marginal likelihood was estimated using the stepping stone sampler method with the settings of using 100 stones and running 10,000 generations for each stone. Finally, we calculated log Bayes factor based on difference between log marginal likelihoods of each fixed character state on the certain node. A value >2 was interpreted as statistically significant according to interpretation of Gilks et al. (1995).

Comparative Studies on Genome Architecture and Genes for Pathogenic Strategies

Genome assembly size, proportion of repetitive elements in genome and gene density of fungal species were obtained from previous literature (Kämper et al. 2006; Cuomo et al. 2007; Hane et al. 2007; Haas et al. 2009; Baxter et al. 2010; Lévesque et al. 2010; Schirawski et al. 2010; Spanu et al. 2010; Amselem et al. 2011; Cantu et al. 2011; Duplessis et al. 2011; Klosterman et al. 2011; Links et al. 2011; Collins et al. 2013; Wibberg et al. 2013; Toome et al. 2014a), otherwise obtained from JGI MycoCosm page (Grigoriev et al. 2014).

Lists of validated organ-specific genes related to U. maydis infection and tumorigenesis, as well as published smut pathogenicity genes were obtained from Skibbe et al. (2010) and Döehlemann et al. (2014), respectively. These sets of genes were used as reference genes for orthology finding across the 15 Ustilaginomycotina species via best-reciprocal blast hit (BBH) method using a customized Perlscript (available upon a request).

The detection, module composition, and family assignment of all CAZymes were performed as described previously (Cantarel et al. 2009; Levasseur et al. 2013; Lombard et al. 2014). Briefly, the method combined BLAST and HMMER searches conducted against sequence libraries and HMM profiles made of the individual functional modules featured in the CAZy database (http://www.cazy.org; last accessed April 20, 2018). Proteins that shared >50% identity over the entire domain length of an entry in the CAZy database were directly assigned to the same family. Proteins with <50% identity to a protein in CAZy were manually inspected (searching for catalytic residues when possible). All positive hits were manually examined for final validation. The PCWDE families were designated based on currently known CAZyme activities and substrates found as a part of the plant cell wall. CAZymes involved in starch/glycogen metabolism were also designated in a similar manner. Published genomes included in the analyses were listed in supplementary file 4, Supplementary Material online. Following family assignment, we identified CAZyme families that differed significantly (Bonferroni corrected p ≤ 0.05) in abundance in Ustilaginomycotina compared with other basidiomycetes using Fisher’s exact test (two-tailed).

Gene Conservation Analysis

To assess gene conservation across lineages in Ustilaginomycotina (from the 15 studied species), the orthologous cluster generated via MCL described earlier was used. Each species was then assigned a taxonomic ID in each rank according to taxonomic assignments from the NCBI taxonomy database (http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi; last accessed April 20, 2018). After that, each cluster from the MCL result was investigated to identify the level of gene conservation (i.e., whether the genes within a given cluster are only found across all members of the same genus, family, order, class, or subphylum). Finally, a number of unique genes and conserved genes in each species were enumerated. Class-specific gene conservation shared by all species of three different classes (Ustilaginomycetes, Exobasidiomycetes, and Malasseziomycetes) was enumerated and presented as a Venn diagram (supplementary fig. S8, Supplementary Material online).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Acknowledgments

We would like to thank Dr Gunther Döehlemann from University of Cologne, Germany for permission to utilize unpublished genomic data of Exobasidium vaccinii. We acknowledge Dr Merje Toome for use of unpublished genomic data of Cystobasidium minutum, Agaricostilbum hyphaenes, Heterogastridium pycnidioideum, and Tritirachium sp. The computing facilities for bioinformatic analyses were granted by Dr Michael Gribskov from Department of Biological Science and Dr Jyothi Thimmapuram from Bioinformatics Core, Purdue University. We thank John F. Klimek and the Aime lab members for molecular support and useful comments on earlier versions of this manuscript. The work conducted by the U.S. Department of Energy Joint Genome Institute, a DOE Office of Science User Facility, was supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Anandamahidol scholarship from Thailand provided a financial support for T.K.’s study at Purdue University; USDA Hatch project 1010662 to M.C.A. provided support as did NSF DEB-1458290 to M.C.A. The research was conducted at Purdue University, US Department of Energy Joint Genome Institute, and Aix-Marseille Université.

Author notes

Teeratas Kijpornyongpan and Stephen J. Mondo contributed equally to this work.