Comparative analysis of genome sequences of the conifer tree pathogen, Heterobasidion annosum s.s.

The causal agent of root and butt rot of conifer trees, Heterobasidion annosum, is widespread in boreal forests and economically responsible for annual loss of approximately 50 million euros to forest industries in Finland alone and much more at European level. In order to further understand the pathobiology of this fungus at the genome level, a Finnish isolate of H. annosum sensu stricto (isolate 03012) was sequenced and analyzed with the genome sequences of 23 white-rot and 13 brown-rot fungi. The draft genome assembly of H. annosum has a size of 31.01 Mb, containing 11,453 predicted genes. Whole genome alignment showed that 84.38% of H. annosum genome sequences were aligned with those of previously sequenced H. irregulare TC 32-1 counterparts. The result is further supported by the protein sequence clustering analysis which revealed that the two genomes share 6719 out of 8647 clusters. When sequencing reads of H. annosum were aligned against the genome sequences of H. irregulare, six single nucleotide polymorphisms were found in every 1 kb, on average. In addition, 98.68% of SNPs were found to be homo-variants, suggesting that the two species have long evolved from different niches. Gene family analysis revealed that most of the white-rot fungi investigated had more gene families involved in lignin degradation or modification, including laccases and peroxidase. Comparative analysis of the two Heterobasidion spp. as well as white-/brown-rot fungi would provide new insights for understanding the pathobiology of the conifer tree pathogen.

Introduction

Heterobasidion annosum is the causal agent of root and butt rot of conifer trees in the forests of northern Europe and North America, causing annual economic loss of 800 million euros in Europe alone [1], [2]. Heterobasidion species complex comprises the European H. annosum, H. abietinum, H. parviporum, and the North American H. irregulare and H. occidentale. They exhibit differences in geographic distribution, lifestyles, interfertility, and host ranges [3]. Heterobasidion spp. have been studied with regard to ecology, control methods, physiology, and evolution [4], [5], [6]. Due to their socioeconomic and biological importance, Heterobasidion species complex has been one of the most studied forest fungi [3].

Recently H. irregulare TC 32-1 was fully sequenced, which was the first sequenced basidiomycete forest phytopathogen [7]. With the aid of this, genome-wide association study was able to identity single nucleotide polymorphisms (SNPs) associated with fungal virulence [8]. Transcriptomic profiling became more informative by using the gene models from the taxonomically close species, vitalizing gene expression-associated studies [9], [10], [11]. However, the sequenced North American isolate H. irregulare is genetically and physiologically different from the European H. annosum [12], [13]. The more genome sequences from various geographic locations would enable more accurate analysis and capture species-specific characteristics as well.

Here, we report the genome sequences of a Finnish isolate of H. annosum sensu stricto (isolate 03012) to further understand the pathobiology of this fungus at genome level. In order to assess the macroscopic similarity between H. annosum and H. irregulare, genome-wide alignment and proteins clustering were conducted. Gene synteny analysis showed species-specific genes, implying genomic rearrangements after speciation. In addition, SNPs and short insertions/deletions (InDels) were identified to elucidate microscopic differences between the two species. With the genome sequences of 35 wood-rotting fungi obtained from JGI MycoCosm [14], comparative gene family analysis was conducted to review the white- and brown-rot classification, focusing on their capability to degrade plant materials. As a result, white-rot fungi had more gene families involved in lignin degradation or modification, including laccases and peroxidase, than brown-rot fungi. Collectively, the results presented in this study may provide new insights for understanding the pathobiology of this conifer tree pathogen.

Results and discussion

Comparisons between the genome and proteome sequences of H. annosum and H. irregulare

The genome of H. annosum sensu stricto (isolate 03012) was assembled into 31,007,627 bp with GC content of 51.55%, representing 66.5-fold coverage. Subsequent gene prediction showed that H. annosum was predicted to have 11,453 protein-coding genes. A series of genomic analyses were carried out so as to identify similarities and differences between the genome sequences of H. annosum and H. irregulare (Fig. 1). As a result, the genome sequences of H. annosum showed similar genomic characteristics with those of H. irregulare (Table 1). Protein clustering results showed that 6719 out of 8647 clusters were shared between H. annosum and H. irregulare. The 6719 clusters encompassed 10,701 and 12,193 proteins from H. annosum and H. irregulare, respectively, which account for 93.43% and 90.96% of the whole protein sequences. Among the shared clusters, 4804 clusters contained a single protein member from each species. A large number of shared clusters clearly showed their close relatedness with each other. Their close connection was also shown in the whole genome alignment. It showed that 84.38% of H. annosum genome sequences were aligned with H. irregulare counterparts. In addition, the aligned regions showed 96.06% of sequence identity across the whole genome, on average (Fig. S1). Meanwhile, structural variations, including SNPs and short InDels, further confirmed the taxonomic differences between the two species. A total of 188,190 SNPs were found when the sequencing reads of H. annosum were mapped against the genome sequences of H. irregulare. The majority of SNPs were found in genic regions (145,640), including exons and introns and untranscribed regions. In spite of macro level similarity between the two species, approximately six SNPs were found in every 1 kb, implying that there were great extent of genome modifications in micro level. In addition, 98.68% of the whole variations were found to be homo-variant, suggesting that there was long enough time for those variations to be fixed. Effects of SNPs and InDels were also pronounced in genic regions, followed by upstream, downstream and intergenic regions (Fig. 2a). Among the SNPs, transitions were far more frequently found than transversions, showing the ratio of transitions to transversions to be 3.57. When investigated the dinucleotide changes, the ratio of two transitions, AT to GC and GC to AT, accounted for 73.90%, on average. Each of the four transversions only accounted for 5.00 to 10.10% of the corresponding types (Fig. 2b). Though 92.11% of InDels showed their length within 10 bp, longer ones, up to 42 bp, were also detected (Fig. 2c). It might also suggest that the two species have gone through long-standing active mutations. This is in the line with the earlier report that the divergence of H. annosum and H. irregulare occurred around 30 mya [6]. In addition, a genome-wide association study by sequencing 23 haploid H. annosum s.s. in low-depth sequencing showed that there were approximately 33,000 SNPs for all individuals [8]. Considering that it was in low coverage, it is in agreement with our results in a way that there are many SNPs to convince they have diverged a long time ago.

Fig. 1

Genome comparisons of the two Heterobasidion spp. (a) Genomic characteristics of H. annosum are displayed, including density of repetitive elements, predicted genes, orthologous genes found in H. irregulare, scaffold representation, and genes involved in degradation of plant materials (aryl-alcohol dehydrogenase in gray, aryl-alcohol oxidase in black, CAZyme predicted by dbCAN in red, laccase in green, peroxidase in blue and quinone reductase in purple). (b) Genomic characteristics of H. irregulare genome and comparisons to H. annosum genome are presented, including SNP density (purple as the lowest and red as the highest), frequency of insertions/deletions (shown in green and red, respectively), density of repetitive elements, predicted genes, orthologous genes found in H. annosum, scaffold representation, and genes involved in degradation of plant materials with the same color codes used in Fig. 1a. The tracks were described from the innermost circle.

Table 1

Summary of H. annosum 03012 genome in comparison with that of H. irregulare TC 32-1.

	H. annosum 03012	H. irregulare TC 32-1
Assembled genome size	31.01 Mb	33.65 Mb
Number of contigs/scaffolds	2415	15
Number of ambiguous bases (Ns)	1574 bp	323 bp
GC contents	51.55%	52.23%
Number of predicted genes	11,453	13,405
Mean gene length (in amino acids)	405.8 aa	377.8 aa
Average number of exons per gene	5.64	5.39

Fig. 2

Summary of single nucleotide changes. (a) Distribution of SNPs by genetic positions, (b) ratio of mutational changes for across genetic positions, and (c) length distribution of insertions and deletions.

H. annosum-specific genes and unique genomic regions in H. irregulare

From the protein clustering result of the two species, we found 735 and 1624 species-specific genes in 616 and 1312 clusters of H. annosum and H. irregulare, respectively. Only 40 of 735 (5.44%) protein sequences, belonging to H. annosum-specific clusters, had protein domains predicted by InterPro scan. It showed huge difference to the fact that 58.38% of whole protein sequences in H. annosum had at least one domain. We hypothesized that species-specific sequences may be more unstructured than the rest of the proteome, since they presumably have emerged in comparatively recent date. In fact, 25.79% and 30.52% of amino acids in the H. annosum and H. irregulare-specific protein sequences were predicted to be global and local disorder, respectively. The rest of protein sequences in H. annosum, however, showed that 15.33% and 18.33% of amino acid residues were predicted to be locally and globally unstructured, respectively. It was also true for H. irregulare-specific proteins. Among the H. irregulare-specific proteins, 32.58% and 39.38% of amino acid residues were predicted as disordered, while 15.81% and 19.08% were marked as disordered in the others. It supports the fact that the proteins coded by species-specific genes tend to be unstructured than those conserved throughout other species [15]. It may address the observation that H. annosum-specific proteins had much less predicted domains compared to the whole proteome.

In order to see distribution of the H. annosum and H. irregulare-specific genes in other species, BLAST search was performed against the 37 rot fungi including the two Heterobasidion spp. A total of 468 and 806 of H. annosum and H. irregulare-specific protein sequences, respectively, did not show any hits in the other rot fungi genomes. In addition, 464 out of the 468 H. annosum sequences did not have any significant hits when searched against the NCBI non-redundant (NR) database, suggesting that orphan genes may encode them. Interestingly, four of the 468 protein sequences only had hits from bacterial species in the BLAST search, of which the hits were found in mainly sphingomonads belonging to the genera Sphingobium and Novosphingobium (Table S1). These bacterial species have been reported to have similar natural wood habitat as H. annosum [16]. Moreover, some isolates of these species are capable of degrading lignin or aromatic compounds like lignin-derived monoaryls and biaryls [17], [18]. Since these genes were not found in H. irregulare, it is speculated that genes encoding these proteins might be acquired from cohabiting sphingomonad bacteria. On the other hand, there were H. irregulare-specific regions retrieved from the dot plot (Fig. S1), which are 10 kb or longer. A total 1.93 Mb covering 35 regions showed that the repetitive sequences were enriched in comparison with the rest of H. irregulare genome sequences (Table S2). These regions encompassed 633 genes, of which 278 genes did not show any significant hits in any species other than H. irregulare from BLASTP searches against the NR database. Furthermore, 190 out of 278 genes were not found in H. annosum (Fig. S2). The results suggest that species-specific evolution is differently occurring in the two species for adapting to their own niches.

Evolutionary implication of gene synteny between the two Heterobasidion spp.

In order to speculate potential driving force underpinning the speciation of the two Heterobasidion spp., macroscopic and microscopic comparisons were conducted as shown above. The macroscopic comparison showed that the majority of contigs of H. annosum genome were mapped against the scaffolds in H. irregulare. According to SNP/InDel identification result, which is a microscopic level of genome comparison, six SNPs were identified in every 1 kb, on average. These comparisons showed that there might have been prolonged genomic rearrangements for adapting to their own niches. To further compare the two genomes at an intermediate scale, genome-wide gene synteny was investigated as well as non-syntenic regions in between. As a result, a total of 528 sets of syntenic genes was identified, encompassing 6894 genes (Table S3). Many sets of genes had small non-syntenic regions, which could be caused by different gene models predicted in two genome sequences. As a result, insertions/deletions of a few kb could be frequently found. For example, around 6 kb of insertion in H. annosum contained five predicted genes (Fig. 3). Interestingly, three of them had the best hit against the same gene in H. irregulare, a cytochrome P450 protein-encoding gene, implying that the inserted region might be caused by duplication. Another example showed that an insertion harboring four genes was found in H. irregulare (Fig. 4). The synteny analysis showed that there are many subtle differences at small and intermediate scale between the two genomes of Heterobasidion spp. A series of insertions/deletions of a few kb might have contributed to the adaptation for each species at different geographic locations.

Fig. 3

An example of non-syntenic regions found in H. annosum. (a) A dot plot presenting the alignment between H. annosum contig98 (64,767 bp) on X-axis and the corresponding region of H. irregulare on Y-axis. (b) Closer outlook of genomic context in H. annosum where an insertion containing the five predicted genes might occurred. (c) Genomic context showing the corresponding region in H. irregulare.

Fig. 4

An example of non-syntenic regions found in H. irregulare. (a) A dot plot presenting the alignment between H. irregulare scaffold_08:340,000–413,000 on X-axis and H. annosum contig120 on Y-axis. (b) Overview of genomic context in H. irregulare where an insertion containing the four predicted genes might occurred. (c) Genomic context showing the corresponding region in H. annosum.

Comparative gene family analysis of 37 rot fungi including Heterobasidion spp.

Many wood-decaying fungi usually have the ability to deconstruct plant materials, including cellulose, hemicellulose and lignin. In particular, white-rot fungi are well known for their ability to delignify wood materials [19]. Enzymes involved in lignin degradation, such as peroxidases and laccases, have been highlighted for their ecological significance and potential applications in industry [20], [21]. Accordingly, the identification of gene families involved in plant material degradation was conducted with 37 genome sequences of rot fungi (Table S4). For the two Heterobasidion spp., the results of identified gene families closely resembled with each other, since their genomes showed high similarity and they cause the same disease. Interestingly, CAZymes, aryl-alcohol oxidases (AAOs), peroxidases and laccases were more frequently found in white-rot fungi than in brown-rot with statistical significance (t-test; P < 0.0078) (Fig. S3). Considering the fact that AAOs, peroxidases and laccases are known to be major players in lignin degradation [22], it might suggest that white-rot fungi would have better performance in lignin degradation.

On the other hand, secretory proteins play important roles in plant-microbe interactions, including colonization of host cells, defense against oxidative stress from plant cells, environmental sensing and acquisition of nutrients [23]. When secretory potential was assessed for the whole protein sequences for each species, white-rot fungi showed higher number of putative secretory proteins than brown-rot (t-test; P = 0.0001). This result may possibly be interpreted that fungi encoding more proteins involved in plant material degradation need higher secretory potential to convey them to the appropriate subcellular locations. In H. annosum, for example, considerable number of laccases, peroxidases, AAOs, and CAZymes were predicted to be secretory (Table S4). However, it has been reported that some basidiomycetes do not show the characteristics of either typical white- or brown-rot [24], [25]. It may underlie the plasticity in fungal life styles and imply that classical white and brown-rot separation has a room for improvement. Availability of more sequenced genomes of wood-decaying fungi may help to elucidate the complexity behind this classification.

Materials and methods

Fungal source and preparation of genomic DNA

The strain used for genome sequencing was Heterobasidion annosum P-type (isolate 03012 provided by Kari Korhonen, METLA, Finland). The strain was grown on malt extract medium. Genomic DNA from H. annosum was prepared by following a published method [26]. Mycelium was grinded in liquid nitrogen, and pulverized and incubated in equal volumes of lysis buffer (0.2 M sodium borate, 30 mM EDTA, 1% SDS, pH 9.0 and phenol) at 60 °C for 5 min. After centrifugation, the supernatant was treated with RNase, and afterwards with an equal volume phenol/chloroform (1:1). Genomic DNA was precipitated with 0.7 volume 2-propanol, washed with 70% ethanol, dried, and resuspended in 100–200 μL of double distilled H2O.

Genome sequencing, assembly and gene prediction

Whole genome shotgun strategy was used in generation of short reads for H. annosum 03012 genome by using Illumina Genome Analyser IIx platform with the insert size of 500 bp. The short reads were assembled by SOAP package including SOAPdenovo (version 2.04) [27] and GapCloser (version 1.12). The generated contigs were further used in generation of scaffolds by SSPACE (version 2.0) [28]. Gene prediction was performed by using MAKER (version 2.28) [29] pipeline, integrating the results from SNAP (version 2006-07-28) [30], AUGUSTUS (version 2.7) [31] and GeneMark-ES (version 2.3e) [32].

Preparation of genome sequences of selected rot fungi

Genome sequences of 36 rot fungi were downloaded from MycoCosm at U.S. Department of Energy Joint Genome Institute website [14]. The data set includes 13 brown-rot and 23 white-rot fungi (Table S5). All the genome data was also archived in the standardized data warehouse of Comparative Fungal Genomics Platform 2.0 (CFGP 2.0; http://cfgp.snu.ac.kr/) for further analysis [33].

Whole genome alignment and detection of structural variations

In order to see genome-wide conservation, contig sequences of H. annosum were aligned against each of the scaffold sequences of H. irregulare by using BLAT [34]. By taking the hits showing ≥ 1 kb matches into consideration, dot plots for each H. irregulare scaffolds were created by using Gepard (version 1.30) [35]. The best hit for each H. annosum contig against the H. irregulare genome was collected for calculation of the net aligned regions between the two genomes. Single nucleotide polymorphisms (SNPs), short insertions and deletions were identified and analyzed by using BWA (version 0.5.9-r16) [36], SAMtools (version 0.1.18) [37] and snpEff (version 3.4i) [38]. The reads were trimmed from 3′ termini if base quality is < 20 in the Phred scale. The SNPs supported by < 4 reads or > 100 reads were discarded for under- and over-representation. The size of upstream and downstream for each transcript were set to be 1 kb each.

Identification of gene synteny between the Heterobasidion spp.

In order to identify syntenic regions between H. annosum and H. irregulare genomes, candidate orthologous gene pairs were defined by using BLASTP [39]. Protein sequences in H. annosum were used as query and searched against those in H. irregulare, followed by filtering the best hits. Syntenic regions were identified by using DAGchainer [40] with default settings, and the micro-syntenic differences were shown with the selected regions.

Gene family annotation

Previously published pipelines were used in identification of putative genes encoding plant cell wall-degrading enzymes (PCWDEs), secretory proteins and peroxidases [23], [33], [41], [42]. To capture the maximum extent of PCWDEs, carbohydrate-active enzymes (CAZymes) were identified by using dbCAN [43] with the E-value cutoff of 1e − 5. And an in-house pipeline was used in laccase and laccase-like multi‑copper oxidase prediction (unpublished). Additional genes involved in lignin degradation, including aryl-alcohol dehydrogenase (AAD), aryl-alcohol oxidase (AAO) and quinone reductase (QR), were identified based on protein domain profiles predicted by InterPro scan [44] (Table S6) The domain profiles were retrieved by using the protein sequences annotated as AAD, AAO and QR in UniProt database [45].

Protein clustering analysis

Whole proteome sequences were subjected to Markov clustering using tribeMCL [46]. All-by-all BLASTP was performed to prepare input data for clustering analysis with the E-value cutoff of 1e − 5 and the inflation factor of 2.0.

Prediction of unstructured protein regions

Whole proteome sequences of the two Heterobasidion spp. were scanned by using IUPred, a software predicting global and short disorder in proteins sequences [47]. Subsequently amino acid residues showing disorder tendency above 0.5 were considered as disordered. In addition, residues predicted as structured were removed from the list of local and global disorder.

Data deposition

The draft genome sequence of Heterobasidion annosum s.s. isolate 03012 was deposited at DDBJ/EMBL/GenBank under accession no. AOSL00000000.

The following are the supplementary data related to this article.

Supplementary material

Table S3

List of the syntenic gene sets between the two genomes of Heterobasidion spp.

Transparency document

Transparency document.