The genome and transcriptome of the pine saprophyte Ophiostoma piceae, and a comparison with the bark beetle-associated pine pathogen Grosmannia clavigera

Genome assembly

Genome features and annotation

We used the Maker annotation pipeline [31] to predict genes. Within the annotated genome of O. piceae, we identified genes and gene families for secondary metabolite processing, cytochrome P450 as well as ABC transporters. We also identified homologous O. piceae and G. clavigera proteins based on reciprocal best BLAST hits. We further characterized the MAT idiomorph that is responsible for the mating type of the sequenced strain.

Although O. piceae and G. clavigera share hosts, cause sap-stain in pine, and are in sister clades in the Ophiostomatales [32] (Figure 1), their genomes showed no large-scale synteny (Additional file 1). This is consistent with synteny being lost with evolutionary time between members of the class Dothideomycetes [33]. Dothideomycetes is a sister clade of the class Sordariomycetes, which include O. piceae and G. clavigera, and we anticipate that similar synteny losses have occurred within this class. A BLAST comparison of the two predicted proteomes showed that 5,450 proteins were reciprocal best hits. These included most of the major metabolic functions. The O. piceae proteins with no significant homolog in the G. clavigera genome were overrepresented by protein kinases (GO:0004672), sequence-specific DNA binding RNA polymerase II transcription factors (GO: 0000981) and zinc ion binding proteins (GO:0008270) (Additional file 2). In addition, proteins involved in transmembrane transport (GO:0055085) were also significantly overrepresented in this group of 3,469 proteins (Additional file 2). Over 40% (1,397) of the O. piceae proteins with no evident homologs in G. clavigera were proteins of unknown function (predicted or hypothetical proteins). None of the six carboxylic ester hydrolases (GO:0052689) in the O. piceae genome had a homolog in the G. clavigera genome.

We searched for genes that may be involved in producing secondary metabolites (SMs). Such genes are typically organized as contiguous genomic clusters and can be identified by tools like SMURF [34], which uses hidden Markov models that consider genomic context and domain content. The first step in fungal SM biosynthesis is usually catalyzed by ‘backbone’ genes like nonribosomal peptide synthases (NRPSs), polyketide synthases (PKSs), hybrid NRPS-PKS enzymes, prenyltransferases, and terpene cyclases [34]. SMURF, which does not identify clusters containing terpene cyclases, identified thirteen backbone genes of which eleven are in SM clusters in O. piceae (Additional file 3), and nineteen genes in fourteen clusters in G. clavigera[35]. All the O. piceae SM genes have homologs in G. clavigera.

Melanin is a secondary metabolite that is produced by O. piceae and related species, but, as in O. piceae, the genes responsible for its production do not always occur in a cluster. Melanin is synthesized through the 1,8–dihydroxynaphthalene (DHN) pathway [36]. In O. piceae we identified a number of genes that were similar to genes that have major roles in the DHN pathway in Ophiostoma, Grosmannia and Ceratocystis species [12, 37, 38]. These genes included a PKS (OPP_00823), two reductases (OPP_02710, OPP_00820) and a scytalone dehydratase (OPP_07153). PKS catalyze both the elongation of five ketide subunits and the cyclization of these units to form the base ring of naphthalene. The first reductase (OPP_02710) converts 1,3,6,8-hydroxynaphthalene to scytalone, while the second (OPP_00820) transforms scytalone to vermelone.

O. piceae is a heterothallic species, so requires two individuals with different mating types for sexual reproduction and production of fertile fruiting bodies. Our genome annotation identified O. piceae’s MAT1-2 idiomorph (OPP_06680). A truncated MAT1-1 gene was next to the MAT1-2 gene, as in Grosmannia and related species [8, 39, 40]. We have produced perithecia for O. piceae by mating UAMH-11346 with UAMH- 11672 (Figure 2), and we have successfully amplified and sequenced the full length of the MAT1-1 loci in this latter strain; the alpha box (~978 bp) was missing in the sequenced strain.

Gene expression patterns

To identify genes that may be critical for the saprophyte O. piceae to grow in the presence of the nutrients and defence chemicals that are characteristic of its natural pine sapwood substrate, we profiled gene expression for the fungus growing on solid agar media supplemented with simple carbon sources (i.e. sugars and lipids), with pine sawdust, or with pine terpenes (see Methods). Mapping the RNA-seq reads to the predicted gene models identified 7,157 genes that had an abundance of at least 10 FPKM (fragments per kilobase of exon per million fragments mapped) in any of the conditions tested.

To select genes that were highly differentially regulated under different growth conditions, we required a gene to have an FPKM abundance that was at least ten times higher in a specific condition, or a related set of conditions, than in all other growth conditions. This approach identified 677 genes whose transcripts were differentially abundant in at least one growth condition (Additional file 4), and 173 genes whose transcripts were differentially abundant in only one condition (Figure 3). By manually comparing the set of 173 genes to functional information in the Gene Ontology (http://www.geneontology.org/) and to reference metabolic pathways KEGG (http://www.genome.jp), we identified pathways that were likely involved in the response of O. piceae to the growth conditions tested in our experiments. We added support for these pathways by manually identifying genes from the 677-gene set (Additional file 4) whose transcripts, while up-regulated, did not pass the stringent 10-fold filter used to identify the set of 173 genes.

We also assessed alternative transcript splicing across the range of growth conditions used for this study (Additional file 5). Splicing appeared not to be an important factor under these conditions.

Growth on mannose, fatty acid (oleic acid) and triglycerides

Mannose is a simple monomeric epimer of glucose and can be readily utilized as a carbon source by O. piceae. We found five genes whose expression was at least ten times higher with mannose than in any other conditions tested (Additional file 4 shown light blue). These included two transporters, one oxidoreductase and two hypothetical proteins. Our data suggests that mannose uptake involves two transporters (OPP_03031, OPP_05665), and a simple isomerisation/ epimerization reaction by an oxidoreductase (OPP_00733) converts it into glucose. The function of two remaining up-regulated genes (OPP_02416, OPP_07274) is unknown.

We grew O. piceae on triglycerides and fatty acids, which are important lipid compounds in lodgepole pine sapwood, and are a major source of carbon for O. piceae[41]. Because most sources of triglycerides contain a small proportion of fatty acids, it was not surprising that most of the 129 genes whose transcripts were differentially abundant between these conditions were highly up-regulated in both of the conditions. Of the 25 up-regulated genes that were significantly induced only in these two conditions (shown in gold in Additional file 4), 18 were predicted to produce secreted proteins. Unexpectedly, the differentially up-regulated genes included no fungal lipases, which are necessary for the hydrolysis of triglycerides (Additional file 4). Twenty-three of the 25 up-regulated genes were predicted to be involved in the breakdown of carbohydrates and sugars; these included eight genes coding for secreted proteins in the glycoside hydrolase family and four genes for secreted proteins involved in carbohydrate and starch binding. We identified a transcription factor (OPP_02429) that showed significant up-regulation in the presence of triglycerides and oleic acid.

One of the genes differentially expressed between olive oil and oleic acid was a cytochrome P450 (OPP_02426) with a significantly higher expression with triglyceride than with fatty acid. Like its G. clavigera homolog (CMQ_5365; CYP630B18) and homologs in several other species including Fusarium graminearum, Aspergillus niger, A. fumigates and others, this gene is in close proximity to genes encoding a myo-inositol transporter, ARCA-like protein and a cytochrome P450 reductase [35].

Growth on pine sapwood, a natural substrate for O. piceae

O. piceae has slower growth rates than G. clavigera on rich media and on wood; for example, here, on MEA, O. piceae grows 5.8 mm/day, while in Wang et al. [21], G. clavigera grew 14 mm/day. Of the treatments used in this study, sawdust obtained by grinding pine sapwood was the closest to the natural substrate. It contains a variety of carbon sources including mannose, triglycerides and fatty acids. In this growth condition, 366 genes were up-regulated, 91 of which were up-regulated only in the presence of sawdust (Additional file 4). The subset of 91 genes was overrepresented in GO terms for transport (GO:0005215, GO:0006810; p < 0.0001) (Figure 4), which may reflect the complexity of the nutrient sources used by O. piceae. The up-regulated transporters included several allantoate, urea, hexose, iron and sugar transporters, and other major facilitator superfamily (MFS) transporters. As well, oxidoreductase genes that code for putative proteins involved in the modification of aromatic compounds, including phenolics, were highly up-regulated (e.g. P450s, dehydrogenases).

One of the above eight up-regulated genomic clusters contained seven genes that were involved in metabolizing quinic acid (OPP_8732 to OPP_8738). These include a quinate permease, two regulatory genes (an activator and a repressor), and the four genes of the quinate/shikamate catabolic pathway [43, 44]. The latter four catabolic genes (OPP_08735 to OPP_08738) suggest that O. piceae uses quinic acid in wood as a carbon source. While this gene cluster has been reported in many fungi, we were unable to find the cluster in G. clavigera. To confirm that O. piceae can use quinic acid, while G. clavigera cannot, we showed that the former, but not the latter, grows on YNB media with quinic acid as the sole carbon source (Additional file 6). Finally, we identified a secreted lipase (OPP_00605) with predicted triglyceride degradation activity, whose transcript abundance was at least 50-fold higher in the mycelium grown in pine sapwood than in the control mannose.

Tolerance of pine tree defence chemicals

O. piceae did not grow when a mixture of monoterpenes (MT) (i.e. (+)-limonene, (+)-3-carene, racemic a-pinene and (−)-ß-pinene at a ratio of 5:3:1:1) was the only carbon source in YNB. However, after a month of incubation in the presence of MT, the inocula resumed normal growth when they were transferred from YNB + MT to MEA. This suggests that O. piceae is able to survive in the presence of very high levels of monoterpenes. When the fungus was inoculated on MEA and treated with different amounts of MT, the growth rate was only significantly affected when at least 100 μl /plate (~ 0.7 g/l) of MT was added (Figure 5). For all MT treatments, the mycelia were more aerial and fluffy, while the asexual reproduction structures (i.e. formation of synemata) were highly inhibited (Additional file 7).

In order to identify genes involved in terpene tolerance, we grew O. piceae on CM and treated it with a mixture of terpenes as previously described in our studies with G. clavigera[8, 21]. While the experiments for the two species were done at different times, we used the same conditions for both growth experiments and the same protocols for RNA extractions. We compared gene expression profiles of O. piceae after 14 h and 40 h treatments to profiles for untreated CM plates at the same time points. At 14 h, 295 genes were differentially abundant, 261 of which were down-regulated. While carbohydrate metabolism (GO:0005975) was associated with the down-regulated genes (p < 0.001), we were unable to identify any GO terms that were associated with the 34 up-regulated genes. After 40 h in the presence of terpenes, 264 genes were differentially abundant, 126 of which were up-regulated. While carbohydrate metabolism was still associated with down-regulated genes (p < 0.001), several transporters (OPP_06758, OPP_05515, OPP_03974, OPP_02103) were significantly up-regulated. In G. clavigera, which is able to utilize terpenes as a carbon source, more than 250 genes were up-regulated by at least 2-fold at 12 h and 36 h in the presence of terpenes [8]. Of the 34 O. piceae genes that were up-regulated at 14 h, 26 had homologs in G. clavigera, and nine of these G. clavigera genes were up-regulated at 12 h. Similarly, of the 126 O. piceae genes that were up-regulated at 40 h, 75 had G. clavigera homologs, twenty of which were up-regulated at 36 h.

We found 26 O. piceae genes that were up-regulated only in the presence of terpenes, at one or both time points. Eighteen of these had G. clavigera homologs. The most highly up-regulated O. piceae gene (OPP_06758) encoded an ABC transporter that was homologous to the G. clavigera transporter (CMQ_4184; GcABC-G1) that confers terpene tolerance to the pathogen [21] (Figure 6). Approximately 1,500 bp upstream of the O. piceae ABC transporter is a gene (OPP_06759) encoding a transcription factor whose expression, like that of the transporter, was up-regulated only in the presence of terpenes (Figure 6). We found that the O. piceae OPP_06758 and the G. clavigera CMQ_4184 (GcABC-G1) ABC transporters were placed in the same clade when we constructed a phylogenetic tree for a subset of the fungal species recently analyzed by Wang et al. [21]. To our knowledge, from currently available sequence data, this clade is unique to these two fungal species (Figure 7); no similar fungal ABC transporter has been reported.

Strain and growth conditions

Genome sequencing

DNA was extracted from fungal hyphae grown on MEA using methods described by Haridas and Gantt [56]. We used two sequencing technologies: Illumina HiSeq 2000, which generated 100 nt reads, and 454 Titanium, which generated reads with a mean length of 318 nt. The libraries for Illumina HiSeq 2000 had two different insert sizes, 200 and 700 nt, while the library for 454 had an insert size of 8000 nt. Illumina HiSeq sequencing was done at the BC Genome Sciences Centre in Vancouver, Canada and 454 sequencing was done at the Plate-forme d’Analyses Génomiques at Laval University in Québec, Canada.

Genome assembly

Reads generated by the two platforms were used with no further processing for genome assembly using ABySS v1.3.0 [25] with a kmer size of 60. This assembler filters the FASTQ sequences based on quality scores. In order to efficiently use the 454 reads for scaffolding, we used a minimum contig size (1000) and read pairs for building scaffolds (2) (SCAFFOLD_OPTIONS = ‘-s1000 -n2’). The assembly was scrubbed and gaps closed with Anchor (v0.3.0; http://www.bcgsc.ca/platform/bioinfo/software/anchor). When Abyss is unable to find overlaps between contigs where paired end data suggests that the contigs should overlap, it joins the contigs with a single lowercase ‘n’. Such overlaps were resolved using transcriptome assembly (described below) or by finding small overlaps at the ends of the contigs using exonerate v2.2.0 [26]. Genome synteny was assessed by MUMmer [57].

Transcriptome assembly

RNA-seq was performed on eight RNA samples extracted from the mycelia of O. piceae hyphae grown under various conditions as shown in Table 5. For each RNA-seq library, we collected samples from three biological replicates, extracted RNA separately from each replicate, and pooled the samples for paired-end sequencing on an Illumina HiSeq (Canada’s Michael Smith Genome Science Center, Vancouver). Multiplexed sequencing with three libraries per lane was done using the Illumina HiSeq platform to obtain 100 bp paired end reads from 250 bp fragments. Reads were analysed using fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and showed read bias and in the first few bases of the reads and poor quality in the last few. Reads with minimum (Phred) quality scores less than 20 were removed and the first six and last four bases of all reads were trimmed using prinseq [58]. Processed RNA-seq reads were assembled using Trinity [28] using the jaccard_clip option to minimize fusion transcripts. The best protein coding transcripts were identified using the included scripts and aligned back to the assembled genome using exonerate v2.2.0 est2genome [26].

Genome annotation

We used the Maker annotation pipeline (v2.26) for genome annotation [31]. In addition to the trinity assembled best candidates, we also used two additional sources of evidence in Maker. The first was transcripts predicted by the Core Eukaryotic Genes Mapping Approach [27], (CEGMA) and the second was coding sequences of transcripts assembled by cufflinks [59] from RNA-seq reads mapped to the assembled genome. Within the Maker framework, we trained SNAP v2006-07-28 [60] using the Trinity assembled transcripts, gene models of Magnaporthe grisea for Augustus (v2.5.5) [61] and an hmm file for Genemark-ES (v2.3) [62] using an independent run. The UniProtKB/Swiss-Prot (release 2012_01) fasta file was provided as protein homology evidence and pred_flank was set to 50 to minimize fusion transcripts. Predicted genes smaller than 100 amino acids were removed unless they were at least 80 amino acids long and had transcript, protein or CEGMA evidence. Selected gene models were manually curated. Functional identification of predicted genes was done using Blast2go (v2.5.1) [63]. tRNA’s were identified using tRNAscan-SE (v1.3.1) [64]. Relative synonymous codon usage (RSCU) was calculated using a local installation of the graphical codon usage analyser [65]. Secretome predictions were made with TargetP [66] and Phobius [67]. A protein was considered to be secreted if either TargetP or Phobius suggested that it was secreted and this result was not in conflict with the other. Identification of secondary metabolism genes and clusters was done using the Secondary Metabolite Unique Regions Finder (SMURF) [34].

RNA-Seq analysis

Quality trimmed RNA-seq reads were aligned to the O. piceae genome using Bowtie (v0.12.7), Tophat (v2.0.4), Cufflinks (v2.0.2) as described by Trapnell et al. [59]. Because mapping the RNA-seq reads to the genome without providing fixed gene models resulted in an unacceptable number of predicted fusion transcripts, reads were mapped using the curated gene models predicted by the Maker pipeline.

Data availability

The sequences reported in this paper are being deposited in NCBI gene bank as assembly and annotations, Project NO. PRJNA182071.