Participants
This study describes 47 affected individuals from 29 families. We describe 26 newly identified individuals from 22 families (families F1-17, F23, F24, F26, F27, and F29) accrued through professional communication as well as through GeneMatcher [19]. They were recruited from hospitals or clinics from the USA and other countries including the UK, Germany, Italy, France, The Netherlands, Canada, and Singapore. Clinical, phenotypic, and KCNK9 genetic variant data were analyzed for each participant. All families providing new or updated data provided informed consent (F1-F19, F23, F24, F26-F27, and F29). Data for families F20-F22, F25 and F28 was abstracted from prior publications only. This study was approved by the Mayo Clinic Institutional Review Board, and local Ethics Committees. Family 9 was recruited to the Genomics England 100,000 Genomes Project. Previously published individuals with KCNK9 variants are also described with summarized published data and/or updated unpublished clinical data [5,6,7,8]. Persons P18.1, P19.1, P20.1, and P21.1 in this study refer to the published Patients 1-4 from Graham et al. [8]. Family F22 refers to the Arab-Israeli kindred described by both Barel et al. and Graham et al. [5, 8]. Person P25.1 refers to the individual reported by Šedivá et al. [6], Family F27 was previously reported without detailed clinical information [20], and person P28.1 refers to BK-227-03 described by Guo et al. [7].
Variant identification
KCNK9 variants were identified through a variety of methods including gene panel testing or exome or genome sequencing through commercial clinical laboratories including Ambry Genetics (Aliso Viejo, CA), Baylor Genetics (Houston, TX), and GeneDx (Gaithersburg, MD), academic clinical laboratories, or through research studies. Additional details are provided in the individual clinical histories in the Additional file 1: Supplementary Note.
Facial analysis
Face2Gene Research application (FDNA Inc., Boston, MA) using DeepGestalt technology (algorithm 19.1.9) [21] was used to evaluate the presence of a distinct facial pattern in individuals with KCNK9 imprinting syndrome. Seventeen frontal photos of the face were obtained from unrelated affected individuals without glasses or black eye bars (individual: age at photo: P3.1: 9y, P6.1: 10y3 m, P7.1: 9y, P8.1: 1y5 m, P9.1: 12y-18y, P12.1: 2y2 m; P13.1: 6y1m, P16.1: 17y, P18.1 (updated photo patient 1) [8]: 6y, P23.1: 1y3 m, P24.1: 3 m, P27.1: 6y1m, P29.1: 8y1m and published photos of patient 2 (P19.1: 11 m), 3 (P20.1: 1y1m), and 4 (P21.1: 3y) from Graham et al. [8] and the individual (P25.1: 17y) described by Šedivá et al. [6] and compared to 17 control images matched for age, sex, and ethnicity provided by Face2Gene. To estimate the power of DeepGestalt in distinguishing affected individuals from controls, a cross validation scheme was used, including a series of binary comparisons between all groups. For these binary comparisons, the data was split randomly multiple times into training sets and test sets. Each such set contained half of the samples from the group, and this random process was repeated 10 times. The results of the comparisons are reported using the receiver operating characteristic (ROC) curve and area under the curve (AUC). The mean AUC = 0.87 and AUC STD = 0.07.
Molecular modeling
No experimentally solved structure currently exists for the protein encoded by KCNK9, TASK3. Therefore, we began molecular modeling from the canonical UniProt sequence, Q9NPC2-1, which corresponds to Ensembl transcript ENST00000303015 to search for existing experimental structures of homologous sequences using Clustal [22, 23] alignment to the Protein Data Bank (PDB) [24]. Using I-tasser [25], we generated homology-based models for the highest homology experimental structures - residues 1-271 from human KCNK1 (PDB: 3ukm) [26] and KCNK4 (PDB: 3um7) [27]. These models were evaluated for quality using online servers and standard metrics including PROCHECK [28], QMEAN [29], QMEANBrane [30], and VADAR [31]. We compared models to one another by calculations of their electrostatic potentials, volumes, and accessible surface areas using DaliLite [32] and APBS [33]. We summarized quality metrics at the residue level to evaluate if differences in quality clustered in 3D. All structure metrics were within expectation for high-confidence membrane proteins. Regions that are embedded within the membrane have high quality scores, comparable to solvent-exposed loops on both intracellular and extracellular-facing sides, which have lower quality scores. We believe that these loops will be more flexible, and therefore, any individual static representation of them will be insufficient. Our use of replicates and comparison to wildtype (WT) will help control for errors in the initial placement of these loops and facilitate interpreting the effects of genomic variants.
After completion of the 3D molecular modeling work described above, an experimental structure for TASK1 (KCNK3) was released [34]. This protein is closer in sequence identity to TASK3 (62%) compared to the TWIK1 (KCNK1) and TRAAK (KCNK4) experimental structures used as templates in this study (27% and 28%, respectively) [26, 27]. In our model, all amino acids through Val243 are concordant with the TASK1 experimental structure. After Val243, our model was disordered and highly dynamic in simulations, concordant with the sequence not conforming to the C-terminal of TWIK1 or TRAAK. In TASK1-based models, Met249 faces into the membrane and we computed Met249Thr to be destabilizing, consistent in effect with our original analysis (see the “Results” section). Thus, our modeling was highly concordant with the TASK1-based model and we have interpreted variants using information from both.
An explicit molecular environment was generated using Visual Molecular Dynamics (VMD) [35] and CHARMM-GUI [36, 37]. An 80x80Å membrane bilayer patch of POPC was generated. The protein was oriented to the patch using sequence-based annotation of trans-membrane helices from UniProt and consensus annotations from MESSA [38]. Water was added to provide a 10 Å distance between TASK3 and environment boundaries. The environment was neutralized at 150 mM KCl by adding 75 K+ and 79 Cl− ions. This initial explicit-environment system for WT TASK3 was used in further modeling, as described below.
Structure-based calculations
We scored genomic variants by 3D structure-based algorithms. First, we generated dimer models of each genomic variant using FoldX [39, 40] version 4 for computational mutagenesis and calculating folding stability changes upon mutation (ΔΔGfold). The explicit-environment WT TASK3 system was modified for each genomic variant using the FoldX-generated models. We also assessed the initial models using Frustratometer [41], to quantify how each genomic variant changes local interactions, from a more integrated perspective that accounts for the many concurrent types of favorable and un-favorable interactions of the mutated residues.
Molecular dynamics simulations
We studied the dynamics of KCNK9 genetic variants on TASK3, and how they affect ionic interactions, using molecular dynamics (MD) simulation as implemented in the NAMD [42] software and the CHARMM27 force field [43]. Each protein model was analyzed in triplicate, beginning with energy minimization for 10,000 steps. Next, harmonic constraints were added to protein non-hydrogen atoms and the system was heated to 300 K over 300 ps using a Langevin thermostat. Protein constraints were gradually released over 1 ns. The entire system was simulated free of restraints for a further 32 ns. In total, we generated 110 ns of MD trajectory per variant. Across variants and replicates, we generated 2.3 μs of MD trajectory to study the effects of the variants in TASK3.
Molecular dynamics analysis
We calculated root mean-squared deviation (RMSD) and root mean-squared fluctuation (RMSF) values using Cα atoms after structurally aligning each trajectory to the initial WT conformation and ignoring the mobile C-terminus. We used principal component analysis (PCA) to summarize the dominant conformational changes across trajectories. PCs were calculated in Cartesian coordinates using Cα atoms. We calculated the radial distribution function (RDF) for ions using VMD. The RDF of K+ describes the normalized probability of observing K+ ions within concentric shells. Electrostatic surfaces were sampled over time by extracting periodic frames from each simulation and assessed them each using APBS [33, 44]. Analyses were carried out using a custom structural bioinformatics workflow and leveraging the bio3d R package [45]. Protein structures and trajectories were visualized using PyMOL [45, 46] and VMD [35]. We generated our time-dependent prediction of effect for each genomic variant using a manual synthesis across multiple structure- and dynamics-based metrics as well as manual inspection of the simulation trajectories.
Mammalian expression plasmids
Human TASK3 (GenbankTM AF212829) cDNA, was cloned into the pcDNA3.1+ vector (Invitrogen, Carlsbad, CA, USA), gifted by Helen Meadows (GlaxoSmithKline, Harlow, UK) or into the pAcGFP1-N1 fluorescent vector (Clontech-Takara Bio Europe). Plasmids encoding for the human M3 muscarinic acetylcholine receptor (GenbankTM AF498917) cDNA were obtained from the UMR cDNA Resource Center (Rollo, MO).
Mutagenesis
Each of the clinically identified KCNK9 variants (Arg131Ser, Arg131His, Arg131Pro, Met132Arg, Phe135del, Met156Val, Met159Ile, Phe164Cys, Thr199Ala, Tyr205Cys, Gly236Arg, Ala237Asp, Met249Thr, Ala320Thr) were introduced by site-directed mutagenesis into human TASK3 cDNA using the QuikChange kit (Agilent, CA, USA) as previously described [14]. Oligonucleotide primers were synthesized by Eurofins MWG Operon, Ebersberg, Germany, and all constructs sequenced by DNA Sequencing & Services, MRC/PPU, University of Dundee, Scotland.
Cell culture
All experiments were performed using a modified human embryonic kidney 293 cell line, tsA201 (European Collection of Authenticated Cell Cultures; Sigma-Aldrich, UK), prepared and maintained as previously described [14]. Once cells reached 80% confluency, they were split and resuspended in media at a density of 7 × 104 and 0.5 mL transferred to a 4-well plate containing a poly-D-lysine-coated coverslip, ready for transfection the following day.
Transfection
Plasmids containing cDNA for either WT or a mutated TASK3 variant and a similar plasmid encoding the cDNA for green fluorescent protein (GFP) were co-transfected at a concentration of 0.5 μg using a modified calcium-phosphate protocol, as previously described [14]. The cells were then incubated at 37 °C in 95% O2 and 5% CO2 for 4–6 h, before being washed with phosphate-buffered saline. The cells were used for experiments 18–24 h later. All variants were expressed as homodimeric channels (each α-subunit of the dimer expresses the incorporated mutation) for all experiments. For experiments measuring the effect of G-protein coupled receptors, M3 receptors were also included in the transfection at a concentration of 0.5 μg.
Whole-cell patch-clamp electrophysiology
Currents were recorded from GFP-fluorescing tsA201 cells expressing the cDNA of interest using whole-cell patch-clamp in a voltage clamp configuration and a step-ramp voltage protocol as previously described [14] using an extracellular solution composed of 145 mM NaCl, 2.5 mM KCl, 3 mM MgCl2, 1 mM CaCl2 and 10 mM HEPES (pH adjusted to 7.4) and an intracellular pipette solution of 150 mM KCl, 3 mM MgCl2, 5 mM EGTA, and 10 mM HEPES (pH adjusted to 7.4). All experiments were conducted at room temperature (20–25 °C), and currents were recorded using an Axopatch 1D patch clamp amplifier (Molecular Devices, Sunnyvale, CA), filtered at 2 kHz, digitized at 5 kHz. Extracellular control solution and modulatory compounds were perfused at a rate of 4–5 mL min−1.
Electrophysiology data analysis and statistics
Data analysis of whole-cell outward current and analysis software was as previously described in Cunningham et al. [47]. Current-voltage graphs were obtained from the voltage ramp (− 120 mV to + 20 mV). Data were expressed as the mean ± 95% confidence intervals (CI), and n represents the number of individual cells recorded. Statistical analyses used were either a one-way ANOVA with a post-hoc Dunnett’s multiple comparisons test or an unpaired/paired Student’s t-test. Data was considered statistically different if P < 0.05 (*), P < 0.01 (**), P < 0.001 (***), and P < 0.0001 (****). Data from variants was compared with matched control data from WT TASK3 recorded either simultaneously or around the same calendar period and cell batch number.
Confocal microscopy
Colocalization studies were performed according to the detailed methods outlined in Cunningham et al. [47]. In brief, WT human TASK3 and human TASK3_Tyr205Cys cDNA’s were subcloned into the pAcGFP1-N1 fluorescent vector (Clontech-Takara Bio Europe) to create a fusion construct with GFP and expressed in tsA201 cells using TurboFect transfection reagent (ThermoFisher, Loughborough, UK). Prior to fixing the cells, the plasma membranes of the cells were stained with CellMask Deep Red (CMR) (ThermoFisher, UK) and the nuclei with Hoechst 33528 (Sigma-Aldrich, UK). Coverslips of cells were mounted with Vectashield anti-fade mounting medium (Vector Laboratories, UK).
Confocal microscopy images were taken using a LSM 880 confocal microscope (Carl Zeiss, Oberkochen, Germany) and processed using Zen Black software (Carl Zeiss). Cells were excited with an argon laser at either 561, 488, or 405 nm for the CMR plasma membrane stain, pAcGFP fused channel and Hoechst 33528 stained nuclei, respectively.
Colocalization was determined using Zen Black software and Pearson’s correlation coefficient (PCC) was used to represent the degree of co-localization.
Chemicals
Muscarine chloride, zinc chloride, dithiothreitol, 5-dithio-bis(2-nitrobenzoic acid), and methanethiosulfonates were all purchased from Sigma-Aldrich, UK. For high potassium (25 mM K) or acidic pH (pH 6.4) experimentation, the extracellular solution was adjusted accordingly. Muscarine chloride was added directly to extracellular solution at a concentration of 0.1 μM prior to use.
Functional study correlation and variant pathogenicity classification
To appropriately weight the strength of functional evidence provided by computational modeling and in vitro electrophysiological studies, we generated a weighting scheme for each. For computational modeling, we primarily used four dynamic features as follows (see Additional file 2: Table S1): The “Δ3-HB bbn RMSD” is the change in protein backbone RMSD among the 3-helix bundles. Since the WT channel exhibits dynamics in simulations (median of 2.8 Å RMSD), we binned this metric into three levels based on how distorted the 3HBs were for each variant compared to WT: ≥ 0.3 ΔRMSD, within 0.3 ΔRMSD, and ≤ − 0.3 ΔRMSD compared to WT simulations. The “SSR of Filter Cα RMSF” is the sum of squared residuals for the selectivity filter alpha-carbon RMSF, between WT and each variant. We binned this metric into three levels compared to WT: greater than 1.0, between 1.0 and 0.5, and less than 0.5. The selectivity filter depends not only on the backbone, but also on sidechain orientations. So, we quantified the “Δside chain heavy atom RMSD TIG (Thr-Ile-Gly) motifs” binned into greater than 0.5, within 0.5, and less than 0.0. Finally, we used “anti-chamber open/close by PCs” to indicate if a single PC feature indicated opening the antechamber, two PC features indicated opening, no change, a single feature of more closed antechamber, and two features of more closed antechamber (Additional file 3: Video S1). For the electrophysiological studies, we used a two-step approach applying more weight to current than the other measured channel attributes. We first weighted each assay result: current amplitude (0–3 corresponding to the level of statistical significance (0: not significant; 1: 0.001 < P < 0.01; 2: 0.0001 < P < 0.001; 3: P < 0.0001) and regulation by GPCR or by pH (0: not significant; 1: P < 0.05) and added weights for each variant (max weight = 5). These were then rescaled to a 0–3 scale (0: sum weight = 0; 1: sum weight = 1; 2: sum weight = 2–3; 3: sum weight > 3). The only exception to this was the Tyr205Cys variant which had no measurable current but would have only received a sum weight of 3 since regulation by GPCR and pH could not be measured without any baseline current, but the damaging evidence was more significant than that weight alone.
Variants were classified according to the 2015 Guidelines by the American College of Medical Genetics and Genomics and Association for Molecular Pathology guidelines [48]. Per current ClinGen guidelines, PM2 was used as a supporting level of evidence [49]. Conservatively, PS2 was applied as a moderate level of evidence due to KCNK9 being a paternally imprinted locus. PP3 was determined using CADD [50], MutationTaster [51], PrimateAI [52], SIFT [53], and Polyphen2 [54] (dbNSFP version 4.1 [55]) applying supporting level of evidence when 3 or more tools predict a damaging impact (Additional file 2: Table S1). The weight for functional evidence (PS3) was applied using only the electrophysiological data by the rescaled weight (none = 0; supporting = 1; moderate = 2; strong = 3). Agreement between computational and in vitro impact of variants was assessed using Spearman correlation and linear regression models.