Transcriptomic response of the mycoparasitic fungus Trichoderma atroviride to the presence of a fungal prey

Characterization of the EST database

We constructed EST libraries of T. atroviride from four different cultivation conditions: (a) confrontation with two plant pathogenic hosts (Botrytis cinerea, Rhizoctonia solani) on agar plates in the dark, using an mRNA extraction time point when their hyphae were only 1-2 mm apart; (b) a respective control consisting of only T. atroviride growing in the dark and not sporulating; (c) cultures illuminated by blue light to induce conidiation [13]; and (d) mycelia growing on plates in the dark and subjected to stress provoked by mechanical injury, a method known to lead to sporulation independent of light [13]. The rationale of this approach was that genes which would be uniquely expressed under (a) could so be distinguished from those required for filamentous growth and sporulation by T. atroviride.

After editing the mean length of the single read (see Materials and Methods), sequences ranged from 0.57 - 0.62 kb. A summary of the properties of the EST database used is given in Additional File 1. The representativeness of the libraries was confirmed by the high degree of transcript diversity (diversity index ranging from 60 to 75%). No ESTs representing rRNA or prokaryotic contaminations were found.

Functional analysis of the EST libraries

Since an annotated genome sequence was available for this study, we developed the following approach for analysis of the ESTs: 27658 ESTs from all four libraries (cf. Additional File 1) were incorporated into a custom-made BLAST database, and the database then queried (BLASTX) with a total of 6889 genes from the T. atroviride genome database for which at least a putative function could be predicted and which were members of one of the 21 KOG groups (5 - 779 protein sequences/KOG group; Additional File 2). The respective protein models were extracted as FASTA files from the T. atroviride genome portal [12]. 9478 ESTs, showed positive matches. A comparison of the ESTs isolated under each of the four conditions illustrated that they were derived to approximately equal portions of ESTs from plate confrontations (2882), mycelial growth (2913), and light induced sporulation (3003). The number of hits from ESTs from plates subjected to mechanical injury were significantly lower (680), possibly due to the stress evoked by this growth conditions, which consequently had a negative impact on RNA quality, and they were thus excluded from most of the comparisons.

The 9478 ESTs represented 2734 genes, thus indicating that 39.7% of the queried genes were indeed detected in our EST collection. We examined to which of the 21 KOG classes the detected genes belonged and found that most KOG classes contained a similar proportion of expressed genes (29 ± 7%). Only three classes (J, translation; U, protein secretion; and Z, cytoskeleton) exhibited a strongly increased percentage of expressed genes (59, 53 and 51%, respectively). EST-density (the number of individual ESTs per average single identified gene) was highest in class J, too (9.1) but also significantly elevated in class C (energy metabolism; 4.2), indicating that these two classes may contain the most abundantly expressed genes (Additional File 2).

Of the total 2734 genes, for which ESTs were found, roughly 60% were found under only a single cultivation condition (Fig. 1). The remaining 1082 genes were retrieved from two or more conditions. Among these, 326 genes were found under all three conditions (mycoparasitism, mycelial growth and light-induced sporulation), and combinations of each two of these conditions were found for 175-192 genes, respectively. Combinations involving the genes expressed under mechanical injury also displayed roughly equal percentages but significantly lower numbers (35-41), which is however consistent with the other data when the lower number of total ESTs from this growth condition is taken into account. ESTs which were found under at least two conditions were particularly found in the KOG groups posttranslational modification, signal transduction, protein synthesis and energy metabolism (Table 2).

Genes involved in sporulation are overexpressed under light-induced sporulation

As a prove of-principle that our dataset allows identification of genes specifically induced under a certain condition, we tested the expression of a set of 15 genes which are known to be involved in sporulation in Neurospora crassa[14], and for which orthologues were identified in T. atroviride (Kubicek et al.; ms in preparation). We found a total of 26 transcripts for 10 of them. These ESTs were distributed between mycoparasistism, mycelial growth, light-induced conidiation and mycelial injury, respectively, in a ratio of 4:4:14:4 (Additional File 3). When normalized to the total number of ESTs/growth condition, this shows that ESTs involved in sporulation (light-induced conidiation; mycelial injury) are ca. 3.5- and 4.3-fold overrepresented under the two conditions leading to sporulation, even though they were represented by different numbers of ESTs. This demonstrates that the number of ESTs in our sample correlates with their expression under the respective condition. The presence of an increased number of ESTs under one condition thus indicates an upregulation of the respective gene.

Constitutively strongly expressed genes of T. atroviride

27 genes made up for 1034 of the 9478 ESTs found, and they were found in comparable numbers under the three major conditions (Additional File 4). The percentage of these strongly and constitutively expressed genes (10.8%) is slightly lower than reported for other Trichoderma EST datasets (15-17% for the 25 most abundantly expressed ESTs; [11, 15]). Some of these genes encoded glycolytic enzymes, hydrophobins or histones that were also found in other studies. Some strongly expressed genes that we detected, however, were unexpected. Intriguingly, the most abundantly expressed gene was cpc1, which encodes the regulator of general amino acid control ("cross pathway control") in fungi [16]. Although its transcript abundance was much higher under conditions of mycoparasitism, its nevertheless significant abundance under the other conditions made us to refrain from calling it "mycoparasitism-specific" (see also below) but we note that it is strongly expressed under this condition.

Genes overexpressed under mycoparasitic conditions

Our EST collection (2734 genes) contained only 28% of all genes present in the T. atroviride genome (11100 genes) and we therefore considered the presence of only single ESTs (albeit clearly proving the expression of the gene under this condition) to be insufficient for claiming specific induction in comparison to other conditions. To identify those genes from our EST collection that are significantly overexpressed under mycoparasitic conditions, we therefore looked for genes which were represented either by at least 3 ESTs under mycoparasitic conditions and by none under the others, or by 4 or more ESTs and thereby accounting for at least 65% of the total number of ESTs from the respective gene. Additional File 5 show s the results from this screening, grouped according to KOG terms: 442 ESTs, corresponding to 66 genes, were found to fulfil this criterion and thus be strongly overexpressed under mycoparasitic conditions. They were distributed over 18 KOG groups, but were most abundant in those comprising posttranslational processing (O) and amino acid metabolism (E): 159 of the total number of ESTs from the KOG group posttranslational processing (183) and 70 out of 159 ESTs assigned to amino acid metabolism belonged to the mycoparasitism-EST library. Conspicuously, several heat shock factors and tRNA synthases, belonging to the posttranslational processing group according to the KOG classification, were present in high numbers and may indicate a significant shift towards response to environmental stress (see Discussion). In addition, lipid metabolism also accounted for a significant portion of these EST (24 from 107). However, apart from these quantitative considerations, a number of genes (e.g. components of signal transduction, proteases) which comprised only a small portion of their KOG groups were also found and will be emphasized in the Discussion.

RT-PCR confirms upregulation of a subset of genes whose ESTs are abundant during mycoparasitism

To test whether the EST abundance in fact reflects gene expression patterns under conditions of confrontation of T. atroviride with its prey, we performed reverse transcriptase PCR (RT-PCR) with 12 genes up-regulated during mycoparasitism and in addition with cpc1, which was according to our EST analysis abundant during mycoparasitism but also during mycelial growth. Two housekeeping genes served as a control (gpd1, tef1). To verify that an upregulation of the selected genes could be attributed to the presence of the host fungus irrespectively of the medium used, plate confrontation assays were set up on potato dextrose agar (PDA), whereas a minimal medium with 0.3% carbon source had been used for generation of the ESTs. Using Rhizoctonia solani as a host fungus, samples were taken before contact (5 mm distance of the mycelia), just at contact, and during overgrowth. It should be noted that for the ESTs mycelia were harvested at a growth stage that was between 'before contact' and 'contact' of this experiment. Expression patterns of the selected genes were monitored. Growth of T. atroviride alone on PDA plates, was used as control growth condition, corresponding to 'mycelial growth' in our EST collection.

In the gene expression analysis (Fig. 2), eight of the twelve putative mycoparasitismspecific genes showed a behaviour consistent with the EST data and were upregulated during mycoparasitism, whereas their transcripts were either absent or of much lower abundance in the control. The four other genes, however, also showed a high level of transcription in the control: hsp26-3, wsc1, pdc and gsy1 (although the latter was still more abundant at the onset of mycoparasitism), indicating that their upregulation is not mycoparasitism-specific. Transcripts of cpc1 were detected under all tested conditions. Therefore, our threshold for detection of mycoparasitism-specific genes appears to have mainly identified genes which are indeed actually upregulated during this process.

Reconstruction and comparative analysis of a metabolic network

In order to obtain a global understanding of the metabolic changes of T. atroviride during the onset of mycoparasitism, we also used the genome sequence data to reconstruct a metabolic network and subsequently attempted to identify pathways that were differently expressed at the onset of mycoparasitism. Only enzymes with standardized EC numbers were thereby considered, which yielded 832 EC numbers in total. From these, a bipartite metabolic network containing 2013 reactions and 2198 metabolites was constructed. Subnetworks containing a portion of these reactions were then extracted from the whole network to address the aforementioned differentially expressed genes and KOG groups in a metabolic context.

The first intriguing finding - was that a significant number of amino acid biosynthetic pathways, especially those for common amino acids were significantly up-regulated (the both number of ESTs was more than two-fold higher than under other conditions). These pathways included (Additional File 6): urea cycle and metabolism of amino groups; glutamate metabolism; methionine metabolism; alanine and aspartate metabolism; valine, leucine and isoleucine biosynthesis; lysine biosynthesis; arginine and proline metabolism; cysteine metabolism. In addition, the biosynthesis of several aminoacyl-tRNAs was also enhanced. Amino acid catabolism (e.g. that of valine, leucine and isoleucine degradation), on the other hand, was obviously downregulated. These data are supportive of our data from EST evaluation and suggest that mycoparasitism is associated with a state of amino acids starvation. We also found an upregulation of pyridoxalphosphate biosynthesis which is consistent with an increased amino acid biosynthesis. Metabolism of sulphur containing compounds was also remarkably up-regulated, and it is important to note that its genes were not expressed under any of the other conditions. Analysis of the lipid metabolism showed that the degradation pathways of lipid metabolism were slightly up-regulated. In addition, biosynthesis of pantothenate was also increased, consistent with its role as a constituent of coenzyme A required for fatty acid catabolism. Another upregulation was found for genes involved in sphingolipid metabolism. Biosynthesis of unsaturated fatty acids pathway was absent. Also, aminosugar catabolism was significantly enhanced under mycoparasitic conditions, which would be in accordance of an attack of T. atroviride on the prey's cell wall.

Cultivation conditions

T. atroviride IMI206040 was grown on Vogel's minimal medium [53] with 0.3% (w/v) glucose as carbon source. For confrontation assays agar plates were overlaid with cellophane and inoculated with an agar plaques of the plant pathogen (Rhizoctonia solani, Botrytis cinerea), and T. atroviride in approximately 5 cm distance from the pathogen. Strains were grown in complete darkness at 28°C and T. atroviride was harvested under red safety light at pre-contact (1-2 mm distance of the mycelia) and immediately frozen in liquid nitrogen. For light induced conidiation, T. atroviride cultures were grown in the dark for 48 h at 27°C on PDA plates and used as pre-inoculum. Mycelial plugs (0.5 cm diam.) were taken from the colony growth front and placed on the centre of plates containing Vogel's medium covered with a cellophane membrane. Cultures were allowed to grow for further 36 h under these conditions, and then photoinduced as described in [54] by exposure to white light for 5 min (fluence rate 27 mmol m^-2 s^-1). For injury-induced conidiation, fungal colonies were grown in total darkness on PDA at 27°C for 72 h, cut in stripes with a scalpel and incubated for an additional 24 h in the dark at 27°C. For mere mycelial growth, T. atroviride was grown on PDA plates in complete darkness at 28°C for 48 hrs.

For RT-PCR experiments strains were grown on PDA, covered with cellophane, in constant light at 25°C and harvested when the mycelia were ca. 5 mm apart (before contact), at contact of the mycelia and after T. atroviride had overgrown the host fungus by ca. 5 mm (after contact). As control T. atroviride was grown alone on plates and, in analogy to the growth conditions for the EST library, the peripheral hyphal zone was harvested.

Nucleic acid isolation and manipulation

Mycelia were ground to a fine powder under liquid nitrogen and total RNA was isolated using the RNeasy kit (Qiagen). cDNA was synthesized using poly A+ RNA and oligo dT primers with 5'-Biotin-GGCGGCCGCACAACTTTGTACAAGAAAGTTGGGT-(T)19-3'. cDNAs were ligated to attB1 adapter (5'-TCGTCGGGGACAACTTTGTACAAAAAAGTTGG-3'), size fractionated using sephacryl S-500 HR columns, and recombined to pDONR222-Lib vector.

Mycelium for the control condition of the RT-PCR was directly extracted from agar. In this case the guanidinium thiocyanate buffer was substituted by a buffer containing 0.6 M NaCl, 10 mM EDTA, 4% SDS, 100 mM Tris, pH8) and after phenol extraction a precipitation step with 8 M LiCl was introduced. The resulting pellet was dissolved in 0.3 ml RNase-free water and subsequently precipitated again with isopropanol. All other steps were carried out identically to the guanidinium thiocyanate method [55].

RT-PCR

For cDNA synthesis, RNA obtained from various cultivations was treated with DNase I (Fermentas, Burlington, Canada) and purified with the RNeasy MinElute Cleanup Kit (Qiagen, Valencia, CA, USA). 5 μg RNA/reaction were reverse transcribed using the the SuperScript™ III Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) and a mixture of random hexamer primer and oligo(dT)₁₈ primer. For transcript analysis cDNA equivalent to 50 ng RNA/reaction, 25 cycles and the primers and annealing temperatures listed in Additional File 7 were used. PCR reactions were carried out in a total volume of 50 μl using the GoTaq system (Promega, Madison, WI, USA) according to the manufacturer's instructions. 40 μl of each PCR reaction were separated on a 1.5% agarose gel containing 0.5 μg/ml ethidium bromide. As negative controls, to ensure the absence of genomic DNA, templates for PCR were prepared by as described above, but no reverse transcriptase was added during the cDNA synthesis step.

cDNA library sequencing

Bacterial colonies containing each T. atroviride cDNA library were plated onto agarose plates (254 mm plates from Teknova, Hollister, CA) at a density of approximately 1000 colonies per plate. Plates were grown at 37 C for 18 hours then individual colonies were picked and each used to inoculate a well containing LB media with appropriate antibiotic in a 384 well plate (Nunc, Rochester, NY). Clones were grown in 384 well plates at 37 C for 18 hours. Contained plasmid DNA for sequencing was produced by rolling circle amplification ([56] Templiphi, GE Healthcare, Piscataway, NJ). Inserts were sequenced from both ends using primers complimentary to the flanking vector sequence and Big Dye terminator chemistry then run on ABI 3730 instruments (ABI, Foster City, CA). The sequencing primers used were the following: for the mycelial injury library (vector: pDONR222; Fw: 5'-GTAAAACGACGGCCAGT, Rv: 5'-AGGAAACAGCTATGACCAT), for the light exposed library and the mycoparasitism library (vector: pEXP-AD502; Fw: 5' CTATTCGATGATGAAGATACC Rv: 5' AGAAGTCCAAAGCTCCACC), and for the dark exposed library (vector: pSPORT1; 5'-GTTTTCCCAGTCACGACGTTGTA, Rv: 5'-AGGAAACAGCTATGACCAT).

EST sequence processing and assembly

ESTs were processed through the JGI EST pipeline (ESTs were generated in pairs, a 5' and 3' end read from each cDNA clone). To trim vector and adaptor sequences, common sequence patterns at the ends of ESTs (Expressed Sequence Tags) were identified and removed using an internally developed software tool. Insertless clones were identified if either of the following criteria were met: >200 bases of vector sequence at the 5' end or less than 100 bases of non-vector sequence remained. ESTs were then trimmed for quality using a sliding window trimmer (window = 11 bases). Once the average quality score in the window was below the threshold (Q15) the EST was split and the longest remaining sequence segment was retained as the trimmed EST. EST sequences with less than 100 bases of high quality sequence were removed. ESTs were evaluated for the presence of polyA or polyT tails (which if present were removed) and the EST re-evaluated for length, removing ESTs with less than 100 bases remaining. ESTs consisting of more than 50% low complexity sequence were also removed from the final set of "good ESTs". In the case of re-sequencing the same EST, the longest high quality EST was retained. Sister ESTs (end pair reads) were categorized as follows: if an EST lacked an insert or was a contaminant then by default the second sister was categorized as the same. However, each sister EST was treated separately for complexity and quality scores. Finally, EST sequences were compared to the GenBank/EMBL/DDBJ nucleotide database in order to identify contaminants; non-desirable ESTs such as those matching non-cellular and rRNA sequences were removed.

For clustering, ESTs were evaluated with MALIGN, a k-mer based alignment tool (Chapman, Unpublished), which clusters ESTs based on sequence overlap (k-mer = 16, seed length requirement = 32, alignment ID ≥ 98%). Clusters of ESTs were further merged based on sister ESTs using double linkage. Double linkage requires that 2 or more matching sister ESTs exist in both clusters to be merged. EST clusters were then each assembled using CAP3 [57] to form consensus sequences. Clusters may have more than one consensus sequence for various reasons to include; the clone has a long insert, clones are splice variants or consensus sequences are erroneously not assembled. Cluster singlets are clusters of one EST, whereas CAP3 singlets are single ESTs which had joined a cluster but during cluster assembly were isolated into a separate singlet consensus sequence. ESTs from each separate cDNA library were clustered and assembled separately and subsequently the entire set of ESTs for all cDNA libraries were clustered and assembled together.

A list of the ESTs and the corresponding accession numbers is given in Additional File 8.

Functional analysis of the EST libraries

Since an annotated genome sequence of T. atroviride IMI206040 was available for this study [12], the following approach was taken for analysis of the 27658 ESTs: they were first used to construct a custom-BLAST database. This database was then queried (BLASTX) with the protein sequences extracted as FASTA files from individual KOG groups available at the T. atroviride genome portal [12]. To look for genes with special functions (e.g. CAZYs, receptors, proteins involved in sporulation etc.), these were manually retrieved from the database, combined to FASTA files and used for the query as above. All proteins contained in the same KOG group were manually verified. In the case of positive matches, the full length amino acid sequences of the respective proteins were further cross-checked by NCBI-BLAST or HMMER [58]. Significant hits were defined as those with 100% similarity between the amino acid sequences. Sequences exhibiting only a few aa mismatches (< 5%), which were speculated to represent sequencing errors, were manually verified by subjecting their ntsequences to TBLASTN against the T. atroviride genome sequence, and considered positive if no other protein showed the same similarity. The same strategy was used for ESTs, whose deduced aa-sequence was < 40 aa (and thus eventually representing a domain shared by several proteins), and identity was assumed when no other protein contained the same sequence. A complete list of identified genes is available in Additional File 9.

Proteins were considered to be orthologues by best bidirectional hits in BLASTP with a cutoff value of 1e-80. In cases of lower cutoff value, phylogenetic trees were constructed by NJ (MEGA 3.1; [59]), and an orthologue confirmed if the best bidirectional hit occurred in the same, strongly supported (bootstrap support > 75%) terminal clade.

Metabolic network construction

We used the EC numbers and the updated KEGG reaction database [60] to build a bipartite metabolic network, which was constructed based on the connection matrix of reactions [61]. Compared to reaction graph or metabolite graph, wherein either reactions or metabolites (called "node") are shown in an interconnected way, the bipartite network is more understandable because, similar to the biochemistry textbook, both the reactions and metabolites are visualized at mean time. The transcriptional data can be easily mapped to the bipartite network by using the visualization attributes of the reaction nodes, representing the expression levels of corresponding genes. Some currency metabolites [60], such as ATP, NADH etc., easily involve in hundreds of biological reactions. In order to reduce the complexity of the network to make the real biological changes more visible, we snipped some connections of reactions mediated by the currency metabolites. The software Cytoscape [62, 63] was used as a layout tool for the metabolic network. We mapped the ESTs to the whole network and its sub networks, by using the opacity of the background of reaction nodes to represent the level of expression which was normalized by the number of ESTs of the corresponding gene expressed in one condition divided by the total number of ESTs of this gene in all compared conditions. A minimal opacity of 50 is set to distinguish them with those not expressed. A series of sub networks, enabling a clearer view on the expression differences between mycoparasitism and mycelial growth, was extracted from the whole network according to the KEGG pathway category of the reactions.