Comparison of transcriptome technologies in the pathogenic fungus Aspergillus fumigatus reveals novel insights into the genome and MpkA dependent gene expression

mRNA-Seq data set summary

To gain a deeper insight into the regulatory circuits of the MAP-kinase MpkA in A. fumigatus we extracted RNA from a null mpkA mutant and from the corresponding wild-type strain. We sequenced a total of six libraries, three biological replicates for the ΔmpkA strain and three for the wild-type strain. We obtained about 263 million paired end reads. Of these, 193 million (74%) were uniquely mapped against the A. fumigatus A1163 genome. Using a stringent RPKM cut-off of 10 reads per gene, 8172 genes (80% of the total annotated genes) were found to be transcribed in the wild-type strain we used. The coverage obtained was similar to that analysed in the dataset published by Gibbons et al.[12] for A. fumigatus (72% of the genes having an RPKM greater than 10). It was also similar to an mRNA-Seq study in A. oryzae[13], where the majority (83.4%) of the genome was covered by at least one read. Requiring coverage of at least 10 overlapping paired-end reads per base-pair, 70.4% (~20 million base pairs) of the genome was found to be expressed. This is remarkable, considering the fact that according to the A. fumigatus CADRE genome database, only ~50% of the genome is potentially transcribable. These findings indicated that the A. fumigatus transcriptome has been greatly underestimated.

UTRs prediction and annotation of new genes

So far, most studies have used mRNA-Seq data to identify transcriptional islands, which are consecutive areas of overlapping reads [13]. We additionally checked these Transcript Active Regions (TARs) for open reading frames. Since this approach is error prone, we employed ab initio gene prediction combined with mRNA-Seq data to improve the quality of gene prediction (Figure 1). This evidence-based gene prediction has the advantage over the classic ab initio prediction of finding gene structures by making use of the mRNA-Seq expression profiles. This combination incorporates intron-junction information derived from intron-spanning reads and allows the prediction of not translated regions (UTR) in a systematic manner.

Gene prediction still requires manual work to find appropriate adjustments for integrating data within a classic gene-predicting framework. By using mRNA-Seq data, we identified 185 new transcripts coding for putative proteins ( Additional file 1: Table S1 and Table S2). However, we could not exclude the possibility that the new transcripts are part of untranslated regions (UTRs) or even part of genes that need to be re-annotated. Among all the identified new transcripts, 44 of them have unknown functions, while 141 have orthologs in the genomes of other fungi.

Incorporating mRNA-Seq data into gene prediction also allowed the identification of UTRs ( Additional file 1: Table S2). However, UTR prediction remains a difficult task as only expressed regions provide evidence for UTR length. The method can therefore result in vague length predictions, especially for genes with low expression, due to insufficient coverage. Nevertheless, the identification of 5'UTR can help to find upstream open reading frames (uORFs) and led in our case to the identification of one potential uORF among the newly identified genes ( Additional file 1: Table S2). However, knowledge of UTR sequences in A. fumigatus remains poor. This is mainly due to the lack of available expressed sequence tag (EST) libraries. Our analysis highlighted the fact that some genes can have multiple or alternative transcriptional starting points. We found an average median length for 5'UTR of 308 and for 3'UTR of 97 nucleotides based on 9912 UTR sequences for each gene considered ( Additional file 2: Figure S1).

Comparison of transcriptome technologies

To better evaluate the performance of microarrays and mRNA-Seq, differentially expressed transcripts were investigated. We compared transcript levels in both the wild-type strain and the ΔmpkA mutant strain by repeating the experiment performed by Jain et al.[8].

Transcriptome data obtained by mRNA-Seq identified 2046 differentially regulated genes (log₂ fold change >2 and likelihood >0.99). For comparison, using the Febit microarray platform, 653 differentially expressed genes were identified (fold change >1.5) (see Additional file 1: Table S3, Table S4, and Table S5). KEGG analysis showed that common gene categories were enriched by the differentially expressed genes identified in the two datasets (Figure 2 and Additional file 1: Table S6, Table S7, and Table S8). The main difference between data sets was that mRNA-Seq gave a more detailed picture. The number of enriched genes for a specific category was substantially higher showing, for example, complete secondary metabolite gene clusters to be differentially regulated (e.g. the gene cluster bordered by genes AFUB_000840 and AFUB_000750). It is noteworthy that genes enriched in the KEGG category 0.1.1., i.e. genes involved in metabolic pathways, were differentially regulated in the mRNA-Seq data, but in contrast, were not enriched in the data set obtained using microarrays. The majority of these genes, in particular genes putatively involved in fatty acid metabolism, were down-regulated in the ΔmpkA mutant. This suggests that the RNA deep sequencing technique is more sensitive in finding differentially expressed genes than the microarray technique. This might be explained by the greater dynamic range covered by the RNA technique [14, 15].

For a more detailed analysis, we computed log fold changes for the ΔmpkA mutant vs the wild-type of the differentially expressed genes based on both microarray and mRNA-Seq data, including two biological replicates for each (Figure 3). Correlation of log fold changes of expressed genes was highest within the same technology. Correlation within mRNA-Seq replicates (Pearson’s correlation coefficient r = 0.96) was higher than that between the two microarray replicates (r = 0.53). This low correlation was mainly due to genes that only slightly changed their expression between conditions and are therefore more affected by the technological noise resulting from the technology used. A possible explanation is the low sensitivity of microarrays for low expressed genes that results in a smaller dynamic range compared to sequence data [15].

When comparing the global transcript fold changes between mRNA-Seq and microarrays the correlation amounted to 0.14-0.19 considering all expressed genes (Figure 3). This was not surprising because most genes are not expected to be differentially expressed, which consequently reduces the correlation. Since these genes are usually not of interest, we selected genes having the highest absolute fold change in both studies (10% quantiles of the total number of genes, intersection from mRNA-Seq and microarrays). The correlation value for this subset resulted in a large increase of both the within-array correlation (from 0.53 to 0.83) and the cross-technology correlation (between mRNA-Seq and microarrays, from 0.19 to 0.41). However, the 0.41 correlation is still lower than those found in other studies (e.g. Marioni et al. [16] found a correlation of 0.73 by comparing log fold changes between Affymetrix arrays and Illumina sequencing in experiments performed on human samples).

With the purpose of having additional values for evaluating the transcriptome data from the two technical platforms, qRT-PCR was performed in order to quantify transcript levels of genes that showed different expression levels in both studies (Figure 4, Additional file 1: Table S9 and Additional file 3 :Figure S2). Among the 15 transcripts investigated, 10 confirmed the direction of expression change detected by both transcriptome analysis techniques (with higher significance for mRNA-Seq data). Surprisingly, qRT-PCR did not reveal significant changes in some genes that showed opposite regulation in the analyses of the two techniques. For instance, qRT-PCR confirmed the mRNA-Seq results for the sidD gene, coding for a non-ribosomal peptide synthetase involved in siderophore production (AFUA_3G03420). Microarray data showed that genes expressed during low iron stress response were down-regulated in the ΔmpkA mutant strain. Previous experiments reversed microarray results, demonstrating that MpkA negatively regulates genes involved in siderophore production, indicating that many genes involved in iron depletion response are up-regulated in the ΔmpkA mutant strain (e.g. sidA or sidD) [8]. Current mRNA-Seq data confirmed our experimental hypothesis, by showing that genes involved in the iron starvation response are, in fact, up-regulated.

In total, the dataset compiled consists of seven wild-type measurements that were compared in a pair-wise manner. A pair wise scatter plot was realised by computing the Pearson and Spearman correlations for all pairs (Figure 5 and Additional file 1: Table S10 and Table S11). Spearman is a rank correlation technique that can capture some non-linear relationships and is less prone to outliers produced by highly abundant genes. As for Marioni et al. [16] and Feng et al. [14] the correlation between replicates for the same experiment ranged from 0.93 for the dataset of Bruns et al.[17] to 0.99 for the both mRNA-Seq studies, indicating high platform-reproducibility for microarrays as well as for sequencing data (data not shown). Therefore, it was sufficient to include only one representative replicate for each platform to obtain concise comparisons.

The transcriptome data reported showed that expression values measured for transcripts were strongly scattered. Cross platform correlation was quite low, ranging from 0.12 ([17]vs[18]) to 0.52 ([4]vs[18]). Correlation between the independent measurements (probes A and B) as calculated for the dataset of Jain et al. [8] was only 0.64, although all probes were exposed to the same samples. This indicates a strong technological bias of microarray data due to variation in hybridisation. At the same time, the correlation between the two Febit datasets for the combined probe measurements was quite high (0.89) pointing to a high conservation of technological bias in microarray platforms. In order to perform a general comparison between technologies, i.e., microarrays vs mRNA-Seq, the correlation of data within the same technology was used as a reference for the correlations between technologies. The highest correlation (0.68) among comparisons of mRNA-Seq data to microarray data was higher than the highest correlation between arrays (0.52, Roche vs. TIGR v1). This points to the fact that RNA sequencing data are less affected by technological bias, which is confirmed by the scatter plot comparing both mRNA-Seq datasets (Figure 5).

Comparison between A. fumigatus transcriptome and proteome

In order to analyse a possible correlation between transcripts and protein abundances as well as fold-changes, mRNA-Seq data, microarray data and the proteome data of the A. fumigatus wild-type and ΔmpkA mutant strain were compared. Looking at differential expression, 48 different proteins were identified by 2D-gel electrophoresis and subsequent MALDI-TOF / TOF, which showed a significant fold change in abundance (employing DIGE) in the ΔmpkA mutant compared to wild type. Gene enrichment and pathway analysis revealed that abundance of many proteins involved in metabolic processes and in oxidation/reduction processes is affected by MpkA ( Additional file 1: Table S12 and Additional file 4: Figure S4).

The regulation of protein levels was also compared with transcriptome data obtained from microarrays developed by Febit and used in the study of Jain et al.[8] and our mRNA-Seq data. As in earlier results (Figure 3, based on 9800 genes rather than 48) the fold change correlation between microarrays and mRNA-Seq was low (r = 0.14, Additional file 5: Figure S3). The measurements for transcriptome fold change often display low reproducibility because of difficulties in capturing the “true” fold change. For this reason, a low correlation between proteome and transcriptome technologies was assumed. However, the correlation between data obtained by microarrays and the proteome data was significantly higher (r = 0.28) than the correlation between mRNA-Seq data vs proteome (r = 0.19). As in the comparative analysis of the transcriptome data, the low correlation applies to all genes with a low fold change resulting in a higher relative impact of technological noise. Furthermore, it has already been shown by studies on human and mouse [19, 20] that the relationship between transcriptome and proteome data is quite complex and protein levels are greatly influenced by post-translational processing and inherent variations in stability. These factors result in low correlations.

In order to study absolute protein abundances, we reanalysed the A. fumigatus 2D proteome map generated by Vödisch et al.[21], resulting in the identification of 312 unique proteins. These proteins were also present in the different transcriptome data mentioned above ( Additional file 1: Table S13 and Table S14). In A. fumigatus, the relative abundances of the proteome map correlated surprisingly well with mRNA-Seq expression levels (r = 0.36 for our data, and r = 0.5 for data from [12], Figure 5) as well as with most microarrays, in particular the Febit arrays (r = 0.47–0.49). It is also noteworthy that the genes on the scatter plots (Figure 6) are more equally distributed for mRNA-Seq vs proteome comparison (left bottom corner) whereas, in particular for Febit arrays, a bias was observed towards genes that are highly expressed on the arrays but with low expression on the proteome-map. This is most probably due to cross-hybridisation. By contrast, this sort of bias was not observed in the NGS data.

To check for genes that are potentially subject to post-transcriptional regulation, we selected genes that were very abundant in the mRNA-Seq data and rare in the proteome-map (coloured green in Figure 6) or the inverse (coloured red). FunCat-analysis using FungiFun tools [22] indicated that the “green” genes (i.e. those genes that are highly expressed according to mRNA-Seq but little expressed according to proteome data) were not significantly enriched within a functional category (even after multiple correction), whereas the red ones (i.e. those genes that are little expressed according to mRNA-Seq but highly expressed according to microarray data) were over-represented within the tricarboxylic-acid pathway and amino acid metabolism. These processes constantly require large amounts of proteins that are therefore less subjected to transcriptional regulation.

cDNA library construction and sequencing

The A. fumigatus strain CEA17^KU80 strain was used for all experiments. The CEA17^KU80 strain was derived from the CEA17 wild-type strain [28].

Total RNA was isolated from wild-type and ΔmpkA strains cultured in Aspergillus minimal medium (AMM) for 16 h at 37°C [8], using the Qiagen RNeasy Plant Mini kit (Qiagen, Germany), according to the manufacturer’s instruction. Three biological replicates for each strain were collected. Total RNA was used for Illumina next-generation sequencing [29]. For library preparation 5 μg of total RNA per sample were processed using Illumina mRNA-Seq sample prep kit (RS-100-0801) following the manufacturer’s instruction. The cDNA libraries were sequenced using a Genome Analyser (GAIIx) in a paired-end approach with 2 × 36 cycles resulting in paired reads with a length of 36 nucleotides per read. Each library was sequenced on a single lane and ended up with around 30–40 million reads per sample. Sequence data were extracted in FastQ format and used for further analysis.

Mapping and normalising of transcriptome reads

All reads were mapped using TopHat 1.2.0 [30] against the A. fumigatus A1163 genome (which is a derivative of the A. fumigatus CEA17 strain) retrieved from Ensemble release 9. Parameters were set according to preliminary investigations, using a minimum intron length of 30 and a maximum of 4000. In agreement with the protocol used, from preliminary mapping a mean inner distance between mate pairs of 307 was found with a standard deviation of 110. To avoid multiple hits for a single read, the maximum multi-hit option was set to one.

Transcript expression levels were normalised by counting the number of reads per kilobase of exon region per million mapped reads (RPKM) [31].

To detect genes differentially expressed in the wild type and the ΔmpkA mutant, we applied the R package baySeq [32] which takes advantage of the three biological replicates for each condition. RPKM and differentially expressed genes were determined using the A. fumigatus A1163 genome annotated by CADRE as a reference.

Differentially expressed genes were categorised using GO enrichment, KEGG enrichment, and FunCat analysis using the FungiFun platform [22].

UTRs prediction

With the purpose of identifying putative new transcripts, we performed new gene annotation employing Augustus ab initio prediction [33] combined with hints generated from the mRNA-Seq reads. We present a new annotation including UTRs for A. fumigatus A1163. These hints consisted of intron junctions that were identified by TopHat as well as exon-parts revealed with the IRanges package which is part of the Bioconductor collection. The workflow for new gene annotation is reported in Figure 1.

Microarray data

The transcriptome and proteome data analysed are summarised in Table 1. The datasets were obtained in a pre-processed state i.e. background-corrected and normalised. Since the different microarray platforms and even the different A. fumigatus genomes have used different Gene-IDs (e.g. accession numbers starting with AFUB for the A1163 strain, and accession numbers starting with AFUA for the Af293 strain), all data were mapped according to the AFUA annotation, using ensemble mapping, to achieve comparability. The data were imported and analysed in the R / Bioconductor software environment employing the packages biomaRt and Biostrings.

Real-Time PCR

The qRT-PCR experiments to determine the amounts of selected putatively differently expressed genes were performed with StepOnePlus Real-Time PCR System (Applied Biosystems), using myTaq HS Mix 2x (Bioline) and Evagreen (Biotium). For each transcript analysed more primers were tested in tandem, in order to obtain primer efficiency close to 100% in all cases. To validate primers, efficiency standard curves were realised considering seven serial dilutions of genomic DNA of A. fumigatus in three replicates. All the primers selected for the qRT-PCR analysis are shown in Additional file 1: Table S15. The A. fumigatus actin1 gene (AFUA_6G04740) was chosen as the normalising gene. The cDNA was generated using RevertAid Premium (Fermentas). As negative controls, a vial was prepared without reverse transcriptase, and another one was prepared without cDNA.

Assessing differential expression

To correlate the mRNA-Seq data with the corresponding presence of proteins, the mRNA-Seq data set was compared with the proteome data published by Vödisch et al.[21]. For our study, the normalised protein spot abundances on the 2D gels were reanalysed using the most recent version of the software Delta 2D 4.3 (Decodon, Germany).

The DIGE (difference in-gel electrophoresis) technique was used to analyse protein samples of wild type and ΔmpkA mutant mycelium cultivated in AMM for 16 h at 37°C. Twelve samples (three independent biological replicates and two technical replicates for each strain as well as one mixed internal standard), were labelled with CyDye minimal dyes according to the manufacturer’s protocol with slight modifications (GE Healthcare Bio-Sciences, Munich Germany). We then labelled 50 μg total protein of each sample with 300pmol of CyDye DIGE fluorophores (dissolved in dimethyl formamide). Samples obtained from the wild-type and from the ΔmpkA mutant strain were labelled with either Cy3 or Cy5. A pool of all six samples (3 × wt and 3 × ΔmpkA) was prepared, labelled with Cy2, and used as a global internal standard. We carried out 2D gel processing, spot screening and MALDI-TOF / TOF analysis with an ultrafleXtreme mass spectrometer (Bruker Daltonics, Germany), as previously described [34] with slight modifications. FlexControl 3.3 software was used for data collection and flexAnalysis 3.3 for spectra analysis and peak list generation (Bruker Daltonics, Germany). Peptide mass fingerprint (PMF) and peptide fragmentation fingerprint (PFF) spectra were submitted to the MASCOT 2.3 server (Matrix Science, U.K.) searching the recent version of the NCPInr database but restricted to the fungi. Mascot search parameters were the following: fixed modification of cysteine thiols to S-carbamidomethyl derivatives, variable methionine oxidation, up to one missed cleavage and a peptide mass tolerance of 100 ppm.

Accession number

The RNA-Seq data were deposited in ArrayExpress (ArrayExpress accession: E-MTAB-1236).