Molecular Dynamics-Based Allosteric Prediction Method to Design Key Residues in Threonine Dehydrogenase for Amino-Acid Production.

Allosteric proteins are considered as one of the most critical targets to design cell factories via synthetic biology approaches. Here, we proposed a molecular dynamics-based allosteric prediction method (MBAP) to screen indirect-binding sites and potential mutations for protein re-engineering. Using this MBAP method, we have predicted new sites to relieve the allosteric regulation of threonine dehydrogenase (TD) by isoleucine. An obtained mutation P441L has been verified with the ability to significantly reduce the allosteric regulation of TD in vitro assays and with the fermentation application in vivo for amino-acid production. These findings have proved the MBAP method as an effective and efficient predicting tool to find new positions of the allosteric enzymes, thus opening a new path to constructing cell factories in synthetic biology.

Introduction

Allostery is a property in certain proteins, which can be regulated by signals from the allosteric site. The allosteric site is in a different position from the catalytical sites. When the effector binds or disturbs the allosteric site, the allosteric signal can be transmitted from the allosteric site to the catalytic site. Through this remote communication, effectors can regulate proteins with allosteric properties.[1−5] It is worth emphasizing that most of the allosteric proteins are rate-limiting steps in the metabolic pathway.[6,7] Moreover, most of the effective residues of key proteins have been protected by patents.[8] Therefore, finding new and effective residues to reduce the allosteric inhibition is one of the most important problems to be solved in the design of synthetic cell factories. However, for most of the interesting enzymes, new positions that can affect allosteric regulation are difficult to discover and identify.

The characteristic of an allosteric protein is that the allosteric binding site and the catalytic site are located at different positions of the protein. Therefore, how to transduce the regulatory signal to the catalytic site is a crucial issue that determines the occurrence of the allosteric action.[9,10] The traditional method to modify allosteric regulation is to screen mutation sites through time and laborious consuming random mutations. In recent years, directed evolution and high-throughput screening methods have greatly increased the speed. However, the uses of directed evolution or high-throughput screening methods are strictly limited by the detection signal of the enzyme, which greatly limited the application of both methods. With the development of computer technology, it is possible to change the regulation performance of allosteric proteins through computer-aided design. In the previous report, using the combination of the multiple sequence alignment method and an SCA algorithm, allosteric proteins, such as aspartate kinase III and DHDPS, have been successfully rationally designed to relieve the feedback inhibitions by the effectors.[11−14] However, these methods are based on sequence alignments derived from differences in protein families, which have a large degree of randomness in sequence selection. Using protein structure data to discover potential sites has become a new direction for reconstructing allosteric proteins. During the allosteric process, the residues of allosteric molecules that directly and indirectly interact with the protein are particularly important. Molecular dynamics (MD) simulation is a computer-based method that can simulate the interaction between allosteric molecules and proteins in the ns range.[15−18] The MM-GBSA method has provided a possibility to decompose the binding free energy of the effector with the protein into interacted sites.[19,20] Therefore, we proposed an MD-based allosteric prediction method (MBAP) to predict indirectly interacted binding sites on allosteric proteins and their potential mutations, which is based on MD simulation and the MM-GBSA method. As a proof of concept, threonine dehydrogenase (TD) has been used as an example for re-engineering.

TD catalyzes the deamination reaction of threonine to α-ketobutyrate, which is a rate-limiting step in the synthase pathway of the leucine family and is allosterically regulated by its final product isoleucine.[21] TD is of interest in the engineering of the biosynthesis factory because it is a key step in the production of high-carbon biofuels (such as n-butanol, n-propanol, etc.).[22] The TD from Escherichia coli consists of two domains: a catalytic domain at the N-terminus with binding sites for threonine and α-ketobutyrate and an “ACT-like” regulatory domain at the C-terminal. Unlike most other ACT-like enzymes, the regulatory domain of each monomeric TD includes two ACT domains, which means there is a complete ligand-binding site in the regulatory domain of each monomer. The ligand binds at the interface between the ACT domains.[23−25] In this work, the MBAP method has been used to predict the indirect-binding sites of TD and potential mutations to relieve allosteric regulation. Then, the predicted mutations have been verified by experimental tests. Finally, the best mutation has been overexpressed in E. coliMG1655 to test its application in amino-acid production.[26−28]

Results and Discussion

We proposed the MBAP method to predict indirectly interacted binding sites on allosteric proteins and their potential mutations. The first advantage of the MBAP method is the efficient residue selection, which is based on MD simulation and the MM-GBSA method. The second feature is the effective prediction of mutation sites by computer-aided selection. The selected residues are saturation mutated by computer and simulated by MD, and then the binding energy of each mutation is used as a selection factor for comparison. In this paper, we selected TD from E. coli as an example to verify the MBAP method.

MD Simulation of TD

TD is allosterically regulated by its final product isoleucine, and isoleucine binds at the interface of ACT domains. The crystal structure of a TD monomer is available from the PDB database (PDB ID: 1tdj), but there is no published record of the complex structure of TD with isoleucine. In this work, the molecular docking method was used to find the binding location of isoleucine in the enzyme structure. The interface of the ACT domains was set as the center of the binding pocket. Then, the conformation with the lowest binding energy in docking serves as the initial structure for the simulation (Figure a). Isoleucine forms hydrogen bonds with residues Y465 and S467 and is located at the binding pocket formed by residues S362, T365, G442, L447, and L448. This position is similar as reported before.[21]

Figure 1

(a) Initial structure of the complex TD with isoleucine (purple) for simulation, where the binding residues are shown in blue, (b) RMSD of the production run of MD simulation, (c) decomposition of the binding energy into the individual residues of TD by the MM-GBSA method, and (d) residues with energy contribution above 0.1 kcal/mol or below −0.1 kcal/mol.

The first step of the MBAP method is finding potential residues related to allosteric regulation by MD simulation. In this study, MD analysis of 100 ns was carried out, and 15–100 ns data were used for prediction. This period is the critical time that the allosteric effector binds to the allosteric protein and begins to transmit the signal to the catalytic site. The TD complex was neutralized by adding sodium counter ions randomly and solvated in a rectangular box of TIP3P water molecules with a solute–wall distance of 12 Å. Molecular dynamics simulations were performed with a periodic boundary condition in the NPT ensemble using Langevin dynamics at 310 K with the damping coefficient of 5.0 ps–1 and the constant pressure of 1 atm.[29−31] The RMSD of the residue backbones has been calculated during the production simulation process to evaluate the stability of the system (Figure b). During the simulation, the binding energy (Gbinding) between isoleucine and wild-type TD has been calculated, as shown in Table , which is −18.29 ± 3.95 kcal/mol. Then, we used the MM-GBSA method to decompose the binding energy into individual residues (Figure c). Here, the value of Gbinding is the sum of all residues involved in the binding process, including both direct- and indirect-binding sites. In another word, residues that indirectly bind to the effectors can also invest energy during the binding process. Therefore, when the binding energy decomposed into different residues, it is possible to find both direct- and indirect-binding sites of the allosteric protein, which has largely expanded the selection of mutated targets. Meanwhile, the binding energy between wild-type TD and isoleucine has been set as a reference for the selection of mutations. The change of energy is an important signal to reflect whether the mutated residues have affected the binding process or not. In this work, 23 residues in TD, which contributed more than 0.1 kcal/mol to the energy change, have been selected as candidates for further investigations, where 22 provided negative energy changes and 1 promoted positive energy changes (Figure d).

Table 1

Values of Energy from the Binding Process of the Wild-Type TD with Isoleucine in the MD Simulation

type of energy	value of energy (kcal/mol)
Gele	24.46 ± 15.98
GvdW	–22.28 ± 2.1
GSA	2.19 ± 16.1
Ggb	–17.7 ± 13.89
Gnonpolar	–20.47 ± 13.9
Gpolar	6.76 ± 4.11
Gbinding	–18.29 ± 3.95

Prediction of the Key Residues by Computer-Aided Selection

The second step of the MBAP method predicts mutations that can affect the allosteric processes from the selected candidates (Figure a). In this work, computer-aided saturation mutagenesis has been applied to all of the 23 candidates to discover potential mutations, respectively. Each candidate has been, respectively, mutated to 19 other types of amino acids. Then, 1 ns MD simulation has been employed on each mutant with isoleucine. The binding energies have been calculated and compared with the reference. Substantial energy changes at the allosteric binding sites can significantly affect the signal transduction process. Mutations, which significantly reduced the binding energy, were predicted as candidates for controlling the allosteric processes. Finally, 10 mutations with the most significant energy changes have been selected as potential candidates, and seven of them have not been published before (Figure b). With this method, it is possible to find different positions to re-engineer allosteric proteins without disturbing the existing patents. It is a method that uses a computer-based rational design to predict mutations and significantly reduces the burden of work for random mutagenesis to find new candidates. Compared to the mutisequence alignment, our MBAP method does not need a vast number of sequences in the protein family, which always have different lengths and need to be manipulated manually. The MBAP method is able to identify indirect-binding positions that the mutisequence alignment-based hotspot methods cannot reach. Most of these positions are perfect targets to change allosteric regulations without disturbing the existing rights. However, the MBAP method cannot predict sites that are unrelated to binding. Because of this, the locations of the signal transmission process after allosteric binding molecules cannot be predicted by this method. Another limitation is that this method cannot determine whether to increase or decrease the allosteric structure, so further experimental verification is needed to be applied. Therefore, the seven unpublished sites were verified by in vitro enzyme assay.

Figure 2

(a) General steps of the MBAP method to find new targets for allosteric processes. (b) Ten mutations predicted by the MBAP method based on energy changes.

In Vitro Test of the Predicted Candidates for Verification

Finally, the activities of seven candidates were tested under 0–10 mM isoleucine to verify the effectiveness of the MBAP method. As shown in Figure , the wild-type TD was strongly inhibited by less than 0.1 mm isoleucine, and the activity loss was more than 95%. Four of the seven mutants significantly reduced the inhibitory effect of isoleucine. In particular, mutant P441L still maintained 80% activity under 2.5 mm isoleucine and 30% activity under 5 mm isoleucine. Moreover, when isoleucine was not present, this mutant has only about 20% activity loss compared with the wild-type TD. Therefore, the MBAP method can predict the new mutation sites for allosteric regulation. Compared to the traditional methods, the MBAP method is more efficient and effective, and compared to other rational re-engineering methods, the MBAP method focused on both the direct and indirect effectors at the binding sites, which greatly expanded the targets for the redesign of the allosteric binding pocket.

Figure 3

In vitro verifications of the predicted candidates. (a) Plasmid schematic of pET28a-ilvA, (b) SDS-PAGE of purified proteins, and (c) relative activities of TDs with 0–10 mM isoleucine. The absolute activities are 70.5 ± 0.2, 18.8 ± 0.1, 15.3 ± 0.1, 44.3 ± 0.1, 55.5 ± 0.1, 7.3 ± 0.1, and 60.1 ± 0.2 mM/(min·mg protein) for WT, G445M, F440T, G445K, P441L, S443A, and L360Y of TDs. All of the data have been tested in triplicates. The relative activities are calculated by the initial slope of the average data with R2 > 0.99.

In Vivo Verification of the Application of the MBAP Method in Amino-Acid Production

Finally, in vivo experiments were carried out in shake flasks to verify the application of the MBAP method in amino-acid production. In this study, three plasmids were constructed to overexpress the leucine family pathway (Figure ): wild-type and mutant p441l were inserted into ptrc99a, denoted pGFI1 and pGFI2, respectively, and another plasmid pBAPA carrying genes PPC, asd, and aspC was constructed based on the pBbS1C plasmid (Figure a). Finally, three strains of E. coliMG1655 (pBAPA–pGFI1), E. coliMG1655 (pBAPA–pGFI2), and E. coliMG1655 (pBbS1C–ptrc99a) were constructed, respectively.

Figure 4

Metabolic pathway of isoleucine, leucine, and valine from glucose. The allosteric regulations are labeled by orange lines. IlvA is allosterically regulated by isoleucine. IlvI/ilvH is allosterically regulated by isoleucine, leucine, and valine. Three enzymes, PPC, aspC, and asd, are overexpressed in plasmid pBAPA to enhance the metabolic flux toward threonine, which is the substrate of ilvA.

Figure 5

In vivo experiments to verify the application of the MBAP method in amino-acid production. The in vivo performance of the three strains was compared by fermentation in shaking flasks for 112 h. Three strains of E. coliMG1655 (pBAPA–pGFI1), E. coliMG1655 (pBAPA–pGFI2), and E. coliMG1655 (pBbS1C–ptrc99a) were compared in the fermentation. (a) Plasmid schematic of pBAPA and pGFI2, (b) OD600 and glucose consumption, and production of isoleucine (c) and leucine (d) are measured by HPLC during the fermentation. All of the data have been tested twice.

The in vivo performance of the three strains was compared by fermentation in shaking flasks for 112 h. The growth rates of the three strains were similar in the fermentation process. Compared with other strains, the strain with an empty plasmid had a faster glucose consumption rate. Both E. coliMG1655 (pBAPA–pGFI1) and E. coliMG1655 (pBAPA–pGFI2) produced isoleucine and leucine during fermentation (Figure b). Compared with wild-type TD protein, the mutant strain produced 2.2 mg/L isoleucine in the fermentation broth, which was 5% higher than the wild type (Figure c). At the same time, E. coliMG1655 (pBAPA–pGFI2) with ilvA-P441L produced 19.3 mg/L leucine, which was more than 400% higher than that of MG1655 (pBAPA–pGFI1) and MG1655 (pBbS1C–ptrc99a) (Figure d). Therefore, the mutation designed by the MBAP method works well in the in vivo system of E. coli. It shifts the metabolic fluxes from glucose to the isoleucine and leucine pathway, especially enhancing the production of leucine. In the isoleucine biosynthesis pathway, ilvI/ilvH is also allosterically regulated by isoleucine. However, to avoid the confusion of affections by these enzymes, ilvI/ilvH were not overexpressed or re-engineered in this work. For this reason, the production of isoleucine only increased 5% with the overexpression of mutated ilvA compared with the wild-type enzyme. In the fermentation, leucine has been measured at the initial point, which may come from the medium. In the experiment, the seeds were cultured in the LB medium, and the fermentation medium contained 1 g/L yeast powder. Leucine could be introduced from both sources. Due to the rapid growth of the bacteria, leucine in the medium was metabolized rapidly.

Mechanism of Mutation of P441L to the Allosteric Effect of TD

To elucidate the mechanism of mutation of TD-P441L to the allosteric effect, a 100 ns MD simulation has been performed on the mutant protein. (Figure ) The simulation results show that the P441 residue is located in a loop of an isoleucine binding pocket. After isoleucine binding, the loop of the wild-type TD seals the binding region and stabilizes the isoleucine binding site (Figure a). When P441 is mutated to leucine, the loop region becomes open and isoleucine cannot bind to the pocket stably (Figure b). The binding energies have proved the difference. The binding energy of the wild-type TD with isoleucine is −18.29 kcal/mol, while that of the mutant protein with isoleucine is only −10.90 kcal/mol. This result also verifies that the sites around the binding sites play a key role in the stability of effector binding sites. The rational design of these sites can effectively change the allosteric regulation.

Figure 6

(a) Structure of TD with isoleucine at 100 ns of MD simulation; the lysine molecule is shown as a blue stick and the residues that bind with lysine are shown as green sticks. (b) Structure of TD-P441L with isoleucine at 100 ns. The lysine molecule is shown as a blue stick and the residues near lysine are shown as purple sticks.

Conclusions

In this work, we have successfully developed an MD-based method to predict the key residues that participate in allosteric regulation. With the MBAP method, we have successfully found a new mutant P441L that can significantly reduce the allosteric regulation of isoleucine, which also worked in vivo, as shown by the fermentation results. Therefore, the mutant P441L is a good candidate for the construction of the isoleucine- or leucine-producing strain. These results also prove that our MBAP method is an effective predicting tool to find more new positions for the key enzymes to construct high-value biofactories.

Materials and Methods

Materials

l-Threonine, l-isoleucine, l-valine, 2-ketobutyric acid, and pyridoxal-5-phosphate hydrate (PLP) were purchased from Sigma-Aldrich. Isopropyl β-d-thiogalactoside (IPTG) was supplied by Promega (Madison). All of the primers were from Genewiz (Suzhou, China). All of the enzymes and other bioreagents were purchased from Takara (Dalian, China).

Generation of the Structure of the TD–Isoleucine Complex

The X-ray crystal structure of TD was obtained from the PDB database (PDB ID: 1TDJ). The complex structure of TD with isoleucine is not available in the database. In this work, the molecular docking method was used to find the binding location of isoleucine in the enzyme structure, and the docking experiments were employed via Autodock Vina software.[32] For the subsequent docking calculations, the structures were converted into the format of pdbqt using Autodock Tools 4.6[33,34] to merge polar hydrogens. The obtained PDB file of the protein was then processed with Autodock Tools 4.6 and converted into the format of pdbqt. We chose coordinates and dimensions along x, y, and z axes of the grid related to the site of presumed pharmacological interest. In particular, a grid box size of 30 × 30 × 30 Å3 and centered at 36(x), −11(y), and −4(z) was set to generate simulation positions. In the docking experiment, the exhaustiveness value was set to 8. Autodock Vina results were analyzed with Autodock Tools 4.6. Illustrations of the three-dimensional (3D) models of all proteins were generated using PyMOL software.[35]

MD Simulations

The binding processes of the wild-type and the mutated TDs have been calculated using MD simulations. The structure of each complex was generated with PyMOL. The system of MD simulation was constructed by AmberTools 18. The ff99SB and GAFF force fields were applied for the calculation of the complex. The TD complex was neutralized by adding sodium counter ions randomly and solvated in a rectangular box of TIP3P water molecules with a solute–wall distance of 12 Å. After the preparation process, 5000 steps of minimization cycles were carried out to minimize the energy of the water molecules and the counter ions with restriction of the complex. Then, 5000 steps of minimization cycles were carried out to minimize the whole system without restriction. After that, the system was heated from 0 to 310 K by a run for 50 000 steps with a time step of 1 fs. A Langevin dynamics approach with a collision frequency of 1 ps–1 has been used for temperature control. Finally, the simulations of wild-type TD were performed employing NAMD[36,37] software with the AMBER[38] force field. The systems were equilibrated at the target temperature 310 K for 15 ns to relax the complex. After that, the production simulation was carried out for 85 ns. The MD simulations were performed with a periodic boundary condition in the NPT ensemble using Langevin dynamics at 310 K with a damping coefficient of 5.0 ps–1 and the constant pressure of 1 atm. The simulation of the nonbonding pair list was updated every 10 steps, and the particle mesh Ewald (PME) method was used to treat long-range electrostatic interactions. A residue-based cutoff of 12 Å was applied to the noncovalent interactions. No constraint was applied to the protein during the MD simulations. A time step of 2 fs was used, and the coordinates of the simulated complexes were saved every 1.0 ps. Analysis of the MD trajectory was conducted on the entire simulation to ensure the dynamic stability of the system. The backbone root-mean-square deviations (RMSDs) were calculated from the trajectory using the first configuration as the reference, and all coordinate frames from the trajectories were first superimposed on the initial conformation to remove any effect of overall translation and rotation. To examine the convergence of the MD simulations, energy, temperature, and pressure were monitored during simulations.

Free-Energy Decomposition

AMBER software was applied to compute the binding free energy (Gbinding)where Gpolar is the polar interaction energy, Gnonpolar is the nonpolar interaction energy, Gele is the electrostatic interaction energy, GGB is the electrostatic solvation free energy, GvdW is the van der Waals interaction energy, and GSA is the nonpolar solvation free energy. The binding free energy was calculated for each residue and then decomposed to involved residues using the MM-GBSA approach,[39] which is implemented in the Amber program. With this method, residues that have not interacted with isoleucine but composed energy changes have also been selected as indirect-binding candidates. Most of the points with energy values below 0.1 kcal/mol have an error range of about 0.1 kcal/mol. Therefore, low energy changes may be caused by system deviations. Finally, 23 residues with energy contribution above 0.1 kcal/mol or below −0.1 kcal/mol were selected as candidates.

Mutation Prediction

Saturation mutageneses have been applied to each candidate. The simulations of mutated TDs lasted 1 ns with the same method. Binding free energies have been calculated for each mutated protein–isoleucine complex. To find mutants that might affect the allosteric processes, the calculated binding free energies have been compared with statics. The calculated binding free energies were then compared with the binding energy of wild-type TD with isoleucine. Ten residues with the most significant energy changes were selected as potential useful mutations and seven unpublished ones have been verified by experimental test in vitro.

Construction of Vectors and Site-Directed Mutagenesis

The ilvA gene was cloned from the genome of E. coli with PCR reactions. The primers used in the PCR are listed in Table S1. The ilvA gene was fused with the backbone of pET28a by the infusion method, and this constructed expression vector is denoted pET28a-ilvA. The mutations were created based on the backbone of pET28a-ilvA using the site-directed mutagenesis kit from Takara. All of the vectors have been sequenced to verify the inserted gene or mutations. All of the vectors constructed in this work are listed in Table S2.

Protein Expression and Enzyme Activity Tests

Enzymes were expressed in E. coliBL21 (DE3) cells (Takara, Dalian, China) using pET derived plasmids. The recombinant cells were first grown in the LB medium supplemented with 50 μg/mL kanamycin at 37 °C until the OD600 reached 0.6 and gene expression was induced by adding 0.1 mM isopropyl β-d-thiogalactopyranoside (IPTG) for an additional 18 h at 20 °C. The harvested cells were washed twice with 20 mM Tris-HCl buffer (pH 7.0) and suspended in a buffer of 50 mM Na2HPO4 (pH 7.0), 0.2 mM EDTA, and 0.1 mM dithiothreitol. Suspended cells were disrupted by sonication and centrifuged at 100 000g for 1 h. The supernatant was purified using a Ni2+-NTA column (GE Healthcare Bio-Sciences, Piscataway, NJ) to obtain samples for the activity assay. The purity of enzymes was checked by SDS-PAGE (Bio-Rad Laboratories, Hercules) and the protein concentrations were quantified using a Bradford protein assay kit (Bio-Rad Laboratories, Hercules).

The in vitro activity of TD was measured following the previous reports. The reaction of each activity assay contained 20 μM PLP, 50 mM potassium phosphate buffer (pH 7.5), 10 mM threonine, an appropriate amount of the enzyme, and varied concentrations of isoleucine. The reaction was carried out at 30 °C and the formation of α-ketobutyrate was measured at 230 nm using a BioTek Microplate Reader. Each measurement has been repeated three times.

In Vivo Measurement by the Shaken Flasks

The ilvA and ilvA-P441L genes were cloned from pET28a-ilvA and pET28a-ilvA-P441L with PCRs. The primers used in the PCR are listed in SI. The ilvA gene was cloned to the backbone of ptrc99a with the enzyme digestion and ligation method. The constructed vectors are denoted pGFI1 and pGFI2.

The recombinant strains were precultured in a 20 mL LB medium with 0.1 mg/mL ampicillin and 0.017 mg/mL chloromycetin at 37 °C for 14–16 h. The seed cultures were then inoculated into a 50 mL MSI medium Hashiguchi et al.[40] containing 100 μg/mL ampicillin, 17 μg/mL chloromycetin, and 0.1 mM IPTG for batch-fermentation in conical flasks (37 °C, 200 rpm). During the fermentation, the cell growth was monitored by measuring the optical density at 600 nm with a spectrophotometer. Glucose concentration was determined using glucose concentration test kits (Jianhe, Nanjing, China). The concentrations of amino acids were measured by high-performance liquid chromatography (HPLC) (Shimadzu, Japan).