The connexin26 human mutation N14K disrupts cytosolic intersubunit interactions and promotes channel opening

Connexins are a family of transmembrane proteins encoded by 21 human genes, found in almost every human tissue, which play key roles in physiology and pathology (Söhl and Willecke, 2003). At the structural level, the basic topology of a connexin protein consists of four transmembrane segments (TM1–4), two extracellular loops (E1 and E2), and one cytoplasmic loop (CL), with both the N- and C-terminal domains facing the cytosol. The assembly of six connexin proteins forms a hemichannel that is trafficked to the plasma membrane, where it can dock with another hemichannel of an adjacent cell to form a gap junction channel (GJC; Goodenough et al., 1996; Gaietta et al., 2002; Sosinsky and Nicholson, 2005). Both unpaired hemichannels and GJCs are permeable to atomic ions and small molecules (Harris, 2001). Hemichannels at the plasma membrane are mostly closed when they are not part of GJCs. This is achieved by negative membrane potential and normal extracellular Ca²⁺ concentrations that together significantly reduce hemichannel open probability (Ebihara et al., 2003; Bukauskas and Verselis, 2004; Fasciani et al., 2013; Verselis and Srinivas, 2013; Bargiello et al., 2018).

The mechanisms by which extracellular Ca²⁺ controls gating in connexin hemichannels is related to interactions of Ca²⁺ ions with negatively charged residues lining the extracellular entrance of the channel pore (near the border of the TM1 and E1 domains). The first published crystals of Cx26 GJC structure and subsequent MD studies indicated a highly conserved electrostatic network located at the extracellular pore entrance, where negatively charged residues at the TM1/E1 region interact with positively charged residues of the E1, TM2, and TM4 domains in the absence of Ca²⁺ (Maeda et al., 2009; Kwon et al., 2012). A more recent crystal structure of Cx26 GJCs in the presence of Ca²⁺ suggests that Ca²⁺ binds at the TM1/E1 region (Bennett et al., 2016), forming an electrostatic ring that restricts cation movement. Functional studies of Cx26 hemichannels reveal that extracellular Ca²⁺ disrupts TM1/E1 electrostatic networks by interacting with negatively charged residues lining the pore, favoring the closed conformation of the channel in an allosteric manner (Lopez et al., 2016). Calcium occupancy of the gating ring, however, did not prevent access of ions or small molecules deeper into the pore, suggesting that the gate itself is intracellular to the Ca²⁺ gating ring (Lopez et al., 2016).

In contrast to Ca²⁺ regulation, the mechanisms governing voltage sensitivity, and specifically the location of the voltage sensor in connexin hemichannels, remain elusive. It is thought that negative voltages interact with the electrostatic network at the extracellular entrances, particularly at the TM1/E1 region, to keep the hemichannels closed. However, the conformational changes that accompany hemichannel closing by this means are not well understood and may involve changes in other parts of the channels, including the cytoplasmic end of the pore (e.g., the N-terminal [NT] domain; Kwon et al., 2013; Bargiello et al., 2018). Because the NT domain is folded into the intracellular pore entrance and forms hydrophobic interactions with residues of TM1, it has been proposed to act as a gate-type “plug” whose movement is dependent on membrane voltage changes (Oshima et al., 2007; Maeda et al., 2009).

Several human connexin mutations that cause disease promote gain of hemichannel function. In Cx26 channels, gain-of-function mutations produce syndromic deafness and skin disorders (Gerido et al., 2007; Lee et al., 2009; Sánchez et al., 2010; Levit et al., 2012; Nielsen et al., 2012). Interestingly, they are all located either in the NT or near the TM1/E1, the latter of which, as described previously, is a critical region for Ca²⁺ and voltage sensing. By investigating the mechanisms of how mutations at TM1/E1 promote gain in hemichannel function, we and others were able to gain valuable information on how Ca²⁺ gates hemichannels (Lopez et al., 2013, 2014, 2016; Sanchez et al., 2013). Recent studies report that mutations at the N-terminal region affect the voltage component of regulation by Ca²⁺ (Sanchez et al., 2016) as well as a distinct voltage-sensitive process involving the NT itself (Sanchez et al., 2016; García et al., 2018). Here, we investigate the molecular mechanisms by which the syndromic deafness N14K mutation promotes gain-of-function in the Cx26 hemichannel. By combining electrophysiological techniques, double mutant cycle analysis, and MD simulations, we report that the N14K mutation increases hemichannel open probability by disrupting the interaction of the NT with residues in the TM2/CL transition in the adjacent connexin subunit. This strongly suggests that NT–TM2/CL interactions promote Cx26 hemichannel closing.

Electrophysiological data were collected using the two-electrode voltage-clamp technique. All recordings were made at room temperature (20–22°C). The recording solutions contained (in mM) 118 NaCl, 2 KCl, and 5 HEPES, pH 7.4, with a range of Ca²⁺ or Ba²⁺ concentration from 0.01 to 10. Currents from oocytes were recorded 1–3 d after cRNA injection using an oocyte clamp (OC-725C; Warner Instruments). Currents were sampled at 2 kHz and low-pass filtered at 0.2 kHz. Microelectrode resistances were between 0.1 and 1.2 MΩ when filled with 3 M KCl.

Single-channel recordings were done using excised inside-out patch clamp from Xenopus laevis oocytes. Oocytes were bathed in internal patch solution containing (in mM) 140 KCl, 1.8 CaCl₂, 1 MgCl₂, 5 HEPES, and 3 EGTA with pH adjusted to 8 with KOH. Patch micropipettes made of glass borosilicate capillaries (World Precision Instruments) were pulled on a horizontal puller (Sutter Instruments), and fire polished using a custom microforge to a resistance ranging between 1 and 10 MΩ when filled with internal patch solution. Patch clamp was performed using an Axopatch 200B amplifier (Axon Instruments, Molecular Devices). Data were digitized at 1 kHz using a 16-bit A/D converter data acquisition system (Digidata1322; Axon Instruments, Molecular Devices) and acquired using Clampex 8 acquisition software (Axon Instruments, Molecular Devices).

Xenopus oocytes express endogenous Cx38 hemichannels, which are activated at extracellular Ca²⁺ concentrations less than ∼0.15 mM. These channels may interfere with the analysis of the heterologously expressed human Cx26 (hCx26) WT and mutant currents. Antisense oligonucleotides against Cx38 were injected in each oocyte to reduce the expression of Cx38 4 h after harvesting the oocytes (1 mg/ml; using the sequence from Ebihara [1996]). The same oocytes were coinjected 1 d later with 9–23 nl complementary RNA (0.5–1 mg/ml) hCx26 or mutants and additional Cx38 (1 mg/ml) antisense. Oocytes injected only with Cx38 antisense were tested at low Ca²⁺ concentration (0.01 mM), and only batches in which no or low levels of endogenous Cx38 currents were observed were used to calculate the Ca²⁺ or Ba²⁺ dose–response curves. Ca²⁺ and Ba²⁺ dose–response measurements were obtained by assessing the peak tail currents after reaching current saturation during a depolarizing pulse from −80 to 0 mV as described by Lopez et al. (2013). Deactivation time constants from tail currents were fitted to the tail current 10 s after reaching steady state, to exponential functions using Clampfit 11 software (Molecular Devices).

Membrane holding potential was set at −80 mV followed by 20-s depolarizing test pulses in 10-mV increments. The peak tail currents were measured for each step and normalized to the maximal current observed at 60 mV. Normalized currents were plotted as function of voltage. The data points from this graph were fitted to a two-state Boltzmann function:

where V is the prepulse potential, and V_1/2 is the potential when the probability for channel opening is 50%.

Double mutant cycle analysis was performed using V_1/2 obtained from the voltage dependence of each channel. These macroscopic parameters have been used previously to explore side-chain interactions and establish coupling coefficients in ligand-gated channels (Kash et al., 2003; Price et al., 2007; Gleitsman et al., 2008). The V_1/2 and z values from the two-state Boltzmann fit were used to estimate changes in Gibbs free energy (ΔG) for channel activation and availability assuming a two-state transition (from closed to open) at test potential V using the equation

Coupling energy (ΔΔG) using the voltage dependence ΔG values was calculated as

If two mutations are functionally independent with respect to their voltage dependence, the coupling energy will be close to 0 kcal/mol. Significant coupling is indicated by any value that deviates from zero, but an accepted cutoff from nonadditivity is 0.5 kcal/mol (Laha and Wagner, 2011). All measurements were obtained at 1.8 mM Ca²⁺.

The WT and N14K simulation systems were based on the protocol of Luo et al. (2016), which included the equilibrated structure of the Cx26 in explicit membrane and solvent. Molecular Operating Environment (MOE; Chemical Computing Group) was used to introduce N14K mutation. CHARMM-GUI (Jo et al., 2008) was used to construct each system, including membrane and solvent. The system was initially neutralized with Cl⁻, followed by addition of K⁺ and Cl⁻ to yield a concentration of 0.1 M KCl. Each system was solvated in a fully hydrated 1-palmitoyl-2-oleoylphosphatidylcholine bilayer with a 15-Å layer of explicit water above and below the bilayer. CHARMM parameter sets were used throughout, including the CHARMM22/CMAP parameters for the protein and ions (MacKerell et al., 1998; Mackerell et al., 2004), the additive C36 all-atom parameters for lipids (Klauda et al., 2010), and TIP3P for water (Jorgensen et al., 1983). All simulations were performed using the PMEMD module of the AMBER16 package (Case et al., 2016). Orthorhombic periodic boundary conditions were used for all simulations in the isobaric-isothermal ensemble using Langevin thermostat and Berendsen barostat. The pressure and temperature were maintained at 1 atm and 303.15°K. Long-range electrostatic interactions were treated using the default particle mesh Ewald method. Nonbonded cutoff was set to be 8 Å. The dynamics were propagated using Langevin dynamics with Langevin damping coefficient of 1 ps⁻¹ and a time step of 2 fs. The SHAKE algorithm was applied to all hydrogen atoms. Both WT and N14K systems were simulated for three independent 300-ns periods starting from different initial velocities, and the last 100 ns were used for data analysis.

Pairwise per-residue-based GB energy decomposition is done using idecomp = 3 (adding one to four interactions to internal energy), igb = 5, saltcon = 0.150 options in MMPBSA.py version 14.0 (Miller et al., 2012). The total pairwise energy is the sum of a van der Waals term, an electrostatic term, a GB polar solvation term, and a nonpolar solvation term. The nonpolar solvation term is calculated using surface areas (gbsa = 2) based on a linear combination of pairwise overlaps algorithm. Each system was simulated for four replicas of 300 ns using different initial velocities. 100 frames (1-ns stride) extracted from the last 100 ns of 300-ns trajectory were analyzed for the interaction between residues 13–15 and 226 residues on the adjacent subunits. The statistical analysis is detailed in the Results section.

The contact correlation maps were computed for each α carbon pair using the last 100 ns with 200-ps stride from each replica using the program CARMA (Glykos, 2006). Pearson correlation coefficients between residue pairs were measured by using

where X_i is the fluctuation vector of residue i, defined as the deviation from the average coordinates. The statistical analysis is detailed in the Results section.

Fig. S1 shows that BAPTA does not eliminate the multiphasic tail current deactivations in Ca²⁺ in N14K channels. Fig. S2 shows currents for Cx26 WT, H100A, N14K, and N14K H100A after a depolarizing step from −80 to 0 mV. Fig. S3 corresponds to single-channel recordings for Cx26 H100A or E101A mutant hemichannels. Fig. S4 shows the wavelet map and RMSD and RMSF plots of WT and N14K simulations. Table S1 is the summary of the pairwise residue interaction energy data for the 29 pairs, with a P value ∼0 from KS test. Table S2 is the summary of the correlation analysis data for the 25 pairs, with a P value <0.005 from KS test.

Because Cx26 hemichannels are gated by changes in voltage and extracellular Ca²⁺, we examined how the N14K mutation affected normal gating functions in response to these stimuli. hCx26 WT and N14K mutant channels were expressed in Xenopus laevis oocytes, and hemichannel currents were assessed using the two-electrode voltage-clamp technique. We used our previously standardized protocols to explore regulation by external Ca²⁺ (Lopez et al., 2013), which focus on changes of peak tail currents and their relaxation kinetics after a 40-s voltage pulse from –80 to 0 mV at different Ca²⁺ concentrations. We did not step to more positive voltages to minimize activation of endogenous currents and the deleterious effects of massive hemichannel opening during the long 40-s pulses required to reach steady-state activation. We have shown that the peak tail currents and relaxation kinetics increase with reduction of external Ca²⁺ (Lopez et al., 2014). The peak tail currents reflect steady-state channel activation at 0 mV at any given Ca²⁺ concentration, and the relaxation of the tail currents reflect the kinetics of the open-to-closed transitions at −80 mV.

Fig. 1 (A and B) shows current traces obtained at 1.8, 0.5, and 0.01 mM extracellular Ca²⁺ from an oocyte expressing either hCx26 WT (A) or N14K (B) hemichannels. For both WT and N14K mutant, reduction of extracellular Ca²⁺ results in an increase of the peak tail current and a slowing of the deactivation kinetics, the latter illustrated in the insets showing currents normalized to the maximal tail current. Increased hemichannel activity as Ca²⁺ is reduced is also reflected in the increased “leak”/holding currents observed before the depolarizing steps at the holding potential of −80 mV (denoted by the red arrows). From the traces, it is clear that relative to WT, the N14K mutant has greater leak current at −80 mV and dramatically slowed deactivation. Interestingly, deactivation kinetics for the mutant N14K does not follow a monoexponential decay; instead it presents a multiphasic deactivation. One possibility is that the complexity of the decay is secondary to Ca²⁺ influx through the greater number of open N14K channels, which could activate endogenous Ca²⁺-activated chloride currents, contaminating the measured macroscopic currents. Previously, we found that such contaminating currents can be eliminated by preinjection of oocytes with ∼120 µM BAPTA (Lopez et al., 2016); all subsequent experiments were performed after preinjection with BAPTA. Strikingly, however, we found that N14K hemichannel deactivation kinetics were unchanged by this treatment (Fig. S1), suggesting an intrinsic multiphasic behavior.

By normalizing the peak tail current at each Ca²⁺ concentration to the maximal activation at the lowest Ca²⁺ concentration, we were able to plot the activation of the hemichannels as a function of the Ca²⁺ concentration (Fig. 1 C). The data then were fitted to a Hill equation of the form

where the fractional current is I/I_max, I is the tail current at a particular Ca²⁺ concentration, I_max is the maximal tail current activation at 0.01 mM Ca²⁺, K_D is the apparent affinity, [Ca²⁺] is the concentration of Ca²⁺ applied to the bath, and n is the Hill coefficient. The best-fit parameter values for K_D and n are 0.35 ± 0.03 mM and 1.43 ± 0.15 for the WT and 0.74 ± 0.05 mM and 1.55 ± 0.12 for N14K, respectively. The mutant’s rightward shift in the K_D without significant change in slope indicates that substitution of a lysine at position N14 disfavors the closed state induced by extracellular Ca²⁺, without effect on the Ca²⁺ sensitivity itself.

Next, we compared changes in voltage-dependent regulation of WT and N14K mutant hemichannels estimated from the current–voltage relationships (Fig. 2). To obtain this relationship, we kept the holding potential at −80 mV and applied 20-s voltage steps from −100 to +70 mV in 10-mV increments. These measurements were done at physiological Ca²⁺ concentration, and peak tail currents were used to obtain voltage dependence. Data were fitted using a two-state Boltzmann equation (Fig. 2 B, solid lines). The best-fit parameter values for V_1/2 and z are 9.6 ± 0.9 and 2.4 for WT (red) and −15.8 ± 0.9 and 0.9 for N14K mutant (blue). Importantly, hCx26 hemichannels activate only at potentials more positive than −20 mV (Fig. 2 B, red circles). In contrast, oocytes injected with N14K show hemichannel currents at potentials as negative as −90 mV (Fig. 2 B, blue circles). Thus, the mutant hemichannel opens at large negative voltages, at which the WT channel is normally closed. The change in slope of the I-V relations suggest that, at a minimum, the mutation affects voltage sensitivity and/or the pathway by which voltage regulates gating.

We previously used analysis of deactivation kinetics and steady-state data to assess the effect of connexin mutations on hemichannel gating (Lopez et al., 2013, 2014, 2016). However, the multiphasic nature of the deactivation currents of the mutant N14K in the presence of extracellular Ca²⁺ made this analysis difficult. We evaluated the effect of several other divalent cations and found that when Ca²⁺ was replaced with Ba²⁺, we were able to retain the effects seen in Ca²⁺ on leak and steady-state currents, but with monotonic deactivation kinetics (Fig. 3 A). The Ba²⁺ dose–response relations for mutant N14K and WT hemichannels (Fig. 3 C, left) showed a shift similar to that observed in Ca²⁺. The calculated values for the apparent K_D and n are 0.407 ± 0.032 mM and 1.00 ± 0.08 for hCx26 WT and 0.795 ± 0.175 mM and 0.71 ± 0.11 for hCx26 N14K mutant, respectively. To assess hemichannel closure, we quantified and compared the deactivation time constants of the tail currents as a function of Ba²⁺ concentration (Fig. 3 C, right). Deactivation kinetics for both hCx26 WT and N14K mutant are accelerated as Ba²⁺ concentrations increase. More importantly, deactivation time constant values for the N14K hemichannels indicate that they close significantly more slowly compared with WT. This suggests that the N14K mutation affects the transition from the open to closed states.

To further quantify this effect, we performed single-channel recordings for both hCx26 WT and N14K. The hemichannel open probabilities at different voltages were determined for hCx26 WT and N14K (Fig. 4 A). At negative potentials, WT Cx26 hemichannels are mostly closed, whereas N14K mutant hemichannels are open. At −40 mV, for example, WT Cx26 hemichannels had Po = 0.02 ± 0.005, whereas mutant N14K channels had Po = 0.8 ± 0.05. Linear regression analysis of the single-channel conductance of mutant N14K hemichannels does not differ from WT Cx26 hemichannels. Our single-channel data demonstrate that the N14K mutation increases open probability, likely by stabilizing the hemichannel open conformation.

Previous work suggested that N14K forms a π-cation interaction with H100 located at the transition between the TM2 and CL domains of the adjacent connexin subunit. This interaction was proposed to be responsible for the stabilization of the hemichannel in the open conformation (Sanchez et al., 2016). We decided to test the suggested interaction between N14K with H100 by following the same experimental design as in Fig. 3. If indeed N14K stabilizes the hemichannel in the open conformation, then replacement of H100 with an Ala should disrupt the interaction, and deactivation kinetics similar to WT should be observed. In Fig. 5, the deactivation kinetics for all three mutants and WT channels are plotted versus Ba²⁺ concentrations (left graph). The normalized deactivation currents for all four constructs are shown at 5 mM Ba²⁺ (right). Strikingly, the deactivation kinetics of the double mutant N14K/H100A were quite similar to those of N14K, ruling out that a π-cation interaction between these residues stabilizes the open channel conformation. In addition, the deactivation kinetics of the double mutant in the presence of Ca²⁺ resemble those of N14K, in further support that N14K does not interact with H100 by a π-cation interaction (Fig. S2).

To further validate our findings, we performed single-channel recordings in double mutant N14K/H100A hemichannels and compared them to those obtained in WT and N14K hemichannels. Single-hemichannel current traces from double mutant N14K/H100A display similar behavior, over a range of voltages, to that observed in single N14K mutant hemichannels. Indeed, open channel probability for N14K/H100A is similar to that found for N14K, but not WT (Fig. 6; compare with Fig. 4). These data confirm that H100 is not critical to stabilize the open channel conformation of N14K.

To investigate the mechanisms by which the mutant N14K increases open channel probability, we first compared pairwise residue interaction energies between WT and N14K using MD simulations. The simulation system was either a WT or N14K hemichannel in explicit membrane and solvent, run for four independent replicas of 300 ns. The pairwise per-residue-based GB energy decompositions between residues 13–15 (in the NT) and rest of the protein residues were derived from the last 100 ns of each trajectory using 1-ns stride (frequency of saving the coordinates). Because sixfold symmetry of the hexamer is not imposed during the MD simulation, the residue–residue interactions between different subunit pairs are not expected to be identical. This is because although the simulations have reached apparent equilibrium (that is, the conformational ensemble is not changing over time), each subunit continues to sample different energies and conformations within the ensemble during the 300-ns simulations. To obtain 100% symmetry across all subunits, much longer simulations would be needed to ensure that each subunit samples the full conformational space so that all the fluctuations would be averaged out. Thus, to avoid misleading results from the motions of a single subunit or subunit pair within the current submicrosecond simulation timescale, we performed an analysis of the pairwise distributions of GB energies of all six subunits, for WT and for N14K mutant hemichannel systems.

To identify the statistically significant interaction energy changes caused by the N14K mutation, we used a nonparametric Kolmogorov–Smirnov (KS) test, which determines the significance of differences in the distributions of energies between two samples and does not assume a normal distribution of energies. To perform the KS test, we first pulled the data from all systems, all replicas, and all six subunits together and filtered out the interaction pairs between neighboring subunits that had a standard deviation <0.1 kcal/mol. Such a narrow distribution of interactions essentially means there was no change between WT and mutant. The 36 pairs with standard deviation >0.1 kcal/mol were subjected to the KS test between WT and N14K mutant using data from all replicas and all subunits. Of these, the energy distributions of the 29 pairs with the most significant P value (<10⁻¹⁶) from the KS test are shown in Fig. 7. The corresponding data, ranked by the absolute value of the mean energy difference, is provided in Table S1. Of the five pairs with the largest mean energy difference, K15–E101, K15–97, and V13–H100 (all between NT and TM2/CL regions) lost interaction energy in the N14K mutant, whereas N14K–D2 and K15–D2 (within the NT) gained interaction energy in N14K relative to WT. Of the 29 most significant changes, the K15–E101 interaction had the largest magnitude of disruption because of the N14K mutation.

To validate the changes in interactions indicated by the GB energies, we performed double mutant cycle analysis using free energy (ΔG) from the transition between the open and the closed states at different voltages. To extract the full ΔΔG from double mutant cycle analysis, mutations must be introduced that fully disrupt the targeted pair interactions. Typically, alanine residues are introduced to achieve this. We tried the N14A mutation to assess interactions with H100, but it gave nonfunctional channels. We therefore used the N14K mutant for these analyses. H100 was replaced by an Ala, so the double mutants were N14K/H100A (Fig. 8). To assess the coupling energy of the K15–E101 pair interaction, we found that none of the residues (Ala, Asn, Met, or Ser) introduced at the K15 position led to functional channels. However, because we had strong indications that the N14K mutation disrupts the K15–E101 interaction (Fig. 7), we used N14K (as a substitute for a K15 mutation), E101A, and N14K/E101A mutants to indirectly evaluate the coupling energy of the K15–E100 pair interaction. We did so knowing that if the MD was incorrect in showing that N14K disrupts the K15–E101 interaction, the mutant cycle analysis would not produce a significant ΔΔG. It is also important to note that our simulations indicate that N14 or N14K and E101 do not physically interact; thus double mutant cycle analysis for N14–E101, while reflecting perturbation of the K15–E101 interaction, might not yield the full value of the perturbation.

Fig. 8 shows the voltage dependence for WT and each of the single and double mutants. By fitting these values to a Boltzmann equation, we were able to calculate the V_1/2 for each of the channels, and a ΔG associated with the voltage dependence was calculated (Table 1 and Fig. 8). From each of these ΔGs, a ΔΔG was obtained (Table 2 and Fig. 7). A ΔΔG value of 0 means that the interaction one is studying has no energetic contribution to the transition between the two assessed conformational states, in this case from open to closed (not that a physical interaction does not exist). A ΔΔG of >1 indicates a significant cooperativity between the two residues in the transition between the open and closed state. In Table 1, the ΔG values calculated using the V_1/2 are shown for each mutant, and the ΔΔG values for each pair of residues is shown in Table 2. Consistent with the simulations and GB energies between N14 and H100, the pair interaction displays a weak coupling energy of −0.59 kcal. Conversely, a highly significant ΔΔG value of −1.5 kcal was obtained for N14–E101. This indicates a strong contribution of the interaction between E101 and the N terminus, likely via K15, in the open-closed reaction of Cx26.

Analysis of correlated atomic motions in proteins offers insight into protein function. Here we use correlation analysis to detect long-range or allosteric change in atomic fluctuations introduced by the N14K mutation. To ensure that changes in the correlation are statistically significant, we applied the KS test (described above) to compare the distributions of pairwise residue correlation coefficients between WT and the N14K mutant using data from all six subunits. Fig. 9 shows the 25 pairs in which the distributions of correlation coefficients between WT and the N14K mutant differed with P < 0.005 (see the legend for detailed explanation). The corresponding raw data are provided in Table S2. It is clear from the distribution plot that all 25 pairs lose correlations in N14K. The top five pairs that showed the largest change in mean difference in correlation coefficient are between residues 19/20 at the NT/TM1 transition and residues 91/94/98 in TM2 of the neighboring subunit. Distinct from the pairwise interaction analysis shown in Fig. 7, the correlated α carbon backbone movement can be caused by direct interaction, concerted movement propagated along the protein backbone, or both, thus providing a broad indicator of the change in protein dynamics. These changes indicate that the rearrangements caused by the N14K mutation are significant enough to be observed even in these submicrosecond simulations.

Human disease-causing connexin mutations, particularly those that have altered function rather than produce functional nulls, are useful tools for understanding gating properties of connexin channels (Lopez et al., 2013, 2014, 2016; Sanchez et al., 2016; García et al., 2018). The recent crystal structures of Cx26 GJCs (Maeda et al., 2009) serve as guides to explore how these mutations affect interactions and domains that control opening and closing of the pore. In the present work, we studied the N14K mutation, located in the cytoplasmic portion of the NT domain. This mutation significantly decreases Ca²⁺ and voltage sensitivity compared with WT Cx26 hemichannels (Sanchez et al., 2016). For example, at extracellular Ca²⁺ concentrations or membrane potentials at which WT hemichannels are normally closed, N14K mutant hemichannels showed significantly greater Po. Single-channel analysis confirms that N14K mutation stabilizes the open channel conformation without change in unitary conductance.

To understand the structural basis of how N14K favors the open state of hemichannels, we performed MD simulations to assess how residue interactions are affected by the mutation. Analysis of pairwise residue decompositions of GB energies in WT and N14K hemichannel simulations revealed that the most statistically significant intersubunit interaction change caused by the N14K mutation is the loss of interaction of residue K15 with E101 of the adjacent subunit. Double mutant cycle analysis agrees with the role of this interaction in favoring the closed state and the ability of the N14K mutation to disrupt it. Correlation analysis of the molecular dynamic simulations indicates that the second half of the NT and the TM2/CL region at the adjacent subunit are functionally coupled, and that the coupling is significantly disrupted in hemichannels containing the N14K mutation. The simplest interpretation for our experimental and computational data is that disruption of pairwise interactions between the NT and the TM2/CL by the N14K mutation enhances occupancy of the open channel conformation. This points to the NT-TM2/CL interaction in WT channels as stabilizing a closed state of the channel.

We observed that the N14K mutation reduced the apparent affinity for Ca²⁺ at least twofold relative to WT hemichannels. A more severe and complex effect is observed on Ca²⁺-dependent deactivation kinetics; although WT hemichannels display monotonic deactivation time constants, the N14K mutants show slower and multiphasic kinetics. The expected source of multiphasic deactivations was contamination by endogenous currents, in particular, Ca²⁺-activated Cl⁻ currents. In studying other gain-of-function mutants, we have eliminated such currents by intracellular injection of 120 µM BAPTA (Lopez et al., 2016; Valdez Capuccino and Contreras, 2016). However, this did not affect the multiphasic nature of the tail currents in N14K (Fig. S1). Substituting Ba²⁺ for Ca²⁺ resulted in slower deactivation kinetics than in Ca²⁺, but they became monoexponential. The origin of the additional kinetic components of deactivation in Ca²⁺ is unclear, but given the resistance to BAPTA, they appear to be intrinsic to the interaction of Ca²⁺ specifically with the N14K hemichannel rather than because of activation of another channel. Using Ba²⁺ as a surrogate for Ca²⁺, we were able to quantitate the effect of N14K on tail current deactivation rates to show that it disfavored transitions to the closed state. The differences between Ca²⁺ and Ba²⁺ on deactivation kinetic behavior could be explained by effects of the mutation on the open-to-closed transition energy landscape unique to Ca²⁺.

At the structural level, there is compelling evidence that the binding site for extracellular divalent cations in connexin hemichannels is located at the external portion of the pore, specifically at the TM1/E1 region (Pfahnl and Dahl, 1999; Bennett et al., 2016; Lopez et al., 2016). The fact that N14K, a cytosolic mutation at the NT, drastically alters channel closure by divalent cations reveals a functionally important allosteric coupling between the extracellular and cytoplasmic sides of the pore, or more specifically between the TM1/E1 and NT. Structural and functional evidence suggests that the NT is part of, or strongly coupled to, a gate in connexin channels (Oshima et al., 2007, 2011; Ek Vitorín et al., 2016; García et al., 2018). This could provide a mechanism by which a mutation in the NT affects hemichannel deactivation in extracellular divalent cations.

The first crystal structure of Cx26 GJCs showed that the NT is folded into the intracellular pore entrance and forms hydrophobic interactions with the first segment of the TM1 domain (Maeda et al., 2009). Unfortunately, a recent high-resolution crystal structure for Cx26 GJCs in the presence of Ca²⁺ could not model the NT domain because the electron density in those regions was not sufficiently well defined (Bennett et al., 2016). Although no strong computational correlation is seen between residue 14 and the Ca²⁺ binding region in WT hemichannels (moving along the same vector simultaneously), our experimental data support the notion that the N14K mutation disrupts conformational changes propagated from Ca²⁺ binding residues to the cytoplasmic domains of the protein, specifically to the interaction between the NT and the TM2/CL region. Further structural and functional studies are necessary to determine how Ca²⁺ bound to the extracellular side of the pore is allosterically coupled with conformational changes at the NT and a gate.

To investigate how the N14K mutation affects gating, we also compared the conductance–voltage relationship in N14K mutant with respect to WT hemichannels. N14K mutant hemichannels display an enhanced opening at the large negative voltages at which WT channels are closed. There is also a decrease in the slope of the voltage dependence, which may suggest a reduction on the intrinsic charges that sense the membrane potential field and/or a less efficient transduction of energy from the voltage sensor to the gating mechanism.

Interestingly, Sanchez et al. (2016) suggest that N14K forms a π-cation interaction with H100 of the adjacent subunit (located at the TM2/CL region) that is critical to stabilize the open hemichannel conformation. However, in our GB energy comparison, N14K causes only a weak decrease of GB energy of positions 14 and 100. By contrast, the K15–E101 energy has a strong and consistent decrease in the N14K mutant across three intersubunit pairs. Thus, it is unlikely that H100 plays a critical role by interacting with N14K. This was supported by the double mutant N14K-H100A, which is unable to form a π-cation interaction. The deactivation kinetics, in Ba²⁺, were essentially the same as for the N14K mutant alone. This suggested that the N14K residue does not form an interaction with the aromatic π-motif of H100 to stabilize the open channel conformation. Further, single-channel studies of the double mutant N14K/H100A hemichannels showed that it has similar open probability to N14K hemichannels in the absence of extracellular Ca²⁺. Thus, a combination of multiple computational and experimental approaches rules out the possibility that a π-cation interaction between these residues stabilizes the open conformation of the channel.

To confirm the functional role of pairwise interactions between the NT and the CL/TM2 region, we performed double mutant cycle analysis. By using the V_1/2 activation for WT and mutants hemichannels, we were able to assess the energetic coupling of pair residues N14–H100 and N14(K15)–E101. The ΔΔG was, then, associated with transitions between open and closed states in response to voltage changes. We designated a ΔΔG of >1 kcal as strong coupling between pair residues; however, many reports use a cutoff for pair interactions of ΔΔG of >0.5 kcal. ΔΔG is an indirect and approximate measure of the free energy of physical interaction between residues; thus the lower the ΔΔG the weaker the interaction between pair residues. Consistent with this idea, ΔΔG values for the N14–H100 pairs interactions (−0.59 kcal) indicate a weak coupling between these residues, consistent with our GB interaction analysis. We could not test directly the K15–E101 coupling interaction because of the inability to find residues that could substitute for K15 and render functional channels. Thus, we tested the coupling energy between residues N14 and E101. Even though our simulation indicates that N14 and E101 are not physically interacting, our rationale to test this pair was that N14K greatly disrupts the electrostatic interaction between K15–E101. Strikingly, the coupling energy obtained when testing the N14–E101 pair interaction is ∼1.29 kcal, which is the strongest coupling among the tested pair residues. Because N14 and E101 are not in direct contact, we interpreted this to indicate that the effects of mutation of these residues reflects disruption of the electrostatic interaction between K15 and E101 (Fig. 7). We also note that because of the potential “indirectness” of the effect of N14K on interactions of K15, the resulting calculated ΔΔG is probably a lower bound of the real coupling energy between K15 and E101.

Analysis of pairwise residue decompositions of GB energies in WT and N14K hemichannel simulations also revealed statistically significant interactions of N14K and K15 with residue D2 in the adjacent subunit. Unfortunately, mutation of D2 to other residues did not yield functional hemichannels. This precluded further analysis to establish whether D2 interactions with N14K and K15 were important to stabilize the open conformation. Because we tried several other residues at different regions, it is also possible that the N14K mutation does not interact with any residue and that just the domain/segment disruption of the NT and TM2/CL interaction is sufficient to stabilize the hemichannel in the open conformation. This idea is further supported by single mutations at the TM2/CL, including H100A and E101A, that themselves can increase open channel probability at negative potentials compared with WT Cx26 hemichannels (Fig. S3). This leads to the notion that interactions between NT and TM2/CL facilitate or stabilize hemichannel closure.

To further understand the extent of the disruptions introduced by the N14K mutation in the Cx26 structure, correlation analysis indicated that the second half of the NT and the TM2/CL region interface of the adjacent subunit are functionally coupled. Statistical analyses reveal that at least 25 pairwise interactions were significantly uncoupled by the N14K mutation (Fig. 9). It is possible, and indeed probable, that we did not detect all the interactions altered by N14K because of the limited duration of the simulations (submicroseconds). This is expected, because gating transitions for N14K mutant hemichannels occur on the scale of seconds to minutes.

Based on our molecular modeling and experimental results, we propose that the interaction between the NT and the TM2/CL is a major contributor to stabilization of the closed state, and that this interaction is disrupted by the N14K mutation, enhancing occupancy of the open channel conformation. A recent study indicates that another human mutation at the NT, G12R, produces gain of hemichannel function. The data are interpreted to suggest that the G12R mutation prevents the NT from functioning as a blocking gating particle (García et al., 2018). The authors suggest that G12R forms an interaction with residue R99 in the adjacent subunit that partially pulls the NT domain away from the pore. In contrast with our observations for the N14K mutant, G12R does not affect the apparent affinity for Ca²⁺ or the overall voltage dependence at physiological Ca²⁺ concentrations. A possible explanation is that the G12R mutation only partially displaces the NT–TM2/CL interactions, whereas these are fully disrupted in N14K hemichannels. Importantly, computational and experimental work indicates that both mutations, N14K and G12R, affect NT-TM2/CL interactions, supporting the critical role of this region in gating.

To date, >100 human mutations in hCx26 have been shown to cause sensorineural deafness and dermatological disorders. Many of these mutations result in defective channel biogenesis or trafficking. However, a significant number produce connexin channels that are functional at the plasma membrane and cause pathology. These mutants likely induce altered gating and/or permeability. Indeed, the lack of associated skin disorders in cases of hCx26-null mutations show that the epidermal disease is not produced by the absence of Cx26 function, but rather gain-of-function changes (Martínez et al., 2009; Xu and Nicholson, 2013). We found that at cellular physiological potentials, when the WT hemichannels are generally closed, a significant fraction of the N14K hemichannels are still open. This could explain the pathology caused by N14K mutation, because severe skin lesions are thought to be caused by the increase of connexin hemichannel activity that disrupts the normal homeostasis necessary for skin proliferation and differentiation (Djalilian et al., 2006). Our work suggests that the molecular basis of the pathology is the disruption of pairwise interactions between the NT and TM2/CL regions, which leads to favoring of the open conformation of the Cx26 hemichannel.

We thank Yu Liu for technical support.

The research reported in this publication was supported by the National Institutes of Health/National Institute of General Medical Sciences (grants RO1-GM099490 to J.E. Contreras and RO1-GM101950 to A.L. Harris and J.E. Contreras); Fondo Nacional de Desarrollo Científico y Tecnológico grant 3150634 and Programa de Atracción e Inserción de Capital Humano Avanzado a la Academia PAI 79170081 to I.E. Garcia. The computational work was supported by National Science Foundation XSEDE research allocation MCB160119. J.E. Contreras is supported by the Health Resources and Services Administration through grant D34HP26020 to New Jersey Medical School Hispanic Center of Excellence.

The authors declare no competing financial interests.

Author contributions: J.M. Valdez Capuccino and J.E. Contreras designed research; Y. Luo designed the MD simulations; J.M Valdez Capuccino, I.E. García, and J.E. Contreras performed research; P. Chatterjee, W.M. Botello-Smith, H. Zhang and Y. Luo performed the MD simulations; A.L. Harris provided guidance for molecular dynamics design and interpretations; J.M. Valdez Capuccino, I.E. García, and J.E. Contreras analyzed experimental data; and J.M. Valdez Capuccino, I.E. García, A.L. Harris, Y. Luo, and J.E. Contreras wrote the paper.

José D. Faraldo-Gómez served as editor.