De novo genome assembly of the soil-borne fungus and tomato pathogen Pyrenochaeta lycopersici

Genome sequencing and assembly

Genomic DNA obtained from the virulent P. lycopersici isolate CRA-PAV_ER 1211 was sequenced using Illumina 100-bp paired-end reads. Two libraries were prepared with mean fragment sizes of 460 and 560 bp respectively and sequenced obtaining 85 millions of fragments for a total of 17Gb, corresponding to approximately 300-fold coverage of the final assembly.

The P. lycopersici genome was assembled de novo using a de Bruijn graph-based (DBG) assembly strategy with a k-mer of 55 (Additional file 1: Figure S1) before reassembly using the overlap-layout-consensus (OLC) method. We chose to use this strategy to take advantage of both the higher efficiency of DBG method in assembling millions of reads and the higher sensitivity of OLC method [22, 23]. In fact, the OLC re-assembly of sequences previously assembled with DBG method reduced the total number of sequences by a 40% (from 11,617 to 7,079) while increasing the average contig length by a 60% (from 4,834 to 7,747 bp) and not reducing significantly the total number of assembled bases (54.8 Mbps vs. 56.2 Mbps).

Gene annotation

Genes were annotated by a combination of ab initio prediction [25] followed by the reference annotation-based transcript (RABT) assembly of sequences obtained from the RNA-Seq analysis of two P. lycopersici samples grown in vitro and during interaction with tomato cv Moneymaker, respectively [26]. We annotated 17,411 genes with 27,275 transcripts, among which 9,553 loci were detected by both approaches and 2,797 only by the analysis of RNA-Seq data (Table 1). Thus, while the number of genes might be slightly overestimated because of ab initio prediction limitations, at least 70.8% of the annotations were supported by experimental evidences. We confirmed that the assembly was able to capture full-length genes by searching the predictions for full open reading frames (ORFs), finding that most of the genes (80.76%) contained start and stop codons. The RNA-Seq data were then assembled de novo and the resulting contigs were mapped onto the draft genome. From the 31,746,550 reads, we obtained 27,982 putative transcripts, 27,574 (98.5%) of which could be mapped onto the assembled genome further validating the comprehensiveness of the gene space represented.

Many of the identified transcripts (75.8%) were conserved in other species as shown by hits against sequences in the NCBI NR protein database (e-value < 1E-0.6) and the Uniprot SwissProt Fungi protein database (e-value < 1E-0.7). This analysis showed that the most represented species were fungal pathogens and the top four species were phytopathogenic fungi (Additional file 1: Figure S2). Based on the identity of the most similar proteins, at least one Gene Ontology term was assigned to 8,507 transcripts. The most represented functional categories are summarized in Additional file 1: Figure S3. We found 3,962 orphan genes (22.8%) with no matches against known proteins or protein domains [27], which is similar to the proportion of orphan genes found in Sordaria macrospora (22%) [28] and Macrophomina phaseolina (29%) [29]. RNA-Seq analysis showed that 1,936 of the 3,962 orphan genes were transcribed and were therefore likely to be functional.

Phylogenetic relationships

The phylogenetic analysis based on whole genome sequences of P. lycopersici and 16 other fungal species with diverse lifestyles (from necrotrophs to biotrophs) is shown in Figure 2. This confirmed that P. lycopersici belongs to the class Dothideomycetes and the order Pleosporales. A neighbor-joining phylogenetic tree, using Rhizopus oryzae as the outgroup, revealed four major clusters representing Saccharomycetes, Leotiomycetes, Sordariomycetes and Dothideomycetes. P. lycopersici appeared to be more closely related to hemibiotrophic and necrotrophic plant pathogens of the genera Leptosphaeria, Pyrenophora and Stagonospora than to biotrophs such as the genera Blumeria and Ustilago.

We performed a pairwise comparison between P. lycopersici and Fusarium oxysporum f. sp. lycopersici genomes to assess the genomic distribution of similarity regions. We choose F. oxysporum among the top ranked species (Additional file 2: Table S2) as a complete genome anchored to chromosomes was available for this fungus. As expected the analysis showed regions of similarity at aminoacidic level throughout all the F. oxysporum chromosomes with the exception of the four F. oxysporum dispensable chromosomes 3, 6, 14 and 15 (Additional file 1: Figure S4).

Vegetative incompatibility

Fungi can propagate both by sexual and vegetative reproduction. The latter is controlled by approximately 50 heterokaryon incompatibility (HET) modules in most pathogenic fungi, whereas there was evidence of 284 modules in P. lycopersici (Figure 3a; Additional file 2: Table S4). Accounting for a possible overestimation of the total number of genes by a 40%, based on the percentage of annotations not supported by RNASeq data, the HET protein modules are anyway expanded by 3.9 times, in average, compared to other fungi. Other proteins that are functionally associated with HET contain NTPase, NB-ARC and NACHT domains, which are involved in the regulation of the immune response and are related to apoptosis/programmed cell death in animals and fungi [30]. We identified 53 NACHT proteins encoded by the P. lycopersici genome, which is slightly more than the number in F. oxysporum (41) and much higher than in other fungi, which range from 0 in Blumeria graminis to 10 in Colletotrichum graminicola. The P. lycopersici genome also encoded a significantly greater number of proteins containing ankyrin (ANK) and tetratricopeptide repeats (TPR), which mediate protein-protein interactions among HET proteins. Altogether, the number of HI related proteins includes 522 domains, e.g. more than double the number encoded by the imperfect fungal pathogen F. oxysporum. About 75% of the modules for all the families related to vegetative incompatibility were supported by RNASeq expression data.

Pathogenesis related genes

The ATP-binding cassette (ABC) transporters and MFS-type membrane transporters were largely represented in the P. lycopersici genome, comprising 109 and 229 modules respectively (Figure 3b and Additional file 2: Table S2). Transcriptome analysis showed that 85% of ABC modules and 75% of MFS-type membrane transporters belonged to transcripts expressed. The P. lycopersici genome also encoded a large number (125) of cytochrome P450 proteins, as shown for other necrotrophic and hemibiotrophic pathogens such as B. cinerea (120), C. graminicola (138), F. oxysporum (157) and S. nodorum (122).

The P. lycopersici predicted genes also included 597 sequences matching 94 different subfamilies of peptidases (Additional file 2: Table S3). Almost all known peptidase families were represented in the P. lycopersici transcriptome, with DmpA aminopeptidase 1 (P1) as the only exception. The most represented clans were the serine peptidase (S clan) and metallopeptidase (M clan), a common feature of fungal pathogens. The comparative analysis of gene families and PFAM domains with several other fungi showed that the number of peptidases in P. lycopersici is comparable to other hemibiotrophic pathogens as C. graminicola and F. oxysporum. Metallopeptidase group of aminopeptidases (M1, M18, M24, M28) were highly represented, with 43 domains in P. lycopersici (shown in yellow in Additional file 2: Table S3) compared to 34 and 30 of C. graminicola and F. oxysporum respectively. Analysis of expression data confirmed that 36 of this modules (85%) corresponded to expressed transcripts. Other metallopeptidase gene families had also undergone expansion in the P. lycopersici genome, including pappalysin, Ste24 and deubiquitinating peptidase (shown in green in Additional file 2: Table S3). A summary is reported in Figure 3d.

Genes involved in carbohydrate degradation (CAZymes)

Carbohydrate degradation is an important component of fungal pathogenicity and virulence, so we examined the P. lycopersici CAZome in detail compared with nine other fungi with complete genome sequences (Additional file 2: Table S6). The P. lycopersici genome encodes 575 CAZyme modules, including 272 glycoside hydrolases, 83 glycosyltransferases, 20 polysaccharide lyases, 149 carbohydrate esterases and 51 carbohydrate-binding modules (CMBs). CE1 and CE10 families of carbohydrate esterases including acetyl xylan esterase (EC 3.1.1.72), cinnamoyl esterase (EC 3.1.1), feruloyl esterase (EC 3.1.1.73), carboxylesterase (EC 3.1.1.1), S-formylglutathione hydrolase (EC 3.1.2.12) and sterol esterases are particularly represented (53 CE1 modules out of which 77% are detected as expressed by RNASeq analysis, 54 modules for CE10 family 73,5% of which are expressed) similarly to other hemibiotrophs such as C. graminicola and F. oxysporum (Additional file 2: Table S6). Cellulose-degrading enzymes were also well represented in P. lycopersici, including seven cellobiohydrolases (GH6), ten β-1,3-glucanases (GH55) and eight GH105 glycoside hydrolases (unknown mechanism). Finally, the complex of carbohydrate cleaving enzymes was completed by 20 genes encoding polysaccharide lyases, including 13 encoding PL3.

Sample and library preparation

Genomic DNA was isolated from a virulent P. lycopersici isolate (CRA-PAV_ER 1211) grown on Potato Dextrose Agar (PDA) medium, by the method described by Cenis as previously described [62] with the following modifications: 200 mg of mycelium was frozen in liquid nitrogen, pulverized, and incubated in 300 μl of lysis buffer (200 mM Tris-HCl pH 8.5, 250 mM NaCl, 25 mM EDTA, 0.5% SDS) for 10 min at 65°C. We then added 150 μl 3 M sodium acetate (pH 5. 2), incubated at –20°C for 10 min and centrifuged at 10000 × g for 30 min. The supernatant was transferred in a fresh tube and the DNA was precipitated by adding an equal volume of isopropanol, centrifuging at 10000 × g for 10 min and washing with 70% ethanol. Finally the DNA was purified using the NucleoSpin Extract II kit (Macherey-Nagel, Düren, Germany). Total RNA from vegetative mycelium grown on PDA medium was extracted using the RNeasy Midi kit (Qiagen, Hilden, Germany). Total RNA from infected tomato roots of cv Moneymaker was extracted using the NucleoSpin RNA Plant 2 kit (Macherey-Nagel) after 8 days post infection (dpi).

Genomic DNA (6 μg) was fragmented by nebulization at 35 psi for 6 min. DNA libraries with insert sizes of 400 bp and 560 bp were prepared from 1 μg of fragmented genomic DNA using the paired-end TruSeq DNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA). The library quality was determined using the High Sensitivity DNA Kit (Agilent, Wokingham, UK).

Total RNA samples were assessed for quality using an RNA 6000 Nano Kit (Agilent) and 2.5-μg aliquots were used to isolate poly(A) mRNA for the preparation of a non-directional Illumina RNA-Seq library using the TruSeq RNA Sample Prep Kit (Illumina). The quality of the library was checked using the High Sensitivity DNA Kit (Agilent).

Sequencing and data preprocessing

Libraries were sequenced with an Illumina GAIIx sequencer generating 100-bp paired-end sequences for DNA libraries and 130-bp paired-end sequences for RNA libraries.

The sequences were pre-processed by removing reads with a number of N >10 or with a read quality <20 using a custom script. Adapters were clipped using Scythe v0.980 [63] and bases on both 3′ ends with a quality <20 were trimmed using Sickle v0.940 [64], eventually entirely removing the fragment if the length was reduced to < 50 bp.

De novo assembly and gene catalog assessment

The genome was assembled de novo using Velvet v1.1.06 [65] with the following parameters: –exp_cov auto (automatic calculation of expected coverage), –scaffolding (scaffolding of contigs with paired-end reads) and –min_contig_lgth 200 (mimimun contig length = 200). The optimal k-mer length was determined by adjusting the k-mer length from 39 to 67 bp in 4-bp increments and using the k-mer for which the N50 and the maximum contig length reached the highest value (Additional file 1: Figure S1). The resulting contigs were then re-assembled with CAP3 v10/15/07 [66] using standard parameters.

Core Eukaryotic Genes (CEGs) were aligned with assembled genome using BLAST and hits were considered significant when the sequence identity was >65%.

RNA-seq reads from in vitro mycelia were assembled using Trinity v r2011-11-26 [67] with standard parameters, jaccard clip on and a minimum contig length of 200 bp. Assembled contigs were mapped using GMAP [68] with standard parameters.

Gene annotation

The final assembly was processed by GeneMark.hmm-ES v2.3e [25] with standard parameters and no a priori information. The genome was masked for repetitive and low-complexity regions with RepeatMasker v open-3.3.0 [69] with standard parameters and a general repeats database. The resulting annotation was refined using TopHat v1.4.1 and Cufflinks v1.2.1 [70] on the two RNA-seq libraries in RABT mode with standard parameters. ORFs were identified for each transcript using CPC v0.9.r2 [71].

Phylogenetic analysis

A phylogenetic tree based on the comparison of whole genomes of P. lycopersici and 16 other fungi was constructed using CVTree v2 [72] at a k-mer length of 7. Rhizopus oryzae was used as the outgroup for building the unscaled tree [73]. The P. tritici-repentis and C. graminicola proteomes were obtained from the Colletotrichum Sequencing Project [74], the L. maculans proteome was obtained from the L. maculans genome project [75], and the B. graminis proteome was obtained from the Blumeria sequencing project [76]. Unless otherwise stated, the remaining proteomes were obtained from the CVtree inbuilt genome database [77].

Functional annotation

Functional annotation was initiated by using each sequence as a BLAST query [78] against the NCBI Non Redundant database retrieved 2012-09-14 [79] (e-value < 1E-0.6) and the Uniprot SwissProt Fungi protein database retrieved 2012-04-30 [80] (e-value < 1E-0.7). The results were analyzed using Blast2GO [81] and integrated with InterPro results [82]. Conserved protein domains in P. lycopersici and all the other fungi considered (AN, BC, BG, CG, FO, LM, NC, PT-R, SN) were identified using HMMer v3.0 [83] to identify homology with proteins in the Pfam-A database (v26.0, 2011-11) [84]. Sequence conservation was considered significant at an e-value threshold <1e-6 for both the entire sequence match and for the independent E-value of the single domain match. CAZymes (v2.0) [85] homology was also inspected using HMMer. Alignments were considered significant at an alignment length > 80 residues, E-value < 1e-5 and HMM profile coverage > 30% or alignment length < 80 residues and E-value < 1e-3 and HMM profile coverage > 30%. BLASTX was used to identify sequences homologous to known pathogenic genes (PHI-base ver 3.2) [86], peptidases (MEROPS ver 9.8) [87], zinc fingers (C2H2 ZNF db, ver. 2007-10-03) [88], MAP kinase sequences from NCBI NR (retrieved 2013-01-30) [79] and membrane transport proteins (TCDB, ver. 2011-July-15) [89] (e-value <1e-10). Sequences with significant hits against membrane transport proteins were also used as BLASTX queries against a G protein-coupled receptors database (GPCRDB, retrieved 2013-01-30) [90] and fungi major facilitator superfamily sequences from NCBI NR [79].

Significance of protein families abundance differences between P. lycopersici and all the other fungal plant pathogens (BC, BG, CG, FO, LM, PT-R, SN) was assessed by a 1-sample t-test. Variance was estimated based on protein families abundance data of all the fungal pathogens taken into account (P. lycopersici excluded). P-values were corrected according to Benjamini and Hochberg [91] on the full dataset of comparisons.

Genome coverage has been estimated by mapping the reads on the assembled genome using BWA v. 0.6.2-r126 [92] using default parameters and calculating coverage on a panel of 140 single copy genes.

Data access

RNASeq reads and transcriptome assemblies have been deposited at the NCBI Sequence Reads Archive (SRA) and NCBI Transcriptome shotgun assembly (TSA) databases respectively and are available under BioProject number PRJNA202292. Genomic reads and genome assembly have been deposited at the NCBI Sequence Reads Archive (SRA) and are available under BioProject number PRJNA202288. Assemblies have been deposited as Whole Genome Shotgun project at DDBJ/EMBL/GenBank under the accession ASRS00000000. The version described in this paper is version ASRS01000000.