J Mol Model (2015) 21: 130 DOI 10.1007/s00894-015-2679-0

ORIGINAL PAPER

Bacterial nitric oxide reductase: a mechanism revisited by an ONIOM (DFT:MM) study Amr A. A. Attia 1 & Radu Silaghi-Dumitrescu 1

Received: 4 November 2014 / Accepted: 13 April 2015 / Published online: 29 April 2015 # Springer-Verlag Berlin Heidelberg 2015

Abstract Bacterial nitric oxide reductase (cNOR) is an important binuclear iron enzyme responsible for the reduction of nitric oxide to nitrous oxide in the catalytic cycle of bacterial respiration. The reaction mechanism of cNOR as well as the key reactive intermediates of the reaction are still under debate. Here, we report a computational study based on ONIOM (DFT:MM) calculations aimed at investigating the reaction mechanism of cNOR. The results suggest that the reaction proceeds via the mono-nitrosyl mechanism which starts off by the binding of an NO molecule to the heme b3 center, N-N hyponitrite bond formation as a result of the reaction with a second NO molecule was found to proceed with an exothermic energy barrier to yield a hyponitrite adduct forming an open (incomplete) ring conformation with the non-heme FeB center (O-N-N-O-FeB). N-O bond cleavage to yield N2O was shown to be the rate-limiting step with an activation barrier of 22.6 kcal mol-1. The dinitrosyl (trans) mechanism, previously proposed by several studies, was also examined and found unfavorable due to high activation barriers of the resulting intermediates. Keywords Bacterial nitric oxide reductase ONIOM . QM:MM . Reaction mechanism Non-heme iron

. cNOR . DFT . . Heme iron .

Electronic supplementary material The online version of this article (doi:10.1007/s00894-015-2679-0) contains supplementary material, which is available to authorized users. * Radu Silaghi-Dumitrescu [email protected] 1

Department of Chemistry, Faculty of Chemistry and Chemical Engineering, Babeș-Bolyai University, Cluj-Napoca, Romania

Introduction Bacterial nitric oxide reductase (cNOR) is a binuclear iron enzyme that catalyzes the 2-electron 2-proton reduction of nitric oxide to nitrous oxide during the anaerobic respiration of bacteria [1]. The crystal structure of cNOR was resolved recently [2] and revealed a catalytic active site that bears a binuclear iron center composed of a heme b3 iron (Feb3) and a non-heme iron (FeB) in addition to a cytochrome c subunit located in close proximity to the active site acting as the electron donor for the catalytic reaction [2]. The 2-electron 2-proton reaction taking place in cNOR is illustrated in Eq. (1). 2NO þ 2Hþ þ 2e− →N2 O þ H2 O

ð1Þ

The presence of two iron centers and the fact that two NO molecules are needed for nitrous oxide generation provide various possibilities for the binding modes of NO to the binuclear active site. Several reaction mechanisms were reported for cNOR over the course of multiple experimental and computational studies [2–13], however, the reaction mechanism of NO reduction by cNOR is still controversial and yet to be understood. Among the proposed reaction mechanisms are the mono-nitrosyl mechanism otherwise known as the cis FeB or the cis Feb3 mechanism which features the binding of the first NO to either heme b3 or the non-heme center followed by the attack of a free NO to form a hyponitrite adduct and, the trans mechanism which involves diiron dinitrosyl as the key intermediate, as illustrated in Fig. 1. The QM study by Blomberg et al. [12] and the QM/MM study by Shiraishi and co-workers [13] are among the most recent computational studies reported for the reaction mechanism of cNOR. Blomberg et al. [12] suggested a mechanism that features a cis Feb3 pathway and reported the formation of a five-membered ring (FeB-O-N-N-O) intermediate that required a total activation barrier of 21.7 kcal mol-1 for N-O

130 Page 2 of 12

J Mol Model (2015) 21: 130

Fig. 1 Proposed mechanisms for NO reduction by cNOR. (I) represents the cis Feb3 mechanism while (II) represents the trans (dinitrosyl) mechanism

bond cleavage to yield nitrous oxide. Conversely, the computational study reported by Shiraishi and coworkers [13] supported an incomplete five memberedring hyponitrous intermediate and proposed N-O bond

Fig. 2 X-ray crystal structure of cNOR (pdb id 3k9z) employed in the present QM:MM study

breaking to form N2O as the rate-limiting step with an activation barrier of 20.9 kcal mol-1 [13]. In the present report, an ONIOM (QM:MM) study is carried out to investigate the reaction mechanism of NO reduction by

J Mol Model (2015) 21: 130

Page 3 of 12 130

Models and theoretical methods

Fig. 3 The QM region employed in the QM:MM and the QM-only calculations comprising the binuclear active site with a total of 130 atoms. The α carbon of each residue lies at the border between and QM and the MM parts whereas for the QM-only calculations, the entire side chain of each residue is considered and the α carbon of each residue is fixed to the x-ray structure coordinates

cNOR and identify the possible key intermediates in the catalytic cycle thus providing additional intuition and insight toward understanding the mechanism of action of cNOR. Fig. 4 Structure (1), the optimized structure of the diferrous state. Bond lengths are in Å. Geometries of the crystal structure are shown in red

The crystal structure of cNOR (PDB id 3k9z) [10] was employed in the present QM:MM investigation and is depicted in Fig. 2. The QM region is comprised of the binuclear active site including the heme center (heme b3 and its axial His93 ligand) and the non-heme center (FeB, and its five first shell ligands Glu68, His64, His43, His29, and water) as depicted in Fig. 3. Several steps were undertaken to prepare the enzyme prior to the QM:MM calculations: hydrogen atoms were added to the crystal structure by utilizing the GROMACS suite of programs [14, 15], the PROPKA [16] toolkit was used for assigning the protonation states of the residues and calculations at pH 7 were taken as a reference, molecular mechanics (MM) relaxation of the enzyme by employing the AMBER [17] force filed was performed in which the backbone of the enzyme and the QM region were fixed in position. ONIOM (QM:MM) calculations were carried out at the (DFT:MM) level of theory by utilizing the Gaussian 09 suite of programs [18] and the Tool to Assist ONIOM Calculations (TAO) toolkit [19]. For the QM layer, the meta generalized gradient approximation (GGA) M06-L functional [20] was used in all calculations. We point out that the M06-L functional has been designed to include medium-range electron correlation effects (dispersion effects) and has been specifically recommended for transition-metal containing systems; M06-

130 Page 4 of 12 Fig. 5 Structure (2), optimized intermediate for the binding of NO to heme b3. Bond lengths are in Å

Fig. 6 Optimized intermediate for the binding of NO to the non-heme (FeB) center. Bond lengths are in Å

J Mol Model (2015) 21: 130

J Mol Model (2015) 21: 130

Page 5 of 12 130

Fig. 7 Structure (3), optimized structure leading to N-N bond formation as proposed for mechanism (I). Bond lengths are in Å

L offered the best agreement of several tested functionals compared to large multireference calculations and several benchmark studies have confirmed its excellent accuracy [21–28]. For iron, the Stuttgart/Dresden effective core potential derived from the SDD basis set was utilized [29] whereas the double ζ 6-31G(d,p) basis set was used for the remainder of the atoms. For the MM layer, all atoms were allowed to move without constraints and the AMBER [17] force field was utilized in all calculations. Thus, geometry optimizations were carried out at ONIOM (M06-L/gen:AMBER)-mechanical embedding (ME) level of theory; subsequent optimizations of the resulting structures were carried out at the ONIOM (M06-L/gen:AMBER)-electronic embedding (EE) level of theory. Gradient-based minimization methods using the unrestricted formulism were used for minimum-energy geometry optimization. Vibrational analyses were carried out for each minimized structure to ensure the absence of imaginary frequencies thus confirming a real minimum. All transition states were characterized by vibrational analysis. All QM:MM calculations were carried out in solvent by employing the conductor-like polarizable continuum model (CPCM) [30] in a dielectric constant of =78.4. As noted by Blomberg et al. [12, 31, 32], the entropy effects resulting from the binding or releasing of di and tri atomic molecules can be significant, such values are assumed to be equal to the translational entropy for the free

molecule and thus estimated at 10.8 and 11.1 kcal mol-1 for NO and N2O respectively. The final energies reported in the present study are free energies in solvent after accounting for enthalpy values, entropy, and zero-point effects. In addition, pure QM calculations were also performed on the QM layer at the same level of theory, the α carbon of the amino acid residues were fixed in position during the geometry optimizations, calculations were carried out in solvent with a dielectric constant of =4 by utilizing the same solvation model. QM energies were reported and compared to those resulted from the QM:MM treatment. All structures in this study were visualized using the UCSF chimera package [33].

Results and discussion In this section, the reaction mechanism as well as key reactive intermediates involved in nitric oxide reduction and nitrous oxide generation in cNOR are investigated and discussed. We start our investigation with the 2-electron reduced (diferrous) state (structure 1). The optimized structure is shown in Fig. 4. The binuclear active site was found to feature a high spin (S=4) ground state where both iron centers (Feb3 and FeB) are in the high spin (S = 2) ferrous state. The

130 Page 6 of 12 Fig. 8 Structure (4), optimized structure of the hyponitrite intermediate. Bond lengths are in Å

Fig. 9 TS4-5, the optimized transition state for N-O bond cleavage and N2O formation. Bond distances are in shown in Å

J Mol Model (2015) 21: 130

J Mol Model (2015) 21: 130 Fig. 10 Structure (5), optimized structure of the totally cleaved N-O bond. Bond lengths are shown in Å

Fig. 11 Structure (6), optimized structure showing the N2O cleavage. Bond lengths are shown in Å

Page 7 of 12 130

130 Page 8 of 12

J Mol Model (2015) 21: 130

Fig. 12 Structure (7), optimized structure of the oxo-bridged intermediate. Bond lengths are shown in Å

antiferromagnetic singlet state (Ms=0) was however almost identical in energy being only 0.06 kcal mol-1 higher in energy than the high spin state indicating that at this point no coupling Fig. 13 The potential energy surface for the reaction mechanism of cNOR as obtained from ONIOM (DFT:MM) calculations

exists between both iron centers — as somewhat expected from a FeB-Feb3 bond distance of 4.72 Å. The geometries featured in the QM layer, as depicted in Fig. 1, are generally

J Mol Model (2015) 21: 130

Page 9 of 12 130

Fig. 14 The optimized structure for the diiron dinitrosyl intermediate. Bond lengths are shown in Å

in line with those reported in previous studies and agree well with geometries of the crystal structure. The first step in the reaction mechanism would constitute the binding of nitric oxide to the binuclear active site. Two possibilities exist for such a step, and NO would in principle have the choice of either binding to the heme b3 center or to the non-heme iron. Both binding modes were investigated and both are depicted below in Figs. 5 and 6. The binding of NO to Feb3 was found exergonic by 11.1 kcal mol-1 relative to structure (1), and, in addition, was found to be more energetically favorable than to FeB by almost 7.5 kcal mol-1. These values were higher by almost 2 kcal mol-1 for QM-only energies. Thus, NO binding to heme b3 is suggested as the first step in the reaction mechanism of cNOR. The optimized structure depicted in Fig. 5 (structure 2) was found to feature an overall S=3/2 ground state, atomic spin densities of 0.29, 4.01, and -1.1 on Feb3, FeB, and NO respectively indicate a high spin FeB center and a low spin Feb3 bound to NO. Excess positive spin on Feb3 and excess negative spin on NO would suggest a slightly oxidized Feb3 center where Fe(III) is starting to form. The next step in the reaction mechanism is proposed to feature the binding of a second NO molecule to the binuclear

active site to form a hyponitrite adduct. For this reaction step, two possible pathways can be followed; in the first pathway, a second NO molecule in a free form would bind directly to the heme-bound NO to form a hyponitrite adduct (see mechanism (I) in Fig. 1). Conversely, a second pathway would feature the binding of the second NO molecule to the non-heme iron center to form a diiron dinitrosyl intermediate (cf. mechanism (II) in Fig. 1). Both scenarios are investigated below. Figure 7 depicts structure (3), the optimized starting structure for the first pathway (the binding of a free NO to the heme-bound NO), a free NO molecule is placed at a distance of 3 Å from the nitrogen atom of the heme-bound NO and the energy required for the hyponitrite N-N bond formation is calculated by a systematic relaxed scan of the N-N bond distance. This optimized structure sits on a (S=2) spin state with 0.72 spin units on the free NO and almost unchanged spins on the rest of the key atoms relative to the preceding structure. NN bond formation was found to be exergonic by 5.8 kcal mol-1 (this values decreased to 2.2 kcal mol-1 for QM-only energy) and proceeded with no energy barrier, thus the optimized structure of the hyponitrite intermediate, structure (4), as depicted in Fig. 8, is exergonic by 2.2 kcal mol-1 for QM energy and 5.8 kcal mol-1 for the total QM:MM energy

130 Page 10 of 12

relative to structure (2), and thus exergonic by 16.9 kcal mol-1 relative to structure (1) in terms of the QM:MM energy. It can be noticed that the terminal oxygen atom of the NO molecule bound to Feb3 is now bound to the non-heme FeB center with a FeB-O bond distance of 2.06 Å. The N-N bond distance is now 1.39 Å and Feb3-Nhyponitrite bond has increased from 1.78 to 1.92 Å. Spin densities of 4.11 on FeB and -0.94 on Feb3 suggest antiferromagnetic coupling between a low spin Fe(III) heme b3 and a high spin non-heme center with a state between ferrous and ferric. Spin population on the formed hyponitrite moiety accounts for 0.40 spin units with almost all the spin localized on the O atom bound to the non-heme center (cf. Fig. 8). N-O bond cleavage to yield N2O bound to heme b3 and oxo-bound non-heme FeB center is proposed to be the next step in the mechanism. The transition state for this reaction step (TS4-5) is depicted in Fig. 9. The transition state (TS4-5) was located at an N-O distance of 1.91 Å and the energy requirement for this step was found to be endothermic by 23 kcal mol-1 for QM-only energy and 22.6 kcal mol-1 for the total QM:MM energy relative to structure (4). Thus, this energy barrier by far constitutes the rate limiting step of this mechanism. FeB-O bond distance of 1.78 Å as well as atomic spin of 3.68 on FeB indicates a Fig. 15 The optimized structure of the trans hyponitrite intermediate. Bond lengths are shown in Å

J Mol Model (2015) 21: 130

Fe(IV)-oxo structure starting to form. Feb3-Nhyponitrite is now increased to 2.02 Å and that of N-N is shortened to 1.21 Å relative to structure (4). The spin population on Feb3 stays unchanged at -0.89 spin units. The optimized structure, (structure 5), for the total cleavage of the N-O bond is depicted in Fig. 10. The QM energy barrier for this reaction step was found to be exergonic by 10.5 kcal mol-1 whereas for the total QM:MM energy, this value was 12.3 kcal mol-1 relative to TS4-5. FeB (IV)-oxo is now clearly visible with FeB-O bond length of 1.66 Å coupled with 3.28 and 0.44 spin units on FeB and O respectively. The spin population on Feb3 is now diminished with only -0.14 spin units indicating that the heme b3 center is now returning to the ferrous state. In order to estimate the energy requirement for the liberation of N2O, a systematic scanning of the Feb3-N bond distance is undertaken and the final structure is depicted in Fig. 11 (structure 6). Cleaving the Feb3-N bond was found to be endergonic by 7.3 kcal mol-1 in terms of QM-only energy, whereas for the total QM:MM energy this value was reduced to 5 kcal mol-1 relative to the preceding structure (5). After the liberation of nitrous oxide, an oxo-bridge is suggested to form between both iron centers. The optimized structure for this intermediate is depicted in Fig. 12 (structure 7) and features an oxo bridge with bond distances of 1.92 and

J Mol Model (2015) 21: 130

2.17 Å for Feb3-O and FeB-O respectively and Feb3-FeB bond distance of 4.22 Å. The formation of this intermediate is found to be energetically favorable being almost 34 kcal mol-1 lower in energy than structure (1) as evident from QM:MM results. The optimized structure now features two ferric centers with an antiferromagnetic singlet ground state which is in line with experimental results reported for this structure [4]. The oxo-bridged structure is then suggested to undergo a two electron reduction and double protonation of the oxo bridged atom to yield water and the diferrous active site. The energy diagram of the reaction mechanism is depicted in Fig. 13. The highest energy barrier was found to be the cleavage of the N-O bond and thus proposed to be the ratelimiting step for the reaction mechanism. As mentioned earlier, an alternative pathway for the reaction mechanism would involve the binding of the second NO to the non-heme FeB center to yield a diiron dinitrosyl intermediate (cf. mechanism (II) in Fig. 1). The optimized structure of this intermediate is shown in Fig. 14. This intermediate was found to feature a S=3/2 ground state with −0.33 and 4.09 spin units on Feb3 and FeB respectively suggesting slightly oxidized diferrous centers. This structure was found to be 8.3 kcal mol-1 higher in energy than structure (1). The optimized structure for the trans hyponitrite intermediate is shown in Fig. 15. Spin populations of -0.99 and 4.38 on Feb3 and FeB respectively now suggest that both centers are in the ferric state; however, this intermediate was found to be 9.1 kcal mol-1 higher in energy relative to the preceding structure. Thus, the combined overall energy requirement for the formation of the dinitrosyl and the trans hyponitrite intermediates is 17.4 kcal mol-1 relative to structure (1) and as much as 34.3 kcal mol-1 higher in energy when compared to structure (4). This high energy demand clearly disfavors the dinitrosyl (trans hyponitrite) mechanism for NO reduction in cNOR. The results acquired in this study are generally in line with those reported by Shiraishi and co-workers [13] concerning the key reactive intermediates, the rate-limiting step of the reaction, as well as the importance of the penta-coordination conformation of the non-heme center in providing energetically accessible barriers for the reactive intermediates. Blomberg et al. have also shown the importance of this ligand arrangement in their earlier study [11]; however, their recent QM results [12] did not support this conformation. On the basis of the results of the present study, the dinitrosyl (trans) mechanism has been disfavored due to high energy barriers. The reaction steps involved in this mechanism are still not completely understood mainly the transfer of electrons and protons during the reaction which are usually difficult to study computationally; there is thus a possibility that an energetically accessible pathway for this mechanism might still exist, thus the dinitrosyl (trans) mechanism should not be completely ruled out.

Page 11 of 12 130

Conclusions To summarize, we carried out an ONIOM (DFT:MM) study to investigate the reaction mechanism of cNOR and identify the possible reactive intermediates of the reaction. The results show that the mechanism proceed via the binding of a NO molecule to heme b3, a hyponitrite bond formation via the binding of a second NO in a free form to the heme-bound NO then takes place with an exothermic activation energy. N-O bond cleavage follows to yield N2O and oxo-bound non-heme center and was found to be the rate limiting step with an activation barrier of almost 22.6 kcal mol-1. The dinitrosyl (trans) mechanism was also explored but found unfavorable due to high energy barriers of the resulting intermediates. Thus, this study provides an insight toward understanding the mechanism of NO reduction by cNOR. Acknowledgments Funding from the Romanian Ministry of Education and Research (Grant PN-II-ID-PCE-2012-4-0488) is gratefully acknowledged.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.

15. 16. 17.

18.

Wasser IM, de Vries S, Meönne-Loccoz P, Schröder I, Karlin KD (2002) Chem Rev 102:1201 Hino T, Matsumoto Y, Nagano S, Sugimoto H, Fuku-mori Y, Murata T, Iwata S, Shiro Y (2010) Science 330:1666 Moenne-Loccoz P, de Vries S (1998) J Am Chem Soc 120:5147 Girsch P, de Vries S (1997) Biochim Biophys Acta 1318:202 Kumita H, Matsuura K, Hino T, Takahashi S, Hori H, Fukumori Y, Morishima I, Shiro Y (2004) J Biol Chem 279:55247 Ye RW, Averill BA, Tiedje JM (1994) Appl Environ Microbiol 60: 1053 Butler CS, Seward HE, Greenwood C, Thomson AJ (1997) Biochemistry 36:16259 Moënne-Loccoz P (2007) Nat Prod Rep 24:610 Collman JP, Dey A, Yang Y, Decréau RA, Ohta T, Solomon EI (2008) J Am Chem Soc 130:16498 Yeung N, Lin YW, Gao YG, Zhao X, Russell BS, Lei L, Miner KD, Robinson H, Lu Y (2009) Nature 462:1079 Blomberg LM, Blomberg MR, Siegbahn PE (2006) Biochim Biophys Acta 1757:240 Blomberg MRA, Siegbahn PEM (2012) Biochemistry 51:5173 Shoji M, Hanaoka K, Kondo D, Sato A, Umeda H, Kamiya K, Shiraishi K (2014) Mol Phys 112:393 Berendsen HJC, van der Spoel D, van Drunen R (1995) GROMACS: a message-passing parallel molecular dynamics implementation. Comp Phys Commun 91:43 Lindahl E, Hess B, van der Spoel D (2001) Gromacs 3.0: a package for molecular simulation and trajectory analysis. J Mol Model 7:306 Rostkowski M, Olsson MHM, Søndergaard CR, Jensen JH (2011) BMC Struct Biol 11:6 Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM Jr, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA (1995) J Am Chem Soc 117:5179 Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K,

130 Page 12 of 12

19. 20. 21. 22.

Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery JA Jr, Peralta JE, Ogliaro F, Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega N, Millam JM, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas Ö, Foresman JB, Ortiz JV, Cioslowski J, Fox DJ (2009) Gaussian 09, Revision C.01. Gaussian, Inc, Wallingford Tao P, Schlegel HB (2010) A toolkit to assist ONIOM calculations. J Comput Chem 31:2363 Zhao Y, Truhlar DG (2006) J Chem Phys 125:194101 Cramer CJ, Gour JR, Kinal A, Wtoch M, Piecuch P, Shahi ARM, Gagliardi L (2008) J Phys Chem A 112:3754–3767 Zheng J, Zhao Y, Truhlar DG (2007) J Chem Theory Comput 3: 569–582

J Mol Model (2015) 21: 130 23. 24. 25. 26. 27. 28.

Torker S, Merki D, Chen P (2008) J Am Chem Soc 130:4808–4814 Zhao Y, Truhlar DG (2008) Theor Chem Accounts 120:215–241 Zhao Y, Truhlar DG (2008) Acc Chem Res 41:157–167 Korth M, Grimme S (2009) J Chem Theory Comput 5:993–1003 Zhao Y, Truhlar DG (2009) J Chem Theory Comput 5:324–333 Attia AA, Lupan A, Silaghi-Dumitrescu R (2013) RSC Adv 3: 26194 29. Dolg M, Wedig U, Stoll H, Preuss H (1987) J Chem Phys 86:866 30. Cossi M, Rega N, Scalmani G, Barone V (2003) J Comput Chem 24:669 31. Blomberg LM, Blomberg MRA, Siegbahn PEM (2004) J Biol Inorg Chem 9:923–935 32. Blomberg LM, Blomberg MRA, Siegbahn PEM, van der Donk WA, Tsai A-L (2003) J Phys Chem B 107:3297–3308 33. Pettersen EG, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, Ferrin TE (2004) UCSF chimera—a visualization system for exploratory research and analysis. J Comput Chem 25:1605