Introduction

YiiP is a bacterial Zn²⁺/H⁺ antiporter and a well characterized representative of the Cation Diffusion Facilitator (CDF) superfamily. Members of this family play important roles in homeostasis of transition metal ions such as Zn²⁺, Mn²⁺, Co²⁺, and Fe²⁺ (Montanini et al., 2007). Mn²⁺ transporters from this superfamily are prevalent in plants, where this metal ion is essential for oxygen generation by photosystem II as well as for a variety of other enzymatic functions (Alejandro et al., 2020). Zn²⁺ CDF transporters are widespread in organisms from all kingdoms of life reflecting a large influence of this ion on cell biology. An estimated 10% of proteins employ Zn²⁺ either as a catalytic co-factor (e.g., carbonic anhydrase and cytochrome C oxidase) or as a structural element (e.g., Zn-finger transcription factors and steroid receptors) (Maret, 2013). In mammals, high concentrations are found in a variety of intracellular vesicles that are key players in the immune system, synaptic transmission, insulin trafficking, and fertilization (Kambe et al., 2015). Zn²⁺ figures in host-pathogen interactions, with the host attempting to either deprive or poison pathogens residing in endosomes (Lonergan and Skaar, 2019). Despite its prevalence, the bulk of intracellular Zn²⁺ is bound to protein, with vanishingly small (<10^-12 M) concentrations of free Zn²⁺ in the cytoplasm (Liang et al., 2016; Outten and O’Halloran, 2001).

Homeostasis is maintained by CDF transporters, which generally export Zn²⁺ from the cytoplasm, as well as members of the Zinc-regulated or Iron-regulated transport Proteins (ZIP), P-type ATPase and ATP-binding cassette (ABC) superfamilies (Yin et al., 2022).

YiiP serves as a model for studying structural and mechanistic properties of CDF transporters. Structures of homologs from both E. coli and S. oneidensis in different conformational states provide a framework for describing the archetypal alternating access mechanism (Lopez-Redondo et al., 2021; Lopez-Redondo et al., 2018; Lu et al., 2009; Lu and Fu, 2007). Molecular Dynamics has been used to study dynamics (Sala et al., 2019), to define the transport pathway (Sharma and Merz, 2022) and, together with in vitro binding assays, to characterize the K_d and pKa of Zn binding sites. Our understanding of CDF transporters is amplified by recent MD simulations of Znt2 (Golan et al., 2018; Golan et al., 2019), cryo-EM structures of Znt8 (Daniels et al., 2020; Xue et al., 2020; Zhang et al., 2023) and Znt7 (Han et al., 2022), as well as X-ray structures of isolated cytoplasmic domains from a number of different species (Udagedara et al., 2020). Together, these studies establish characteristics of Zn²⁺ binding sites, a conserved homodimeric architecture, and the nature of inward-facing (IF) and outward-facing (OF) states that presumably characterize all CDF transporters.

Despite this extensive work, there is uncertainty about key mechanistic questions that we have sought to address in the current work. Our primary goal was to obtain experimental evidence for the individual roles of the three Zn²⁺ binding sites found on each YiiP protomer (Cotrim et al., 2019). Site A is within the transmembrane domain (TMD) and features three Asp and one His residues; it is alternately exposed to the cytoplasm or to the periplasm in IF and OF states, respectively, and is thus directly responsible for transport of ions across the membrane. Although each molecule composing the dimer appears to have an independent transport pathway, there is unresolved potential for cooperativity. Site C within the C-terminal domain (CTD) is a binuclear site featuring four His and two Asp residues; one of these sites (C1) is well conserved whereas the other (C2) is not (Parsons et al., 2018). Site B is on the loop between transmembrane helices 2 and 3 (TM2-TM3) featuring two His and one Asp residue. Despite a lack of sequence homology for sites B and C2, the structure of Znt8 showed Zn ions at similar locations of the tertiary fold, but bound at different locations in the linear sequence, suggesting that these non-conserved, auxiliary sites may have functional significance (Xue et al., 2020).

Another goal was to address the role of protons in the transport process. There is ample evidence supporting Zn²⁺/H⁺ antiport. In vitro assays with Znt1, Znt2, ZitB and CzcD showed pH dependence of transport (Chao and Fu, 2004a; Cotrim et al., 2021; Guffanti et al., 2002; Shusterman et al., 2014). Binding studies with Znt8 (Zhang et al., 2023) as well as computational studies with Znt2 (Golan et al., 2019) and YiiP (Sharma and Merz, 2022) showed an interplay between Zn²⁺ and H⁺ binding at site A. However, there are discrepancies in stoichiometry, with work on ZitB (Chao and Fu, 2004a) and YiiP (Chao and Fu, 2004b) supporting a 1:1 exchange of Zn²⁺ and H⁺, whereas studies of CzcD (Guffanti et al., 2002) and Znt2 suggest 1:2 (Golan et al., 2019). There is also recent evidence for Zn²⁺/Ca²⁺ antiport by Znt1 in neurons (Gottesman et al., 2022).

For the current study, we have generated mutants of YiiP and have used cryo-EM and microscale thermophoresis (MST) to measure the pH dependence and structural effects of Zn²⁺ binding at each of the three sites. In addition, we have used molecular dynamics (MD) simulations together with the experimental MST data to deduce the pKa of residues at these sites and to address the stoichiometry of transport. From these data, we conclude that Zn²⁺ binding at site C is responsible for the integrity of the homodimer. Release of Zn²⁺ from site B triggers a conformational change in which the transport site A becomes occluded, suggesting a potential relay of Zn²⁺ between these two sites. Occlusion of only one protomer breaks the symmetry of the dimer and suggests that the transport process is not cooperative. In addition, binding affinity at site A displays a dramatic pH dependence, which can be explained by the protonation of 2 or 3 of the residues comprising this site. A corresponding Zn²⁺/H⁺ antiport stoichiometry of 1:2 or 1:3, depending on pH, is consistent with energetic coupling of Zn²⁺ export to the proton-motive force in a physiological setting.

Results

Mutants used to study the structural effects of Zn²⁺ binding

In order to study effects of Zn²⁺ binding at individual sites, we produced mutants of YiiP from Shewanella oneidensis. The D51A mutation was introduced to preclude binding at site A and D70A to preclude binding at site B. For site C we mutated Asp287, because it bridges the two Zn²⁺ ions at that site and we anticipated that D287A would therefore eliminate both ions.

The wild-type (WT) protein studied in previous work (Coudray et al., 2013; Lopez-Redondo et al., 2021; Lopez-Redondo et al., 2018) served as a positive control. For each of these constructs, protein was expressed in E. coli, solubilized in decyl-β-D-maltoside and purified by affinity chromatography. To ensure binding of Zn²⁺ rather than other metal ions picked up during growth and initial purification (e.g., Ni²⁺), we incubated the protein overnight with chelators and then added 250 μM Zn₂SO₄ prior to the final purification step. As in previous work (Lopez- Redondo et al., 2021), we used an Fab antibody fragment to produce a larger complex amenable to cryo-EM analysis: 165 kDa comprising the YiiP homodimer (65 kDa) bound to two Fab’s (50 kDa each).

Site A: D51A mutation

The structure of the D51A complex was solved at 3.6 Å resolution (Fig. 1b, Suppl. Fig. 1, Table 1) revealing an architecture very similar to untreated, WT YiiP (PDB code 7KZZ) that represents the IF, holo state (Lopez-Redondo et al., 2021). Specifically, YiiP formed a two-fold symmetric homodimer with Fab molecules bound near the C-terminus of each CTD. Density for Zn²⁺ ions is clearly visible at sites B and C (Suppl. Fig. 2e and f). At site A, however, this structure lacked density (Suppl. Fig. 2d), confirming that the D51A mutation effectively eliminated binding at this site.

Figure 1.

Table 1.

For comparison, we solved the structure from the WT construct that was similarly loaded with Zn²⁺. As expected, this 3.8 Å resolution structure (Fig. 1a, Suppl. Fig. 1) is indistinguishable from 7KZZ, which is also from the WT construct but not treated with chelators or explicitly loaded with Zn²⁺ (RMSD 1.07 Å for all 562 Cα atoms), Despite overall similarity, comparison with the D51A structure showed movements of TM1, TM4 and TM5 (Fig. 2a), which resulted in a notably elevated RMSD for this region: 4.7 Å based on the corresponding 200 Cα atoms compared to 0.69 Å for the 352 residues composing the rest of the YiiP dimer. The cytoplasmic end of TM5 had the largest differences. Specifically, a kink is introduced near His155 and the TM4-TM5 loop is disordered (Suppl. Fig. 2d), suggesting a flexibility of this region in the absence of Zn²⁺ at site A.

Figure 2.

Site B: D70A mutation

An initial structure of the D70A construct revealed a global conformational change reminiscent of the previously solved WT, apo state (PDB 7KZX) produced by removal of Zn²⁺ with chelators (Lopez-Redondo et al., 2021). The most salient feature of this conformation is a kink between TMD and CTD, which disrupts the overall two-fold symmetry of the complex.

During image processing, it became evident that there were in fact multiple conformations in the dataset from D70A, which ultimately generated two distinct structures at 4.0 and 3.9 Å resolution (Fig. 1c, Suppl. Fig. 3). Unlike the WT, apo structure (7KZX), both of these new structures revealed densities at sites A and C (Suppl. Fig. 2g,i,j, and l), consistent with Zn²⁺ binding to intact sites. However, the two structures differed both in the conformation of the TMD and in the orientation of the CTD. One of these structures, termed D70A_sym, has a TMD displaying local two-fold symmetry, in which both protomers adopt the IF state (Fig. 2b). Like the WT, apo structure, the TM2/TM3 loops, which harbor Asp70 and site B, are disordered. In the other structure, termed D70A_asym, one TMD adopts the IF state with a disordered TM2/TM3 loop, whereas the other TMD adopts a novel conformation in which the TM2/TM3 loop extends ∼17 Å away from the membrane to interact with H1 and β1 elements of the CTD from the opposing protomer (Fig. 2c, Suppl. Fig. 5b). In addition to reconfiguring this loop, there are movements in TM1,2,4,5 (Fig. 2d) that close off site A from the cytoplasm. This occlusion is documented by the considerably smaller radius of the cavity leading to this site (Suppl. Fig. 13).

Despite these changes, local dimeric elements from the D70A structures superimpose closely onto the corresponding elements from the WT, holo structure. In particular, alignment of CTD’s produces low RMSD’s of 0.82 Å and 1.08 Å relative to D70A_sym and D70A_asym, respectively (168 Cα atoms). The core helices mediating the TMD dimer interface (TM3/6) also align closely with WT with RMSD’s of 0.79 Å and 1.22 Å for D70A_sym and D70A_asym (106 Cα atoms, Fig. 2b, d). For D70A_sym, the entire TMD is generally consistent with the WT, holo structure (Fig. 2b), though slight displacement of TM1,4,5 and bending of the cytoplasmic end of TM2 lead to a somewhat elevated RMSD of 2.1 Å for 366 Cα atoms from both protomers after alignment based on TM3/6. For D70A_asym, the TMD in the IF conformation has a similar RMSD of 2.6 Å relative to WT, holo (183 Cα atoms from protomer B, Fig. 2d), whereas the occluded TMD has an elevated RMSD of 5.8 Å (197 Cα atoms from protomer A after aligning TM3/6) due to large differences in the cytoplasmic ends of TM1,2,4,5.

The observed structural effects of the D70A mutation are broadly consistent with MD simulations based on the WT structure. In particular, when site B was empty, the TM2/TM3 loop exhibited greatly increased mobility documented by >2-fold increase in root mean square fluctuations (RMSF, Fig. 3a); this result is consistent with disorder of this loop in D70A cryo-EM structures. The CTD also exhibited increased mobility as documented both by the broader distribution of the angles relative to the TMD (Fig. 3d) and by per residue RMSF, in particular when the TMD was used as reference for alignment (Fig. 3b), indicating that the CTD moved as a fairly rigid body relative to the TMD. In structures with Zn²⁺ bound at site B, the TM2/TM3 loop from one protomer is close to the linker between TM6 and the CTD of the opposing protomer (Fig. 3e,f, Suppl. Fig. 2b,e, and n). We therefore used the distance between Cα atoms of Asp72 in the TM2/TM3 loop and Arg210 in the TM6/CTD linker from the opposing chain as a collective variable for tracking CTD movement. The corresponding distance distributions (Fig. 3c) are broader and extend to larger values when site B is empty. In simulations of the D72A mutant with Zn²⁺ bound at site B there is a marked increase in CTD/TMD angle (Fig. 4d) as well as a more modest increase in the distance distribution (Fig. 4c). These results, together with structures of D70A_sym and WT, apo, suggest that a salt bridge between Asp72-Arg210 helps maintain the symmetry seen in the holo state. Release of Zn²⁺ from site B causes disordering of the TM2/TM3 loops leading to loss of the salt bridge and thus greater mobility of the CTD. However, the D70A_asym structure indicates that this disordering is transient, as the loop adopts the extended conformation and associates with different elements of the CTD (Fig. 2c). Interestingly, this novel change in the TM2/TM3 loop is coupled with movements of TM1,2,4,5 that close site A off from the cytoplasm (Fig. 3i-j, Suppl. Fig. 13). Mapping of electrostatic charge shows that whereas the access channel to site A in the IF protomer is negatively charged, positive charge dominates the cytoplasmic surface of the protomer in the occluded conformation (Fig. 3k,l).

Figure 3.

Figure 4.

Table 2.

Site C: D287A and D287A/H263A mutations

The D287A mutant was chosen to disrupt site C since it coordinates both Zn²⁺ ions at this site. However, the corresponding structure (3.7 Å resolution, Fig 1d, Suppl. Fig. 4) shows no structural changes relative to the WT structure (1.18 Ä RMSD for all 562 Cα atoms) and Zn²⁺- related densities are still clearly visible at site C (Suppl. Fig. 2o). We therefore introduced a second mutation, D287A/H263A, and the corresponding particles had an unusual appearance consistent with formation of higher order oligomers (Suppl. Fig. 4b). This observation suggested that site C may be responsible for the integrity of the homodimer. Indeed, SEC elution profiles of the YiiP-Fab complex show that the main peak from D287A/H263A is shifted relative to the other mutants (from 11.1 to 10.1 ml, c.f., Suppl. Figs. 1a,b, 3a, 4a,b). Elution profiles from complexes formed with D70A and D287A also have an earlier peak at ∼10 ml and, in both cases, image processing revealed a relatively small subpopulation of particles forming a dimer of dimers (Suppl. Fig. 3c & 4e). However, for D287A/H263A, this dimer of dimers represents the main peak and the only class of particles that could be refined to high resolution (3.5 Å) (Suppl. Fig. 4h). For D287A/H263A, higher order oligomers are apparent during 2D and 3D classification (Suppl. Fig. 4d & f) and are likely to explain the earlier elution peak at 8-9 ml.

The dimer of dimers from D287A/H263A is stabilized by an unprecedented domain swap in which CTD’s and TMD’s affiliate with different protomers (Fig. 2g). This domain swap is enabled by a structural transition of the linker between the TMD and the CTD. In particular, the normally unstructured loop between TM6 and the first helix of the CTD (H1) is reconfigured into a very long, continuous α-helix extending from the N-terminus of M6 (Trp178) to the C-terminus of H1 (Glu226). Site C itself is completely disordered and the CTD’s are splayed apart (Figs. 2g & Suppl. Fig. 5c): i.e., the distance between Cα atoms of Arg237 and Glu281 is 22.7 Å compared to ∼12 Å in structures with an intact site C (Lopez-Redondo et al., 2021). This domain swap was not observed in the dimer-of-dimers formed by D70A and D287A, where the CTD and site C were unperturbed and Fab molecules were uniquely responsible for inter-dimer interactions. Although it was not possible to refine a structure from the higher-order oligomers from D287A/H263A, due to heterogeneity and preferred orientation, they appeared to comprise a linear chain of molecules propagated via this domain swap. These surprising changes in the CTD indicate that site C is indeed crucial in maintaining the integrity of the native homo-dimer.

Despite the domain swap, however, TMD’s from the D287A/H263A complex are quite congruent with the WT structure (RMSD 1.35 Å for 346 Cα atoms excluding the TM2/TM3 loop) and density for Zn²⁺ ions is clearly visible at site A (Suppl. Fig. 2p). Weak density is visible for the TM2/TM3 loops, indicating that they adopt an extended configuration that disrupts site B and allows them to interact with the CTD; it does not appear that Zn²⁺ is bound at site B in this domain-swapped complex (Suppl. Figs. 2q and 5c,d).

Figure 5.

Zn²⁺ binding affinity

To assess the binding affinity of individual Zn²⁺ binding sites, we used MST to analyze a series of mutants designed to isolate the individual sites (Fig. 4 and Suppl. Fig. 6). Specifically, the triple mutant D70A/D287A/H263A was used to study site A, D51A/D287A/H263A to study site B, and D51A/D70A to study site C. To maintain accurate and reproducible Zn²⁺ concentrations ranging from nanomolar to micromolar, we used nitrilotriacetic acid (NTA, K_d = 14 nM) to buffer Zn²⁺ for sites with high affinity, and citrate (K_d = 12 μM) for sites with lower affinity. At pH 7, sites A and C displayed relatively high affinity (K_d = 16 and 33 nM, respectively, Table 2, Fig. 4), whereas site B had considerably lower affinity (K_d = 1.2 μM). Site C is binuclear, and the apparent affinity was reduced when further mutations were introduced to isolate individual C1 (D51A/D70A/H263A with K_d = 153 nM) and C2 (D51A/D70A/H234A with K_d = 223 nM) sites, suggesting cooperative binding of ions at C1 and C2 as implied by the coordination geometry of this site (Fig. 1d inset).

To explore the basis for coupling of Zn²⁺ transport to the proton-motive force, we measured Zn²⁺ binding affinity at pH values from 5.6 - 7.4. We found that affinity at site A changed by five orders of magnitude: K_d ranging from 1 nM at pH 7.4 - consistent with the cytosol - to 302 μM at pH 5.6 (Fig. 4a, Table 2). Site B had only modest pH dependence ranging from 1-16.6 μM, whereas Site C changed by two orders of magnitude from 0.033 - 6.7 μM. Cooperativity at site C is consistent with a high Hill coefficient at pH 7 (n=2.9), which fell below 1 at lower pH’s. However, Hill coefficients obtained for other sites were quite variable (Table 3), including those for sites C1 and C2, making conclusions about cooperativity inconclusive.

Table 3.

The domain swap seen in cryo-EM structures of D287A/H263A raises concern that site C mutations might affect affinity measured at other sites. However, in the absence of Fab, all of the mutants used for MST studies eluted from SEC columns at the same volume as WT (Suppl. Fig. 6a-d) indicating that they all adopted the native dimeric structure. This observation suggests that the domain-swapping seen in the structure of D287A/H263A was induced by the Fab molecules. This conclusion is supported by direct comparison of preparations with site C mutations before and after addition of Fab (Suppl. Fig. 6 h & i). In the absence of Fab, they both elute from SEC as a single peak at the expected volume of 12 ml. In the presence of Fab, there is a shift of the main peak consistent with formation of a higher molecular weight complex and, additionally, appearance of a second peak. In the case of D287A-Fab, the main peak at 11.1 ml is consistent with the YiiP-Fab complex seen for the other mutants (Suppl. Fig. 1 & 3) and the second peak at 9.9 ml is consistent with the dimer of native homodimers seen during image processing. In the case of D287A/H263A, the main peak in the presence of Fab is at 10.1 ml and thus consistent with a dimer of domain-swapped dimers and the secondary peak at 8.6 ml is consistent with higher order oligomers seen during cryo-EM processing (Suppl. Fig. 4). In addition, we performed MST analysis with the D70A/D287A mutation as an alternative for measuring affinity at site A. The cryo-EM structure of D287A shows that the native homodimeric assembly is retained and although Zn²⁺ was observed at site C, we expect it will bind at lower affinity and thus not overlap with higher affinity binding at site A. Indeed, this construct produced Kd= ∼1 nM at pH 7 and ∼6 μM at pH 6, which is consistent with results from D70A/D287A/H263A (Table 2), thus confirming the pH dependence of site A.

We also used MD simulations in conjunction with the experimental MST data to address pH dependence of Zn²⁺ binding and to evaluate the contribution of individual residues. To start, we estimated pKa values for each titratable residue at site A (Asp47, Asp51, His155 and Asp159) and site B (Asp70, His73 and His 77) using constant pH MD (CpHMD) simulations with pH replica exchange (Huang et al., 2021) in the absence of Zn²⁺. It was not possible obtain data for site C because the Zn²⁺-free CTD proved unstable in CpHMD simulations. These simulations generated probabilities for each protonation state as a function of pH; these probabilities were fit with the Hill-Langmuir equation to derive per-residue pK_a’s (Suppl. Figs. 7-9, Table 4). For site A residues, His155 had the highest predicted pK_a (8.0), followed by Asp159 (4.9), Asp47 (3.8) and Asp51 (2.9). Predicted populations of the individual microstates (s0-s15, Suppl. Figs. 10a & b) indicated that two residues are likely protonated at pH 5.5 with His155 getting protonated first followed by Asp159. Although these simulations were conducted with the IF conformation, we do not expect significant differences for the OF state, given the similarity of site A seen in the X- ray structure (Lu et al., 2009) of the OF state (RMSD of 0.49 Å for 17 Cα atoms on TM2 and TM5 and 1.01 Å for all atoms composing the four residues that coordinate the Zn²⁺ ion). For site B, His73 and His77 have almost identical pK_a’s of 8.1, and the dominant state at neutral pH consists of both His73 and His77 protonated (Suppl. Figs. 11a & b).

Table 4.

We went on to use a thermodynamic model incorporating both protonated and Zn²⁺- bound states to assess the stoichiometry of transport (Fig. 4h). We used this model to generate Zn²⁺ binding curves based on pK_a’s of individual residues and binding free energies for Zn² derived from the experimental MST data (Table 2). Initially, pK_a’s of each residue in the absence of Zn²⁺ were taken from the Hill-Langmuir fits to the CpHMD data (Table 4, Suppl. Fig. 9).

However, this CpHMD-based model did not agree with the experimental data (Figs. 4d & e), likely due to inaccuracy in microstate pK_a values. We therefore employed a novel modeling method based on the inverse Multibind approach (see Methods) to iteratively refine the CpHMD pK_a’s and the MST Zn²⁺ binding free energies using the MST data points as a target for Monte Carlo (MC) minimization. This analysis, referred to henceforth as MST inference, allowed us to estimate probabilities of all the microstates as a function of pH and Zn²⁺ concentration (Suppl. Figs. 10 & 11). After refinement of parameters by MST inference, the predicted Zn²⁺ binding curves closely reproduced the observed pH-dependence of Zn²⁺ binding at both sites A and B (Figs. 4f & g). The dramatically improved agreement is a consequence of substantial changes in the pK_a’s of individual residues based on the MST inference algorithm (Table 4).

A more detailed look at the population of microstates predicted by MST inference allows us to explore the interplay of Zn²⁺ and H⁺ binding under physiological conditions. For site A in the absence of Zn²⁺, the dominant state at cytosolic pH of 7.5 consists of 2 protons bound by His155 and Asp159 (Suppl. Figs. 10c & d). Based on the microscopic pK_a values, we calculated the “coupling energy” (Ullmann, 2003) and thus deduced that these two residues form a coupled system (see Methods). As a result, protonation of His155 and Asp159 is highly cooperative and produces a very steep binding curve (Suppl. Fig. 10e); the singly protonated state is strongly suppressed and the effective pK_a is 7.8 for both residues (Suppl. Fig 10d). At the lower pH of the periplasm, a third proton is recruited by Asp51, whose pK_a increased to 6.5 during refinement by MST inference (Table 4). Regardless of the pH, as Zn²⁺ binds at site A, all the residues are deprotonated (Suppl. Fig. 10f), supporting the idea that Zn²⁺ transport can be coupled with a proton gradient across the plasma membrane. For site B, the dominant state at cytosolic pH in the absence of Zn²⁺ has one protonated residue: either His73 or His77 (Suppl. Fig. 11c & d). Although the microstate pK_a’s are equivalent, these site B residues behave quite differently from Site A residues and generate an anti-cooperative binding curve with effective pK_a’s of 5.6 and 12.4 (Suppl. Figs. 11e & f). Zn²⁺ binding at cytosolic pH displaces all protons (Suppl. Fig. 11g). Even though protons may co-exist with Zn²⁺ at lower pH’s, this is not a physiologically relevant condition given that site B is uniquely exposed to the cytoplasm.

Discussion

In this study, we have focused on the properties of individual Zn²⁺ binding sites of YiiP from S. oneidensis. Previous studies have shown that this protein forms a homodimer in the IF state with Zn²⁺ ions constitutively bound at three sites (Lopez-Redondo et al., 2018). Treatment with EDTA to remove Zn²⁺ from all three sites induced a conformational change to produce an occluded state, a necessary precursor to the OF state (Lopez-Redondo et al., 2021). For the current work, we generated mutations at each site to measure their binding affinities and to assess their respective roles in generating conformational change. We found that site A, considered to be the transport site in the middle of the TMD, has nanomolar affinity at a cytoplasmic pH of 7.4, which is reduced 10⁵-fold at the more acidic pH of 5.6. Release of Zn²⁺ from site A induces only modest movement of the peripheral transmembrane helices, but has no global effect on the architecture of the homodimer. Site B, on the short loop between TM2 and TM3, has much lower affinity in the micromolar range with relatively little pH dependence. When site B lacks Zn²⁺, the loop becomes disordered and the protein undergoes a global conformational change leading to an occluded conformation in one of the protomers. Site C, a binuclear site in the CTD that engages residues from both protomers, has relatively high affinity and intermediate pH dependence. When site C was fully disrupted with a double mutation, the homodimer became destabilized such that the Fab molecules used for cryo-EM imaging induced a domain swap of the CTD leading to higher order oligomers.

Elements stabilizing the homodimer

Dimerization appears to be a universal feature of Cation Diffusion Facilitators. All structures so far show similar elements contributing to a common dimer interface (Suppl. Fig. 12). In particular, a conserved salt bridge exists at the cytoplasmic surface of the membrane (Lys79-Asp209 for soYiiP); the cytoplasmic end of the CTD generally includes one or more Zn²⁺ sites, and extensive hydrophobic interactions involving TM3 helices mediate TMD interactions (Lopez-Redondo et al., 2021). In early work, the salt bridge was proposed to act as a fulcrum for alternating, scissor-like movements in TMD and CTD driven by relay of Zn²⁺ between sites A and C (Lu et al., 2009). Furthermore, the original X-ray structure of ecYiiP featured a V-shaped architecture in which TMD’s were completely disengaged (Lu et al., 2009; Lu and Fu, 2007).

However, subsequent work showed that stabilization of TMD interactions by cysteine crosslinking did not inhibit transport activity (Lopez-Redondo et al., 2018) and compact dimeric TMD interfaces were observed in cryo-EM structures of soYiiP, Znt8 and Znt7 (Han et al., 2022; Lopez-Redondo et al., 2018; Xue et al., 2020). Indeed, the TMD interface is a remarkably conserved feature in our new structures of D70A and D287A/H263A mutants, despite large- scale conformational changes involving the TM2/TM3 loop and domain swapping of the CTD. We therefore conclude it to be unlikely that the TMD’s undergo large-scale scissoring movements as part of the transport cycle.

For isolated CTD’s, metal ion dependent scissor-like movements have been well documented using truncated constructs from a variety of species (Cherezov et al., 2008; Udagedara et al., 2020; Zeytuni et al., 2014). We observed similar movements in full-length soYiiP: i.e., the CTD’s move apart when Zn²⁺ is removed from site C either by chelation or mutation. However, other elements must contribute, given that isolated CTD’s dimerize even in the absence of metal ions and Zn²⁺ ions were not observed at all in the CTD of Znt7 (Han et al., 2022). Nevertheless, the Fab-induced domain swap seen in D287A/H263A indicates that, for soYiiP, CTD interactions are weakened in the absence of Zn²⁺. Conversely, it seems likely that dimer stabilization, via TMD and salt-bridge interactions, lead to enhanced Zn²⁺ binding in the CTD, as evidenced by the considerably higher Zn²⁺ affinity at site C of full-length soYiiP (K_d of 33 nM at pH 7) compared to those measured from isolated CTD’s (K_d in micromolar range (Udagedara et al., 2020; Zeytuni et al., 2014)). In any case, the relatively high affinity of site C at pH 7 suggests that this site will remain occupied and will serve to ensure dimer stability under physiological conditions.

Role of the site B

Structural effects of the D70A mutation support our previous hypothesis that site B is responsible for inducing a global conformational change. The most conspicuous change is the kink between TMD and CTD, which breaks the global two-fold symmetry of structures that retain Zn²⁺ binding at site B, in which local two-fold axes of CTD and TMD are aligned (Suppl. Fig. 12). Although stable, kinked conformations were not achieved over the time-scale of MD simulations, increased movement of the CTD relative to the TMD is seen both in the previously reported apo state simulations (Lopez-Redondo et al., 2021) and with an empty site B (Fig. 3), indicating a role of the Zn²⁺-bound TM2/TM3 loop in stabilizing the position of the CTD. Indeed, van der Waals interactions are observed between the structured loop from one protomer and the TM6- CTD linker of the opposing protomer on both sides of the dimer (Fig. 3e,f), and MD simulations indicate that a salt bridge between Asp72 and Arg210 may reinforce these interactions (Fig. 3d, c.f., site B in Suppl. Fig. 2).

Both cryo-EM structures and MD simulations show that the TM2/TM3 loop becomes disordered when Zn²⁺ is released from site B and this is a likely explanation for increased mobility of the CTD. Two structures were obtained from the D70A mutation in which the CTD becomes progressively more tilted and twisted. In one structure, the TMD’s are symmetrical and retain the IF conformation (D70A_sym), whereas in the other structure with more extreme CTD movement (D70A_asym), the TMD transitions to an occluded state. We speculate that these structures represent a sequence of conformational change leading to the occluded state.

Although both TM2/TM3 loops are disordered in the D70A_sym structure, one of these loops reforms into an extended conformation in the occluded protomer of the D70A_asym structure (Suppl. Fig. 5). This change results in a novel interaction between the loop and the CTD, which may drive the displacement and bending of membrane helices thus leading to occlusion of the transport site A (Suppl. Fig. 13). This interaction would also prevent the TM2/TM3 loop from rebinding Zn²⁺ until the protein returned to the IF state. On the other side of the dimer, the CTD has lost contact with the TMD, allowing the TM2/TM3 loop to remain disordered and the TMD to remain in the IF state. Interestingly, the TM2/TM3 loops are also extended in the D287A/H263A domain-swapped structure (Suppl. Fig. 5), but in this case, the interactions of this loop with the reconfigured CTD’s are quite different and the TMD’s remain in the IF state, indicating that a specific structural constraint is required to instigate occlusion.

Occlusion of the transport site A is a result of bending of TM5 and tilting of TM1 and 4, thus narrowing the gap between these helices and closing off access to the cytosol. These movements also generate a dramatic difference in the electrostatic surface at the cytoplasmic side of the TMD. In the IF state, the open cavity leading toward site A is negatively charged, as also shown for Znt7 (Han et al., 2022), thus serving to attract Zn²⁺ toward this site. After transition to the occluded state, not only is the cavity closed (Suppl. Fig. 13), but the surface becomes positively charged thus repelling Zn²⁺ (Fig. 3k,l). Previous work has identified a “hydrophobic gate” consisting of two residues at the cytoplasmic ends of TM5 and TM6: Leu154 and Leu199 in soYiiP (Gupta et al., 2014). Our previous comparison of holo and apo states showed that these two residues do come closer together in the occluded state (Lopez-Redondo et al., 2021) and homologous residues in Znt7 are close in the OF state and separated in the IF state (Han et al., 2022). Interaction of Leu154 and Leu199 is also seen in the occluded protomer in the D70A_asym structure (Fig. 2d), but these residues are at the periphery of the cavity leading to site A. This observation suggests that interaction of these residues may play an important role in stabilizing the occluded state, but that bending at the cytoplasmic end of TM5 (Val148-Ala151) and tilting of TM1 and TM4 may be more directly responsible for blocking access to the transport site.

For soYiiP, the IF state appears to be a low-energy, ground state when Zn²⁺ is present at all three sites. Indeed, the IF state has been seen for WT protein both in lipid-based helical crystals (Lopez-Redondo et al., 2018) and detergent micelles (Lopez-Redondo et al., 2021), as well as in D287A and D287A/H263A mutant structures, all with RMSD’s <1.5 Å. Removal of Zn²⁺ from site B can be viewed as a source of energy for inducing conformational change and likely leading to transport. Interestingly, the occlusion only occurs in one protomer, suggesting a lack of synchrony across the dimer axis. A mixture of conformational states was also observed for Znt8, in which one protomer adopted the OF state while the other was in the IF state. These observations imply that a given protomer may be able to undergo transport while the other remains inactive, as illustrated in Fig. 5.

The proposed sequence of conformational change is consistent with the notion that Zn²⁺ is initially recruited to site B and that transfer to site A would act as a trigger for the first step of the transport cycle. Such a two-step mechanism might be necessary if Zn²⁺ were delivered to YiiP as a chelate – e.g., with metallothionein or glutathione – that would not have direct access to site A. In the case of Cu⁺, metallochaperones are thought to play a role in delivery (Robinson and Winge, 2010), but so far there is very limited evidence for Zn²⁺ metallochaperones (Chandrangsu et al., 2019). Nevertheless, essentially all of the Zn²⁺ in the cytoplasm exists in chelated form with vanishingly small concentrations of free Zn²⁺ available for direct binding (Choi and Koh, 1998; Outten and O’Halloran, 2001). The accessibility of site B makes it a plausible acceptor as an initial recruitment site. We postulate that site B would only be available for binding Zn²⁺ in the IF apo state (Fig. 5d). Once bound, two features would facilitate transfer to site A: the drastic difference in affinity between sites A and B would generate a large free- energy gradient and the negative charge of the cavity leading to site A, which lies almost directly below site B in the WT, holo structures, would direct the ion toward the transport site. Although the amino acid sequence of site B is not conserved amongst CDF transporters, a Zn²⁺ ion was observed in an analogous juxtamembrane position on Znt8 (Xue et al., 2020) and many CDF transporters have histidine-rich loops, typically between TM4 and TM5, which might fulfill a similar role in recruiting Zn²⁺ and orchestrating conformational change. Indeed, for Znt7 this His- rich loop was shown to bind two Zn²⁺ ions leading these authors also to speculate that it plays a role in recruitment and shuttling to the transport site (Han et al., 2022).

Stoichiometry of transport and energy coupling

The proton-motive force is generally postulated to drive the antiport mechanism utilized by CDF transporters. In the case of YiiP, this antiport involves exchange of Zn²⁺ from the cytoplasm for H⁺ in the periplasm and the stoichiometry of this process is key to energy coupling. In previous work, the Zn²⁺:stoichiometry has been measured to be 1:1 for ZitB (Chao and Fu, 2004a) and deduced by cellular or computational studies to be 1:2 for CzcD (Guffanti et al., 2002) and Znt2 (Golan et al., 2019), respectively. In addition, isothermal calorimetry was used to deduce a binding stoichiometry for Cd²⁺ of 1:1 for ecYiiP (Chao and Fu, 2004b). For the current work, we refined CpHMD simulation data with experimental MST data using the inverse Multibind approach (Kenney and Beckstein, 2023) to infer the prevalent protonation states for site A and thus to address stoichiometry. Our analysis is consistent with Zn²⁺ binding to a fully unprotonated site A regardless of the pH. After Zn²⁺ release, the transport site becomes either doubly or triply protonated, depending on the pH. Due to strong coupling between His155 and Asp159, the singly protonated state is essentially not seen. Thus, this model is consistent with a stoichiometry of at least 1:2, possibly 1:3 in more acidic environments (pH < 6).

Our analysis predicts significant shifts in pK_a for Zn-binding residues. Such shifts reflect strong Coulomb interactions due to both the low dielectric of the membrane and the clustering of several titratable residues at the binding site. Effects of the membrane environment have been studied both by experimental (Gayen et al., 2016; Isom et al., 2010; Morrison et al., 2015) and computational (Henderson et al., 2020; Panahi and Brooks, 2015) methods, showing pK_a values up to 9 for acidic residues. Coupling between nearby residues has also been shown to shift pK_a values (Yue et al., 2017), with estimated pK_a values of 3 and 11 for a di-aspartyl pH sensor in a pH-sensitive calcium channel serving as a dramatic example (Chang et al., 2014). Together with the current work, these examples illustrate the importance of tuning the local environment to harness the energy of the proton motive force.

The Zn^2+:H⁺ stoichiometry has crucial physiological consequences for the energetics of transport, which are governed by the Nernst equation describing electrochemical potential. The overall transport cycle can be described as

with the free energy calculated as

where m is the number of protons, n is the number of Zn²⁺ ions, z is the net charge for the reaction and V_m is the membrane potential. Based on this equation, the higher stoichiometries of 1:2 or 1:3 provide the cell with increased leverage over Zn²⁺ transport. In the absence of membrane potential, for example, a 10-fold proton gradient (ΔpH of 1) can produce Zn²⁺ gradients at equilibrium of up to 10² or 10³, respectively, instead of simply 10 for a 1:1 stoichiometry (more generally, 10^m/n). Membrane potential represents an important component of the proton motive force that will influence electrogenic transport, where z is non-zero. In particular, a stoichiometry of 1:1 would produce net positive charge transfer out of the cell and thus would require working against this membrane potential. Given a potential of -80 mV, relatively normal in E. coli (Felle et al., 1980), a 1:1 stoichiometry would produce an unfavorable energy term that would overcome a chemical gradient of one pH unit. In contrast, a 1:2 stoichiometry would be electroneutral and would therefore be unaffected by membrane potential thus generating a Zn²⁺ gradient of 10² as described above. Finally, a stoichiometry of 1:3 would benefit both from the membrane potential as well as from the increased number of protons to theoretically generate a gradient of 10^4.4.

These calculations are based on equilibrium thermodynamics, but the pH dependence of Zn²⁺ binding also has implications for the kinetics of transport. In particular, the dramatic difference in K_d at cytoplasmic pH of 7.4 (1 nM) compared to more acidic pH’s in the periplasm (e.g., 6 μM at pH 6) implies that the on-rate of Zn²⁺ is dominant in the cytoplasm but that the off- rate is enhanced in the periplasm, assuming that the binding site in the OF state has similar pH dependence. This assumption seems plausible given the similarity of transport site geometry of YiiP in the OF state (Lu et al., 2009). Thus, although the final equilibrium concentrations of Zn²⁺ are not affected by the pH dependence of the transport site, the rate of transport and the rate of equilibration will be greatly enhanced (Tanford, 1983).

YiiP has three aspartates and one histidine at site A, whereas many CDF transporters including the mammalian Znt’s substitute Asp47 for a second histidine. Previous work has focused on changes in ion specificity associated with this change (Hoch et al., 2012), but it also seems likely to affect stoichiometry. At a minimum, two histidines would make the 1:2 stoichiometry almost certain and a 1:3 stoichiometry highly plausible. The resulting electrogenicity could be desirable for eukaryotic cells in which these transporters operate predominantly in intracellular organelles such as insulin secreting granules, synaptosomes, golgi or zinc-o-somes that are responsible for the dramatic Zn²⁺ sparks during oocyte fertilization (Chu, 2018; Hara et al., 2017). Since the pH gradient of these organelles is modest, additional driving force from the membrane potential might be key in producing the high internal concentrations that are sometimes required for function.

MATERIALS and METHODS

Protein expression and purification

YiiP was expressed in E. coli (BL21(DE3)-CodonPlus-RIPL) from a pET vector that included an N-terminal decahistidine tag. Cells were grown in LB media supplemented with 30 μg/ml kanamycin at 37°C until they reached an OD₆₀₀ of 0.8. After cooling the media to 20°C, expression was induced by addition of 0.5 mM isopropyl-β-D-thiogalactoside followed by overnight incubation at 20°C. Cells were harvested by centrifugation at 4,000xg for 1 h, resuspended in lysis buffer (20mM HEPES, pH 7.5, 100mM NaCl, 10% glycerol, and 500 μM tris(2-carboxyethyl)phosphine) - 100 ml of buffer per 20 g of cells - and then lysed with a high- pressure homogenizer (Emulsiflex-C3; Avestin, Inc. Ottawa Canada). Protein was extracted from the membrane by adding 1.5 g dodecyl-β-D-maltoside per 100 ml of cell lysate followed by 2h of incubation at 4°C. Insoluble material was removed by centrifugation at 100,000xg for 30 min. The supernatant was loaded onto a Ni-NTA affinity column pre-equilibrated in buffer A (20 mM HEPES, pH 7.5, 100 mM NaCl, 10% glycerol, and 0.05% dodecyl-β-D-maltoside). The column was washed by addition of buffer A supplemented with 20 mM imidazole and protein was then eluted using a gradient of imidazole ranging from 20 to 500 mM. Peak fractions were combined, supplemented with tobacco etch virus (TEV) protease (1:10 weight ratio of TEV:YiiP) to cleave the decahistidine tag, and dialyzed overnight at 4°C against buffer A. TEV protease was removed by loading the dialysate onto an Ni-NTA column and collecting the flow-through fractions. After concentration, a final purification was done with a Superdex 200 size-exclusion chromatography (SEC) column (GE Healthcare, Chicago, Illinois) equilibrated with SEC buffer (20 mM HEPES, pH 7.5,150 mM NaCl, 0.2% n-decyl-β-D-maltoside, and 1 mM tris(2- carboxyethyl)phosphine).

Fab selection, modification, expression, and purification has been described previously (Lopez-Redondo et al., 2021). Briefly, the construct designated Fab2r was expressed in E. coli strain 55244 from a freshly prepared transformation. Cells were cultured for ∼24 h at 30°C with constitutive expression behind an innately leaky T4 promoter. Cell pellets were harvested by centrifugation, lysed with a high-pressure homogenizer and Fab was purified with a 5-ml HiTrap Protein G HP column (GE Healthcare). Pooled fractions were dialyzed against sodium carbonate buffer (pH 5.0) overnight at 4°C and further purified with Resource-S cation exchange column (GE Healthcare). Finally, pure Fab protein was pooled and dialyzed against SEC buffer.

Cryo-EM sample preparation and structural analyses

To ensure that YiiP was fully loaded with Zn²⁺, the purified protein was initially incubated with metal ion chelators (0.5 mM EDTA and 0.5 mM N,N,N′,N′-tetrakis(2-pyridinylmethyl)-1,2- ethanediamine [TPEN]) for 24hrs at 4°C to remove ions that co-purified with the protein, as seen in previous work (Lopez-Redondo et al., 2021). The sample was then dialyzed against YiiP SEC buffer containing 0.5 mM EDTA for 10 h to eliminate chelator-metal complexes. Finally, YiiP was loaded with Zn²⁺ by dialysis (10 h with four buffer exchanges) against SEC buffer supplemented with 0.25 mM ZnSO₄. YiiP was then incubated with Fab2r at a 1:1 molar ratio for 1 h at 20°C to form the YiiP-Fab2r complex. This complex was purified by SEC using a Superdex 200 column equilibrated with SEC buffer supplemented with 0.25 mM ZnSO₄. Peak fractions at 3-5 mg/ml were used immediately for preparation of cryo-EM samples. Specifically, 3–4 μl were added to glow-discharged grids (C-Flat 1.2/1.3-4Cu-50; Protochips, Inc.) that were blotted under 100% humidity at 4°C and plunged into liquid ethane using a Vitrobot (Thermo Fisher Scientific, Inc. Bridgewater NJ).

The YiiP-Fab2r complex was imaged with a Titan Krios G3i electron microscope (Thermo Fisher Scientific, Inc) equipped a Bioquantum energy filter and K2 or K3 direct electron detector (Gatan, Inc. Pleasanton CA) with a pixel size of ∼1 Å and a total dose of ∼50 electrons/Ų. Micrographs containing crystalline ice, excessive contamination, or imaging artifacts were removed and the resulting micrographs were imported into cryoSPARC v2.15 (Punjani et al., 2017) for analysis. Particles were picked based on templates generated in previous work and an initial set of particles were subjected to two rounds of 2-D classification to remove false positives. The resulting particles were then subjected to successive rounds of ab initio reconstruction with C1 symmetry and a resolution cutoff starting at 12 Å and declining to 8 Å; at each step, the best of two output classes was carried forward to the next round. The resulting particles were used for heterogeneous refinement against two or three reference structures derived from the ab initio jobs, still with C1 symmetry. A final selection of particles was then used for non-uniform refinement using both C1 and C2 symmetry. Postprocessing steps included calculation of local resolution and evaluation of 3-D variability (Punjani and Fleet, 2020). For the D70A construct, RELION was used to process a class of particles corresponding to a dimer of dimers (Suppl. Fig. 3). After exporting these particles from cryoSPARC, steps of 2D classification and ab initio reconstruction were repeated in RELION, followed by 3D refinement of ∼117,000 dimer-of-dimer particles with C2 symmetry. A mask encapsulating one dimer was created with Chimera (Pettersen et al., 2004) and symmetry expansion followed by signal subtraction was used to generate a new set with ∼234,000 dimeric particles. This expanded particle set was imported back into cryoSPARC and combined with the class of isolated dimer particles. This combined particle set was used for hetero-refinement to segregate dimer complexes with symmetrical and asymmetrical TMD’s. These segregated particle sets were then used for final non-uniform refinement.

For model building of symmetrical structures from WT, D287A, D51A, we started with the deposited coordinates from wild type (PDB, 7KZZ). This model was docked as a rigid body to the map and adjusted manually in Coot (Emsley et al., 2010) as a starting structure for refinement using PHENIX (Adams et al., 2010). For D70A and D287A/H263A mutants, the CTD and the TMD from either 7KZZ or 7KZX were separately docked to the map. After adjustment and crude rebuilding in Coot, these models were submitted to NAMDINATOR for Molecular Dynamics Flexible Fitting (Kidmose et al., 2019). The resulting models were used as starting points for PHENIX refinement. This refinement consisted of multiple rounds that alternated between real-space refinement and manual adjustment using Coot. Atomic models were then displayed using CHIMERA (Pettersen et al., 2004) and the PyMOL Molecular Graphics System (Schrödinger, LLC, New York NY).

For characterization of cavities running from site A to the cytoplasm, we used Caver Analyst 2.0 Beta (Jurcik et al., 2018). The starting point was defined set near D51 in models for WT and D70A_asym and we selected the shortest tunnel that ran to the cytoplasm. Settings used for the calculations were as follows. Approximation: 12, Minimum probe radius: 0.9, Clustering threshold: 3.5, Shell depth: 2.0, Shell radius: 3.0. Zn²⁺ ions were excluded from the calculation.

MST experiments

YiiP mutants were labelled with Alexa fluor 488 fluorescent dye (Invitrogen life Technologies, Carlsbad CA) by adding 2.5 μL of dye from a 16 mM stock solution in DMSO to 200 μL of protein at 1-2 mg/ml in SEC buffer supplemented with 10% glycerol at pH 7. The reactive group of this dye was N-hydroxysuccinimide, which at pH 7 preferentially targets the N- terminus of the polypeptide. Although labeling of lysine residues is also possible, we believe this was minimal due to the low labeling stoichiometries of ∼1:1 used for our experiments. Zn²⁺ was removed by adding 0.5 mM EDTA and 0.5 mM TPEN followed by overnight incubation at 4°C. Excess dye and chelated Zn²⁺ were removed using several cycles of dilution with SEC buffer and concentration with a 50 kDa cutoff concentrator (AMICON, Millipore Sigma, Burlington MA). For titration, Zn²⁺ was buffered either by 0.5 mM sodium citrate or 0.2 mM NTA; the total amount of added ZnSO₄ was varied to achieve the desired concentration of free Zn²⁺, according to the program MAXCHELATOR (Bers et al., 1994). Protein concentration during the titration varied from 8 to 100 nM. After 1:1 mixing of protein and Zn²⁺ solutions, samples were centrifuged for 5 min at 2000xg, then loaded into standard treated capillaries for measurement with a Monolith NT.115 MST instrument (NanoTemper Technologies, South San Francisco, CA); measurements were taken at 37 °C with LED power ranges from 20-60% and medium MST power. Data from three independent titrations were analyzed with the MO.Affinity Analysis software v2.3 using MST on-time of 15 s. For determination of K_d, data was fitted with a curve based on the law of mass action,

where F(C_Zn) is the fraction bound, C_Zn is concentration of free zinc, C_P is the concentration of YiiP protein, F_b and F_u refer to the normalized fluorescence in the bound and unbound state, and K_d is the affinity constant. Alternatively, the data was fit with the Hill equation to assess cooperativity,

where EC₅₀ is the half-maximal effective concentration (akin to K_d) and n is the Hill coefficient.

Overview of MD simulations

Two types of all-atom, explicit solvent MD simulations were conducted to assess the effects of proton and zinc binding on the structure of the YiiP dimer in a lipid bilayer. We used fixed charge equilibrium MD simulations (using a force field with unchanging parameters) to investigate the interactions between Zn²⁺ ions and YiiP. We also carried out replica exchange constant pH MD (CpHMD) simulations to calculate the pK_a values of all titratable residues and assess the microscopic protonation states of the ion binding sites in the absence of Zn²⁺. Unlike CpHMD simulations, fixed charge simulations do not allow for dynamic protonation or deprotonation of residues. Despite this limitation, they are our preferred approach to study protein-zinc ion interactions and conformational changes because they can be run more efficiently than CpHMD and were thus used to quantify the effect of zinc ions on the structure of YiiP (see Table 5 for a summary of these simulations). The zinc ions in the system were simulated with the non-bonded dummy model described in our previous work (Lopez-Redondo et al., 2021); parameter files are available as package 2934 in the Ligandbook repository (Domański et al., 2017) (https://ligandbook.org/package/2934). Briefly, parameters governing this model were refined based on experimental hydration free energy, ion-oxygen distance and coordination number of the water in the first hydration shell, ultimately reproducing these values with errors of 1% and 0.3%, respectively. The model was further validated using simulations of known zinc-binding proteins, in which stability of the protein and the geometry of the binding complexes were well maintained, though coordination distances were slightly longer than experimental values, especially for sulfur atoms from cysteine residues (Lopez-Redondo et al., 2021). Given that YiiP employs aspartate or histidine residues for Zn²⁺ coordination, the model is well-suited for the current work.

Table 5.

Fixed charge equilibrium MD simulations

MD simulations in the holo (with Zn²⁺ in A, B, and C sites) and apo (no Zn²⁺) state were taken from our previous study (Lopez-Redondo et al., 2021). To investigate the influence of the binding site B, we generated an “empty site B” structure by removing the Zn²⁺ ions from site B in both protomers of PDB ID 5VRF, the same starting structure used previously. In order to study the influence of the salt-bridge D72-R210, simulations were also conducted after applying the D72A mutation to the holo structure (PDB ID 5VRF). In these D72A simulations, Zn²⁺ ions were bound in all sites (A, B, C1, C2), i.e., D72A simulations differed from the holo simulations only in the mutation.

The membrane-protein systems were modeled by embedding the YiiP dimer into a 4:1 palmitoyloleoylphosphatidylethanolamine:palmitoyloleoylphosphatidylglycerol (POPE:POPG) lipid bilayer, which approximates the composition of the plasma membrane from E. coli (Raetz, 1986), and solvating the system with water as well as sodium and chloride ions corresponding to a concentration of 100 mM using CHARMM-GUI v1.7 (Jo et al., 2008; Jo et al., 2009; Lee et al., 2016). We used GROMACS 2021.1 (Abraham et al., 2015) with the CHARMM36 force field, the CMAP correction for proteins (MacKerell et al., 1998; Mackerell et al., 2004), CHARMM36 lipids (Klauda et al., 2010), and the CHARMM TIP3P water model. The empty site B system contained 117,512 atoms, the apo system contained 117,394 atoms, and the D72A system contained 117,418 atoms. All systems were constructed in a hexagonal simulation cell with initial dimensions 101 Å × 101 Å × 135 Å. Default protonation states of all ionizable residues were used based on the experimental pH of 7. The neutral HSD tautomer (proton on the N_δ) was selected to model all histidines except H73 and H155, which were modeled with HSE (proton on the N_ε) based on their orientation relative to Zn²⁺ ions in the cryo-EM structure PDB ID 5VRF.

The systems were first energy minimized and underwent a 3.75-ns six-stage equilibration procedure with position restraints on protein and lipids, following the CHARMM- GUI protocol (Jo et al., 2008). Three copies of 1-μs production simulations were carried out with empty site B, starting from the same initial system conformation but with different initial velocities. All simulations were performed under periodic boundary conditions at constant temperature (T = 303.15 K) and pressure (P = 1 bar). The velocity rescaling thermostat (Bussi et al., 2007) was used to maintain the temperature with a time constant of 1 ps and separate temperature-coupling groups for protein, lipids, and solvent. A semi-isotropic pressure coupling scheme was implemented using the Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with a time constant of 5 ps, a compressibility of 4.6 × 10⁻⁵ bar⁻¹. Long-range electrostatics were calculated with the smooth particle mesh Ewald method (Essmann et al., 1995) under tinfoil boundary conditions with an initial cutoff of 1.2 nm, which was optimized during the simulation, and interactions beyond the cutoff were calculated in reciprocal space with a fast-Fourier transform on a grid with spacing 0.12-nm and fourth-order spline interpolation. The van der Waals interactions were switched smoothly to 0 between 1.0 nm and 1.2 nm, and the interactions were shifted over the whole range and reduced to 0 at the cutoff. The Verlet neighbor list was updated dynamically by GROMACS for optimized performance with a buffer tolerance of 0.005 kJ/mol/ps. Bonds to hydrogen atoms were treated as rigid holonomic constraints with the P-LINCS algorithm (Hess, 2008) with an expansion order of four and two LINCS iterations; alternatively, SETTLE (Miyamoto and Kollman, 1992) was used for water molecules. The classical equations of motion were integrated with the leapfrog algorithm with a time step of 2 fs.

CpHMD simulations with pH replica exchange

In conventional MD simulations, the protonation states of titratable groups in the system are fixed. To investigate the role of protons in the Zn²⁺ transport of YiiP, we performed membrane-enabled hybrid-solvent continuous constant pH MD simulations (Huang et al., 2021). The currently available implementation of the membrane-enabled CpHMD method does not take into account the direct effect of ions on the titration of nearby residues so we only ran CpHMD simulations for the apo system.

CpHMD simulations were initialized with the apo model (PDB ID 5VRF with Zn²⁺ ions removed). An initial 250-ns production simulation was run for the CpHMD apo system to fully relax the membrane, using the same GROMACS 2021.1 protocol described above for fixed charge equilibrium simulations.

The membrane-enabled hybrid-solvent continuous CpHMD simulations were performed using the CHARMM program version c42a2 (Brooks et al., 2009) with the PHMD module (Khandogin and Brooks, 2005; Lee et al., 2004) and the pH replica-exchange (REPDSTR) module (Wallace and Shen, 2011). The conformations of YiiP were sampled using conventional all-atom simulations with the CHARMM22/CMAP all-atom force field (MacKerell et al., 1998; Mackerell et al., 2004), the CHARMM36 lipid force field (Klauda et al., 2010), and the CHARMM modified TIP3P water model. The titration coordinates were propagated using the membrane GBSW implicit-solvent model (Im et al., 2003) with the GBSW input radii for the protein taken from Chen et al. (Chen et al., 2006). Based on the average distance between the C2 atoms of the lipids in the cytoplasmic- and periplasmic-facing leaflets, the thickness of the implicit bilayer was set to 40 Å with a switching distance of 5 Å for the transition between the low dielectric slab and bulk solvent. The implicit membrane was excluded by two cylinders with a radius of 14 Å placed at the center of mass of each protomer of YiiP. The radius was selected to maximize the coverage of the interior of the protein with minimal overlapping of the implicit membrane.

The final snapshot of the equilibrium simulation was used as the initial structure for the CpHMD simulation. Dummy hydrogen atoms were added to the carboxylate groups of acidic residues following the documentation of the PHMD module (Wallace and Shen, 2011) in CHARMM (Brooks et al., 2009) using the HBUILD facility. The system was then equilibrated with energy minimization using 50 steps of steepest descent followed by 50 steps of adopted basis Newton-Raphson algorithms and CpHMD at pH 7 for 1 ns, whereby the harmonic restraints on the protein heavy atoms were reduced from 1 kcal·mol⁻¹Å⁻¹ to zero.

The production simulation was then performed using hybrid-solvent CpHMD with the pH replica-exchange protocol (Huang et al., 2021; Wallace and Shen, 2011), using 30 replicas with pH ranging from 1.5 to 11.5. The specific pH conditions were 1.5, 1.75, 2, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.25, 4.5, 4.75, 5, 5.25, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 8.75, 9, 9.5, 10, 10.25, 10.5, 11, 11.5, chosen to ensure that the exchange rate between the nearby replicas was higher than 0.2. Each replica was simulated under periodic boundary conditions at constant temperature (T = 303.15 K), pressure (P = 1 bar), and specified pH. A modified Hoover thermostat method (Hoover, 1985) was used to maintain the temperature, while pressure was maintained using the Langevin piston coupling method (Feller et al., 1995) with a piston mass of 2500 amu. Long-range electrostatics were evaluated with the particle mesh Ewald method (Darden et al., 1993) with a real-space cutoff of 1.2 nm, and interactions beyond the cutoff were calculated in reciprocal space with a fast-Fourier transform on a grid with 0.09-nm spacing and sixth-order spline interpolation. The Lennard–Jones forces were switched smoothly to 0 between 1.0 and 1.2 nm, and the potential was shifted over the whole range and reduced to 0 at the cutoff. Bonds to hydrogen atoms were constrained with the SHAKE algorithm to allow a 2 fs time step. To avoid a spike in potential energy due to a lack of solvent relaxation (Wallace and Shen, 2011), a GBSW calculation was executed every 5 MD steps to update the titration coordinates. An attempt to exchange adjacent pH replicas was made every 1000 MD steps (corresponding to 2 ps). Each replica simulation lasted 12 ns for a total aggregate sampling time of 360 ns. Many replicas exchanged across a large fraction of available pH space (Supplementary Figure 7), indicating that 12 ns per replica were sufficient for sampling the degrees of freedom near the protonation sites.

Analysis of fixed charge MD simulations

Analysis of the trajectories was carried out with Python scripts based on MDAnalysis (Gowers et al., 2016). RMSDs of C_α atoms of the whole protein, TMD, and CTD were calculated using the qcprot algorithm (Liu et al., 2010) after optimally superimposing the structure on the same C_α atoms of the cryo-EM structure. Similarly, root mean square fluctuation (RMSF) of C_α atoms of the whole protein, TMD, and CTD were calculated. To assess the relative motion between the two domains, a CTD-TMD rotation angle was calculated from the rotation matrix which minimized the RMSD of the CTD after superimposing the protein on the TMD domain of the reference structure (PDB ID 5VRF), as in our previous work (Lopez-Redondo et al., 2021). The existence of the salt-bridge Asp72-Arg210 was quantified using the shortest distance between O_δ atoms of Asp72 and hydrogen atoms on the side chain of Arg210.

Analysis of CpHMD simulations

The titration coordinates were extracted from CpHMD output files as time series S(t) (Supplementary Figure 8) with the CpHMD-Analysis scripts (https://github.com/Hendejac/CpHMD-Analysis). The primary purpose of these simulations is to obtain microscopic pK_a values for titratable residues. We first describe the conventional approach to obtain per-residue pK_a values using the heuristic generalized Hill equation and then in the next section demonstrate an alternative inference approach based directly on statistical mechanics.

The deprotonation fraction S of a titratable site was calculated from the titration time series (Supplementary Figure 8) as the number of trajectory frames N of the residue in the deprotonated and protonated states as

where the site is identified as deprotonated when the CpHMD titration coordinate λ is greater than 0.8 and protonated when the titration coordinate is less than 0.2. Individual residue pK_a’s were obtained from a fit to the generalized Henderson-Hasselbalch equation (Hill equation) (Henderson et al., 2020; Huang et al., 2016; Huang et al., 2021) for the deprotonated fraction S as a function of pH,

where n is the Hill coefficient, which represents the slope of the response curve (Supplementary Figure 9). The mean pK_a of a residue was calculated for each pair of equivalent residues in protomer A and protomer B because these residues were sampled independently in the CpHMD simulation. The statistical error was estimated as the absolute difference between the pK_a’s of protomer A and B residues and their mean value.

The CpHMD simulations track the protonation of every single titratable residue and thus provide detailed microscopic information on the exact protonation state for each binding site.

The distribution of these microscopic states named S0 to S15 for site A (Supplementary Figure 10, inset table) and S0 to S7 for site B (Supplementary Figure 11, inset table) forms the basis for our alternative calculation of microscopic pK_a values. In order to obtain these CpHMD protonation state distributions, we treated titration coordinate data S(t) from protomer A and B as independent; thus, by concatenating them we effectively doubled our sampling. Using the same criterion for bound/unbound protons and the definition of microstates for site A (Supplementary Figure 10) and site B (Supplementary Figure 11) we generated separate microstate time series for sites A and B. These time series were then histogrammed to derive the distributions of the microscopic protonation states of site A (Supplementary Figure 10a) and B (Supplementary Figure 11a).

Calculation of microscopic, state-dependent pKa values with the Multibind method

The Multibind method generates thermodynamically consistent models for systems of coupled reactions by using a maximum likelihood approach to combine kinetic or thermodynamic measurements from different sources (Kenney and Beckstein, 2023). Given free energy differences between states (either from simulations or experiments), Multibind generates a complete set of free energies for all states in the form of a potential graph that obeys path- independence of free energies and detailed balance while being maximally consistent with the input data. From this set of free energies, all macroscopic thermodynamic observables can be calculated. Without such an approach, thermodynamically inconsistent models arise due to random errors in the input measurements. For Zn²⁺ and proton binding, the free energy differences in the potential graph are calculated as functions of the external parameters, namely the free Zn²⁺ concentration [X] and the pH, and with pK_a and Zn²⁺ standard state binding free energy values as input. The binding free energy of a reaction A + X ⇌ A: X is calculated as

where (k_B is Boltzmann’s constant), K_D is the dissociation constant, and the binding free energy is written as a sum of ΔG⁰_bind, the binding free energy at the standard state concentration c₀=1 M, and a term depending on the ligand concentration [X]. For proton binding, the free energy difference is expressed equivalently as

where pK_a is the acid dissociation constant.

We used the Python implementation of Multibind (https://github.com/Becksteinlab/multibind) to construct thermodynamic models for binding sites A and B, with the assumption that all binding sites across both protomers are independent.

Each model describes the transitions between all possible states as either binding of a proton or binding of a Zn²⁺ ion. With four titratable residues, site A has M=2⁴=16 Zn²⁺-free protonation microstates (states S0 to S15 in Supplementary Figure 10, inset table) and 16 Zn²⁺-bound protonation microstates for a model with 32 states in total (Fig. 4h). Site B has three titratable residues and thus M=8 for a model with 16 states in total (states S0 to S7 in Supplementary Figure 11, inset table). Because only free energy differences between states are measurable, we can arbitrarily specify one state as a reference with zero free energy. Here we chose the 0- proton unbound state S0 as the reference state with G₀=0, but calculated observables are independent of this particular choice. The probability of observing the system in the microstate i is

where G_i is the free energy of any state i and z is the partition function

which contains the partition functions of the unbound and bound states. We note that the microstate probabilities P_i can also be directly obtained from the CpHMD simulations, which enables us to use Multibind as an inverse method to infer microscopic pK_a from CpHMD data (without resorting to the generalized Hill equation fit), as described below.

The fraction of YiiP protomers that have a Zn²⁺ ion bound to the binding site in question, the bound fraction, is

The bound fraction can be obtained experimentally from MST measurements.

With the ability to calculate microstate probability or bound fraction at specific pH or Zn²⁺ concentration, we devised an inverse approach to infer microscopic, state-dependent pK_a’s and standard binding free energies (ΔG⁰_bind) from either CpHMD data (to avoid using the Hill equation with the deprotonated fractions, which are aggregated over multiple states and may therefore mask coupling between residues) or from experimental MST measurements.

The inverse Multibind approach generates successive microscopic models with a Monte Carlo (MC) scheme and compares the calculated observables to the target values. Observables correspond to microstate probability distributions P_i(pH) for CpHMD or bound fraction <X> as a function of pH and Zn²⁺ concentration for MST. Once convergence is reached, the microscopic model contains the inferred microscopic state-dependent pK_a values and Zn²⁺ binding free energies. We implemented a simple finite-temperature MC algorithm that minimizes the RMSD between the given and computed macroscopic features. The model was initialized with pK_a values from CpHMD (obtained with the Hill equation) or ΔG⁰ binding free energies calculated from the experimental K_d from the MST analysis. For each MC step, the target observable was calculated for a range of pH values and Zn²⁺ concentrations. For as long as the RMSD between the computed and the target observable remained above a cutoff (see Table 6), a new set of pK_a or ΔG⁰ was generated by adding a random value drawn from a uniform distribution ranging from –0.2 to 0.2 pK_a units (for pK_a) or –0.2 to 0.2 k_BT (for binding free energies) to all microscopic parameters. A new set of changes was accepted with probability (Metropolis criterion)

where ΔRMSD is the difference between the RMSDs between the target and the computed observable for the old set of pK_a’s or ΔG⁰ ’s and the new set

The fictional temperature in the Metropolis criterion, r, was set to 0.0001 to be on the order of typical per-step changes. The RMSD cutoffs (Table 6) were chosen from initial runs using a zero-temperature MC algorithm, which accepted a new set of changes only when ΔRMSD < 0. These initial runs were stopped after 2000 steps and after no changes had been accepted in the last 100 or more steps. The zero-temperature MC ensured that the final set generated the smallest RMSD, while the finite-temperature MC used in the production runs allowed larger parameter space to be sampled. Fifty runs were performed for each target.

Table 6.

Averaged pK_a values and corresponding standard deviations were calculated from the 50 independent MC runs. Although these standard variations of our parameter estimates weakly depend on the somewhat arbitrary choice of r, we chose to report them instead of a rigorous statistical error (which is unavailable with the current approach) to provide a sense of variability of the estimates.

For site B, the CpHMD simulations suggested that H73 and H77 behaved identically and had almost the same microscopic pK_a values. We calculated the “coupling energy” (Ullmann, 2003) between H73 and H77,

(expressed in pK_a units, i.e. energy divided by k_BT ln10), from the microscopic CpHMD state pK_a’s (see Supplementary Data). The coupling energy measures how the protonation of a specific residue (e.g. H73) depends on the protonation of another residue (here, H77). The coupling energy is the difference in pK_a when H77 remains deprotonated while H73 binds a proton (S0→S2, see Supplementary Figure 11) compared to the situation when H77 is already protonated (S1→S3). The H73-H77 coupling energy for the CpHMD simulation was W=+1.0 for D70 deprotonated and W=+0.2 for D70 protonated. Because W>0, binding of a proton to one of the histidines decreases the pK_a of the other, thus decreasing proton binding via anti- cooperative coupling. Preliminary calculations with the MST data, however, showed that small initial differences in the starting values could lead to large pK_a shifts in H73 and H77, with the directions determined by the initial ordering of pK_a values. Because we had no specific evidence that the two histidines should behave differently, we considered the behavior of the initial unconstrained MST inference calculations to be problematic. This problem indicated that the inverse approach can be sensitive to initial values and that there may not always be sufficient target data to constrain the microscopic model. We therefore treated the two residues as symmetrical, i.e., both should be behaving in the same way, and imposed a constraint in the Multibind approach so that both were assigned the same pK_a at each step, thus effectively imposing the coupling observed in the CpHMD simulations. The interaction energies between H73 and H77 for the MST-inference data are W=+6.2 (D70 deprotonated) and W=–0.70 (D70 protonated). The dominant state is the one with D70 deprotonated due to the low pK_a of D70 and hence the overall behavior of H73 and H77 remains strongly anti-cooperative in the MST- inference model.

For Site A, the initial CpHMD simulations did not indicate any degeneracy in pK_a values such as the one for H73/H77 and so no additional constraints were applied for any of the Site A MST-inference calculations. We also calculated the coupling energy for H155 and D159 in site A from the microscopic pK_a values (see Supplementary Data) for the CpHMD and the MST- inference models for (1) no protons bound to D51 and D47 (W_MST=–5.05, W_CpHMD=+1.26), (2) D51 protonated (W_MST=–7.39, W_CpHMD=+0.59), (3) D51 protonated (W_MST=–5.26, W_CpHMD=+1.02), and (4) both D47 and D51 protonated (W_MST=–2.39, W_CpHMD=+1.79).Thus, refinement against the experimental MST data changes the model based on the CpHMD results from anti- cooperative binding (W>0) to strong cooperative binding (W<0).

The inverse Multibind approach has some limitations when applied to binding curves such as the MST data. The binding isotherms, even when covering a range of pH values, do not contain enough data to determine the microscopic model fully. Therefore, it was necessary to use the CpHMD pK_a values to initialize the model instead of arbitrary starting values. Thus, in the current implementation, the MST inference approach should be viewed as a refinement process for the simulation-derived parameters, restricted or guided by the experimental data.

The resulting thermodynamic model allows the calculation of state probabilities (and any other properties) at pH values outside the experimental range but it should be noted that these calculated quantities are extrapolations that may not be accurate due to the lack of experimental data to constrain the model; for instance, conformational changes may occur that could drastically alter the interactions but these changes would not have been captured in a pure binding model such as the one in Fig. 4h.

An effective pK_a under the Zn²⁺-free condition was estimated for each residue by fitting the deprotonation fractions (generated from the microscopic Multibind model) to the Hill- Langmuir equation. The deprotonation fraction of a specific residue was calculated by summing up the probabilities of microscopic states where the residue was deprotonated. The Zn²⁺ concentration was set to 10^-20 M to approximate the Zn²⁺-free condition. For the symmetrized residues H73 and H77, the Hill-Langmuir equation did not produce a satisfactory fit. In order to compute macroscopic pK_a’s for such a coupled system, we followed previous work (Henderson et al., 2020; Ullmann, 2003) and fitted the total number of protons bound to the two residues to the so-called “coupled titration model”

where N_prot is the total number of protons, S₁ and S₂ are the deprotonation fractions of the two residues, and pK₁ and pK₂ are the two coupled pK_a’s describing the binding of the first proton and the second proton to the coupled titrating sites. These two pK_a’s can be interpreted as effective pK_a values of two uncoupled residues.

Acknowledgements

We are grateful for help with initial MST experiments from Brian Kloss and Renato Bruni in the Center on Membrane Protein Production and Analysis at the New York Structural Biology Center, which receives funding from NIH grant P41GM116799. Electron microscopy was performed at the Cryo-EM core facility at NYU Langone Health. The authors acknowledge High Performance Computing core facility at NYU Langone Health as well as Research Computing at Arizona State University for providing computational and storage resources that have contributed to the research results reported within this paper. We also thank Bjørn Panella Pederson for critical reading of the manuscript.

Funding for this work was provided by NIH grants R01GM125081 (D.L. Stokes) and R35GM144109 (D.L. Stokes). A portion of this research was also supported by NIH grant U24GM129547 and performed at the PNCC at OHSU and accessed through EMSL (grid.436923.9), a DOE Office of Science User Facility sponsored by the Office of Biological and Environmental Research. MD simulations were performed using PSC Bridges at the Pittsburgh Supercomputing Center (allocation TG-MCB130177), a resource of the Extreme Science and Engineering Discovery Environment, which is supported by National Science Foundation grant ACI-1548562.

Competing Interests

None.