Molecular dynamics simulation of GJC2 mutants reveal pathogenic mechanisms of PMLD1 and SPG44

GJC2 is the gene-encoding connexin 47 (Cx47), a gap junction (GJ) protein expressed in oligodendrocytes as well as venous (Munger et al., 2016) and lymphatic endothelial cells (Ferrell et al., 2010). In the central nervous system, Cx47 provides for coupling between oligodendrocytes; Cx47 also participates in oligodendrocyte–astrocyte coupling by forming heterotypic channels with astrocytic Cx43 (Maglione et al., 2010; Wasseff and Scherer, 2011). Autosomal recessive pathogenic variants in GJC2 lead to at least three distinct disorders in the central nervous system. These include Pelizaeus–Merzbacher-like disease 1 (PMLD1; HLD2) (Henneke et al., 2008; Uhlenberg et al., 2004), Spastic paraplegia 44 (SPG44) (Orthmann-Murphy et al., 2009), and a minimal leukodystrophy (Abrams et al., 2014). Cx47, like almost every other connexin studied, forms intercellular accumulations (plaques), which are organized aggregates of cell–cell channels. We have previously shown that several pathogenic variants associated with PMLD1 lead to abnormal trafficking of Cx47, with ER-localized immunoreactivity for Cx47, and loss of plaque formation and coupling (Flores-Obando et al., 2022; Orthmann-Murphy et al., 2007a). On the other hand, the SPG44-associated mutant p.I33M appears to traffic normally, forming WT-like plaques; however, cell–cell channels formed by this mutant are predicted to be closed under normal physiologic conditions (Orthmann-Murphy et al., 2009). The p.R98L mutation, associated with minimal leukodystrophy, also shows Cx47 staining in a pattern indistinguishable from that produced by the WT but induces reduced levels of GJ coupling between apposed cells expressing the mutant (Abrams et al., 2014). In this communication, we describe the clinical and cellular phenotypes seen in two families with novel GJC2 mutations. In the first case, the proband and two siblings show a complicated spastic paraparesis phenotype caused by homozygous pathogenic variants at codon 40 leading to a p.G40S mutation. In the second patient, compound heterozygous variants, p.R244P and p.G40S, lead to a much more severe, classic PMLD1 phenotype.

To gain further insight into the mechanisms that distinguish p.G40S-induced pathology from those caused by the p.R244P mutant, we conducted MD simulations of WT Cx47 and the two mutants. Our trajectory analyses demonstrate that p.G40S introduces a secondary structure disturbance in the pore-lining α-helix, TM1. This disrupts the interactions of E41 with the N-terminal amide and those of multiple residue pairs in the electrostatic network, which are strongly conserved across the connexin family. In contrast, p.R244P nullifies electrostatic interactions between R244 and lipid headgroups in the plasma membrane and inhibits β-sheet formation for multiple residues in Cx47’s extracellular domains. These two disease-causing mutations are spatially distinct, suggesting that the impeded functions of Cx47 are unrelated; however, MD simulations demonstrate participation of R244 in the electrostatic network, suggesting that there may be overlap in their underlying biophysical mechanisms of disease.

Cx47 WT, Cx43 WT, and Cx47 mutants were cloned into pIRES2-EGFP (Clontech) or pIRES2-DsRed as described (Orthmann-Murphy et al., 2007a, 2007b, 2009). For recording, confluent Neuro2a cells obtained directly from ATCC were transiently transfected with 500 ng of WT or mutant Cx47 constructs using LTX and Plus Reagent (Life Technologies), as per the manufacturer’s instructions and as previously described (Abrams et al., 2013). For cells expressing both mutants, 250 ng of each construct was mixed prior to mixing with the transfection reagents. Heterotypic cell pairing used the method of Orthmann et al. (Orthmann-Murphy et al., 2007b), mixing cells expressing a pIRES2-EGFP construct with cells expressing a pIRES2-DsRed construct in a 1:1 ratio. Coupling was assessed by dual whole-cell patch clamping of cell pairs 6–48 h after replating as previously described (Abrams et al., 2014). The following recording solutions (in mM) were used: pipette solution, 145 CsCl, 5 EGTA, 1.4 CaCl₂, and 5.0 HEPES, pH 7.2; bath solution, 150 NaCl, 4 KCl, 1 MgCl₂, 2 CaCl₂, 5 dextrose, 2 pyruvate, and 10 HEPES, pH 7.4. Heterotypic pairings between cells are shown as “connexin expressed in cell 2/connexin expressed in cell 1,” where cell 1 is the cell being pulsed. Junctional conductances (g_j) were determined by measuring instantaneous junctional current (I_j) responses to ±40 mV junctional voltage (V_j) pulses. Additional ±100 mV pulses were evaluated to ensure that activating conductances resulting from shifts in the voltage dependance of gating would be detected. Cytoplasmic bridges were excluded by demonstrating sensitivity of junctional conductance to application of octanol-containing bath solution.

Neuro2a cells were transiently transfected as described above. The cultures were incubated for 24 h before being processed for immunofluorescence labeling, and immune-labeled of Cx47 using rabbit anti-Cx47 (gift of Steven Scherer, University of Pennsylvania, Philadelphia, PA, USA, Cat# Rb/Cx47 935, RRID: AB_2819034, 1:1,500) in conjunction with Alexa 594 (1:400) chicken anti-rabbit secondary and counterstained with 4′,6-diamidino-2-phenylindole as previously described (Orthmann-Murphy et al., 2007b). For presentation, slides were z-stack imaged using an Olympus Bx61 microscope with a 100× oil immersion lens (NA: 1.30) with identical exposures. Z-stacks were merged via maximal intensity projection, producing individual images for display, which were montaged in Fiji (Schindelin et al., 2012). To quantify plaques, 10 distinct regions were imaged per condition under a 40× lens (NA: 0.75). Both plaques and pairs of apposing, anti-Cx47–labeled cells were counted to determine plaque formation frequency per cell pair.

Statistical testing was done in GraphPad Prism 10.0 (Dotmatics) using parametric or nonparametric inferential statistics as noted in the text or figure legends.

A Protein Data Bank (PDB) structure of Cx47 WT was constructed via homology modeling in SWISS-MODEL (Waterhouse et al., 2018) using the extensively equilibrated Cx26 hemichannel as the template (Kwon et al., 2011). Sequence similarity of Cx47 to Cx26, excluding the disordered cytoplasmic domains, is 66.1%. CHARMM-GUI’s Membrane Builder was used to assemble the final systems (Lee et al., 2016). Amino acid substitutions were introduced via the PDB manipulation module. Three disulfide bridges per protomer, identified in the Cx26 crystal structure and predicted to be widely conserved across the connexin family (Bai et al., 2018; Flores et al., 2020), were explicitly included between residues 53 and 253, 60 and 247, and 64 and 242. Each hemichannel was then embedded in a POPC lipid membrane bilayer with the channel pore aligned along the z axis. This protein–membrane system was solvated in 0.15 M KCl with TIP3 waters. The final atom counts for each tetragonal, solvated, protein–membrane system are Cx47 WT, 341,841 atoms; p.G40S, 342,105 atoms; and p.R244P, 341,752 atoms.

MD simulation was carried out with the pmemd.cuda 20.0 module of Amber20 (Case et al., 2020) employing the CHARMM36m force field (Huang et al., 2017) with periodic boundary conditions. Each system was initially minimized for 2,500 cycles of steepest descent, followed by 2,500 cycles of conjugate gradient descent, then equilibrated for 1.875 ns using the default initial equilibration scheme of CHARMM-GUI’s Membrane Builder, consisting of six steps of successive constraint relaxations. Finally, three independent NPAT ensembles per condition were simulated for a minimum of 100 ns each, extending one simulation per condition to 500 ns to confirm system stability. Constant pressure was maintained at 1.0 atm with a Monte Carlo barostat, constant temperature was maintained at 303.15 K via Langevin dynamics, hydrogen bonds were constrained with SHAKE, and long-range electrostatics were calculated using the Particle Mesh Ewald method. System equilibration of the constraint-free dynamics was gauged via RMSD of the whole protein, the channel-forming domains, and the minimum channel pore width (See Figs. S1 and S2).

For each condition, three constraint-free dynamics replicates were simulated at a 2-fs time step, retaining trajectory frames at every 10 ps for a minimum of 100 ns, with one replicate for each condition extended to 500 ns. In addition to these three closed system replicates, we simulated eight additional constraint-free systems per condition, initialized from the same equilibrated structures with an additional external electric field applied along the channel pore axis using the Amber mdp parameter efz to validate the conductive nature of our models under physiological transmembrane potentials. The applied external electric fields spanned ±50 and ±100 mV, with two replicates simulated for at least 100 ns per voltage.

All trajectory statistical analyses were done sampling every 10 frames (0.1 ns) and without the first 10 ns of simulation to mitigate artifacts from early equilibration. CPPTRAJ 6.18.1 was used for trajectory analysis to obtain atomic positions, dynamics, and interactions of interest; secondary structure dynamics as defined by the DSSP algorithm (Roe and Cheatham, 2013). VMD 1.9.4a51 Salt Bridges Plugin 1.1 was used to quantify salt bridges across the trajectories (Humphrey et al., 1996). HOLE 2.2.005 program was used to quantify structural properties of the connexin pore (Smart et al., 1996). Generally, a cutoff of 4 Å was taken to define an interaction between residues unless noted otherwise. Custom Python scripts were written to count K/Cl ion passaging through the channel pore in simulations conducted with an applied external electric field. Data after processing and plotting were performed in Python 3.11.4, and molecular visualizations were produced using VMD 1.9.4a57.

The Office for the Protection of Research Subjects of the University of Illinois Chicago College of Medicine reviewed this study and assigned a determination of Not Human Research. Informed consent was obtained in each country of patient encounter. Per Dr. El-Hattab, patients from family 1 provided consent at Tawam Hospital in accordance to the Helsinki Declaration. The patient from family 2 provided consent through the Myelin Disorders Bioregistry Project (IRB 14–011236; MDBP) at the Children’s Hospital of Philadelphia. Dr. Orthmann-Murphy at the University of Pennsylvania has a reliance agreement with the MDBP.

The proband is a 41-year-old man and the product of a consanguineous marriage (first cousins). His relevant medical history is as follows: He began to experience difficulty walking at age 3 and was first evaluated at age 10 when his 17-year-old sister was also evaluated for gait abnormalities. By age 18, on joining the military, he was noted to have gait abnormalities and some difficulty with using his hands and was assigned to clerical duties. At age 29 he began to require crutches and by age 41 was wheelchair bound. On exam at age 41, he was noted to have breakdown of smooth pursuit eye movements but no nystagmus, normal strength in the upper extremities, 4/5 weakness in iliopsoas muscles, bilateral dysmetria, hyperactive deep tendon reflexes with ankle clonus, and Babinski responses bilaterally. His parents were examined; his father was healthy, and his mother walked with a cane. He had two children who were reportedly healthy. He also has a 42-year-old brother with a milder phenotype (gait disorder beginning at age 30) who works in a cognitively demanding clerical role. The neurological exam of the brother demonstrated nystagmus on lateral gaze, hyperreflexia, mild iliopsoas weakness, and difficulty with tandem gait. The 48-year-old sister of the proband had a progressive gait disorder beginning in childhood, requiring use of crutches by age 18 due to progressive weakness, as well as spasticity and urinary urgency. Her exam demonstrated evidence of scanning speech, bilateral dysmetria, and bilateral lower extremity spasticity. MRI imaging (Fig. 1) showed thinning of the corpus callosum and diffuse moderately abnormal white matter signal. All three siblings appeared cognitively normal; however, formal neurocognitive testing was not performed. Genetic testing of the proband (Centogene Neuro panel, including 1,038 genes) was conducted and showed homozygosity for GJC2 pathogenic variant c.118G>A (p.G40S), and genetic testing confirmed the same homozygous pathogenic variant in the GJC2 gene, c.118G>A (p.G40S) in both siblings.

Limited past medical history is available, but the patient suffered from a very early age from visual impairment, seizures, developmental regression, ataxia, myoclonus, and dysphagia. The patient is nonambulatory. There was also a history of cleft palate repair. MRI at 6 mo of age (Fig. 2) showed profound hypomyelination which persisted on studies at 2, 3, and 5 years of age (not shown). Genetic testing showed compound heterozygous mutations in GJC2, c.118G>A (paternal), and c.731G>C (maternal), resulting in p.G40S and p.R244P, respectively.

Analyses conducted to evaluate the quality of the MD simulation replicates, as well as the statistics and visualizations describing all residue side-chain pairs, are included as supplementary material. Fig. S1 shows Cx47 WT and variant simulation RMSDs across time relative to the first frame of simulation. Fig. S2 shows minimum pore radius, determined by the HOLE2 program, across simulation time for three replicate simulations for Cx47 WT, Cx47 p.G40S, and Cx47 p.R244P. Fig. S3 shows a volumetric map-derived water density isosurface superimposed on our Cx47 WT structure. Fig. S4 shows the side-by-side alignment of our Cx47 homology model (left) with the AlphaFold2-predicted model (right). Fig. S5 shows the MD simulation calculated I-V curve fitted piecewise for positive and negative transmembrane voltages to partially account for possible channel rectification. Fig. S6 shows the side-chain center of mass distances between oppositely charged residue pairs that demonstrate salt bridge formation. Table S1 shows the proportion of WT Cx47 protomers that, at any point along the MD trajectories, demonstrate salt bridge interactions between the noted residue pairs.

To begin our investigations into the mechanisms by which expression of the p.G40S and p.R244P Cx47 mutants lead to disease, we expressed each in Neuro2a cells to examine their ability to form GJ plaques and intercellular communication channels. As shown in Fig. 3, p.G40S shows a pattern of staining very similar to that of WT with areas of intense staining at appositions between expressing cells. In contrast, the p.R244P mutant showed staining in a pattern suggestive of ER accumulation, similar to what we have shown with other PMLD mutants (Orthmann-Murphy et al., 2007a, 2009). To model the compound heterozygous state, we co-expressed Cx47 containing the p.R244P mutation with Cx47 expressing the p.G40S mutation. In this case, plaques are seen, albeit at a significant reduction compared with the cells expressing the p.G40S mutant alone. This reduction is roughly proportional to the amount of p.G40S DNA used to transfect each of these sets of cells (500 ng/well for p.G40S alone versus 250 ng/well for each construct in the p.G40S+p.R244P condition), suggesting that the p.R244P mutant does not interact with the p.G40S mutant in a way that reduces the ability of p.G40S to form junctional plaques.

We then used dual whole-cell patch clamping to evaluate the ability of these two mutants to form cell–cell channels. As shown in Fig. 4 A, cell pairs expressing Cx47 WT efficiently formed intercellular channels, leading to high levels of intercellular coupling; on the other hand, when expressed in apposed cells, neither the p.R244P nor the p.G40S mutant induces intercellular coupling statistically above the levels seen in the negative control cells. Furthermore, pairs of cells co-expressing both the p.R244P and the p.G40S mutant forms of Cx47 also failed to induce conductances statistically above those seen in negative controls. As noted above, Cx47 is expressed in oligodendrocytes, while astrocytes express Cx43. Thus, oligodendrocytes can couple with one another via homotypic Cx47 cell–cell channels, but they are also able to couple with astrocytes via heterotypic channel formation between Cx47 and Cx43. For this reason, we evaluated whether the p.R244P and the p.G40S mutants could induce functional coupling when paired with cells expressing WT Cx43.

As shown in Fig. 4 B, neither the p.R244P nor the p.G40S mutant induces intercellular coupling statistically above the levels seen in the negative control cells paired with Cx43. Furthermore, when paired with cells expressing Cx43, cells co-expressing both the p.R244P and the p.G40S mutants also failed to induce conductances statistically above those see in negative controls. The extremely low levels of activating currents seen in heterotypic pairings of mutants with Cx43 are qualitatively and quantitatively similar to those seen for endogenous connexin expression cells paired with Cx43. Fig. 4 C shows representative current traces for pairs used for evaluation of coupling in Fig. 4, A and B.

Connexins are the family of transmembrane proteins that serve as the building blocks for GJ channels. Subcellular location divides a connexin protein into nine distinct motifs: an α-helical N-terminal domain (NT-helix) that resides in the channel pore, a four-helix transmembrane bundle (TM1–4), two extracellular stretches (E1 and E2, which connect TM1 to TM2 and TM3 to TM4, respectively) that form antiparallel β-pleated sheets, a cytoplasmic loop connecting TM2 to TM3, and a cytoplasmic C terminus. These secondary and tertiary structural features are conserved across the entire protein family, evidenced by topology analyses (Milks et al., 1988; Willecke et al., 1991; Yeager and Gilula, 1992), by high sequence identity within these domains via multiple sequence alignments (Söhl and Willecke, 2004), and more recently, by the growing number of experimental structures which confirm this patterning. The majority of sequence variability among connexins is localized to the cytoplasmic loop and the C-terminal tail. Six connexin protomers aggregate annularly, forming a connexin hemichannel (Fig. 5 A); two hemichannels, one each from apposing cell membranes, join at their extracellular domains to form a dodecameric GJ channel.

The average equilibrated Cx47 hemichannel, shown in Fig. 5, A and B, has a minimum pore radius of 5.1 Å, determined by the N-terminal methionines, and reaches a maximum of ∼10 Å at the intra- and extracellular orifices. A volumetric map derived from bulk water average density confirms channel patency and hydration (Fig. S3). Comparison to the AlphaFold model (AF-Q5T442-F1; Fig. S4) demonstrates highly corresponding secondary and tertiary structural organization of the transmembrane and extracellular domain; however, the disordered cytoplasmic domains of the AlphaFold homology model display clearly nonphysical residue placements (stretches of disordered protein where the plasma membrane would necessarily be). Simulations with applied transmembrane voltages were carried out to assess the functional nature of the equilibrated Cx47 homology model: two replicates each at −100, −50, 50, and 100 mV (mirroring the range of transmembrane conditions of a patch-clamp experiment). An estimate of ion conductance derived from counting the net charge passage through the channel over time is 224 pS, within a reasonable, twofold range for MD simulations of an open state Cx47 hemichannel. Given a measured Cx47 cell–cell conductance of 55 pS, determined previously (Orthmann-Murphy et al., 2007b), one would predict a hemichannel conductance of ∼110 pS. However, it is not unexpected for connexin channel simulations to demonstrate variability in ion permeation counts when constrained to relatively short simulation times (Kwon et al., 2011). Furthermore, MD simulations employing the standard TIP3P water model have been shown to overestimate experimental diffusion coefficients by a factor of 2–3 (Aksimentiev and Schulten, 2005; Vega and Abascal, 2011), bringing our simulation conductance to a range of 75–112 pS, strongly agreeing with experimental results. At the very least, our simulations confirm the conductive capability of our Cx47 model, the alternative being an impermeable system. The transmembrane voltage simulations exhibit slight channel rectification. Piecewise fitting of positive and negative transmembrane voltages piecewise predicts conductances spanning 182–269 pS (Fig. S5). The inward rectification is qualitatively similar to that seen for the highly homologous Cx45 hemichannel (Valiunas, 2002) (Cx45 primary sequence is highly homologous to Cx47 in the NT, TM1, TM2, and E1 domains, the predicted contributors to the pore lining segments of the hemichannel).

Our Cx47 MD model conforms to the structural patterning that has been observed in experimental structures of homologous connexins. All experimentally determined connexin structures demonstrate three distinctive physical features in addition to those already mentioned: (1) a network of intra- and interprotomeric electrostatic interactions formed by the residues within the extracellular mouth of the hemichannel, (2) a nexus of hydrophobic side chains that binds the four transmembrane α-helices together, lying adjacent to the electrostatic network, and (3) three charged amino acids distributed among TM1, TM3, and TM4, at the level of the NT-helix (Fig. 5, C and D).

Using a distance cutoff of ≤4 Å between charged side chains as defining a salt bridge interaction, 14 residues from each Cx47 protomer contribute to the electrostatic network of the Cx47 hemichannel. In our WT simulations, the three most prominent interprotomeric salt bridges of this network are between D66–R257, R75–E260, and K50–D61, and the salient intraprotomeric salt bridges are between D46–K50, D66–R244, E235–R237, E47–K261, R56–D61, and R257–E260. Fig. 5 D depicts how these electrostatic interactions are organized along the Cx47 hemichannel membrane topology. Some of these salt bridges are dynamic, breaking and forming over time, so there are transient interactions of lesser frequency, which may serve the distinct functional roles of each residue. A full list of electrostatic interactions observed in at least one monomer across all simulations is presented in the Table S1; the dynamics of all salient interactions are presented in Fig. S6.

Hydrophobic cores, observed in experimentally determined connexin structures, are present here as well. One core, located within the electrostatic network, is centered around residue Y44 (homologous to W44 observed in the Cx26 crystal structure; PDB ID 2ZW3), along with L43 and M268. The second hydrophobic core is located near the geometric center of the four-helix bundle (W77, F76, and L265) and forms a nexus, around which hydrophobic side chains of neighboring residues from all four helices are organized. This stable center may serve as a structural anchor, resisting denaturation induced by gating conformational change.

At the lateral level of the NT-helix, adjacent to the W77 electrostatic cluster, TM1 is strongly electrostatically bound to TM3 via R32–E221. This salt bridge is observed in all experimentally determined connexin structures. Then, one step further toward cytoplasmic face, TM1, TM3, and TM4 contain K22, R217, and E282, respectively. These three residues’ side chains demonstrate relatively greater conformational flexibility across our simulations, and their electrostatic pairing pattern varies (for example, E282 shows relatively similar preferences for K22 or R217). Interestingly, in other experimental and theoretical connexin structures, pairings among the side chains of members of this electrostatic triad vary. The average equilibrated Cx26 (Kwon et al., 2011) and the homology modeled Cx32 structures (Abrams et al., 2013) (our unpublished data) demonstrate pairings of the homologous residues in TM1–TM3; the Cx46 cryo-EM GJ (Flores et al., 2020), a MD model of Cx46 (Abrams et al., 2018), and the presented homology modeled Cx47 structures demonstrate pairings of the homologous residues in TM3–TM4; and both of the putatively closed Cx43 structures show an absence of either interaction (Lee et al., 2023; Qi et al., 2023). Given their proximity to the NT, further study of the dynamics of this charged residue triad may provide insight into the mechanism of fast voltage gating via NT-helix pore-occlusion (Verselis et al., 1994).

To further examine the effects of the p.G40S and p.R244P mutations on Cx47 structure and gain insight into the potential effects on channel function, we employed in silico mutagenesis combined with MD simulations.

p.G40S exhibits marked instability around the voltage-sensing apparatus of Cx47. Analysis of electrostatic interactions reveal that E41 is the only negatively charged residue within at least 10 Å of the positively charged N terminus (Fig. 6 A). Simulations of p.G40S reveal that introduction of a serine induces a secondary structure shift such that the serine side chain forms a hydrogen bond with the backbone of I39. This restricted conformation disrupts the flexibility of the neighboring V38, G39, and E41 as shown in the DSSP visualizations (Fig. 6 B). By inhibiting E41 from forming its native open-state interaction with the positively charged N terminus, this mutation may disrupt the voltage-sensing capabilities of the assembled connexon. It is interesting to note that the α-helical content of the adjacent upstream residues, V34 through A37, is increased, suggesting that the flexibility of this portion of TM1 is likely reduced, which may also participate in reducing the desirable interaction between E41 and the N terminus. In our WT simulations, the average distance between the NT amine and E41 of a protomer is 5.53 ± 3.04 Å. In contrast, the same distance measured in the p.G40S simulations is 7.52 ± 3.04 Å. On the other hand, p.R244P simulations are more similar to the WT, with distance measured being 5.83 ± 3.70 Å (Fig. 6 C).

In addition, as discussed above, Y44 serves as a hydrophobic stabilization center. The residues in this neighborhood form the TM1/E1 boundary, which has been suggested to play a role in voltage sensing via “loop gating” (Verselis et al., 1994). Y44 and its neighboring residues I43, S45, and D46 lose significant α-helical content.

Further downstream are D46 and E47, members of the network of charged residues that are strongly conserved across the connexin family. Previous experimental and simulation studies implicate this domain’s central role in mediating calcium-induced hemichannel closure. Upon p.G40S mutation, the intraprotomeric E47–R75 interactions appear slightly weakened, whereas the interprotomeric R75–E260 interactions become strongly clamped (Fig. 6 C). These observations, in tandem with the weakening of the E41–NT interaction and alterations at residues 38 and 39 noted above, suggest that p.G40S assumes a closed conformation at equilibrium.

p.R244P alters local lipid dynamics and perturbs the native conformation of the extracellular loops. Two significant properties are altered upon p.R244P mutation. First, the positively charged side chain is lost. Second, proline introduces a strict condition on the backbone conformation at that position within the second extracellular loop. The extracellular loops are known to be strongly conserved across connexins and are thought to determine the coupling propensity between hemichannels. Disruption of essential intercellular side chain interactions would be expected to impede the formation of GJ plaques.

It is apparent through our simulations that R244 consistently forms stable electrostatic bonds with the negatively charged phosphate group of lipid heads as shown in Fig. 7 A. At most time points, R244 from two of the six protomers is closely associated with a polar lipid headgroup, leading to an average probability of association of ∼0.33 for each R244 residue. This probability of interaction at residue 244 is significantly higher for both WT and p.G40S than for the p.R244P mutant. The tendency for lipid molecules to dwell in the milieu of the residue at 244 suggests that the p.R244P mutation disrupts normal protein–lipid dynamics and subsequently has a negative impact on hemichannel formation and plaque aggregation.

p.R244P also appears to restrict secondary structure formation at the extracellular loops, which are known to take the form of β-sheets (Maeda et al., 2009). DSSP analysis of the WT (Fig. 7 B) reveals that residues 243, 244, and 247–250 are able to intermittently assume β-strand conformations; however, the introduction of proline at 244 prevents that from occurring.

The electrostatic network is also perturbed by this mutation. Obviously, D66–R244 is lengthened, as the arginine side chain is far longer than that of proline. More interestingly, the interprotomeric R75–E260 interaction appears far weaker, in the very opposite direction as induced by p.G40S. In fact, p.R244P appears to weaken salt bridges without compensatory effect, unlike p.G40S for which strengthening for at least one salt bridge is apparent.

Inspection of the original Cx26 crystal structure revealed that the predominant intercellular interactions between the extracellular loops of apposed hemichannels were N54–L56, Q57–Q57, N176–K168, N176–D179, and N176–T177 (Maeda et al., 2009). Comparing the homologous residues (as determined by sequence alignment) in Cx47 (See Table 1) there is no obvious difference in their equilibrium conformations. However, there is a noticeable distortion of the center β-strand, which contains H249 and D252. Since H249 forms three of the five homologous hydrogen bonds found in Cx26, it may be the case that the restricted conformation of p.R244P contributes to inhibition of proper interfacing between apposed hemichannels.

In this communication, we present data outlining the consequences of two pathogenic variants in GJC2 that lead to either mild (SPG) or more severe phenotype (PMLD). In the homozygous configuration, p.G40S is associated with a relatively mild later-onset complex SPG phenotype, similar to what we previously reported for a family with SPG who had a homozygous pathogenic variant in GJC2 (p.I33M) (Orthmann-Murphy et al., 2009). p.G40S appears three times in the Genome Aggregation Database (gnomAD v2.0), giving an estimated population prevalence of 1.2 × 10⁻⁵, a reasonable frequency for a rare allele causing a rare recessive disorder seen in a consanguineous situation. Like Cx47 WT, when expressed exogenously, the p.G40S mutant traffics to the cell surface and forms GJ plaques but does not form functional channels between apposed cells. We previously showed that in the SPG44 variant p.I33M, subcellular localization was similar to WT and channel function was altered in a way predicted to lead to complete loss of coupling between oligodendrocytes and near complete or complete loss of coupling between oligodendrocytes and astrocytes. Together, these assays support the idea that GJ channel formation without establishment of functional GJs leads to a milder neurologic phenotype than when GJ formation is impacted.

We also present the case of a patient with compound heterozygous pathogenic variants, a combination of the above variant (p.G40S) with p.R244P, a variant which causes a severe conatal form of PMLD. This variant is found only once in gnomAD 2.0, giving a population frequency of 4.0 × 10⁻⁶. When this mutant is expressed in Neuro2a cells, no GJ plaques are formed, and the staining pattern is suggestive of a failure to exit the ER, similar to what we have reported for several other PMLD1 mutants (Orthmann-Murphy et al., 2007a). In addition, no GJ coupling is detected in either the homotypic or heterotypic (with Cx43) configurations. When p.R244P is co-expressed in Neuro2a cells with p.G40S, GJ plaque formation still occurs, but no coupling is seen in any tested configuration (homotypic or heterotypic with Cx43). The simplest interpretation of this observation is that expression of p.R244P protomers do not exert a dominant negative effect on the ability of p.G40S to traffic to the cell surface.

However, the reason for the much more severe clinical phenotype in the patient with p.G40S and p.R244P compound heterozygous variants, compared with the patients homozygous for p.G40S alone, is not immediately evident. One possibility is that both p.G40S and p.I33M, two examples of GJC2-associated hereditary spastic paraparesis, retain some function of Cx47 that is dependent on trafficking to the cell surface, such as ZO-1 binding and localization of the transcription factor ZONAB to the cell membrane (Li et al., 2008). In this model, PMLD1-associated pathogenic variants (including p.R244P) lead to complete loss of Cx47 function, with lack of proper cell surface localization and failure to form functional GJs. Our data, however, would seem to argue against this idea, since they suggest that a compound heterozygote p.G40S+p.R244P patient would form plaques. However, we cannot rule out the possibility that the p.R244P mutant may have a dominant negative effect on p.G40S plaque formation in vivo or that a quantitative reduction in GJ plaques may be sufficient to cause the more severe phenotype. Another possibility is that the p.R244P mutant confers a toxic gain of function, which in combination with a loss of function allele (p.G40S) leads to a more severe phenotype than in patients homozygous for the loss of function mutant (p.G40S). While we have no evidence for p.R244P conferring a toxic gain of function, we (Flores-Obando et al., 2022) and others (Chen et al., 2017) have recently shown that several other PMLD1 mutants cause activation of the CHOP-mediated arm of the unfolded protein response. A third possibility is that there are other genetic loci that may modulate the severity of PMLD1, driving the variability of prognoses across patients. To date, at least 50 variants in Cx47 have been described in the literature; however, at this point only a limited number have been evaluated for trafficking abnormalities or ability to promote cell–cell coupling. Thus, further investigation may still be required to identify the unique determinants of an SPG44 versus a PMLD1 phenotype in the homozygous versus compound heterozygous state. This understanding is critically important in the design of potential gene therapy interventions for PMLD (Georgiou et al., 2017).

To better understand the functional implications of these pathogenic variants on Cx47 structure, we performed MD simulations to evaluate biomolecular systems at the atomic level by leveraging physical laws. This complements experimental approaches which cannot routinely probe such small spatial resolutions and timescales. We conducted MD experiments on Cx47 and two of its pathogenic variants to probe for the biophysical determinants, which might distinguish SPG44 and PMLD1 causing variants. The p.G40S and p.R244P mutants are predicted to have distinct, non-WT behaviors when modeled in silico. In particular, each mutant induces significant secondary structure disturbances in the spatial neighborhood around their residue positions. The p.G40S simulations predict a closed conformation due to deviations of the N-terminal helix from its native “open” WT position. In the case of p.R244P simulations, the extracellular β-sheet domains show a pronounced decrease in integrity, coincident with a loss of a potentially important protein–lipid interaction, suggesting that conformational distortions in Cx47’s extracellular domains nullify proper GJ docking and assembly.

Immunofluorescent labeling shows that GJ plaques composed of Cx47 p.G40S do form between adjacent cells. However, our electrophysiological experiments demonstrate that GJs are unable to conduct ions. Therefore, the connexin localization mechanism remains intact, and the dysfunction of p.G40S mutants is most likely to be rooted in the Cx47 gating mechanism. The N-terminal α-helix has long been thought to govern connexin channel closure, and recent cryo-EM structures of Cx43 reveal a stable conformation in which the six NT helices intertwine to form a wall at the intracellular orifice (Lee et al., 2023; Qi et al., 2023). In our simulations, we predict that the p.G40S mutant introduces destabilization in the NT helices in TM1. This finding is consistent with the suggestion that TM1 has a role as a core regulator of the gating mechanism for connexins (Bargiello et al., 2018).

Cells expressing Cx47 p.R244P do not form GJ plaques. Immunofluorescent labeling indicates that translated connexin gets trapped intracellularly, likely in the ER. In agreement with this, we also observe the absence of GJ-mediated ion conductance. Thus, experimentally, p.R244P mutants are translated but fail to properly localize. Possible explanations for this may include misfolding, mis-aggregation, loss of transport signaling, or overactivation of protein degradation pathways.

In our simulations, introduction of the p.R244P mutation into Cx47 predicts that the six residues 243, 244, and 247–250 no longer adopt the β-sheet conformation seen in the WT. This β-sheet extracellular domain is strongly conserved across connexins and is thought to define docking compatibility across the connexin family (Bai et al., 2018). Such a conformational disturbance likely has consequences on Cx47’s ability to form channels. Looking more closely, one of the affected residues, C247, is involved in a disulfide bond that is universally conserved for GJ-forming connexins (Foote et al., 1998). The conformational constraints imposed by a proline at position 244 may prevent proper alignment of the cysteine residues that pair together to stabilize the extracellular domain. This may have consequences not only on channel docking but more fundamentally on the folding of individual protomers.

A significant question we wish to pose is: “Is it possible to differentiate amino acid variants that are disease-causing from those that are benign?” Table 2 lists the connexin variants that are present in ClinVar and occur in disease-associated connexins, at positions homologous to 40 and 244 of Cx47 (Landrum et al., 2018).

None of the variants that occur at these two positions is classified as benign or likely benign, so we predict that the residues located at 40 and 244 or their homologous positions are critical to WT connexin functions. If mutant-induced conformational dynamics correlate strongly with disease phenotype, then MD simulations may serve as a reliable prognostic predictor of arbitrary connexin variants and may help adjudicate variants of uncertain significance in GJC2.

Our computational setup presumes that we may garner insight about the macroscopic consequences of amino acid substitutions from nanosecond scale simulations. However, it is important to note that our connexin homology models were produced using a template in a putative open-channel conformation. This has allowed us to observe consistent differences between their subsequent dynamics on a nanosecond timescale. Secondary structure perturbations are particularly clear, and potentially important protein–nonprotein interactions are highlighted. However, we are only capturing the stability (or instability) of the “open-state.”

Connexin proteins are multifunctional. Hemichannel assembly, protein localization, GJ docking/plaque formation, slow and fast gating, and possible allosteric interactions are phenomena that cover a wide range of timescales from ns to many minutes, and their dysfunctions lead to distinct disease phenotypes. All-atom computational studies of Cx47 provide access to the sub-microsecond timescale, which can tell us the local stability of a system. However, phenomena such as plaque formation and connexin channel gating, which operate at the µs or longer timescales require far more computational resources to fully sample. Future studies leveraging coarse-grained approaches and enhanced sampling techniques will be required for a holistic computational investigation of these phenomena.

All data underlying the study are summarized by the figures. Data sheets and raw data underlying the experimental and computational work are available upon reasonable request to the corresponding author.

Crina M. Nimigean served as editor.

We thank the patients and their representatives for granting permission to publish this information.

Author contributions: D. Gong: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, and writing—original draft, review, and editing. J.L. Orthmann-Murphy: conceptualization and writing—original draft, review, and editing. D. Kumar: formal analysis and methodology. G.D. Dungan: investigation and writing—review and editing. A.W. El-Hattab: conceptualization, data curation, investigation, resources, and writing—review and editing. N. Schiess: conceptualization and writing—original draft. Y.L. Luo: conceptualization, data curation, formal analysis, project administration, supervision, and writing—review and editing. M.M. Freidin: investigation, project administration, and writing—review and editing. C.K. Abrams: conceptualization, formal analysis, investigation, methodology, resources, supervision, and writing—review and editing.