Electric Field Measurements Reveal the Pivotal Role of Cofactor–Substrate Interaction in Dihydrofolate Reductase Catalysis

The contribution of ligand–ligand electrostatic interaction to transition state formation during enzyme catalysis has remained unexplored, even though electrostatic forces are known to play a major role in protein functions and have been investigated by the vibrational Stark effect (VSE). To monitor electrostatic changes along important steps during catalysis, we used a nitrile probe (T46C-CN) inserted proximal to the reaction center of three dihydrofolate reductases (DHFRs) with different biophysical properties, Escherichia coli DHFR (EcDHFR), its conformationally impaired variant (EcDHFR-S148P), and Geobacillus stearothermophilus DHFR (BsDHFR). Our combined experimental and computational approach revealed that the electric field projected by the substrate toward the probe negates those exerted by the cofactor when both are bound within the enzymes. This indicates that compared to previous models that focus exclusively on subdomain reorganization and protein–ligand contacts, ligand–ligand interactions are the key driving force to generate electrostatic environments conducive for catalysis.


Circular dichroism analysis
Structural characterisation of the enzymes was carried out with Applied Photophysics circular dichroism spectroscopy. Comparison of the molar ellipticity of the wild type enzymes with the variants reveals that introduction of the nitrile probe into the active site have negligible effects on their structures. Figure S1. (A) CD spectra at 20 °C of cysteine free EcDHFR (blue), EcDHFR S148P/T46C (green) and EcDHFR S148P/T46C-CN (red) and (B) their thermal melting curves. CD spectra were recorded in 10 mM potassium phosphate buffer (pH 7.0) at a protein concentration of 10 μM.   BsDHFR T46C-CN 12.1 ± 0.9 123 ± 6 1.9 ± 0.6 82 ± 5 * data obtained from reference 3 Figure S10. Graphical steps of spectra treatment employed for extracting the maximum vibrational frequency of labeled DHFRs. Reported peaks were averages of the three spectra treatments; minimum (2 nd derivative), maximum (Gaussian fit), and inflexion point (1 st derivative) from at least two independent measurements. Error is reported as the standard deviation of the treatments.

NMR spectra analysis
The NMR spectroscopy samples were prepared according to the previous report with 2.0 mM of S 13 CN labeled protein, 10 mM ligands with 2.5 mM of 3-(trimethylsilyl)-1-propane sulfonic acid sodium salt as an external standard.

Gas phase calculations
Geometry optimization and frequency calculations of the nitrile bond of MeCN and MeSCN were carried out with Gaussian09 package of programs in the gas phase with three semiempirical methods and 18 density functionals ( Figure S12). 5 Frequency calculations were carried out with the harmonic and the anharmonic approximations, as reported in Table S4.

System Preparation
The calculations of the closed Michaelis complex, E:NADP + :folate of EcDHFR were carried out from the initial X-ray structures as deposited in the PDB with code 4P66. This structure contains a nitrile probe (XCN) attached at position 46, as well as C85A, D37N, and C152S mutations. NADP + and MTX (methotrexate) were in the active site. MTX was manually modified to FOL (NA4 mutated to O4, and CM removed). The missing side chain for Glu17, Asn18, Asn23, Arg33, Arg52, Gln108 and Glu118 where added.
The protonation states of titratable residues were determined with PROPKA ver. 3.0 3. 6 Accordingly, His149 was protonated in δ position, His124 in ε position and His 45, 114 and 141 in both δ and ε positions. The hydrogen atoms were then added considering pH 7 with the utilities of fDYNAMO library. 7 Finally, the system was solvated by placing it in a 100 ´ 80 ´ 80 Å 3 pre-equilibrated box of water molecules. Any water with an oxygen atom lying in a radius of 2.8 Å from a heavy atom of the protein was deleted. The geometry of the system was optimized applying a series steepest descent conjugated gradient minimization (gradient tolerance of 5 kJ·mol -1 ·Å -1 ) followed by a series of L-BFGS-B (gradient tolerance of 0.1 kJ·mol -1 ·Å -1 ). A truncation function (shift function) was applied to treat the non-bonded interactions, with an internal cut-off of 14.5 Å and an external of 16 Å. All the optimization algorithms were implemented in the fDYNAMO library. Then, 5 ns of classical MD simulations were run to equilibrate the system with AMBER force field, 8 as implemented in NAMD software. 9 Parameters for both substrates obtained with Antechamber 10 at AM1 level. Time evolution of the RMSD on the 5 ns of classical MM molecular dynamics simulation in the Michaelis complex (E:NADP + :folate) is displayed in Figure S16B, while B-factor of the backbone after the 5 ns of classical MM molecular dynamics simulation is shown in Figure S17.
In the case of the preparation of the occluded product ternary complex, E:NADP + :THF, the initial structure was obtained from PDB 1RX6. In this case, the three mutations have been introduced: C85A/C152S/T46SCN + D37N (names of last atoms in SCN SG-CSNC). The active site contains NADP + and DDF (5,10-dideazatetrahydrofolic acid). This later was modified to THF (C10 mutated to N10; C5 mutated to N5) while the missing ring was added in the NADP + . Same as in the closed conformation, the protonation state of the titratable residues of the EcDHFR was determined with PROPKA ver. 3.0 3. Accordingly, His149 was protonated in δ position, His124 in ε position and His 45, 114 and 141 in both δ and ε positions. Once the full protein was prepared, it was solvated using the same procedure as in the closed conformation and, after the same series of optimization algorithms, 5 ns of classical MM MD simulation was carried out. The time evolution of the RMSD on this occluded conformation is displayed in Figure S16C, while B-factor of the backbone after the 5 ns of classical MM molecular dynamics simulation is shown in Figure S17.
Finally, the holoenzyme binary complex, E: NADPH, was prepared from the closed Michaelis complex, E:NADP + :folate, by removing the substrate and running the corresponding 5 ns of classical MD simulations for equilibration. The time-dependent evolution of the RMSD is displayed in Figure S16A.

Nitrile frequency calculations.
In order to carry out the frequency determination of the nitrile bond, 500 structures were randomly selected from the last 500 ps of the classical MD simulations. Then, a hybrid quantum mechanics/molecular mechanics (QM/MM) scheme was used, where the probe was described by BVP86 functional with the standard 6-311++g(d,p) basis set, and the rest of the protein and water molecules with OPLS-AA 11 and TIP3P 12 classical force fields, respectively. A detail of the active site showing the QM/MM partition is presented in Figure S18.

Quantification of the hydrogen bond interactions
The contributions of H-bond interactions between the nitrile group of the probe and the active site of the EcDHFR in both states were computed according to the geometrical distance and angle function of equation S1, as proposed by Jensen and co-workers. 13

Equation (S1)
Distance contribution:   Table S5. QM/MM theoretical calculation of the frequency of the nitrile bond (n), total electric field projected on the nitrile bond ( ⃗ !"" ) and its contributions coming from the cofactor ( ⃗ #$%& ). Values averaged over structures selected from the last 51 ps of the full classical MD simulations. The hydrogen bond interactions (HB) were quantified on representative structures using equation S1.