The voltage sensor (VS) domain in Hv1 proton channels mediates a voltage-dependent and H+-selective “aqueous” conductance (GAQ) that is potently modulated by extracellular Zn2+. Although two conserved His residues are required for Zn2+ effects on GAQ gating, the atomic structure of the Zn2+ coordination site and mechanism by which extracellular Zn2+ stabilizes a closed-state conformation remain unknown. Here we use His mutagenesis to identify residues that increase Zn2+ potency and are therefore likely to participate in first solvation shell interactions with Zn2+. Experimental Zn2+-mapping data were then used to constrain the structure of a new resting-state Hv1 model (Hv1 F). Molecular dynamics (MD) simulations show how protein and water atoms directly contribute to octahedral Zn2+ coordination spheres in Zn2+-bound and -unbound Hv1 F models. During MD simulations, we observed correlated movements of Zn2+-interacting side chains and residues in a highly conserved intracellular Coulombic network (ICN) that contains highly conserved Arg “gating charges” in S4 as well as acidic “counter-charges” in S2 and S3 and is known to control VS activation, suggesting that occupancy of the extracellular Zn2+ site is conformationally coupled to reorganization of the ICN. To test this hypothesis, we neutralized an ICN Glu residue (E153) and show that in addition to shifting GAQ activation to more negative voltages, E153A also decreases Zn2+ potency. We speculate that extracellular gating-modifier toxins and other ligands may use a generally similar long-range conformational coupling mechanism to modulate VS activation in related ion channel proteins.

Divalent metal ions such as Zn2+ commonly function as protein cofactors by directly interacting with oxygen, sulfur, and nitrogen atoms in His, Cys, Glu, and Asp side chains (Stote and Karplus, 1995; Dudev et al., 2003; Laitaoja et al., 2013). Increasing the concentration of extracellular ZnCl2 causes the apparent open probability (POPEN)–voltage relation for the intrinsic (i.e., “aqueous”) H+-selective conductance (GAQ) in native and expressed Hv1 channels to shift positively (Cherny and DeCoursey, 1999; Ramsey et al., 2006; Sasaki et al., 2006; Musset et al., 2010; Takeshita et al., 2014; Qiu et al., 2016). The effect of ZnCl2 to modulate Hv1 gating is inferred to mean that zinc(II), which is a divalent cation (Zn2+) in aqueous solution, stabilizes a closed-channel conformation and therefore attenuates H+ current amplitude (Cherny and DeCoursey, 1999; Ramsey et al., 2006; Sasaki et al., 2006; Musset et al., 2010; Takeshita et al., 2014; Qiu et al., 2016). The relatively high concentration of Zn2+ in human seminal fluid is hypothesized to prevent spermatozoon activation in the male reproductive tract in vivo by preventing Hv1 channel opening, indicating that Zn2+ modulation of GAQ gating is physiologically meaningful (Lishko and Kirichok, 2010). Although recent studies have provided important insights into the mechanism of Zn2+ effects on Hv1 gating (Ramsey et al., 2006; Takeshita et al., 2014; Qiu et al., 2016), essential mechanistic details remain unknown.

The effects of Zn2+ and other divalent metal cations on GAQ gating are concentration dependent (IC50 ≈ 0.2–2 µM) and sensitive to changes in pHO (apparent pKa ≈ 6.4), suggesting that one or more His residues are required for Zn2+ coordination (Meech and Thomas, 1987; Byerly and Suen, 1989; Cherny and DeCoursey, 1999; Ramsey et al., 2006; Sasaki et al., 2006; Lishko et al., 2010; Musset et al., 2010). The availability of human, mouse, and Ciona intestinalis HVCN1 cDNAs encoding bona fide voltage-gated H+ channel proteins (Ramsey et al., 2006; Sasaki et al., 2006) led to the identification of His residues in the predicted S2 helical segment (H140/H2.40) and S3–S4 linker (H193/H3.71) in the human Hv1 voltage sensor (VS) domain that are necessary for Zn2+ effects on GAQ gating (Fig. 1 A; Ramsey et al., 2006). The H140A-H193A double mutant is essentially insensitive to submillimolar [Zn2+], but single H140A and H193A mutations decrease Zn2+ potency only 10- to 30-fold (Ramsey et al., 2006). Hv1 orthologues from C. intestinalis and Coccolithus pelagicus also lack H3.71 (Fig. 1 A) and exhibit similar Zn2+ sensitivity to human Hv1 H193A (Ramsey et al., 2006; Sasaki et al., 2006; Taylor et al., 2011; Qiu et al., 2016), indicating that H3.71 is not required for the effect of Zn2+ on GAQ gating. Electron density attributed to Zn2+ is observed near E1.58, D1.62, and H2.40 in the anomalous difference map of a mouse Hv1-based protein chimera (mHv1cc) x-ray structure solved to 3.45 Å (Takeshita et al., 2014), but the position of H3.71 is not resolved, and the low overall resolution precludes unambiguous determination of Zn2+ coordination sphere geometry or identification of liganded protein and water atoms.

Although E1.58 is close to Zn2+ in the mHv1cc x-ray structure, single Ser mutations at E115/E1.58 or D119/D1.62 do not measurably alter Zn2+ sensitivity (Takeshita et al., 2014). Mutation of another phylogenetically conserved acidic side chain (D/E3.61) that is predicted to be near E1.58 (Fig. 1 D) alters Zn2+ sensitivity in C. intestinalis Hv1. D1.51 is required for exquisite H+ selectivity in Hv1, but the large effects of mutations at this site on GAQ gating and channel expression (Rada and Leto, 2008; Ramsey et al., 2010; Qiu et al., 2016) could confound the interpretation of D1.51 mutant effects on Zn2+ potency. The differential sensitivity of H+ currents and fluorescence signals in Alexa Fluor 488–maleimide–labeled C. intestinalis Hv1 S242C (S3.70C) mutant channels to extracellular Zn2+ may indicate that more than one Hv1 conformational state is capable of binding Zn2+ (Qiu et al., 2016), but additional studies are needed to correlate structural changes and fluorescence signal (Randolph et al., 2016) changes with functional readouts such as GAQ (Randolph et al., 2016) and GSH (Randolph et al., 2016) gating. Collectively, data presented by Qiu et al. (2016) suggest that E1.58, H2.40, and D3.61 are necessary for Zn2+ modulation of H+ channel gating in C. intestinalis Hv1. However, putative Zn2+-coordinating oxygen atoms in a C. intestinalis Hv1 model structure are 4.2 Å (E1.58) and 3.7 Å (D3.61) away from the metal ion (Qiu et al., 2016), which is too far for these ligands to contribute to first solvation shell interactions with Zn2+ (Stote and Karplus, 1995; Dudev et al., 2003; Laitaoja et al., 2013). Additional studies are therefore needed to elucidate the atomic structure of the extracellular Zn2+ binding sites and to determine how Zn2+ occupancy is linked to changes in GAQ gating.

### Electrophysiology

FlpIn293-TRex cells (Thermo Fisher Scientific) were transfected with cDNAs encoding an N-terminal Venus-tagged human Hv1 fusion protein (Hv1) and subjected to whole-cell voltage clamp electrophysiology 24–48 h later as previously described (Ramsey et al., 2010; Villalba-Galea, 2014; Randolph et al., 2016). Briefly, pipette and bath solutions contained (in mM) ∼80 TMA-MeSO3, 100 Bis-Tris, and 8 HCl, pH 6.5, at room temperature (∼21°C), 310–320 mOsm. An acidified 1-M stock solution of ZnCl2 was serially diluted in bath recording solution to the indicated concentrations, and solution changes were accomplished by a gravity-fed superfusion of the bath recording chamber. Cells were superfused for at least 2 min between solution changes to ensure that the superfused [Zn2+] had reached steady state. Currents were measured using a model 2400 integrating patch clamp amplifier (A-M Systems) and digitized for storage and analysis using a USB-2551 A/D interface controlled by Lab View software (National Instruments) as described previously (Villalba-Galea, 2014) [REF submitted manuscript]. Signals were low-pass filtered (2 kHz) online and digitally sampled at 20 kHz. Current records were analyzed using Clampfit 9.2 (Molecular Devices), Origin 6.0 (Microcal), and custom software that is available on request (Villalba-Galea, 2014).

Conductance during the applied voltage step (GSTEP) was calculated from the current measured at the end of the voltage step (ISTEP) by GSTEP = ISTEP/V−VREV, where VREV is the measured reversal potential of the ISTEP-V relation, as previously described [REF submitted manuscript]. Data were fit to a Boltzmann function of the form

$GSTEP=(GSTEPmax)-(GSTEPmin)1+eV-V0.5dx+GSTEPmin,$
(1)

where V0.5 is the voltage at which 50% of the maximum GSTEP is reached, dx is a slope factor, and GSTEPmax and GSTEPmin represent the maximum and minimum GSTEP amplitudes, respectively. After subtraction of GSTEPmin, GSTEP data are normalized to the fitted value of GSTEPmax. V0.5 values in Table 1 represent means ± SEM from n determinations in separate cells; fitted curves and V0.5 and dx values reported in Figs. 2, 3, 4, and S1 represent Boltzmann fits to the mean data shown in plots.

Various mutations tested here have different effects on the position of the GSTEP-V relation; to measure Zn2+ effects at similar apparent POPEN for all mutants, the value of VSTEP used to activate GAQ in the presence and absence of Zn2+ was varied commensurately with the V0.5 for that mutant (Fig. S1 and Table 1). The midpoint of Zn2+ concentration responses were determined by calculating GSTEP in the absence (GSTEPControl) and presence (GSTEPZn2+) of the indicated [Zn2+] and normalizing the data to GSTEPControl. Data are plotted as the fractional decrease in GSTEPZn2+ [Zn2+ effect = 1 − (GSTEPZn2+/GSTEPControl)]. Concentration–response curves are fit to a Hill function of the form,

(2)

where IC50 is the [Zn2+] at which 50% effect is reached and nH is the Hill coefficient. IC50 and nH values in Table 1 represent means ± SEM from n determinations of Zn2+ potency measured in separate cells; fitted curves, IC50 values, and nH values reported in Figs. 2, 3, 4, and S1 represent Hill fits to the mean data shown in plots. To compare Zn2+ effects at equivalent apparent POPEN in WT and mutant Hv1 channels, we measured currents elicited by voltage steps to the voltages indicated in Table 1. pIC50 values are calculated as log (IC50); for determination of statistical significance of ΔpIC50s, individual pIC50 values determined from data in n different cells are compared using Student’s unpaired t test with significance set at P < 0.05.

### Homology modeling and simulation

The previously described Hv1 D resting-state human Hv1 VS domain model (Randolph et al., 2016) was used as the template for subsequent modeling (Modeller 9.16; Webb and Sali, 2016) to produce a starting structure (Hv1 F) for subsequent MD simulations. We generated 100 models in each of four consecutive rounds of refinement; the best model from each round was selected based on molpdf and DOPE_HR scoring functions (Modeller 9.16) and agreement to experimentally defined structural restraints (Table S1) for continued refinement in the next stage. In the first round of refinement, the x-ray crystal structure of the Hv1 coiled-coil (CC) motif (PDB accession no. 3A2A; Li et al., 2010) was used as a template for residues 226–266 with α-helical restraints applied to S4 residues 197–266. In the second round of refinement, the dopehr_loopmodel feature (Modeller 9.16) was used to remodel the backbone of residues 185–194 and 215–225, and α-helical restraints were applied to 165–184, 216–221, and 223–226. To recapitulate the electrostatic networks determined from whole-cell patch clamp experiments (Ramsey et al., 2010), distance restraints were manually applied to Cα atoms in selected residues (Table S1). The resultant model (preHv1 F) was subjected to a third round of refinement in which dopehr_loopmodel was used to remodel the backbone structure of residues 124–139 (S1–S2 loop) and 184–186 with α-helical restraints applied to residues 132–158 (S2 helix), 165–185 (S3 helix), 216–221 (S4-CC linker helix), and 223–226 (N terminus of CC motif). For the final refinement, dopehr_loopmodel was used to remodel residue ranges 124–139 and 184–186 with α-helical restraints applied to residue ranges 132–158 (S2 helix), 165–184 (S3 helix), 216–221 (S4–CC linker helix), and 223–226 (CC). Additional experimentally derived distance restraints were used to constrain the position of amino acid side chains (Table S1). To focus on the VS domain structure, the C-terminal CC motif structure (residues 227–266) was subsequently removed from the model. In all models, the imidazolium groups of His residues were defined such that Nε1 is protonated and Nδ1 is deprotonated (HSE).

The final model (Hv1 F) was embedded into a presolvated (12,057 TIP3P waters; Jorgensen et al., 1983), preequilibrated POPC lipid bilayer (163 POPC molecules) in a simulation box of dimensions 88 Å × 92 Å × 93 Å, and POPC lipids occupying overlapping space with protein (any atom–atom distance <0.8 Å or heavy atom–heavy atom distance <1.3 Å) were deleted using the psfgen plugin in VMD 1.9.2 (http://www.ks.uiuc.edu/). The AutoIonize plugin (VMD 1.9.2) was used to add either Na+ (34 Na+ cations) or Zn2+ (34 Zn2+ cations; 100 mM) and sufficient Cl to achieve charge neutrality (Hv1 F + Na+ system: 34 Na+, 35 Cl, 11,885 waters, 55,718 total atoms; Hv1 F + Zn2+ system: 34 Zn2+, 69 Cl, 11,954 waters, 55,959 total atoms). A single Zn2+ (Hv1 F + Zn2+) or Na+ (Hv1 F + Na+) ion was manually moved into proximity of Nδ1 atoms in H140/H2.40 and H193/H3.71 in the Hv1 F model structure to allow for spontaneous metal–protein interactions. Energy minimization was achieved by 5,000 steps of conjugate gradient energy minimization. After successful energy minimization, the POPC lipids were equilibrated around Hv1 F. After a final round of equilibration with harmonic restraints to equilibrate the overall system (harmonic restraints: Cα atoms, 5 kcal/mol/Å), an all-atom MD simulation was performed using the CHARMM36 force field (MacKerell et al., 1998, 2004; Klauda et al., 2010, 2012; Best et al., 2012; Huang et al., 2017; http://mackerell.umaryland.edu/charmm_ff.shtml). To determine equilibration, simulations were assessed for stabilization of RMSD, all energy terms calculated by NAMD 2.9, and periodic boundary condition cell size. Energy minimization and equilibration were simulated using a graphical processing unit (GPU) build of NAMD 2.9. MD simulations were performed using either a GPU build of NAMD 2.9 or a central processing unit build of NAMD 2.10. Tcl forces (NAMD 2.9 or 2.10) were used to keep water out of the membrane during equilibration and MD simulations. Unless otherwise indicated, energy minimization and MD simulations were conducted under standard conditions (temperature, 300°K; pressure, 1 atm; time step, 2 fs; particle mesh Ewald electrostatics grid size, 125 Å × 125 Å × 125 Å; cutoff distance, 12 Å; switching distance, 10 Å; constant temperature control, Langevin dynamics; constant pressure control, Langevin piston).

We then conducted an MD simulation in the presence of Zn2+ (Hv1 F + Zn2+), but in contrast to our expectations from experimental data (Ramsey et al., 2006), we found that Zn2+ does not stably interact with Nδ1 atoms in H140/H2.40 or H193/H3.71, resulting in a rapid increase and the distance between these residues during the first 2 ns of the MD trajectory (Fig. 5 F). This result suggests that the geometry of Zn2+ coordination by His Nδ1 and carboxylate and/or water oxygen atoms is not sufficiently optimized during model equilibration to yield a stable coordination sphere. To address this limitation, we conducted a second MD simulation (Hv1 F·Zn2+) in which harmonic constraints (initially 5 kcal/mol/Å) were applied to bonds between Zn2+ and Nδ1 atoms of H140/H2.40 and H193/H3.71 to mimic stable first solvation shell interactions between nitrogen atoms and Zn2+ (Ramsey et al., 2006). After 50 ns of Hv1 F·Zn2+ simulation time, harmonic constraints were decreased to 3 kcal/mol/Å; 20 ns later, harmonic constraints were removed (Fig. 5 G). After removal of Nδ1-Zn2+ harmonic constraints, additional MD simulation of the Hv1 F·Zn2+ system (130 ns) was conducted identically to the Hv1 F + Zn2+ systems. We also introduced H140A or H193A mutations into the Hv1 F model structure and conducted MD simulations (Hv1 F H140A·Zn2+ and Hv1 F H193A·Zn2+) in which harmonic constraints (5 kcal/mol/Å) were applied to bonds between Zn2+ and the Nδ1 atom of H140/H2.40 (in H193A) or between Zn2+ and the Nδ1 atom of H193/H3.71 (in H140A), respectively (Fig. S4). A third MD simulation (Hv1 F + Na+) was conducted in the presence of Na+ instead of Zn2+ without harmonic constraint between Nδ1 atoms and Na+ (Fig. S5 and Tables S2, S3, S4, and S5).

Images of atomic structures were rendered using Tachyon ray tracing in VMD 1.9.2 or 1.9.3 (Humphrey et al., 1996; Stone, 1998). Salt bridges (O–N distance cutoff ≤3.2 Å) and H-bonds (donor–acceptor distance ≤3.0 Å; angle cutoff, 20°) in Hv1 F and Hv1 F·Zn2+ MD simulation trajectories were identified using plugins in VMD 1.9.2 or 1.9.3 (Humphrey et al., 1996). Model structures were compared after structural alignment using the MultiSeq STAMP structural alignment plugin in VMD 1.9.2 or 1.9.3 (Humphrey et al., 1996; Eargle et al., 2006; Roberts et al., 2006). Atomic distance changes during MD simulation trajectories were measured using VMD 1.9.2 or 1.9.3 (Humphrey et al., 1996) by manually picking selected atoms. To ensure proper treatment of Zn2+ by the CHARMM36 force field (MacKerell et al., 1998, 2004; Klauda et al., 2010, 2012; Best et al., 2012), the radial distribution function (RDF) plugin in VMD 1.9.2 (Humphrey et al., 1996) was used to calculate the Zn2+ RDF for every frame for the Hv1 F·Zn2+ and Hv1 F + Zn2+ systems, as illustrated for selected Zn2+–ligand interactions during MD simulations (Figs. 6 and S7) that reproduce experimentally observed Zn2+ RDFs (Hitoshi et al., 1976; Johansson, 1992; Stote and Karplus, 1995). VMD 1.9.2 (Humphrey et al., 1996) was used to calculate the distance of interactions between Zn2+ and H2.40-Nδ1 or H3.71-Nδ1 atoms for every frame in Hv1 F·Zn2+ and Hv1 F + Zn2+ systems. Distance histograms (bin width = 0.25 Å) were constructed using OriginPro 8.1. The Salt Bridges plugin of VMD 1.9.2 (Humphrey et al., 1996) was used to determine the salt bridges for every frame in the Hv1 F·Zn2+, Hv1 F + Zn2+, and Hv1 F + Na+ systems. To determine if changes in salt bridges observed during MD simulations were correlated to distances between Zn2+ and H2.40-Nδ1 or H3.71-Nδ1, the Pearson’s correlation coefficient (PCC) was calculated using custom Python 2.7 scripts (http://www.python.org). H-bond occupancy was calculated using the VMD 1.9.2 (Humphrey et al., 1996) Hydrogen Bonds plugin over every frame in the Hv1 F·Zn2+, Hv1 F + Zn2+, and Hv1 F + Na+ systems using default settings described previously.

### Online supplemental material

Fig. S1 shows Zn2+ concentration–response curves for E119H at various voltages, scatter plots of pIC50 or τACT values in function of V0.5, and Zn2+ concentration–response curves D112N, D112H, and E153A. Fig. S2 shows a comparison of C. intestinalis Hv1, Hv1 FL, mHv1cc (PDB accession no. 3WKV), and Hv1 D model structures. Fig. S3 shows Zn2+ coordination geometries in Hv1 F·Zn2+ and Hv1 F + Zn2+; water occupancy in Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c0, and Hv1 F + Zn2+; and ZN13 to D1.51-Cδ1 and D1.51-Cδ1 distance histograms for Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c0, and Hv1 F + Zn2+ MD trajectories. Fig. S4 shows representations and distance plots for Hv1 F·Zn2+ H140A and H1933 model MD trajectories. Fig. S5 shows comparisons of Zn2+ and Na+ positions in snapshots taken from Hv1 F + Zn2+ and Hv1 F + Na+ MD trajectories and distance plots for selected atoms in the Hv1 F + Na+ MD trajectory. Fig. S6 shows plots of distances between selected terminal side-chain atoms in the Hv1 F·Zn2+ c0 MD trajectory, snapshots showing differences in Zn2+ position before (Hv1 F·Zn2+ c5) and after (Hv1 F·Zn2+ c0) ZN13 dissociation from H2.40-Nδ1 and H3.71-Nδ1, and PCC for distance between terminal atoms of selected salt bridge pairs and distance between ZN13 and H2.40-Nδ1 or H3.71-Nδ1 during Hv1 F·Zn2+ c0 and Hv1 F·Zn2+ c5 MD trajectories. Fig. S7 shows salt bridge distance plots and amplitude histograms for selected residues in the Hv1 F·Zn2+ c0 MD trajectory. Fig. S8 shows distance plots for selected atom pairs during Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c3, and Hv1 F·Zn2+ c0 MD trajectories.

Table S1 shows distance restraints for the last three rounds of Hv1 F model refinement. Table S2 shows mean Cα – Cα distances in Hv1 F + Zn2+, Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c0, and Hv1 F + Na+ MD trajectories. Table S3 shows RMSD calculations for Hv1 F + Zn2+, Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c0, and Hv1 F + Na+ MD trajectories. Table S4 shows mean salt bridge distances in Hv1 F + Zn2+, Hv1 F·Zn2+ c5, and Hv1 F·Zn2+ c0 MD trajectories. Table S5 shows occupancy of selected H-bonds in Hv1 F + Zn2+, Hv1 F·Zn2+ c5, Hv1 F·Zn2+ c0, and Hv1 F + Na+ MD trajectories.

The locations of candidate Zn2+-coordinating side chains in a previous resting-state Hv1 model structure (Hv1 D; Randolph et al., 2016) and the mHv1cc crystal structure (Takeshita et al., 2014) are shown in Fig. 1. In Hv1 D, Cγ atoms in H2.40 and H3.71 are separated by ∼20 Å (Fig. 1 B), and His nitrogen atoms are therefore too distant to simultaneously substitute for water oxygen atoms in the first solvation shell (distance ≤2.5 Å) of a liganded Zn2+ ion (Stote and Karplus, 1995; Dudev et al., 2003; Laitaoja et al., 2013). We remodeled extracellular loops in Hv1 D to generate a new Hv1 resting-state model (Hv1 F) in which the Cγ atoms in H2.40 and H3.71 are separated by only 7.3 Å and imidazolium Nδ1 atoms are close enough to potentially participate in first-shell interactions with Zn2+ (Fig. 1 C). Except for an added S4–S5 linker helix (5L) in Hv1 F, the positions of backbone atoms in S1–S4 helical segments in Hv1 D and Hv1 F are highly similar, indicating that the remodeling procedure did not substantively alter the overall structure (Fig. 1 E). Over the final 20 ns of Hv1 D and initial 20 ns of Hv1 F + Zn2+ MD simulations, the overall backbone RMSD is 4.2 Å, and RMSDs in individual helical segments are 2.5 Å (S1), 3.0 Å (S2), 2.7 Å (S3), and 1.7 Å (S4). Hv1 F thus represents a new candidate Zn2+-liganded conformation of the Hv1 VS domain. To identify side chains that are close to the coordinated Zn2+ ion, we mutated candidate extracellular residues to Ala or His and measured their effects on the potency of Zn2+ to modulate GAQ gating. If the residue contributes to stabilizing a divalent metal-liganded conformation, Ala mutations should decrease, and His substitutions increase, Zn2+ potency. Because mutations may also alter voltage-dependent gating (Ramsey et al., 2010; Randolph et al., 2016), we first determine the midpoint (V0.5) for voltage-dependent activation of GAQ from GSTEP-V relations and adjust the voltage used to measure Zn2+ sensitivity such that each mutant is tested at similar apparent POPEN (Figs. 2 B, 3 B, and S1, and Table 1).

Consistent with a previous study (Ramsey et al., 2006), Zn2+ potently modulates GAQ gating in WT Hv1 (Fig. 2, A and C). H140A and H193A mutations, alone or together, cause only small shifts in V0.5 (Fig. 2 B), and GAQ activation kinetics are also similar in WT, H140A, H193A, and H140A-H193A (Table 1), indicating that H2.40 and H3.71 do not directly control VS activation gating. As expected (Ramsey et al., 2006), Zn2+ potency is markedly decreased in H140A and H193A (H140A: IC50 = 38.9 ± 9.4 µM; H193A: IC50 = 87.8 ± 23.8 µM) compared with WT Hv1 (WT: IC50 = 3.8 ± 1.7 µM), and Zn2+ sensitivity in the H140A-H193A double mutant is shifted into the millimolar range (H140A-H193A: IC50 = 1.7 ± 0.3 mM; Fig. 2 C and Table 1). Relative changes in Zn2+ potency between selected mutants are reported by ΔpIC50 values: mutations that decrease Zn2+ potency yield positive ΔpIC50s (Table 2). For example, H140A-H193A causes a 2.7 log unit decrease in Zn2+ potency compared with WT (ΔpIC50 = −2.7; Table 2). E119A also decreases Zn2+ potency when engineered into the background of H193A (ΔpIC50 = −0.5; Fig. 4 A and Table 2), consistent with the hypothesis that H2.40 and E1.58 might interact to coordinate Zn2+ (Takeshita et al., 2014; Qiu et al., 2016). Zn2+ potency is also decreased in D112H and E153A, but E119A, D123H, D130H, and E196H do not produce statistically significant changes (Fig. 3 C; Fig. 4, A–D; Fig. S1 E; and Table 2). The inability of His substitutions at D1.51, D1.62, D1.69, and E4.38 to increase Zn2+ sensitivity suggests that the side chains of these residues are likely to be too far away from Zn2+ to directly contribute to its coordination sphere. The effect of E153A to reduce Zn2+ potency is surprising because this residue is part of the intracellular Coulombic network (ICN) and thus distant the from the extracellular Zn2+ coordination site in Hv1 F (E2.53-Cα and ZN13 are separated by 23 Å in Hv1 F; Figs. 1 C and 5 A).

Mutations in Hv1 can independently affect Zn2+ potency and GAQ gating. For example, V0.5 differs by only 11.2 mV in H140A-H193A and E119H-H140A-H193A (Table 1), but Zn2+ IC50 increases 2.4-fold (from 1.7 mM in H140A-H193A to 7.3 µM in E119H-H140A-H193A; Tables 1 and 2). The converse is observed in D185H, where V0.5 shifts positively (+67.6 mV relative to WT Hv1), but Zn2+ potency (IC50 = 4.2 µM) is not different from WT (Fig. 4 D and Table 1). In E153A, V0.5 is shifted −57 mV (Table 1) and Zn2+ potency decreases substantially (ΔpIC50 = −1.29; Fig. S1 F and Table 2) compared with WT Hv1 (Tables 1 and 2). A plot of pIC50 versus V0.5 values for all mutants tested here reveals that the effects of mutations on V0.5 and pIC50 are poorly correlated (R2 = 0.08; Fig. S1 C). The time constant for GAQ activation is also poorly correlated with V0.5ACT; R2 < 0.001; Fig. S1 D). D112H causes V0.5 to shift +55 mV, consistent with our previous work showing that D1.51 and D3.61 likely interact with R1 (R4.47) to stabilize an activated-state conformation of the Hv1 VS domain (Randolph et al., 2016). However, His substitutions of other extracellular acidic side chains (E1.58, D1.62, and D1.69, but not E4.38) also shift the GAQ-V relation positively. The effect of D185H on V0.5 is also strongly dependent on the mutant background (Table 1). Together, the data suggest that His mutations of extracellular acidic residues may perturb Coulombic interactions that serve to stabilize an activated-state VS conformation. In summary, D1.51, E1.58, D1.62, D1.69, and D3.61 mutations produce the same effect on GAQ gating (positive shift) as extracellular Zn2+ (Cherny and DeCoursey, 1999), whereas E2.53 mutations shift GAQ gating negatively (Ramsey et al., 2010).

In contrast to E119H, D123H and D130H single mutants do not significantly alter Zn2+ potency (Fig. 3 C; Fig. 4, A–C; and Table 2). The lack of effect in D123H is somewhat surprising given that in the mHv1cc x-ray structure, the D1.62 side chain is close to the putative Zn2+ binding site (Fig. 1 D; Takeshita et al., 2014). However, when H2.40 is neutralized, His substitution at D1.62 does appear to partially restore Zn2+ potency (D123H-H140A vs. H140A: ΔpIC50 = 0.47; Fig. 4 B and Table 2). Consistent with a previous study showing that suggesting that D1.62 neutralization does not alter Zn2+ sensitivity (Takeshita et al., 2014), our data suggest that although this side chain probably does not directly contribute to Zn2+ coordination in WT Hv1, it could participate in the coordination sphere in H140A and H193A mutant channels (Fig. S4, B and C). In contrast to D1.62, introducing His at D1.69 dramatically increases Zn2+ sensitivity the background of H193A (D130H-H193A vs. H193A: ΔpIC50 = 1.44; Fig. 4 C and Table 2). Zn2+ is nearly equipotent in D130H-H193A (IC50 = 3.7 ± 2.7 µM; Table 1) and WT Hv1 (IC50 = 3.8 ± 1.7 µM; Table 1). His substitution of D1.69 suffices for Zn2+ coordination in the absence of H3.71 but not in the absence of H2.40 (Table 2). We find that several other His mutant combinations do not significantly alter Zn2+ potency (Tables 1 and 2), suggesting that the gain-of-function effect depends on local structural features that differ in the various mutant proteins. The observation that WT-like Zn2+ potency is measured in E119H-H140A, E119H-H193A, E119H-H140A-H193A, and D130H-H193A suggests that Zn2+ coordination spheres are likely to be structurally distinct in various mutant proteins, highlighting a potential limitation of mutagenesis strategies for understanding the mechanism of Zn2+ effects in WT Hv1.

Among the mutants tested here, E119H appears to be special in producing a gain-of-function phenotype with respect to Zn2+ potency. E119H decreases Zn2+ IC50 from 3.8 µM (WT) to 0.6 µM (ΔpIC50 = 0.8; Fig. 3 C and Tables 1 and 2), consistent with studies indicating that this side chain is close to the liganded Zn2+ ion (Takeshita et al., 2014; Qiu et al., 2016). The effect of E119H is magnified in the background of H140A and/or H193A mutants; in E119H-H140A-H193A, Zn2+ potency increases 2.4 log units (Fig. 4 A and Table 2). Zn2+ modulates Hv1 gating more potently in E119H and E119H-H140A mutants than in WT, and five of the six largest positive ΔpIC50s are measured in mutants containing E119H (Fig. 4 A and Table 2). In summary, we find that mutations of both extracellular (E1.58) and intracellular (E2.53) acidic residues affect Zn2+ potency. One possibility, which we explore later (Fig. 6), is that extracellular Zn2+-coordinating and ICN structures are functionally coupled in Hv1, similar to allosteric modulation in G-protein–coupled receptors (GPCRs; Wacker et al., 2017).

To gain insight into the structural basis for Zn2+ coordination, we conducted an MD simulation of the Hv1 F model structure. Prior to MD, a Zn2+ (ZN13) was manually positioned between E1.58, H2.40, and H3.71 to encourage spontaneous coordination by His-Nδ1 and Glu carboxylate oxygen atoms (Figs. 1 C and 5 C). The Hv1 F + Zn2+ model is stable during MD simulation (Fig. 5 E), but we do not observe stable first solvation shell coordination of ZN13 by E1.58, H2.40, or H3.71 side chains (Fig. 5, D and F). Instead, H2.40 and H3.71 migrate to positions that are separated by >15 Å, similar to the Hv1 D template (Figs. 1 B and 5 F). To optimize the geometry of the Zn2+ coordination sphere in Hv1 F, we placed harmonic constraints (5 kcal/mol/Å) on ZN13 to H140-Nδ1 and ZN13 to H193-Nδ1 bonds and ran a separate “constrained” MD simulation (Hv1F·Zn2+ c5). As expected, the distance (∼2.1 Å) between ZN13 and Nδ1 atoms in H2.40 and H3.71 remains stable in Hv1F·Zn2+ c5 (Fig. 5, D, G, and H). The Hv1F·Zn2+ c5 system also exhibits less backbone dynamic fluctuation than Hv1 F + Zn2+, as evidenced by differences in Cα atom RMSD values (Fig. 5 G and Table S3).

As expected for Zn2+ in solution (Stote and Karplus, 1995), we observe that ions are persistently coordinated in octahedral geometry during MD simulations. In Hv1F·Zn2+ c5, H2.40-Nδ1, H3.71-Nδ1, E1.58-Oε1, E1.58-Oε2, and two water (W213 and W1080) oxygens form a stable His2-Glu Zn2+ coordination sphere, and Zn2+ coordination is maintained when His Nδ1-Zn2+ harmonic constraints are decreased (3 kcal/mol/Å in Hv1F·Zn2+ c3; Fig. 5, D and H; and Fig. S3 A). In the absence of external harmonic constraints (Hv1 F + Zn2+), Zn2+ ions are octahedrally coordinated by six water oxygens (Figs. 5 B and S3 B). Histograms of Zn2+-ligand distances in both Hv1 F + Zn2+ and Hv1F·Zn2+ simulations show a pattern of peaks that is consistent with the expected RDF for Zn2+ interactions (Figs. 5 D and S7 B; Stote and Karplus, 1995). In a control simulation conducted with Na+ instead of Zn2+ (Hv1 F + Na+), we also observe constitutive octahedral coordination of Na+ by six waters (Fig. S5). Comparing Hv1 F + Zn2+ to Hv1 F + Na+ shows that Zn2+ and Na+ ion stably occupy similar positions in the hydrated extracellular vestibule, but in neither of these systems do we observe proteins atoms engaging in first solvation shell interactions with the metal ion (Fig. S5). Although we observe longer-distance interactions between Zn2+ and potential ligand atoms in MD simulations (Figs. 5 D and S7), we do not analyze such interactions further here because the experimental data (Ramsey et al., 2006; Musset et al., 2010; Takeshita et al., 2014; Qiu et al., 2016) indicate that first solvation shell interactions with H2.40, H3.71, and E1.58 are likely to be the main determinants of potent modulation of GAQ gating by Zn2+, potentially explaining why Na+ does not substitute for divalent metal cations (Cherny and DeCoursey, 1999; Ramsey et al., 2006; Musset et al., 2010; Takeshita et al., 2014; Qiu et al., 2016). Our results are consistent with MD simulations and x-ray structures of class A GPCRs that are modulated by Na+ via interactions with a highly conserved acidic side chain (Katritch et al., 2014).

Importantly, His2-Glu coordination of ZN13 persists for >40 ns after harmonic constraints are removed (Hv1F·Zn2+ c0; Fig. 5 H). The lack of spontaneous first solvation shell interactions between Zn2+ and H140-Nδ1 or H193-Nδ1 in Hv1 F + Zn2+ (Fig. 5 F) is therefore explained by insufficient optimization of ZN13 coordination geometry before MD simulation. In contrast to Hv1F·Zn2+, the first solvation shells of Zn2+ ions are constitutively occupied by six waters in Hv1 F + Zn2+ (Fig. S3, A and B). The stability of the His2-Glu coordination in the absence of externally applied forces directly demonstrates that an empirical force field can be used to simulate Zn2+-protein interactions in silico. Furthermore, Zn2+ ions consistently interact with protein atoms at discrete distances, consistent with the expectations for the Zn2+ RDF (Figs. 5 D, 6 C, and S7; see Materials and methods; Hitoshi et al., 1976; Johansson, 1992; Stote and Karplus, 1995). Although H2.40-Nδ1 and H3.71-Nδ1 leave the first solvation shell of ZN13 after harmonic constraints are removed, the coordinating nitrogen atoms remain in close proximity (Fig. 5 H) even though local backbone protein structure is thermodynamically stable before and after the ZN13 dissociation, as indicated by Cα atom RMSDs at H140 (0.6 Å) and H193 (0.9 Å) calculated for snapshots at t = 70 ns and t = 130 ns of the Hv1F·Zn2+ c0 MD trajectory (Fig. 5 H).

Zn2+ dissociation from H2.40 and H3.71 Nδ1 atoms during Hv1F·Zn2+ c0 is correlated with the reorganization of H-bonds and Coulombic interactions (i.e., salt bridges), and changes in protein backbone structure result. Although ZN13 remains octahedrally coordinated throughout the Hv1F·Zn2+ trajectory, the identities of coordinating ligands change (Fig. 5, G and H). After ∼40 ns, first H3.71-Nδ1 and later H2.40-Nδ1 are seen to dissociate from the first solvation shell of ZN13; exchange of first-shell waters (W213 and W1080) occurs within 10 ns of the change in H3.71-Nδ1 coordination (Fig. 5 H). At the end of the Hv1F·Zn2+ c0 simulation, ZN13 is octahedrally coordinated by E1.58-Oε1, E1.58-Oε2, and four different water oxygens. ZN13 also migrates deeper into the extracellular vestibule, toward D112/D1.51 and F150/F2.50 (Fig. S6). The inward migration of Zn2+ observed during Hv1F·Zn2+ c0 after ZN13 dissociates from H2.40-Nδ1 and H3.71-Nδ1 is reminiscent of that in a previous study (Qiu et al., 2016). However, D1.51-Oδ1 does not participate in first-shell interactions with Zn2+ in any of our MD simulations. In Hv1 F + Zn2+, two hydrated Zn2+ ions occupy the extracellular vestibule, close to the location of ZN13 in the Hv1F·Zn2+ c5 and Hv1F·Zn2+ c0 systems (Fig. 5, B and C). In summary, we observe that Zn2+ is persistently coordinated in an octahedral geometry by waters and side-chain atoms of E1.58, H2.40, and H3.71.

Comparing structural differences between Hv1F·Zn2+ c5, Hv1F·Zn2+ c0, and Hv1 F + Zn2+ systems provides insight into the possible mechanism by which Zn2+ exerts an effect on GAQ gating. In Hv1F·Zn2+ c5, we observe stable Coulombic interactions between terminal atoms of residues that participate in salt bridges, but marked changes in the distances between salt-bridged side chains occur during the Hv1F·Zn2+ c3 and Hv1F·Zn2+ c0 simulations (Fig. 5 G). For example, the D3.50-R3 distance decreases by ∼5 Å between Hv1F·Zn2+ c0 and c5 (Fig. 5 G), indicating that the salt bridge is substantially strengthened. In contrast, the D1.51-R2 salt bridge begins to weaken even before ZN13 dissociates from H2.40-Nδ1 and H3.71-Nδ1 (Fig. 5 G). The D1.51-R1 interaction remains strong throughout the MD trajectory, resulting in a net movement of the D1.51-R1 pair away from other residues in the ICN (Fig. 5 G). Together, side chain and backbone movements observed in constrained and unconstrained Hv1F MD simulations indicate that reorganization of Coulombic interactions and Zn2+ occupancy of the His2-Glu coordination sphere exhibit correlated movements that are suggestive of allosteric coupling.

Snapshots taken from Hv1F·Zn2+ and Hv1 F + Zn2+ also show that although the positions of D1.51, F2.50, and other nearby residues are similar, the organization of salt bridges between conserved S4 Arg “gating charge” side chains and acidic “counter charges” in S2 and S3 is characteristically different (Fig. 6, A and B). During the full Hv1F·Zn2+ MD simulation, both short- (<4 Å) and longer- (>4 Å) range Coulombic interactions are measured (Fig. 6, C–E; and Fig. S7). For example, stable short-distance salt bridges are observed between E2.53-R2 and D3.50-R3 pairs in Hv1F·Zn2+ c0 (Fig. 6, A, D, and E; and Fig. S7). However, in Hv1 F + Zn2+, E2.53 and D3.50 salt bridge partners are swapped, and E2.53-R3 and D3.50-R2 now form stable pairs (Fig. 6, B, D, and E). A stable D1.51-R1 salt bridge is observed in both Hv1F·Zn2+ and Hv1 F + Zn2+ simulations (Fig. 6, A and B), consistent with experimental data showing that D112 is not necessary for Zn2+ interactions (Tables 1 and 2 and Fig. S1 E). Occupancy of H-bonds between donor and acceptor atoms in ICN side chains shows a similar pattern to salt bridges: E2.53-R2 and D3.50-R3 H-bonds are common in Hv1F·Zn2+ c0 but absent from Hv1 F + Zn2+ (Fig. 6 F). In contrast, D3.50-R2 H-bonds exhibit high occupancy in Hv1 F + Zn2+, but are rare in Hv1F·Zn2+ c0 (Fig. 6 F). Overall, it appears that the number and strength of interactions between S4 Arg and S2/S3 acidic groups are increased when Zn2+ interactions are constrained.

To test the hypothesis that dissociation of ZN13 from its coordination site is correlated with changes in ICN salt bridge distances, we measured atomic distances in pairs of interacting groups in the presence (Hv1F·Zn2+ c5) and absence (Hv1F·Zn2+ c0) of ZN13-Nδ1 harmonic constraints. We measured distances between terminal atoms in identified salt-bridged residue pairs and distances between terminal nitrogen and oxygen atoms in each frame of the Hv1F·Zn2+ c5 and Hv1F·Zn2+ c0 MD trajectories and analyzed covariance and correlations (see Materials and methods). Correlated changes in salt bridge distance and ZN13 to H2.40-Nδ1 or ZN13 to H3.71-Nδ1 distance are reported by PCC for each pair (Table 3). Several interaction pairs exhibit positive or negative PCC values in Hv1F·Zn2+ c0 that are not observed in the Hv1F·Zn2+ c5 MD simulation (Fig. 6 G, Fig. S6 G, and Table 3). For example, the increase in D1.51-R2 distance is positively correlated (PCC > 0.5) with an increase in both the ZN13 to H2.40-Nδ1 and ZN13 to H3.71-Nδ1 distances (Fig. 6 G and Table 3), consistent with our observation that the D1.51-R1 pair moves away from the ICN after ZN13 dissociation from H2.40-Nδ1 and H3.71-Nδ1 (Fig. 5 G). Also, consistent with previously described changes in ICN architecture (Fig. 6, A–E), we observe a strong negative correlation between the E2.53-R2 and D3.50-R3 salt bridge distance and ZN13 to H2.40-Nδ1 and ZN13 to H3.71-Nδ1 distances (Fig. 6 G). In summary, changes in both salt bridge and H-bond networks in the ICN are strongly correlated with changes in Zn2+ occupancy of the primary His2-Glu coordination sphere.

Zn2+-dependent changes in Coulombic and H-bond network structures are also reflected by dynamic reorganization of the protein backbone. Changes in the distance between Cα atoms of key side chains are correlated with the relief of ZN13-Nδ1 harmonic constraints or dissociation of Zn2+ from its coordination site in Hv1F·Zn2+ c0 (Fig. S8). For example, D1.51-Cα moves ∼4 Å closer to D3.61-Cα and ∼4 Å away from D3.50 and R2 Cα atoms (Fig. S8 A). D1.51-Cα also gets ∼6 Å closer to ZN13 in Hv1F·Zn2+ c0, partly because D1.51 moves extracellularly (Fig. S8 G) and partly because ZN13 moves intracellularly (Fig. S6). Reorganization of the Zn2+-liganded E1.58 side chain also results in an increase in the R1-Cα to E1.58-Cα distance (Fig. S8 E), whereas the D1.51-Cα to E1.58-Cα distance remains almost unchanged in the unconstrained simulations (Fig. S8 A). R1-Cα and R2-Cα atoms move closer to F2.50-Cα (Fig. S8, E and F), indicating that S4 undergoes a small (2.0–2.5 Å) outward translation during the Hv1F·Zn2+ c0 simulation. In contrast to Hv1F·Zn2+ c0, Cα atom distances are relatively constant in the Hv1 F + Zn2+ (Fig. S8 H) and Hv1 F + Na+ control simulations (Fig. S5 G). Changes in backbone structure demonstrate that alterations in the patterns of Coulombic interactions and H-bonds do not merely result from dynamic side-chain switching in structurally similar conformations. Overall, our data demonstrate that side-chain interactions among highly conserved ICN residues are reorganized in a Zn2+-dependent fashion. Extracellular Zn2+ coordination and voltage-dependent gating therefore appear to be allosterically coupled by long-range conformational rearrangements within the Hv1 VS domain.

Here we show that His mutations of extracellular residues in Hv1 are sufficient to confer a gain of function by increasing Zn2+ potency effects. Among the mutations tested here, E119H has the most remarkable effect. Substituting His at position E1.58 lowers the Zn2+ IC50 by 0.8 log units, and WT-like Zn2+ potency is maintained even in the absence of H2.40 and H3.71 (Tables 1 and 2), indicating that E119H fully substitutes for the loss of these His residues. A straightforward interpretation of the experimental data is that an imidazole nitrogen atom in the introduced His side chain of E119H contributes to the Zn2+ coordination sphere via a first solvation shell interaction. This interpretation is reinforced by data showing that E119H-H140A-H193A responds to Zn2+ equipotently with WT Hv1. The ability of H140A, H193A, and E119H-H140A-H193A mutants (and nonmammalian Hv1 orthologues containing H3.71 substitutions) argues that the effect of Zn2+ to modulate GAQ gating is unlikely to require a H2.40- or H3.71-dependent “metal ion lock” that shifts voltage dependence merely by restricting the mobility of S4.

The relatively potent effects of Zn2+ in D123H, D130H, and D185H suggest that additional extracellular acidic side chains can function as Zn2+ ligands in Hv1, as previously suggested (Takeshita et al., 2014; Qiu et al., 2016; Fig. 1 and Tables 1 and 2). However, single His mutants at D1.62, D1.69, and D3.61 are insufficient to lower Zn2+ IC50 like E119H, and the effects of His substitutions at these positions are strongly sensitive to the mutant background, suggesting that conformational heterogeneity of the Zn2+ coordination site structure could account for moderately potent (i.e., IC50 = 10–100 µM) responses in mutant and nonmammalian Hv1 channels (Sasaki et al., 2006; Qiu et al., 2016). The apparent ability of Hv1 to coordinate Zn2+ in more than one conformation suggests that unidentified side chains or conformational changes may therefore be necessary for Zn2+ coordination in mutant and nonmammalian Hv1 channels, and caution is therefore warranted in the interpretation of neutralizing mutagenesis studies designed to identify metal ion binding sites (Takeshita et al., 2014; Qiu et al., 2016). In contrast to D1.62, D1.69, and D3.61, His substitutions at D1.51 and E4.38 decrease Zn2+ potency, suggesting that these residues are too distant from the primary coordination site to participate in first solvation shell interactions with the divalent metal cation.

Our experimental data are consistent with previous studies in that they support the hypothesis that H2.40 and H3.71 form a stable Zn2+ coordination sphere together with E1.58 (Ramsey et al., 2006; Takeshita et al., 2014; Qiu et al., 2016) and provide valuable spatial constraints for generating a refined structural model of the human Hv1 resting-state VS domain in a Zn2+-liganded conformation (Hv1 F). Hv1 F is the only reported Hv1 VS domain structure in which oxygen and nitrogen atoms in E1.58, H2.40, and H3.71 side chains may simultaneously serve as Zn2+ ligands (Figs. 1 C, 5 A, S2, and S3 A; Chamberlin et al., 2014; Takeshita et al., 2014; Li et al., 2015; Qiu et al., 2016; Randolph et al., 2016). Optimal Zn2+ coordination geometry in Hv1 F is observed after MD simulation in the presence of harmonic constraints between ZN13 and H2.40-Nδ1 and H3.71-Nδ1 (Fig. 5 G). In the absence of ZN13-Nδ1 harmonic constraints, the Hv1 F model structure readily adopts a conformation that is similar to the previously reported (Hv1 D; Randolph et al., 2016) template during MD simulation (Figs. 1 B and 5 F), indicating that Hv1 D represents a reasonable model of the Zn2+-free resting-state Hv1 VS domain conformation. The position of R1 is not markedly different in Hv1 D and Hv1 F (Fig. 1 E), demonstrating that the remodeled Hv1 F structure remains compatible with the structural constraints imposed by the resting-state GSH in R1H mutant Hv1 channels (Randolph et al., 2016).

Additional MD simulations conducted after ZN13-Nδ1 harmonic constraints were relieved (Hv1F·Zn2+ c0) reveal the structure of a Zn2+ coordination site in the Hv1 VS domain in atomic detail. A His2-Glu coordination sphere is observed during constrained MD simulations (Hv1F·Zn2+ c5 and Hv1F·Zn2+ c3), and this atomic architecture persists for >40 ns in the absence of externally applied forces (Fig. 5 G). During the Hv1F·Zn2+ c0 MD trajectory (t ≈ 110 ns), Zn2+ coordination by H2.40-Nδ1 and H3.71-Nδ1 is replaced by water oxygen atoms, which is unexpected if the atomic architecture of the first Zn2+ solvation shell represents a thermodynamically stable conformation. We speculate that the empirical force field used here may underestimate the strength of electronic interactions between Zn2+ and imidazole nitrogen atoms in H140 and H193, allowing water oxygens to substitute for Nδ1 atoms in the first solvation shell during Hv1F·Zn2+ c0 and preventing spontaneous coordination during Hv1F + Zn2+ MD simulations. Thus, although the backbone structure of the protein near the Zn2+ binding site is thermodynamically stable, the atomic architecture of the first Zn2+ solvation shell is evidently not, and His2-Glu coordination of Zn2+ does not persist indefinitely under our MD simulation conditions.

The aforementioned result highlights an inherent limitation in the use of an empirical force fields to simulate electronic interactions between Zn2+ and His imidazole nitrogen atoms (Stote and Karplus, 1995), and preclude accurate estimation of Zn2+-dependent free energy changes or apparent Zn2+ affinity from the MD data. On the other hand, Zn2+-ligand distance measurements (Figs. 5 D and S7) show that the CHARMM36 force field used here is appropriately parameterized to yield good estimates of coordination spheres involving protein atoms and biologically import transition metal ions, including Zn2+, and in this respect our results compare well to results obtained from semiempirical quantum mechanical methods (Yu et al., 2018). As expected for Zn2+ in aqueous solution (Burgess, 1978; Stote and Karplus, 1995), we observe that the first Zn2+ solvation shell contains six ligand atoms in an octahedral coordination geometry almost exclusively over >400 ns of total MD simulation time. Consistent with the expected Zn2+ RDF (Hitoshi et al., 1976; Stote and Karplus, 1995), we also observe second- and third-shell interactions at the expected distances from Zn2+ ions (Figs. 5 D and S7 B). In contrast to a previous study in which Zn2+ position was constrained using an umbrella sampling protocol (Qiu et al., 2016), we observe that Zn2+ interacts exclusively with protein side-chain and water oxygen atoms, suggesting that applying constraints to specific atoms based on experimental data can yield a more refined picture of the metal ion coordination geometry.

The in silico approach used here offers new insights into the structural determinants of Zn2+ coordination in Hv1. Although D1.51 (D160) mutation in C. intestinalis Hv1 alters voltage-dependent fluorescence changes (Qiu et al., 2016), we show that D112N does not decrease and D112H does not increase Zn2+ potency (Tables 1 and 2 and Fig. S1 E), and we do not observe first solvation shell interactions between Zn2+ and the D1.51 side chain in any of our MD simulations (Figs. S3 F and S6 F). Hv1 F is distinct from a C. intestinalis Hv1 resting-state model in which Zn2+ participates in first solvation shell interactions with a D233/D3.61 carboxylate oxygen (d = 2.7 Å), the hydroxyl oxygen of S229/S3.57 (d = 2.1 Å), a backbone oxygen (d = 1.9 Å), and two or three water oxygens at the putative “site 2” (Qiu et al., 2016). Although our experimental and computational results do not support the hypothesis that Zn2+ interactions with or near D1.51 are required for the observed effects on GAQ gating, we do observe that Zn2+ migrates inward toward the previously defined “site 2” after Zn2+-Nδ1 harmonic constraints are relieved (Fig. S6, C–E). In our study, Zn2+ (or Na+) ions that occupy the vicinity of “site 2” (Qiu et al., 2016) do not make first solvation shell interactions with the H2.40 or H3.71 side chains that are necessary for micromolar Zn2+ sensitivity (Fig. 2; Fig. 3; Fig. S5, A–F; Fig. S6, E and F; and Tables 1 and 2; Ramsey et al., 2006). Additional studies will be needed to determine whether Zn2+ affects GAQ gating by occupying a “deeper” site, closer to D1.51, in mutants such as E119H-H140A-H193A.

We hypothesize that Hv1F·Zn2+ represents a biologically relevant Zn2+-liganded Hv1 resting-state structure, and that structural insights revealed in our MD simulations may be generally useful for understanding Zn2+ coordination at solvent-accessible sites in metalloproteins. In Hv1F·Zn2+ MD simulations, we observe His2-Glu-(H20)2 (bidentate Glu Oε1 and Oε2 in Zn2+ first solvation shell) or His2-Glu-(H20)3 (monodentate Glu Oε1 or Oε2 interaction with Zn2+) coordination spheres that are distinct from the tetrahedral geometry commonly observed for Zn2+ in high resolution in NMR and x-ray structures (Dudev et al., 2003; Laitaoja et al., 2013). Whereas Zn2+ binding sites in solution NMR and x-ray structures appear to be relatively inaccessible to bulk solvent (Dudev et al., 2003; Laitaoja et al., 2013), Zn2+ and its ligands are located in the highly hydrated VS central crevice in our MD simulations (Fig. S3, C–E). The rapid reversibility and pHO dependence of Zn2+ effects on GAQ gating strongly argue that the binding site is readily solvent accessible (Cherny and DeCoursey, 1999; Ramsey et al., 2006; Musset et al., 2010; Takeshita et al., 2014; Qiu et al., 2016). Consistent with the solvent accessibility of the coordinated ZN13, we directly observe exchange of Zn2+ ligands (H2.40-Nδ1, H3.71-Nδ1, and W213 and W1080 oxygens) for other water oxygens in the first Zn2+ solvation shell in the Hv1F·Zn2+ c0 MD trajectory after relieving ZN13-Nδ1 harmonic constraints (Fig. 5 H; and Fig. S3, A and B). Together, our experimental and computational data strongly argue that Zn2+ interacts with Hv1 residues in a well-hydrated environment, indicating that high apparent affinity and solvent accessibility are not mutually exclusive.

Rather than acting as a cofactor for enzymatic catalysis or to stabilize protein structure, we hypothesize that Zn2+ might function analogously to Na+ in class A GPCRs, where the metal ion allosterically tunes agonist affinity and effector interactions by biasing receptor conformation (Gerwert et al., 2014; Katritch et al., 2014; Massink et al., 2015). In high-resolution x-ray structures of inactive class A GPCRs (β2AR, PDB accession no. 4BVN; A2AR, PDB accession no. 4EIY; δOR, PDB accession no. 4N6H), Na+ exhibits trigonal bipyramidal coordination geometry. The spontaneous conversion from octahedral to trigonal bipyramidal coordination seen in microsecond MD simulations of GPCRs (Shang et al., 2014) is consistent with the hypothesis that dynamic changes in the positions of Na+ and its coordinating protein atoms are correlated with conformational changes that alter ligand affinity and effector coupling (Gerwert et al., 2014; Katritch et al., 2014). Additional experimental and computational studies are needed to investigate how subtle differences in Zn2+ position or coordination geometries influence Hv1 function, and our results provide a foundation for such future work.

Changes in Zn2+ occupancy at the primary His2-Glu coordination site observed during MD simulation are correlated with structural changes in other regions of the Hv1 VS domain, suggesting the existence of structural pathways for long-range allosteric coupling. Prior to the replacement of H2.40-Nδ1 and H3.71-Nδ1 in the first solvation shell with water oxygens (at t = 110 ns in Hv1F·Zn2+ c0; Fig. 5 H), we find dynamic changes in both side-chain (H-bonds and salt bridges) and backbone structures (Fig. 5, G and H; and Fig. S8, A–G). Additional functional and computational studies will be needed to further examine Zn2+-dependent conformational changes and determine the molecular pathways by which energy is transduced in the Hv1 VS domain. Similar to other VS domain proteins (Swartz, 2008; Lacroix et al., 2014), neutralization of ICN residues in Hv1 (i.e., E153A or D174A) causes the position of the GAQ-V relation in Hv1 to shift negatively (Ramsey et al., 2010; Chamberlin et al., 2014), indicating that these side chains normally help to stabilize a resting-state conformation. Comparing MD simulations of Zn2+-liganded and unliganded conformations illustrates that the pattern of Coulombic interactions and H-bonds between highly conserved acidic “counter-charge” residues (E2.53 and D3.50) and S4 Arg “gating charge” (R2/R4.50 and R3/R4.53) in the ICN are sensitive to Zn2+ occupancy (Figs. 6, S6, and S8, and Tables S4 and S5). Although R1 and D1.51 form a stable pair in both Hv1 F·Zn2+ and Hv1 F + Zn2+ MD simulations, D3.50 exhibits an R3/R2 and E2.53 selectively makes a stable salt bridge with R3 when Zn2+ is constrained (Fig. 6). Differences in ICN structure are consistent with the hypothesis that Zn2+ occupancy is associated with an overall strengthening of interactions that stabilize a GAQ-closed, resting-state conformation in Hv1 (Fig. 6 and Tables S4 and S5).

An analysis of coupled salt bridge and Zn2+-ligand distance changes in Zn2+-constrained (Hv1 F·Zn2+ c5) and unconstrained (Hv1 F·Zn2+ c0) MD simulations further supports the idea that Zn2+ occupancy of the primary extracellular His2-Glu binding site is allosterically coupled to changes in ICN structure. We therefore hypothesize that MD simulations elaborate a conformational coupling mechanism for Zn2+ modulation of GAQ gating in Hv1. The effect of E2.53 neutralization to decrease Zn2+ potency shows that changes in ICN structure reciprocally alter Zn2+ coordination, further supporting the conformational coupling hypothesis. The allosteric interaction between extracellular Zn2+ and ICN residues occurs over a long distance: E2.53-Cα is 19.9 ± 0.3 Å or 24.1 ± 0.6 Å (mean ± SD during 120 ns Hv1 F·Zn2+ c0 MD simulation) away from H2.40-Cα or H3.71-Cα, respectively. Long-range allosteric coupling provides a structural explanation for a biophysical phenomenon that has remained unexplained since Zn2+-dependent modulation of native voltage-gated H+ currents was first measured (Thomas and Meech, 1982; Barish and Baud, 1984; Cherny and DeCoursey, 1999). Additional studies are needed to more precisely understand what causes the ICN conformational differences reported here, but we hypothesize that local structural changes at the His2-Glu site and ICN may be readily propagated to remote sites via relatively subtle changes in side-chain chemistry (i.e., pKa shifts at ICN side-chain atoms) that alter the strength of H-bonds and Coulombic interactions in an extended electrostatic network. Long-range electrostatic networks have been previously reported for GPCRs (Isom and Dohlman, 2015; Miao et al., 2015). Finally, we speculate that extracellular “gating modifier toxins” may elicit changes in ICN structure to alter VS activation gating in other VS domain proteins (Rogers et al., 1996; Alabi et al., 2007; Milescu et al., 2009; Salari et al., 2016) by a conformational coupling mechanism that is similar to the effect of Zn2+ on GAQ gating in Hv1 postulated here.

The authors wish to thank Philip Mosier and Carlos A. Villalba-Galea for useful discussions and technical assistance.

This work was supported by a Consejo Nacional de Ciencia y Tecnología (CONACyT) postdoctoral fellowship 251376 to V. De La Rosa and National Institutes of Health grant GM092908 to I.S. Ramsey.

The authors declare no competing financial interests.

Author contributions: V. De La Rosa designed, performed, and analyzed experimental data; A.L. Bennett designed, conducted, and analyzed molecular modeling and simulations; I.S. Ramsey designed experiments, analyzed electrophysiological and computational data, prepared figures, and wrote the manuscript.

Richard W. Aldrich served as editor.

Alabi
,
A.A.
,
M.I.
Bahamonde
,
H.J.
Jung
,
J.I.
Kim
, and
K.J.
Swartz
.
2007
.
Portability of paddle motif function and pharmacology in voltage sensors
.
Nature.
450
:
370
375
.
Barish
,
M.E.
, and
C.
Baud
.
1984
.
A voltage-gated hydrogen ion current in the oocyte membrane of the axolotl, Ambystoma
.
J. Physiol.
352
:
243
263
.
Best
,
R.B.
,
X.
Zhu
,
J.
Shim
,
P.E.
Lopes
,
J.
Mittal
,
M.
Feig
, and
A.D.
Mackerell
Jr
.
2012
.
Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles
.
J. Chem. Theory Comput.
8
:
3257
3273
.
Burgess
,
J.
1978
.
Metal ions in solution.
Ellis Horwood
,
Chichester
.
481
pp.
Byerly
,
L.
, and
Y.
Suen
.
1989
.
Characterization of proton currents in neurones of the snail, Lymnaea stagnalis
.
J. Physiol.
413
:
75
89
.
Chamberlin
,
A.
,
F.
Qiu
,
S.
Rebolledo
,
Y.
Wang
,
S.Y.
Noskov
, and
H.P.
.
2014
.
Hydrophobic plug functions as a gate in voltage-gated proton channels
.
111
:
E273
E282
.
Cherny
,
V.V.
, and
T.E.
DeCoursey
.
1999
.
pH-dependent inhibition of voltage-gated H(+) currents in rat alveolar epithelial cells by Zn2+ and other divalent cations
.
J. Gen. Physiol.
114
:
819
838
.
Dudev
,
T.
,
Y.L.
Lin
,
M.
Dudev
, and
C.
Lim
.
2003
.
First-second shell interactions in metal binding sites in proteins: a PDB survey and DFT/CDM calculations
.
J. Am. Chem. Soc.
125
:
3168
3180
.
Eargle
,
J.
,
D.
Wright
, and
Z.
Luthey-Schulten
.
2006
.
Multiple Alignment of protein structures and sequences for VMD
.
Bioinformatics.
22
:
504
506
.
Gerwert
,
K.
,
E.
Freier
, and
S.
Wolf
.
2014
.
The role of protein-bound water molecules in microbial rhodopsins
.
Biochim. Biophys. Acta.
1837
:
606
613
.
Hitoshi
,
O.
,
Y.
Toshio
, and
M.
Masunobu
.
1976
.
X-Ray Diffraction Studies of the Structures of Hydrated Divalent Transition-Metal Ions in Aqueous Solution
.
Bull. Chem. Soc. Jpn.
49
:
701
708
.
Huang
,
J.
,
S.
Rauscher
,
G.
Nawrocki
,
T.
Ran
,
M.
Feig
,
B.L.
de Groot
,
H.
Grubmüller
, and
A.D.
MacKerell
Jr
.
2017
.
CHARMM36m: an improved force field for folded and intrinsically disordered proteins
.
Nat. Methods.
14
:
71
73
.
Humphrey
,
W.
,
A.
Dalke
, and
K.
Schulten
.
1996
. VMD: visual molecular dynamics. J. Mol. Graph. 14:33-8, 27-8.
Isom
,
D.G.
, and
H.G.
Dohlman
.
2015
.
Buried ionizable networks are an ancient hallmark of G protein-coupled receptor activation
.
112
:
5702
5707
.
Johansson
,
G.
1992
.
Structures of Complexes in Solution Derived from X-Ray Diffraction Measurements
.
39
:
159
232
.
Jorgensen
,
W.L.
,
J.
Chandrasekhar
,
J.D.
,
R.W.
Impey
, and
M.L.
Klein
.
1983
.
Comparison of simple potential functions for simulating liquid water
.
J. Chem. Phys.
79
:
926
935
.
Katritch
,
V.
,
G.
Fenalti
,
E.E.
Abola
,
B.L.
Roth
,
V.
Cherezov
, and
R.C.
Stevens
.
2014
.
Allosteric sodium in class A GPCR signaling
.
Trends Biochem. Sci.
39
:
233
244
.
Klauda
,
J.B.
,
R.M.
Venable
,
J.A.
Freites
,
J.W.
O’Connor
,
D.J.
Tobias
,
C.
Mondragon-Ramirez
,
I.
Vorobyov
,
A.D.
MacKerell
Jr
., and
R.W.
Pastor
.
2010
.
Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types
.
J. Phys. Chem. B.
114
:
7830
7843
.
Klauda
,
J.B.
,
V.
Monje
,
T.
Kim
, and
W.
Im
.
2012
.
Improving the CHARMM force field for polyunsaturated fatty acid chains
.
J. Phys. Chem. B.
116
:
9424
9431
.
Lacroix
,
J.J.
,
H.C.
Hyde
,
F.V.
Campos
, and
F.
Bezanilla
.
2014
.
Moving gating charges through the gating pore in a Kv channel voltage sensor
.
111
:
E1950
E1959
.
Laitaoja
,
M.
,
J.
Valjakka
, and
J.
Jänis
.
2013
.
Zinc coordination spheres in protein structures
.
Inorg. Chem.
52
:
10983
10991
.
Li
,
Q.
,
R.
Shen
,
J.S.
Treger
,
S.S.
Wanderling
,
W.
Milewski
,
K.
Siwowska
,
F.
Bezanilla
, and
E.
Perozo
.
2015
.
Resting state of the human proton channel dimer in a lipid bilayer
.
112
:
E5926
E5935
.
Li
,
S.J.
,
Q.
Zhao
,
Q.
Zhou
,
H.
Unno
,
Y.
Zhai
, and
F.
Sun
.
2010
.
The role and structure of the carboxyl-terminal domain of the human voltage-gated proton channel Hv1
.
J. Biol. Chem.
285
:
12047
12054
.
Lishko
,
P.V.
, and
Y.
Kirichok
.
2010
.
The role of Hv1 and CatSper channels in sperm activation
.
J. Physiol.
588
:
4667
4672
.
Lishko
,
P.V.
,
I.L.
Botchkina
,
A.
Fedorenko
, and
Y.
Kirichok
.
2010
.
Acid extrusion from human spermatozoa is mediated by flagellar voltage-gated proton channel
.
Cell.
140
:
327
337
.
MacKerell
,
A.D.
,
D.
Bashford
,
M.
Bellott
,
R.L.
Dunbrack
,
J.D.
Evanseck
,
M.J.
Field
,
S.
Fischer
,
J.
Gao
,
H.
Guo
,
S.
Ha
, et al
1998
.
All-atom empirical potential for molecular modeling and dynamics studies of proteins
.
J. Phys. Chem. B.
102
:
3586
3616
.
MacKerell
,
A.D.
Jr
.,
M.
Feig
, and
C.L.
Brooks
III
.
2004
.
Improved treatment of the protein backbone in empirical force fields
.
J. Am. Chem. Soc.
126
:
698
699
.
Massink
,
A.
,
H.
Gutiérrez-de-Terán
,
E.B.
,
N.V.
Ortiz Zacarías
,
L.
Xia
,
L.H.
Heitman
,
V.
Katritch
,
R.C.
Stevens
, and
A.P.
IJzerman
.
2015
.
Sodium ion binding pocket mutations and adenosine A2A receptor function
.
Mol. Pharmacol.
87
:
305
313
.
Meech
,
R.W.
, and
R.C.
Thomas
.
1987
.
Voltage-dependent intracellular pH in Helix aspersa neurones
.
J. Physiol.
390
:
433
452
.
Miao
,
Y.
,
A.D.
Caliman
, and
J.A.
McCammon
.
2015
.
Allosteric effects of sodium ion binding on activation of the m3 muscarinic g-protein-coupled receptor
.
Biophys. J.
108
:
1796
1806
.
Milescu
,
M.
,
F.
Bosmans
,
S.
Lee
,
A.A.
Alabi
,
J.I.
Kim
, and
K.J.
Swartz
.
2009
.
Interactions between lipids and voltage sensor paddles detected with tarantula toxins
.
Nat. Struct. Mol. Biol.
16
:
1080
1085
.
Musset
,
B.
,
S.M.
Smith
,
S.
Rajan
,
V.V.
Cherny
,
S.
Sujai
,
D.
Morgan
, and
T.E.
DeCoursey
.
2010
.
Zinc inhibition of monomeric and dimeric proton channels suggests cooperative gating
.
J. Physiol.
588
:
1435
1449
.
Qiu
,
F.
,
A.
Chamberlin
,
B.M.
Watkins
,
A.
Ionescu
,
M.E.
Perez
,
R.
Barro-Soria
,
C.
González
,
S.Y.
Noskov
, and
H.P.
.
2016
.
Molecular mechanism of Zn2+ inhibition of a voltage-gated proton channel
.
113
:
E5962
E5971
.
,
B.
, and
T.L.
Leto
.
2008
.
Oxidative innate immune defenses by Nox/Duox family NADPH oxidases
.
Contrib. Microbiol.
15
:
164
187
.
Ramsey
,
I.S.
,
M.M.
Moran
,
J.A.
Chong
, and
D.E.
Clapham
.
2006
.
A voltage-gated proton-selective channel lacking the pore domain
.
Nature.
440
:
1213
1216
.
Ramsey
,
I.S.
,
Y.
Mokrab
,
I.
Carvacho
,
Z.A.
Sands
,
M.S.P.
Sansom
, and
D.E.
Clapham
.
2010
.
An aqueous H+ permeation pathway in the voltage-gated proton channel Hv1
.
Nat. Struct. Mol. Biol.
17
:
869
875
.
Randolph
,
A.L.
,
Y.
Mokrab
,
A.L.
Bennett
,
M.S.
Sansom
, and
I.S.
Ramsey
.
2016
.
Proton currents constrain structural models of voltage sensor activation
.
eLife.
5
:
e18017
.
Roberts
,
E.
,
J.
Eargle
,
D.
Wright
, and
Z.
Luthey-Schulten
.
2006
.
MultiSeq: unifying sequence and structure data for evolutionary analysis
.
BMC Bioinformatics.
7
:
382
.
Rogers
,
J.C.
,
Y.
Qu
,
T.N.
,
T.
Scheuer
, and
W.A.
Catterall
.
1996
.
Molecular determinants of high affinity binding of alpha-scorpion toxin and sea anemone toxin in the S3-S4 extracellular loop in domain IV of the Na+ channel alpha subunit
.
J. Biol. Chem.
271
:
15950
15962
.
Salari
,
A.
,
B.S.
Vega
,
L.S.
Milescu
, and
M.
Milescu
.
2016
.
Molecular Interactions between Tarantula Toxins and Low-Voltage-Activated Calcium Channels
.
Sci. Rep.
6
:
23894
.
Sasaki
,
M.
,
M.
Takagi
, and
Y.
Okamura
.
2006
.
A voltage sensor-domain protein is a voltage-gated proton channel
.
Science.
312
:
589
592
.
Shang
,
Y.
,
V.
LeRouzic
,
S.
Schneider
,
P.
Bisignano
,
G.W.
Pasternak
, and
M.
Filizola
.
2014
.
Mechanistic insights into the allosteric modulation of opioid receptors by sodium ions
.
Biochemistry.
53
:
5140
5149
.
Stone
,
J.
1998
. An efficient library for parallel ray tracing and animation. MS thesis. University of Missouri, Rolla, MO. 79 pp.
Stote
,
R.H.
, and
M.
Karplus
.
1995
.
Zinc binding in proteins and solution: a simple but accurate nonbonded representation
.
Proteins.
23
:
12
31
.
Swartz
,
K.J.
2008
.
Sensing voltage across lipid membranes
.
Nature.
456
:
891
897
.
Takeshita
,
K.
,
S.
Sakata
,
E.
Yamashita
,
Y.
Fujiwara
,
A.
Kawanabe
,
T.
Kurokawa
,
Y.
Okochi
,
M.
Matsuda
,
H.
Narita
,
Y.
Okamura
, and
A.
Nakagawa
.
2014
.
X-ray crystal structure of voltage-gated proton channel
.
Nat. Struct. Mol. Biol.
21
:
352
357
.
Taylor
,
A.R.
,
A.
Chrachri
,
G.
Wheeler
,
H.
Goddard
, and
C.
Brownlee
.
2011
.
A voltage-gated H+ channel underlying pH homeostasis in calcifying coccolithophores
.
PLoS Biol.
9
:
e1001085
.
Thomas
,
R.C.
, and
R.W.
Meech
.
1982
.
Hydrogen ion currents and intracellular pH in depolarized voltage-clamped snail neurones
.
Nature.
299
:
826
828
.
Villalba-Galea
,
C.A.
2014
.
Hv1 proton channel opening is preceded by a voltage-independent transition
.
Biophys. J.
107
:
1564
1572
.
Wacker
,
D.
,
R.C.
Stevens
, and
B.L.
Roth
.
2017
.
How Ligands Illuminate GPCR Molecular Pharmacology
.
Cell.
170
:
414
427
.
Webb
,
B.
, and
A.
Sali
.
2016
. Comparative Protein Structure Modeling Using MODELLER. Curr. Protoc. Bioinformatics. 54:5.6.1–5.6.37.
Yu
,
Z.
,
P.
Li
, and
K.M.
Merz
Jr.
2018
.
Extended Zinc AMBER Force Field (EZAFF)
.
J. Chem. Theory Comput.
14
:
242
254
.

## Author notes

*

V. De La Rosa and A.L. Bennett contributed equally to this paper.

V. De La Rosa’s present address is Dept. of Cellular and Integrative Physiology, University of Texas Health, San Antonio, TX.