## 3 Results and Discussion  
The QM/MM optimized geometries of each conformer are shown in Figure 5 for both Fe\({}^{2+}\) (gray) and Fe\({}^{3+}\) (orange) superimposed on the oxidized structure determined with X-ray crystallography at 1.3 A resolution [1]. In general, the QM/MM optimized geometries reproduce well the structures resolved with X-ray crystallography [1]. Bond lengths and relative orientations of atoms are well preserved. The positions of oxygen atoms from crystal water w365 and w589 are also preserved in each model, suggesting that the methodology is suited to describe the overall electronic environment. Bond lengths of atoms coordinated to a nickel or iron atom are reported in Tables 1 and 2 for [Ni\({}^{2+}\)Fe\({}^{2+}\)Se] and [Ni\({}^{2+}\)Fe\({}^{3+}\)Se], respectively; for comparison, the bond lengths obtained from the crystal structure (PDB ID: 5 JSH from [PERSON] et al. [1]) are listed in bold print. Ni-S lengths for the coordinated cysteine thiolates are well reproduced, but Fe-Ni distances in all three conformers are overestimated by 0.1-0.2 A compared to the crystal structure. Particularly, conformer III demonstrates an overestimation in the Fe-Ni bond length for both Fe\({}^{2+}\) and Fe\({}^{3+}\) oxidation states (2.61 A and 2.76 A, respectively, versus the experimental value of 2.41 A). This Fe-Ni separation is accompanied by a rotation of the side chain of Sec489, which leads to a 1 A shift in the position of the selenium atom (see arrow pointing to Sec489 in Figure 5C). This overestimation in bond lengths is also observed for Ni-S and Ni-Se bonds for conformers I and II. For the oxidized species of [Ni\({}^{2+}\)Fe\({}^{3+}\)Se], reasonable agreement is observed between the computed Fe-Ni bond length of conformer II (2.45 A) and the experimentally reported value of 2.41 A (Table 2). Conformers I and III overestimate the Fe-Ni distance by more than 0.1 A and 0.34 A, respectively, compared to the experimentally measured values [1]. The discrepancies between the computed and measured bond lengths suggest that either the assumed electronic state of the catalytic site may not accurately reflect the state what is captured in the crystal structure or that the assignment of atomic positions may be ambiguous due to insufficient experimental resolution. Here, we have treated the electronic state of nickel as 2+ since [NiFeSe] hydrogenases are known to be EPR silent, in contrast to standard Ni-Fe hydrogenases [8]. Another factor that may lead to discrepancy between computed and experimental geometries may be the existence of a bridging ligand (such as a hydride) that may not be resolved from electron density. Figure 5: The QM/MM optimized geometry of the atoms included in the QM region of conformers I, II, and III for Fe\({}^{2+}\) (atoms shown in gray) and Fe\({}^{3+}\) (atoms shown in orange) are compared; the optimized geometries are superimposed on the geometries obtained from the crystal structure (PDB ID 5 JSH [1]) (atoms shown in color: nickel atom is depicted in light blue, iron atom in green, sulfur and selenium atoms in yellow, carbon atoms in light blue, oxygen atoms in red, and nitrogen atoms in dark blue; hydrogen atoms are omitted for clarity). The bond lengths between nonproteic ligands (CO, CN1, and CN2) coordinated to Fe show reasonable agreement (0.05-0.14 A). In our models, we have adopted the assignment of positions from the crystal structure of [PERSON] et al. [1], namely CO pointing toward Leu425 and CN-ligands pointing toward Pro421 and Ser445. Nonetheless, the experimental resolution (1.3 A) does not allow for unambiguous assignment of ligand position. In other words, if the CO ligand were exchanged with either of the CN-ligands, the electronic structure of the hydrogen bonding network to neighboring amino acid side chains and crystal water molecules would change (see Figure 3 for up-close view of active site). These alternate two scenarios have not been tested yet but would be worthwhile extensions of the current investigation.