Genome wide analysis of protein production load in Trichoderma reesei

The filamentous fungus Trichoderma reesei (teleomorph Hypocrea jecorina) is a widely used industrial host organism for protein production. In industrial cultivations, it can produce over 100 g/l concentrations of extracellular protein [1]. Interestingly, it exhibits a low growth rate protein production phenotype [2], i.e. the highest rates of cellulase and hemicellulase secretions are observed at low growth rate [3–6].

In industrial environments such conditions can exist for a prolonged period during the feeding stage of a fed-batch cultivation. However, many industrial processes are batch cultivations with constantly changing growth rates. In both conditions T. reesei efficiently produces cellulases and also heterologous proteins under cellulase promoters, given a cellulase expression inducing carbon source. Cellulase expression is well induced by sophorose (for review see [7, 8]), but also by cellobiose [9], lactose [10] and cellobiono-1,5-lactone [11]. Recent transcriptome analyses of cellulase producing conditions, have revealed various responses to protein production: partial induction of the secretion system, induction and carbon source specific modulation of cellulases, but also complex metabolic responses whose importance and regulation remains to be discovered [6, 12–18]. The major regulators of cellulases are transcription factors cre1 [19] and xyr1 [20] (for review see [8]), but also ace1 [21], ace2 [22], ace3 [15], clr2 [23] and other Zn2Cys6 zinc cluster transcription factors [15], GCN5-like histone acetyl transferases [6] and histone methyltransferases [24] have been implicated.

The level of endoplasmic reticulum (ER) stress in protein producing conditions is typically unclear in published T. reesei transcriptome analysis. Induction of the unfolded protein response (UPR), the conserved eukaryotic stress response to high ER load, can be assessed by non-conventional splicing of the hac1-intron [25] and hence is not detectable from microarray analysis. Major factors like foldase bip1 [26] and pdi1 [27] are typically induced. However, in Saccharomyces cerevisiae, the expression of BIP1 and PDI1 has been shown to be a poor indicator of UPR [28].

With the secretion pathway factors (for review see [29]) and cellulase transcription factors extensively studied, new approaches for improving protein production in T. reesei are needed. Genome scale stoichiometric metabolic modelling has been very successfully applied to improvement of production of a single metabolite (for review see [30]). In this field, the stoichiometric matrix, i.e. a mathematical description of most biochemical reactions of the cell, is analyzed to select enzyme genes to be removed, regulated or added, or to optimize cultivation strategies in order to increase metabolite production. In contrast, successful application of stoichiometric modelling to improve protein production has been rare, although various other strategies for engineering metabolite supply for protein production have been used (for review see [31]). A notable exception is analysis of super oxide dismutase production with stoichiometric modelling in Komagataella phaffii (Pichia pastoris) [32]. Filamentous fungi, such as T. reesei, can reach far higher production levels and are also in many other respect distinct organism from yeasts. Stoichiometric modelling of T. reesei has not been, to our knowledge, reported previously. Also, successful stoichiometric modelling typically requires growth conditions rarely used in published T. reesei work i.e bioreactor cultivations, extensive culture monitoring and a defined carbon source.

In this paper we carry out a well controlled bioreactor experiment to study the effect of variation of protein production load to the physiology of T. reesei. Our RNA-sequencing transcriptomics detects the concentration of carbon source as the most important determinant of transcriptome. We extensively analyze the growth medium during cultivation and find novel metabolic responses such as production of glycerol and a cellotriose-like compound. We then use this cultivation data for flux balance analysis of the metabolism of T. reesei and demonstrate for the first time the use of genome wide stoichiometric metabolic modelling for T. reesei. We show that our model can predict protein production rate and provides novel insight into the metabolism of protein production. We also provide this unprecedented cultivation and transcriptomics data set for future modelling efforts

In order to dissect the effect of variable protein production load to the secretory machinery and physiology of T. reesei in batch cultivations, sufficient variation in protein production is required. In our data set, a 1.7 times higher protein amount was produced by the highest producing strain LipPr4d, i.e. the lipase producer, as compared to the lowest producing strain CutCBHd. Together with the above presented results this clearly shows that sufficient variation was reached. The range of protein amount produced is typical for T. reesei cultures on minimal medium [58] and the specific protein production rate reaches levels reported for industrial protocols [16, 59]. As a carbon source we selected cellobiose. Cellulose or other more complex lignocellulosic material would have been likely to lead to stronger induction of protein production. However, in order to be able to quantify the carbon source uptake rate to enable stoichiometric modelling, we needed a defined carbon source.

A number of variables, such as growth phase, growth rate or the specific product proteins itself, can confound the effects of the load of protein production on transcriptome and physiology. We try to overcome these problems firstly by including variation of these variables in the data set and secondly by data analysis methods. All of our strains have different specific growth rates and produce a different protein mixture, some with heterologous product proteins and some with different sets of endogenous proteins. Hence, we are able, to a great extent, to control for these factors and detect responses that are not merely results of such confounding factors. On the data analysis side we use a combination of correlation, differential expression and clustering methods to overcome biases in data and shortcomings of individual methods. For example gene expression cluster 5 (Additional file 1: Figure S4) shows a peak of expression at 40  h in all strains except for the lipase (LipPr4d and LipPr4dCt) strains. The cluster is enriched in genes which are expressed at a significantly lower level in the high producing strains (Cel4dCt, LipP4rd and LipPr4dCt) than in the low producing strains (Cel4d, CutCBHd and CutCBHdCt) at 40  h (Additional file 1: Figure S6). However, the responses of the lipase production strain (LipPr4dCt) and cellulase deletion control strain (Cel4dCt) (i.e. the strains producing protein best) at 40h are quite opposite even though their protein production rate (Fig. 1f) and amount of protein produced at this time point are almost identical (Fig. 1b). Hence, expression of the genes in cluster 5 is not effected by the protein production load but by some other factors. Without including sufficient variation in the experimental design this response could be easily misinterpreted.

In order to calculate meaningful uptake and secretion rates, i.e. first order derivatives, from the concentration data of batch cultivations, the measurement data has to be either interpolated, smoothed or modeled. Typically this is done for example with higher order polynomials, smoothing splines [60] or a Monod model [61]. We propose to use heteroscedastic Gaussian processes [38] for this task. Gaussian processes infer a distribution of all interpolating functions that match the observation, instead of a single interpolant. The distribution mean characterises the most likely interpolation function, which are shown in Figs. 1 and 2 (left). We employ time-dependent, i.e. heteroscedastic, observation noises that can realistically model the differing observation variances at different measurement times. Finally, we derive the derivative distribution as a second Gaussian process [62], whose mean indicates the most likely uptake and secretion rates over time, which are shown on the right side column of Figs. 1 and 2. At the first and last time point of the time series these rate estimates might slightly suffer from lack of preceding and subsequent data, respectively.

We detected a cellotriose-like compound, glucose, glycerol and ethanol from the growth media. For the moment we can only speculate about the source of these compounds, but clearly, based on their dynamics, they are produced by T. reesei cells or by secreted enzymes. To our knowledge this is the first time that such metabolism has been reported for T. reesei. The cellotriose-like compound and glucose are likely to be produced by extracellular enzymes from the carbon source cellobiose. Indeed, cellotriose production from cellulose has been described for an extracellular extract of Streptomyces by Elwyn T. Reese [63]. Furthermore, it has been proposed that cellotriose and cellotetraose act as cellulase inducers in Phanerochaete chrysosporium [64] and that production of transglycosylation products from cellobiose is important for cellulase induction of T. reesei [9]. The formation of the cellotriose-like compound we detect, could be a side activity of CAZymes expressed at low basal level even in non-inducing conditions [65], which after detection by T. reesei enables full induction of the cellulytic system as has been previously suggested (for review see [8]). The fact that most of this cellotriose-like compound is produced by the cellulase deletion strain Cel4d (Fig. 2) suggests that it is made by an auxiliary enzyme, whose production is increased by cellulase deletion. Transport system for such small molecule inducers have been described [66, 67] and accordingly 5 out of our 18 top extracellular protein production rate correlated genes in continuous and batch conditions are transporters (Table 5). In addition the putative lactose and/or cellobiose permease crt1 (3405) [13, 67] lies beside the homologue of Neurospora crassa clr2 (26163) [15, 23] found in cluster 24 and having a FIRE predicted Xyr1 element on its promoter region. The rlog2 expression values of crt1 and clr2 are correlated with a Pearson correlation of 0.94 in our experiment, hence they are likely to be commonly regulated also in general. Co-expression of clr2 and crt1 under cellulase inducing conditions have also been reported previously [15].

Glucose, like the cellotriose-like compound, can be expected to be a product of extracellular enzymes hydrolysing the carbon source cellobiose or the cellotriose-like product. However, why would it accumulate in conditions where carbon source starts to be consumed and there is ample living biomass? The answer could lie in the discrepancy between cellobiose consumption rate and growth rate (Figs. 1e, g, 3b). Although these two variables correlate strongly, the dynamic range of the growth rate is smaller and in particular the cellobiose uptake rate decelerates faster than the growth rate at 40–64 h. Hence, while the cells decelerate cellobiose uptake, the already secreted cellulotic system is not regulated and therefore residual glucose is temporarily accumulated.

Glycerol and ethanol are likely to be products of intracellular metabolism secreted to the growth medium. Glycerol is produced by Saccharomyces cerevisiae under anaerobic conditions, i.e. while it is producing ethanol, to compensate for cellular reactions that produce NADH and therefore to balance its redox state [68]. In addition, glycerol is produced as a response to a number of stresses in S. cerevisiae such as osmotic and oxidative stress [69, 70].

In our study, glycerol production was detected in all the cultures at the time when the biomass amount in the cultures was the highest (i.e. at 60  h), whereas ethanol was produced only in the cultures of the strain producing cutinase (CutCBHd) and its control strain (CutCBHdCt) at a time period preceding the glycerol production phase (i.e. from 16–24  h). It is possible that there would be transient oxygen limitation during the cultivation phase when biomass reaches the maximum level and starts to decline, i.e. at the same time when glycerol production was detected. However, due to the differences in the production time points, a direct link between glycerol and ethanol production in T. reesei is less plausible. At the time when glycerol production is detected, at 40–64  h, respectively, most protein production has already taken place. Possible stress conditions at that stage could be due to accumulation of vesicles and oxidative stress from protein folding [28, 71, 72].

The gene expression clusters whose expression correlates well with the rate of glycerol production are enriched in genes that encode proteins acting in the later steps of protein secretion after ER. These genes tend to be expressed at a higher level in the strains producing high amount of protein at 64  h as compared to the strains producing less protein (Fig. 5b). Hence, although there is no direct link between protein production and glycerol production, the 64h induction of these genes could be a response to accumulated protein secretion stress.

Although the ethanol production we detect is low and transient, ethanol production has been described for T. reesei also previously [73, 74]. We provide the first transcriptomics data set for such conditions which could be useful for development of a consolidated bioprocess of ethanol production.

We measured the concentrations of intracellular amino acids for time points 24, 40 and 64  h and found that the case-control strain pairs resemble each other, indicating that the strain history and other properties of the strains are more important explanatory factors for the intracellular amino acid content than the protein production load. Given these circumstances we did not try to use the intracellular amino acids for modelling, however most of the amino acids, especially vacuolar storage amino acids arginine and glutamine [75, 76], exhibit a clear downward trend over the cultivation time (Additional file 1: Figure S3). This apparent use of storage amino acids could explain why our flux balance analysis in many cases predicts lower protein production rates than are measured, in particular for early the time points at 16 and 24  h (Fig. 7). The cells could contain resources stored in the early phase of the cultivation which are not covered by our modelling, but which are in reality used in protein production.

It has been shown that protein production rates can be predicted for Komagataella phaffii (Pichia pastoris) based on online measurements of OUR and CER by a stoichiometric model [77]. Also in our data set such a relationship exists (Fig. 3f), but it breaks down in the samples where the protein production rate is the highest. In these sampling points the production rates are approximately 1000 times higher than in [77]. In general, our protein production rates are 10 to 100 times higher than those used in previous stoichiometric modelling of fungal protein production systems [32, 77–79]. Hence, the combination of T. reesei metabolic model and our cultivation data set provides a set up of modelling protein production in far more relevant conditions for industrial protein production. This could also indicate that at protein production rates typical for T. reesei yet undiscovered metabolic interactions exist.

It has been argued that the metabolic activity of cells, i.e. the flux through metabolic reactions, is controlled by the rate of carbon source uptake or by the rate of following key reactions [80–82]. As flux and transcript levels do not correlate in general in eukaryotes [83, 84] at least in the central carbon metabolism, what then controls transcript levels? In our experiment, based on the correlation of gene expression levels and cultivation variables (Table 2), the concentration of biomass and carbon source are the most important determinants of the transcriptome, effecting a similar set of genes. Although it has been shown that cell density has an effect on the trancriptome of T. reesei, the carbon source availability is a far more fundamental question for an organism and hence it is likely that the carbon source concentration is the main determinant of transcriptome. The fact that the carbon source uptake rate and carbon source concentration both exert control on different cellular levels highlights the importance of assessing both factors in order to understand the cell.

At the transcriptome level we find two major responses in the secretion machinery and three major responses in amino acid metabolism. Regarding the amino acid metabolism, the most prominent response is the general downward regulation of primary metabolism genes with declining carbon source concentration (mainly gene expression cluster 31, Fig. 4). Secondly, there is the induction of starvation related amino acid metabolism genes that have a negative correlation with the carbon source concentration (cluster 1). Thirdly, a special sub set of a amino acid metabolism genes is correlated with the protein production rate and differentially expressed between low and high producing strains (cluster 24, peak at 16  h and cluster 28 peak at 24  h) and putatively regulated by Xyr1. However, it is not obvious from our data why these individual genes are regulated.

The development of a stoichiometric genome wide metabolic model is an incremental project as exemplified by the progress from first published S. cerevisiae model [85] to the latest S. cerevisiae consensus model Yeast 7 [86]. Our T. reesei model was constructed without manual refinement though the CoReCo pipeline [55]. We have shown that the performance of models created using CoReCo is comparable to early S. cerevisiae models [55] in simulations. Our ability to predict protein production rates in this experiment with the T. reesei model shows its functionality and usefulness (Fig. 7). Overall the model predicted protein production rates of almost half the speed than were measured. Cellulases are heavily glycosylated proteins, a fact which our current model does not take this into account, but rather tries to produce all protein mass from amino acids. Amino acids are energetically very costly to produce for the cellular metabolism. Hence, the actual use of stored resources by the cells and omission of glycosylation could explain the prediction of lower than measured protein production rates.

The model predicts particularly high protein production rates at 40  h in comparison to the measured ones. As the transcriptome is mainly determined by the carbon source concentration (Table 3), which at 40  h has dropped 20–50  % from 16  h (Fig. 1c), this discrepancy is likely a result of transcriptional repression of main cellulases transcripts and/or the actual protein secretion machinery (Figs. 4, 6 clusters 24 and 28), yet again another factor not taken into account by our stoichiometric model.

The metabolism of cysteine and methionine (a precursor of cysteine) is not only highlighted by transcriptomics data, but also by our flux modelling (Fig. 9). In contrast to most metabolic reactions that correlate with the carbon source uptake rate as expected, the flux through the methionine salvage cycle exhibits a higher correlation with the protein production rate. In parallel, the gene expression for a cysteine synthase 76018 correlates with the protein production rate, but the flux through this reaction does not. Hence, it could be that cellulase protein production requires elevated cysteine and consequently methionine metabolisms and possibly flux to cysteine. Overall asparagine and cysteine are the two most over represented amino acids (43 and 37  % more, respectively) when comparing the relative amino acid contents of sequences of T. reesei proteins found to be secreted based on 2D-gel analysis [6] and other T. reesei proteins. More fundamentally, the capability to assimilate sulphur could be the limiting factor. A regulatory link between sulphur metabolism and cellulase expression has been shown earlier [87]. The induction of taurine catabolism genes at 24 and 40  h in high producing strains (LipPr4d, LipPr4dCt and Cel4dCt) also highlights the relevance of sulphur metabolism. In T. reesei’s natural environment taurine is a major source of sulphur [88]. Hence, the induction of these genes could be a natural, yet in industrial conditions futile, response to lack of sulphur.

The protein secretion machinery responds in our experiment by ER related genes that are correlated to protein production rate and are putatively controlled by Xyr1 (Fig. 4 cluster 24) and by genes that are involved in the later stages of the secretory pathway, after ER and which are correlated to the glycerol production rate (Fig. 4, cluster 22). Such a regulation of the secretion machinery and amino acid metabolism by Xyr1 has not been previously implicated based on expression data, although genome sequence analysis has proposed that Xyr1 would regulate these and numerous other functional categories [50]. As no hac1 splicing was detected, it seems that in these conditions mere Xyr1 controlled induction of the secretion machinery is sufficient to cope with the protein production load.

We found that the expression of main cellulases followed the protein production rate (Fig. 6) and that their lack was compensated by expression of other cellulases. The induction of other cellulases implies that the repression under secretion stress (RESS) response was alleviated in the Cel4d strain [26].

In order to assess our results in the light of previous experiments we compared our analysis to an analysis of a large set of T. reesei microarray data [51] and found that genes involved in the major responses discussed above (gene expression clusters 22, 28 and 31, Fig. 4) are also co-regulated in the microarray data set. Furthermore, gene expression cluster 24 contains four genes (clr2 26163, ace3 77513, 108381, 70351) previously implicated as regulators of cellulase gene expression [15]. The overexpression of ace3 77513 and 108381 has been shown to have an effect on cellulase gene expression [15, 89, 90]. We have previously shown that although the level of paralogy in T. reesei genome is very low in general [91, 92], the genes of central carbon metabolism enzymes nevertheless have paralogous gene pairs with opposite regulation in response to protein production. Out of seven pairs of paralogues reported in [6] we detected five pairs with opposite regulation in response to protein production (LPD1a 67699, LPD1b 77373; LSC1a 22910, LSC1b 2223; LSC2a 103451, LSC2b 80881; PDB1a 76744, PDB1b 122745; TKL1-2a 120635, TKL1-2b 2211; TPI1a 68606, TPI1b 121789;) in this work. We also compared the correlation of gene expression to protein production rate between this experiment and [6] and found no similarity between the two studies. This is not surprising given the completely different physiological conditions in a chemostat [6] and a batch cultivation experiment. In this experiment we find that the carbon source concentration is the main determinant of transcriptome. However, in a chemostat the carbon source concentration is zero as all available carbon is immediately taken up. Therefore, in these batch experiments the growth rate appears as only a very minor determinant of the transcriptome while in a chemostat it emerges as dominant [6]. Nevertheless, central carbon metabolism paralogues and known key secretion pathway factors like pdi1 (gene expression correlation to protein production rate is here 0.69 and 0.80 in [6] ), dpm2 (0.56 and 0.77, respectively) emerge from both experiments. Hence, genes correlating with the protein production rate in these completely different cultivation regimes could reveal novel key factors (Table 5).

It has been argued that the flux of carbon source uptake can actively control metabolic activities of the cell. In our experiment the carbon source concentration appears as the main controller of the transcriptome. These parallel control systems would allow the cell to integrate information at different time scales: i.e. short-term responses (such as changes in metabolism in less than seconds) and long-term responses (such as transcriptional changes over hours). Their parallel existence also resolves the surprising lack of correlation between flux and transcriptome in eukaryotes.

Various approaches have been taken to improve protein production by filamentous fungi and T. reesei specifically. Manipulation of transcription factors has proven to be an efficient method to increase cellulase protein expression in many cases [15, 19]. Our modelling proposes that stoichiometry would allow production of more cellulases than measured, but the transcriptome becomes repressed due to depletion of the carbon source. Thus, the adjustment of metabolic control and the transcriptional regulation offers an interesting target for modification when aiming at improved protein production, e.g. by alleviating such a repression.

As a completely novel target for improvement of cellulase production our data analysis and modelling proposes processes of sulphur assimilation and cysteine metabolism. As amino acid metabolism is essential and tightly controlled at many levels detailed modelling will be required to select the correct targets and their manipulations.