A QM/MM Study of Nitrite Binding Modes in a Three-Domain Heme-Cu Nitrite Reductase

Copper-containing nitrite reductases (CuNiRs) play a key role in the global nitrogen cycle by reducing nitrite (NO2−) to nitric oxide, a reaction that involves one electron and two protons. In typical two-domain CuNiRs, the electron is acquired from an external electron-donating partner. The recently characterised Rastonia picketti (RpNiR) system is a three-domain CuNiR, where the cupredoxin domain is tethered to a heme c domain that can function as the electron donor. The nitrite reduction starts with the binding of NO2− to the T2Cu centre, but very little is known about how NO2− binds to native RpNiR. A recent crystallographic study of an RpNiR mutant suggests that NO2− may bind via nitrogen rather than through the bidentate oxygen mode typically observed in two-domain CuNiRs. In this work we have used combined quantum mechanical/molecular mechanical (QM/MM) methods to model the binding mode of NO2− with native RpNiR in order to determine whether the N-bound or O-bound orientation is preferred. Our results indicate that binding via nitrogen or oxygen is possible for the oxidised Cu(II) state of the T2Cu centre, but in the reduced Cu(I) state the N-binding mode is energetically preferred.


Introduction
Microbial copper-containing nitrite reductases (CuNiRs) are enzymes that catalyse the reduction of NO 2 − to NO, a key denitrification step in the global nitrogen cycle [1]. The most studied and well-characterized are the two-domain CuNiRs, which are homotrimeric proteins with each monomer consisting of two cupredoxin-like domains containing an electron-transfer type 1 Cu (T1Cu) site and a catalytic type 2 Cu (T2Cu) site [2,3]. The T1Cu is typically coordinated with Cys-Met-His 2 and the T2Cu is coordinated with His 3 -H 2 O in the resting phase. The two sites are separated by approx.  There is a striking similarity between the core cupredoxin domain of RpNiR and the two-domain AcNiR (Figure 1a,b). However, the additional heme domain necessitates certain rearrangements in the core structure. In two-domain CuNiRs, X-ray studies revealed the presence of two solvent access channels with the potential to deliver protons from the bulk solvent to the T2Cu site during the reduction of NO2 − , one of which also functions as the substrate access channel [10,11]. Structural, spectroscopic and computational studies have revealed the importance of the substrate access channel in controlling the hydration of the active site [12][13][14], the orientation of the nitrite and selectivity for small ligands [15,16], and determining the rate of turnover [12]. In RpNiR, the linker that tethers the heme c domain to the cupredoxin domain is poised in such a way that the sidechain of one of its residues, Tyr323, resides close to the T2Cu site, thereby blocking this substrate channel (Figure 1c). The position of the Tyr323 is critical and it is hypothesized to have mechanistic There is a striking similarity between the core cupredoxin domain of RpNiR and the two-domain AcNiR (Figure 1a,b). However, the additional heme domain necessitates certain rearrangements in the core structure. In two-domain CuNiRs, X-ray studies revealed the presence of two solvent access channels with the potential to deliver protons from the bulk solvent to the T2Cu site during the reduction of NO 2 − , one of which also functions as the substrate access channel [10,11]. Structural, spectroscopic and computational studies have revealed the importance of the substrate access channel in controlling the hydration of the active site [12][13][14], the orientation of the nitrite and selectivity for small ligands [15,16], and determining the rate of turnover [12]. In RpNiR, the linker that tethers the heme c domain to the cupredoxin domain is poised in such a way that the sidechain of one of its residues, Tyr323, resides close to the T2Cu site, thereby blocking this substrate channel ( Figure 1c). The position of the Tyr323 is critical and it is hypothesized to have mechanistic implications. Unusually, the three-domain RpNiR has a remarkably low affinity for binding NO 2 − in the oxidized T2Cu state [7]. Another striking similarity among all CuNiRs is the presence of conserved Asp and His residues in the T2Cu active site, which are believed to participate in the proton transfer that is necessary for reduction [16][17][18][19]. In RpNiR the corresponding Asp and His are Asp97 and His240, respectively. A mutant of RpNiR in which the catalytically important Asp97 was mutated to Asn (D97N) [20] revealed that binding of NO  Figure 2 in [21]), which is thought to occur as a result of reduction of the T2Cu. In addition, the mode of nitrite binding is important for mechanistic understanding of the reduction of nitrite in two-domain CuNiRs, which may proceed according to either an ordered [24][25][26] or random-sequential mechanism [27,28] and there is similar speculation about the mechanism in three-domain RpNiR [7,20]. As binding via either O or N is possible for the nitrite ligand, it is essential to identify the initial binding mode of nitrite in the active site to elucidate the mechanism for RpNiR. Until now no direct structural information of the binding mode of NO 2 − in native RpNiR is available, although binding via nitrogen was observed for the D97N mutant [20]. We have therefore used combined quantum mechanical and molecular mechanical (QM/MM) methods to determine whether binding of NO 2 − occurs via the Nor O-atom in native RpNiR and have compared the calculated structures to the available experimental data for the D97N mutant as well as the available structural data on two-domain CuNiRs.

Molecular Dynamics of Native and D97N RpNiR in the Resting State
Initial classical MD simulations were carried out to investigate the positioning of the crucial active site residues such as Tyr323 and provide snapshots for subsequent QM/MM calculations. Analogously to other CuNiRs, the conserved active site Asp residue, Asp97 in RpNiR is likely a key player for the proton transfer required during catalytic reduction at the T2Cu. We therefore ran simulations of the resting states of three RpNiR trimer systems: the two protonation states of Asp97 (which we label "D97" and "D97p" for the unprotonated and protonated states respectively) and an in silico mutated structure "D97N" where Asp97 is replaced with Asn. Conformational variations within the T2Cu active sites of the three independent subunits (chains A, B and C) of the trimer were analysed over the unrestrained last 75 ns of a total of 80 ns of MD trajectories for each of these three systems (see the methodology section for further details).
As the Tyr323 residue sits in the substrate access channel it interacts with the hydration network connecting Asp97, T2Cu and His240. Figure 2a depicts the dynamic hydrogen bond distance between the phenolic H atom of Tyr323 and the nearest O atom of Asp97/nearest N-or O-atom of Asn97. In the case of D97, a strong hydrogen bond is present between the negatively charged carboxylic acid group of Asp97 and Tyr323. Protonation of Asp97 weakens this interaction and as a result a transient break in this interaction is frequently observed in the D97p system. A similar situation is observed for the D97N mutant with the transient breaks becoming more frequent. Interestingly, there is a complete breaking of the TyrO-H···O-Asp97 interaction at~69 ns in one subunit of the D97 system (chain A in Figure 2a). This is due to movement of the Tyr323 away from the T2Cu pocket (vide infra). To understand the relative position of these two amino acids with respect to the T2Cu ion, the time evolution of the distances between the T2Cu and the centre of mass (COM) of the side chains of Tyr323 and Asp97 during MD were measured and are plotted in Figure 2b, while the averaged values over the entire trajectory for each subunit are given in Table 1. The average distances between Asp97 and T2Cu are similar to those in the crystal structure. The average distance of Asn97-T2Cu for the mutant D97N simulations is larger than that measured in the recent D97N crystal structure [20] (4.30 Å). The other distance, T2Cu and COM of the Tyr323 side chain, in the native crystal geometry is 7.16 Å. The average values for the T2Cu-Tyr323 distance in D97, D97p and D97N systems also agree well with the crystal structure (Table 1). (chain A in Figure 2a). This is due to movement of the Tyr323 away from the T2Cu pocket (vide infra).
To understand the relative position of these two amino acids with respect to the T2Cu ion, the time evolution of the distances between the T2Cu and the centre of mass (COM) of the side chains of Tyr323 and Asp97 during MD were measured and are plotted in Figure 2b, while the averaged values over the entire trajectory for each subunit are given in Table 1. The average distances between Asp97 and T2Cu are similar to those in the crystal structure. The average distance of Asn97-T2Cu for the mutant D97N simulations is larger than that measured in the recent D97N crystal structure [20] (4.30 Å). The other distance, T2Cu and COM of the Tyr323 side chain, in the native crystal geometry is 7.16 Å. The average values for the T2Cu-Tyr323 distance in D97, D97p and D97N systems also agree well with the crystal structure (Table 1).  The exception, as noted above, is in one subunit of D97 (chain A in Figure 2b and Table 1) where the Tyr323 becomes displaced from the T2Cu site by 2.5-3.0 Å after~69 ns. Since the relative positions of Asp97 and T2Cu are unchanged in the same MD timeframe, this implies that it is the flexibility of the linker chain that enables displacement of the Tyr323. Although the orientations of the Tyr323 sidechain in the two cases are different, the displacement seen by MD corroborates the flexibility observed in the D97N crystal structure ( Figure S1), where the sidechain is rotated~90 • away from its position in the native crystal structure (Figure 2a in [20]).
In the D97p and the D97N systems ( Figure 2a, middle and right panels) the local H-bond network between Tyr323 and Asp/Asn97 is frequently disrupted. However, within the timescale of measurement Tyr323 remains 'locked' in its position in the substrate channel and so would sterically hinder the entry of NO 2 − . In the D97 system Tyr323 is held strongly in H-bond interaction by the negatively charged Asp97 but in the case of chain A the Tyr323 nevertheless breaks away towards the end of the simulation (Figure 2b, left panel), becoming displaced in the channel by 2.5-3.0 Å from its original position relative to the T2Cu ( Figure S1). The distance of Tyr323 from T2Cu in chain A (Table 1) at~69 ns is similar to that found by Dong et al. in the crystal structure of NO-and nitrite-bound D97N RpNiR [20].

QM/MM Optimization of NO 2 − Bound Structures
Snapshots of structures from the MD simulations were taken as the basis for studying the binding mode of nitrite to the T2Cu site. The snapshots were selected to be representative of the conformational variations observed in the MD trajectories described in the previous section, taking them from 25 ns onwards to ensure that the systems were completely equilibrated after removal of backbone constraints. The nitrite ion was added to the systems as described in the methodology section. QM/MM optimizations were carried out for both the Cu(II) and Cu(I) states of the T2Cu ion, starting from both 'top-hat' and 'N-bound' orientations of NO 2 − . For each optimisation, the geometry and relative energy of the relaxed structure were assessed. We categorised the relaxed geometries of the T2Cu site in two ways. First, by the overall geometry of the T2Cu-nitrite complex, which can be either tetragonal or trigonal. This was determined by measuring the twisting angle ϕ between two planes: one defined by N(NO 2 − )-T2Cu-N(His1) and the other defined by N(His2)-T2Cu-N(His3), following the procedure used by Ryde et al. for blue and green Cu proteins [29,30]. His1 is the axial ligand, defined as the ligand with the largest N(NO 2 − )-T2Cu-N(His) angle, with His2 and His3 the remaining equatorial ligands. For a tetragonal structure ϕ will be approximately 0 • or 180 • , while for trigonal structures ϕ will be approximately 90 • . The second assessment was of the orientation of the NO 2 − ion bound to the T2Cu centre.
The nitrite orientation can be classified as symmetrical N-bound, L-shaped, bidentate top-hat, monodentate top-hat, or side-on. To assign these labels we define a coordinate system in which the centre of mass of NO 2 − is placed at the origin with the nitrogen atom aligned along the positive X-axis, and the oxygen atoms lie on the XY plane, as shown in Figure 3. Two angular measurements describing the position of Cu in relation to this coordinate system can then be used to classify the different orientations. We define an altitudinal angle θ as the angle between the projection of Cu onto the XZ plane and the X-axis. This angle can be thought of as the 'pitch' of the NO 2 − ion and describes a change in its orientation from N-bound (0 • ) to side-on (90 • ) to top-hat (180 • ). We also define an azimuthal angle ψ as the angle between the projection of Cu onto the XY plane and the X-axis. This angle can be thought of as the 'yaw' of the NO 2 − ion which circles from N-bound (0 • ) to L-shaped (~30 • ) to perpendicular (90 • , not observed in practice) to monodentate top-hat (~150 • ) to bidentate top-hat (180 • ). The thresholds we use for labelling the nitrite ion orientation in each structure are shown in Table 2.
Molecules 2018, 23, x FOR PEER REVIEW 6 of 16 The second assessment was of the orientation of the NO2 − ion bound to the T2Cu centre. The nitrite orientation can be classified as symmetrical N-bound, L-shaped, bidentate top-hat, monodentate top-hat, or side-on. To assign these labels we define a coordinate system in which the centre of mass of NO2 − is placed at the origin with the nitrogen atom aligned along the positive Xaxis, and the oxygen atoms lie on the XY plane, as shown in Figure 3. Two angular measurements describing the position of Cu in relation to this coordinate system can then be used to classify the different orientations. We define an altitudinal angle θ as the angle between the projection of Cu onto the XZ plane and the X-axis. This angle can be thought of as the 'pitch' of the NO2 − ion and describes a change in its orientation from N-bound (0°) to side-on (90°) to top-hat (180°). We also define an azimuthal angle ψ as the angle between the projection of Cu onto the XY plane and the X-axis. This angle can be thought of as the 'yaw' of the NO2 − ion which circles from N-bound (0°) to L-shaped (~30°) to perpendicular (90°, not observed in practice) to monodentate top-hat (~150°) to bidentate top-hat (180°). The thresholds we use for labelling the nitrite ion orientation in each structure are shown in Table 2. The energies of the relaxed conformers were assessed to determine whether there was any preference for one binding mode over another. Absolute energies between the conformers cannot be compared since the total number of atoms of the conformers is not conserved as they are taken from random time steps and a shell of solvent is used. Instead a relative energy was calculated from the difference in energy between relaxed top-hat and N-bound structures (ΔE = Etop-hat -EN-bound) such that a positive ΔE indicates that N-bound is favoured while a negative value indicates top-hat is favoured. Figure 2 indicates that within the timescale of the MD simulation the majority of the Tyr323 conformations are hydrogen bonded to Asp97 and positioned well into the substrate channel. However, a fraction of the conformations involves a disruption of the Tyr323-Asp97 hydrogen bond and a major displacement of Tyr323 from its position in the substrate channel. We therefore selected two snapshots to represent these possibilities. The first of these was taken at 72.53 ns where Tyr323 maintained an H-bond interaction with Asp97 in all three subunits. It represents the major conformational orientation observed from MD in this system (Figure 2a, left panel). The second snapshot was taken at 79.60 ns where the Tyr323 is displaced from the active site (Figure 2a   The energies of the relaxed conformers were assessed to determine whether there was any preference for one binding mode over another. Absolute energies between the conformers cannot be compared since the total number of atoms of the conformers is not conserved as they are taken from random time steps and a shell of solvent is used. Instead a relative energy was calculated from the difference in energy between relaxed top-hat and N-bound structures (∆E = E top-hat -E N-bound ) such that a positive ∆E indicates that N-bound is favoured while a negative value indicates top-hat is favoured. Figure 2 indicates that within the timescale of the MD simulation the majority of the Tyr323 conformations are hydrogen bonded to Asp97 and positioned well into the substrate channel. However, a fraction of the conformations involves a disruption of the Tyr323-Asp97 hydrogen bond and a major displacement of Tyr323 from its position in the substrate channel. We therefore selected two snapshots to represent these possibilities. The first of these was taken at 72.53 ns where Tyr323 maintained an H-bond interaction with Asp97 in all three subunits. It represents the major conformational orientation observed from MD in this system (Figure 2a, left panel). The second snapshot was taken at 79.60 ns where the Tyr323 is displaced from the active site (Figure 2a H-bond interaction with Asp97. As detailed in the Methods section optimizations were carried out separately for each monomeric unit, resulting in six independent starting structures. Relative energies and geometric descriptions of these are presented in Table S1. 1) Geometries of optimized Cu(II) and Cu(I) states starting from top-hat and N-bound orientations of NO 2 − In the case of Cu(II) starting from a top-hat orientation, all the relaxed structures retain the top-hat orientation, with both θ and ψ within 170-180 • . The nitrite oxygen atoms are roughly equidistant from Cu, with distances differing by~0.3 Å (Table S1). When reduced to Cu(I), the relaxed coordination changes to monodentate top-hat orientation with asymmetric Cu-O distances. θ lies within the range from 160 • to 180 • , with the majority of structures in the 170-180 • range, while ψ is around 145-165 • , a signature of monodentate NO 2 − in top-hat orientation. Here Cu is bound to O, hence N is necessarily displaced, resulting in the azimuthal angle of~150 • . In one conformation, θ is 159.1 • , at the lower limit of the observed range. Starting from N-bound orientations, there is more variation in the relaxed structures. All the final structures represent N-bound geometry with θ < 10 • . However, both symmetric (ψ < 10 • ) and L-shaped (ψ > 10 • ) orientations are observed. This is also reflected in the two Cu-O distances: for L-shaped N-bound structures, the difference between the two Cu-O distances are >0.3 Å. On reduction to Cu(I) we observe θ < 10 • as in Cu(II), but the N-bound structures are more symmetrical in nature, where ψ < 3 • and in many cases close to 0 • , with the oxygen atoms therefore equidistant from Cu. The range of θ and ψ observed for this system across all the 6 structures are given in Table 3. The angle ϕ is in the range of 80-90 • and six angles around Cu are in the range of 90-135 • , but in some cases, one of the N(NO 2 − )-T2Cu-N(His) angles is larger and fell in the range of 140-162 • (Table S1). Overall the structures can be classified as tetrahedrally distorted trigonal geometries.

2) Energetics of Cu(II) and Cu(I) states
In the case of Cu(II), for snapshots where the Tyr323 is engaged in a H-bond interaction with Asp97, the structures are mostly close in energy, and in one structure the bidentate orientation is favoured over the L-shaped NO 2 − orientation (Table S1). In the structure where Tyr323 is displaced from the T2Cu site, an extended H-bond network exists linking Tyr323 to Asp97 and the L-shaped N-bound structure has more favourable energy than the bidentate top-hat structure. The complex active site H-bond pattern as well as the interactions of the other active site residues is susceptible to the orientation of NO 2 − at the T2Cu site and likely contributes to the difference in energy observed in these structures. In Table 4 Molecules 2018, 23

D97p RpNiR
In contrast to the D97 system, protonation of Asp97 causes frequent transient disruption of the Tyr323-Asp97 H-bond (Figure 2a, middle panel), while in the timescale of the MD simulation no significant movement of Tyr323 from its position in the substrate channel is observed (Figure 2b, middle panel). Three representative snapshots were taken which differed by the number of monomeric units that featured the presence of a Tyr323-Asp97 H-bond: in all monomeric units (74.52 ns), in two of the three monomeric units (44.31 ns), and in one monomeric unit (48.40 ns) respectively. Optimization of these snapshots was carried out for all monomeric units, providing 9 independent conformations. Two structures in Cu(II) state resulted in orientations with very short Cu-O-Asp97 distance, which is not observed experimentally. These structures were discarded, resulting in 7 independent starting structures. Relative energies and geometric descriptions of all conformations are presented in Table S1.

−
In the case of Cu(II) starting from top-hat orientation, all the relaxed structures optimized in a bidentate top-hat orientation, with both θ and ψ within the range 170-180 • . The oxygens are approximately equidistant from Cu, with distances differing by~0.3 Å. When reduced to Cu(I), the relaxed orientation changes to monodentate top-hat. In the majority of the structures θ is in the range of 170-180 • . In two conformations, however, θ is 153.2 • and 145.0 • , respectively falling outside the range given in Table 2. Such cases fall in the borderline between the top-hat and side-on categories. The azimuthal angle ψ lies in the range 145-165 • , the same as observed for the D97 system, which is a signature for monodentate top-hat orientation.
Similarly to the D97 system, the Cu(II) structures starting from an N-bound orientation show more structural variation than its Cu(II) top-hat counterpart. All the final structures feature an N-bound geometry with θ < 10 • ; however both symmetric (ψ < 10 • ) and L-shaped (ψ > 10 • ) orientations are observed. This is also reflected in the two Cu-O distances: for the L-shaped N-bound structures the difference between the two Cu-O distances are >0.3 Å (Table S1). On reduction to Cu(I), all N-bound structures have θ < 20 • as shown in Table 2; with the majority within θ < 10 • . The L-shaped orientation is not observed, rather all of them fall in the category of symmetrical N-bound orientation with equidistant O-atoms and ψ < 6 • . In the majority of structures, the variation in ψ is even less; ψ < 3 • (Table S1). The distribution of θ and ψ observed for this system are also provided in Table 3. The angle ϕ is in the range of 78-90 • and six angles around Cu are in the range of 90-135 • but in a few cases one of the N(NO 2 − )-T2Cu-N(His) angles is larger and fell in the range of 140-145 • (Table S1). Like D97, the geometries can overall be classified as tetrahedrally distorted trigonal structures.

2) Energetics of Cu(II) and Cu(I) states
In the Cu(II) case the energy difference between individual top-hat and N-bound structures is usually less than 5 kcal/mol (see Table S1) which could arise from small changes in the active site H-bond network. There is no clear evidence of energetic differences arising from the orientation of the NO 2 − ion. Unlike the D97 system, there is a small trend towards symmetrical N-bound structures over monodentate structures. Table 4 shows the average ∆E for both Cu(I) and Cu(II) states. For Cu(II) the structures are effectively isoenergetic with a small standard deviation, whereas for Cu(I) the average ∆E suggests that N-bound is favoured, but with a large standard deviation showing there is significant variations over the snapshots.

D97N Mutant RpNiR
As with the D97p system, a mutation of Asp97 to Asn weakens the hydrogen bond between Asn97 and Tyr323 causing frequent transient disruption of the H-bond interaction (Figure 2a Table S1.  Table 2. The oxygen atoms are approximately equidistant from Cu, with distances differing by~0.3 Å. When reduced to Cu(I), the relaxed coordination changes to a monodentate top-hat orientation. θ ranges from 160 • to 180 • , with the majority of structures within 170-180 • . Unlike for D97 and D97p, in one conformation a side-on orientation of NO 2 − is observed (θ and ψ~88 • ). This is consistent with data for the two-domain CuNiR, whose T2Cu site lacks the linker residue [21]. Another conformation falls in the borderline of top-hat and side on with θ = 149.4 • . In most of the structures the value of ψ is around 145-165 • as observed for D97 and D97p systems and implies a monodentate top-hat orientation. There is another outlier, where ψ is 172.1 • and is of bidentate nature. Again, as with the D97 and D97p systems more structural variations were observed starting from the N-bound NO 2 − orientation. All the final Cu(II) structures feature an N-bound geometry with θ <10 • ; however symmetrical N-bound structures (ψ < 10 • ) are more abundant than the L-shaped N-bound (ψ > 10 • ) ones. In one N-bound case a penta-coordination around T2Cu is observed with symmetrical N-bound NO 2 − and a water. On reduction to Cu(I), the structures relaxed to symmetrical N-bound orientation with θ < 10 • , and in the majority of structures ψ < 3 • . For two structures ψ values of 4.3 • and 6.0 • were observed, but this is still well within the thresholds defined in Table 2. The distribution of these two geometrical parameters is provided in Table 3 (Table S1), hence, the overall structures are classified as tetrahedrally distorted trigonal structures.

2) Energetics of Cu(II) and Cu(I) states
In the Cu(II) case the individual energy differences between top-hat and N-bound structures are slightly higher than observed for the D97p systems; the maximum difference ranging up to 8 kcal/mol (see Table S1). The distribution of N-bound orientations in Cu(I) state slightly exceeds the top-hat orientations. In this system too, there is no clear correlation between the orientation of the NO 2 − at the active site and the energy variation. The average ∆E for both Cu(I) and Cu(II) states is provided in Table 4.
Here we see that for Cu(II) the average energy difference suggests that top-hat may be marginally favoured overall but with significant variation between snapshots as indicated by the large standard deviation. For Cu(I) the result is clearer cut with N-bound favoured over top-hat overall.

Discussion
The MD simulations in this work reveal an important consequence of the flexibility of the linker region in RpNiR. They suggest that Tyr323 does not rigidly block the substrate access channel to the active site, rather it is flexible and can be displaced by the natural water dynamics in the resting state of the enzyme, causing the channel to open, which would allow NO or NO 2 − to enter or leave the catalytic site. The MD corroborates the findings of Dong et al. for the D97N mutant [20], who observed displacement of Tyr323 for the NO-bound structure, resulting in an opening up of the channel. It is well known that sampling of the full configurational space is not possible in the time limit of typical all-atom MD simulations, and future work will use enhanced sampling methods to explore the flexibility of the Tyr323 residue and substrate entry to the T2Cu. A striking difference between the crystal structure of mutant RpNiR and two-domain CuNiRs is the mode of NO 2 − binding; in two-domain NiRs NO 2 − binds in a bidentate top-hat orientation to T2Cu whereas in the RpNiR mutant an N-bound orientation was observed [20]. In this work we have investigated the orientation of nitrite binding to T2Cu for both oxidized Cu(II) and reduced Cu(I) states of native and mutant RpNiR. To categorize the different modes by which NO 2 − can bind to the T2Cu we introduced two geometrical parameters, θ and ψ, and observed that the conformations obtained on QM/MM optimization of the nitrite bound structures clustered into particular ranges. Figure 4 shows representative structures of the four clusters of geometries as classified by θ and ψ, which we label bidentate top-hat, monodentate top-hat, L-shaped N-bound and symmetric N-bound. There were however a few structures which adopted other orientations: one Cu(I) conformer in D97N optimized to a side-on NO 2 − orientation, and one Cu(I) conformer in D97N and two Cu(I) conformers in D97p have θ values of approximately 150 • , which lies in-between top-hat and side-on orientation. The side-on orientation was seen in MSOX data for the reduced Cu(I) state for AcNiR [21], and hence could be a viable ligand orientation for binding to reduced RpNiR. The twist angle, ϕ, though originally used to define the type I Cu coordination geometry, can also be extrapolated to define the geometry around the T2Cu site [29]. In the case of RpNiR, for all the systems studied ϕ is in the range of 80-90 • , inclining towards trigonal structures. Considering the His-Cu-NO 2 − and His-Cu-His angles along with ϕ, the structures are best classified as tetrahedrally In our calculations on RpNiR, however, we only find clustering of the former category.
In the oxidized Cu(II) state the relative stabilities of the N-bound and top-hat states is not clear cut. However, for the reduced Cu(I) state there is a slight preference for the N-bound state in both the D97p and D97N systems and a clear preference for it in the D97 system (  There were however a few structures which adopted other orientations: one Cu(I) conformer in D97N optimized to a side-on NO2 − orientation, and one Cu(I) conformer in D97N and two Cu(I) conformers in D97p have θ values of approximately 150°, which lies in-between top-hat and side-on orientation. The side-on orientation was seen in MSOX data for the reduced Cu(I) state for AcNiR [21], and hence could be a viable ligand orientation for binding to reduced RpNiR.
The twist angle, φ, though originally used to define the type I Cu coordination geometry, can also be extrapolated to define the geometry around the T2Cu site [29]. In the case of RpNiR, for all the systems studied φ is in the range of 80-90°, inclining towards trigonal structures. Considering the His-Cu-NO2 − and His-Cu-His angles along with φ, the structures are best classified as tetrahedrally distorted trigonal structures. In the work of Karllot et al. two distinct clustering of φ with maximum N(His)-Cu-N(His) angle was observed in the reported crystal structures of all two-domain NiRs [29]. One had φ in the range of 82-89° and maximum N(His)-Cu-N(His) angles around 111-125° whereas in the other φ was in the range of 65-77° and the maximum N(His)-Cu-N(His) angles around 134-  Figure 5 shows the overlay of the X-ray structures of the nitrite-bound D97N mutant [20] and two simulated conformers, from the D97 and D97N systems, with similar L-shaped NO 2 − orientations bound to Cu(II). In both systems, the Tyr323 was oriented in the substrate access channel as in the native RpNiR crystal structure, rather than being displaced as observed experimentally for D97N. Apart from the difference in the Tyr323 position, the orientation of NO 2 − , the coordinated histidines and the other active site residues are all similar to the observed X-ray structure. In the crystal structure θ is 4.9 • and ψ is 45.2 • and the corresponding angles from the optimized structures are θ: 5.8 • , 3.5 • and ψ: 11.3 • , 16.3 • . Our findings indicate that the X-ray structure is consistent with the L-shaped N-bound NO 2 − orientations observed for the Cu(II) state in RpNiR and not the reduced Cu(I) state (SI, Figure S2), owing to the symmetrical N-bound orientation we observe in the reduced state. Such an L-shaped orientation is not primarily observed for two-domain CuNiRs in the Cu(II) state where a bidentate top-hat orientation of NO 2 − is mostly prevalent.
θ is 4.9° and ψ is 45.2° and the corresponding angles from the optimized structures are θ: 5.8°, 3.5° and ψ: 11.3°, 16.3°. Our findings indicate that the X-ray structure is consistent with the L-shaped Nbound NO2 − orientations observed for the Cu(II) state in RpNiR and not the reduced Cu(I) state (SI, Figure S2), owing to the symmetrical N-bound orientation we observe in the reduced state. Such an L-shaped orientation is not primarily observed for two-domain CuNiRs in the Cu(II) state where a bidentate top-hat orientation of NO2 − is mostly prevalent.

Materials and Methods
Nitrite binding to native RpNiR has not been structurally characterised experimentally to date and only very recently has nitrite binding in a mutant (D97N) been reported [20]. We used a native structure of RpNiR in its resting state as a template to study the mode of nitrite binding to the T2Cu site. The native crystal structure at 1.01 Å resolution [9] was taken from the PDB databank (PDB ID: 3ZIY). The structure was first 'cleaned' for MD simulations by removing all double occupancy mainchain and sidechain atoms and partial water molecules, keeping those with greater fractional occupancy or lower B-factors. The homotrimeric biological unit was generated by symmetry operations using the clean version. Hydrogen atoms were added and the protonation states of the titratable residues apart from those at the active site were adjusted to be consistent with neutral pH using the propKa module of the PDB2PQR suite [33], followed by visual inspection of the local sidechain environments. The protonation states of the catalytically important residues, Asp97 (D97) and His240, were treated carefully. The residue His240 was singly protonated as HSD but the residue Asp97, the main residue involved in proton transfer, was considered in both its deprotonated and protonated forms in the native state. An in silico mutation to obtain D97N was performed on the

Materials and Methods
Nitrite binding to native RpNiR has not been structurally characterised experimentally to date and only very recently has nitrite binding in a mutant (D97N) been reported [20]. We used a native structure of RpNiR in its resting state as a template to study the mode of nitrite binding to the T2Cu site. The native crystal structure at 1.01 Å resolution [9] was taken from the PDB databank (PDB ID: 3ZIY). The structure was first 'cleaned' for MD simulations by removing all double occupancy mainchain and sidechain atoms and partial water molecules, keeping those with greater fractional occupancy or lower B-factors. The homotrimeric biological unit was generated by symmetry operations using the clean version. Hydrogen atoms were added and the protonation states of the titratable residues apart from those at the active site were adjusted to be consistent with neutral pH using the propKa module of the PDB2PQR suite [33], followed by visual inspection of the local side-chain environments. The protonation states of the catalytically important residues, Asp97 (D97) and His240, were treated carefully. The residue His240 was singly protonated as HSD but the residue Asp97, the main residue involved in proton transfer, was considered in both its deprotonated and protonated forms in the native state. An in silico mutation to obtain D97N was performed on the native RpNiR structure, so that Tyr323 retained the same initial orientation in the substrate channel/T2Cu pocket in all three systems being investigated: D97, D97p and D97N.
The prepared structures corresponding to the three systems were then solvated with a 15 Å layer of TIP3P water [34]. The electroneutrality of the D97, D97p and D97N systems were maintained by adding 6, 3, and 3 sodium counterions, respectively. Explicit all-atom MD simulations were performed on these systems using NAMD 2.9 [35] with the CHARMM36 force field [36]. These simulations employed Langevin dynamics with periodic boundary conditions at 293 K. Long range electrostatics were treated by the Particle Mesh Ewald method. In the NPT simulations the pressure was maintained with the Langevin piston method. Both the systems were initially subjected to 5000 steps of conjugate gradient (CG) minimization to eliminate any unphysical contacts. Next, the water and ions were equilibrated using an NVT ensemble, keeping the protein fixed for 1 ns. This was followed by 5000 steps of CG minimization and 5 ns equilibration under an NPT ensemble, keeping the backbone harmonically restrained (5 kcal −1 mol −1 Å 2 ) and the coordination spheres of both the T1 and T2 sites (Cu(His) 2 (Met)(Cys − ) and Cu(His) 3 (H 2 O), respectively) constrained at their crystallographic coordinates. The simulation was continued for another 75 ns after removing the backbone restraints and allowing the coordinated water to move freely. The trajectories obtained from the MD simulations were analysed using VMD [37].
Snapshots were taken from the equilibrated MD trajectory as a starting point for QM/MM optimizations. The location of Tyr323 in the active site and the H-bond interaction of Tyr323 with Asp/Asn97 were used as the essential criteria for selecting the snapshots of both the native and mutant systems as detailed in the results section for each system. To obtain NO 2 − bound structures, these selected snapshots were overlaid with the analogous "top-hat" NO 2 − coordinates from the nitrite bound AcNiR structure [38]. The water(s) coordinated to Cu and any nearby water molecules that were in unphysically close steric interaction with the overlaid NO 2 − ion were then deleted. The nitrite coordinates were flipped to create analogous starting structures for the N-bound optimisations. In the absence of experimental data of the binding mode, this approach ensured that there was no bias in the initial starting structures towards either the top-hat or N-bound orientations. The forces acting on the monomers along the simulation diverge, thus, the conformational arrangement at the catalytic sites which are housed at the interphase of two monomeric units are independent of one another, providing three independent starting structures in each snapshot. All these structures were subjected to QM/MM optimization of the nitrite bound T2Cu site following the standard QM/MM partition scheme. The QM region consists of the three histidine residues coordinating T2Cu (His99, His134 and His289), the NO 2 − ligand, the three important active site residues Asp/Asn97, His240, and Ile242, the Tyr323 side chain, and any water molecules immediately H-bonded to Asp/Asn97, His240, Tyr323 and the NO 2 − ligand. The residues were cut at the neutral Cα-Cβ bond and were capped with H-atoms to satisfy the valency for the QM calculations. Atoms within 7 Å of the QM region remained unconstrained during QM/MM geometry optimizations while the remaining atoms were frozen. The QM/MM optimizations were carried out for both the oxidized Cu(II) and reduced Cu(I) states of Cu starting from two orientations of NO 2 − : the η 2− OO (bidentate) 'top-hat' and η 1 -N (monodentate) 'N-bound', to ensure that the potential energy surface was sufficiently explored. QM/MM calculations were performed with the Tcl-based version of the ChemShell computational chemistry environment [39,40], using ORCA [41] and DL_POLY [42] for the density functional theory and MM calculations respectively, and the DL-FIND module for geometry optimisations [43]. The electrostatic embedding scheme with charge shift correction was used to represent the surrounding MM partial charge distribution. The density functional B3LYP [44] with the DFT-D3 dispersion correction [45] was used for the QM atoms during geometry optimization and def2 − SVP basis sets were used for all QM atoms except Cu, which was treated with the def2 − TZVP basis set [46]. The CHARMM36 force field was used for the MM region.
Supplementary Materials: The supplementary materials are available online, Figure S1: Overlap of the NO 2 − bound X-ray structure of D97N and the frame at 79.6 ns of D97 MD in its resting state. Chain A, where the Tyr is displaced (Figure 3 in main text) is highlighted here. The coordination around the T2Cu site and active site residues Asp97 and Tyr323 illustrate the similarity in the position of Asp97 and displacement of Tyr323 in MD that structurally matches the mutant crystal structure; Figure S2: Overlap of the NO 2 − bound X-ray structure of D97N and two symmetrical N-bound conformers obtained from an initially N-bound NO 2 − in the Cu(I) state of the D97 and D97N systems; Table S1: Geometric parameters and energies for the optimized conformers obtained from the D97, D97p and D97N systems. The distances are Cu-N, Cu-O1 and Cu-O2, respectively, where O1 and O2 are the two equivalent O-atoms of NO 2 − .