Genome analysis of the ubiquitous boxwood pathogen Pseudonectria foliicola

Fungal isolate and nucleic acid isolation

An axenic culture of P. foliicola isolate ATCC13545^® (also known as isolate A.R. 2711) was used for genome sequencing. This fungal isolate was originally cultured from B. sempervirens in Maryland, US. The isolate was grown on potato dextrose agar (BD Difco™, Sparks, MD, USA) for 5-days under 12-h white light photoperiod and then transferred to yeast extract potato dextrose liquid media at 25 °C for 2-days under continuous light. Genomic DNA was extracted from hyphal tissue harvested from liquid media using the OmniPrep DNA kit (G-Biosciences, St. Louis, MO, USA) according to manufacturer’s instructions, and subsequently purified using the Zymo Genomic DNA Clean and Concentrator kit (Zymo Research, Irvine, CA, USA).

Whole genome sequencing and de novo assembly

A genomic DNA library was constructed using the TruSeq Nano DNA Library Prep kit (Illumina, Inc., San Diego, CA, USA) and quantified using the Qubit 2.0 fluorometer (Life Technologies, Grand Island, NY) and the LabChipXT DNA 750 (Caliper Life Sciences, Hopkinton, MA, USA). The library was sequenced on an Illumina MiSeq in two independent runs using paired-end 600-cycle reagent cartridge v.3 (Illumina, Inc.). Reads were processed and assembled using CLC Genomics Workbench v.7.5.1 (CLC Bio, Boston, MA, USA) with a k-mer of 24 and a minimum contig length of 500 bp. Illumina adapters were trimmed and low quality reads (Phred score < 0.05) were removed. Summary statistics for the draft genome were generated using CLC Genomics Workbench and QUAST (Gurevich et al., 2013).

Completeness of the P. foliicola draft genome assembly was evaluated using BUSCO v.1.1b1 (Simão et al., 2015). The genome assembly of P. foliicola used in this study was deposited in NCBI GenBank under accession LMTV00000000, and datasets are also available at the US National Agricultural Library on AgData Commons (http://dx.doi.org/10.15482/USDA.ADC/1408094).

Nuclear genome annotation

Ab initio gene predictions for the draft genome assembly of P. foliicola were performed using the MAKER2 v.2.31.6 annotation pipeline (Holt & Yandell, 2011). Gene training was performed after running three rounds of the program according to the program documentation. Gene boundaries were assigned using protein homology evidence from Fusarium graminearum strain PH-1 (NCBI BioProject Acc: PRJNA13839; Cuomo et al., 2007). Additional ab initio gene predictions were made using the program SNAP (http://korflab.ucdavis.edu/software.html) and AUGUSTUS v.3.2.1 (Stanke et al., 2004) with F. graminearum set as the prediction species model organism. New gene predictions using MAKER2 were also performed for the previously published genome assembly of D. macrodidyma (Malapi-Wight et al., 2015) using the same parameters as above.

Identification of transposable elements and repeat-induced mutations

The presence of transposable elements (TEs) was evaluated from the P. foliicola and D. macrodidyma genomes using the REPET v2.5 (Flutre et al., 2011) pipeline, along with supporting databases Repbase v.20.05 (Kapitonov & Jurka, 2008; Bao, Kojima & Kohany, 2015) and Pfam v.27.0, and run according to the accompanying program documentation (https://urgi.versailles.inra.fr/Tools/REPET). The TEdenovo stage was used first to produce a database of four de novo identified TEs: an incomplete helitron, a MITE (miniature inverted-repeat TE), a TRIM (terminal-repeat retrotransposon in miniature), and one uncategorized TE sequence. These sequences were then used as initial input for the first of two runs of the TEannot pipeline stage, which produced a final genome-wide GFF annotation file using the TE classification scheme described by Wicker et al. (2007). A custom script was used to tabulate counts of TE classes, orders and superfamilies and filtered out fragments less than 80-bp in length according to Wicker et al. (2007) recommendations to avoid misclassification. Results from TEannot based on tblastx and blastx (BLAST+ vr. 2.2.31; Camacho et al., 2009) searches against Repbase were also filtered to remove hits sharing <70% sequence identity.

Individual TE families with ten or more sequences identified, including at least one sequence ≥ 300-bp in length, were checked for signatures of repeat-induced point (RIP) mutation activity using the RIPCAL vr. 2 program (Hane & Oliver, 2008; Hane, 2015). The 300-bp length cutoff was selected based on RIPCAL’s default setting for scanning subsequences and the program was run using both the alignment mode consensus model with ClustalW (Larkin et al., 2007) and di-nucleotide frequency-based methods. Evidence of the signature of RIP mutations was present if di-nucleotide frequencies matched the indexes: (TpA / ApT) ≥ 0.89 and (CpA + TpG) / (ApC + GpT) ≤ 1.03 (Margolin et al., 1998) and visual inspection of RIPCAL alignments showed one or more peaks for (CA ←→TA) + (TG ← →TA) mutations (Hane & Oliver, 2008).

Mitochondrial genome

The P. foliicola and D. macrodidyma mitogenomes were identified by performing tBLASTx searches of the complete genome assemblies using the 95.7 kb F. graminearum mitochondrial genome as a query (Al-Reedy et al., 2012), with the genetic code set to four (mold mitochondrial). Mitogenomes were annotated using MITOS (Bernt, Donath & Jühling, 2013). Comparative analysis of genes, rRNA and tRNA was performed in the program SimpleSynteny (Veltri, Malapi-Wight & Crouch, 2016) using the P. foliicola CDS file to annotate all genomes (e-value cutoff 1e-5, minimum query cutoff 10%), with circular genome mode implementation and genomes organized to minimize Euclidean line distance.

Identification of mating type idiomorphs

We assessed the presence of the MAT1-1 and MAT1-2 idiomorphs in the P. foliicola genome by performing a local BLASTn search (e-value cutoff 1e-5) against a MAT gene database. The database contained nucleotide sequences for the highly conserved alpha domain DNA binding motif (MAT1-1-1) and the high mobility group (HMG) box DNA binding motif (MAT1-2-1) for 13 different filamentous fungal species, retrieved from the NCBI GenBank database (Files S1 and S2).

Phylogenetic reconstruction

A phylogenetic analysis was used to illustrate the relationship between P. foliicola and 14 other fungal species. For this analysis, the predicted proteomes of Aspergillus nidulans FGSC A4 (ASM114v1; Galagan et al., 2005), Botrytis cinerea BcDW1 (Assembly GCA000349525; Blanco-Ulate et al., 2013), F. graminearum PH-1 (GCA000240135; Cuomo et al., 2007), Macrophomina phaseolina MS6 (GCA000302655; Islam et al., 2012), Magnaporthe oryzae 70-15 (MG8; Kim et al., 2010), Neurospora crassa (GCA000786625; Baker et al., 2015), Penicillium oxalicum 114-2 (GCA000346795; Liu et al., 2013), Pyrenophora tritici-repentis (GCA000149985; Manning et al., 2013), Sclerotinia sclerotiorum 1980 UF-70 (ASM1469v1; Amselem et al., 2011), Trichoderma reesei RUT C-30 (GCA000513815; Koike et al., 2013), Ustilago maydis 521 (UM1; Kämper et al., 2006), Verticillium dahliae JR2 (GCA000400815; De Jonge et al., 2012) and Yarrowia lipolytica CLIB122 (GCA000002525; Dujon et al., 2004) were downloaded from either the EnsemblFungi database (https://fungi.ensembl.org/index.html) or the NCBI Genbank Genome database (https://www.ncbi.nlm.nih.gov/genbank/). The predicted proteome of D. macrodidyma generated in our study was also included in the analysis. All proteomes were searched against each other using BLASTp (e-value 0.001) and clustered in orthologous gene sets using OrthoMCL v1.4 in the CYVERSE Discovery Environment (https://de.cyverse.org/de/). Single copy genes (clusters with exactly one member per species) found in all 15 fungal proteomes were extracted from the orthologous dataset and the amino acid alignments were performed using MUSCLE v3.8.31 (Edgar, 2004). Gblocks v.0.91b was used to remove ambiguously aligned regions using relaxed selection parameters following Talavera & Castresana (2007). Maximum likelihood (ML) phylogenetic analyses were performed in RAxML (Stamatakis, 2006), using the RAxML GUI v. 1.5b1 (Silvestro & Michalak, 2012). Phylogenetic trees were constructed using the JTT matrix-based model (Jones, Taylor & Thornton, 1992) and 1,000 bootstrap replicates. Ustilago maydis 521 was used as outgroup in the phylogenetic analyses.

Estimation of evolutionary divergence times

To obtain approximate information on divergence events, the phylogenetic tree was timed using RelTime (Tamura et al., 2012) as implemented in MEGA 7 (Kumar, Stecher & Tamura, 2016). RelTime estimated the relative divergence times for each node of the ML tree using the same outgroup taxa. We used seven confidence intervals obtained from the TimeTree database (http://timetree.org; Hedges et al., 2015) as minimum and maximum times to convert the relative times into absolute times (Table S4). Time estimates were performed using maximum-likelihood branch length, local clocks, a JTT matrix-based model (Jones, Taylor & Thornton, 1992) and a discrete Gamma distribution among sites (five categories).

Comparative genomic analyses

The genome sequences of D. macrodidyma JAC15-245 (NCBI GenBank accession JYGD00000000) and F. graminearum PH-1 were downloaded from NCBI GenBank and Ensembl-Fungi database respectively, and used for comparative analysis against the P. foliicola draft genome assembly generated in this study. The predicted proteome generated from this study was used for D. macrodidyma, while the published proteome dataset for F. graminearum was downloaded from the Ensembl-Fungi database. Genome-wide identification, comparison and visualization of orthologous gene clusters among P. foliicola, D. macrodidyma and F. graminearum were performed using the web platform OrthoVenn (http://www.bioinfogenome.net/OrthoVenn/; Wang et al., 2015). The program uses a modified OrthoMCL algorithm to identify orthologous gene clusters from the UniProt/Swiss-Prot database. Proteins potentially involved in carbohydrate metabolism were annotated for each proteome by searching against the database of automated Carbohydrate-active enzyme ANnotation (dbCAN; http://csbl.bmb.uga.edu/dbCAN/annotate.php; Yin et al., 2012). A X² test of independence, similar to that used in Martinez et al. (2008), was used to identify statistically significant differences across the CAZyme repertoire. Identification of putative enzymes related to the biosynthesis of secondary metabolites was performed from the predicted proteins using the web-based program AntiSMASH (http://antismash.secondarymetabolites.org; Medema et al., 2011; Weber et al., 2015). The secretome was predicted by screening the predicted proteins for different features using a bundle of eight different prediction tools implemented in the web-based program SECRETOOL (Cortázar et al., 2014). Pathogenicity associated genes were identified by performing a local BLASTp search of the predicted proteome against the curated pathogen-host interaction database (PHI-base v.4.0; http://www-phi4.phibase.org/; Winnenburg et al., 2008).

Nuclear genome assembly and annotation

A de novo genome assembly of P. foliicola ATCC13545^® was generated from 19.5 million paired end, 300-bp reads, comprising a total of 5.4 Gb of raw sequence data. The resulting genome assembly for P. foliicola was 28.7 Mb, organized in 425 scaffolds (≥500-bp), with an average read depth coverage of 188-fold (Table 1). The N50 scaffold length is 184,987 bp and the three longest scaffolds are 619,696 bp, 551,775 bp, and 524,658 bp. Analysis of the P. foliicola genome assembly using BUSCO identified 420 complete genes out of 429 conserved eukaryotic ortholog dataset and 1,416 complete genes out of 1,438 conserved fungal ortholog dataset (genome completeness scores of 97% and 98%, respectively). Ab initio gene predictions performed with the MAKER2 pipeline for the P. foliicola genome assembly identified 9,272 protein-coding genes with an average predicted protein length of 484 amino acids (Table 1). Meanwhile, gene predictions for the genome of D. macrodidyma identified 16,959 protein-coding genes with an average length of 478 amino acids. Of the three genomes evaluated in this study, gene number predictions were highest for D. macrodidyma, which also had the largest estimated genome size at 58 Mb when compared to P. foliicola and F. graminearum (13,313 genes, 36.1 Mb; Cuomo et al., 2007).

Transposable elements, repetitive DNA and repeat induced mutations

Only 0.7% (196,205-bp) of the P. foliicola nuclear genome assembly contained TEs based on the REPET pipeline analysis (Table 1). Annotation of TEs across the P. foliicola genome assembly identified 191 TE matches. Of these, 141 were Class I (retrotransposons) TEs comprising ∼0.4% of genome and 45 were Class II (DNA) TEs making up ∼0.3% of the genome. Long terminal-repeats (LTRs) were found to be the most abundant Class I order (∼0.3% of genome) and included matches in the Copia, Gypsy and BEL/Pao retrotransposon superfamilies. Almost all Class II TEs were identified as Helitrons, with only ∼0.001% of the genome identified with TIR or “unknown” TEs (Wicker code: DXX). Evidence of RIP mutation affecting P. foliicola TEs was only found for Helitrons based on alignment of 37 sequences (3.7-kb in length) and di-nucleotide indexes: TpA/ApT = 1.8 and (CpA + TpG)/(ApC + GpT) = 0.3.

The draft genome of D. macrodidyma genome as originally published did not include information about TEs. Here, we also analyzed that assembly for the presence of TEs and RIP mutations using the same parameters used for P. foliicola for comparative purposes. As REPET identified over 2,000 D. macrodidyma TEs of “unknown” superfamily, we performed an additional tblastx search (e-value: 0.01, percent query coverage per hsp: ≥70%, minimum percent identity of hits: ≥70%) against Repbase and reclassified five TEs based on the top match if it shared the same class and order. REPET analysis identified ∼6.5% (3,794,887-bp) of the D. macrodidyma genome to be made up of TEs (2,178 elements). Of these, 1,489 elements were Class I TEs and comprised ∼5.3% of the genome and 689 were Class II TEs and comprised ∼1.2% of the genome. The most common Class I retrotransposon orders were: LTRs (∼3.2%), “Unknown RXX” (∼1.1%), DIRS-like (∼0.9%), LINE (∼0.1%) and SINE (0.001%). For Class II DNA elements, the most common orders were: TIR (∼0.7%), Helitron (∼0.4%) and “Unknown DXX” (∼0.2%). Evidence for RIP mutation affecting D. macrodidyma TEs using alignment and di-nucleotide counts was identified in the sets of 91 Helitron and 146 DIRS TEs. For the Helitrons, RIP indexes were calculated as: TpA/ApT = 1.3 and (CpA + TpG)/(ApC + GpT) = 0.9, while for DIRS: TpA/ApT = 1.4 and (CpA + TpG)/(ApC + GpT) = 0.3.

Mitochondrial genome size and gene content

The mitogenomes of P. foliicola and D. macrodidyma were each contained within a single scaffold, measuring 58.3 kb and 44.2 kb respectively (scaffold 56, scaffold 68). Although the mitogenome-containing scaffolds for these two organisms were considerably smaller than that of F. graminearum (95.7 kb), a full complement of protein-coding genes was contained: apocytochrome b (cob), ATP-synthase subunits (atp6, atp8, atp9), cytochrome oxidase subunits (cox1, cox2, cox3), and NADH subunits (nad1, nad2, nad3, nad4, nad4L, nad5, nad6). This set of 14 protein-coding mitochondrial genes is highly conserved among fungi, and shared with animal mtDNA (Bullerwell & Lang, 2005). A full complement of tRNAs necessary for translation was encoded in the mitogenomes (25 total), as were small and large subunit rRNAs (one and three copies, respectively). Neither P. foliicola nor D. macrodidyma mitogenomes contained copies of the four large unidentified open-reading frames that are encoded by the F. graminearum mitogenome (Al-Reedy et al., 2012), accounting for the observed difference in sizes between these species.

Mitochondrial genome synteny

The organization of the mitogenomes of P. foliicola, D. macrodidyma and F. graminearum was highly conserved at the gene level. As observed in the mitochondria of most ascomycete fungi, all genes, rRNAs and tRNAs were oriented in the same direction. The three mitogenomes contained a nearly identical ordering of shared genes (Fig. 2). Only a single tRNA (tRNA-G) was positioned differently in the P. foliicola and D. macrodidyma mitogenomes relative to F. graminearum. Two tRNAs (tRNA-R and tRNA-L) were positioned differently between D. macrodidyma and P. foliicola/F. graminearum. This conserved ordering of the three mitogenomes was consistent with previous observations that the mitochondrial DNA of fungi in the order Hypocreales is highly conserved (Al-Reedy et al., 2012; Pantou, Kouvelis & Typas, 2008).

Identification of mating type idiomorphs in P. foliicola

The presence of the conserved alpha domain indicative of the MAT1-1 mating type idiomorph was identified in P. foliicola scaffold 102. In the same scaffold, the APN2 (encoding DNA lyase) and SLA2 (encoding cytoskeletal protein) genes were found flanking the mating type idiomorph, as well as the MAT-locus associated gene COX13 (cytochrome c oxidase subunit VIa homolog). The HMG-box region indicative of the MAT1-2 mating type idiomorph was not found within the P. foliicola genome. Based on these data, P. foliicola appears to be a heterothallic fungus, requiring a partner of the alternate mating type in order to initiate the sexual life cycle.

Phylogeny and divergence estimation

The phylogenetic relationships and divergence times among the fungal genomes studied is displayed in Fig. 3. Fourteen publicly available fungal genomes were used to examine the phylogenetic placement of P. foliicola through the analysis of single copy orthologous genes. The program OrthoMCL identified 16,356 gene clusters, from which 1,884 orthologous genes were shared across all 15 fungal species. From these shared gene clusters, 1,511 orthologous genes present as single copies were used for the phylogenetic analysis. The final dataset after removal of ambiguously aligned regions consisted of 388.7 Mb. The ML analysis identified, with high bootstrap support (>70%), five major clusters representing the ascomycete classes Sordariomycetes, Leotiomycetes, Eurotiomycetes, Dothideomycetes, and Saccharomycetes (Fig. 3). Within the Sordariomycetes, P. foliicola was more closely related, albeit appearing basal, to D. macrodidyma and F. graminearum, all of which belonged to the Nectriaceae in the order Hypocreales. These phylogenetic relationships were consistent with previous reports (Lombard et al., 2015; Yin et al., 2015). Using RelTime methods and a JTT matrix-based model, the estimated log likelihood value was −142633.0660. The divergence of P. foliicola from the Nectriaceae species D. macrodidyma and F. graminearum was estimated to have occurred ∼132 Mya. The overall estimates of divergence times in the tree are in agreement with divergence times reported for the order Hypocreales (Sung, Poinar & Spatafora, 2008).

Comparative genomic analysis of P. foliicola and other fungi in the Nectriaceae

Gene orthology

Orthologous gene clusters were identified for P. foliicola, D. macrodidyma and F. graminearum using OrthoVenn (Fig. 4). Proteins from these three fungal species formed 10,403 orthologous clusters, of which 7,135 clusters were shared among all three species. The top three Swiss-Prot annotations among the core shared clusters were the Acyl-CoA-binding domain-containing protein (44 proteins), pleiotropic drug resistance protein 4 (30 proteins) and a short-chain dehydrogenase TIC 32, chloroplastic (12 proteins). Sixteen clusters representing 48 predicted proteins were unique to P. foliicola, while 610 and 124 species-specific protein clusters were identified in D. macrodidyma and F. graminearum, respectively. Only 0.2% of the P. foliicola proteome was unique, a low percentage relative to 1.3% for the F. graminearum and 6.0% for the D. macrodydima proteomes. Although the majority of the 16 clusters unique to P. foliicola had no annotations based on the UniProt/Swiss-Prot and Gene Ontology (GO) databases, three of the unique clusters were identified as (1) vegetative cell wall protein Gp1/structural constituent of cell wall (GO:0005199), (2) thioredoxin/protein disulfide oxidoreductase activity (GO:0009507), and (3) leucine-rich repeat extensin-like protein 3/structural constituent of cell wall (GO:0005618). Biological processes and molecular functions annotated by the GO database for the species-unique gene clusters were most abundant in D. macrodidyma (66 biological processes, 28 molecular functions), and least abundant in P. foliicola (five biological processes, two molecular functions; Table S1). Three GO categories were found enriched (hypergeometric test on OrthoVenn, p-value <0.05) in the P. foliicola unique clusters: (1) a glycerol ether metabolic process, (2) structural constituent of cell wall and, (3) a protein disulfide oxidoreductase activity. Despite the large number of unique gene clusters found in D. macrodidyma, only one GO category, an oxidoreductase activity acting on single donors with incorporation of molecular oxygen, was found to be enriched.

The CAZyme repertoire

Annotation of the P. foliicola proteome identified 448 CAZyme modules, including domains encoding 179 glycoside hydrolases (GH), 88 glycosyl-transferases (GT), 18 polysaccharide lyases (PL), 64 carbohydrate esterases (CE), 46 carbohydrate-binding modules (CBM), and 53 enzymes with auxiliary activities (AA) (summarized in Table 2, Fig. 5 and Table S2). The D. macrodidyma and F. graminearum genomes encoded 1,086 and 707 CAZyme modules, respectively. The CAZyme encoding genes represented between 4.8 and 6.6% of the predicted proteome for the three fungal species (Table 2).

Table 2:

Summary of the carbohydrate-active enzyme (CAZyme) modules identified from the predicted proteome of Pseudonectria foliicola, Dactylonectria macrodidyma and Fusarium graminearum.

RF%: Relative frequency of CAZyme modules over the total number of predicted proteins for the corresponding genome.

Name	Pseudonectria foliicola	RF%	Dactylonectria macrodidyma	RF%	Fusarium graminearum	RF%
AA	53	0.52	147	0.84	110	0.79
CBM	46	0.53	116	0.73	80	0.65
CE	64	0.83	210	1.37	127	1.08
GH	179	1.94	447	2.59	268	2.01
GT	88	1.01	125	0.78	100	0.82
PL	18	0.17	41	0.20	22	0.16
Total	448	4.83%	1,086	6.40%	707	5.31%

DOI: 10.7717/peerj.5401/table-2

Notes:

AA

Auxiliary activity families

CBM

Carbohydrate-binding modules

CE

Carbohydrate esterase families

GH

Glycoside hydrolase families

GT

Glycosyltransferase families

PL

Polysaccharide lyase families

An increased number of CAZyme modules was identified from the D. macrodidyma proteome when compared to P. foliicola, F. graminearum and a database of 91 other plant pathogenic, facultative pathogenic, saprophytic and symbiotic fungi (Table S2; Zhao et al., 2013). Although not significant in X² tests against P. foliicola and F. graminearum, D. macrodidyma has an increased number of GH and CE enzymes. This is largely due to the GH3 family, which is associated with cellulose degrading activities, the GH28 family, with pectinase activities, as well as the less characterized GH78 and GH109 families. Within the CE and PL module groups, increased numbers in the D. macrodidyma genome were observed in the families CE3, CE5, CE10 and PL1 (Fig. 5). The CE3 and CE5 have acetyl xylan esterase and cutinase activities and are common modules with a higher representation in Ascomycetes relative to Basidiomycetes (Zhao et al., 2013). CE10 families have carboxylesterases activities but can also act on non-carbohydrate substrates (Cantarel et al., 2009). The PL1 family is the most commonly found among fungi, particularly in plant pathogenic species (Zhao et al., 2013). The number of CE families in D. macrodidyma (233) is comparable to that of the pea root pathogen Fusarium solani, previously regarded as having the most CEs with 223 (Zhao et al., 2013).