In silico engineering of Pseudomonas metabolism reveals new biomarkers for increased biosurfactant production

Reconstruction of the engineered constraint-based genome-scale metabolic model

To reconstruct a metabolically engineered model of P. putida for rhamnolipids production, following Wittgens et al. (2011) and Tiso et al. (2016), we introduced two pathways for rhamnolipids biosynthesis by collating the full list of reactions catalyzed by RhlA and RhlB to an existing genome-scale model of P. putida, iJP962 (Oberhardt et al., 2011). The RhlA and RhlB genes produce rhamnolipids by three sequential reactions (Fig. 1A). RhlA is involved in the synthesis of HAA (Déziel et al., 2003), and is loosely bound to the inner membrane (Rahim et al., 2001). The next reaction is catalyzed by the membrane-bound RhlB rhamnosyltransferase and uses dTDP-L-rhamnose and an HAA as precursors, yielding mono-rhamnolipids (Rahim et al., 2001). The RhlA and RhlB genes are clustered with rhlR and rhlI, which encode proteins involved in their transcriptional regulation through the quorum-sensing (QS) response and they are arranged as an operon. RhlI and LasI synthesize the QS autoinducer molecules butanoyl-homoserine-lactone (C4-HSL) and 3-oxo-dodecanoyl-homoserine-lactone (3-oxo-C12-HSL), respectively. When their concentration reaches a threshold, they bind to the regulator proteins and induce the expression of the Rhl-genes (Wittgens et al., 2017). QS response regulates the production of rhamnolipids (Van Delden & Iglewski, 1998), as well as hundreds of additional genes (Hentzer et al., 2003; Schuster et al., 2003; Wagner et al., 2003).

The full set of known biochemical reactions for rhamnolipids biosynthesis were added from the P. aeruginosa model (Oberhardt et al., 2008) to the P. putida model iJP962. Where appropriate, stoichiometrically balanced reactions of the rhamnose pathway from KEGG (Kanehisa & Goto, 2000), MetaCyc (Caspi et al., 2016) and BRENDA (Schomburg et al., 2013) were added. Table 1 shows the reactions for rhamnolipids biosynthesis that were added to the P. putida model. Reactions RHLA, RHLB and RHLC represent the rhamnosyltransferase chain A, rhamnosyltransferase chain B and rhamnosyltransferase 2, respectively. Reactions 3H3H and PHAC are involved in the poly (3-hydroxyalkanoic acid) synthase (Oberhardt et al., 2008). The reaction flux across inner and outer membranes was carried out by transport reactions, which were modelled as reactions converting intracellular into extracellular compounds. A transport reaction was also added for the export of rhamnolipids across the cell membrane. For the full list of reactions for rhamnolipids synthesis, exchange and transport see Additional File 1.

The model was then manually curated to establish that the new reactions were fully integrated. This was achieved by evaluating metabolite specificity and metabolite charges accordingly, as well as reaction directionality to confirm that each reaction carried a flux. Gene-protein-reaction rules of the added reactions were also curated from literature (for the final metabolic model see Additional File 2). To run the model, a linear optimization for rhamnolipid production was then carried out, and flux balance analysis (FBA) (see the following subsections for a detailed description) was used to analyze the newly reconstructed engineered model.

Geometric flux balance analysis

Flux balance analysis is a widely used mathematical approach for modelling large-scale metabolic networks (Orth, Thiele & Palsson, 2010). Because FBA assumes the homeostasis of a system, it does not require knowledge of metabolite concentrations and enzyme kinetics. This differentiates FBA from other modelling techniques that require kinetic parameters, usually difficult to obtain. In FBA, the set of biochemical reactions is represented mathematically in the form of a stoichiometric matrix (S) with dimensions of m × n, where the m metabolites are represented a rows and the n reactions are represented as columns. The stoichiometric matrix is a numerical matrix of stoichiometric coefficients for each metabolite participating in a reaction. The stoichiometric coefficient for every metabolite consumed and produced in the system has a negative and positive coefficient respectively. A zero stoichiometric coefficient is given for every metabolite that does not take part in a given reaction. It is assumed that the system is at a pseudo-steady state S · v = 0 that holds for internal metabolites, i.e. those reactants and products of the chemical reactions constituting the model that cannot be imported or exported directly. The vector v represents the flux distribution of the n reactions. Exchange metabolites can be imported and exported from the system, so they do not satisfy the steady state assumption. This flux distribution v therefore represents a feasible flux of metabolites through the reaction network, where under the principle of mass conservation the total amounts of internal metabolite consumed and internal metabolite produced are the same, and the derivative of their concentration is therefore equal to zero.

Constraints such as directionality and capacity (based on enzyme activity, Gibbs free energy change and uptake rates from the literature) are placed on individual reactions by defining the upper (V_max) and lower (V_min) bounds on the range of values that the flux of each reaction can hold (V_min ≤ v ≤ V_max). These constraints define the space of allowable flux distributions at which every metabolite is consumed or produced by each reaction in the system. Despite these constraints, the system is still underdetermined (there are more unknowns than equations), and therefore infinite possible solutions exist. A flux distribution can be obtained by defining an objective function that is a scalar product of the vector of flux rates v, and a vector of weights c, measuring how each component in the network contributes to the production of a biologically desirable phenotype. The set of all possible solutions to the FBA problem is given by the equation and constraints: (1) $$\begin{array}{l}\max c\cdot v,\text{\hspace{0.17em}such that}\hfill \\ S\cdot v=\dot{x}\hfill \\ {\dot{x}}_i=0\text{\hspace{0.17em}if\hspace{0.17em}}{M}_i\in \text{internal metabolites}\hfill \\ {\dot{x}}_i\in ℝ\text{\hspace{0.17em}if\hspace{0.17em}}{M}_i\in \text{exchange metabolites}\hfill \\ V_{\text{min}}\le v\le V_{\text{max}}\hfill \end{array}$$ where in our case the vector c allows to select either the biomass or rhamnolipids as the objective function.

In our pipeline, we use the geometric flux balance approach to define a unique flux balance solution (Smallbone & Simeonidis, 2009). Geometric FBA is based on a geometric representation of a FBA problem. In particular, every FBA problem defines a polyhedron which can be naturally decomposed as the sum of a convex hull and a pointed cone; FBA solutions are to be found within the hull. Since the vertices of the hull and the rays of the cone are uniquely defined, the center of the solution hull (i.e., the final FBA solution) is uniquely defined. Using the geometric FBA algorithm allows us to choose a unique and well-defined flux from the space of all possible solutions. The solution provided also satisfies a number of additional constraints. Indeed, the model assumes that flux correlates with enzyme levels, which is equivalent to the cell minimizing the amount of enzyme required to satisfy this objective. Moreover, the algorithm removes any fluxes representing thermodynamically infeasible internal cycles and selects the solution required to satisfy the given objective from the remaining set of solutions. Hence, the chosen unique solution flux is in a sense “central” and can be considered unbiasedly representative of all possible FBA solutions.

Objective function and uptake rates for optimal rhamnolipids synthesis

Consistent with the reference model iJP962, we used the uptake of glucose at 10 millimoles per gram dry weight per hour (mmol/gDW/h) as a control growth condition. To determine the best carbon source for optimal rhamnolipids synthesis, we investigated alternative carbon sources separately: fructose, sucrose, glycerol, benzoate and myristic acid. We simulated the growth medium with a single carbon source by setting to 10 mmol/gDW/h the uptake of the carbon source under investigation and to zero the uptake of the other sources.

Our P. putida model was optimized to maximize the production and export of rhamnolipids. Hence, we used maximum rhamnolipids production as the objective function in our engineered model. Geometric FBA was used to calculate the optimal flux distribution that maximizes the objective function. Simulations were carried out in MATLAB (version R2018a) using the COBRA toolbox (Schellenberger et al., 2011), with the linear programming solver GLPK (the Matlab script is provided as Additional File 3).

Using gene expression data to build condition-specific metabolic models

Understanding how the transcriptomic alterations change the metabolic phenotype can provide an effective method for data interpretation and analysis (Vijayakumar et al., 2018; Stephens et al., 2015). To this end, we also used the P. putida metabolic model to investigate the transcriptomic effects on different pathways and reactions by including gene expression data into the proposed model.

GEMsplice (Angione, 2018) was used to merge gene expression data with the P. putida metabolic model. The main idea is to create a profile-specific metabolic model for each single gene expression profile. This is done by defining the constraints on fluxes in Eq. (1) as (2) $$V_{\text{min}}\text{φ}\left(\text{θ}\right)\le v\le V_{\text{max}}\text{φ}\left(\text{θ}\right)$$ where the function φ maps the expression level θ of each gene to a coefficient for the lower- and upper-bounds of the corresponding reactions, and is defined as (3) $$\text{φ}\left(\text{θ}\right)={\left[1+\text{γ}\left|\log \text{θ}\right|\right]}^{sgn\left(\text{θ}-1\right)}$$ where the sgn operator returns a vector of ±1 (signs of θ − 1). The constant γ sets the weight of the gene set expression level as an indicator of the rate of production of the associated enzyme (Angione, 2018). We ran our model with γ = 1 in order to ensure a logarithmic effect of the transcriptomic value on the flux bounds of the metabolic model.

We used the integrated model to investigate the relation between gene-expression data and rhamnolipids production in P. putida. We downloaded the expression data of P. putida from GEO (accession number: GSE28257). The dataset provides the expression levels θ of 5,547 genes for 40 samples of the P. putida wild type and 40 samples of P. putida Tn5 mutants. For each sample, a condition specific model was created by using Eq. (2) as constraint in the geometric FBA problem. We used maximum rhamnolipids production as the objective function in our engineered condition-specific models in order to maximize the production and export of rhamnolipids.

Elastic-net regression identifies key genes driving metabolic alterations

After running the condition-specific models, we compared the predicted flux rates (i.e., the FBA solution vector v) of the two groups (wild type and Tn5 mutant) to identify a set of differentially active reactions (DARs), that is, reactions with an adjusted p-value < 0.05. The identified reactions belong to disrupted metabolic pathways that carry a significantly different flux between the wild type samples and the Tn5 mutant samples when the rhamnolipids production is maximized. To further investigate those disrupted metabolic pathways and identify the genes contributing to the flux rates of the DARs in the two groups, we applied the variable selection regression method described below. The idea is to identify the genes that are highly predictive of the flux rate of each DAR when optimizing the rhamnolipids production rate.

Let t be the number of observations (samples) with p predictors (genes). Let y = (y₁,…, y_t)^T be the response (the FBA solution vector) and X = (x₁|…|x_p) be the model matrix (gene expression matrix), where x_j = (x_1j,…, x_tj)^T, j = 1,…, p are the predictors. For any fixed non-negative λ₁ and λ₂, we use the elastic-net regularization criterion (Zou & Hastie, 2005), namely a linear combination of lasso and ridge regression penalties: (4) $$L({\lambda }_1,{\lambda }_2\text{,\hspace{0.17em}β) = }\parallel \mathbf{\text{y}}-\mathbf{\text{Xβ}}{\parallel }_2+{\text{λ}}_1\parallel \mathbf{\text{β}}{\parallel }_1+{\lambda }_2\parallel \mathbf{\text{β}}{\parallel }_2$$ where β = (β₁,…,β_p) is the vector of coefficients to be estimated, $\parallel \mathbf{\text{β}}{\parallel }_1=\sum _{j=1}^p\left|{\text{β}}_j\right|$ and $\parallel \mathbf{\text{β}}{\parallel }_2=\sum _{j=1}^p{\text{β}}_j^2$. The elastic-net estimator $\widehat{\beta }$ is the minimizer of Eq. (4): (5) $$\widehat{\beta }={\text{argmin}}_{\text{β}}\left\{L({\lambda }_1,{\lambda }_2,\text{β})\right\}$$ Let α = λ₂/(λ₁ + λ₂) and λ = λ₁ + λ₂; then solving $\widehat{\beta }$ in Eq. (5) is equivalent to the optimization problem (6) $$\widehat{\beta }={\text{argmin}}_{\text{β}}\parallel y-X\beta {\parallel }_2+P_{\text{α},\text{λ}}\left(\beta \right)$$ where P_α,λ(β) is the elastic-net penalty function defined as (7) $$P_{\text{α,λ}}\left(\mathbf{\text{β}}\right)=\text{λ}\left[\left(1-\text{α}\right)\parallel \mathbf{\text{β}}{\parallel }_1+\text{α}\parallel \mathbf{\text{β}}{\parallel }_2\right]$$

In our analysis, the model matrix X was set equal to the normalized gene expression matrix, with t = 80 observations (expression profiles) and p = 5,547 predictors (genes). The response variable y was set equal to the vector of flux rates of the DAR to be analyzed (y is a vector with dimension t × 1). Hence, y is the vector with the flux rates of a given DAR resulting from running each of the t = 80 FBA condition-specific models. For each DAR, we set the regularization parameter α = 0.5 to achieve a balance between lasso and ridge regression. We used a 10-fold cross validation to identify the optimal λ. Simulations were carried out in R version 3.5.1 using the glmnet package 2.0-16 (Friedman, Hastie & Tibshirani, 2010).

Rhamnolipids synthesis by the engineered model of P. putida

To determine the maximum rhamnolipids production by the engineered in silico model, we investigated several carbon sources and evaluated the metabolic network comprehensively (Fig. 1B). When the P. putida model was optimized for biomass and rhamnolipids production, the amount of rhamnolipids production increased with the uptake of each metabolite. Metabolism of myristic acid (C-14), followed by fructose and sucrose/glucose, provided the best condition for optimal rhamnolipids synthesis. As expected, rhamnolipids synthesis increased in a linear relationship with the increased uptake of various carbon sources (Fig. 1B). Our results also suggest that most of the rhamnolipids production derives from the rhamnose pathway rather than from the FA degradation pathway (Figs. 1A and 1C). Interestingly, when myristic acid was supplied as a carbon source, both pathways contributed to rhamnolipids production to a similar degree. Figure 1D shows rhamnolipids synthesis and biomass by the engineered model under each carbon source.

To determine whether a mixture of metabolites increased rhamnolipids production, we increased the uptake of mixed metabolites simultaneously. When glucose and glycerol or glucose and myristic acid were supplied as combinations of metabolites simultaneously, rhamnolipids production increased to 2.19 mmol/gDW/h and 4.50 mmol/gDW/h, respectively, compared to when each metabolite was supplied individually. In our previous study, a mineral salt medium used for growing biosurfactant producers was initially supplemented with two g/L glucose to initiate biomass production. This was followed by the addition of glycerol to test their influence on biosurfactant production. Pseudomonas aeruginosa DS10-129 produced a maximum of 1.77 g/L rhamnolipid with glycerol at 288 h (Rahman et al., 2002).

Comparison between P. putida and P. aeruginosa for production of rhamnolipids

Figure 1E shows the comparison between our model and the Pseudomonas aeruginosa PAO1 model (Oberhardt et al., 2008) in terms of rhamnolipids production under six different carbon sources (glucose, fructose, sucrose, glycerol, benzoate and myristic acid). The transport reactions for the export of rhamnolipids across the cell membrane were added to the P. aeruginosa model. In order to compare the two models under the same carbon sources, the transport reactions across inner and outer membranes for sucrose and benzoate were also included in the P. aeruginosa model (see Additional File 4 for the full list of reactions).

To analyze the different production rates of rhamnolipids, we investigated the alternative carbon sources separately by setting an uptake rate of 10 mmol/gDW/h for the carbon source under investigation, and zero uptake for the other sources. Both models show a high production rate of rhamnolipids when either glucose, fructose or sucrose is provided, consistent with previous results (Bahia et al., 2018). Glycerol provides enough nutrients for P. putida and P. aeruginosa for the production of rhamnolipids in accordance with Rahman et al. (2002) and Silva et al. (2010). When benzoate or myristic acid were provided as sole carbon source, the production rate of rhamnolipids was 1.967 mmol/gDW/h in P. putida and null in P. aeruginosa, which might be due to the unrelated genome codon index and codon adaptation index profiles of the two bacteria (Weinel et al., 2002). However, if we used the uptake of glucose at 10 mmol/gDW/h as a control growth condition (Oberhardt et al., 2011), the production rate of rhamnolipids was 1.818 mmol/gDW/h in both models.

Principal component analysis reveals biomarkers of rhamnolipids production in P. putida

Principal component analysis (PCA), a form of unsupervised machine learning, identifies data similarities from multidimensional biological datasets (Brunk et al., 2016). More specifically, PCA is a statistical technique that uses a multi-dimensional space to convert a set of correlated variables into linear uncorrelated latent variables called principal components. In our case, it is based on the singular value decomposition of the matrix of flux rates, and is therefore equivalent to finding the system of axes in the space of flux rates such that the covariance matrix is diagonal.

We investigated the individual reactions and identified the key components that drive change in the growth and rhamnolipids production in the engineered model. We applied PCA on our observed flux dataset generated under different growth media; glucose, fructose, sucrose, glycerol, benzoate and myristic acid. To characterize the unique features of individual reactions and variables in the observed flux datasets, we plotted the first two singular vectors of PCA (Figs. 2A and 2B). We found that the first two eigenvectors sum to 88.7% of the variance in the observed flux. These findings suggest that changes in reaction fluxes correlate with the availability of various carbon sources for growth and rhamnolipids production. Figure 2A shows the variable correlation plot of each variable and the contribution for the corresponding carbon source. We found that the first component correlates highly with the variables fructose and sucrose, while the second component correlates with the variables myristic acid and benzoate. Table 2A shows the detailed contribution of each carbon source on the principal components. These variations are driven by changes in the amount of carbon sources used for growth, indicating the network adaptation, particularly in the rate of core metabolic reactions.

Figure 2B shows the top-30 individual reactions, with the highest mean scores on the components, mapped on the first two principal components (for the full list of contributions see Additional File 5). ATP synthase (reaction id: RR08593) and cytochrome-c oxidase (reaction id: IR10022), together with the uptake of oxygen, H₂O and H₂O transport (reaction ids: EX_EC0001, EX_EC0007 and RR08674), scored highly with the first component, indicating the energy demand for growth and rhamnolipids synthesis (Table 2B). The utilization of cytochrome-c oxidase is a common feature of several proteobacteria (Osamura et al., 2017); it is involved in the production of ATP via the respiratory electron transport chain and contributes to the production of the necessary enzymes subsequently used for ATP production by the ATP synthases (Tremblay & Déziel, 2010). One of the reactions that scored highly with the second component is succinate dehydrogenase (reaction id: RR04368). This reaction is involved both in the tricarboxylic acid cycle (TCA) and in respiration via the electron transport chain linked to rhamnolipid production (Wittgens et al., 2011). Figure 2C reports the correlation matrix of the six variables under investigation (the distribution of each variable, the absolute value of the correlation, the result of the correlation test and the bivariate scatterplots with a fitted line). The plot shows that the results of this preliminary analysis are in accordance with the results reported in the PCA variables factor map (Fig. 2A).

In order to analyze the quality of our PCA analysis, we report the scree plot (Fig. 2D) and the cos² correlation map (Fig. 2E). The first two components retain 88.7% of the information (variances) contained in the data, which allows us to focus only them for the statistical analysis of the model. Moreover, the correlation plot of cos² (Fig. 2E) indicates a good representation of the variables on the first two principal components. This also explains the position of the six variables in Fig. 2A (they are close to the circumference of the correlation circle).

In conclusion, the PCA analysis shows that our results are in agreement with those obtained by Wittgens et al. (2011) and Tiso et al. (2016). Indeed, the high cos² value of glucose and fructose shows that they both play a key role in the metabolic pathway of rhamnolipids synthesis chosen as FBA objective. Hence, fluxes from the degradation of glucose and fructose generate rhamnolipids via the rhamnose pathway.

Regression analysis identifies disrupted pathways and genes

We integrated gene expression profiles into the proposed metabolic model of P. putida to investigate disrupted metabolic reactions and pathways using GEMsplice (Angione, 2018). We compared 40 samples of the P. putida wild type with 40 samples of P. putida Tn5 mutants. We simulated the growth medium using the single carbon source that allowed the highest production of rhamnolipids, that is, myristc acid (Fig. 1C). Hence, we set the uptake of myristic acid equal to 10 mmol/gDW/h while the uptake of the other sources was set equal to zero. Table 3 reports the list of the pathways associated with the top-5 DARs, that is, reactions with an adjusted p-value < 0.05 (for the full list of 15 DARs see Additional File 6).

Purine metabolism and FA biosynthesis are the two pathways associated with the top-2 DARs with adjusted p-value of 0.0042 and 0.0043, respectively. Both pathways play a key role in the rhamnolipids production (Rehm, Mitsky & Steinbüchel, 2001) and bacterial membrane biogenesis (Zhang & Rock, 2012). The benzoate degradation via hydroxylation pathway (adjusted p-value = 0.0066) has also been previously linked to the rhamnolipids pathways (Procópio et al., 2012). Indeed, the genes encoding enzymes involved in the rhamnolipids productions also encode enzymes for the benzoate degradation via hydroxylation pathway. It is noteworthy that the identified DARs reflect the disruption of metabolic pathways from the interaction between gene expression profiles (integrated through GEMsplice) and metabolic networks (represented by the P. putida metabolic model). As a consequence, these results are complementary to the outcomes obtained through metabolic network analysis alone, which does not take into account specific transcriptomic profiles.

Regression analysis (elastic-net, see Methods) was then applied to identify the key genes contributing to the flux rate of the DARs. Figure 3 shows the distribution of the top-10 genes (genes with the highest |β|) c the most disrupted DAR (in the purine metabolism pathway) for the wild type and Tn5 mutant samples. By analyzing these distributions it is possible to characterize the metabolic diversity of the different samples and predict their behavior under different conditions. For example, the different distribution of the gene PP2431 might reveal a different cellular adaptation (Fernández et al., 2013). Moreover, the gene PP4355 has been identified as a gene involved in the encoding process of diverse flagellar components in Tn5 mutants samples, which might explain the different distributions in the two types of samples (Sharma et al., 2014). We stress that the procedure proposed here is a single case study, and it can be adapted and extended to identify or compare any two different types of P. putida samples.