Detecting O2 binding sites in protein cavities

Internal cavities are important elements in protein structure, dynamics, stability and function. Here we use NMR spectroscopy to investigate the binding of molecular oxygen (O2) to cavities in a well-studied model for ligand binding, the L99A mutant of T4 lysozyme. On increasing the O2 concentration to 8.9 mM, changes in 1H, 15N, and 13C chemical shifts and signal broadening were observed specifically for backbone amide and side chain methyl groups located around the two hydrophobic cavities of the protein. O2-induced longitudinal relaxation enhancements for amide and methyl protons could be adequately accounted for by paramagnetic dipolar relaxation. These data provide the first experimental demonstration that O2 binds specifically to the hydrophobic, and not the hydrophilic cavities, in a protein. Molecular dynamics simulations visualized the rotational and translational motions of O2 in the cavities, as well as the binding and egress of O2, suggesting that the channel consisting of helices D, E, G, H, and J could be the potential gateway for ligand binding to the protein. Due to strong paramagnetic relaxation effects, O2 gas-pressure NMR measurements can detect hydrophobic cavities when populated to as little as 1%, and thereby provide a general and highly sensitive method for detecting oxygen binding in proteins.


Oxygen-induced spectral changes.
N chemical shifts are about 0.1-1.0 ppm. In addition, cross-peaks for Y88, I100, and L118 became weaker or disappeared with increasing O 2 concentration. In contrast, N 2 and Ar gas did not induce perceptible changes in chemical shifts and cross-peak intensities in the same ranges of gas concentrations (N 2 ~3.3 mM, Ar ~7 mM; see Supplementary Fig. S2 online). Figure 1b shows the region for methyl group signals in 1 H/ 13 C constant time (CT) HSQC spectra of 15 N/ 13 C-labeled L99A. O 2 -induced chemical shift changes and/or loss of signal intensities are significant for the methyl groups of I78δ 1 , I78γ 2 , L84δ 1 , L84δ 2 , V87γ 2 , A99β , M102ε , V103γ 2 , V111γ 1 , V111γ 2 , L118δ 1 , L118δ 2 , L121δ 1 , A129β , A130β , L133δ 1 , and I150γ 2 . In contrast, chemical shift changes were not observed when Ar was increased to 5 bar (~7 mM) or N 2 was increased to 7 bar (~4. 6 Figure 3 shows the O 2 -induced chemical shift changes observed for amide nitrogens, methyl carbons, and methyl proton nuclei of the residues around the enlarged hydrophobic cavity (cavity 4). The O 2 association constant can be estimated from the concentration dependence of peak positions, assuming exchange between O 2 -bound and free states. Global fitting to all chemical shift changes was performed, using the following equation (1),  Supplementary Table S1 online. Figure 6 shows O 2 -induced Δ R 1 for methyl protons of the C-terminal domain (i.e., residues 71-160). Methyl protons exhibiting marked PREs (Δ R 1 ≥ 4 s −1 ) were also located around the two hydrophobic cavities, as shown in Fig. 5b. Note that R 1 values of several methyl protons around the two cavities could not be obtained at 3.8 mM of O 2 concentration, because their cross-peaks were severely broadened or disappeared (see section: The paramagnetic effect leads to line broadening). These PRE data closely match the O 2 -induced changes in chemical shifts and peak intensities (Fig. 1). The PRE arises from dipolar interactions between a nucleus and unpaired electrons of the paramagnet and spin relaxation contributions show < r −6 > distance dependence between the paramagnetic center and the nucleus of interest undergoing rotational motion, as described by the Solomon-Bloembergen equation 35,36 . O 2 -induced Δ R 1 for amide protons predicted by 1/r 6 -weighted distance analysis is shown in Fig. 4a 13 . The crystal structure of L99A at 8 atm of xenon pressure (PDB ID, 1c6k) has three xenon atoms in cavity 4 but none in cavity 3 10 . Therefore, we added two xenon molecules (the maximum number of Xe that can be accommodated) to cavity 3 and minimized the total energy. Δ R 1 was estimated from 1/r 6 distance dependent PRE contribution from each xenon site, using equation ( where r 1-3 are the distances to xenon binding sites 1-3, respectively, in cavity 4 and r 4-5 are the distances (Å) to the xenon binding sites 4 and 5, respectively, in cavity 3, and a, b, c, d, e, and f are fitting parameters. Hydrogen atoms were added to the crystal structure by using the WHATIF server 37 . A linear combination of predicted PREs from the five xenon-binding sites matches the observed Δ R 1 pattern well (Fig. 4a) Table 1. Chemical shift changes of representative amide nitrogen, methyl protons, and methyl carbons for the oxygen binding to L99A. a Δ δ are obtained by a global fitting for changes in chemical shifts using eq. 1.
shows the relative O 2 occupancy at each xenon site. O 2 occupancy in cavity 3 is about 5% of that in cavity 4. In addition, extremely large Δ R 1 values were predicted for the amide protons of Y88 and L118, which indeed show severe line broadening with increasing O 2 concentration. The R 2 -value of the correlation between observed and predicted Δ R 1 was 0.82. Although O 2 occupancies at each xenon site are estimated above, contributions of each binding site to Δ R 1 of each amide proton are expected to be different, as borne out by Fig. 4b. For instance, the amide protons of residues 129, 130, 153, and 154 exhibit large PREs from site 5 due to their proximity to bound oxygen molecules, even though the O 2 occupancy at site 5 is much smaller than that for sites 1-3. The contributions to Δ R 1 from the individual binding sites as a function of distance are given in Supplementary Fig. S5 online. The data can be adequately modeled with a 1/r 6 distance dependence. O 2 -induced Δ R 1 for methyl protons of residues 71-160 was also predicted by the 1/r 6 -weighted distance analysis in Fig. 6. The R 2 -value of the correlation between observed and predicted Δ R 1 was 0.80. These statistically substantial correlations indicate that O 2 molecules associate with cavities 3 and 4. The parameters a, b, c, d, e and f obtained were 1.8, 0.0, 3.1, 0.0, 0.046 (Å 6 /s), and 0.79 (s −1 ) respectively. Standard error of the estimate was 3.0 (s −1 ). The predicted O 2 occupancy in cavity 3 is much smaller than that of cavity 4, suggesting that O 2 binding to cavity 3 is weaker than that to cavity 4. Furthermore, O 2 occupancy at site 2 appears to be lower than that at sites 1 and 3. Although these tendencies are consistent with the results of amide protons, the ratios of parameters a-e between amide and methyl protons are different. We considered two explanations. First, the number of the data of methyl protons we used in the fit is smaller than that for amide protons. Even at a few bars of O 2 pressure, line broadening by PREs is sufficiently large to prevent a correct estimate of Δ R 1 (see section The paramagnetic effect leads to line broadening). Indeed, the standard error of estimate in the case of methyl protons is four times larger than that of amide protons. Alternatively, the values of the predicted Δ R 1 increase very steeply, as the distance is less than 3 Å to the xenon site. For instance, even within the same methyl group, the values for each proton are largely different. Thus, in order to obtain a more quantitative prediction, we need to know the exact location and probability of O 2 molecules at these binding sites much more precisely to calculate a correct estimate.
According to a calculation by Teng and Bryant, the O 2 -induced proton relaxation rate constant for a proton in van der Waals contact with O 2 is about 6200 s −1 if O 2 is bound at all times 12 . Accordingly, an observed value of 10 s −1 for Δ R 1 suggests that the binding probability is ~0.2% or less. However, the O 2 -bound state probability    Supplementary Table S1. Severe line-broadening prohibited quantitative evaluation of Δ R 1 for L84δ 2 , A99β , L118 δ 2 , and L121 δ 1 (asterisks). The crystal structure of L99A at 8 atm of xenon pressure possesses three xenon molecules in cavity 4. We added two xenon molecules in cavity 3 and energy minimized. Δ R 1 were estimated from 1/r 6 weighted distance dependence from each xenon site, using equation (2 of the protein seems to be about 5-10 times greater than the predicted binding probability according to our estimation of K d (i.e., 1.3% at atmospheric pressure). This difference indicates that the effective distance between the O 2 and the protein proton is larger than the van der Waals contact. We speculate that O 2 does not tightly contact the protein proton and is allowed to move in the cavity space. We discuss the effective distance and dynamics of O 2 in hydrophobic cavities in the section: Rotational and translational diffusion of O 2 in hydrophobic cavities. Small, but substantial, Δ R 1 values were observed at the residues in the N-terminal domain (i.e., residues 1-70), corresponding to the parameter f, and originates from O 2 diffusing around the protein surface and proton spin diffusion, as discussed previously 38 .
The paramagnetic effect leads to line broadening. We investigated the line-widths (Δ ν 1/2 ) of resonance lines in the 15 N and 1 H dimensions of the refocused-HSQC as a function of O 2 concentration. Because the refocused-HSQC 39 provides line widths that are directly proportional to the transverse relaxation rate of in-phase nitrogen coherences, R 2 N , without contributions from 1 H relaxation, it is suitable to estimate the contributions of conformational exchange and PRE contributions to the 15 N line-width 31 . Figure S6 shows 1 H and 15 N line-widths for residues 85, 88, 89, 99, 100, and 118 as a function of dissolved oxygen concentration. To within experimental error, the 15 N line widths with increasing O 2 concentration are constant for residues around the hydrophobic cavities, while a strong increase in 1 H line widths for the same residues is observed. This observation supports the notion that dipolar PRE is the primary cause of line broadening, as this effect is proportional to the square of the gyromagnetic ratio of the nucleus involved. The observation that all of the amide and methyl protons showing severe line broadening and loss of signal intensities are located less than 6 Å from the closest xenon binding site (Table S1) further supports this. These results show that the PRE to the transverse relaxation rates lead to the line broadening.

Origin of O 2 -induced chemical shift change.
We showed that changes in 15 N, 13 C, and 1 H chemical shifts were specific to O 2 and not observed for diamagnetic gases N 2 and Ar. We next sought to understand whether O 2 -induced changes result mainly from the paramagnetic property of O 2 . In general, localized unpaired electrons of a paramagnet couple to the surrounding nuclei (i.e., hyperfine coupling) and may induce chemical shift changes through spin polarization and delocalization conveyed through the molecular orbitals of the molecule (contact shifts) or through the magnetic field emanating directly from the paramagnetic center (pseudocontact shifts, PCSs). In the case of O 2 , PCSs result from the anisotropic g-tensor of the unpaired electrons. PCSs depend on the distance between the paramagnet center and the nucleus of interest and the orientation with respect to the principal axes of the magnetic susceptibility tensor (i.e., χ-tensor). If PCSs were significant, the additional magnetic field would be sensed to a similar degree by the nuclei in 15 N-1 H and 13 C-1 H bonds and therefore lead to diagonal displacements of signals in 1 H/ 15 N and 1 H/ 13 C HSQC spectra, where the magnitude of change would be about the same for the bonded nuclei when measured in ppm. Such peak movement was not observed with increasing O 2 concentration. To the contrary, 15 N and 13 C chemical shift changes upon O 2 binding were much larger than 1 H shift changes. In addition, if oxygen exchanges rapidly with the protein interior without preferred bound orientation or rapidly reorients itself with respect to the protein in the bound state, PCSs would average to zero. Instead, the O 2 molecules in the hydrophobic cavities may have frequent collisions with nuclei. The collisions of O 2 might allow delocalization of unpaired electron spins to 13  Do N 2 and Ar interact with the cavities of the protein? Similarities of the properties of the gases, such as mole fraction solubility, van der Waals radius and polarizability 42,43 , imply that N 2 and Ar could associate to the hydrophobic cavities of L99A (Further discussions are in Supplementary Information). Indeed, binding of Ar at cavity 4 was observed in crystal structures of L99A at 8 to 32 bar of Ar pressure 10 . These results indicate that chemical shift changes by binding of noble gases, which can be attributed to changes in structure and conformational equilibrium, are in general much smaller than those by paramagnetic shifts of O 2 . Based on these observations and considerations, we conclude that O 2 -induced changes in chemical shifts resulted primarily from changes in contact shifts rather than due to the PCSs of O 2 (see following section) and structural changes and conformational equilibria of the protein.

Rotational and translational diffusion of O 2 in hydrophobic cavities.
In order to investigate rotational and translational diffusion of O 2 in cavities 3 and 4, molecular dynamics (MD) simulations of 100 nanoseconds were performed five times. An O 2 molecule was inserted in both cavities 3 and 4 of the crystal structure of L99A (PDB ID; 1c6k) from which the pre-existing three xenon molecules were removed. One of the MD simulations is shown in Supplementary Movie S1 online. Note that the movie consists of 500 snapshots taken every 0.2 nanoseconds. O 2 frequently moves around each hydrophobic cavity and rotates many times within 100 ns. Similar results were obtained in the four other MD simulations. Figure 7a shows xenon binding sites 1-5 in L99A, and Fig. 7b shows a density map of O 2 molecules in cavities 3 and 4, obtained by the 100 nanoseconds MD simulation. While O 2 samples almost all spaces in cavity 4, sampling frequencies that were more than 4 times higher than the average were observed only at sites 1 and 3 (Fig. 7b). These results suggest that O 2 density is substantially higher at xenon binding sites 1 and 3 than at site 2, which is qualitatively consistent with the prediction by Δ R 1 for amide and methyl protons. In the small hydrophobic cavity 3, O 2 seems to mostly populate site 5. These results are consistent with the previous one nanosecond MD simulation by Mann and Hermans 44 .
Interestingly  Supplementary Fig. S7 online. Although it is known that L99A allows xenon and benzene to bind to cavity 4, the pathway of ligand access and egress is unknown. The present results provide the first insights in the potential pathway of ligand binding and unbinding to cavity 4, as well as egress from cavity 3.
We separately performed additional one nanosecond MD simulations to understand more details of the rotational diffusion of O 2 in cavity 4. A time dependent relaxation of the rotational correlation function was more reasonably fitted to a bi-exponential function than a single-exponential one, as shown in Supplementary  Fig. S8 online. Rotational correlation times estimated by a bi-exponential function were 0.164 ± 0.006 ps and 1.41 ± 0.02 ps. Such a bimodal rotational correlation is well known for small molecules in many solvents 45,46 . The fast and slow components would be related to the inertial characteristics of the small molecule and diffusive solvent motions, respectively 45,46 . Therefore, we speculate that the fast and slow components in the present case are related to the inertial characteristic of O 2 and protein internal dynamics taking place around the cavity, respectively. A different explanation is anisotropy of rotational motion of O 2 in the cavity. In any case, it appears that the rotational motions taking place at sub-nanosecond will average the orientation with respect to the principal axes of the magnetic susceptibility tensor and reduce the PCSs. On the other hand, chemical shifts can be changed with an increase in a population of "the O 2 -bound state". Because the observed NMR signals are the ensemble average of exchanging conformations, O 2 -induced Δ R 1 of a particular proton would depend on the averaged-distance from the O 2 molecule, which would be larger than the van der Waals contact. These results strongly support our discussions on PREs and PCSs. O 2 associates selectively and preferentially to hydrophobic and not to hydrophilic internal protein cavities. Although the size of hydrophilic cavities (cavities 1 and 2) is considered to be enough for O 2 binding, we could not detect O 2 -induced changes in spectral parameters, such as chemical shifts, peak line-widths, and longitudinal relaxation rate constants, for residues around these cavities. Interestingly, X-ray crystallography found electron density in the hydrophobic and hydrophilic cavities (cavities 1, 2, and 4), which are attributed to water molecules. In hydrophilic cavity 1, two well-ordered water molecules were seen, while a single well-ordered water molecule was seen in hydrophilic cavity 2. In hydrophobic cavity 4, weak electron density was distributed around the cavity 47 . In contrast, no gas molecules or water molecules have previously been detected in cavity 3, as far as we are aware. O 2 binding to cavities is generally considered to be in competition with water binding. Water molecules associate to the hydrophilic cavities more than O 2 probably due to its dipolar property and higher concentration (i.e., ~55.6 M) in solution. Therefore, the observation of O 2 penetration into the hydrophilic cavities of proteins is expected to be difficult. The present data show that O 2 penetrates into the protein interior and selectively and preferentially associates with the two hydrophobic cavities of the protein. The preferential partitioning of O 2 into hydrophobic regions is consistent with earlier observation of O 2 binding to ribonuclease A 12 and lipid bilayer 18 , for example. O 2 association with the hydrophobic cavities and dynamic motion of the protein. Penetration of O 2 into the protein interior requires spaces greater than its molecular size, such as a cavity, and transient conformational fluctuations, which provide pathways for penetration. Nucleus-electron dipolar interactions may therefore be modulated by conformational fluctuations of the protein and depend on the solubility of O 2 in the protein interior. Our previous work by high-pressure NMR spectroscopy revealed that T4 lysozyme L99A has at least two types of conformational fluctuations 6 ; one takes place within the ground state ensemble, which is limited to the C-terminal domain. These fluctuations occur more rapidly than a millisecond and provide heterogeneous conformations around the hydrophobic cavities. The conformational fluctuations within the ground state ensemble may allow a penetration of gas molecules into the protein interior. Indeed, MD simulations showed that O 2 molecules frequently move around cavities 3 and 4 and may go back and forth between the cavities and the outside of the protein during 100 ns. A second motion present for T4 lysozyme L99A is a conformational fluctuation between the ground state and a transiently formed high-energy "excited" state of the protein, which takes place on the millisecond time scale (average 0.7 ms) 30 . Because the aromatic side chain of F114 flips into the enlarged cavity in the excited state, O 2 association will compete with the F114 flip-in in this excited state. Finally, O 2 -induced changes in chemical shifts and signal intensities were also observed for some methyl groups around cavity 4 (39 Å 3 ) of cysteine-free wild-type (C54T/C97A; WT*) T4 lysozyme (see Supplementary Fig. S9 online). These observations indicate that O 2 molecules associated to the cavity 4 of WT* in the absence of a conformational fluctuation between ground and transiently formed excited states. Taken together, these data suggest that O 2 association is facilitated by the conformational fluctuation taking place within the ground state ensemble, rather than between the ground and excited states of L99A, and that the transiently formed excited state is not involved in gas binding to the enlarged cavity. The fact that quenching of the tryptophan fluorescence by O 2 occurs on nanosecond time scale for several native proteins 16,17 agrees with this conclusion.

Conclusion
Gas-pressure NMR spectroscopy using O 2 has been used to explore dynamic cavities in T4 lysozyme L99A. We have come to the following conclusions: (1) O 2 preferentially interacts with hydrophobic cavities and induces significant changes in NMR spectra, such as increased peak-widths and longitudinal relaxation rate constants, due to its paramagnetic property. (2) O 2 associates to the two hydrophobic cavities in T4 lysozyme L99A. So far, no gas or water molecules have been detected in cavity 3. The present study provides the first evidence of ligand binding to cavity 3. (3) O 2 -induced relaxation enhancements could be adequately accounted for by the paramagnetic dipolar relaxation, assuming 1/r 6 -weighted contributions from five sites, where r is the distance to the paramagnet. (4) The dissociation constant for O 2 binding to cavity 4 of the protein is 21 mM, indicating that about 1% of the protein contains O 2 molecules in the dynamic hydrophobic cavity at ambient pressure. (5) According to MD simulations, O 2 molecules in the hydrophobic cavities of the protein frequently move and rotate on the picosecond to nanosecond time scale. The cleft between helices D and G and the channel consisting of helices D, E, G, H, and J could be the potential gateway for ligand binding to cavity 4. O 2 association with the hydrophobic cavities would be facilitated by the conformational fluctuations taking place within the ground state ensemble, rather than between the conformational ground and excited states of the protein. (6) The rotational and translational motions of O 2 in the hydrophobic cavities may effectively reduce potential pseudocontact shift contributions to nuclear shielding.
The combination of NMR and MD simulation provides static and dynamic aspects of O 2 binding to hydrophobic cavities. This approach might also be useful to probe the permeation pathway of ions or small molecules, such as channel-blocking molecules in membrane proteins 48 and hydrophobic binding pockets for ligands, including drug compounds. Knowledge of protein permeation by oxygen is also highly relevant for optical spectroscopy and microscopy, where O 2 dissolved in the protein matrix leads to quenching and bleaching, and knowledge of oxygen association pockets may facilitate the elimination of oxygen-free cavities through protein engineering. This strategy has the potential to greatly improve our understanding of the role played by protein cavities in biologically relevant processes.

Methods
Sample preparation. T4 lysozyme L99A was prepared from the recombinant cysteine-free T4 lysozyme (WT*, C54T/C97A) 31 . Uniformly 15 N-labeled or 15 N/ 13 C-labeled L99A was produced in M9 media with 15 NH 4 Cl and 13 C 6 glucose as the sole nitrogen and carbon sources, following established protocols 31 . The purified protein sample was dialyzed in 50 mM phosphate buffer including 25 mM NaCl at pH 5.5. Sample concentration was measured by UV absorption at 280 nm and was calculated with a molar extinction coefficient of 25440 M −1 cm −1 at 280 nm. NMR experiments. We used 1 H 500 MHz (Bruker BioSpin Co. AVANCE ΙΙΙ) or 600 MHz (Bruker BioSpin Co. AVANCE) NMR spectrometers. In order to study the binding of oxygen (O 2 ), nitrogen (N 2 ), and argon (Ar) to the protein, we used a pressure resistance NMR tube (528-QPV-7, Wilmad-Lab Glass Co.) connected to a gas cylinder by PTFE tubing. Gas pressure was applied to adjust their concentrations in the protein solution. Mole fraction solubility of O 2 , N 2 , and Ar in water are 2.3 × 10 −5 , 1.2 × 10 −5 , and 2.5 × 10 −5 , respectively, at 298 K 43 . In this article, we use absolute pressure (gauge pressure + atmospheric pressure). 1 H-NMR, 1 H/ 15 N refocused-HSQC 39 , and 1 H/ 13 C CT-HSQC spectra were obtained for 0.50 mM uniformly 15 N-labeled or 1.0 mM 15 N/ 13 C-labeled T4 lysozyme L99A solution at 298 K at different gas pressures. 1 H longitudinal relaxation enhancements for amide and methyl protons were obtained from 1 H/ 15 N HSQC and 1 H/ 13 C CT HSQC spectra using saturation-recovery, achieved with proton x and y purge pulses followed by a relaxation delay before each scan 49 . Seven to ten relaxation delays ranging from 0.003 s to 1.5 s were used. Spectral analysis was performed using NMRPipe 50 and Sparky 51 . Molecular dynamics simulation. Molecular dynamics (MD) simulations of 100 nanoseconds were performed five times using GROMACS 4.6.4 simulator 52 . The system contained a T4 lysozyme L99A (PDB ID; 1c6k), two O 2 molecules, 8 chloride ions, and about 15,000 water molecules. Three xenon molecules in cavity 4 of the protein were removed from the structure, and one O 2 molecule was inserted in each hydrophobic cavity (i.e., cavities 3 and 4). The OPLSLL force field 53 was used for the protein, and the TIP4P model was used for water 54 . Potential parameters for O 2 and chloride ions were as described in the literature 55,56 . MD simulations were conducted with the NPT ensemble (300 K, 1 bar) in a truncated dodecahedron box with dimensions of 25.8 Å. Temperature was controlled using a Langevin thermostat with a viscosity of 0.5 ps −1 . Pressure was controlled by a Parrinello− Rahman barostat with relaxation times of 2.0 ps 57 . Electrostatics were treated using the particle mesh Ewald (PME) method with a 10.0 Å cutoff distance 58 . The van der Waals interactions were expressed using the twin-range cutoff method with 10.0 and 12.0 Å distances. Covalent bonds for hydrogen atoms in the polypeptide were constrained using the linear constraint solver (LINCS) 59 . Covalent bonds in water were constrained using the SETTLE algorithm 60 . The integration time step was 2 femtoseconds.
In order to estimate the rotational correlation times of O 2 in cavity 4, we performed separate one nanosecond MD simulations. To avoid artifacts from the thermostat and barostat, we carried out the MD simulations in the NVE ensemble with the structure after 100 ns simulation in the NPT ensemble 61 . Snapshots were recorded every 0.01 picoseconds. Rotational correlation times were calculated using 1 ns trajectories of the O 2 molecule in cavity 4. We defined a direction vector between the two oxygen atoms relative to the orientation of the protein.