Genome-wide transcriptional response of Trichoderma reesei to lignocellulose using RNA sequencing and comparison with Aspergillus niger

The wheat straw –induced transcriptome of T. reesei QM6a

The T. reesei genome is 33.9 Mb in size with 9,126 predicted genes [20]. Transcriptomes were sequenced from replicated independent cultures under 3 different sets of conditions: 1) after growth from conidia for 48 h in the presence of glucose as sole carbon source, a monosaccharide which represses expression of many genes involved in plant cell wall degradation, 2) 24 h after transfer of washed mycelia from 1) into media containing ground wheat straw as the sole carbon source to monitor the induction of genes involved in polysaccharide deconstruction and 3) 5 h after addition of glucose to the straw cultures from 2) to determine genes responsive to carbon catabolite repression. The ball milled wheat straw used in this study contained 37% cellulose, 32% hemicelluloses, 22 ± 0.1% lignin and was 25% crystalline [3]. Statistical tests [21–23] were applied to enable us to identify all genes which were significantly differentially expressed (p-value of <0.001 for all three tests) between the three conditions studied (see Additional file 1). RPKM values were calculated for each of the biological replicates as well as for the combined mapping of all replicates at 48 h, 24 h straw and 5 h glucose (see Additional file 1). The results shown in this study are from the combined mapping scores and only inductions showing a significant score in all statistical tests are discussed.

Expression of CAZy genes in T. reesei and comparison with A. niger

The degradation of plant cell wall carbohydrates is mediated by enzymes of four different classes: the carbohydrate esterases (CEs), the polysaccharide lyases (PLs), the glycoside hydrolases (GHs) and the auxiliary activities (AAs). These enzymes are classed, based on their primary amino acid sequence and related activity, into families in the Carbohydrate Active Enzyme database (CAZy) (http://www.cazy.org) [24]. Analysis of the T. reesei QM6a genome identified 22 CE-encoding genes representing 8 families, 5 PL-encoding genes representing 3 families, 195 GH-encoding genes, representing 49 families and 6 AA-encoding genes, representing 1 family [18]. The A. niger ATCC 1015 genome contains 25 CEs representing 9 families, 8 PLs representing 2 families, 239 predicted GHs representing 50 families and 7 AAs representing 1 family [25]. There are differences in the families of CEs, PLs and GHs encoded by the genomes of both fungi [18, 20, 25]. PLs are not as important as GHs and CEs for wheat straw degradation, as PLs mainly target pectin, a structure which is also degraded by enzymes of many GH families, including GH family 28 [26]. The family of AAs encoded by the genomes of T. reesei and A. niger (AA family 9) were formerly known as GH61s but were shown to be copper-dependent oxidases and have a different catalytic mechanism to the GHs [27]. Enzymes of AA family 9 play important accessory roles in enhancing lignocellulose degradation [28].

Carbohydrate esterases which play a role in lignocellulose degradation and which are encoded by T. reesei but not by A. niger belong to CE family 15. The genome of T. reesei encodes one CE family 15 glucuronoyl esterase [JGI:123940], also known as CIP2, which contains a cellulose binding module (CBM1) and which plays an important role in dissociating lignin from hemicelluloses through targeting the ester bonds between the aromatic alcohols of lignin and the glucuronic acid residues from the xylose backbone in hemicelluloses [29]. CE families 8 and 12 are present in A. niger but not in T. reesei[18, 25] and the genome of A. niger encodes 3 CE family 8 pectin methylesterases involved in the de-esterification of pectin [30], and two CE family 12 rhamnogalacturonan acetylesterases involved in the deconstruction of plant cell wall pectin [31].

Glycoside hydrolases involved in lignocellulose deconstruction and encoded by T. reesei and not A. niger belong to GH families 39, 115 and 45 and assist in the degradation of xylan (GHs 39, 115) and cellulose (GH 45) [18]. In A. niger, proteins from GH families 26 and 51, which are not encoded by T. reesei, are involved in the degradation of hemicellulosic mannan and arabinan residues [26]. Furthermore, the genome of A. niger encodes proteins of GH families 53 and 88, which are involved in the degradation of pectin [25, 26, 32], inulinases and invertases belonging to GH family 32 and which degrade polysaccharides containing fructose and sucrose [25].

After 48 h growth in glucose, CAZy gene mRNA represented 1.14% of total RNA in T. reesei (c.f. 3% in A. niger, Figure 1), with proteins from GH families 16, 18 and 72 (glucanases, chitinases and glucanosyltransferases) representing approximately half (45%) of the total CAZy mRNA (Figure 2). Thus in T. reesei, low levels of mRNA from genes encoding enzymes involved in the degradation of complex carbohydrates, including hemicellulose and chitin, are present when the fungus is cultivated in glucose-based medium. In this medium many of these enzymes are likely to be involved in cell wall remodelling during hyphal extension as high growth rates are achieved in the presence of glucose in T. reesei[33]. In contrast, in A. niger, transcripts from the glucoamylase glaA gene accounted for over 65% of total CAZy mRNA in glucose medium [3]. Induction of this gene in the presence of glucose and glucose-containing polysaccharides such as such starch has previously been described [34].

Although both fungi seem to use a similar array of GHs, they up-regulate the expression of different CE-encoding genes. Whereas transcripts from CE family 1 were most abundant in A. niger (~10% of total CAZy mRNA), transcript levels of genes encoding proteins in CE family 5 were highest in T. reesei (~5% of total CAZy mRNA). This is an agreement with previous microarray studies [18]. The genome of T. reesei encodes 3 acetyl xylan esterases and 1 cutinase, all belonging to CE family 5 [18]. The highest expression values were recorded for 2 acetyl xylan esterases [JGI:73632, JGI:54219] one of which contains a CBM1 module [JGI:73632]. Acetyl xylan esterases remove acetyl groups at O-2 and O-3 positions of the xylose chain in arabinoglucuronoxylans [5], a process which has been shown to significantly enhance subsequent hydrolysis of xylans and cellulose [37]. The genome of A. niger encodes 3 CE family 1 members: one acetyl xylan esterase, one feruloyl esterase and one un-defined esterase [25]. Expression values of the acetyl xylan esterase and the feruloyl esterase were very high in the presence of straw [3]. Feruloyl esterases cleave ferulic acid groups which are esterified to the 5′-OH of arabinofuranosyl groups (arabinose residues linked to O-2 or O-3 of xylose) and which can be covalently linked to lignin or other ferulic acid groups in xylans [5]. It appears that the enzyme mix secreted by A. niger aids in loosening the lignin-hemicellulose structure in addition to de-acetylating the xylan backbone in order to allow access of other CAZymes to the underlying hemicellulose and cellulose polysaccharides. Thus the bulk of GH and AA enzymes used to degrade straw are from the same GH and AA families in both organisms, whereas different CE family members are used suggesting that both fungi specialised also in the cleavage of different bonds found within plant cell walls.

Not surprisingly, transcript abundance of PL-encoding genes was very low in both fungi (~0.012% of total CAZy mRNA in T. reesei and ~0.5% of total CAZy mRNA in A. niger), confirming that PLs do not play an important role in wheat straw degradation (as mentioned above). Similar results were obtained through previous microarray studies [18].

T. reesei also highly induced the transcription of genes encoding proteins other than hydrolases, which have been proposed to be involved in enhancing cellulose degradation (Table 1). One such enzyme is the expansin-like swollenin, swo1 (Table 1), thought to play a role in the loosening of the plant cell wall by disrupting hydrogen bonds between plant polysaccharides, thus increasing cell wall area and access of hydrolytic enzymes (such as cellulases) to the underlying polymers [38]. Another enzyme, CIP1 (Table 1), which contains a CBM belonging to family 1, is thought to enhance cellulose hydrolysis [9]. Transcript levels of cip1 were also detected in the presence of sophorose and regulation of this gene is the same as for other well characterised cellulases (e.g. cbh1), indicating a potential role for this protein in cellulose degradation [39]. Our results show that genes encoding both enzymes are highly induced in the presence of straw and that they are regulated similarly to CAZy enzymes-encoding genes, suggesting that these enzymes could be key players in wheat straw degradation. Genes encoding CIP1 and swollenins are absent from the genome of A. niger.

Addition of glucose to the straw cultures exerted strong carbon catabolite repression of the CAZy-encoding genes, and CAZy transcript abundance decreased to 0.82% in T. reesei (1.4% in A. niger, Figure 1) of the total cellular mRNA, with members from GH families 16, 17, 67 and 72 (glucosidases, glucuronidase and glucanosyltransferases) (and GH15 in A. niger) being the most expressed CAZy genes under this condition in T. reesei (Figure 2).

Expression of non-CAZy genes in T. reesei and comparison with A. niger

Antisense transcription

Regulation of gene expression can also occur at the post-transcriptional level through regulatory RNAs. Natural antisense transcripts (NATs) are non-protein coding, fully processed mRNAs, which can partially overlap the protein-coding transcripts and which have many regulatory roles [51]. NATs have been found in several fungi including A. niger and N. crassa[3, 52]. To calculate the levels within our T. reesei transcriptomes, the number of reads corresponding to the non-coding strand was counted for each gene and AS (antisense) RPKM values were calculated in each condition. Approximately 1.82%, 1.47% and 2.79% of all reads in T. reesei were antisense reads when mycelia were grown in 48 h glucose, 24 h straw and 24 h straw + 5 h glucose respectively. Thus, AS transcription was detected in T. reesei and 630 genes had an AS RPKM > 1 in at least one condition (see Additional file 3). The 630 genes encoded proteins involved in a wide variety of cellular functions. The% of AS transcripts was similar in T. reesei when compared to A. niger[3].

The ratio of antisense:sense expression under glucose 48 h and straw 24 h conditions was calculated for these 630 genes in order to find genes with AS transcripts which change between the two conditions. Most of these genes have a low AS/S ratio on 24 h straw indicating that sense transcripts dominate over AS reads in this condition.

Confirmation of the presence of NATs in T. reesei was achieved by strand-specific PCR analysis for a gene [JGI:76852] with NATs (Figure 3). This gene is predicted to encode a secreted β-glucuronidase, belonging to GH family 2. Strand-specific RT-PCR confirmed the presence of spliced and non-spliced S and AS transcripts of different sizes (Figure 3B and C) in all three conditions as was previously reported for a gene containing NATs in A. niger[3]. Regulation at the post-transcriptional level presents an interesting area for further research.

Strains and growth conditions

T. reesei QM6a [53, 54] was used throughout this project. Conidia were produced from glycerol stocks of T. reesei grown on potato dextrose agar medium (PDA: 4.0 g/L potato extract, 15.0 g/L agar, 20.0 g/L) at 28°C. Conidia were harvested with 2 ml 0.01% (w/v) Tween 80. Liquid batch cultures were inoculated at a concentration of 10⁵ spores/mL.

Strains were cultured in 100 ml Trichoderma Minimal Media [TMM: 15 g/L KH₂PO₄, 5 g/L (NH₄)₂SO₄, 10 g/L of the respective carbon source, 0.005 g/L FeSO₄.7H₂0, 0.0016 g/L MnSO₄.H₂0, 0.0014 g/L Zn.SO₄.H₂0, 0.0037 g/L CoCl₁₂.6H₂0, 0.6 g/L MgSO₄, 0.6 g CaCl₂] in 250 mL Erlenmeyer flasks at 28°C, shaking at 150 rpm. Mycelia were grown for 48 h in 1% (w/v) glucose, after which they were removed by filtration through Miracloth (Merck), washed with double-distilled water (ddH₂O), and transferred to fresh media supplemented with the relevant carbon source at 1% w/v. Three different sets of conditions were distinguished: 1) growth from conidia for 48 h in the presence of glucose as the sole carbon source (48 h glucose), 2) 24 h after transfer of washed mycelia from 1) into media containing ground wheat straw as the sole carbon source (24 h straw) and 3) 5 h after addition of glucose to the straw cultures from 2) (5 h glucose). Transcriptomes were analysed from triplicate cultures of 48 h glucose and 24 h straw and from duplicate cultures of 5 h glucose.

Ball milling, sugar, lignin and crystallinity analysis of the wheat straw used in this study can be found elsewhere [3].

RNA extraction

Mycelia from each condition were snap-frozen and ground to a fine powder under liquid N₂ using a mortar and pestle. 100 mg of mycelial powder was used for RNA extraction, the procedure of which was described elsewhere [3]. Briefly, total RNA was extracted from mycelial powder using TriZol reagent (Invitrogen) according to manufacturer’s instructions. Extracted RNA was purified using the Qiagen RNeasy Mini Kit following the manufacturer’s instructions of the RNA clean-up protocol with on-column DNA digestion.

RT-PCR and qRT-PCR

The synthesis of cDNA was carried out as described by [3]. PCR reactions were performed using RedTaq DNA Polymerase (Sigma) and 1 μL of cDNA in a 20 μL reaction. PCR conditions were 30 cycles of denaturation at 94°C for 1 min, annealing at 60°C for 2 min and extension at 72°C for 3 min. Strand-specific PCRs (ssPCRs; as described in [3]) were carried out at an annealing temperature of 64°C using the primer pairs indicated in the Additional file 4.

qRT-PCR reactions were run using the same system as described by [3] and carried out for 40 cycles with denaturation at 95°C for 30 s and annealing at 64°C for 30 s and extension at 60°C for 60 s. Briefly, qRT-PCR reactions were run in triplicates per gene in each condition in a total reaction volume of 20 μL containing 11.0 μl Fast SYBR green master mix (Applied Biosystems), 0.11 μL/primer (175 nM final concentration), 2.2 μL of cDNA and 8.58 μL ddH₂O. Gene expression values were calculated against a standard curve of known genomic DNA concentrations. All primer pair sequences are listed in the Additional file 4.

RNA-seq and subsequent data analysis

Ribosomal RNA was degraded in 10 μg of total RNA using the Ribominus Eukaryotic kit (Invitrogen). SOLiD whole transcriptome libraries were made according to the SOLiD Whole transcriptome kit protocol (Applied Biosystems) and the library concentrations were measured with the Quant-it HS dsDNA assay kit (Invitrogen). Libraries were pooled to equimolar amounts (Invitrogen) and gel purified using 2% size-select E-gels to 200–300 bp (Invitrogen). Emulsion PCR (0.5 M final concentration of pooled libraries) and bead-based enrichment was done according to the SOLiD 4 Templated bead preparation guide containing library. Sequencing was carried out on a SOLiD 4 ABi sequencer platform according to manufacturer’s instructions to generate single fragment 50 bp colour space reads.

Using the BioScope version 1.3 Whole Transcriptome Pipeline (LifeTechnologies), reads from each SOLiD 4 libraries were initially filtered against sequencing library adaptors and other sequencing artefacts. Reads were then mapped independently to the entire unmasked and masked versions of the annotated genome assembly of T. reesei, and to defined gene transcript sequences (JGI T. reesei assembly version 2, annotation filtered models version 2.0). For mapping against the genome assemblies, it was also possible to align reads against a library of exon junctions derived from the exon coordinates detailed in the genome annotation. This allowed reads spanning exon boundaries to be included in alignment results and be recorded in BAM format. Subsequently, the reads corresponding to rRNA encoding gene regions were removed from the BAM records before further downstream analysis. The resulting BAM file of mapped reads from each sample against the unmasked and masked genomes were processed with HTseq-count [3] to generate read counts per gene from uniquely aligned reads. These counts were determined for both the sense strand only and to both strands (unstranded). Summary metrics for these results are shown in Additional file 5.

When comparing the read alignment metrics between the unmasked and masked genome assemblies a negligible difference of less than 1% is seen for total mapped reads, uniquely mapped reads and for reads mapped within annotated genes. This suggests that the masked regions of the genome do not correspond to transcribed genomic regions. The percentage of reads that were uniquely mapped to the unmasked assembly, as a proportion of total mapped reads, was 63.60% for all 48 h glucose replicates, 69.64% for all 24 h straw replicates and 52.95% for the 5 h glucose duplicates (Additional file 5).

For transcript sequence read mapping, less than half the total number of mapped reads compared to the mapping against the annotated genome was observed for all conditions. This can be explained by the difference in reference sequence space used in mapping. Transcript sequences were derived from the annotated exon coordinates within the genome sequence; therefore it would not be possible to map reads corresponding to transcribed regions outside of these defined coordinates. These unmapped reads could be attributed to incomplete annotation of the genome and/or to reads that are not mRNA coding. This is supported by a similar number of uniquely mapped reads to the transcripts within gene coordinates of the genome mapping when compared to the transcripts. The percentage of reads mapped uniquely to transcripts sequences, as a proportion of total transcript mapped reads, was 38.49%, 54.93% and 33.83% for the conditions 48 h glucose, 24 h straw and 5 h glucose respectively (Additional file 5).

For gene expression analysis the mapped read counts per gene calculated against the unmasked genome. Counts were determined for both sense strand only and to both strands (unstranded), as described previously. Antisense read counts per gene were calculated by subtracting the sense counts from the other. Read counts were then normalized to RPKM (Reads Per Kilobase per Million mapped reads) expression values for each gene [3] and visualised with the Integrative Genome Viewer (IGV 1.5) programme [3].

Differential expression values were determined using DEGseq [23] using sense read counts per gene for each experimental condition. DEGseq implements three independent statistical tests (Fisher’s Exact Test [21], Likelihood Ratio Test [22] and an MA-plot-based method with the Random Sampling method [23]).

Availability of supporting data

The raw and processed RNA-sequencing data sets supporting the results of this article are available in the NCBI’s Gene Expression Omnibus [55] repository under GEO accession number GSE44648 [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE44648].