Coregulation of Terpenoid Pathway Genes and Prediction of Isoprene Production in Bacillus subtilis Using Transcriptomics

Abstract

The isoprenoid pathway converts pyruvate to isoprene and related isoprenoid compounds in plants and some bacteria. Currently, this pathway is of great interest because of the critical role that isoprenoids play in basic cellular processes, as well as the industrial value of metabolites such as isoprene. Although the regulation of several pathway genes has been described, there is a paucity of information regarding system level regulation and control of the pathway. To address these limitations, we examined Bacillus subtilis grown under multiple conditions and determined the relationship between altered isoprene production and gene expression patterns. We found that with respect to the amount of isoprene produced, terpenoid genes fall into two distinct subsets with opposing correlations. The group whose expression levels positively correlated with isoprene production included dxs, which is responsible for the commitment step in the pathway, ispD, and two genes that participate in the mevalonate pathway, yhfS and pksG. The subset of terpenoid genes that inversely correlated with isoprene production included ispH, ispF, hepS, uppS, ispE, and dxr. A genome-wide partial least squares regression model was created to identify other genes or pathways that contribute to isoprene production. These analyses showed that a subset of 213 regulated genes was sufficient to create a predictive model of isoprene production under different conditions and showed correlations at the transcriptional level. We conclude that gene expression levels alone are sufficiently informative about the metabolic state of a cell that produces increased isoprene and can be used to build a model that accurately predicts production of this secondary metabolite across many simulated environmental conditions.

Figure 1. Quantitation of isoprene production under diverse experimental perturbations inB. subtilis.

Flow chart of the perturbation and sampling protocols. Details are provided in Materials and Methods. B, The isoprene concentrations in the headspace were determined by GC-MS. “OX” represents genetic perturbations in which the listed gene was overexpressed. Bars indicate the mean and error bars indicate standard deviation from the mean.

Figure 2. Terpenoid pathway genes inB. subtilis.

Metabolic pathway layout of terpenoid genes. The 1-deoxy-D-xylulose-5-phosphate (DXP) pathway enzymes (shaded) are: Dxs, 1-Deoxy-D-xylulose-5-phosphate synthase; Dxr, 1-deoxy-D-xylulose-5-phosphate reductoisomerase; IspD, 4-diphosphocytidyl-2-C-methyl-D-erythritol synthase; IspE, 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase; IspF, 2C-methyl-D-erythritol 2,4-cyclodiphosphate synthase; IspG, 1-hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate synthase; IspH, 1-hydroxy-2-methyl-butenyl 4-diphosphate reductase; Idi (encoded by fni), isopentenyl pyrophosphate isomerase; IspA, farnesyl diphosphate synthase; HepS, heptaprenyl diphosphate synthase component I; UppS, undecaprenyl pyrophosphate synthetase; IspS, putative isoprene synthase. The mevalonate (MVA) pathway enzymes (white) are: MmgA, degradative acetoacetyl-CoA thiolase; YhfS, hydroxymethylglutaryl CoA synthase; PksG, 3-hydroxy-3-methylglutaryl-ACP synthase. Metabolite abbreviations: G3P, glyceraldehyde-3-phosphate; IPP, isopentenyl diphosphate; DMAPP, dimethylallyl diphosphate.B, The genome position, length, and sense of each of the terpenoid genes in the B. subtilis genome are represented in the table and graphically.

Figure 3. The terpenoid genes cluster into three main groups with respect to coexpression and into two main groups with respect to perturbations.

Two-dimensional hierarchical clustering of the terpenoid genes against isoprene production in the investigated environmental and genetic perturbations in B. subtilis.

Figure 4. The terpenoid genes cluster into their main groups with respect to coexpression and isoprene production; covariance is dependent upon the type of perturbation.

Genes were clustered using k-means. The normalized isoprene production level is represented by the dashed blue line and the mean expression for each cluster is represented by the solid black line in each cluster. The legend in each of the three panels indicates the name of the genes in each cluster. The mRNA levels are normalized (z-score) for each gene. The x-axis is representative of the perturbations as follows (1 through 12): wild type control, acetic acid, ethanol, lactic acid, indole, 0.005% H₂O₂, 0.02% H₂O₂, OX-dxs, OX-fni, DMSO, empty vector, and OX-ispA.

Figure 5. dxs gene expression has the strongest correlation with isoprene production.

A, Spearman’s correlation coefficients were calculated between gene expression and isoprene production using all 12 perturbations; wild type control, acetic acid, ethanol, lactic acid, indole, 0.005% H₂O₂, 0.02% H₂O₂, OX-dxs, OX-fni, DMSO, empty vector, and OX-ispA. Refer to Figure 2A for gene names. B, Spearman’s correlation coefficients were calculated for the subset of five perturbations: wild type control, acetic acid, ethanol, lactic acid, and indole. Refer to Figure 2A for gene names. C, Genes that positively correlate with isoprene are shown in red, genes that inversely correlate with isoprene are shown in blue, genes with little correlation are not shaded, and genes with an unknown relationship are shown in yellow.

Figure 6. Clustering of the reduced transcriptome identifies similar expression patterns in acetic acid, lactic acid, and H₂0₂ perturbations.

Genes and conditions were clustered using squared Euclidean distance metric. The genes in each cluster are listed in Tables S2, S3, S4.

Figure 7. Isoprene production induced by perturbations can be predicted by a PLSR model based on a reduced transcriptome.

The model was created using 213 genes and trained using ten perturbations with cross validation. The green line is the fit for the training data set; the red line is the fit for the testing data set. The R² value for prediction of the test set is 0.64. Closed circles are representative of training set values that are labeled A through J; the red triangles are the test set conditions, the test set conditions that are identical to the training set conditions are labeled A’ through J’. Unique conditions in the training set are labeled numerically 1 through 12. Perturbation abbreviations: A and A’, 2% acetic acid; B and B’, 2% lactic acid; C and C’, 1% ethanol; D and D’, 0.2 mg/mL indole; E, DMSO; F, OX-ispA; G and G’, OX-fni, H, empty vector; I and I’, wild type control; J and J’, OX-dxs; 1, 0.2% acetic acid; 2, 0.2% lactic acid; 3, 0.02 mg/mL indole; 4, 0.1% ethanol; 5, 0.2% xylose; 6, 0.2% mannose; 7, 0.2% glucose; 8, 0.3 M NaCl; 9, OX-dxsdxr; 10, OX-dxsfni; 11, 0.005% H₂O₂; 12, 0.02% H₂O₂.

Figure 8. Gene expression levels are predictive of isoprene production in a PLSR model.

A total of 213 genes were used in the analyses. Genes (squares), terpenoid genes (red squares), training set perturbations (red triangles), and test set perturbations (closed circles) are indicated in the plot. LV1 and LV2 indicate the latent variable values for each axis, with the percent variance in isoprene captured by the model indicated in parentheses. Perturbation abbreviations for the training set are labeled A through J and two perturbations for the testing data set are labeled as 11 and 12. Abbreviations: A, 2% acetic acid; B, 2% lactic acid; C, 1% ethanol; D, 0.2 mg/mL indole; E, DMSO; F, OX-ispA; G, OX-fni, H, empty vector; I, wild type control; J, OX-dxs; 11, 0.005% H₂O₂; 12, 0.02% H₂O₂.