Conserved patterns across ion channels correlate with variant pathogenicity and clinical phenotypes

Abstract

Clinically identified genetic variants in ion channels can be benign or cause disease by increasing or decreasing the protein function. As a consequence, therapeutic decision-making is challenging without molecular testing of each variant. Our biophysical knowledge of ion-channel structures and function is just emerging, and it is currently not well understood which amino acid residues cause disease when mutated. We sought to systematically identify biological properties associated with variant pathogenicity across all major voltage and ligand-gated ion-channel families. We collected and curated 3049 pathogenic variants from hundreds of neurodevelopmental and other disorders and 12 546 population variants for 30 ion channel or channel subunits for which a high-quality protein structure was available. Using a wide range of bioinformatics approaches, we computed 163 structural features and tested them for pathogenic variant enrichment. We developed a novel 3D spatial distance scoring approach that enables comparisons of pathogenic and population variant distribution across protein structures. We discovered and independently replicated that several pore residue properties and proximity to the pore axis were most significantly enriched for pathogenic variants compared to population variants. Using our 3D scoring approach, we showed that the strongest pathogenic variant enrichment was observed for pore-lining residues and alpha-helix residues within 5Å distance from the pore axis centre and not involved in gating. Within the subset of residues located at the pore, the hydrophobicity of the pore was the feature most strongly associated with variant pathogenicity. We also found an association between the identified properties and both clinical phenotypes and functional in vitro assays for voltage-gated sodium channels (SCN1A, SCN2A, SCN8A) and N-methyl-D-aspartate receptor (GRIN1, GRIN2A, GRIN2B) encoding genes. In an independent expert-curated dataset of 1422 neurodevelopmental disorder pathogenic patient variants and 679 electrophysiological experiments, we show that pore axis distance is associated with seizure age of onset and cognitive performance as well as differential gain versus loss-of-channel function. In summary, we identified biological properties associated with ion-channel malfunction and show that these are correlated with in vitro functional readouts and clinical phenotypes in patients with neurodevelopmental disorders. Our results suggest that clinical decision support algorithms that predict variant pathogenicity and function are feasible in the future.

Keywords: bioinformatics; epilepsy; genetics; ion channel; neurodevelopmental disorder.

Figure 1

Graphical summary of the study. (A) Dataset. First, a subset of ion channels was selected from IUPHAR and screened for available protein structures. Second, missense variants were assembled from three different databases: gnomAD, ClinVar and HGMD. Third, structural features were annotated to describe the location and secondary structure of all residues comprising the ion channels. (B) Analyses. Sets of residues described by a single feature or a combination of features that were enriched for pathogenic variants were identified across all ion channels. A subset of these features was identified to show enrichment for pathogenic variants across all considered channel families. On the basis of these shared features a two-dimensional reference framework was developed for describing the location of a residue with respect to the distance from the membrane and the pore. Using this framework, the highest variant burden was observed in pore residues that were located inside the membrane. Finally, within these pore residues, the impact of the biophysical pore environment on the variant burden was investigated. (C) Case studies. Our framework identified correlations between functional and clinical phenotypes and the localization of missense variants in the protein structure. GRIN = GRIN1, GRIN2A, GRIN2B genes; SCN = SCN1A, SCN2A, SCN8A genes; Mixed = electrophysiological readouts showed conflicting functional changes; PCA = Principal component analysis; WT = wild-type.

Figure 2

Identification of single features and feature combinations that are enriched for pathogenic variants. (A) Scatterplot of the odds ratio (x-axis) versus the P-value (y-axis) of the variant burden analysis of pathogenic versus population variants in each of the 163 different sets of residues described by one feature or a feature combination in the discovery cohort (Supplementary Fig. 1). Residue sets with significant enrichment of pathogenic variants after multiple-testing corrections (Fisher’s exact test, OR > 1, P < 0.05) are displayed above the dotted line. (B) Scatterplot of variants in the validation cohort. Only the subset of the 55 features and feature combinations that were found to be enriched for pathogenic variants in the discovery cohort is shown. Residue sets above the dotted line are enriched for pathogenic variants in both the discovery and validation cohort.

Figure 3

Residue distance from pore and membrane correlates with the pathogenic burden, harbouring the highest burden at membrane-spanning pore residues. (A) Graphical representation of the framework that describes the localization of each residue in ion-channel proteins on two dimensions (2D). The localization of each residue was depicted by the gene-wise normalized distance from the pore axis and the distance from the membrane centre, thereby allowing a comparison of the variant distribution across ion channels. (B) Table of the combined enrichment or depletion of pathogenic variants across all ion-channel proteins summarized over 200 different 2D regions. Each bin represents a distinct localization in the protein structure that is described by the normalized distance from the pore axis (x-axis) and the distance from the membrane centre (y-axis). Significant enrichments or depletions of pathogenic variants are indicated with a star. (C) The variant densities of pathogenic variants and population variants are displayed along with the normalized distance from the pore axis, separately for each channel class (n = 8). (D) Table of the enrichment or depletion of pathogenic variants in different sets of pore residues. Each set of pore residues was defined by their location (membrane, extracellular, intracellular) and secondary structure (x-axis), together with a description of which residues were considered as pore residues (y-axis). Significant values are indicated by a star. (E) Table as in D but based on a subset of channel genes that are constrained for variants in the general population (n = 13 genes).

Figure 4

Biophysical pore properties discriminate pathogenic and population variants and identify residues that are most likely to harbour pathogenic variants in the hydrophobic pore sections. [A(i)] Graphical representation of the different biophysical pore properties for each pore section along the pore axis. [A(ii)] Each pore residue is assigned to the biophysical properties of its closest pore section. (B) Contribution of the seven considered biophysical pore properties to the first (x-axis, PC1) and second (y-axis, PC2) dimension of the PCA that was performed on the pore properties together for all ion-channel proteins. (C) Scatterplot along the first two dimensions of the PCA (PC1 and PC2). Each dot represents a pore residue where a pathogenic variant (ClinVar/HGMD) or population variant (gnomAD) was observed. Pore residues where population and pathogenic variant were observed are not shown.

Figure 5

Distances to membrane and pore correlate with the clinical representation in sodium channels and NMDA receptors. (A) Nav1.2 protein structure (PDB-ID: 6j8e) encoded by SCN2A showing patient variants and population variants observed in SCN1A, SCN2A and SCN8A (SCN genes) that were aligned to SCN2A. Pathogenic SCN variants were curated from the SCN Portal, population variants from gnomAD. (B) Boxplot of the distribution of patients’ variants grouped by their clinical phenotypes along with the normalized distance from the pore centre. Boxes represent patient variants that were associated with an early seizure onset (seizure onset < median onset) or a late seizure onset (seizure onset > median onset) or patient variants associated with intellectual disability (ID) or developmental delay (DD). (C) Boxplots showing the same groups as in C, but along the distance to the membrane centre. (D) Heterotetrameric protein complex consisting of two Glu1N and two Glu2NA subregions (PDB-ID: 6ira) encoded by GRIN1 and GRIN2A, respectively. Patient and population variants observed in GRIN1, GRIN2A and GRIN2B (GRIN genes) are visualized on the structure. Patient variants were recruited from our internal variant database, population variants from gnomAD. Variants in GRIN2B were aligned to GRIN2A and were visualized on the GluN2A subregions. (E and F) Boxplots as in B and C of the clinical phenotypes observed in the GRIN patient cohort.

Figure 6

Distances to membrane and pore correlate with the molecular effect in voltage-gated sodium channels and NMDA receptors. (A) Scatterplot of variants associated with a molecular effect plotted with respect to the normalized distance to the pore axis (x-axis) and the distance from the membrane centre (y-axis). All readouts were assigned to the complementary SCN2A protein position based on a multiple sequence alignment. Variant with a mixed effect had contrary effects in different electrophysiological measurements. (B) Box plot of the molecular effects scattered across the normalized distance from the pore axis in SCN genes. (C) Box plot of the molecular effects scattered across the distance from the membrane centre in SCN genes. (D) Scatterplot as in A visualizing variants associated with a molecular effect in GRIN genes. (E and F) Box plots as in B and C show the variants classified with a molecular effect along with the normalized distance from the pore axis and membrane centre in GRIN genes.