Evolutionary adaptation from hydrolytic to oxygenolytic catalysis at the α/β-hydrolase fold

Protein fold adaptation to novel enzymatic reactions is a fundamental evolutionary process. Cofactor-independent oxygenases degrading N-heteroaromatic substrates belong to the α/β-hydrolase (ABH) fold superfamily that typically does not catalyze oxygenation reactions. Here, we have integrated crystallographic analyses under normoxic and hyperoxic conditions with molecular dynamics and quantum mechanical calculations to investigate its prototypic 1-H-3-hydroxy-4-oxoquinaldine 2,4-dioxygenase (HOD) member. O2 localization to the “oxyanion hole”, where catalysis occurs, is an unfavorable event and the direct competition between dioxygen and water for this site is modulated by the “nucleophilic elbow” residue. A hydrophobic pocket that overlaps with the organic substrate binding site can act as a proximal dioxygen reservoir. Freeze-trap pressurization allowed the structure of the ternary complex with a substrate analogue and O2 bound at the oxyanion hole to be determined. Theoretical calculations reveal that O2 orientation is coupled to the charge of the bound organic ligand. When 1-H-3-hydroxy-4-oxoquinaldine is uncharged, O2 binds with its molecular axis along the ligand's C2–C4 direction in full agreement with the crystal structure. Substrate activation triggered by deprotonation of its 3-OH group by the His-Asp dyad, rotates O2 by approximately 60°. This geometry maximizes the charge transfer between the substrate and O2, thus weakening the double bond of the latter. Electron density transfer to the O2(π*) orbital promotes the formation of the peroxide intermediate via intersystem crossing that is rate-determining. Our work provides a detailed picture of how evolution has repurposed the ABH-fold architecture and its simple catalytic machinery to accomplish metal-independent oxygenation.


DNA techniques, protein expression, and protein purification
For protein overexpression in E. coli, the target gene was inserted into the pQE30 vector (Qiagen) at the BamHI and SalI restriction sites.An uncleavable N-terminal His6 fusion tag of sequence MRGSHHHHHHGS was added to the gene product to facilitate protein purification.
As wild-type HOD tends to dimerize due to oxidative formation of an intermolecular disulfide bridge, we took advantage of the C69S substitution that abrogates this problem. 1 The HOD C69S variant (HOD C69S ) has catalytic properties that are identical to those of wild-type HOD.
All work described here was carried out in this background.For simplicity, throughout the text, we refer to His6-HOD C69S simply as HOD.Site-directed mutagenesis to generate all HOD variants was performed using standard techniques followed by gene sequencing for validation.

X-ray crystallography
Crystals of HOD, HOD H251A , HOD S101A exhibiting variable degree of twinning and belonging to space group P212121 or P41 were obtained by the vapor-diffusion method at 18 °C by mixing equal volumes of protein solution (150 mg/ml in 20 mM Tris-HCl,100 mM NaCl, 2 mM EDTA, 1mM DTT, pH 7.5) and a reservoir solution containing 1.65 M Na/K tartrate and 100 mM Hepes, pH 7.0.Complexes with the organic compounds were obtained by soaking crystals in a reservoir solution enriched by a 2 mM solution of the appropriate ligand for typically around six hours.Noble gas derivatization was achieved by pressurizing at 30 bar a capillarymounted HOD crystal in pure Xe atmosphere.Pressure was maintained for the entire duration of the room-temperature diffraction experiment.For dioxygen pressurization experiments, we employed the methodology dubbed 'soak-and-freeze'. 2 Crystals were mounted in specific sample high pressure-supports (SPINE standard pins with a capillary) and exposed to 40 bar O2 for two minutes at ambient temperature prior to cryo-cooling in liquid N2 whilst still under pressure.Except for the Xe pressurization measurement, all data collections were performed at 100K with crystals cryo-protected using the reservoir condition spiked with 20% glycerol (v/v).
Complete data sets were measured at beam lines I02, I04-1, I24 of Diamond Light Source (Didcot, United Kingdom) or ID29, BM30A of the European Synchrotron Radiation Facility (Grenoble, France).Data processing was carried out using the xia2 pipeline 3 and crystal structures solved either by the molecular replacement technique using Phaser 4 or by difference Fourier methods.Model refinement was achieved with the Refmac5 5 package.Model building was performed with COOT. 6A summary of data collection and refinement statistics are shown in Table S1.

In crystallo turnover
In solution, the H251A replacement renders the enzyme almost inactive (kcat of 0.0006 s -1 at pH 6.5 compared to 16.3 s -1 for wtHOD) as it lacks the functional H251-D126 dyad required for the fast substrate activation step in which its 3OH group is deprotonated. 7,8 reviously, we have obtained the X-ray structure of the HOD H251A -QND complex under normoxic conditions. 7e structure shows that under these conditions the substrate binds in the active site without any appreciable turnover.Despite the very limited activity of HOD H251A , O2 pressurization of HOD H251A -QND crystals (40 bar for 2 minutes) results in complete turnover and electron density maps at 2.0 Å show the bound N-acetylanthranilate product (Figure S8A).This experiment clearly demonstrates that O2 can reach the R-site in the preformed E-S complex promoting QND conversion into the NAA product during pressurization.Reactivity must therefore be attributed to a small quantity of the QND anion present at the equilibrium.In crystallo UV/Vis spectra are consistent product formation without substantial accumulation of the substrate anion (Figure S8B).The structure suggests that proton abstraction is mediated by solvent molecules that connect the 3OH group of the substrate with D126 (Figure S8A).
UV/Vis spectra shown in (Figure S8B) were recorded using the in crystallo spectroscopic facility at ESRF. 9

Standard molecular dynamics (MD) simulations
Initial coordinates for substrate-free simulations were taken from the relevant crystallographic structure (PDB code 2WJ3).Hydrogen atoms were added by the antechamber module of the AMBER package using geometric rules. 10The protonation state of histidine residues and other amino acids were checked using the PROPKA tool. 11,12 he protonation state of residue H251 was chosen in relation to the charge of the substrate.The position of hydrogens on the histidine ring (delta/epsilon) were determined by visual inspection of the high-resolution structure looking at possible hydrogen bonds with nearest residues (epsilon: 7, 40, 100, 102, 154, 164, 171, 238, 264; delta: 38, 64, 172, 217, 251 when neutral).The enzyme was solvated by in a cubic box of 85´85´85 Å 3 , with water molecules placed initially 5 Å away from the protein surface whilst retaining all crystallographic water molecules.Fifteen sodium ions were added to neutralize the protein charge.The system was initially relaxed using the NAMD package 13 while production runs were performed with ACEMD. 14For the first 2 ns, restraints were applied to both Ca and Cb protein carbon atoms and to crystallographic waters.At T = 300 K in the NPT ensemble, all restraints were gradually removed over 4 ns.After 2 ns in the NPT ensemble at 300 K and P = 1 atm, we calculated the average size of the box and we used this value for simulating the system in the NVT ensemble (L=73.67Å, 39081 atoms).Long-range electrostatic interactions were treated with the Soft-Particle-Mesh-Ewald algorithm 15 (64 grid points with direct cutoff at 9.0 Å, as suggested by the ACEMD software; the same cutoff holds for van der Waals interactions).We employed the Amber99SB ILDN forcefield 16 and TIP3P 17 for the protein and water molecules, respectively.The GAFF force-field parameters 18 were used to describe the organic molecules.Partial atomic charges were evaluated according to the RESP approach: 19 the molecule was first optimized at the HF/6-31G(d) level, up to a convergence in energy of 10 -5 au, using the Gaussian software. 20Atomic RESP charges were derived from the electrostatic potential using the antechamber module of the AMBER package. 10O2 was parametrized with the LJ Amber parametrization, by adding a small charge to have a dipole of 0.121 eÅ.
Two independent 1 µs-long replica simulations were carried out without or with explicit O2 (WAT-MD and OXY-MD, respectively).For OXY-MD, ten O2 molecules were added to the system in the cubic box with L = 73.67Å. Normalized dioxygen distribution was calculated by sampling O2 positions every 50 ps onto a matrix of 0.6 ´ 0.6 ´ 0.6 Å 3 elementary cubes and dividing the total atom count in each volume by the total number of frames.Initial coordinates for substrate-bound simulations were taken from the anaerobic HOD-QND complex (PDB code 2WJ4) and the system was setup as described above.QND was initially considered in its protonated form at O3 (QNDH) with neutral H251 protonated on ND1.The force field for the substrate was obtained as described previously. 21The protonation of all other histidine residues was chosen to preserve the H-bond pattern observed in the X-ray structure.At 300 K in the NPT ensemble, we removed protein restraints over 4 ns and transferred the proton from the substrate's O3 atom to H251(NE2) as this represents the mechanistically-relevant complex. 7Weak restraints on the substrate positioning were maintained throughout to prevent its escape from the active site during the runs.To study O2 entry in the HOD-QND complex we performed a total of 80 independent 200-ns MD runs by adding to the system either ten O2 molecules in the L = 73.67Å box (OXY10-S-MD runs) or one O2 molecule confined in a sphere with a 17-Å radius of centered on S101 (OXY1-S-MD runs).Of the 40 OXY10-S-MD runs, 25 were carried out with a harmonic restrains on S101 to maintain its trans rotamer and 15 had no restraints on S101.Of the 40 OXY1-S-MD runs, 20 were carried out with a harmonic restrains on S101 to maintain its trans rotamer and 20 had no restraints on S101.

Metadynamics Simulations
To investigate the oxygen diffusion within HOD, we also performed metadynamics simulations employing a bias on the O2 position. 22As coordinates to be biased we selected the position of the dioxygen center of mass in the 3D Cartesian space.As consequence of the use of absolute coordinates, we kept the protein frame in place by restraining the position of a set of Ca atoms (residues 31, 101, 112, 243) to the reference structure obtained after relaxation.We performed well-tempered metadynamics, 23 DT= 3000 K, with bias added as Gaussians with a frequency of 1 ps.As the R-site cavity is rather small, we decided to work with narrow Gaussians (width=0.3Å).For this reason, we started with a high value for the initial height (1 kcal/mol) and this value is smoothly decreased by the well-tempered algorithm during long runs.
Dioxygen was confined to a sphere of 17-Å radius centered on S101.We performed three simulations of 1.8 µs each.In two replica runs (identical conditions except for initial velocities chosen randomly) S101 was restrained to its trans rotamer.No restraints were employed in the third one.We calculated the error on metadynamics by reweighting the positions sampled during the runs with the corresponding bias term and comparing the free energy surface with the one obtained integrating the bias potential.

Analysis of O2 and H2O escape times from the R-site
Preliminary tests indicated that O2 residence times at the R-site are in the order of ns, corresponding to an energy barrier of few kcal/mol, achievable with thermal motions (kT = 0.6 kcal/mol at 300 K).Thus, the use of enhanced sampling techniques such as metadynamics 24 or TPS-derived methods 25 was not required for this analysis.As our starting point we took advantage of a protein state found in our OXY1-S-MD run in which O2 was found at the Rsite.This was subjected to 200 steps of conjugate gradient minimization and then thermalised for 1 ns at 300 K before assigning new random velocities to 50 independent standard MD simulations.For each trajectory we evaluated the time spent by O2 at the R-site before it abandons this location.The escape process can be seen as a diffusion process over a barrier dividing two basins.Although there are uncertainties in the estimation of the precise time when a molecule goes over the unknown transition state for calculating the residence time, 25 we can associate to the transition over a barrier the abrupt change in a collective variable/parameter.By plotting the distance of oxygen from the substrate geometric center (calculated using heavy atoms only), we recognized that around 6 Å there is an abrupt change of position.Thus, we defined the residence time as the time to reach in the sampled trajectories that distance, without considering the possibility to recross back the transition state or considering the transmission coefficient equal to one.Harmonic restraints were employed to preserve S101 rotamers and the substrate position (on average) during the simulations.
Escape times (koff -1 ) are expected to follow the exponential distribution: belonging to the family of Poisson statistics.From this distribution we can derive the probability that an event occurs before time t: Numerically we can fit F(t) with the cumulative histogram of escape times to obtain the single parameter lambda that characterizes this distribution.We can also determine the average escape time using block averages.Fitting the cumulative distribution is however more appropriate because it is difficult to sample long escape times.The cumulative distributions were obtained by accumulating data with a bin width of 1 ns.

QM calculations of the isolated ligand-O2 complexes
DFT calculations to determine the equilibrium geometries of the ligand-O2 complexes were carried out at the M062X/6-31+G(d,p) level of theory using QChem 5.2 26 using the D3(0) dispersion correction. 27Since even for simple systems, such as toluene-O2 there are many conformers possible, 28 we have reduced the dimensionality of the problem by restraining the sum of the (C2-O1) and (C4-O2) ligand-O2 distances to their crystallographic values.We find however that around the equilibrium geometry, results are quite insensitive to small changes in the restrain value (Figure S10).This restraint was not applied in the QM/MM calculations, as the protein scaffold prevents O2 motion away from its pocket.Energy decomposition analysis (EDA) 29 as implemented in QChem 5.2 was employed to calculate the interaction energy between O2 and QND/MQO, and its contributing terms.Delocalization indices 30 (d) using the partitioning Mulliken scheme 31 were calculated using the NDELOC code, 32 interfaced with Gaussian16. 20.They are a measure of 'electron sharing' between sets of atomic centers, and their 2-order magnitudes allowed us to assign "bond orders".To avoid a wrong partition within the Hilbert space, diffuse functions were not included in d values calculations.

QM/MM energy minimizations
The HOD-MQO/QND-O2 ternary complexes used as starting point for the QM/MM calculations were prepared from the relevant crystallographic structures using the CHARMM software 33 and the CHARMM36 forcefield. 34MD simulations were performed with NAMD 2.13. 13The coordinates of the hydrogen atoms were generated using the standard protonation states for the titrable residues.H251 was protonated in the HOD-QND system, as this residue abstracts the 3OH proton from QND.The protein complexes were placed in the center of a cubic TIP3 water box, large enough to include the protein and at least 10 Å of solvent on all sides.KCl at a concentration of 0.15 M was added to neutralize the system, and the Particle Mesh Ewald method 35 was used to treat the electrostatics of the periodic boundary conditions.
The system was subjected to a classical MD equilibration of 10 ns at constant temperature (303.15K) and pressure (1 atm), where to avoid O2 diffusion outside the active site, the position of QND/MQO and O2 were restrained throughout the equilibration.After equilibration, the system was trimmed to a sphere of 26 Å centered around O2, where the atoms further away than 18 Å were kept frozen during all QM/MM calculations.Full electrostatic embedding was adopted, and hydrogen link atoms were used to treat the QM/MM boundaries, when necessary.
QM calculations were carried out at the M062X/6-31+G(d,p) level of theory using Grimme's dispersion correction as in the gas-phase calculations.In a first set of calculations, the reactants were optimized and, subsequently, the (C4-C2-O2-O1) dihedral angle (F) was scanned through geometry restrained minimizations.For this study, the QM region included the substrate, O2, the sidechain of H251 (protonated for QND-O2 simulation), H102, G103, the backbone of G35, W36, C37, and three water molecules close to either O2 or the histidine sidechain.No restraints were applied the (C2-O1) and (C4-O2) distances.In a second set of calculations, we calculated the energy profile along the reaction path.For the CO-release stage, we used as a reaction were carried out using two different QM regions, one including just the substrate (QND or MQO) and O2, and another also including the sidechain of H100, S101, and H251 as well as a water molecule nearby.The results were found to be relative insensitive to the size of the QM region.To describe the intersystem crossing (ISC) stage, the energy difference between the singlet and triplet states was chosen as a reaction coordinate, 36 so the energy of the Minimum Energy Crossing Point (MECP) is that obtained when the energy difference is zero.The results are shown in (Figure S13).For this system, the energy difference between the triplet and singlet states was strongly coupled with the approaching of O2 to QND, so it was also possible to represent the reaction pathway as a function of the composite reaction coordinate α.This is done in Figure 6 of the main text, where all the geometries for the intersystem crossing stage appear for α < -0.4 Å.For the ISC stage we only considered the small QM region, and for the sake of simplicity only results with the small QM region are shown in Figure 6.

Calculation of the energy barrier for the intersystem crossing (ISC) step
To compare the effective energy barriers between for two reactions, one spin-allowed (adiabatic) and another spin-forbidden, the barrier height of the latter should be corrected to account for the limited swapping probability between the electronic states involved.According to the Non-Adiabatic Transition State Theory, the rate coefficient of a spin-forbidden transition should be calculated as: 37 where kNon-Adiabat is the actual rate coefficient of the process, kAdiabat would be the rate coefficient of a hypothetical analogous spin-allowed reaction (just one electronic state), and ápshñ is the probability of spin-flipping, in this case the probability of hopping from the triplet to the singlet state at the MECP.ápshñ can be calculated from the Landau-Zener theory. 38,39 sing the double passage version of the spin-diabatic Landau-Zener formula for a given energy is, where  -. is given by: where E is the energy, EMECP is the energy of the MECP, µ is the reduced mass of the normal mode orthogonal to the crossing seam, | DF | is the norm of the difference between the gradients of the two Potential Energy Surfaces at the crossing point, and SOC is the effective spin-orbit coupling calculated at the MECP.For E < EMECP,  -. = 0.
SOC was calculated at the CASPT2 (22,38)/cc-pvdz level of theory using Open-MOLCAS, 40 with the active space formed by all the valence orbitals of O1, and O2 and those related to the C2-O2 and O1-O2 bonds.The value obtained (SOC= 22 cm -1 ) was found to be barely sensitive to small changes in the active space.We obtained ápshñ = 0.005 which leads to an increase of 3.1 kcal/mol to the barrier height.
According to our calculations the energy of the MECP was 10.2 kcal/mol above the energy of the triplet QND -O2 complex.At the MECP, rC2-O2 =1.62 Å, and the QND moiety is no longer planar.The geometry at the MECP is similar to that obtained by Minaev 41 at the multireference CASSCF (16,11) level of theory.In Ref. 42 , Silva computed the MECP for four different (1H)-3-hydroxy-4-oxoquinolines-O2 complexes, and found considerably different rC2-O2 depending on the substituent, ranging from 2.23 Å for QND to 1.56 Å when the methyl group is replaced by a fluoride.The energy of the MECP for the F substituted oxoquinoline was the smallest, 9.2 kcal/mol, 1 kcal/mol lower than in our calculations, also associated with a small rC2-O2 distance.Differences between our results and those in Ref. 42 can be associated to the different functional used (B3LYP vs M06-2X) or different simulation of the interactions with the solvent (PCM continuum model using water as a solvent vs QM/MM model).To check the effect of rC2-O2 on the energy of the MECP, we calculated the minimum of the crossing seam for different given rC2-O2 distances and we obtained that for rC2-O2 = 2.0 and 2.3 Å, the energies of the crossing point were only 1.8 kcal/mol and 2.7 kcal/mol above the MECP, respectively.
With respect to the SOC value, a value of 75 cm -1 was obtained for other similar spinforbidden cofactor-less enzymatic reactions, 36 leading to an increase of just 1.6 kcal/mol to the barrier, 1.5 kcal/mol lower than for the reaction between QND and O2.Assuming that this value is the upper-limit for a cofactor-less reaction involving O2, other points of the crossing seam (and not necessarily the MECP) could act as effective transition state if their energy is no more than 1.5 kcal/mol above the energy of the MECP and show a SOC value closer to 75 cm -1 .
Even if that is the case, the effective barrier of the intersystem-crossing stage would be higher than for the CO release stage.
áI/s(I)ñ    Whilst the His-Asp dyad is totally conserved, the nucleophile is not and, in several organisms, the latter residues is replaced by an alanine.The map displayed in magenta is the NCS-averaged anomalous difference map contoured at +5.5σ.We note that the anomalous signal for this experiment is weak as the wavelength used for data collection was far from the f" peak to minimize radiation damage at room temperature.Also, data at 2.9 Å-resolution exhibit non-merohedral twinning (twinning fractions 0.52/0.48)which negatively affects the anomalous signal.Despite these limitations the interpretation of this site as Xe is unambiguous.Hydrophobic residues stabilizing the Xe atom are highlighted in gold and shown as sticks.The catalytic triad (S101/H251/D126) is also shown.).In the case of (E) fitting is not possible as there are too few escape events with 50 ns.Pressurization enables O2 access to the R-site in the pre-formed HOD H251A -QND complex promoting turnover without buildup of the substrate anion.The structure suggests that the marginal residual catalytic activity of the HOD H251A variant (kcat = 0.0006 s -1 at pH 6.5 compared to 16.3 s -1 for wtHOD) relies on solvent-mediated substrate deprotonation.UV/Vis spectra were recorded using the online microspec at the BM30A beamline at the ESRF.Around the equilibrium geometries, the energy is relatively insensitive to small changes in the restraint value.The two minima for the MQO-O2 complex at Φ=0 and 180º represents the same structure, reflecting that O1 and O2 are equivalent.Figure S15.Shape of the half-occupied orbitals for the QND-O2 complex at the reactants asymptote (left panels), and at the MECP (right panels).The orbitals were localized following the procedure described in Ref. 43 At the reactants asymptote, half-occupied orbitals correspond to the two p* orbitals of O2.The p* perpendicular to the p cloud of QND (px*), is not mixed with any QND orbital.However, the p* orbital that stacks with the p cloud of QND (py*) is significantly mixed with the p cloud of QND.As QND and O2 approach each other, mixing is more substantial, leading to an occupied orbital with bonding character along the C2-O bond, and a higher energy orbital that shows an antibonding character both along O1-O2 and C2-O bonds.The latter orbital is singly occupied in the triplet state but it is a virtual orbital in the singlet state, explaining why the singlet state is stabilized along the reaction path while triplet state is destabilized.
coordinate α = d(O1-O2) + d(C3-C4) + d(C2-C3) -d(O1-C4) -d(O2-C2), which includes all bond distances that are formed and broken during the reaction.Calculations with and without the additional restraint of d(C3-C4) = d(C2-C3) were run, showing that including this restraint aids to reach the product state, without modifying the barrier heights.To ensure convergence, geometry restrained minimizations were carried out forwards and backwards.Calculations

Figure S1 .
Figure S1.ABH fold.(A) General ABH fold topology and (B) ABH three-dimensional architecture exemplified by the crystal structure of Arthrobacter nitroguajacolicus Rü61a 1-H-3-hydroxy-4oxoquinaldine 2,4-dioxygenase (HOD) (PDB code 2WJ3).HOD is the focus of the present study.The core domain of ABHs (here shown in light grey) is an eight-stranded, mostly parallel, β-sheet surrounded by α-helices.This is very often embellished by a cap or lid domain (here shown in dark grey) inserted between β6 and the short αD that provides structural and functional variation.In the case of HOD the cap domain is constituted by four α-helices.A very conserved feature of ABHs is the Ser (nucleophile), His, and Asp (acidic residue, sometimes replaced by Asn) catalytic triad represented by black dots in the topology diagram.S101, H251 and D126 form the catalytic triad in HOD.These residues are shown as sticks in panel B with nitrogen and oxygen atoms in blue and red, respectively.The nucleophile is hosted by the 'nucleophilic elbow', a tight turn connecting β5 and αC.In some ABHfold dioxygenases the nucleophile is replaced by an alanine.The acidic residue, generally positioned after β7, is relocated after β6 in HOD.The catalytic histidine is near the C-terminus after the β8 strand.N and C indicate the protein N-terminus and C-terminus, respectively.

Figure S2 .
Figure S2.Sequence alignment of various predicted and validated bacterial ABH-fold cofactor-independent 2,4-dioxygenases exhibiting high similarity to A. nitroguajacolicus Rü61a 1-H-3-hydroxy-4-oxoquinaldine 2,4-dioxygenase (HOD).Residues are highlighted in different shades of green according to their sequence identity.Secondary structure elements for HOD are also shown for reference.Residues of the nucleophile-histidine-acidic catalytic triad are indicated.Whilst the His-Asp dyad is totally conserved, the nucleophile is not and, in several organisms, the latter residues is replaced by an alanine.

Figure S3 .
Figure S3.Root mean square deviation (r.m.s.d.) from the starting equilibrated crystallographic structure for HOD during 1 µs-long replicas Molecular Dynamics (MD) simulations.OXY-MD runs (A) were carried out in the presence of ten O2 molecules in a solvated box of approximately 74´74´74 Å 3 .WAT-MD did not include O2 molecules in the simulations.Sampling was performed every 50 ps of each run.Thick lines represent the moving average calculated using a 5 ns window.

Figure S4 .
Figure S4.Root mean square fluctuations (r.m.s.f.) from the average structure for HOD during WAT-MD (A) and OXY-MD (B) simulations.OXY-MD runs were carried out in the presence of ten O2 molecules to a solvated box of approximately 74´74´74 Å 3 .WAT-MD did not include O2 molecules in the simulations.Two replica 1-µs simulations were run for each system and data are presented as average r.m.s.f. with standard deviations.(C) is an overlay of A and B with asterisks highlighting the regions with different mobility.

Figure S5 .
Figure S5.Validation of Xe binding at the B-site.(A) Interpretation of the Xe site as a water molecule results in significant residual density in Fourier difference maps.Post-refinement 2mFo-DFc and mFo-DFc electron density maps are shown at the +1.0σ (blue) and +3.0σ (green) levels, respectively.(B)Interpretation as a Xe is satisfactory and occupancies refine between 0.4 and 0.6 in the four different molecules in the a.u.The 2mFo-DFc electron density map is shown at the +1.0σ level in blue.The map displayed in magenta is the NCS-averaged anomalous difference map contoured at +5.5σ.We note that the anomalous signal for this experiment is weak as the wavelength used for data collection was far from the f" peak to minimize radiation damage at room temperature.Also, data at 2.9 Å-resolution exhibit non-merohedral twinning (twinning fractions 0.52/0.48)which negatively affects the anomalous signal.Despite these limitations the interpretation of this site as Xe is unambiguous.Hydrophobic residues stabilizing the Xe atom are highlighted in gold and shown as sticks.The catalytic triad (S101/H251/D126) is also shown.

Figure S6 .
Figure S6.Molecular Dynamics (MD) simulation of the HOD-QND complex.(A) Root mean square deviation (r.m.s.d.) from the starting equilibrated crystallographic structure for the HOD-QND complex (PDB code 2WJ4), during a 1 µs-long MD simulation.Sampling was performed every 50 ps of each run.The thick yellow line lines represent the moving average calculated using a 5 ns window; (B) root mean square fluctuation (r.m.s.f.) from the average structure during the simulation.

Figure S7 .
Figure S7.Fitting of residence times for O2 (A-C) or water (D-F) at the R-site with S101 restrained to either the trans rotamer (A,D), or to the plus rotamer (B,E), or, alternatively, with S101 is replaced by an alanine (C,F).The cumulative distributions were fitted to the equation (2) to determine the residence times ( = koff -1 ).In the case of (E) fitting is not possible as there are too few escape events with 50 ns.

Figure S8 .
Figure S8.In crystallo turnover.(A) Active site of the HOD H251A -QND complex following O2 pressurization at 40 bar for 2 minutes.QND is fully converted to N-acetylanthranilate (NAA), the product of the HOD-catalyzed reaction.Electron density maps (2mFo-DFc) are shown at the +1.0σ level for NAA, the side chains of S101, D126, A251 (catalytic triad), the water molecule at the R-site (W1) and the solvent network (highlighted by the red dotted line) that connects D126 to the organic ligand.(B) in crystallo UV-visible spectra of HOD H251A -QND under normoxic conditions (thick orange line), after 1-minute O2 pressurization (medium thickness orange line), after 2-minutes O2 pressurization (thon orange line).The in crystallo UV-visible spectrum of the anaerobic HOD-QND complex (dotted green line) is also shown for reference.A small offset was added to the green trace to improve clarity.Pressurization enables O2 access to the R-site in the pre-formed HOD H251A -QND complex promoting turnover without buildup of the substrate anion.The structure suggests that the marginal residual catalytic activity of the HOD H251A variant (kcat = 0.0006 s -1 at pH 6.5 compared to 16.3 s -1 for wtHOD) relies on solvent-mediated substrate deprotonation.UV/Vis spectra were recorded using the online microspec at the BM30A beamline at the ESRF.

Figure S9 .
Figure S9.Representative 2mFo-DFc electron density for the HOD S101A -MQO complex under normoxic conditions.Electron density maps are shown at the +1.0σ level as chicken-wire representation for MQO (blue) and a water molecule (W1) at the R-site.The latter is displayed in red for clarity.The view is the same as in Figure 5B of the main text.

Figure S10 .
Figure S10.3D contour maps for the energy profiles for the O2-ligand (MQO on the left and QND on the right) complexes as a function of the dihedral angle Φ = (C4-C2-O2-O1), and the sum of (C2-O1) and (C4-O2) distances.The latter was used as a restraint in the calculations for the isolated complexes.Around the equilibrium geometries, the energy is relatively insensitive to small changes in the restraint value.The two minima for the MQO-O2 complex at Φ=0 and 180º represents the same structure, reflecting that O1 and O2 are equivalent.

Figure S12 .
Figure S12.QM geometry restrained optimization and delocalization indices for QNDH-O2 (A-E), MQO-O2 (F-J) and QND-O2 (K-O) complexes as a function of the dihedral angle F = (C4-C2-O2-O1).Panels (F-O) are identical to panels (C-L) in Figure 5 of the main text and are shown here for reference.In panels (D,I,N) delocalization indices for isolated O2 (dashed red line) and the aromatic C-C bond in benzene (dashed grey line) are given for reference.Protonation of QND (QNDH) results in an overall neutral system QNDH-O2 system that behaves like MQO-O2.

Figure S13 .
Figure S13.QM/MM geometry restrained optimization and delocalization indices for the MQO-O2 (A-D) and QND-O2 (E-H) complexes as a function of the dihedral angle F = (C4-C2-O2-O1).(A,E) are stick representations of the minimized structures.(B,F) show the total energy of the complexes as a function of F. (C,G) display the delocalization index (d) for the O-O bond in the complex.For comparison, the dashed red and grey lines represent the delocalization indices for O-O and C-C bonds in isolated O2 and in benzene, respectively.(D,H) show the delocalization index for the interaction between the substrate and O2.

Figure S14 .
Figure S14.Energy profile for the transition from the triplet to the singlet state and formation of the peroxide intermediate.The energy is presented as a function of the energy difference between the triplet and the singlet state, which was chosen as a reaction coordinate.The relevant structures at the points indicated are shown as stick (ball-and-stick for O2) representations.Structure 2 represents the minimum energy crossing point (MECP).

Table S1 .
Data collection and refinement statistics.