Theoretical Insights into the Aerobic Hydrogenase Activity of Molybdenum–Copper CO Dehydrogenase

Abstract

The Mo/Cu-dependent CO dehydrogenase from O. carboxidovorans is an enzyme that is able to catalyse CO oxidation to CO${}_2$; moreover, it also expresses hydrogenase activity, as it is able to oxidize H${}_2$. Here, we have studied the dihydrogen oxidation catalysis by this enzyme using QM/MM calculations. Our results indicate that the equatorial oxo ligand of Mo is the best suited base for catalysis. Moreover, extraction of the first proton from H${}_2$ by means of this basic centre leads to the formation of a Mo–OH–Cu${}^{\mathrm{I}}$H hydride that allows for the stabilization of the copper hydride, otherwise known to be very unstable. In light of our results, two mechanisms for the hydrogenase activity of the enzyme are proposed. The first reactive channel depends on protonation of the sulphur atom of a Cu-bound cysteine residues, which appears to favour the binding and activation of the substrate. The second reactive channel involves a frustrated Lewis pair, formed by the equatorial oxo group bound to Mo and by the copper centre. In this case, no binding of the hydrogen molecule to the Cu center is observed but once H${}_2$ enters into the active site, it can be split following a low-energy path.

1. Introduction

Carbon monoxide dehydrogenases (CODHs) have proved to be an essential component for the biogeochemical carbon monoxide (CO) consumption. They contribute to maintenance of sub-toxic concentration of CO in the lower atmosphere by processing approximately 2 × 10${}^8$ tons of it annually [1]. To date, only two enzymes have been found to be able to use CO as carbon and energy source [2]. The first type is represented by the oxygen-sensitive Ni,Fe-dependent CO dehydrogenase enzyme (NiFe-CODH), found in anaerobic bacteria and archea, whereas in aerobic carboxido-bacteria, such as Oligotropha carboxidovorans, the same function is performed by the oxygen-tolerant MoCu-dependent CO dehydrogenase (MoCu-CODH) enzyme [3]. These bacterial enzymes catalyze the oxidation of CO to CO${}_2$ following the reaction:

CO + H₂O → CO₂ + 2H⁺ + 2e⁻

In addition to this reaction, both enzymes show the ability to catalyze other reactions, although at lower activity. Reduction of CO${}_2$, following the reverse mechanism, is explicated by NiFe-CODH, while the ability of oxidizing H${}_2$ to protons has been reported for MoCu-CODH [3,4].

The latter belongs to the xanthine oxidase family of enzymes, whose members are usually characterised by a mononuclear redox-active site, LMo${}^{\mathrm{VI}}$OS(OH), with a square-pyramidal coordination geometry [5]. The oxo (O${}^{2-}$) group is an apical ligand and the equatorial plane has a dithiolene ligand from a molybdopterin–cytosine dinucleotide (MCD) cofactor, a catalytically labile hydroxyl group (OH${}^{-}$) and a sulphide ion (S${}^{2-}$) [6]. MoCu-CODH presents a noncanonical binuclear active site, in which the sulphido ligand bridges to a second metal, –S–Cu–S-Cys and the equatorial hydroxyl group is deprotonated to form another oxo group. The copper ion is found in a +1 oxidation state and it is also coordinated by a weakly-bound water molecule [7].

The protein is a ($\mathsf{\alpha }\mathsf{\beta }\mathsf{\gamma }$)${}_2$ heterodimer. The active site is located in the large subunit (CoxL; 89 kDa), the medium subunit (CoxM; 30 kDa) contains a FAD cofactor, whereas the small one (CoxS; 18 kDa) possesses two [2Fe–2S] iron–sulfur clusters [7].

The protein resting state is characterised by a Mo${}^{\mathrm{VI}}$(=O)SCu${}^{\mathrm{I}}$ core. Binding of either CO or H${}_2$ is believed to occur at the Cu${}^{\mathrm{I}}$ centre, allowing the proper placement and activation of the substrate [8,9]. The subsequent steps lead to the transfer of two electrons to the Mo ion. This process is favoured by the presence of a highly delocalised redox-active orbital over the Mo–$\mathsf{\mu }$S–Cu moiety [10]. The catalytic activity of MoCu-CODH is also promoted by the presence of highly conserved amino acids in the metal second coordination sphere [11]. A glutamate residue (Glu763) positioned in proximity of the equatorial oxygen ligand of molybdenum is believed to promote deprotonation events [9,12]. In addition, a phenylalanine (Phe390), situated in a flexible loop in front of the copper centre, has been proposed to contribute to the proper orientation of the Cu-bound ligand [13].

Notwithstanding the fact that the general features of the catalytic reaction cycle of the enzyme has been clarified thanks to previous experimental and theoretical efforts, a detailed understanding of the reaction mechanisms for the oxidation of both CO and H${}_2$ has not yet been reached. As for H${}_2$, a mechanism for the dihydrogen molecule splitting was proposed by Wilcoxen and Hille (see Figure 1) [9]. According to electron paramagnetic resonance (EPR) spectroscopy studies, H${}_2$ is believed to coordinate Cu${}^{\mathrm{I}}$ in a side-on fashion. Once the hydrogen is bound, it would be activated by the metal and could easily be split by extraction of a proton thanks to a nearby base. Subsequently, a second proton would leave the active site and the two electrons would be transferred to the Mo centre. To date, no theoretical data have been reported in support of the experimental proposal. Only one previous theoretical work investigating the hydrogenase activity has been published so far [14]. In this study, it was unexpectedly found that the copper ion could not bind the H${}_2$ substrate. Therefore, the authors suggested that changes in the active site are needed, for example, by protonation events, before the catalytic cycle can start. In such context, it is interesting to notice that, to the best of our knowledge, side-on binding of H${}_2$ to molecular species containing a Cu(I) ion has been previously reported only in a couple of gas-phase case studies, that is, H${}_2$–CuCl [15] and H${}_2$–CuF [16]. Somehow, the interaction between H${}_2$ and the enzyme appears to find a closer parallel with the case of porous materials containing electron-rich copper ions, such as certain zeolites proposed for H${}_2$ storage applications [17,18,19] and Prussian-blue analogues [20].

In the present work we present a plausible mechanistic picture for the MoCu-CODH hydrogenase activity, taking into account the effects of the protein matrix by means of hybrid quantum mechanical/molecular mechanical (QM/MM) enzyme models.

2. Results and Discussion

2.1. Study of H${}_2$ Binding Modes to the Copper Centre

We started with an investigation of the most plausible binding sites for the H${}_2$ substrate. Considering the resting state of the protein, two types of coordination to the copper centre were considered (1HR and 2HR in Figure 2). We find that H${}_2$ can bind to the Cu ion both when the latter is bi-coordinated ($\mathsf{\mu }$S–Cu–S-Cys388) and also when the Cu coordination sphere includes a weakly bound water molecule, as proposed in previous investigations of the enzyme resting state [7,13]. The hydrogen molecule is found to be slightly activated, with a H–H distance of 0.77 Å in both cases (the equilibrium distance of H${}_2$ at BP86/def2-TZVP level is 0.75 Å).

For both binding modes, we estimated the binding energy of H${}_2$ to Cu${}^{\mathrm{I}}$ by comparing the energy of the coordinated structure to the corresponding one in which H${}_2$ is not linked to copper but it is already present in the active site (1HNoB and 2HNoB in Figure 3). The resulting QM/MM total energy differences were found to be +46.9 and +46.4 kJ/mol, in the absence (1H) or presence of the coordinated water molecule (2H), respectively (see

Table 1. Relative energies (kJ/mol) calculated using the quantum mechanical/molecular mechanical (QM/MM) approach at B3LYP-D3(BJ)/def2-TZVPD  level.

	NoB	R	P1	TS	P	TS2	P2
1H	−46.9	0.0	57.7		−2.9	31.8	−130.1
2H	−46.4	0.0	83.3
3H	−10.9	0.0		29.3	28.0	56.4	−90.8
1H-FLP	0.8	0.0		45.2	44.4	79.1	−82.8

).

Thus, formation of these reactants was found to be strongly unfavourable from an energetic point of view. Therefore, we contemplated the possibility that a sterically hindered Lewis acid–Lewis base pair (a so-called frustrated Lewis pair, FLP), may form in the active site, which would eventually lead to heterolytic H${}_2$ cleavage. Indeed, combination of Lewis acid and bases have been demonstrated to activate a wide range of small molecules [21]. Moreover, the plausible involvement of an FLPs in enzymatic catalysis has been recently proposed for the H${}_2$ splitting by [NiFe]-hydrogenase [22] and an FLP might actually be active also in the case of [FeFe]-hydrogenases [23,24,25].

We found a variant of 1HR (1H-FLP_R in Figure 2), in which H${}_2$ assumes a different position in the active site, such that the Mo–O${}_{\mathrm{eq}}$ atom and the copper centre appear to be able to act as a FLP. The 1H-FLP_R reactant is characterised by a H${}_2$ molecule positioned in front of the copper atom, with long Cu–H${}_{{\mathrm{H}}_2}$ distances (2.80 and 2.93 Å) and a O${}_{\mathrm{eq}}$–H${}_{{\mathrm{H}}_2}$ distance of 2.33 Å. 1H-FLP_R was found to be 0.8 kJ/mol lower in energy than the 1HNoB species in which the hydrogen molecule is separated from the bimetallic centre.

Based on previous literature, which proposed that the active site could be protonated prior to substrate binding [14,26], we also evaluated the influence of protonation on the coordination of H${}_2$ to the metal. Protonation of the bridging sulphide ion turned out not to favour dihydrogen binding, as no H${}_2$-bound geometry was found on the QM/MM potential energy surface. On the other hand, protonation of the sulphur atom of Cys388 favours coordination of H${}_2$ to Cu. The latter process would lead to the formation of a H${}_2$-coordinated geometry (3HR in Figure 2), in which dihydrogen is activated to a H–H bond of 0.78 Å and the computed binding energy is +10.9 kJ/mol.

It is noticeable that energy minimizations of the one-electron reduced counterparts of 1HR, 2HR and 3HR lead to the detachment of the hydrogen molecule from the copper centre. This result is not compatible with the experimental outcomes that identified a H${}_2$-bound Mo${}^{\mathrm{V}}$Cu${}^{\mathrm{I}}$ species by means of EPR experiments [9]. However, since one-electron reduced counterparts do not represent plausible intermediates in the catalytic cycle—there are no evidences that H${}_2$ binding is followed by Mo reduction, prior to substrate oxidation—we have not investigated this aspect any further. The aim of the work is focused on the mechanistic characteristics of the H${}_2$ oxidation reaction, the details of which are described in the subsequent sections.

2.2. Exploring Basic Residues in the Active Site

Next, we investigated which functional group may represent the base for the abstraction of the first proton from the activated hydrogen molecule. We considered the Glu763 residue, which has been suggested to take part into deprotonation events (vide infra) and the equatorial oxo ligand of molybdenum. As for the reactants, we considered not only 1H-FLP_R and 3HR but also 1HR and 2HR. In fact, even though H${}_2$ binding in the latter two cases was computed to be disfavoured, one cannot completely rule out the hypothesis that such adducts have some minor role in dihydrogen oxidation catalysis, a possibility that would be enabled in case the subsequent reaction steps are associated with an overall smooth energy profile.

The two bases would lead to the formation of two copper hydrides (see Figure 4), namely either GluOOH–Cu${}^{\mathrm{I}}$H (P1) or Mo–OH–Cu${}^{\mathrm{I}}$H (P). For the 1HR and 2HR states (see Figure 2), the oxygen of the glutamate residue is closest (the O${}_{\mathrm{Glu}}$–H_H2 distance is 1.87 Å and 1.89 Å for 1HR and 2HR), while the Mo–oxo ligand is 3.09 and 3.28 Å from the closest H atom of H${}_2$, respectively. However, for 3HR, O${}_{\mathrm{eq}}$ is the closest base (the O–H distance is 1.91 Å) while the oxygen of Glu763 is 3.83 Å away.

The search for the GluOOH–Cu${}^{\mathrm{I}}$H products was found to be challenging, as the glutamic acid sidechain in the guess geometries turned out to release a proton to reform a Cu-bound H${}_2$ along optimizations. Therefore, to optimize the intermediates 1HP1 and 2HP1, (see Figure 5) we had to fix the O${}_{\mathrm{Glu}}$–H distance; this was done in the hypothesis of the existence of a barrier for proton transfer, although so small that convergence on the desired minimum could be hampered by the choices made on the starting input structures to be optimized. The total energy differences of these structures, with respect to the corresponding reagents, were found to be 57.7 and 83.3 kJ/mol for 1H and 2H, respectively. In the case of 3H, no GluOOH–Cu${}^{\mathrm{I}}$H (3HP1) product was found.

The product formed by extraction of the proton by the equatorial ligand of Mo (P) was not found in the case of 2H, whereas it corresponds to structures with a total energy difference of −2.9 and 28.0 kJ/mol with respect to the reagents for 1H and 3H, respectively (1H-FLP_P in Figure 6 and 3HP in Figure 7). 1H-FLP_P can also represent the product for the H${}_2$ splitting reaction starting from the 1H-FLP_R structure. In this case, the resulting Cu–hydride was found to be 44.4 kJ/mol less stable than the reagent.

2.3. Plausible Activation Mechanisms for H${}_2$ Splitting

Based on the above results on the relative energies of the intermediates potentially involved in the catalysis, we proceeded with the computation of the activation barriers in order to estimate which reaction mechanism would represent the most energetically favourable path for the hydrogenase activity of MoCu-CODH.

Starting from the H${}_2$-bound conformation 1HR, we have just shown that the equatorial Mo–oxo ligand represents the most favourable basic group for the first proton abstraction, leading to a rather stable product (1H-FLP_P). However, despite extensive efforts, we were not able to locate any transition state for the reaction 1HR → 1H-FLP_P. This issue is probably due to the large distance between the labile proton and the base. We also considered the 1HR → 1HP1 reaction because, even if the GluOOH–Cu${}^{\mathrm{I}}$H product is not as stable as the Mo–OH–Cu${}^{\mathrm{I}}$H product, the whole catalysis may be driven by subsequent exoenergetic reactions. However, again, no transition state was found.

On the other hand, the 1H-FLP_R species can be linked to the 1H-FLP_P product through a low-lying transition state (1H-FLP_TS), reported in Figure 6. The reaction 1H-FLP_R→ 1H-FLP_P was found to be endergonic ($\Delta$E = 44.4 kJ/mol) with an associated activation barrier of 45.2 kJ/mol. The latter results support the hypothesis that Cu and the totally oxidised equatorial ligand of molybdenum could play the role of a FLP.

For the complex with a water molecule at the copper centre (2HR), we could not locate any transition state geometry that would allow the deprotonation of the hydrogen molecule by the anionic Glu763 residue.

Finally, we considered protonation of Cys388 and found that such protonation can be functional for the hydrogenase activity. In fact, in this case the activation energy for H${}_2$ cleavage involving the equatorial oxo-ligand is only 29.3 kJ/mol (3HTS in Figure 7), that is, 15.9 kJ/mol lower than for 1H-FLP.

Independently of the protonation state of Cys388, our results support a mechanism in which the equatorial oxo group of Mo extracts the first proton. Such process is characterised by a rather low activation barrier that may allow a relatively rapid proton exchange, in accordance with experimental evidence [9]. However, differently from what was hypothesised in the literature [9], we propose that this rapid proton exchange is promoted by the oxo ligand on Mo, rather than by Glu763. Moreover, the equilibrium of this step favours the H${}_2$ reactant complex over copper hydride in both cases, as discussed by Wilcoxen and Hille in light of the pH independence of the reaction. Finally, the intermediates 3HP or 1H-FLP_P should be converted to the fully reduced form of the enzyme by the transfer of the Cu-bound hydride to the (now protonated) O${}_{\mathrm{eq}}$, giving Mo${}^{\mathrm{IV}}$O(OH${}_2$)Cu${}^{\mathrm{I}}$ [8,13,27] (P2; see Figure 8). In both cases, the resulting reduced intermediates were found to be very stable, −90.8 and −82.8 kJ/mol for 3HP2 and 1H-FLP_P2, respectively (see Table 1 and Figure 6 and Figure 7). Moreover, the reaction barriers for the above two reactions are relatively low (28.4 and 34.7 kJ/mol, respectively). Finally, the reverse reactions are energetically impeded since the reverse activation barriers exceed 140 kJ/mol in both cases (TSs energies are and 147.2 and 161.9 kJ/mol respectively). Again, this is in good agreement with experiments, showing no production of H${}_2$ by the enzyme [9].

3. Methods

3.1. The Protein

All calculations were based on the crystal structure of MoCu-CODH hydrogenase in its oxidised form (PDB ID: 1N5W) [7]. Only the large subunit (CoxL) of one monomer, containing the active site, was considered in this study. The enzyme was set up in the same way as in our previous study [28]. The protonation state of all the residues was determined based on calculations with PROPKA [29] and on studies of the hydrogen-bond pattern, of the solvent accessibility and of the possible formation of ionic pairs. All Arg, Lys, Asp and Glu residues were assumed to be charged, with exception of Glu29 and Glu488 that were protonated on OE2, whereas Asp684 was protonated on OD1. Cysteine ligands coordinating to metals were deprotonated. Among the His residues, His61, 339, 766 and 793 were protonated on the ND1 atom, His177, 178, 210, 213, 243, 700, 753, 754 and 788 were assumed to be protonated on NE2 atom, whereas the other His residues were assumed to be doubly protonated. The protein was solvated with water molecules, forming a sphere with a radius of 60 Å around the geometric centre of the protein. The added protons and water molecules were then optimised by a 1-ns simulated-annealing molecular-dynamics simulation, followed by energy minimization [28] (for details on the adopted force field, see the next subsection).

3.2. QM/MM Calculations

The QM/MM calculations were performed with the ComQum software [30,31]. According to this approach, the protein and the solvent are split into three subsystems—System 1 corresponds to the QM region and it was relaxed by QM methods. It consisted of the molybdenum ion and its first coordination sphere (two O${}^{2-}$, the bridging sulphide and the MCD cofactor, truncated to exclude the phosphate and cytosine moieties), the copper ion, the ligand molecules H${}_2$, H${}_2$O, as well as the sidechains of residues Cys388 and Glu763 (see Figure 9).

System 2 consists of all residues with any atom within 6 Å of any atom in System 1, whereas System 3 contains the remaining part of the protein and the water molecules. The latter two systems were kept fixed at the crystallographic coordinates in the QM/MM calculations.

In the QM calculations, System 1 was represented by a wavefunction whereas all the other atoms were represented by an array of partial point charges.

When there is a covalent bond between the QM and MM systems, the QM system was truncated using the hydrogen link-atom approach [32]. The QM system was capped with hydrogen atoms (hydrogen link atoms, HL), the position of which are linearly related to the corresponding carbon atoms (carbon link atoms, CL) in the full system [30]. All atoms were included in the point-charge model, except the CL atoms. ComQum employs a subtractive scheme with electrostatic embedding and van der Waals link-atom corrections [33]. The total QM/MM energy is calculated as

$$E_{\mathrm{QM}/\mathrm{MM}}=E_{\mathrm{QM}1+\mathrm{ptch}23}^{\mathrm{HL}}+E_{\mathrm{MM}123,q_1=0}^{\mathrm{CL}}-E_{\mathrm{MM}1,q_1=0}^{\mathrm{HL}}$$

(1)

where E${}_{\mathrm{QM}1+\mathrm{ptch}23}^{\mathrm{HL}}$ is the QM energy of System 1 truncated by HL atoms and embedded in the set of point charges modeling Systems 2 and 3. E${}_{\mathrm{MM}1,q_1=0}^{\mathrm{HL}}$ is the MM energy of System 1, truncated by HL atoms, without any electrostatic interactions. E${}_{\mathrm{MM}123,q_1=0}^{\mathrm{CL}}$ is the classical energy of the whole system, with CL atoms and with the charges in System 1 set to zero, to avoid double counting of the electrostatic interactions. The QM calculations were carried out using TURBOMOLE 7.1 software [34]. Geometry optimizations and TS searches (the latter based on minimizations with geometric constraints imposed step-wise to selected atoms, along the putative reactive path) were performed using the BP86 functional [35,36] in combination with def2-TZVP basis set [37]. All calculations included Grimme’s dispersion correction with Becke–Johnson damping (D3(BJ)) [38]. The resolution-of-identity technique was employed to accelerate the calculations [39].

The MM calculations were performed with the Amber software [40], using the Amber ff14SB force field for the protein [41] and the general Amber force field [42] with restrained electrostatic potential (RESP) charges [43] for H${}_2$ and MCD. The two Fe${}_2$S${}_2$ clusters were described with RESP charges and a non-bonded model (they are kept fixed in the calculations).

4. Conclusions

We have studied dihydrogen oxidation by MoCu-CODH using QM/MM calculations. Our results indicate that the equatorial oxo ligand of Mo is a better base than Glu763 during catalysis. Moreover, extraction of the first proton from H${}_2$ by means of this basic centre leads to the formation of a Mo-OH–Cu${}^{\mathrm{I}}$H hydride that allows for the stabilization of the copper hydride, otherwise known to be very unstable [44]. It is intriguing to notice that our proposal of a direct involvement of the Mo-bound oxo group in H${}_2$ splitting finds a conceptually similar case in a recently published mechanistic study on the [Fe]-hydrogenase. In the latter enzyme, a deprotonared OH group that belongs to the iron–guanylylpyridinol cofactor and that is in $\mathsf{\gamma }$-position with respect to the metal-activated H${}_2$, was suggested to be involved in substrate splitting [45]. Such mechanistic picture on the [Fe]-hydrogenase was the result of crystallographic studies that took advantage of the fact that the latter enzyme exists in two forms—an open substrate-accessible form and a closed catalytically active form. This is not the case for MoCu-CODH; however, key details on the proton transfer events following H${}_2$ binding might still be gained at experimental level, most likely by means of application of high-sensitivity infrared spectroscopy techniques (see for example, Reference [46]).

In light of our results, two plausible mechanisms for the hydrogenase activity of the enzyme can be proposed, as reported in Figure 8. The first reactive channel (cycle A in Figure 8) depends on protonation of the sulphur atom of Cys388, which appears to favour the binding and activation of the hydrogen molecule. The second reactive channel involves a frustrated Lewis pair, formed by the Mo–O${}_{eq}$ oxo group and the copper centre (cycle B in Figure 8). In this case, no binding of the hydrogen molecule to the Cu center is observed but once H${}_2$ enters into the active site, it can be split following a low-energy path. All in all, we cannot exclude that both pathways may be active in MoCu-CODH, thus enhancing the overall catalytic turnover. The two mechanisms for H${}_2$ oxidation are characterized by energy profiles that were found to be quite similar (see Figure 10) and although the FLP-based mechanism features a ~20 kJ/mol higher rate-determining barrier, it does not require any preliminary active-site protonation. However, the H${}^{+}$ source for protonation of Cys388 might be the H${}_2$ molecule itself, as a result of a possible proton transfer following initial H${}_2$ splitting at the active site by means of the FLP mechanism. Alternatively, protonation of Cys388 may occur as a parallel, concerted process at the onset of H${}_2$ oxidation catalysis. Notably, the existence of a pair of distinct proton pathways has been recently suggested for [FeFe]-hydrogenases. There, a ”regulatory” protonation of the [4Fe–4S] portion of the hexanuclear active site was proposed to occur in concomitance with a ”catalytic” protonation step occurring at the di-iron portion of the same active site [47]. The results presented in the latter study and—hopefully—those discussed in the present paper are likely to provide useful information for future efforts to deepen insights into proton transfer mechanisms from and to the MoCu-CODH active site.