Genome sequencing and transcriptome analysis of Trichoderma reesei QM9978 strain reveals a distal chromosome translocation to be responsible for loss of vib1 expression and loss of cellulase induction

Genome analysis of QM9978

Trichoderma reesei QM9123 was, until now, commonly believed to be the parental strain of QM9978. Since the genome sequence of QM9123 was published in 2010, we performed a comparison to QM9978 and showed that none of the 20 mutations identified in QM9123 is shared by QM9978 [16]. We therefore exclude the possibility that QM9123 is the parental strain of QM9978 and present an up-dated version of the T. reesei strain lineage in Fig. 1. No detailed record on the generation of the strain QM9978 and the closely related strain QM9979, which is also cellulase negative, exists, but we assume that the high number of mutations is a result of two or more subsequent mutagenesis events. Sequence comparison of QM9978 to QM9979 showed common mutations between those two strains indicating that they originate from the same parental strain (F.B., personal communication). Compared to high cellulase producer strains like T. reesei NG14 (one InDel, 78 SNVs, and nine large structural variations) and Rut-C30 (one deletion, two InDels, and 106 SNVs additional to NG14), the number of mutations is comparably low, probably because of alternative mutagenesis methods, and number and severity of treatments [14, 16, 17, 22].

From the 19 genes where a mutation occurred in the coding region, only the role of the blue-light receptor BLR1 in cellulase regulation has been studied in detail. Although BLR1 and BLR2 are the key regulators of the blue-light response in fungi and conidiation, it was shown that blr1 is not essential for cellulase expression [23] (Table 1).

Our research focused on the identified deletions and those mutations were subsequently verified by PCR analysis using primers flanking the deleted fragment. For two predicted deletions (on scaffold 1 and scaffold 16), we could not confirm the results obtained by Illumina sequencing as no PCR product was obtained. The recent genome assembly of T. reesei using chromosomal contact data in a computational approach identified scaffolds 1 and 16 to be located on chromosomes number five and seven, respectively [24–26]. Upon close examination of our sequencing data using the Integrative genome viewer program (http://software.broadinstitute.org/software/igv/) combined with the information on the Joint Genome Institute homepage (http://genome.jgi.doe.gov), it appeared that what we originally believed to be a deletion event could in reality be a translocation resulting from a chromosome break 1095 bp upstream of the ATG start codon of Trire2:54675 located on chromosome VII (Fig. 3a). To verify this hypothesis, we designed oligonucleotides flanking the assumed DNA break. Only PCRs with the primer combination 80028_MutaC_R and 54675_MutaC_R binding on the two different chromosomes resulted in an amplicon (Fig. 3b). Sequencing of this amplicon identified the region where the DNA breaks occurred. As a result from the translocation, the gene Trire2:37262 (a putative ubiquitin-protein ligase) upstream from Trire2:54675 on chromosome VII is replaced by Trire2:80028, a putative ABC transporter on chromosome V (Fig. 3b).

Gene expression profiling in T. reesei QM9978

Changes in gene expression between QM9978 and QM6a were analyzed in order to define the impact of mutations at the transcriptional level. Strains were grown on d-glucose for approximately 24 h until the carbon source was depleted from the medium (less than 3 g/L glucose remaining). Then, a lactose pulse of 10 g/L was added to the medium as lactose is known to be a strong inducer of cellulases [27]. Total RNA of two biological replicates was extracted at three different time points: 0 h (after carbon source depletion and before lactose addition), 24 and 48 h after the lactose pulse. Statistical treatments and differential analyses were performed using DESeq 1.8.3. A threshold of two for the log2 ratios with an adjusted p value <0.001 was used to identify genes differentially expressed (DE). Gene functions were identified by the aid of a manually annotated T. reesei genome database proprietary to IFPEN.

Comparison of the gene expression in QM9978 to QM6a identified a total of 785 genes differentially expressed (DE) in at least one time point. The number of down-regulated genes slightly increases during the time course, whereas the number of up-regulated genes remains stable. Interestingly, DE genes are mainly down-regulated, meaning their expression is lower in QM9978 then in QM6a. We compared the gene expression data with our sequencing results and identified five DE genes, which were affected by mutations. The genes were characterized by homology search using the NCBI BLASTp tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and the Joint Genome Institute genomics resource database (http://genome.jgi.doe.gov). The annotation was additionally verified manually using the FUNGIpath database (http://embg.igmors.u-psud.fr/fungipath/). This way we identified ID51893—an orthologue of the lipid body protein ppoC (A. nidulans) involved in coordinating sexual and asexual sporulation [28], ID54675—a homologue of the transcription factor vib-1 involved in vegetative incompatibility and cellulose utilization in N. crassa [29], ID68803—a predicted RTA1-domain protein, and ID108418 and ID122376—two hypothetical proteins of unknown function.

From the total of 785 genes, only 153 genes (23%) were DE at all time points. In this gene set, we find 14 glycoside hydrolases (including cbh2, xyn2, and agl1), eight putative transporters (including crt1), and four transcription factors including vib1 (ID54675) and azf1/sah-1 (ID61055) as identified by homology search as described above [30]. The N. crassa homologues of these two transcription factors have been shown to be involved in glucose signaling and to be up-regulated during growth on the silvergrass Miscanthus [31]. Of these, 153 genes two are impacted by mutation: vib1 and the putative RTA1-domain protein (ID68803).

Among the four genes surrounding the translocation on chromosomes V and VII identified during our sequence analysis, only vib1 and Trire2:80028, a putative ABC transporter whose expression is diminished in QM9978, show altered expression levels (Fig. 3).

Blastp analysis showed that T. reesei vib1 shows 49% similarity to N. crassa vib-1 [29, 32]. Since N. crassa vib-1 has been shown to be essential for cellulase production in N. crassa, we investigated whether a link could be established between the translocation, vib1 expression, and the loss of cellulase production in the QM9978 strain.

Analysis of the role of vib1 in cellulase induction

Constitutive expression of vib1 restores cellulose utilization in QM9978

To determine whether vib1 is involved in cellulose deconstruction by T. reesei, we constructed a vib1 overexpression QM9978 strain (QM9978 +vib1) by ectopic integration of the gene under the control of the strong constitutive promoter of T. reesei glyceraldehyde-3-phosphate dehydrogenase (gpd) [33]. For two independent clones, we monitored enzyme production during 4 days using a fast screening assay based on azurine-crosslinked cellulose (AZCL-HE-cellulose) that releases a soluble blue dye into the agar medium after enzymatic hydrolysis by endoglucanases, in which the size of the blue halo on an agar plate corresponds to the amount of secreted cellulases. Plates supplemented with 10 g/L glucose or lactose (for cellulase induction) and 1 g/L AZCL-HE-cellulose were used and the results of one clone per strain are representatively shown in Fig. 4 and Additional file 1: Figure S1. Strain QM9978 failed to secrete a sufficient amount of cellulases on lactose to decompose the AZCL-HE polymer, whereas the agar of the QM9978 +vib1 strain was stained blue due to enzymatic hydrolysis. To verify that no residual VIB1 activity remains in QM9978, we constructed the vib1 deletion strain QM9978 −vib1. Like QM9978, this strain was unable to decompose AZCL-HE-cellulose, indicating that the mutation in QM9978 equals a null phenotype (Fig. 3).

Since the data from the plate assay suggested a strong impact of vib1 on cellulase production, we produced liquid cultures grown on 10 g/L lactose/Solka-floc (cellulase induction) and on 10 g/L glucose (repression) to determine the protein content in the supernatant and the volumetric enzyme activity (see “Methods” for details on the protocol). Both the extracellular protein concentration and the filter paper activity were found to be significantly increased in strain QM9978 +vib1 compared to QM9978 (p < 0.0001 and p < 0.01, respectively) (Fig. 5).

The translocation in QM9978 occurred 1095 bp upstream of the vib1 ATG start codon. Therefore, the indicated promoter of approximately 1 kb remains unaffected [34, 35]. To test whether the abolished vib1 expression in QM9978 is due to the missing 5′ upstream region, we constructed a QM9978 strain expressing vib1 under the control of 850 bp of its native promoter, ectopically integrated. This allowed us to investigate the role of the 5′ region absent in QM9978, for example, if activator-binding sites were affected. The resulting strain QM9978_vib1np had an enzyme secretion capacity comparable to QM9978 +vib1 on plate test using AZCL-HE-cellulose (Fig. 4). These results were confirmed by real-time qPCR analysis, where we detected a significant accumulation of cbh1 (cel7a) transcript, both for the overexpression strain QM9978 +vib1 and for QM9978_vib1np (Fig. 6a). Expression levels in both strains remained lower compared to strain QM9414, probably due to the fact that strain QM9414 has undergone a different mutagenesis treatment during strain engineering leading to higher expression of chb1 than a cellulase-restored QM9978 strain.

Analysis of vib1 transcription in QM9978 +vib1 and QM9978_vib1np confirmed that vib1 transcripts were present at levels equivalent to QM9414, a strain with unimpaired vib1 expression (Fig. 6b).

Fungal strains, culture conditions, and measurement of cellulase formation

Trichoderma reesei wild-type strain QM6a (ATCC 13631), the early cellulase producing mutant T. reesei QM9414 (ATCC 26921), the cellulase-negative mutant QM9978 (obtained originally from Dr. K.O’ Donnell, U.S. Department of Agriculture, Peoria, IL), and the hypercellulosic strain Rut-C30 (ATCC 56765) were used throughout this study and kept on potato dextrose agar (Sigma, St. Louis, MO) at 30 °C (Additional file 2: Table S1). T. reesei was grown in liquid culture in 250-mL Erlenmeyer flasks on a rotary shaker (250 rpm) at 30 °C in 50 mL minimal medium as described [60] using frozen spores as inoculum. Carbon sources used are specified in the legends to the respective figures, and were used as 1% (w/v) in batch cultures. To assess the cellulase formation on different carbon sources using a colorimetric assay based on AZCL-HE-cellulose, strains were grown on minimal medium plates supplemented with 5 g/L of the respective carbon source in addition to AZCL-HE-cellulose (Megazyme International, Ireland) to a concentration of 1 g/L. Determination of protein content in culture supernatants was assessed by the modified Lowry protein assay as described using three biological and at least two technical replicates [61]. Global cellulase activity was measured using the IUPAC standard Filter Paper Assay [62], after a 15-fold miniaturization similar to [63] which allowed working in 2-mL Eppendorf tubes. In order to surround the desired 4% conversion yield, four enzyme dilutions were tested for each sample. This miniaturized protocol was validated by comparison with the standard protocol. Statistical analysis of the results was performed with the Prism GraphPad software (https://www.graphpad.com/scientific-software/prism/).

For transcriptomic studies 50 mL of minimal medium (5 g/L (NH₄)₂SO₄, 5 g/L KH₂PO₄, 11.7 g/L trisodium citrate, 0.6 g/L MgSO₄, 0.6 g CaCl₂, 20 g/L glucose) were prepared in a 250-mL Erlenmeyer flask. The medium was supplemented with the following trace elements: 0.5 mL/L FeSO4·7H₂O, 2 mL/L CoCl₂, 0.16 mL/L MnSO₄, and 0.14 mL/L ZnSO₄. Prior to sterilization, the pH of the medium was adjusted to 5.4. Strains QM6a and QM9978 were inoculated in duplicate with approx 10⁶ conidia/mL. Cultivation was carried out at 30 °C at 150 rpm. When the glucose was close to depletion (<3 g/L), 5 mL of a lactose solution was added (25 g/L lactose, 6 g/L (NH₄)₂SO₄, 5 g/L trisodium citrate). Glucose concentration was assessed by enzymatic reaction using an analox glucose analyzer GM10 (Imlab). Extracellular protein concentration was measured against a BSA standard (0–1.5 g/L range with second-order regression) using the Quick Start Bradford Protein Assay kit (Bio-Rad).

To monitor gene expression levels via real-time qPCR, the replacement technique described by Sternberg and Mandels [64] was used: mycelia, grown for 24 h on 1% (w/v) glycerol as the carbon source under otherwise similar conditions as described above, were cautiously filtered through Miracloth (Calbiochem, EMD Biosciences, La Jolla, CA), washed with tap water, and resuspended in minimal medium supplemented with the respective carbon source at 1% (w/v) and the incubation continued for 5, 16, and 24 h. To guarantee equal treatment of the cultures, they were done in parallel and results of three technical replicates were obtained.

Escherichia coli NEB 10-beta competent cells (New England Biolabs, France) were used for the propagation of vector molecules and DNA manipulations.

Vector construction for generation of gene deletion and overexpression mutants

Deletion cassettes consisting of 1–1.5 kb fragments of the gene specific flanking regions interrupted by the hygromycin resistance cassette were assembled by yeast recombinational cloning [50, 65]. Oligonucleotides (10 mM) 5F + 3R and 3F + 3R were used for amplification of the individual flanking regions from genomic T. reesei DNA using Phusion polymerase (Thermo Fisher Scientific). Approximately 19 bp that overlap with the pRS426 (URA +) yeast shuttle vector or the hph gene were introduced by PCR at each flanking end to allow homologous recombination. The hph marker cassette was amplified from pLHhph1 with oligonucleotides hphF and hphR [65]. A PCR touchdown program ranging from 62 to 58 °C for annealing was used for amplification. Oligonucleotide sequences are shown in Additional file 3: Table S2. Yeast transformation was performed as described [66] using the lithium chloride/polyethylene glycol procedure. Yeast strain WW-YH10 (ATCC 208405) was transformed with both flanking regions, the hph marker cassette and an EcoRI/XbaI-digested plasmid pRS426 [67]. Transformants were selected on synthetic complete dropout medium (SC-URA with uracil dropped out). Following total DNA isolation from liquid SC-URA medium, plasmids were introduced into chemically competent NEB 10-beta E. coli cells (New England Biolabs) [68]. The outside primers 5F + 3R were used to synthesize the complete deletion cassette from plasmids isolated from S. cerevisiae.

For construction of overexpressing strains, the coding region of vib1 was amplified using oligonucleotides VIB1 for and VIB1 rev and placed under the control of the strong constitutive gpd promoter (oligonucleotides promGPD for and PromGPD rev) followed by the gpd terminator (termGPD for and termGPD rev) and the hph resistance cassette (hph for and hph rev, see Additional file 3: Table S2) all inserted into the XbaI/EcoRI-linearized plasmid pUC19 yielding plasmid pUC_vib1_OE. Oligonucleotides K7 VIB1 for and K7 VIB1 rev were used to amplify the complete overexpression cassette from pUC_vib1_OE.

For reintroduction of vib1 under the control of the native promoter into the strain QM9978, the coding region plus approximately 850 bp of the promoter region (after the translocation breakpoint) together with 500 bp of the terminator region were amplified using oligonucleotides VIB1_F and VIB1_R (Additional file 3: Table S2) and introduced in the SacI/BamHI-linearized plasmid pLHhph1 using the Gibson Assembly Cloning Kit (New England Biolabs) [65]. Five micrograms of the circular plasmid were used for protoplast transformation of T. reesei.

Transformation of T. reesei and analysis of transformants

Trichoderma reesei QM9414, QM9978, and Rut-C30 were used for the generation of vib1 deletion and overexpression strains. Protoplast preparation and transformation were performed as previously described [69]. The deletion cassettes were purified after PCR (QIAquick PCR Purification kit, QIAGEN) and concentrations were determined (Nanodrop Spectrophotometer, Peqlab). After transformation protoplasts were stabilized and regenerated on selection medium supplemented with 100 mg/mL hygromycin B. For sporulation, transformants were transferred to 12-well plates and purified by plating conidiospores onto plates with 0.1% Triton X-100 as colony restrictor. Single colonies were transferred to selective medium and screened for correct integration of the deletion or overexpression cassettes by PCR analysis (primer abbreviated with “ch” in Additional file 3: Table S2). Insertion of the gpd-vib1 overexpression cassette was verified with primer pairs TR46-verif 5′vib1 for + TR47-verif 5′vib1 rev, TR48-verif 3′ vib1 for + TR49-verif 3′vib1 rev, and TR50-verif hph vib1 for + TR51-verif hph vib1 rev amplifying the 5′or 3′ region and the hph marker gene, respectively. Three or four independent clones were isolated for each strain construction and analyzed as biological replicates.

Gene expression analysis by real-time qPCR

DNase treated (DNase I, RNase free; Fermentas) RNA (5 mg) was reverse transcribed with the RevertAidTM First Strand cDNA Kit (Fermentas, Germany) according to the manufacturer’s protocol using a combination (1:1) of the provided oligo-dT and random hexamer primers. All assays were carried out in 96-well plates, which were covered with optical tape, as previously described [70]. Primers, amplification efficiency, and R-square values are given in Additional file 4: Table S3. The amplification efficiency was determined using triplicate reactions from a dilution series of cDNA, and then calculated from the given slopes in the IQ5 Optical system Software v2.0. Expression ratios were calculated using the REST 2009 Software (Quiagen, Germany). All samples were analyzed in two independent experiments with three replicates in each run.

RNA sample preparation

At each time point, 5 mL culture medium containing mycelia was filtrated on 1.2-µL GF/C glass microfiber membranes and flash frozen in liquid nitrogen. Frozen mycelia were disrupted and homogenized in Lysing Matrix C (Fast RNA Pro Red Kit, MP Bio) in the presence of 700 µL Lysis buffer RLT (RNeasy Mini Kit, Qiagen) complemented with 7 µL of β-mercaptoethanol on a rotor–stator homogenizer (FastPrep, MP Bio, Santa Ana, CA, USA). After centrifugation, the lysate was transferred on a Qiashredder spin column (Qiagen) and centrifuged for 2 min to eliminate all cell debris. The flow-through was then mixed with 0.5 volumes of pure ethanol and transferred to an RNeasy spin column. All further steps, including an on-column DNAse digestion, were realized following the manufacturer’s instructions. RNA quantification was performed with Qubit 2.0 (Life Technologies, Carlsbad, CA, USA). The quantity and quality of the total RNA was determined using a NanoDrop ND-1000 spectrophotometer (Nanodrop Technologies, Wilmington, USA) and by electrophoresis on 1% agarose gel. An additional quality control was performed on the Bioanalyzer 2100 system (Agilent, Santa Clara, USA).

RNA-seq library preparation and analysis

Library preparation and Illumina sequencing were performed at the Ecole normale superieure genomic platform (Paris, France). Messenger (polyA+) RNAs were purified from 1 µg of total RNA using oligo(dT) primer. Libraries were prepared using the strand non-specific RNA-Seq library preparation TruSeq RNA Sample Prep Kits v1 (Illumina). Libraries were multiplexed by four on one single flowcell lane and subjected to 50-bp paired-end read sequencing on a HiSeq 1500 device. A mean of 52 ± 10 million passing Illumina quality filter reads was obtained for each of all the samples. We performed RNA-seq analysis using the Eoulsan pipeline [71]. Only read 1 was considered for the further analyses. Before mapping, poly N read tails were trimmed, reads ≤11 bases were removed, and reads with quality mean ≤12 were discarded. Reads were then aligned against the T. reesei genome (version 2 from the JGI website) using the Bowtie mapper (version 0.12.7) [72]. Alignments from reads matching more than once on the reference genome were removed using Java version of samtools [73]. To compute gene expression, the T. reesei genome annotation from JGI (version 2) was used. All overlapping regions between alignments and referenced exons were counted. The sample counts were normalized using DESeq 1.8.3 [74]. Statistical treatments and differential analyses were also performed using DESeq 1.8.3. A threshold of two for the log2 ratios with an adjusted p value <0.001 was used to identify genes differentially expressed (DE). Gene functions were identified by the aid of a manually annotated T. reesei genome database proprietary to IFPEN.

The RNA-Seq gene expression data and raw fastq files are available at the GEO repository (http://www.ncbi.nlm.nih.gov.gate1.inist.fr/geo/) under Accession Number GSE89199.

Illumina genome sequencing of T. reesei QM9978

Chromosomal DNA from T. reesei QM9978 was prepared as previously described [14]. Fragment libraries were prepared according to the TruSeq^® DNA Sample Preparation Guide from Illumina (http://www.illumina.com). These libraries were then loaded onto the cluster generation station for single-molecule bridge amplification using the Standard Cluster Generation kit from Illumina. The slide with amplified clusters was then subjected to sequencing on the Illumina Genome Analyser I (GAI) for single reads using the 36 cycle Sequencing Kit version 1 from Illumina.

The whole genome sequence of T. reesei QM9978 will be deposited at the Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra).

Sequence alignment and analyses

Illumina short reads from QM9978 were mapped onto the T. reesei genome (http://genome.jgi-psf.org/Trire2/), using the Maq 0.6.6 software solution [75]. Mapping was done with two maximum mismatches. InDels and SNVs were identified using Maq 0.6.6. Homozygous mutations were selected and filtered on genomic context (complexity  =  1, uniqueness  > 15.8, GC percentage between 0.31 and 0.74). Each mutation was manually checked using the Integrative Genome Viewer and compared to the sequence of the Rut lineage to remove the false mutations coming from initial QM6a error sequencing. Large structural variations were searched using the BreakDancer software and filtered on the reads number covering the genomic variations [76].

We evaluated the location of SNVs and deletions according to gene annotations using the “filtered models” from the JGI website. From this annotation, we calculated the position of intron, promoter (using an 800 bp upstream region), and terminator (using a 200-bp downstream region).

Functions of genes were predicted by similarity to orthologous genes when found in other fungal taxa using the FUNGIpath database (http://embg.igmors.u-psud.fr/fungipath/). If no orthologous gene function could be identified but a gene with a homology e-value superior to e-10 was identified, the product of this gene was annotated as “similar to.” Conserved domains were identified either by interproscan (http://www.ebi.ac.uk/interpro/search/sequence-search) or by NCBI BLASTp search (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

To validate the translocation between chromosomes V and VII, we amplified the region around genes encoded by the protein IDs 80028 and 54675 by PCR using oligonucleotides up and downstream of the coding region (Primer 54675_MutaC_F + 54675_MutaC_R and 80028_MutaC_F + 80028_MutaC_R Fig. 3b; Additional file 3: Table S2). PCR amplicons were sequenced with adequate primers using Sanger cycle sequencing/capillary electrophoresis (Microsynth AG, Austria).