All fungal genomes were downloaded from the NCBI Genome Database, with the exception of the O. ulmi W9 genome (downloaded from http://www.moseslab.csb.utoronto.ca/o.ulmi/). The descriptions and accession numbers are displayed in Supplementary Table 1 The predicted proteins/protein annotation for O. novo-ulmi H327 is available as supporting information from Comeau et al. (2015), and their proposed nomenclature was used throughout the text. To assess the completeness and quality of genome sequence assembly/predicted proteins, BUSCO version 2.0 (Simao et al., 2015) was used with the Sordariomyceta dataset and with species parameters set to the nearest represented species (Neurospora crassa). This quality assessment was applied to all genomes from the Ophiostomataceae family used during this work and its results are displayed in Supplementary Table 1. Putative BGCs in the O. ulmi W9 and O. novo-ulmi H327 genomes were identified with antiSMASH 3.0 (using the genome sequence assembly as input) (Weber et al., 2015) and SMIPS 2016-07-26 (using the predicted proteins as input) (Wolf et al., 2015) algorithms. Previously reported results were also examined (Comeau et al., 2015). Notably, BGCs and their fragments can remain in the genome without generating functional products (Campbell et al., 2012). Some genes of these apparently non-functional BGCs may be subject to pseudogenization, driving the synthesis of truncated products (Campbell et al., 2012). To validate the backbone gene of these predicted BGCs and obtain insight about possible non-functional BGCs, we analyzed RNA-seq data from O. novo-ulmi H327 (NCBI’s BioProject access: PRJNA325932) (Nigg et al., 2015). Briefly, in the experiment conducted by Nigg et al. (2015), O. novo-ulmi H327 was grown in three conditions and both fungal lifestyle structures were analyzed: (1) yeast in liquid culture medium; (2) mycelium in static liquid culture medium; and (3) mycelium in solid culture medium. The detailed RNA-seq experimental procedure, sequencing, data management, and statistics have been previously described (Nigg et al., 2015). Expression data for each gene and differentially expressed genes were extracted from supplementary material (http://www.g3journal.org/content/5/11/2487.supplemental; Supplementary Tables S4–S6). Furthermore, the borders of each cluster that were initially delimited based on antiSMASH 3.0 prediction were subsequently confirmed using CASSIS 2016-10-20 (Wolf et al., 2015) or compared with conserved orthologous gene clusters of other filamentous fungi species using MultiGeneBlast 1.1.14 (Medema et al., 2013).

The conservation of predicted BGCs among Ophiostomataceae species (Ophiostoma piceae strain UAMH 11346, Sporothrix schenckii strain 1099-18, Sporothrix brasiliensis strain 5110, Sporothrix pallida strain SPA8, Grosmannia clavigera strain kw1407, Graphilbum fragrans strain CBS 138720, Leptographium lundbergii strain CBS 138716, and Leptographium procerum strain CMW34542) was assessed using MultiGeneBlast 1.1.14 (analyses of genomic loci with predicted proteins) (Medema et al., 2013), based primarily on backbone gene conservation (amino acid sequences; e-value < 1 × 10^-5, query coverage ≥ 50%, and identity ≥ 50%). The same approach was used for vicinity analysis of the OpPKS8 locus. Afterward, BLASTP¹ (non-redundant protein sequence database) (Altschul et al., 1990) was used to search for and curate orthologous gene clusters among other genomes of filamentous fungi and to select putative orthologous backbone genes for the phylogenetic analysis of the OpPKS8 backbone gene (OphioH327gp7312) (three cutoffs were defined: e-value < 1 × 10^-5, query coverage ≥ 50%, and identity ≥ 45%). To further confirm the orthology of collected genes, amino acid sequences of all proteins per genome of species harboring putative OpPKS8 backbone gene orthologs (13 species with annotated genomes) were subjected to OrthoMCL clustering (clustering thresholds: e-value < 1 ×-10^-5 and identity ≥ 30%) (Supplementary Table 1) (Li et al., 2003). Additionally, several fungal genome sequences have been deposited at NCBI as raw or incomplete assemblies from projects that generally employed whole genome shotgun (WGS) strategies. These projects do not have predicted proteins/genes/mRNAs deposited at NCBI, and are therefore inaccessible through routine BLASTP against the non-redundant protein sequences (nr) database and BLASTN against the nucleotide collection (nr/nt) database. Thus, initial screening using BLASTN against the WGS database revealed several unannotated genomes (Purpureocillium lilacinum strain TERIBC 1, Biatriospora mackinnonii strain E5202H, Phialocephala scopiformis strain CBS 120377, Cairneyella variabilis strain VPRI 42388, Shiraia spp. strain Slf14, Endocalyx cinctus strain JCM 7946, Paecilomyces hepiali strain FENG, Pyrenochaeta lycopersici strain CRA-PAV ER-1211, Talaromyces purpureogenus strain MYA38, Aureobasidium pullulans strain AY4, and Scedosporium aurantiacum strain WM 09.24) (Chan et al., 2012; Vandeputte et al., 2014; Yang et al., 2014; Malapi-Wight et al., 2015; Pérez-Bercoff et al., 2015; Walker et al., 2016; Yu et al., 2016) that harbor putative OpPKS8 backbone gene orthologs. The content of these genomes was downloaded and further explored for the presence of OpPKS8 gene cluster orthologs (through MultiGeneBlast 1.1.14), which were annotated using FGENESH (gene-finding parameters for the closest species presented in the algorithm) (Solovyev et al., 2006). Corresponding amino acid sequences were aligned with the predicted OpPKS8 backbone gene products from O. ulmi W9 and O. novo-ulmi H327. Genes that satisfied the previously fixed cutoffs were included in further analyses.

Since the phylogeny of PKS genes can be problematic, particularly for ortholog definition, we included collected entries comprising putative orthologs of the OpPKS8 backbone gene and all characterized PKS from MIBiG, a database of characterized BGCs (Medema et al., 2015), for phylogenetic reconstruction (Sbaraini et al., 2016). This tree helps to determine if the collected entries are true orthologs or artifacts. The amino acid alignment was built using PRANK v.100701, without manual curation (Loytynoja and Goldman, 2010). The best-fit evolutionary model was estimated using Prottest 3.4 (Darriba et al., 2011) (Supplementary Table 2), and phylogenetic reconstruction was conducted using two different methods (both results are available in Supplementary Figure 1): Maximum Parsimony, as implemented in MEGA 6 (Tamura et al., 2013) with 1,000 bootstrap replicates, and Bayesian Inference using MrBayes (Altekar et al., 2004; Ayres et al., 2012; Ronquist et al., 2012). The latter method was executed for up to 10⁷ Markov chain Monte Carlo (MCMC) generations (sampled every 100 steps), applying an average standard deviation of split frequencies < 0.01, calculated every 1000th generation, as the main convergence criterion. The first 25% of samples, corresponding to the initial stages in MCMC sampling where likelihood values increase rapidly, were discarded as burn-in, summarizing parameters and trees afterward. After this phylogenetic confirmation of orthology, amino acid sequences (BLAST-collected entries for the OpPKS8 backbone gene product in addition to FGENESH-annotated entries) were manually trimmed and only the keto-synthase, acyltransferase, and dehydratase protein domains sequences were aligned. Alignment reliability was assessed with GUIDANCE 2.0 (alignment is shown in Supplementary Data Sheet 1), using PRANK (Loytynoja and Goldman, 2010) for sequence alignment with 100 bootstrap replicates and variable gap penalties and a GUIDANCE2 score cutoff of 0.90 for site removal (Landan and Graur, 2008). The best-fit evolutionary model was estimated as described above, and trees for this dataset were constructed using PhyML 3.1 (Guindon et al., 2010) with aLRT SH-like (approximate likelihood ratio test Shimodaira–Hasegawa) branch support estimation (Anisimova and Gascuel, 2006; Anisimova et al., 2011) and MrBayes, as described above.

A species tree based on 18S ribosomal RNA gene (partial), internal transcribed spacer 1, 5.8S ribosomal RNA gene, internal transcribed spacer 2, and 28S ribosomal RNA gene (partial) sequences was constructed to clarify the relationships between the Ophiostomataceae species used in comparative genomic analysis (Carlier et al., 2006; Romon et al., 2014). Additionally, a Neurospora crassa entry was added as an outgroup. The best-fit evolutionary model was assessed using jModeltest 2.1.9 (Darriba et al., 2012), and a phylogenetic tree was constructed using PhyML 3.1 with 100 bootstrap replicates (alignment shown in Supplementary Data Sheet 2). A logical diagram describing the step-by-step process and the connections between methodologies applied in this section is shown in Supplementary Figure 2.

A genome survey predicted 19 putative BGCs/backbone genes in O. ulmi W9 and O. novo-ulmi H327 genomes. Interestingly, the metabolic potential of O. ulmi W9 and O. novo-ulmi H327 seemed identical. All BGCs found in O. ulmi W9 were also found in O. novo-ulmi H327, and vice versa, and included 1 NRPS, 10 PKS, 5 terpenes (TERP), and three BGCs that were classified as “OTHER” by antiSMASH, a generic class of gene clusters grouping BGCs with unusual backbone genes (e.g., NRPS-like and PKS-like) (Figure 1 and Supplementary Data Sheet 3).

To validate the predicted backbone genes of DED pathogens and obtain insight about possible non-functional BGCs, RNA-seq data from O. novo-ulmi H327 was analyzed as it was expected that possible non-functional BGCs would not have mapped reads. All of the predicted backbone genes exhibited mapped reads in at least one of the three conditions (yeast in liquid culture medium, mycelium in static liquid culture medium, or mycelium in solid culture medium), indicating gene functionality (Supplementary Table 3). Additionally, the expression of OpPKS10 BGC and OpOTHER1 BGC appeared to be regulated by fungal lifestyle structures, since the backbone genes (OphioH327g7900 and OphioH327g5165) of these clusters are up-regulated in mycelium (mycelium in static liquid culture medium × yeast in liquid culture medium [logFC: 6.71; FDR ≤ 1%] and mycelium in solid culture medium × yeast in liquid culture medium [logFC: 5.46; FDR ≤ 1%]) and yeast (yeast in liquid culture medium × mycelium in static liquid culture medium [logFC: 1.2; FDR ≤ 1%] and yeast in liquid culture medium × mycelium in solid culture medium [logFC: 1.9; FDR ≤ 1%]), respectively (Nigg et al., 2015) (Supplementary Table 3).

To identify SM gene clusters unique to DED pathogens, we assessed BGC conservation among eight Ophiostomataceae family species. The proximity of these species to O. ulmi W9 and O. novo-ulmi H327 was confirmed by phylogeny (Supplementary Figure 3), and completeness and quality of genome sequence assembly/predicted proteins was assessed using BUSCO v2 (all genomes had >89% of the benchmarking universal single-copy orthologs expected to be conserved in Sordariomycetes; Supplementary Table 1). The majority of BGCs (>89%) found in DED pathogens were well-conserved in the saprophytic O. piceae UAMH 11346 (18 conserved BGCs) and in species from the closely related Sporothrix genus (17 conserved BGCs in S. schenckii 1099-18, S. brasiliensis 5110, and S. pallida SPA8) (Figure 1). Interestingly, more than 50% of the predicted BGCs detected in DED pathogens were also conserved in more distantly related species such as G. fragrans CBS 138720 (14 conserved BGCs), G. clavigera kw1407 (11 conserved BGCs), L. lundbergii CBS 138716 (11 conserved BGCs), and L. procerum CMW34542 (11 conserved BGCs) (Figure 1). These results indicate that the SM potential is remarkably similar between Ophiostomataceae family members. Notably, two BGCs were found only in the Ophiostoma genus: OpOTHER3 (OphioH327gp1586) was identified in O. piceae UAMH 11346, O. ulmi W9, and O. novo-ulmi H327, while OpPKS8 is exclusive to DED pathogens. To further confirm the absence of OpPKS8 in O. piceae UAMH 11346, the vicinities of the OpPKS8 locus were explored. The results confirmed the absence of the OpPKS8 gene cluster in O. piceae UAMH 11346, as several genes in the vicinity were conserved despite the absence of OpPKS8 (Supplementary Figure 4).

Considering that no putative BGCs have been functionally characterized in DED pathogens to date, we employed comparative genomics to determine the likely final products of these BGCs. These comparisons revealed three interesting BGCs that exhibited previously characterized orthologous gene clusters in other species and which were putatively linked with the biosynthesis of a fujikurin-like compound (OpPKS8), a siderophore compound (OpNRPS1), and pyrones and resorcylic acids (OpPKS10, a type III PKS).

Orthologs for the OpPKS10 backbone gene were found in all evaluated Ophiostomataceae species (Figure 1). The predicted backbone gene protein for this cluster (OphioH327gp7900) exhibited 63% identity with the protein sequence of a type III PKS involved in the biosynthesis of triketide pyrones, tetraketide pyrones, tetraketide resorcylic acids, and pentaketide resorcylic acids characterized in N. crassa (Goyal et al., 2008; Rubin-Pitel et al., 2008). N. crassa and the Ophiostomataceae family species are grouped in the same class, which may indicate that proteins with similar functions to those of N. crassa could represent the final product of the OpPKS10 backbone gene, and that these similar proteins could act in the biosynthesis of similar metabolites. Synteny comparison between the type III PKS locus in DED pathogens and in N. crassa did not reveal extensive conservation (Supplementary Figure 4), suggesting that it might not be grouped in a biosynthetic cluster. Reinforcing this hypothesis, CASSIS was unable to discriminate homologous regulatory sequences between type III PKS of DED pathogens and neighboring genes.

OpNRPS1 BGC shows similarity to characterized NRPS involved in siderophore biosynthesis (Supplementary Figure 4). Iron is an essential micronutrient for almost all organisms, including fungi (Silva-Bailao et al., 2014). Thus, these microorganisms have developed highly efficient iron-acquisition systems. Siderophores are chelators synthesized by microbes to sequester iron (Haas, 2014). Siderophore biosynthesis has been widely explored in model fungi such as Aspergillus fumigatus, Paracoccidioides brasiliensis, and Metarhizium anisopliae (Haas, 2014; Silva-Bailao et al., 2014; Giuliano Garisto Donzelli et al., 2015; Sbaraini et al., 2016). In these model species, siderophores are essential for successful infection (Schrettl et al., 2007; Giuliano Garisto Donzelli et al., 2015). Primary screening (BLASTP cutoffs: e-value < 1 × 10^-5, query coverage ≥ 50%, and identity ≥ 45%; searching against characterized gene clusters/backbone genes deposited in MIBiG) did not reveal genes for siderophore biosynthesis in DED pathogens. However, in view of the importance of siderophores for the infection of several fungal species, we deepened our analysis by decreasing the previously fixed cut-offs (BLASTP cut-offs: e-value < 1 × 10^-5, query coverage ≥ 50%, and identity ≥ 30%). Additionally, since siderophore biosynthesis in filamentous fungi is usually catalyzed by NRPS and DED pathogens have a single NRPS gene cluster (OpNRPS1), we focused on this BGC. Secondary screening revealed that the protein sequence of the OpNRPS1 backbone gene from O. ulmi W9 exhibited 38% identity with the backbone gene of the intracellular siderophore ferricrocin (KFG82816.1), which has been characterized in M. robertsii and is conserved in M. anisopliae (Giuliano Garisto Donzelli et al., 2015; Sbaraini et al., 2016). Additional analysis of OpNRPS1 gene cluster revealed an L-ornithine N5-oxygenase gene (OphioH327gp2447) in the BGC (Supplementary Figure 4). The product of this gene exhibited 48% identity with an L-ornithine N5-oxygenase (SidA) from Aspergillus nidulans (Eisendle et al., 2003). It has been hypothesized that, in A. nidulans, the first biosynthetic step for both intra- and extracelluar siderophores is catalyzed by ornithine oxygenates such as SidA (Plattner and Diekmann, 1994). In agreement with this hypothesis, mutational inactivation of SidA abolishes siderophore biosynthesis, completely attenuating the virulence of A. nidulans to mice (Eisendle et al., 2003). This clustering of the ferricrocin ortholog backbone gene and the L-ornithine N5-oxygenase gene suggests, therefore, that DED pathogens utilize a siderophore biosynthesis pathway that is similar to other model fungi. Furthermore, we hypothesize that the products of OpNRPS1 BGC modulate iron dynamics in DED pathogens.

The OpPKS8 gene cluster matched the characterized fujikurin BGC from Fusarium fujikuroi IMI 58289, and showed considerable identity (36–68% identity in protein-by-protein comparisons of BGC constituent genes), albeit a different cluster configuration (Figure 2) (Wiemann et al., 2013; von Bargen et al., 2015). von Bargen et al. (2015) have proposed that the fujikurin gene cluster (also called PKS19) of F. fujikuroi IMI 58289 is composed of five genes: a highly reducing PKS (CCT72377) exhibiting 58% identity with the OpPKS8 backbone gene (OphioH327gp7312); a protein with a predicted enoyl reductase domain (CCT72378) exhibiting 63% identity with the ortholog present in OpPKS8 BGC (OphioH327gp7311); a protein with a predicted hydrolase domain (CCT72379) exhibiting 50% identity with the ortholog present in OpPKS8 BGC (OphioH327gp7310); a putative transcription factor (CCT72381) exhibiting 36% identity with the ortholog present in OpPKS8 BGC (OphioH327gp7309); and a putative transporter (CCT72380) that does not have an ortholog in OpPKS8 BGC (von Bargen et al., 2015). Thus, the basic components of the fujikurin biosynthetic route appear to be conserved in DED pathogens. Additionally, comparative genomics and phylogeny have revealed that the OpPKS8 gene cluster is phylogenetically related to putative BGCs found in Scedosporium apiospermum IHEM 14462 (66–84% identity in protein-by-protein comparisons of constituent genes of these BGCs), Aureobasidium pullulans AY4 (58–80% identity) and Setosphaeria turcica Et28A (66–70% identity) (Figures 2, 3). Some genes conserved between the OpPKS8 gene cluster and the BGCs from S. apiospermum, A. pullulans, and S. turcica coded for proteins with known implications in SM biosynthesis. This led to the hypothesis that other genes may be involved in the biosynthesis of the final compound of the OpPKS8 gene cluster (unconserved genes in F. fujikuroi), namely, a fujikurin-like compound. In addition to the previously cited genes, the proposed OpPKS8 BGC contains the following: a putative O-methyltransferase protein (OphioH327gp7306) which exhibits 77% identity with an ortholog from S. apiospermum KEZ41288; absent in A. pullulans and S. turcica); a protein with a predicted hydrolase domain (OphioH327gp7305) exhibiting 79% identity with an S. apiospermum ortholog (KEZ41289); a putative MFS transporter protein (OphioH327gp7308) 75% identical to an ortholog from S. apiospermum (KEZ41294); and a putative monooxygenase protein (OphioH327gp7307) exhibiting 66% identity with an ortholog in S. apiospermum (KEZ41295; absent in S. turcica) (Figure 2).

Moreover, it is important to note that S. apiospermum belongs to a different order (Microascales), and that A. pullulans and S. turcica belong to a different class (Dothideomycetes), and are therefore considerably distant from DED pathogens (Ophiostomales, Sordariomycetes). Evolutionary distances are best displayed in the schematic Ascomycota tree of life based on the work of Schoch et al. (2009) (Supplementary Figure 5), which contrasts with the OpPKS8 backbone gene tree (Figure 3) (Schoch et al., 2009). Additionally, the OpPKS8 gene cluster is the only BGC exclusively found in DED pathogens (among the evaluated species), raising the possibility that HGT events could have shaped the evolution of OpPKS8 orthologs. Alternatively, an ancestor of the Sordariomycetes and Dothideomycetes classes may have contained a fujikurin-like compound gene cluster, which remained in the genome of some species. This hypothesis, however, would have required the loss of constituent genes across several Sordariomycetes and Dothideomycetes species. Additionally, species from Leotiomycetes (Hymenoscyphus fraxineus, Cairneyella variabilis, and Phialocephala scopiformis) and Eurotiomycetes (Talaromyces cellulolyticus, Talaromyces purpureogenus, and Penicillium griseofulvum) classes also have orthologs, making this vertical descent hypothesis unlikely.

Elucidation of the structure of fujikurin has revealed that this compound was previously described as CR377, an SM with antifungal activity isolated from an endophytic Fusarium species (Brady and Clardy, 2000). Furthermore, the F. fujikuroi IMI 58289 BGC (PKS19) involved in fujikurin production was up-regulated during the infection of rice plants, highlighting a possible role of this metabolite as a phytopathogenic virulence determinant (Wiemann et al., 2013). We sought to analyze the lifestyles (pathogenic or not) of those fungi harboring OpPKS8 backbone gene orthologs. Although cluster configuration changes in these species (Figure 2), lifestyle analyses could aid in the prediction of a role for the metabolic product of the OpPKS8 BGC. We were particularly interested in whether other phytopathogenic fungi harbor fujikurin or fujikurin-like compound gene clusters, in view of the putative role of fujikurin in the infection of Fusarium species. Regarding BGCs closely related to OpPKS8 (shaded in gray in Figure 3), five gene clusters have been found in species with a known phytopathogenic trait; three in widely known opportunistic pathogens (although both A. pullulans and Scedosporium spp. are associated with plants); one in Endocalyx cinctus (a poorly documented species linked with saprophytic growth in dead palms); one in the ericoid mycorrhizal fungus Cairneyella variabilis; and one in Paecilomyces hepiali (a species with nematophagous/entomopathogenic traits that may be associated with plants). Regarding BGCs more distantly related to the OpPKS8 gene cluster: 16 have been found in phytopathogens, with several Fusarium representatives; three in endophytic fungi species; two in saprophytic species; two in opportunistic pathogens; two in widely known mycopathogens; and one in a species with nematophagous/entomopathogenic traits (Figure 3). Clearly, several species that possess orthologous genes have a plant-associated lifestyle, leading us to speculate that fujikurin-like SMs allow several fungal species to interact with plants, not only as phytopatogens or endophytes but also in a nematophagous/entomopathogenic or opportunistic manner.