New gene models and alternative splicing in the maize pathogen Colletotrichum graminicola revealed by RNA-Seq analysis

Quality analysis of RNA-Seq data

We totally obtained 1.31 billion reads of a length of 100 nucleotides. Of these, approx. 51 million quality-processed reads with an average length of 82 nucleotides mapped to the genome of C. graminicola. This significant difference can be explained by the relatively low amount of fungal biomass that is present in the host, especially during the early stages of infection. Since about 67% of the C. graminicola genome contains protein-coding sequence information, the average coverage of these regions was approx. 122-fold.

Analysis of the localization and distribution of the obtained reads in the genomic context revealed that among the matched reads 93% were uniquely mapped to the genome with their full read length allowing two mismatches. About 7% mapped to multiple locations in the genome, of which 96% matched 2 to 10 different locations, and 4% to more than 10 different locations. RNA-Seq data showed that 82.5% of the reads that matched unique regions could be mapped to known coding/UTR regions, 15.9% to intergenic/unknown UTR sequences, and 1.6% mapped to genomic regions encoding rRNA and tRNA (Figure  1A). Of all reads that mapped to multiple locations, 56.8% mapped to genomic regions encoding rRNA and tRNA, 23.3% mapped to known UTR/protein coding regions, and 19.9% mapped to intergenic regions and yet unidentified UTRs (Figure  1B).

An indispensable prerequisite for reliable RNA-Seq data is the sufficient coverage of transcripts, preferably by a large number of reads that are evenly distributed across the target sequences [17]. Although 3-prime ends of transcripts were found to be represented by a somewhat higher quantity of reads, the overall coverage shows a largely even distribution and sufficient coverage of the coding transcript sequences by RNA-Seq reads (Figure  1C).

Altogether, transcripts corresponding to 11,732 (97.7%) of the 12,006 genes predicted in the current annotation of the C. graminicola genome could be detected. This number contrasts the 2,766 genes that had been previously validated by EST contigs [9]. About 9,675 (80%) of all predicted genes were found to be covered in full-length by RNA-Seq reads. Reads corresponding to about 87% of C. graminicola genes were detected at approx. 2.5 million uniquely mapped RNA-Seq reads. Subsequently, gene coverage increased only moderately despite the increasing sequencing depth (Figure  1D).

Further quality evaluation was performed in order to assess the reproducibility of biological repeats. Merged data sets of biological repeats 1 vs. 2 and 2 vs. 3 show consistent distribution of data points in the scatter plots (Figure  1E, 1F).

These results clearly demonstrate the good quality of the obtained RNA-Seq data as proven by the extensive coverage of hitherto identified genes, sufficient depth of the sequencing approach and the broad diversity of RNAs transcribed during pathogenic development of C. graminicola, which is a prerequisite for both the identification of a large number of nTARs and improvement of the annotated gene models.

Identification of novel gene models and non-coding RNAs

In C. graminicola, a total of 12,006 genes, encoding proteins of more than 100 amino acids, have already been predicted by gene-finding software tools [9]. However, 15.9% of the mapped reads fell within intergenic/unidentified UTR regions. A substantial part of these reads correspond presumably to non-coding RNAs and transcripts of as yet unrecognized genes or exons of incorrectly annotated genes.

Altogether, we found 819 nTARs, of which 430 harbor introns indicating that they could be genes. In addition, about 60% of all nTARs were longer than 500 bp (Additional file 1). Conceptional 6-frame translations employed to search the non-redundant protein database of NCBI (ftp://ftp.ncbi.nlm.nih.gov/blast/db/) and the protein domain database of the Sanger-Institute (ftp://ftp.sanger.ac.uk/pub/databases/Pfam) revealed the presence of coding sequences with similarity to known proteins in 146 nTARs (E-value <10^-8), of which 31 displayed a conserved protein domain (E-value <10^-3, Figure  2A, Additional file 2). Notably, 74% of these putative proteins were found to have orthologs in other Colletotrichum species (Figure  2B).

The majority of nTARs could not be assigned to particular genes, probably due to the fact that they either lack conserved orthologous genes in other species or these genes have not been annotated so far. Alternatively, these TARs could correspond to non-coding or regulatory RNAs. In order to identify non-coding/regulatory RNAs we compared database entries of known non-coding RNAs with nucleotide sequences of the C. graminicola nTARs. Despite the enrichment of mRNA sequences carrying terminal poly(A)-tails during sample preparation we found 5 nTARs with similarity to non-coding RNAs, i. e. 2 transfer RNAs, 1 spliceosomal RNA,1 small nuclear RNA and 1 ribosomal RNA (Figure  2A, Additional file 3). None of the nTARs shared similarity with database entries of known regulatory RNAs. One reason for this could be that regulatory RNAs are species-specific, which would prevent identification of similar RNAs in public databases containing predominantly entries of mammals and Drosophila melanogaster.

Though the function of the majority of nTARs remains elusive, 106 of the nTARs that were found to have orthologs in other fungal species display intact full-length ORFs (Additional file 2). These results convincingly show that a substantial proportion of nTARs are most likely transcripts of protein-encoding genes rather than pseudogenes or regulatory RNAs. As a consequence, a significant number of genes, some of which with putative functions in virulence or pathogenicity, have not been recorded in the current version of the genome annotation. For example, an nTAR designated COGR_14733.1 (Figure  2C) encodes a protein with assumed effector function as evidenced by the occurrence of an ortholog with putative effector function in C. higginsianum [18]. Transcripts of COGR_14733.1 were detected exclusively between 24-120 hpi (Figure  2D), which could be indicative of a role in the switch from the biotrophic to the necrotrophic lifestyle.

Identification of incorrect gene models

Following mapping and assembly of the RNA-Seq reads, we conducted a thorough manual inspection of the currently annotated genes using the integrative genome viewer [19] that allowed simultaneous survey of mapped reads, genome annotation of the Broad Institute and a genome annotation based on the algorithms of AUGUSTUS, a gene prediction software dedicated to the analysis of fungal genomes that has been trained using the supplied Fusarium graminearum data set [20]. Based on the coverage of the genomic sequence by RNA-Seq reads the following categories of incorrectly annotated features have been identified:

• I.

  Incorrect intron-exon boundaries: This type of incorrect annotation was identified either by reads that mapped partially to annotated intron regions or partial missing coverage of annotated exon regions. Prerequisite for the correction of these errors was a sufficient read coverage of the genes to be re-annotated. In comparison with the current annotation we found 434 introns that were too long and 45 that were too short (Figure  3A,B; Additional file 4).

• II.

  Incorrect exons/introns: The procedure that led to the identification of this type of error was similar to the approach chosen to identify errors of type I with the exception that only absent or redundant full-length introns were captured. At least one intron and one exon needed to be added to the sequence of 111 genes whereas one intron and one exon was removed in 192 genes (Figure 3A,C; Additional file 4).

• III.

  Incorrectly split or joined genes: This category refers to adjacent genes that were incorrectly joined in the annotation of the Broad Institute or individual genes that were incorrectly split due to misinterpretation of the genomic sequence data. Incorrectly joined genes were identified by the lack of read coverage in the sequence of the putative gene which could not be explained by the occurrence of an intron. Vice versa, incorrectly split genes could be found by a largely uniform read coverage throughout the sequence of the incorrectly annotated gene. In both cases, newly annotated genes were subjected to Blast search in order to confirm the identity of the new gene by sequence comparisons with orthologous genes of closely-related fungal species. In total, 31 incorrectly joined genes and 52 incorrectly split genes could be identified (Figure 3A, 3D; Additional file 4).

• IV.

  Incorrect genomic sequence: By comparing the genomic nucleotide sequence provided by the Broad Institute with our RNA-Seq and genomic sequence data we were able to detect sequencing errors or point mutations that led to incorrect gene models (Figure 3A). Among the 81 aberrations identified, 40 did not alter the genetic information of coding sequences, i. e. they were silent nucleotide substitutions, insertions or deletions within intergenic regions and nucleotide substitutions located either upstream or downstream of coding sequences or within TARs (Additional file 5). The gene structure or sequence of the translation product was affected in 41 genes by nucleotide substitutions leading to missense, loss of the stop codon and insertions or deletions of nucleotides in coding sequences (Additional file 4).

In some cases, these errors had fatal effects with respect to the quality of the algorithm-based gene annotation. For example, an insertion of adenine at nucleotide position 948 of gene GLRG_08997 encoding a cyclohexanone monooxygenase resulted in a frameshift and hence stop of translation. Automated annotation by the Broad Institute led to the incorrect insertion of a short intron in order to solve the conflict. The correct ORF was restored by removing the adenine (Figure  3E) as proven by blastp homology search against the NCBI database. Altogether, incorrect annotations of 41 genes originating from substitutions, deletions or insertions of single nucleotides could be improved by manual revision of gene models.

Using RNA-Seq data, a total of 906 incorrect gene models identified in the current genome annotation were successfully revised, which will significantly contribute to both increased quality and reliability of the genomic data of C. graminicola.

Identification of alternative splicing (AS)

Alternative splicing describes the regulated processing of mRNAs that leads to the expression of multiple proteins encoded by a single gene, which is an inevitable prerequisite for diversity and functional complexity of eukaryotic proteomes [21–23]. In order to assess alternative splicing in C. graminicola we inspected the re-annotated genome with special emphasis on reads that match intronic regions. Candidate genes for alternative splicing were identified by a significant read coverage within intron sequences as compared to the read coverage of the surrounding exons. We found 75 putative AS events in 69 genes, including exon skipping, intron retention and alternative 5′ or 3′ splice sites (Figure  4A, Additional file 6).

Intron retention appears to be the most prevalent form of AS in C. graminicola, accounting for about 56% of all AS events. However, analysis of the alternate ORFs revealed premature termination of translation in most of the candidate genes irrespective of the AS mechanism. There was no evidence for a stage-specific exchange of entire domains, a characteristic effect of AS in higher organisms that leads to functionally diverse proteins.

The majority of the AS candidate genes that have been characterized are involved in transcription or translation. Further putative functions include maintenance of DNA, signal perception or transmission and cellular structure or transport (Figure  4B).

In total, about 0.6% of all C. graminicola genes were estimated to undergo AS, though it must be emphasized that delayed processing of premature mRNAs may mimic alternative splicing [24, 25]. In line with this, detection of alternatively spliced RNAs by RT-PCR, which has been frequently used, does not appear to be an appropriate technique providing sufficient evidence for alternative splicing.

Identification of untranslated regions (UTRs)

We identified the 5′- and 3′-boundaries of transcripts by searching for a sharp decline of RNA-Seq reads signals at both ends of the annotated transcripts. UTRs of genes that overlap with UTRs of adjacent genes were excluded from further analysis. Altogether, 5′- and 3′-UTRs of 4378 genes could be retrieved (Additional file 7). The average length of UTRs in the genome of C. graminicola was found to be 272 nucleotides. UTRs of the 5′- and 3′-regions have an average length of 229 and 316 nucleotides, respectively. Among the identified UTRs, 893 harbor intronic sequences, of which 535 are located in the 5′-UTRs and 363 in the 3′-UTRs.

Screening for RNA editing

Comparing the genomic sequence and RNA-Seq data obtained from C. graminicola strain M2 offered the opportunity to identify nucleotide modifications that occurred due to RNA editing.

A genome-wide search for SNPs, deletions and insertions (INDELS) of single nucleotides yielded 81 of these modifications, however, careful manual examination revealed that the majority of putative RNA editing-based substitutions or insertions were restricted to repetitive sequences of a particular nucleotide or located near intron splice sites. Therefore, we concluded that single nucleotides were erroneously introduced or deleted by DNA polymerases used in the sequencing reactions rather than by RNA editing mechanisms. In addition, several SNPs appeared to be caused by misalignment of RNA-Seq reads to the genomic DNA sequence. In order to verify the nucleotide sequence of regions harboring putative SNPs we performed local re-sequencing of the genomic DNA. The obtained results strongly suggest that no RNA editing occurs in C. graminicola.

Colletotrichum graminicolastrain, culture conditions and maize infection

The wildtype strain CgM2 of C. graminicola (Cesati) Wilson [teleomorph Glomerella graminicola (Politis)] was provided by R. L. Nicholson, Purdue University, IN, USA. CgM2 was cultured on oat-meal-agar plates at 23°C under near UV light (Philips TLD36W/08, Hamburg, Germany).

Conidia were harvested from 2-week-old plates in distilled water and washed three times. The final concentration of conidia was adjusted to 10⁶/ml in 0.01% Tween 20 prior to inoculation of maize plants or leaf segments.

For RNA-Seq, maize plants (Zea maize cv. Golden Jubilee) were grown in individual pots for 14 days at 23°C. The 3^rd leaf was carefully fixated in a horizontal position and inoculated with 10-μl droplets containing 10⁴ conidia in 0.01% (v/v) Tween20. After 24 h in a moisture chamber at 25°C, plants were transferred to an environmentally-controlled greenhouse cabinet and further incubated at 25°C. Samples corresponding to three biological replicates were taken at 0, 12, 24, 48, 72 and 120 hours postinoculation using a cork borer (8 mm diameter) and immediately frozen in liquid nitrogen. The infection progress was observed microscopically and anthracnose lesions were photographed at 5 days postinoculation.

Preparation of RNA and Illumina sequencing

Frozen leaf discs (500 mg per sample) were homogenized in 2 ml reaction tubes using pre-cooled steel beads (3 mm) and a tissue lyser (Qiagen, Hilden, Germany) at 30 Hz/60 s. Four hundred and fifty μl buffer RLT (Qiagen, Hilden, Germany) was added to the frozen powder. RNA extraction was performed using the Plant RNA Kit (Peqlab, Erlangen, Germany) according to the manufacturer’s instructions. Quality and quantity of the extracted RNA was examined using Nanodrop 1000 photometer (Peqlab, Erlangen, Germany) and Bioanalyzer 2100 (Agilent, Böblingen, Germany). RNA samples were further processed and subjected to high-throughput sequencing by Illumina HiSeq^TN2000 at ServiceXS (Leiden, Netherlands).

Processing of RNA-Seq reads and mapping

The C. graminicola CgM2 genome and gene information (version 11/2012, bioproject PRJNA37879) were downloaded from http://www.broadinstitute.org/annotation/genome/colletotrichum_group/MultiDownloads.htm. Reads containing sequencing adapters and reads of low quality were identified by fastqc (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc). Quality-processing of raw reads included Phred score-guided trimming (Phred score >20) of reads using the software fastx_trimmer of the FASTX toolkit (http://hannonlab.cshl.edu/fastx_toolkit). Elimination of adapter sequences was achieved by removing 12 nucleotides of the 5′-ends of all reads using TrimGalore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore). All further analysis was done with quality-processed reads.

RNA-Seq reads were mapped against the genome of C. graminicola with Tophat version 2.0.8 [39] using default options except for the minimum intron length which was set to 30 nucleotides. Each time point and biological replicate of the maize infection approach was mapped independently, and individual datasets were subsequently joined prior to further analysis.

In order to evaluate the reproducibility of RNA-Seq data obtained from biological repeats we compared the read number per gene of samples belonging to three independent infection series (0-120hpi) using SAMtools [40], htseq-count [41] and R [42].

The average coverage of all cds by RNA-Seq reads was determined using RSeQC [43]. Genomic features were analysed using featureCounts [44] of the software package Subread version 1.4.5. Coverage and distribution of transcriptional active regions was analysed using BEDTools [45] together with custom bash scripts.

A comprehensive scheme illustrating sample preparation and workflow of the in silico analysis is provided in the supplement (Additional file 8).

Transcript assembly and identification of incorrect gene models

Initial transcript assembling was done as described [46]. Additional transcript assembling and gene prediction using introns as hints was performed by AUGUSTUS [20]. Predicted gene models were manually inspected using the Integrative Genomics Viewer (IGV) [19] that enabled simultaneous display of gene models predicted by AUGUSTUS, Cufflinks and Broad Institute together with RNA-Seq reads. For each gene model it was decided whether one of the above annotations was supported by RNA-Seq reads. A manual curation of gene models was performed if none of the gene models was in agreement with our RNA-Seq data.

Identification of novel transcriptionally active regions

Novel transcriptionally active regions (nTARs) within intergenic regions of the current genome annotation were detected essentially as described in Trapnell and associates [46]. Briefly, nTARs were identified by comparing the current broad annotation with the annotation obtained from the cufflinks pipeline (filtering class code u) using cuffcompare (option -r).

Further characterization of nTARs was performed by blastx search against the non-redundant protein database hosted by NCBI. Subjects with an E-value <10e^-8 were considered to be putative orthologs.

Nucleotide sequences were translated using TransDecoder (http://transdecoder.sourceforge.net) and analyzed with respect to the occurrence of functional protein domains in the Pfam-A database [48] using HMMER 3.1b1 (threshold parameter -cut_ga) [47], defining an E-value <10e^-3 as significant hit.

For further characterization, databases fRNA [49] and NONCODE v3.0 [50] both harboring collections of non-coding RNAs were searched using nTAR sequences as input.

Analysis of alternative splicing (AS)

The improved genome annotation was used to screen for AS in C. graminicola. Significant coverage of intron sequences by RNA-Seq reads or locally decreased read numbers in exons were considered to be indicative of AS. Alternative splice sites were required to be supported by a read number of at least 25% uniquely mapped reads compared to the read number of the surrounding exons. Putative AS events were grouped into 5′- or 3′- AS, intron retention and exon skipping. Transcripts displaying insufficient read coverage (read number <30) at junction sites and putative AS events that would lead to very short proteins due to the occurrence of an early stop codon in the alternate ORF were excluded from further analysis.

Annotation of UTRs

UTRs were annotated based on the predicted gene models of the current genome annotation and RNA-Seq sequence information using BEDTools [45]. To improve signal-to- noise ratio the number of aligned reads was doubled. Positions covered by 0, less than 8 and less than 15 reads had been removed in three independent calculations. The resulting fragments were annotated on the basis of the current genome annotation and fragments containing sequences of adjacent annotated genes were removed. UTRs were considered to be correct if two out of three calculated UTR models resulted in a similar UTR length.

Screening for RNA editing

Comparative analysis of the genomic sequence of the reference strain CgM2 of C. graminicola and our RNA-Seq data was performed in order to identify nucleotide modifications based on RNA editing. Read duplicates in the RNA-Seq data set and in the genomic sequences were removed using PicardTools (picard.sourceforge.net). Nucleotide substitutions and insertions/deletions (INDELS) were detected using SAMTools [40] and snpEff [51]. Candidate sites were required to be covered by more than 30 reads with at least 25% of reads possessing the alternate allele. Each genomic position harboring a putative nucleotide substitution or INDEL was carefully examined with respect to sequence errors or misalignments. Insertions or deletions of one nucleotide in oligo-nucleotide repeats in RNA-Seq data were considered as amplification bias and excluded from further analysis. Differences in the nucleotide sequence between RNA-Seq data and genomic sequence data that could not be explained by misalignments were further analyzed by re-sequencing the corresponding genomic DNA.

Data deposition

Illumina sequencing data were deposited in the sequence read archive at NCBI (http://www.ncbi.nlm.nih.gov/sra) with accession number SRX516617. The annotation file can be downloaded at http://www.landw.uni-halle.de/download/.