Research article | Open | Published:
Comparative genomics of Coniophora olivacea reveals different patterns of genome expansion in Boletales
BMC Genomicsvolume 18, Article number: 883 (2017)
Coniophora olivacea is a basidiomycete fungus belonging to the order Boletales that produces brown-rot decay on dead wood of conifers. The Boletales order comprises a diverse group of species including saprotrophs and ectomycorrhizal fungi that show important differences in genome size.
In this study we report the 39.07-megabase (Mb) draft genome assembly and annotation of C. olivacea. A total of 14,928 genes were annotated, including 470 putatively secreted proteins enriched in functions involved in lignocellulose degradation. Using similarity clustering and protein structure prediction we identified a new family of 10 putative lytic polysaccharide monooxygenase genes. This family is conserved in basidiomycota and lacks of previous functional annotation. Further analyses showed that C. olivacea has a low repetitive genome, with 2.91% of repeats and a restrained content of transposable elements (TEs). The annotation of TEs in four related Boletales yielded important differences in repeat content, ranging from 3.94 to 41.17% of the genome size. The distribution of insertion ages of LTR-retrotransposons showed that differential expansions of these repetitive elements have shaped the genome architecture of Boletales over the last 60 million years.
Coniophora olivacea has a small, compact genome that shows macrosynteny with Coniophora puteana. The functional annotation revealed the enzymatic signature of a canonical brown-rot. The annotation and comparative genomics of transposable elements uncovered their particular contraction in the Coniophora genera, highlighting their role in the differential genome expansions found in Boletales species.
Coniophora olivacea is a basidiomycete fungus belonging to the order Boletales. C. olivacea produces brown-rot decay on dead wood of conifers (softwood) and, less frequently, on hardwood species. In addition, C. olivacea also damages wood buildings or construction materials. The genome sequence of its sister species C. puteana was made public in 2012  and contributed to the understanding of genomic differences between brown and white-rot fungi. White-rot fungi are efficient lignin degraders, whereas brown-rot fungi attack cell wall carbohydrates leaving lignin undigested. The main responsible of this behavior are lignin-degrader peroxidases, which are abundant in white-rot species and particularly contracted in brown-rot and mycorrhizal fungi . The Boletales order comprises a diverse group of species including saprotrophs and ectomycorrhizal species such as Suillus sp. or Pisolithus sp. During the last 6 years, up to 12 Boletales genomes have been sequenced and annotated [1, 3, 4]. Information that emerged from these studies showed important differences in genomic characteristics between the species belonging to this group, whose predicted common ancestor was dated 84 million years ago. Evolution from this boletales ancestor (supposed to be a brown-rot saprotroph) lead to the diversification and the appearance of ectomycorrhizae, which shows a particular contraction of the number of plant cell wall-degrading enzymes coding genes (PCWDE) [4, 5]. In addition, Boletales show important differences in their genome size and gene content. For example, the smallest assembled Boletales genome spans 38.2 Mb and has 13,270 annotated genes (Hydnomerulius pinastri), but the largest (Pisolithus tinctorius) spans 71.0 Mb and has 22,701 genes . Previous studies in saprophytic basidiomycetes have shown that species with higher genome sizes tend to have more transposable elements . Also, it has been described that species associated with plants (pathogenic and symbiotic) have genomes with expanded TE families [1, 7], although this trend varies between the three basidiomycete phyla . In this paper, we describe the draft genome sequence and annotation of the brown-rot C. olivacea, and we compare it with the genomes of C. puteana as well as with that of three other Boletales showing important differences in genome sizes (Serpula lacrymans, Pisolithus tinctorius and Hydnomerulius pinastri). The results show that C. olivacea displays enzymatic machinery characteristic of brown-rot fungi encoded in a compact genome, carrying a small number of repetitive sequences. The comparative analysis with other Boletales shows that both ancient and modern LTR-retrotransposon amplification events have greatly contributed to the genome expansion along the evolution of Boletales.
Fungal strains and culture conditions
Coniophora olivacea MUCL 20566 was obtained from the Spanish Type Culture Collection and was cultured in SMY submerged fermentation (10 g of sucrose, 10 g of malt extract and 4 g of yeast extract per litre).
Nucleic acid extraction
Mycelia were harvested, frozen, and ground in a sterile mortar in the presence of liquid nitrogen. High molecular weight DNA was extracted using the phenol-chloroform protocol described previously . DNA sample concentrations were measured using a Qubit® 2.0 Fluorometer (Life Technologies, Madrid, Spain), and DNA purity was measured using a NanoDrop™ 2000 (Thermo-Scientific, Wilmington, DE, USA). DNA quality was verified by electrophoresis in 0.7% agarose gels. Total RNA was extracted from 200 mg of deep-frozen tissue using Fungal RNA E.Z.N.A Kit (Omega Bio-Tek, Norcross, GA, USA), and its integrity was verified using the Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA).
Genome and transcriptome sequencing and assembly
A detailed description is provided in Additional file 1: Text S1. Briefly, the C. olivacea MUCL 20566 genome was sequenced using Illumina HiSeq-1 TB Regular 2 × 151 bp 0.309 kb. Sequenced reads were QC filtered for artifact contamination using BBDuk from the BBMap package (https://sourceforge.net/projects/bbmap/) and subsequently assembled with Velvet 1.2.07 . The result -pair library with an insert size of 3000 +/− 300 bp in silico that was then assembled together with the original Illumina library with AllPathsLG . Raw sequences were deposited in SRA (Sequence Read Archive) NCBI database under accession number SRP086489. Strand-specific RNASeq libraries were created and quantified by qPCR. Sequencing was performed using an Illumina HiSeq-2500 instrument. Reads were filtered and trimmed to remove artifacts and low quality regions using BBDuk. Transcriptome was de novo assembled using Trinity  and used to assist annotation and assess the completeness of the corresponding genome assembly using alignments of at least 90% identity and 85% coverage.
The genome assemblies of C. olivacea MUCL 20566 and C. puteana (http://genome.jgi.doe.gov/Conpu1/Conpu1.home.html) were aligned using the Promer tool from the MUMmer 3.0 package . Genome rearrangements were identified in the alignment with dnadiff tool from the same package.
The annotation of the C. olivacea MUCL 20566 assembly was performed using the Joint Genome Institute pipeline  to predict and functionally annotate protein-coding genes and other features such as tRNAs or putative microRNA precursors. The SECRETOOL pipeline  was used to identify putatively secreted proteins, considering the presence of signal peptides, cleavage sites, transmembrane domains and the GPI (glycosylphosphatidylinositol) membrane anchor. Carbohydrate-active enzymes (CAZys) were annotated based on BLAST  and HMMER  searches against sequence libraries and HMM (Hidden Markov Models) profiles of the CAZy database  functional modules. Protein structure predictions were carried out with Phyre2 . Raw sequencing reads, genome assembly, transcriptome assembly, gene predictions and functional annotations are publicly available in the C. olivacea genome portal of Mycocosm database (http://genome.jgi.doe.gov/Conol1/Conol1.home.html).
Annotation of transposable elements
Transposable elements (TEs) were identified and annotated in the C. olivacea assembly using REPET package [20, 21], as well as in the following boletales assemblies available in Mycocosm database (http://genome.jgi.doe.gov/programs/fungi/index.jsf): Coniophora puteana v1.0 (ID: Conpu1), Hydnomerulius pinastri v2.0 (ID: Hydpi2), Serpula lacrymans S7.3 v2.0 (ID: SerlaS7_3_2), Pisolithus tinctorius Marx 270 v1.0 (ID: Pisti1). Briefly, de novo TE detection was carried out with the TEdenovo pipeline  and the elements were classified with PASTEC . The resulting TE library was fed into TEannot pipeline  in two consecutive iterations: the first one with the full library, and the second with an improved library consisting on consensus elements carrying at least one full-length copy after manually discarding false positives (i.e., C. olivacea genes).
Insertion age of LTR-retrotransposons
Full-length LTR-retrotransposons were identified using LTRharvest  followed by BLASTX against Repbase . Long Terminal Repeats were extracted and aligned with MUSCLE . Alignments were trimmed using trimAl  and used to calculate Kimura’s 2P distances. The insertion age was calculated following the approach described in  using the fungal substitution rate of 1.05 × 10−9 nucleotides per site per year [6, 28].
Identification of gene families
All-by-all BLASTP followed by MCL (Markov Cluster Algorithm) clustering  was carried out with C. olivacea protein models using a threshold value of e−5 and an inflation value of 2. We considered gene families carrying four or more genes for further analyses.
The predicted proteomes of the following species were downloaded from Mycocosm database (Mycocosm ID in parenthesis):
Agaricus bisporus var. bisporus H97 v2.0 (Agabi_varbisH97_2), Boletus edulis v1.0 (Boled1), Coniophora olivacea MUCL 20566 v1.0 (Conol1), Coniophora puteana v1.0 (Conpu1), Cryptococcus neoformans var. grubii H99 (Cryne_H99_1), Fomitopsis pinicola FP-58527 SS1 v3.0 (Fompi3), Gyrodon lividus BX v1.0 (Gyrli1), Hydnomerulius pinastri v2.0 (Hydpi2), Leucogyrophana mollusca KUC20120723A-06 v1.0 (Leumo1), Paxillus involutus ATCC 200175 v1.0 (Paxin1), Phanerochaete chrysosporium RP-78 v2.2 (Phchr2), Pisolithus tinctorius Marx 270 v1.0 (Pisti1), Pleurotus ostreatus PC15 v2.0 (PleosPC15_2), Rhizopogon vinicolor AM-OR11–026 v1.0 (Rhivi1), Scleroderma citrinum Foug A v1.0 (Sclci1), Serpula lacrymans S7.3 v2.0 (SerlaS7_3_2), Suillus luteus UH-Slu-Lm8-n1 v2.0 (Suilu3), Trametes versicolor v1.0 (Trave1). Species phylogeny was constructed as follows: all-by-all BLASTP followed by MCL clustering was carried out with a dataset containing the proteomes of all the species. The clusters carrying only one protein per species were identified, and the proteins were aligned using MAFFT . The alignments were concatenated after discarding poorly aligned positions with Gblocks . The phylogeny was constructed using RaxML  with 100 rapid bootstraps under PROTGAMMAWAGF substitution model. Phylogenetic reconstruction of Gypsy reverse-transcriptases was carried out as follows: Reverse transcriptase RV1 domains were extracted from LTR-retrotransposons of the TE consensus library using Exonerate  and aligned with MUSCLE. The alignments were trimmed using trimAl with the default parameters, and an approximate maximum likelihood tree was constructed using FastTree .
C. olivacea assembly and annotation
The nuclear genome of C. olivacea was sequenced with 137 X coverage and assembled into 863 scaffolds accounting for 39.07 Mb, 90.3% of the genome size estimation based on k-mer spectrum (43.28 Mb). The mitochondrial genome was assembled into two contigs accounting for 78.54 kb. The assembly completeness was 99.78% according to the Core Eukaryotic Genes Mapping Approach (CEGMA ), with only one missing accession (KOG1322, GDP-mannose pyrophosphorylase). We assembled 66,567 transcripts (mean lenght = 2,744 nt, median = 2,154 nt) of which 97.8% could be mapped to the genome. The C. olivacea assembled genome was more fragmented than its close relative C. puteana (Table 1). The total repeat content was 2.91% of which 2.15% corresponded to transposable elements, 0.64% to simple repeats, and 0.12% to low complexity regions. The estimation of repeat content from low-coverage Illumina data (3.8X) yielded 6% of the genome size covered by transposable elements (Additional file 2: Table S1). We used transcriptomic information, ab initio predictions and similarity searches to predict a total of 14,928 genes—84.5% of them having a strong transcriptome support (spanning more than 75% of the gene length). In addition, 88.3% of the annotated genes had significant similarity to proteins from the NCBI nr database and 46.6% to the manually curated proteins from the Swiss-Prot database (cutoff e−05) . A total of 7,841 predicted proteins (52.3%) carried Pfam domains and 1,471 (9.8%) carried signal peptide, of which 470 were predicted to be secreted using the more stringent SECRETOOL pipeline.
The multigene phylogeny based on 1,677 conserved single copy genes displayed different classes, orders and families in branches congruent with previous phylogenetic data  and with very high support. C. olivacea was placed in a branch next to its sequenced closer species C. puteana representing the Coniophoraceae family in the order Boletales (Fig. 1).
The whole-genome protein-based alignment between the two Coniophoraceae species spanned 52.7% of the C. olivacea and 48.0% of C. puteana assemblies. It shows evidence of macrosynteny between the two species (Fig. 2a, Additional file 3: Fig. S1), with an average similarity of 78.4% in the aligned regions (Fig. 2b) and numerous inversions (1,027 regions). The good conservation between both genomes in protein coding regions was evidenced by the amount of orthologous genes obtained using the reciprocal best hit approach (7,468 genes with more than 70% identity over 50% of protein sequences) and by the number of C. olivacea proteins yielding significant tBLASTN hits against the C. puteana genome (13,572 genes, cutoff e-5, Fig. 2c). For the remaining 1,352 C. olivacea-specific (orphan) genes, only 48 could be functionally annotated based on KOG (Eukaryotic Orthologous Groups), KEGG (Kyoto Encyclopedia of Genes and Genomes), GO (Gene Ontology) or InterPro databases.
Carbohydrate-active enzymes of C. olivacea
The annotated proteome was screened for the presence of carbohydrate-active enzymes (CAZy). A total of 397 proteins were annotated and classified into different CAZy classes and associated modules. The CAZyme profile of C. olivacea was very similar to that of C. puteana although small differences were found in the glycoside hydrolases (GH, Additional file 4: Table S2). Some families such as GH5, GH18 or GH31 were smaller than in C. puteana. Similar to other brown-rot basidiomycetes, C. olivacea lacked Class II peroxidases (Auxiliar Activities AA2) and displayed a reduced set of other cellulolytic enzymes such GH6 (1), GH7 (1) and CBM1 (2) and AA9 (6).
Functional characteristics of C. olivacea predicted secretome
Using SECRETOOL pipeline we predicted 470 putatively secreted proteins in C. olivacea and 504 in C. puteana. An enrichment analysis of gene ontology (GO) terms was performed to determine what gene functions were over-represented in the secreted proteins. Thirty GO terms were significantly enriched including 24 corresponding to molecular functions, four to biological processes and two to cellular components (Table 2). The most enriched molecular function was “feruloyl esterase activity,” which is responsible for plant cell-wall degradation. “Polysaccharide catabolic process” was the most enriched GO term within the biological processes, and “extracellular region” within the cellular components (Table 2).
Analysis of putatively secreted multigene families
Using all-by-all BLASTP followed by MCL we clustered by similarity the 1,471 proteins carrying signal peptides in C. olivacea. We used all proteins carrying signal peptides rather than only SECRETOOL predictions in order to obtain larger protein clusters. Up to 60% of the 1,471 proteins grouped in clusters were formed by 2 to 59 genes (Additional file 5: Table S3), showing the same distribution as the whole proteome (p = 0.6032, Wilcoxon test, 61% of the 14,928 predicted genes were found in clusters containing 2 to 157 members). For further analysis of the secreted genes found in clusters, we focused on the 70 clusters (families) formed by four or more gene members. Using the KOG, KEGG, InterPro and GO databases, we could assign functions to 45 out of the 70 gene families (Table 3). Cytochrome P450, hydrophobins and aspartic-peptidases were the largest gene families. In addition, 17 CAZys clusters were found including glycoside hydrolases (GH), carbohydrate esterases (CE), carbohydrate-binding modules (CBMs) and redox enzymes classified as auxiliary activities (AA). 25 clusters lacked functional annotation, and some of them had a high number of genes (clusters 2, 6 and 7 in Table 3). All of these genes belonging to families with unknown function were further analyzed with Phyre2 to predict their protein structure and used for PSI-BLAST (Position-Specific Iterated BLAST) analysis. Using this approach, two gene families were functionally annotated with high confidence (96.3–97.4% confidence for individual protein predictions): one as a copper-dependent lytic polysaccharide monooxygenase (LPMO, also known as AA9; cluster 16), and the other as thaumatin-lyke xylanase inhibitor (tlxi, cluster 48). The Cluster16 containing putative LPMOs was particularly interesting. This was formed by 10 genes coding for small proteins ranging from 130 to 162 amino acids with three exons (with the exception of protein ID839457 that shows only two). All these genes coded for proteins that have a signal peptide but lack of known conserved functional domains. Six were confidently annotated as LPMOs by Phyre2, and four of them were predicted to be secreted by SECRETOOL. In addition, this family of unknown proteins is conserved in all the agaricomycetes shown in Fig. 1. Interestingly, four members of this family appear as a tandem located in C. olivacea scaffold_124 (scaffold_426:4800–12,000).
Impact of repeat content on C. olivacea genome size and other Boletales
To study the role that TEs have played in the evolution of the Boletales genomes, we annotated and quantified the TE content in five species showing important differences in genome size: C. olivacea (39.1 Mb), C. puteana (42.9 Mb) , Hydnomerulius pinastri (38.2 Mb) , Serpula lacrymans (47.0 Mb)  and Pisolithus tinctorius (71.0 Mb)  (Additional file 6: Dataset S1, Additional file 7: Dataset S2, Additional file 8: Dataset S3, Additional file 9: Dataset S4, Additional file 10: Dataset S5). TEs were de novo identified and annotated using pipelines of the REPET package. The results yielded major differences in TE content between the five species, with C. olivacea, C. puteana and H. pinastri having low TE content (2.15%, 3.94% and 6.54% of their corresponding genome sizes), and S. lacrymans and P. tinctorius having up to 29.45% and 41.17% of their genomes occupied by TEs, respectively (Fig. 3, Table 4). In addition to higher TE content, species with larger genome assembly size showed higher TE diversity as reflected by the higher number of TE families, which ranged between 43 in C. olivacea to 432 in P. tinctorius.
The TEs found belong to seven out of the nine TE orders described by Wicker et al : LTR, DIRS (Dictyostelium Intermediate Repeat Sequences), PLE (Penelope-like Elements), LINE (Long Interspersed Nuclear Elements), SINE (Small Interspersed Nuclear Elements), TIR (Terminal Inverted Repeats) and Helitrons. Two of the orders (LTR and TIRS, which contain long terminal repeats or terminal inverted repeats, respectively) were present in the five species. Class I TEs were primarily responsible for the observed genome size differences—especially the elements belonging to LTR in the Gypsy superfamily, which accounted for more than 15% of the assembly in S. lacrymans and P. tinctorius, but less than 3% in H. pinastri, C. olivacea and C. puteana. Of all the LTR/Gypsy families detected by TEdenovo, we observed that those elements belonging to the Chromoviridae group (carrying a Chromatin organization domain, PF00385, in the N-terminal region after the integrase, Fig. 4) were the most abundant LTR-retrotransposons in these five species, ranging from 44 to 83% of the total Gypsy coverage. LTR-retrotransposons in the Copia superfamily were also particularly abundant in S. lacrymans and P. tinctorius (accounting for 2.4–6% of the total assembly size). Remarkably, non-coding LTR-retrotransposons such as TRIM (Terminal-repeat Retrotransposons In Miniature) and LARD (Large Retrotransposon Derivatives) were also found in three out of the five genomes, but in lower amounts (<1% of the genome, Table 4).
LINE, SINE, DIRS and PLE elements were also found in low copy numbers, but none of these were present in the five species. Regarding Class II transposons, TIR order was the most important in terms of abundance and copy number with elements encoding DDE transposases present in the five species. The second most important were MITEs (Miniature Inverted–repeat Transposable Elements) and other non-coding elements carrying structural features (classified as TIR/unknown in Table 1). Rolling-circle helitrons were found in H. pinastri, S. lacrymans and P. tinctorius, while putative Mavericks were present only in this latter one.
Phylogenetic reconstruction of the LTR reverse-transcriptases
To understand the phylogenetic relationship between the LTR-retrotransposon familes in the five analyzed genomes, we inferred a maximum likelihood phylogeny of the LTR reverse-transcriptases of the Gypsy consensus sequences (Fig. 5). Three main clades were obtained (A, B and C). Clades A and B were formed, almost exclusively, by families found in the P. tinctorius genome. Moreover, while clade B is formed mostly by distantly related families, the profile of clade A suggests that an important fraction of the families underwent recent diversification. All LTR families found in the other four species grouped in clade C along with the remaining families of P. tinctorius. This clade contained several retrotransposon sub-clades sharing closely related families from three to five species.
Age of the LTR-retrotransposon amplification bursts in the Boletales
LTR-retrotransposons carrying conserved domains as well as intact Long Terminal Repeats (putative autonomous elements) were subjected to further study to investigate their amplification dynamics over the course of evolution. Based on the nucleotide divergence between the two LTRs, we estimated the time of insertion of each element using a substitution rate of 1.05 × 10−9 nucleotide substitutions per site per year. The number of intact, putative autonomous LTR-retrotransposons varied greatly in the five species ranging from 26 elements in C. olivacea to 944 in P. tinctorius. The LTR profiles of C. olivacea, C. puteana and S. lacrymans showed recent peaks of amplification with insertion dates at 0–5 million years (MY). LTR amplification in H. pinastri showed a peak at 10–15 MY ago, whereas the profile of P. tinctorium pointed to a much older amplification burst showing a maximum peak at 25–30 MY ago and few recent retrotransposition events (Fig. 6).
Genomic and proteomic characteristics of C. olivacea
We report the 39.07 Mb draft genome assembly and annotation of brown-rot basidiomycete C. olivacea. In terms of genome size, this species is slightly smaller than C. puteana, but it falls in the range of other brown-rot basidiomycetes such as Hydnomerulius pinastri (38.3 Mb)  or Serpuyla lacrymans (47.0 Mb). As expected for closely related species, C. olivacea and C. puteana show macrosynteny, although due to the short scaffold lengths it is impossible to establish comparisons at a chromosome scale. We found very good conservation of protein-coding genes, although C. olivacea has up to 1,352 orphan genes—most of these are supported by structure and RNA evidence (i.e., no homology to any other known gene). In this sense, the higher number of annotated genes in C. olivacea relative to C. puteana is probably related to the higher amount of assembled RNA contigs used to assist the annotation of the former (resulting from the higher RNAseq depth). The presence of about 10% of orphan genes is common in fungal genomes, and these genes often lack an in silico functional annotation as we found for C. olivacea [39, 40].
Wood-decaying species require a complex enzymatic machinery to degrade lignin and obtain nutrients. According to the CAZy enzymes identified in the genome, the C. olivacea proteome carries the main signatures of canonical brown-rot: (i) it completely lacks Class II peroxidases—enzymes primarily involved in lignin degradation , and (ii) it carries a reduced set of enzymes involved in degradation of crystalline cellulose. In fact, its profile is very similar to that of C. puteana, displaying only minor differences in several enzyme groups. As previously seen in other wood-degrading fungi, the in silico secretome of C. olivacea is enriched in functions related to lignocellulose degradation . Our analysis showed that most intracellular and secreted proteins are members of multi-gene families of diverse size originating from gene duplications. The number of gene families that could not be functionally annotated by standard similarity-based methods was high, a phenomenon that is frequently observed in fungi.
To overcome this drawback, we used an alternative approach that combines similarity with structural information (Phyre-2). We then assigned a putative function to two multi-gene families conserved across the basidiomycete phylogeny but for which a putative function had not been previously proposed. Of special interest is the newly identified family of putative copper-dependent lytic polysaccharide monooxygenases (AA9, LPMO). The LPMOs are recently discovered enzymes used by microbes to digest crystalline polysaccharides . They increase the saccharification yield of commercial enzyme cocktails . Nevertheless, despite the promising results obtained in silico, experimental assays will be necessary to confirm the function of the members of this newly described gene family.
Impact of TEs in the evolution of Boletales genomes
The results of TE annotation in the five Boletales showed how different patterns of LTR-retrotransposon amplifications have shaped the architecture of their genomes. The expansion of LTR/Gypsy retrotransposons belonging to Chromoviridae occurred mainly in the species with large genomes, whereas the smaller genomes have a small amount of these families (ie, three families in C. olivacea and C. puteana). Chromoviruses are the most common LTR-retrotransposons in fungi , and the key to their success might be the presence of a chromo-integrase, which is thought to guide the integration of these elements into heterochromatic regions . Heterochromatin is gene-poor, and it is silenced by epigenetic mechanisms such as DNA methylation and RNAi . Thus, integration of these elements in such regions would allow them to skip purifying selection and increase their probability to persist in the genome. In fact, this could be the reason for the longer prevalence of Gypsy over Copia LTR-retrotransposons in most fungal species—the latter tend to integrate at random locations including euchromatic regions where transposon fixation is more difficult . The LTR-retrotransposon amplification bursts of the Boletales indicate that elements from both Coniophora species are young and thus putatively active, and the profile of S. lacrymans also indicates a very strong activity of young copies with a progressive decrease in the amplification signals of older elements. Our findings suggest that the latter three species are currently in a period of genome expansion. Despite the different profile of H. pinastri and P. tinctorius we cannot rule out the same hypothesis, as both assemblies contain high gap content (7.7% and 13.3%, respectively). This fact usually leads to an underestimation in the amount of young retrotransposons , as they are difficult to assemble due to their repetitive nature and high sequence identity. In fact, we show that due to this reason the assembly-based TE quantification underestimated LTR content in C. olivacea in comparison to non-assembly based quantification (Additional file 2: Table S1). The profile of P. tinctorius is intriguing. This ectomycorrhizal (ECM) species undergoes a massive expansion of LTR-retrotransposons in the Gypsy superfamily (similar to that found for other symbiotic species in Agaricomycotina [7, 49]; however, the majority of elements are very old (20–40 MY) and still carry structural and coding domains necessary for transposition. The phylogeny of Gypsy reverse-transcriptases suggests that many P. tinctorius-specific families are distantly related to the other four species. In fact, its impressive retrotransposon content might be partially explained by the amplification and diversification of ancestral families (giving rise to clades A and B in Fig. 5). Our phylogenetic reconstruction suggests that such ancestral families were also present in other boletales but didn’t proliferate in the genome (ie, H. pinastri or C. puteana). Whether genome defense mechanisms or lifestyle constraints are responsible of this phenomenon is still to be demonstrated. In this regards, it is interesting to note that the LTR-mediated genome amplification of P. tinctorius roughly coincides with the estimated origins of ECM symbiosis in Boletales . Of the four Class I TE orders found, only the LTR elements were present in the five species. The most plausible scenario is that the elements from the other three orders (DIRS, LINE, and PLE) were lost by random drift in some of the species. Alternatively, they might be present in some genomes but in the form of very ancient and degenerated copies that are not detectable. Similarly, this patchy distribution was also found in class II elements (ie, helitrons were absent in the Coniophora genus and present in the remaining three species). Previous studies have shown that besides the conserved presence of LTR and TIR orders, the remaining TE groups tend to be present in variable amounts in basidiomycetes .
In this study we present the draft genome sequence and annotation of the brown-rot fungi Coniophora olivacea, along with a comparative analysis with C. puteana and other members of the Boletales order. Our results show evidence of macrosynteny and conservation in the protein coding genes of the two species. The functional analysis of C. olivacea secretome showed that it displays the main signatures of a canonical brown-rot, and uncovered a new family of putative LPMOs widely conserved in basidiomycota. The annotation of transposable elements revealed a particular contraction in these two species in comparison to other Boletales, mainly due to the differential expansion of Chromoviridae LTR-retrotransposons. By analyzing the distribution of insertion ages and phylogenetic relationships of these elements we show that these LTR-retrotransposons have played a key role in the genome expansion experienced by certain species in the Boletales order.
Core Eukaryotic Genes Mapping Approach
Dictyostelium intermediate repeat sequence
Hidden Markov Models
Kyoto Encyclopedia of Genes and Genomes
Eukaryotic Orthologous Groups
Large retrotransposon derivative
Long interspersed nuclear elements
Lytic polysaccharide monooxygenases
Long Terminal Repeats
Miniature inverted-repeat transposable elements
Plant cell wall-degrading enzymes
Reciprocal best hit
Sucrose, malt, yeast
Sequence Read Archive
Terminal inverted repeats
Terminal-repeat retrotransposon in miniature
Floudas D, Binder M, Riley R, Barry K, Blanchette RA, Henrissat B, et al. The paleozoic origin of enzymatic lignin decomposition reconstructed from 31 fungal genomes. Science. 2012;336:1715–9.
Riley R, Salamov AA, Brown DW, Nagy LG, Floudas D, Held BW, et al. Extensive sampling of basidiomycete genomes demonstrates inadequacy of the white-rot/brown-rot paradigm for wood decay fungi. Proc Natl Acad Sci U S A. 2014;111:9923–8.
Eastwood DC, Floudas D, Binder M, Majcherczyk A, Schneider P, Aerts A, et al. The plant cell wall-decomposing machinery underlies the functional diversity of forest fungi. Science. 2011;333:762–5.
Kohler A, Kuo A, Nagy LG, Morin E, Barry KW, Buscot F, et al. Convergent losses of decay mechanisms and rapid turnover of symbiosis genes in mycorrhizal mutualists. Nat Genet. 2015;47:410–5.
Martin F, Kohler A, Murat C, Veneault-Fourrey C, Hibbett DS. Unearthing the roots of ectomycorrhizal symbioses. Nat Rev Microbiol. 2016;14:760–73.
Castanera R, López-Varas L, Borgognone A, LaButti K, Lapidus A, Schmutz J, et al. Transposable elements versus the fungal genome: impact on whole-genome architecture and transcriptional profiles. PLoS Genet. 2016;12
Hess J, Skrede I, Wolfe BE, LaButti K, Ohm RA, Grigoriev IV, et al. Transposable element dynamics among asymbiotic and ectomycorrhizal Amanita fungi. Genome Biol Evol. 2014;6:1564–78.
Castanera R, Borgognone A, Pisabarro AG, Ramírez L. Biology, dynamics, and applications of transposable elements in basidiomycete fungi. Appl Microbiol Biotechnol. 2017:1–14.
Larraya LM, Perez G, Penas MM, Baars JJP, Mikosch TSP, Pisabarro AG, et al. Molecular Karyotype of the white rot fungus Pleurotus ostreatus. Appl Envir Microbiol. 1999;65:3413–7.
Zerbino DR, Birney E. Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18:821–9.
Gnerre S, MacCallum I, Przybylski D, Ribeiro FJ, Burton JN, Walker BJ, et al. High-quality draft assemblies of mammalian genomes from massively parallel sequence data. Proc Natl Acad Sci. 2010;108:1513–8.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.
Kurtz S, Phillippy A, Delcher AL, Smoot M, Shumway M, Antonescu C, et al. Versatile and open software for comparing large genomes. Genome Biol. 2004;5:R12.
Grigoriev IV, Nikitin R, Haridas S, Kuo A, Ohm R, Otillar R, et al. MycoCosm portal: gearing up for 1000 fungal genomes. Nucleic Acids Res. 2014;42:D699–704.
Cortázar AR, Aransay AM, Alfaro M, Oguiza JA, Lavín JLSECRETOOL. Integrated secretome analysis tool for fungi. Amino Acids. 2014;46:471–3.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
Johnson LS, Eddy SR, Portugaly E. Hidden Markov model speed heuristic and iterative HMM search procedure. BMC Bioinformatics. 2010;11:431.
Cantarel BL, Coutinho PM, Rancurel C, Bernard T, Lombard V, Henrissat B. The carbohydrate-active EnZymes databse (CAZY): an expert resource for Glycogenomics. Nucleic Acids Res. 2009;37(Database issue):D233–8.
Kelley LA, Mezulis S, Yates CM, Wass MN, Sternberg MJE. The Phyre2 web portal for protein modeling, prediction and analysis. Nat Protoc. 2015;10:845–58.
Quesneville H, Bergman CM, Andrieu O, Autard D, Nouaud D, Ashburner M, et al. Combined evidence annotation of transposable elements in genome sequences. PLoS Comput Biol. 2005;1:166–75.
Flutre T, Duprat E, Feuillet C, Quesneville H. Considering transposable element diversification in de novo annotation approaches. PLoS One. 2011;6:e16526.
Hoede C, Arnoux S, Moisset M, Chaumier T, Inizan O, Jamilloux V, et al. PASTEC: an automatic transposable element classification tool. PLoS One. 2014;9:e91929.
Ellinghaus D, Kurtz S, Willhoeft U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinformatics. 2008;9:18.
Jurka J. Repbase update - a database and an electronic journal of repetitive elements. Trends Genet. 2000;16:418–20.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.
Capella-Gutierrez S, Silla-Martinez JM, Gabaldon T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009;25:1972–3.
SanMiguel P, Gaut BS, Tikhonov A, Nakajima Y, Bennetzen JL. The paleontology of intergene retrotransposons of maize. Nat Genet. 1998;20:43–5.
Dhillon B, Gill N, Hamelin RC, Goodwin SB. The landscape of transposable elements in the finished genome of the fungal wheat pathogen Mycosphaerella graminicola. BMC Genomics. 2014;15:1132.
Enright AJ, Van Dongen S, Ouzounis CA. An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res. 2002;30:1575–84.
Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.
Talavera G, Castresana J. Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst Biol. 2007;56:564–77.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
Slater GS, Birney E. Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics. 2005;6:31.
Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol. 2009;26:1641–50.
Parra G, Bradnam K, Korf I. CEGMA: a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics. 2007;23:1061–7.
Bairoch A, Boeckmann B. The SWISS-PROT protein sequence data bank. Nucleic Acids Res. 1991;19(Suppl):2247–9.
Hibbett DS, Binder M, Bischoff JF, Blackwell M, Cannon PF, Eriksson OE, et al. A higher-level phylogenetic classification of the fungi. Mycol Res. 2007;111:509–47.
Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, et al. A unified classification system for eukaryotic transposable elements. Nat Rev Genet. 2007;8:973–82.
Grandaubert J, Bhattacharyya A, Stukenbrock EH. RNA-seq-based gene annotation and comparative genomics of four fungal grass pathogens in the genus Zymoseptoria identify novel orphan genes and species-specific invasions of transposable elements. G3 (Bethesda). 2015;5:1323–33.
Nagy LG, Riley R, Tritt A, Adam C, Daum C, Floudas D, et al. Comparative genomics of early-diverging mushroom-forming fungi provides insights into the origins of Lignocellulose decay capabilities. Mol Biol Evol. 2015;3:msv337.
Fernández-Fueyo E, Ruiz-Dueñas FJ, Martínez MJ, Romero A, Hammel KE, Medrano FJ, et al. Ligninolytic peroxidase genes in the oyster mushroom genome: heterologous expression, molecular structure, catalytic and stability properties, and lignin-degrading ability. Biotechnol Biofuels [Internet]. 2014;7:2.
Alfaro M, Oguiza JA, Ramírez L, Pisabarro AG. Comparative analysis of secretomes in basidiomycete fungi. J Proteome. 2014;102:28–43.
Vaaje-Kolstad G, Westereng B, Horn SJ, Liu Z, Zhai H, Sorlie M, et al. An oxidative enzyme boosting the enzymatic conversion of recalcitrant polysaccharides. Science. 2010;330:219–22.
Müller G, Várnai A, Johansen KS, Eijsink VGH, Horn SJ. Harnessing the potential of LPMO-containing cellulase cocktails poses new demands on processing conditions. Biotechnol Biofuels. 2015;8:187.
Muszewska A, Hoffman-Sommer M, Grynberg M. LTR retrotransposons in fungi. PLoS One. 2011;6:e29425.
Gao X, Hou Y, Ebina H, Levin HL, Voytas DF. Chromodomains direct integration of retrotransposons to heterochromatin. Genome Res. 2008;18:359–69.
Lippman Z, Martienssen R. The role of RNA interference in heterochromatic silencing. Nature. 2004;431:364–70.
Pereira V. Insertion bias and purifying selection of retrotransposons in the Arabidopsisthaliana genome. Genome Biol. 2004;5:R79.
Labbe J, Murat C, Morin E, Tuskan GA, Le Tacon F, Martin F. Characterization of transposable elements in the ectomycorrhizal fungus Laccaria bicolor. PLoS One. 2012;7:e40197.
The authors want to thank the helpful advices of the URGI team and especially Véronique Jamilloux on the use of REPET package, as well as Francis Martin, Joey Spatafora and In-Geol Choi for allowing the use of unpublished genome data and Bernard Henrissat for the annotation of CAZYmes.
Availability of data and material
Datasets generated and analyzed during the current study are available in the NCBI-SRA repository under accession number SRP086489, and in Mycocosm database http://genome.jgi.doe.gov/Conol1/Conol1.home.html.
This work was supported by Spanish National Research Plan (Project AGL2014–55971-R) and FEDER funds; Public University of Navarre. The work of RC is supported by an FPI grant from the Spanish Ministry of Economy and Competitiveness. The work conducted by the U.S. Department of Energy Joint Genome Institute, a DOE Office of Science User Facility, is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Funders had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Coniophora olivacea strain MUCL 20566 was obtained from the Spanish Type Culture Collection.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Text S1. Supplementary methods. (DOCX 16 kb)
Comparison of TE content estimation from REPET and Repeatexplorer. (XLSX 9 kb)
Snapshot of synteny dot plot between C. olivacea and C. puteana. (TIFF 582 kb)
Comparison of CAZy proteins annotated in C. olivacea and C. puteana. (XLSX 7 kb)
Protein IDs of genes belonging to the 70 gene families with more than four members. (XLSX 9 kb)
Dataset S1. TE annotation in C. olivacea. Classification information at the order level is included in the output format of PASTEC. (GFF3 1342 kb)
Dataset S2. TE annotation in C. puteana. Classification information at the order level is included in the output format of PASTEC. (GFF3 1273 kb)
Dataset S3. TE annotation in H. pinastri. Classification information at the order level is included in the output format of PASTEC. (GFF3 1674 kb)
Dataset S4. TE annotation in S. lacrymans. Classification information at the order level is included in the output format of PASTEC. (GFF3 8198 kb)
Dataset S5. TE annotation in P. tinctorius. Classification information at the order level is included in the output format of PASTEC. (GFF3 12,724 kb)