Monthly Archives: September 2013

Variant analysis of different isolates and fruiting bodies of Chalara fraxinea

Contributors

Christine Sambles, David Studholme. University of Exeter, Devon.

Introduction

To identify the extent of genetic variation in the sequenced samples, we undertook variant analysis from seven isolates (KW1, FERA 105, FERA 232, FERA 233, FERA 88, FERA 93 & FERA 94), two fruiting body samples (MFB1 & PFB1) and four mixed material samples (HP1, UB1, AT1 and AT2).

Materials

Raw reads:

C. fraxinea:  KW1

Fruiting bodies: MFB1, PFB1

Mixed material: AT1, AT2, Upton, Holt

 

Genome assembly

C. fraxinea scaffolds:  KW1

Methods

Qualities of raw reads from the thirteen samples were assessed with FASTQC.   Adapter- and quality- trimming was performed with Trim galore, a wrapper script using FASTQC and cutadapt (phred cutoff:20, error rate:0.1, adapter overlap: 1bp, min. length: 20bp, paired read length cut-off: 35bp). FASTQC was automatically run on all trimmed files to confirm trimmed read quality. Raw and trimmed metrics and GC content of raw reads can be seen in Table 1.

Trimmed reads were aligned to the pre-assembled KW1 genome from OADB using the splice-aware aligner Tophat (Trapnell, et al., 2009). The resulting assembly BAM file was used to create a pileup file using MPILEUP from Samtools.

The variant detection software, VarScan2 pileup2snp (or mpileup2snp) was used ( –p-value 0.05 –min-coverage 10 –output-vcf -min-var-freq 0.95) to call SNPs (Koboldt, et al., 2012). The numbers of SNPs were normalised to the number of bases with >=10X coverage to take into account the different depths of sequencing.

Results

The Fera samples contain fewer sequencing reads compared to the fruiting body samples, mixed material samples and KW1 (Table 1). Sample KW1 covers 32% of the genome positions with a coverage of >=10X compared to between 0.17-0.27% for the Fera samples. The remaining samples range from 13% (AT1) to 32% (Upton) (Table 1). This suggests that limited SNP calling is possible in the FERA samples. The KW1 genome assembly has a size of 65Mb and transcriptome assembly has a size of ~35Mb suggesting that spliced transcripts represent 55% of the genome size.

Sample

Sample type

Raw reads x2

GC%

Filtered x2

% 10X KW1 genome positions covered

KW1

C. fraxinea

41,036,951

49

38,594,279

32%

FERA 105

C. fraxinea

1,442,015

48

1,440,015

0.26%

FERA 88

C. fraxinea

1,731,461

47

1,729,869

0.20%

FERA 93

C. fraxinea

1,355,421

47

1,353,849

0.24%

FERA 94

C. fraxinea

1,642,453

47

1,639,096

0.17%

FERA 232

C. fraxinea

1,919,713

48

1,918,416

0.27%

FERA 233

C. fraxinea

1,838,409

47

1,837,097

0.26%

Mature fruiting body (MFB1)

C. fraxinea

15,033,876

49

14,524,487

27%

Primordial fruiting body (PFB1)

C. fraxinea

7,298,316

49

7,187,891

16%

AT1

Mixed material

31,059,879

48

30,388,033

29%

AT2

Mixed material

23,096,021

44

21,924,824

13%

Holt

Mixed material

36,960,651

49

36.231.834

21%

Upton

Mixed material

43,099,662

48

41,367,071

32%

Table 1:  Read number pre- and post-filtering, G+C content of raw data and percentage of KW1 genome positions covered at a depth of >=10X.

SNP calling using VarScan revealed a high level of genetic variation between the samples (Table 2).  The mature fruiting body (MFB1) and primordial fruiting body (PFB1) samples showed a high level of variation (~16,000 SNPs) when aligned to the KW1 genome. Samples AT1 (~6,000), AT2 (~3,000), Holt (~14,000) and Upton (~12,000) also suggested a high level of genetic variation although this could be due to the alignment of non-Chalara reads present in the mixed sample, which have a high degree of similarity to regions of the C. fraxinea genome. This highlights the importance of the availability of multiple genomes due to such high levels of genetic variation.

Sample Sample type # Unambiguous bases # SNPs VarScan Total

KW11

C. fraxinea

19,974,663

59

NA

FERA 88

C. fraxinea

167,378

40

204

FERA 93

C. fraxinea

127,827

20

FERA 94

C. fraxinea

150,981

47

FERA 105

C. fraxinea

104,667

22

FERA 232

C. fraxinea

169,295

53

FERA 233

C. fraxinea

165,989

35

MFB1

C. fraxinea

16,519,209

11,852

16,396

PFB1

C. fraxinea

9,869,668

9,545

AT12

Mixed material

18,108,777

6,029

NA

AT22

Mixed material

8,016,701

3,138

Holt

Mixed material

12,975,549

13,654

Upton

Mixed material

19,789,944

12,259

Table 2: VarScan SNP calling of C. fraxinea isolates, fruiting body and mixed material samples.  1MAT1-1 mating type; 2MAT1-2 mating type; UAP: unambiguous positions.

We identified respectively 40, 20, 47, 22, 53, 35 single-nucleotide differences (10X coverage, 95% consensus) for samples FERA88, FERA93, FERA94, FERA105, FERA232, FERA233 when compared to the KW1 genome. In total, 204 unique SNPs were identified in the FERA samples when multi-sample calling was used. P-values were calculated using a Fisher’s Exact Test on the read counts supporting reference and variant alleles. Of the 204 SNPs, 165 were located within an annotated gene, 155 of these within annotated exons. These 155 SNPs were distributed across 70 genes which are we are currently annotating (BLAST & PFAM).

One of these genes (CHAFR746836.1.1_0031990) has a sequence that is similar (62% identity) to cerato-platanin from the basidiomycete Trametes versicolor (GenBank: EIW62259.1) and was first identified in the ascomycete, Ceratocystis fimbriata, the causal agent of “canker stain disease” in Platanus x acerifolia in Europe (Pazzagli, et al., 1999; Pazzagli, et al., 2006). It belongs to a family of cerato-platanin phytotoxic proteins which are found in the cell wall of the fungus (Boddi, et al., 2004) and are involved in the host-pathogen interaction (Pazzagli, et al., 1999). Six SNPs were identified in this protein. All the FERA samples except FERA88 and FERA105 differ from KW1 at all six SNPs; FERA 88 and FERA105 are identical to KW1 at all six SNPs (Fig 1). At these six SNPs, the primordial fruiting body sample resembles FERA 105 and FERA 88 and differs from FERA 90, FERA 232, FERA233, FERA 94 and the mature fruiting body sample. These two fruiting body samples originate from different sources. In the mature fruiting body sample, which is likely to contain both dikaryotic hyphae and haploid ascospores, all 6 SNPs appear to be heterokaryotic. In all other samples the SNPs are homokaryotic/homozygous.

preview_1215165

Fig 1: (CLick for larger view) IGV view of putative cerato-platanin gene in C. fraxinea. The six SNPs called by VarScan are shown in the top track. The next four tracks show the homokaryotic/homozygous SNPs which are the same in FERA93, FERA232, FERA233 and FERA94, the sixth track shows possible heterokaryotic SNPs in the mature fruiting body (MFB1) and the last three tracks show that the homokaryotic/homozygous SNPs which resemble the KW1 genome sequence are also present in FERA88, FERA105 and the primordial fruiting body (PFB1).

 

Normalisation of the identified SNPs shows that the highest variation occurs in fruiting body samples, Holt and Upton when compared to the KW1 genome (Fig 1). Further analysis is required to understand the significance of this difference in variation and whether the Upton and Holt variation is due to the presence of other fungi, such as Phytothphora sp. and Togninia sp. in the mixed samples interfering with the SNP calling.

 

preview_1215163

Fig 2: Normalised number of SNPs against bases with >=10X . Click for larger view

References

Boddi, S., et al. (2004) Cerato-platanin protein is located in the cell walls of ascospores, conidia and hyphae of Ceratocystis fimbriata f. sp. platani, FEMS Microbiology Letters, 233, 341-346.

Koboldt, D.C., et al. (2012) VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing., Genome research, 22, 568-576.

Pazzagli, L., et al. (1999) Purification, characterization, and amino acid sequence of cerato-platanin, a new phytotoxic protein from Ceratocystis fimbriata f. sp. platani., The Journal of biological chemistry, 274, 24959-24964.

Pazzagli, L., et al. (2006) Cerato-platanin, the first member of a new fungal protein family: cloning, expression, and characterization., Cell biochemistry and biophysics, 44, 512-521.

Trapnell, C., Pachter, L. and Salzberg, S.L. (2009) TopHat: discovering splice junctions with RNA-Seq., Bioinformatics (Oxford, England), 25, 1105-1111.

 

Metagenomic analysis of Fraxinus excelsior, Chalara fraxinea and infected material.

Contributors

Christine Sambles, David Studholme. University of Exeter, Devon.

Introduction

To check their species composition, we performed a metagenomic analysis on several datasets that had been generated from samples described as Fraxinus excelsior, Chalara fraxinea or as mixed infected material. We used MEGAN (v4, MEtaGenome Analyzer (Huson, et al., 2011)) and the assembled transcripts from F. excelsior and C. fraxinea to identify the taxonomic groups to which uninfected sample transcripts are allocated, to use as a reference database for binning. We identified many transcripts specific to the infected samples (not in uninfected Fraxinus) that fell outside of the Chalara and Fraxinus bins; these may indicate additional microbial species present during infection. These might include species acting synergistically with C. fraxinea during infection or opportunists present as a secondary consequence of infection. The uninfected F. excelsior sample identifies species that are part of the ‘normal’ or ‘healthy’ tree microbiota and could therefore be excluded from the list of infection-related species.

Material

Transcriptome assemblies:

F. excelsior: ATU1

C. fraxinea:  KW1

Mixed material: AT1, AT2, Upton, Holt

BLASTX against GenBank:

F. excelsior: ATU1

C. fraxinea: KW1

Mixed material: AT1, AT2, Upton, Holt

Methods

We identified sequence similarity between assembled transcripts and GenBank protein sequences using BLASTX; we used as queries the transcripts from uninfected F. excelsior (ATU1), a C. fraxinea isolate (KW1) and four mixed material samples (AT1, AT2, Upton and Holt). We loaded the output from BLASTX into MEGAN and performed taxonomic binning using a minimum support value of 35, a minimum BLAST score of 50 and only retaining hits whose bit scores lie within 10% of the best score. The analyses were normalised, compared and rendered within MEGAN.

Results

For each of the six samples, we identified the numbers of transcripts assigned to Helotiales (the order in which C. fraxinea is classified) and to Viridiplantae (green plants). The results are summarised in Table 1. Transcripts from the (nominally) C. fraxinea isolate KW1, binned into the classes of Dothideomycetes, Eurotiomycetes, Leotiomycetes and Sordariomycetes, which all reside within the subphylum of Pezizomycotina. This result is consistent with the sequenced sample being pure Chalara. As expected, all of the bin-able transcripts from the F. excelsior ATU1 transcripts fell within the Viridiplantae kingdom, specifically within the group of flowering plants (Magnoliophyta).

Sample

Sample Type

Helotiales

Viridiplantae

Non H/V Percent

Normalised reads

Percent

Normalised reads

Percent

ATU1

F. excelsior

0

0%

79,438

79%

21%

KW1

C. fraxinea

32,350

32%

0

0%

68%

AT1

Mixed material

15,853

16%

47,202

47%

37%

AT2

Mixed material

8,889

8.9%

54,434

54%

37%

Upton

Mixed material

12,619

13%

12,798

13%

74%

Holt

Mixed material

6,588

6.6%

31,741

32%

61%

Table 1: Reads binned to Helotiales and Viridiplantae in normalised comparison for each sample.

In the data from Upton mixed material, 74% of the transcripts were not binned within the Helotiales or Viridiplantae, which is where C. fraxinea and F. excelsior transcripts are expected to fall, based on the results from pure isolate of C. fraxinea and the uninfected F. excelsior. In the Upton data, 34% of the total number of transcripts was assigned to Oomycetes; specifically 33% to Phytophthora spp.. Additionally, 13% are not assigned to any taxon and a further 11% had no significant similarity to proteins in the GenBank database, detectable by BLAST. The presence of Phytophthora spp. might be attributed to cross-lane contamination during sequencing, since the Norwich laboratory handling the Upton data also work with Phytophthora infestans. A similar contamination had also been reported for Fera samples with Maize Chlorotic Mottle Virus (MCMV) and Sugarcane Mosaic Virus (SMV) sequences being present. This is a common problem in Illumina sequencing which has led to the incorporation of a taxonomic binning step into sequencing pipelines including at our own sequencing facility at the University of Exeter. These contamination issues highlight the importance of confirming the taxonomic distribution of sequence data in addition to quality checks before performing any downstream analyses. Once identified, the contaminant reads can be removed from the dataset.

The Holt mixed material sample analysis showed 1.5% transcripts binned to Togninia minima, an ascomycete in the order Calosphaeriales. T. minima is a pathogen of grapevines and Prunus spp., however, the closely related T. fraxinopennsylvanica (anamorph: Phaeoacremonium mortoniae) has been observed in dead vascular tissue of declining ash tree branches (Fraxinus latifolia) in California (Eskalen, et al., 2005a; Eskalen, et al., 2005b). It may be that T. fraxinopennsylvanica is present in the Holt material and that the reads were assigned to the species T. minima because that is the most closely related species for which extensive sequence data is available.

For the nominally pure sample of C. fraxinea isolate KW1, only 32% of reads were assigned to Helotiales. However, 66% of reads were assigned to Fungi with 16% not assigned to a taxa and 17% had no significant similarity to proteins in the GenBank BLAST database. This is likely to be due to insufficient sequence data in the GenBank database from Chalara and closely related species.

 

preview_1215170

Fig 1: Click for larger view. Metagenomic analysis of Fraxinus excelsior, Chalara fraxinea and four infected material samples (AT1, AT2, Upton and Holt).

Further analysis using alignments to the F. excelsior and C. fraxinea genome will help interpret whether or not other taxa such as Togninia sp., indicated to be present by the MEGAN analysis, are present or whether they are mis-assigned reads due to the lack of Chalara- and Fraxinus-related proteins in the database.

References

Eskalen, A., Rooney-Latham, S. and Gubler, W.D. (2005a) First Report of Perithecia of Phaeoacremonium viticola on Grapevine ( Vitis vinifera ) and Ash Tree ( Fraxinus latifolia ) in California, Plant Disease, 89, 686-686.

Eskalen, A., Rooney-Latham, S. and Gubler, W.D. (2005b) Occurrence of Togninia fraxinopennsylvanica on Esca-Diseased Grapevines ( Vitis vinifera ) and Declining Ash Trees ( Fraxinus latifolia ) in California, Plant Disease, 89, 528-528.

Huson, D.H., et al. (2011) Integrative analysis of environmental sequences using MEGAN4., Genome research, 21, 1552-1560.

 

Analysis of the Chalara genome suggests that it is AT repeat rich

The Contributor

Dan MacLean

The Analysis

I analysed the TGAC KW1 assembly of Hymenoscyphus pseudoalbidus / Chalara fraxinea  to display various features of the assembly. Click for larger image

 

preview_1194882

  1. (outer black ring) – scaffolds of length >= 10 kbp
  2.  (blue stacks) – gene models
  3.  (blue heatmap) – gene density in 10 kbp windows
  4.  (line plot) – aligned distance between paired end read for library with 196 bp insert size (< 31 bp  – 5th percentile – rendered in orange, > 351 bp – 95th percentile – rendered in blue)
  5.  (line plot) – aligned distance between paired end read for library with 570 bp insert size (< 265 bp – 5th percentile – rendered in orange, > 2313 bp – 95th percentile – rendered in blue) Thresholds calculated as per (https://github.com/danmaclean/h_pseu_analysis/blob/master/circos/scripts/assess_insert_size_distributions.md)
  6.  (pink line plot) – uniquely mapped read covereage. Maximum plotted = 300, minimum plotted = 80
  7.  (pink line plot) – GC percent. Black line = 50% GC, light grey area = 50% – 30 % GC, dark grey area = > 30% GC.
  8. Maximum plotted = 60% GC, minimum plotted = 20% GC
  9.  (green stacks) – Repeat Masker matches.

 

http://figshare.com/articles/Genome_analysis_of_the_Ash_Dieback_Fungus_suggests_a_repeat_rich_genome/791640

The interpretation

The genome assembly contains numerous low GC, high coverage, repeat rich regions, with low overall gene density. This suggests an assembly with collapsed repeats of a genome that is overall rich in repeats.

Cite this:

http://dx.doi.org/10.6084/m9.figshare.791640

Lignin degrading enzymes in the Chalara genome assembly

The contributors

Dan MacLean

The analysis

To extend the analysis of Chalara’s ability to degrade wood and thus invade a tree directly I examined the lignin-decaying gene complement in the genome. The previous work on identifying decay related enzymes in [Floudas]: http://dx.doi.org/10.1126/science.1221748  was used as a base.

The list of proteins in the TGAC 1.1 KW 1 genome assembly of Chalara fraxinea and the list of wood decaying proteins in Floudas  were used as input to BLAST searches to identify proteins with strong sequence identity in the Chalara (see [protocol]: https://github.com/danmaclean/h_pseu_analysis/blob/master/interesting_wood_gene_analysis.md ).

For each group in the Floudas decay-related list Chalara proteins were assayed to identify proteins in the with >50% sequence identity to the representative protein in at least 10 members. Counts of the number of Chalara proteins passing this threshold were used as the estimate of the number of members of each group in the Chalara genome.

 

Click for larger imagepreview_1194872

The interpretation

The analysis shows that the Chalara genome is in fact very poor in the sorts of enzymes that are typically used by wood rotting fungi to decay the wood.

 

Assessing the origin of the UK Ash dieback pathogen

The contributors

Kentaro Yoshida, Dan MacLean, Daniel Bunting, Diane Saunders at The Sainsbury Laboratory

The analysis

To clarify the origin of Ash-dieback pathogen, we constructed a phylogenetic tree of UK, French, and Japanese isolates. Based on high quality SNPs and alignments of short reads from RNA sequencing of infected Ash tree samples from seven locations across Norfolk and Suffolk of UK and Japanese samples, we reconstructed the consensus sequences of genes for these tested samples. We searched single copy genes in Fusarium graminearum, Sclerotinia sclerotiorum, Botrytis cinera, Hymenoscyphus albidus, and H. pseudoalbidus using OrhoMCL (Li et al. 2003 Genome Res. 13:2178-89.) The 2,964 single copy genes were identified and used for the construction of the tree. We built a maximum likelihood phylogenetic tree based on third codon positions of the genes using RAxML software (Stamatakis et al. 2005 Bioinformatics 21:456-463).

3rd_codon_trees_large

MLtreeB

Maximum likelihood trees based on 3rd codon positions of 2964 genes.
The number on braches shows bootstrap probability.
(A) All tested samples.
(B) The enlarged tree in Hymenoscyphus species.

The interpretation

The tree showed H. albidus and H. pseudoalbidus; were clearly separated. Japanese isolates were separated from a single common ancestor of European isolates, suggesting that Japanese isolates can be ancestral. This observation supports the idea that ash-dieback pathogen was originated from Asian isolates (Queloz et al. 2011 For. Path. 41:133-142, Pautasso et al. 2013 Biological Conservation 158:37-49.)

Acknowledgements

Tsuyoshi Hosoya from National Museum of Nature and Science, Japan kindly provided Japanese samples.
Steve Collin from Norfolk Wildlife Trust kindly collected some of the UK samples.

 

 

The top 100 ranked C. fraxinea candidate effector tribes

The contributors

Daniel Bunting (Nuffield student), Kentaro Yoshida and Diane Saunders at TSL.

The material

We used the potential C. fraxinea KW1 effector candidates identified in (http://oadb.tsl.ac.uk/?m=20130910).

The analysis

To select proteins with high likelihood as candidate effectors we focused on tribes that ranked highly in our scoring system. We used hierarchical clustering of the top 100 ranked C. fraxinea effector candidate tribes to group them further based on shared features.

 

Chalara_circos

Figure 1. The top 100 ranked protein tribes containing putative effectors. A. Combined score used to rank tribes based on their content of effector features. B. Score for number of members classified as secreted. C. Score for number of members identified as expressed during infection. D. Score for number of members with similarity to known fungal AVRs. E. Score for number of members with effector motifs or nuclear localisation signals (NLS). F. Score for number of members classified as repeat containing. G. Score for number of members classified as small and cysteine rich. H. Score for number of members encoded by genes with at least one flanking intergenic region > 10Kb. I. Score for number of members not annotated by PFAM domain searches. Red star indicates tribe that contains potential Nep1-like protein (http://oadb.tsl.ac.uk/?p=235).

Mining for putative effectors in the Chalara fraxinea KW1 genome

The contributors

Daniel Bunting (Nuffield student), Kentaro Yoshida and Diane Saunders at TSL.

The material

We used the C. fraxinea KW1 predicted proteome (data/ash_dieback/chalara_fraxinea/Kenninghall_wood_KW1/annotations/Gene_predictions/TGAC_Chalara_fraxinea_ass_s1v1_ann_v1.1) as a basis to mine for candidate effectors.

The analysis

  1. First, the predicted proteome of C. fraxinea KW1 was searched for potential secreted proteins using SignalP2 with parameters described in [1]. Transmembrane domain containing proteins and proteins with mitochondrial signal peptides were removed using TMHMM [2] and TargetP [3], respectively.
  2. We then clustered all proteins using TribeMCL [4], following the methods described in [5]. We identified clusters of proteins (known as tribes) that contained at least one secreted protein. These 593 tribes were then used for all further analysis.
  3. Next, we annotated the protein tribes for known effector features as described in [6].
  4. Finally, we assigned an e-value to each feature within a tribe using the method described in [7] in order to rank tribes based on their likelihood of containing effector proteins.

A spreadsheet that contains the above analysis is available at (data/ash_dieback/chalara_fraxinea/Kenninghall_wood_KW1/annotations/Effector_mining).

flow_diagram

Figure 1. Pipeline used to mine for potential effector proteins in C. fraxinea KW1 isolate. Programs are indicated in red.

References
1. Torto TA, Li S, Styer A, Huitema E, Testa A, Gow NAR, Van West P, Kamoun S: EST mining and functional expression assays identify extracellular effector proteins from the plant pathogen Phytophthora. Genome Res 2003, 13(7):1675–1685.
2. Krogh A, Larsson BÈ, Von Heijne G, Sonnhammer ELL: Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. J Mol Biol 2001, 305(3):567–580.
3. Emanuelsson O, Brunak S, von Heijne G, Nielsen H: Locating proteins in the cell using TargetP, SignalP and related tools. Nat Protoc 2007, 2(4):953–971.
4. Enright AJ, Van Dongen S, Ouzounis CA: An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res 2002, 30(7):1575–1584.
5. Haas BJ, Kamoun S, Zody MC, Jiang RHY, Handsaker RE, Cano LM, Grabherr M, Kodira CD, Raffaele S, Torto-Alalibo T, et al: Genome sequence and analysis of the Irish potato famine pathogen Phytophthora infestans. Nature 2009, 461(7262):393–398.
6. Cantu D, Seqovia V, MacLean D, Bayles R, Chen X, Kamoun S, Dubcovsky J, Saunders DGO, Uauy C: Genome analyses of the wheat yellow (stripe) rust pathogen Puccinia striiformis f. sp. tritici reveal polymorphic and haustorial expressed secreted proteins as candidate effectors. BMC Genomics 2013, 14:270.
7. Saunders DGO, Win J, Cano LM, Szabo LJ, Kamoun S, Raffaele S: Using hierarchical clustering of secreted protein families to classify and rank candidate effectors of rust fungi. PLoS One 2012, 7(1):e29847.