|Home | About | Journals | Submit | Contact Us | Français|
The respiratory cytochrome bo3 ubiquinol oxidase from E. coli has a high affinity ubiquinone binding site that stabilizes the one-electron reduced ubisemiquinone (SQH), which is a transient intermediate during the electron mediated reduction of O2 to water. It is known that SQH is stabilized by two strong hydrogen bonds from R71 and D75 to the ubiquinone carbonyl oxygen O1, and weak hydrogen bonds from H98 and Q101 to O4. In the current work, SQH was investigated with orientation selective Q-band (~34 GHz) pulsed 1H ENDOR spectroscopy on fully deuterated cyt bo3 in an H2O solvent so that only exchangeable protons contribute to the observed ENDOR spectra. Simulations of the experimental ENDOR spectra provided the principal values and directions of the hyperfine (hfi) tensors for the two strongly coupled H-bond protons (H1 and H2). For H1, the largest principal component of the proton anisotropic hfi tensor Tz′ = 11.8 MHz, whereas for H2 Tz7prime; = 8.6 MHz. Remarkably, the data show that the direction of the H1 H-bond is nearly perpendicular to the quinone plane (~70° out of plane). The orientation of the second strong hydrogen bond, H2, is out of plane by about 25°. Equilibrium molecular dynamics (MD) simulations on a membrane-embedded model of the cyt bo3 QH site show that these H-bond orientations are plausible but do not distinguish which H-bond, from R71 or D75, is nearly perpendicular to the quinone ring. Density functional theory (DFT) calculations support that the distances and geometries of the H-bonds to the ubiquinone carbonyl oxygens, along with the measured proton anisotropic hfi couplings, are most compatible with an anionic (deprotonated) ubisemiquinone.
Cytochrome bo3 ubiquinol oxidase (cyt bo3) from Escherichia coli is a member of the heme-copper oxidoreductase superfamily of enzymes, which includes the mitochondrial cyt c oxidase and most prokaryotic respiratory oxidases (oxygen reductases).1–4 Cyt bo3 is located in the cytoplasmic membrane where it catalyzes the reduction of molecular oxygen to water using ubiquinol-8 (UQ8H2) as the electron donor and also functions as a proton pump, conserving much of the energy available from the redox reaction as the proton motive force.5–7 Whereas most of the heme-Cu oxygen reductases utilize cyt c as a substrate, cyt bo3 is in a subgroup that utilizes quinol as a substrate.2,3 Depending on the detergent utilized, cyt bo3 can be co-purified with or without ubiquinone-8 (UQ8).8 The enzyme purified using the detergent n-dodecyl β-D-maltoside (DDM) contains one equivalent of UQ8. This “tightly bound” quinone, located at the high-affinity QH site, is not displaced from the protein even during catalytic turnover using the soluble substrate ubiquinol-1 (UQ1H2).9,10 This has provided support for the existence of two separate quinone binding sites in which the high affinity site (QH) binds a quinone that acts as a cofactor and the low affinity site (QL) binds to the ubiquinol substrate. 9
The QH site stabilizes the one-electron-reduced semiquinone (SQH) detected by EPR spectroscopy.11,12 This ubisemiquinone can be trapped as a kinetic intermediate during the electron transfer process from the substrate ubiquinol to the low-spin heme b component of the enzyme, or can be generated by redox poising the enzyme-quinone complex.10,13–15 The structure of cyt bo3 published at 3.5 Å resolution does not contain quinone.16 The quinone binding site in the protein was identified by locating conserved residues as likely candidates within the structure and was then confirmed by site-directed mutagenesis.16 Residues R71, D75, H98, and Q101 from subunit I were proposed to interact with the bound UQ8 (Figure S1).16, 17 Mutations in each of the four proposed QH-site residues severely reduce the quinol oxidase activity and eliminate the semiquinone EPR signal, with a notable exception being the D75H mutant.17 This mutant stabilizes a SQH radical with a midpoint potential similar to that of the wildtype (WT) enzyme.17 Hence, a protein environment stabilizing the SQ radical is necessary but not sufficient for its proper function. A precise spatial arrangement of the SQ and the surrounding residues at the QH site is crucial for the efficient electron transfer process. The interactions of SQH with the protein environment were extensively studied using X-band 1D and 2D ESEEM (reviewed in18).
Auxotroph strains were developed to perform selective 15N labeling of different nitrogens in Arg, His and Gln residues in the cyt bo3 protein.19,20 2D ESEEM studies have identified Nε of R71 in cyt bo3 as the H-bond donor carrying the most unpaired spin density transferred from SQH.20 In addition, weak hfi couplings with the side-chain and peptide nitrogens from R71, H98, and Q101 were resolved, thus providing a full view of the spin density transfer from SQH to the nearest residues.19,20 These data are supported by the ~10–11 MHz isotropic coupling with the 5′-methyl substituent protons (about twice the value for the SQ in alcohol solutions),21–25 and 13C couplings in 5′-methyl26 and carbonyls,24.all which indicate significant asymmetry in the distribution of the unpaired spin density.
Two strongly coupled exchangeable protons with hfi anisotropic components of 6.3 and 4.2 MHz were found near SQH using X-band 2D ESEEM in conjunction with deuteration of the solvent.25 The 6.3 MHz anisotropic coupling is equal to the value reported for the monoprotonated benzosemiquinone radical generated in frozen alcohol.27 This coupling, along with the isotropic hfi of the methyl protons, were used to support the conclusion that the SQ in the QH site of cyt bo3 is a neutral radical.25 However, the DFT calculations utilizing different simplified models of either the ubiquinone anion-radical or neutral radical H-bonding failed to reproduce all experimentally determined 1H (H-bonds, 5′-methyl), 13C (carbonyls), and 14N (nitrogen H-bond donors) couplings simultaneously.20,28–30 Hence, the SQH protonation state, as well as the general hydrogen bond network, remain uncertain. What is required to clarify the situation is a more precise, experimentally characterized model of the SQH binding in cyt bo3 including the orientations of the H-bonds and the corresponding donor for each strongly coupled proton. This is the goal of the current work.
In order to overcome the existing uncertainties in the structure of SQH we describe a pulsed Q-band 1H ENDOR study on fully deuterated cyt bo3 in an H2O solvent to characterize the hfi tensors of the exchangeable protons. The full deuteration of the enzyme and quinone molecule limits the observed spectra to these protons only. Previous X-, Q-, and W-band 1H ENDOR studies of SQs in proteins and model systems are reviewed in detail.18,31–35 Earlier Q-band ENDOR has been successfully applied for the characterization of the hydrogen bond network around SQA and SQB in deuterated bacterial reaction centers36,37 and the anion and neutral radicals of deuterated p-benzoquinone in water or alcohols.27,38 As in these studies, our Q-band ENDOR work defines the geometries of the H-bonds around the SQH from the analysis of the hfi tensors that is supported by DFT calculations allowing us to directly compare our results with those from the previous studies.
The X-band EPR spectrum of a SQ in frozen solutions is a single line with unresolved g-tensor anisotropy. The spectral width is comparable to the excitation width by microwave pulses, so in the X-band experiment, pulses can be considered as giving a complete excitation of the EPR spectrum. Therefore, at this microwave frequency the ESEEM and ENDOR powder spectra exhibit nuclear frequencies from all orientations of the applied magnetic field relative to the hfi tensor principal axes. On the other hand, at Q-band, with a magnetic field ~1.2 Tesla, the principal components of the SQH g-tensor are sufficiently well resolved, allowing for orientation selective measurements. In these types of experiments, particular sections of the EPR spectrum that correspond to molecules with well-defined orientations with respect to the magnetic field are excited one at a time.27,36–38 The principal values of the 1H hfi tensors and their orientations relative to the SQH g-tensor reference frame for two exchangeable protons are obtained by simulations of the orientation selective Q-band 1H ENDOR spectra recorded at different points on the EPR spectrum. 1H tensors obtained from these measurements were combined with all other EPR parameters referenced above and used as constraints to model the SQH spatial conformation within cyt bo3 and its electronic structure using MD and DFT calculations.
WT cyt bo3 was overexpressed by IPTG induction in E. coli strain C43(DE3) transformed with plasmid pETcyo as previously described.39 Full deuteration of this enzyme was achieved by growing the bacteria in deuterated water with uniformly deuterated glucose as the sole carbon source. To ensure proper growth, cells were inoculated stepwise into the minimal medium with 80%, 90%, and finally 100% deuterated water. This complex was purified with a NTA-Ni affinity column after solubilizing the isolated membrane fraction in 1% DDM (dodecyl β-D-maltoside). After eluting the protein from the column, five successive buffer exchanges with a 100-kDa Amicon Ultra-0.5 mL (EMD Millipore) were carried out to change the buffer to 50 mM potassium phosphate pH 8.3, 10 mM EDTA, 5% glycerol and ~0.05% DDM. The final concentration of cyt bo3 was approximately 500 μM. Sodium ascorbate (10 mM) was added from a freshly prepared 1 M stock. The sample was transferred to a Q-band EPR tube topped with argon, which was further incubated on ice for 3 hours for the SQ signal to develop.9 The tube was then frozen in liquid nitrogen for storage.
The instrumentation and pulse sequences for X-band 1-D 4-pulse ESEEM, HYSCORE, and Q-band pulsed Davies ENDOR have been previously described.25,37,40 All EPR measurements were performed at 60 K. Only the SQH EPR signal is observed in pulsed EPR experiments at this temperature as demonstrated previously.25 The repetition rate in the pulsed EPR experiments was set to 1,000 Hz to avoid saturation of the echo signals. Orientation-selective Davies ENDOR spectra were acquired in stochastic mode with an average of 200 traces.
Q-band 1H ENDOR simulations were performed in Matlab R2014b with EasySpin v5.0.9 (http://www.easyspin.org). The principal values of each rhombic 1H hfi tensor [a−T(1+δ), a−T(1−δ), a+2T] were defined using the isotropic constant a, anisotropic coupling T, and rhombicity parameter δ (δ=0 for axial hfi tensor). The orientation of the hfi tensor with respect the g-tensor reference frame was described using the Euler angles α, β, and γ. These angles are the series of rotations that bring the g-tensor into the hfi tensor eigenframe: (i) counterclockwise rotation around the g-tensor Z axis by α, (ii) counterclockwise rotation around its new Y axis by β, (iii) counterclockwise rotation around its new Z axis by γ to arrive at the hfi tensor.
For orientation selective spectra, the effective excitation bandwidth at each field position can be estimated from the EPR broadening and the length of the microwave pulses. The EPR broadening was modeled using the EasySpin H-Strain parameter in simulations of the Q-band EPR spectrum (Figure 1). The pulse excitation bandwidth was approximated by multiplying the inverse of the initial microwave π-pulse by two for Davies ENDOR. Additionally, the characteristic suppression of small couplings in Davies ENDOR was taken into account by multiplying the ENDOR simulations with an upside-down Lorentzian function with a width (full width at half-maximum) of 1/2tp = 2.1 MHz, where tp is the length of the first microwave π-pulse (240 ns).41 All other parameters were the same as those used in the experiments.
The parameters to be optimized in the simulations included those describing the principal values and directions of the hfi tensors (a, T, δ, α, β, and γ) and the ENDOR linewidth broadening parameter (lwENDOR). Optimizations were carried out by a least-squares minimization using a Nelder-Mead style simplex method.42 Only data below 51.1 MHz and above 52.4 MHz were included in the least-squares fitting. The errors were estimated by shifting the parameters away from the optimized solution until the least-squares deviation rose above a set threshold. This threshold was predetermined by visually inspecting at which point the simulations clearly did not reproduce the general features of the experimental data.
Hyperfine couplings were translated into inter-atom distance constraints based on past DFT calculations, EPR characterization and corresponding crystal structures. Euler angles, on the other hand, were converted into angle/dihedral constraints directly with trigonometry (Table S1).
In conjunction with the EPR experiments, MD simulations were conducted on a membrane-embedded model of cyt bo3 with bound UQ8 to assess how the geometric restraints derived from EPR manifest themselves in binding orientations of the quinone to the oxidase. The starting structure was based on chains A–D from the ubiquinol oxidase structure (PDB code: 1FFT).16 The missing ubiquinone ligand was placed roughly in the middle of the high affinity binding pocket. After adding missing side chains and protons using the VMD43 plugin PSFGEN, the resulting structure was embedded into a bilayer with mixed lipid composition (2:1:1:1 ratio of PMPE:PYPE:POPE:PMPG; a composition that approximately reflects the dominant lipid species in E. coli membranes)44 using CHARMM-GUI. 45,46 The bilayer and inserted protein were solvated using explicit TIP347 water molecules, and ionized to a 150 mM concentration of NaCl using the AUTOIONIZE plugin of VMD.43 The final size of the system was 111×111×129 Å3. CHARMM36 was used for classical force field parameters for the lipid,48 protein,49 hemes, and ionic50 components of the system. Parameters for ubiquinone were obtained from a previous parameterization study.51 Due to the functionalization of one of the heme groups and a copper-containing complex, additional topological components derived by analogy to existing parameters that describe the interactions within the protein were also included (provided in Supporting Information).
In addition to the Hamiltonian formed by the classical force field, additional terms were added such that the sampled configurations are consistent with the EPR experiments, as the crystal structure lacks a quinone ligand. The restraint potentials are composed of half-harmonic potential terms applied to specific bonds and angles within the system (details in Supporting Information). For the strong H-bonds determined between O1 of UQ and Nε of R71 as well as protonated D75 (Figure S1), the distance between the H atom of the donor and the O of the acceptor is harmonically (k = 1000 kcal/mol/Å2) constrained to lie between 1.2 and 2.1 Å. The linearity of the H-bond is maintained by applying a harmonic restraining potential (k = 1 kcal/mol/degree2) if the angle between the hydrogen and its heavy atom is more than 30 degrees away from the line between the donor and acceptor heavy atoms. Weak H-bonds between O4 of UQ and H98 and Q101 are restrained only by the distance between the donor and acceptor heavy atoms: 1.5–4.0 Å for H98 (protonated Nε) H-bond, and 2.2–4.0 Å for Q101 H-bond. In addition, the angle formed by the H of the donor with O1 and C1 of UQ is also restrained (k = 1 kcal/mol/degree2) to a certain range determined by the EPR simulations.
With this setup, NAMD 2.1052 was used to conduct 10 independent simulations. Simulations 6–10 began from the crystallographic state, while simulations 1–5 started after rotating the R71 sidechain to start on the opposite side of the quinone plane. Each copy was equilibrated for 5 ns, during which protein heavy atoms more than 8 away from the binding site were harmonically restrained (k = 1 kcal/mol/Å2) to their initial position. From this equilibrated starting point, data were collected over 50 ns of production simulations where the EPR restraints were active, but the protein was released.
All simulations were performed in an NPT ensemble, using a Langevin thermostat with a damping coefficient of 1/ps to maintain the temperature at 310 K and a Nosé-Hoover Langevin piston barostat53 with period and decay times of 200 fs to maintain the pressure at 1 atm. The barostat acted in a semi-isotropic manner, where dimensions of the periodic cell were held in constant ratio in the membrane plane. The particle mesh Ewald (PME) method54,55 with 1.2 Å grid spacing was used to account for long-range electrostatic contributions during each of the 2 fs time steps.
All density functional calculations, including geometry optimization and hyperfine couplings, were performed using the ORCA package56 and the B3LYP functional and the EPR-II basis set.57,58 Additionally, the optimizations included the third generation (D3) semi empirical van der Waals corrections with Becke-Johnson damping.59,60 The conductor like screening model (COSMO) with a dielectric constant of 8.0 was used for all calculations to replicate the protein environment.61 Geometry optimizations were performed for the three frames selected from the MD trajectory, with all heavy atom coordinates constrained except those belonging to the quinone, in order to optimize the quinone position relative to the residues’ positions in that frame.
A multifrequency CW EPR study of SQH in cyt bo3 was previously performed at X (~9.5 GHz), Q (~34 GHz) and W (~94 GHz)–bands. High-frequency Q- and W-band spectra resolve the g-tensor anisotropy of SQH. Numerical simulations of the W-band spectrum indicate a slightly rhombic g tensor with principal values gx =2.00593, gy=2.00543, gz= 2.00220, giso=2.00452 (error ±0.00005).24
Q-band CW EPR and two-pulse field-swept ESE spectra of SQH in fully deuterated cyt bo3 in an H2O solvent are shown in derivative mode in Figure 1. The overall lineshape of these spectra is predominantly axial with well-separated gx,y() and gz(||) components. In addition, the gz component of the spectra possesses a well-pronounced splitting produced by an exchangeable proton. Two-pulse ESE relaxation rates were found to change with field position, so the CW spectrum was used for simulations. Least-square fitting yields the following parameters for the g and H-Strain values: g = [2.00593 2.00543 2.00228], H-Strain = [11.7 14.8 7.8] (MHz). The simulated g-tensor components obtained here are nearly identical to the values reported previously.24 The splitting at gz was found to be reproduced when the hfi tensor of the hydrogen bond proton with largest components A(1H) = [−6.6 −5.3 11.7] (MHz) (see Table 1) was included in the simulations. The simulated spectrum with the optimized parameters is shown along with the experimental data in Figure 1.
The low-frequency region of the HYSCORE spectrum (Figure 2) of SQH in deuterated cyt bo3 is dominated by a deuterium matrix peak Dwc, which comes from weakly coupled nonexchangeable 2H nuclei of the quinone substituents and amino acid residues. Cross-peaks 1D, located symmetrically on the anti-diagonal relative to the deuterium matrix line, identify strongly coupled deuterium nuclei interacting with the electron spin of SQH. Coordinates of the peak maximum (3.5, 5.0) MHz define a hyperfine coupling ~1.56 MHz for 2H, corresponding to a 10.16 MHz coupling scaled for 1H. This value is consistent with the coupling ~10–11 MHz assigned previously to the 5′-methyl protons based on the HYSCORE and ENDOR spectra of nonexchangeable protons interacting with the SQH spin.21–25
In addition to the deuterium lines, the HYSCORE spectrum contains cross peaks 1N with coordinates (5.1, 3.3) MHz previously observed for SQH in native (not perdeuterated) cyt bo3 and assigned to the Nε of R71.19 These intensive cross-peaks correlate (double-quantum) nuclear transitions of the highest frequency of 14Nε from opposite electron spin manifolds. This nitrogen is involved in the H-bond with SQH and possesses the largest hfi coupling A14N ~1.8 MHz from the delocalization of unpaired spin density from SQH onto the s-orbital of the Nε. On the other hand, the appearance of deep 2H ESEEM in the fully deuterated sample produces intensive peaks in the 2D spectrum that leads to cross-suppression of weaker 14N cross-features and 1H cross-ridges of low intensity in the HYSCORE spectra.
Four-pulse ESEEM provides an alternate way to detect the exchangeable protons through observation of sum combination peaks.63 The 1H sum combination spectrum of SQH in deuterated cyt bo3 in H2O buffer contains three well-resolved lines in the region of the double Zeeman frequency 2ν1H as shown in Figure 3. In comparison with the previously reported spectrum of SQH in a fully protonated background, the intensity of the proton matrix peak 2ν1H is greatly diminished (Figure 3) and is now comparable in intensity with two other peaks shifted from 2ν1H to higher frequencies by ~0.7 and ~1.4 MHz. The shifts of these lines, assigned to strongly coupled exchangeable protons, are the same as in protonated cyt bo3 and correspond to anisotropic hyperfine couplings of ~6.0 and ~4.2 MHz.64 Thus, 1D and 2D ESEEM spectra of SQH in deuterated cyt bo3 with deuterated quinone in a protonated buffer indicate that the hydrogen bonding network around SQH is not significantly perturbed by the deuteration of the protein and quinone molecule.
Orientation selective pulsed Davies 1H ENDOR measurements of SQH were performed at thirteen evenly spaced points in steps of 0.2 mT on the field-swept ESE spectrum (Figure 1). Accumulated ENDOR spectra in absorption and derivative presentation are shown in Figures 4 and and5.5. The sample consists of fully deuterated protein and quinone in an H2O buffer, so only exchangeable protons contribute to the observed spectra. The experimental ENDOR traces exhibit various features of high and low intensity that are resolved from the central 51.72–51.85 MHz 1H matrix region (Figure 4).
Previous analysis of the X-band HYSCORE spectra of SQH in cyt bo3 reported three sets of cross-ridges from exchangeable protons H1–H3 with principal values of the hfi tensors in an axial approximation A = a − T: −7.0, −5.4, 6.3 MHz and A|| = a + 2T: 11.9, 7.2, 1.2 MHz, respectively (signs are relative).25 In agreement with these results, the ENDOR spectra at the gz edge of the EPR spectrum (traces B10–B13) shows a feature of low intensity with a splitting ~11.5 MHz. The splitting and lineshape allows us to assign this feature to the proton H1 with the largest A|| principal value and conclude that the A|| principal axis for this proton is closely collinear with the gz principal axis. This A|| feature resembles a “single-crystal-like-shape” in trace B10. There are also two pairs of intensive features with resolved maxima of A<6 MHz at the gz edge. Those two pairs of lines with the splittings 4.1 and 2.7 MHz possess the best resolution and most symmetrical lineshape in trace B10. The ENDOR spectra recorded at other parts of the EPR spectrum exhibit more complex shape and poorer peak resolution, though up to three maxima are clearly seen in traces B4–B8. The width of this central pattern resulting from the overlap of the most intensive areas around the A features from different protons is increased up to ~9 MHz at the gx,y side of the EPR spectrum, preventing identification of any other remaining A|| lines. The peaks with the splitting 2.3 MHz at the gz edge are present in all spectra except for at the gx,y edge where it is 1.6 MHz. For the simulation of the orientation selective ENDOR spectra we used the representation in first derivative mode (Figure 5), because it better resolves minor variations of the broad features in the traces B1–B9.
The coupling constants a and T determined for the three exchangeable protons from X-band HYSCORE measurements25 were used as initial parameters in the Q-band ENDOR spectral simulations. The inclusion of rhombicity (δ ≠ 0) into all three hfi tensors was found to be necessary to reproduce the experimental line shapes. For protons H1 and H2, approximately axial tensors characteristic of hydrogen bonds fit the ENDOR traces well. On the other hand, a single H3 with previously determined parameters was insufficient to reproduce the weakly hfi coupled ENDOR region (49.5 MHz to 54 MHz). We believe that the ridge previously thought to be solely produced by H3 actually results from contributions from more than one proton, each with T less than 2.9 MHz.25 Indeed, simulations with two protons possessing different anisotropic couplings (H3 and H4) can yield a satisfactory fit to the experimental data in the region around 1H Zeeman frequency. The final ENDOR simulations overlaying the experimental traces in Figure 5 (derivative mode) are shown in red, with the simulation parameters reported in Table 1. The simulation parameters are also demonstrated to reproduce the absorption mode spectra as shown in Figure S2. The hfi coupling parameters agree well with previous reports. Most importantly, the Euler angles of these protons are also determined, and these are crucial to locate the protons with respect to the SQH radical (Figure 6).
With the higher resolution of the hfi tensors of H-bond protons afforded by the full deuteration of the protein and orientation-selective Q-band ENDOR, we are now able to consider their assignments. We will focus on the protons H1 and H2, which possess the largest anisotropic couplings. Results from previous FTIR65 and detailed 15N HYSCORE20 studies allow us to suggest that the two strongly coupled protons belong to H-bonds between the carbonyl O1 and donor atoms of D75 and R71 (Figure S1).
The principal values and principal directions of the hfi tensors provide insight into the locations of the hydrogen bonded protons. Specifically, the Euler angles α and β in Table 1 describe the in-plane and out-of-plane orientations of the TZ′(||) axis with respect to the SQH g-tensor reference frame (see Figure 6 for a pictorial representation of the Euler angles). In contrast to α and β, the Euler angle γ describes the directions of the perpendicular components of the hfi tensor (a − T(1 + δ) and a − T(1 − δ)), and because of the very low values of δ found for H1 and H2, is not considered in this discussion.
When TZ′(||) is assumed to lie along the O···H direction, α and β provide a direct means of determining the direction of the H-bond and, therefore, locating the proton with respect to the carbonyl oxygen. This assumption about the 1H tensor orientation is valid to a good approximation when the dominant contribution to the anisotropic hfi comes from the magnetic dipolar interaction between the unpaired π-electron spin density on the carbonyl oxygen and its H-bonded proton, and is confirmed by our DFT calculations for values up to T ~ 5 MHz. DFT provides a more realistic model of the unpaired electron spin density distribution over the quinone molecule, resulting in a more accurate prediction of the anisotropic hfi tensors. Calculations on the optimized structures of SQA and SQB in bacterial RCs from Rb. sphaeroides were found to be in strong agreement with the ENDOR defined Euler angles.37
Table 2 summarizes the principal values and orientations of the principal axes of the hfi tensors determined by the same Q-band ENDOR approach for the protons of the H-bonds with carbonyl oxygens of the anion-radical of BQ-d4 in H2O,38 the monoprotonated benzosemiquinone BQH•-d4 in 2-propanol (CH3)2CHOD,27 and semiquinones of deuterated UQ10 in deuterated RCs from Rb. sphaeroides in H2O.36,37 Notably, experiments performed with BQ•−-d4 in H2O (and alcohols) defining β ~ 90° and α ~ −54° allowed the conclusions that the geometry of H-bonds is largely in plane with the quinone ring along the lone pair orbitals of the sp2 hybridized oxygen.38 The H-bond protons are characterized by almost purely anisotropic hfi couplings, with T ~ 3 MHz. These experimentally observed values are well supported by DFT calculations.38,66 A hfi coupling of T ~ 3 MHz is consistent with a proton participating in a planar hydrogen bond, forming an angle ±60° with the C=O bond (Euler angle α, refer to Figure 6) and an O···H bond length ~ 1.8 Å.
A Q-band ENDOR study on the monoprotonated benzosemiquinone (BQH•-d4) produced by UV illumination of BQ dissolved in 2-propanol at cryogenic temperatures reported a similar in plane orientation (β=90°) of the TZ′ axis of the anisotropic hfi tensor for the proton in the covalent O-H bond.27 This axis deviated equally (α = −56°) from the C=O bond in comparison to the similar axis for the H-bonded proton. DFT calculations performed on different structures of BQH• coordinated by four, three, or zero 2-propanol molecules have also found the O–H proton lying in the quinone plane and approximately along the direction of the lone pair orbital of the oxygen. The O–H bond length was found to be around 1.0 Å. The largest component of the anisotropic hfi tensor is parallel to the O–H direction (α = −43° to −50°) similarly to the H-bonded protons. The DFT results show very good agreement with the experimental anisotropic hfi components but predict a larger absolute value for the negative isotropic coupling (−8.0 MHz vs − 6.2 MHz from experiment).27
In contrast to the results from the model systems in organic solvent, the angle β characterizing the orientation of the largest component of the anisotropic hfi tensor for the protons of the H-bonds with SQA and SQB in RCs from Rb. sphaeroides is 50–73° (i.e. deviation from the quinone plane is 17–40°). In proteins, one factor influencing the geometry of hydrogen bonds and proton hfi couplings is the structure of the quinone-binding site, particularly the location of suitable H-bond partners for the SQ oxygens. As a result, the hydrogen bonds are likely forced either above or below the ring plane. In addition, the components of the anisotropic hfi tensor are usually larger than the value of T ~ 3 MHz for H-bonded protons with water and alcohol that suggests shorter H-bonds in proteins.18 The extent of the deviation of the hydrogen bond from the quinone plane can influence the isotropic and anisotropic coupling of the hydrogen-bonded proton and should be verified by DFT calculations.30,67
In light of this discussion, the hfi tensor H2 in cyt bo3 possesses characteristics typical for an H-bonded proton between the SQ and protein residues: T ~ 4.3 MHz, a small isotropic coupling, and an out of plane deviation of the TZ′ principal axis (and O···H direction) of about 25°. The in-plane deviation of the H-bond from the lone pair direction is close to zero (lone pair direction is 60°, while the α for H2 is 55°) though it is noted that the angle α is determined with significantly lower accuracy than angle β (error margins are 2 to 3 times larger in Table 1). In contrast to H2, the components of the anisotropic hfi tensor for H1 are closer to the tensor of the proton in the covalent O-H bond. Together with the large isotropic coupling for methyl protons these data best fit the description of a monoprotonated, neutral radical for SQH.25 However, there remains a need for additional modeling of SQH to support this conclusion, particularly due to the significant difference in the isotropic constant with the value reported for BQH•, as well as the unusual orientation of the O···H bond in the direction almost normal to the quinone plane (β ~ 20°). The lack of bound quinone in the existing X-ray structure prevents DFT computations of the SQH electronic structure. Therefore, at the first stage of our theoretical analysis we have performed equilibrium MD simulations on a membrane-embedded model of cyt bo3 ubiquinol oxidase with bound ubiquinone.
The EPR measurements place restraints on allowable binding geometries of SQH and the oxidase. MD simulations of the quinone and oxidase together places this geometric information of protons into the context of the protein pocket to obtain structural insights of SQH binding. As discussed in Experimental Procedures, MD simulations of the bo3 oxidase with a quinone bound to the “high affinity” pocket are carried out with additional harmonic constraints enforcing structural features derived from pulsed EPR. These biased simulations are designed to mimic the situation of SQ binding, and set up the initial geometries for subsequent DFT calculation. Additionally, the flexibility of working in silico can permit different hypotheses to be tested as to the identity of the H1 and H2 protons. It is known that the strongly H-bonded protons H1 and H2 are on the O1 side of the SQH, however, their assignment to residues D75 and R71 is ambiguous.18,20 Parallel MD runs were set up with proton H1 being assigned either to the β-carboxyl of D75 or to the Nε of R71 to test both possibilities.
By the end of each 50 ns-long MD trajectory, the protein RMSD is commensurate with the crystal structure resolution in all 10 trajectories (Figure S3). To facilitate the evaluation of the quinone conformation, frames from the trajectories were grouped into 4 clusters based on the structural similarity of the quinone binding site under restraint (Figure S4). The representative structures from each group are shown in Figure S5. It is clear that these structures are quite different from each other, suggesting that accessible quinone conformational space is quite large. These restrained conformations cover a wide distribution of the relevant bond lengths and angles (Figures S6 and S7) defining the geometry of SQH within the protein binding pocket. From the MD trajectories, some general structural features of SQH binding are observed. For instance, the 2-methoxy group tends to interact with D75, and the two H-bond donors R71 and D75 most likely are located on opposite sides of the SQH plane. The MD runs used previously determined parameters for the fully oxidized quinone,51 which cannot alone correctly represent SQH. This prevents us from using frames directly from MD to evaluate the SQH binding conformation. We have further performed structural optimization using DFT on selected candidate structures. These structures were chosen since they simultaneously satisfied all bond lengths, bond angle and dihedral angle constraints derived from EPR. Three frames closely satisfying all the geometrical parameters determined by the EPR simulations are shown in Figure 7. Representative distances and angles are listed in Table 3. The side chains of D75 and R71 are always located on opposite sides (above or below) of the quinone plane and the H-bonds from these two residues are quite short (~2 Å O···H distance). In contrast, the H-bonds between H98/Q101 Nε and the O4 carbonyl are considerably longer and do not exhibit a particular spatial pattern.
Initial models for DFT calculations were generated using the three frames selected from the MD trajectories with hydrogen bonds from the carboxylic acid OH group of D75 and the NεH guanidium group of R71 to O1, and with the NεH of the imidazole group of H98 to O4. By comparing different models, the calculated values can provide valuable insight into understanding how structural variations between the three different frames will affect the resulting EPR parameters. Table 4 shows the calculated principal values for the anisotropic hfi tensor and isotropic coupling for two H-bonded protons from residues R71 and D75, which are the main focus of this study. Comparison of the values for the three frames indicates that calculated hfi tensors for both H-bonded protons are consistent with experimentally determined tensors H1 and H2 for MD2 only. The largest calculated principal value 11.3 MHz for the H-bonded proton of D75 and 8.4 MHz for the proton of R71 are in a very good agreement with the experimental values 11.8 MHz (H1) and 8.4 MHz (H2), respectively. In addition, the calculated tensors possess low hfi tensor rhombicity δ=0.06 (D75) and 0.15 (R71) comparable with that estimated from spectral simulations δ=0.1 for H1 and H2.
In the DFT optimized model the H2 proton is significantly out of the quinone plane and located practically under the O1 atom, i.e. β = 180° (Table 3). This angle is substantially different from the EPR spectral analysis (β ~ 115°). We have investigated the dependence of the calculated hfi couplings on this angle. As shown in Table S2 rotation of this angle for the optimized DFT model leads to only minor changes in calculated hfi couplings for the hydrogen bonded protons. In contrast, a more planar H2 proton (β = 170°, 160°, 150°) possesses progressively smaller 1H and 13C hfi couplings for the 5′-methyl group, yielding a better agreement with experimental values. However, it is impossible to evaluate the conformation dictated by EPR simulation (β ~ 115°) here because at 115° the R71 sidechain has severe steric clashes with the 2-methoxy group which is pointing down in this DFT model. In summary, judging from the DFT calculations so far, it appears that H1 belongs to D75.
The small values of anisotropic hfi couplings for protons H3 and H4 (~1 MHz) indicate weakly coupled exchangeable protons in the SQH environment. The calculations performed in our previous work20 exploring idealized small models of the SQH–protein interactions with geometry optimization have shown couplings of this order for the Nη protons of the R71. Segments of O-H directions defined by Euler angle uncertainties for these protons are consistent with the location of NηH2 in the MD2 (Figure 7).
Calculated hfi couplings can also be compared with previous experimental determinations of 1H (10–11 MHz)21–25 and 13C (−6.1 MHz)26 isotropic constants for the 5′-methyl group and 13C anisotropic hfi tensors and isotropic constants for C1 (−8.9, −17.3, 26.1) MHz, a=4.7 MHz and C4 (−7.9, −11.3, 19.3) MHz, a = 0.9 MHz carbonyl atoms in SQH, respectively.24 All these characteristic hfi couplings have been measured in several studies of SQs and indeed are the principal indicator that the spin density distribution in the SQH is highly asymmetric compared with a practically symmetric distribution in the anion radical generated in water or alcohol solutions.18
DFT calculations for the 1H and 13C isotropic couplings of the 5′-methyl group of SQH, listed in Table 4, are in good agreement with experimentally determined values for MD1, but are elevated (~2.5–3.0 MHz for 1H and ~1.5 MHz for 13C) for MD2 and MD3. In contrast, the fit is poorer for 13C carbonyl tensors. It is noted that the experimental hfi tensors for 13C1 or 13C4 labeled quinone were measured from Q- and W-band EPR spectra in frozen samples, where the hfi splittings are poorly resolved only at the high-field component with gz~2.0022.24 No any other splittings are resolved. The complete tensors were obtained indirectly from simulation of the axial spectral lineshape. Such analysis provides, at best, a reasonable estimate of the largest component of the hfi tensor only. Particularly, it clearly indicates that Azz of C1 is larger than that of C4.
Our 13C, C1 and C4 computations for frames MD1–MD3 show greater asymmetry than observed experimentally for SQH. For example, if we consider the computed Az=a+Tzz components of the hfi tensors ~39–58 MHz for C1, and in particular the very low |Az| ~0.9–12 MHz for the C4 position, the experimentally reported values are ~30 MHz for C1 and ~20 MHz for C4 in SQH.24 The better agreement of the DFT calculations with the hydrogen bonded and methyl group protons described above would indicate that the inferior 13C agreement for C1 and C4 is probably due to computational overestimation of the CO polarization which historically has been a difficult to predict using DFT calculations. Previous reports28,30 indicated relatively better agreement between experiment and computation for 13C1 and 13C4 hfi tensors. These models used, principally, neutral water molecules as H-bond donors to both O1 and O4. This type of neutral and symmetrical H-bond donation to both oxygens generally leads to less polarization of the CO groups compared with our charged one-sided models obtained from the MD simulations.
Results of our calculations may be understood from the Mülliken unpaired spin populations for MD2 (Figure S8), which places relatively large spin density on O4 (0.30) but very low spin density on C4 (0.04). Difference between spin populations of O1 (0.13) and C1 (0.23) is not so significant. The experimental 17O couplings for the carbonyl oxygens in SQH are thus far not available. However, the difference between the calculated Az of 17O1 (~ −60 MHz) and 17O4 (~ − 110 MHz) in MD2 exceeds the largest difference between 17O1 (~ −74 MHz) and 17O4 (~ −105 MHz) calculated previously for the anion-radical SQH model with three H-bonds to O1 (3/0, 1HO-1HO-1HN) among other different models.28 All demonstrate the highly asymmetric nature of the spin density distribution on the SQH which was originally described and analyzed in detail in ref. 20.
Our previous assignment of SQH as a neutral radical was based on (i) The high anisotropic coupling 6.3 MHz for one of the exchangeable protons corresponding to a short, in plane O···H bond and (ii) An isotropic coupling constant ~10–11 MHz for the protons of the 5′-methyl substituent, which is two-fold larger than that of the anion-radical with symmetrical hydrogen bonds.25 In the current work, orientation selective Q-band ENDOR spectra of exchangeable protons near SQH allowed us to determine the orientation of the anisotropic hfi tensor for the H1 proton in the g-tensor coordinate system. The principal axis of its largest component, coincident with high accuracy with the O···H1 direction, forms an angle β ~ 20° with the normal to the quinone ring plane. The O···H1 bond direction is nearly perpendicular to the plane of the quinone ring, unlikely to be a covalent bond. The O···H2 bond direction is also out of plane, but only by about 25°.
From the combined MD and DFT computational efforts, the H1 proton is assigned to the H-bond between the SQH carbonyl O1 and D75 and the H2 proton is assigned to the H-bond between the same O1 and R71 Nε. The O···H1 and O···H2 distances in the DFT optimized structure of MD2 are equal to 1.61 and 1.74 Å, respectively. These distances are within the typical range of H-bonds to an anionic SQ and significantly exceed the length of the covalent O-H bond.18 The distances and angles above allow us to conclude that the SQH stabilized in cyt bo3 is the anionic SQ. In addition, the significant asymmetry of the spin density distribution in SQH revealed by 1H and 13C 5′-methyl couplings21–26 and the 13C hfi tensors for C1 and C424 are consistent with this H-bond geometry. The unusual characteristics of H1 result from the strong ~70° out of plane H-bond deviation. The asymmetry of the spin density distribution is likely caused by the additive influence of both H-bonds between carbonyl O1 and D75 and R71 Nε, because DFT calculations give similar values for 1H and 13C 5′-methyl and 13C carbonyl couplings for all Frames.
The assignment of SQH as an anionic semiquinone also agrees with the conclusion based on the optical absoption spectrum of the formation of SQH following pulsed radiolysis of cyt bo3.15 The pH-dependence of the potentiometric titration of SQH indicates that the pKa of the semiquinone Q•−•/QH is about pH 7.3,12 also consistent with the observation of the anionic form of SQH in the current work, performed at pH 8.3.
This research was supported by Grants DE-FG02-08ER15960 (S.A.D.) and DE-FG02-87ER13716 (R.B.G.) from Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Sciences, US DOE, and NCRR/NIH Grant S10-RR15878 and S10-RR025438 for pulsed EPR instrumentation. The computational work at the UIUC was supported by the grants from the NIH to E.T. (U54-GM087519 and P41-RR05969). PJOM acknowledges the use of computer resources granted by the EPSRC UK national service for computational chemistry software (NSCCS). A.T.T. and J.V.V. gratefully acknowledge past support as NIH trainees of the Molecular Biophysics Training Program (5T32-GM008276). J.V.V. also acknowledges support from the Sandia National Laboratories Campus Executive Program, which is funded by the Laboratory Directed Research and Development (LDRD) Program at Sandia National Laboratories. Sandia is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the US Department of Energy’s National Nuclear Security Administration under Contract No. DE-AC04-94AL85000.