Identifying mutation hotspots reveals pathogenetic mechanisms of KCNQ2 epileptic encephalopathy

Kv7 channels are enriched at the axonal plasma membrane where their voltage-dependent potassium currents suppress neuronal excitability. Mutations in Kv7.2 and Kv7.3 subunits cause epileptic encephalopathy (EE), yet the underlying pathogenetic mechanism is unclear. Here, we used novel statistical algorithms and structural modeling to identify EE mutation hotspots in key functional domains of Kv7.2 including voltage sensing S4, the pore loop and S6 in the pore domain, and intracellular calmodulin-binding helix B and helix B-C linker. Characterization of selected EE mutations from these hotspots revealed that L203P at S4 induces a large depolarizing shift in voltage dependence of Kv7.2 channels and L268F at the pore decreases their current densities. While L268F severely reduces expression of heteromeric channels in hippocampal neurons without affecting internalization, K552T and R553L mutations at distal helix B decrease calmodulin-binding and axonal enrichment. Importantly, L268F, K552T, and R553L mutations disrupt current potentiation by increasing phosphatidylinositol 4,5-bisphosphate (PIP2), and our molecular dynamics simulation suggests PIP2 interaction with these residues. Together, these findings demonstrate that each EE variant causes a unique combination of defects in Kv7 channel function and neuronal expression, and suggest a critical need for both prediction algorithms and experimental interrogations to understand pathophysiology of Kv7-associated EE.


MHF algorithm identifies EE mutation clusters in K
). These variants include 10 submicroscopic and partial gene deletions, 17 splice site mutations, 10 nonsense mutations, 25 frameshift mutations, 2 non-initiation mutations, 126 missense mutations that lead to single amino acid substitutions, and 4 mutations that result in single amino acid deletions (Fig. 1a). These mutations were classified into three groups according to the severities of their clinical outcomes described in the RIKEE database (Fig. 1a, Supplementary Table S1). The "mild or BFNE" mutations lead to seizures but not developmental delay in patients. The "severe or EE" mutations cause neonatal encephalopathy, seizures, and developmental delays. The "uncertain severity" mutations are associated with both benign seizures and EE or have limited clinical information. In addition, 130 silent mutations (with highest allele frequency from 0.0017% to 19%) and 25 relatively common nonpathogenic missense mutations of K v 7.2 (with allele frequency ≥ 0.01%) were identified from the Exome Aggregation Consortium (ExAC) database that collected protein-coding genetic variations from 60,706 humans (http://exac.broadinstitute.org).
In contrast to the evenly distributed silent mutations, pathogenic single amino acid mutations are concentrated at transmembrane segments S1 to S6, the pore loop, and intracellular helices A and B of K v 7.2 (Fig. 1b). To test if this trend was statistically significant, we developed a resampling algorithm titled Mutation Hotspot Finder (MHF). This algorithm was applied under the null hypothesis that pathogenic mutations are equally observed at every residue of a functional domain in the full-length K v 7.2 protein when there is no further association between the mutations and the domains. Because our MHF examines the association between the pathogenic variants and the functional domains, we used 130 single amino acid mutations and excluded nonsense and frameshift mutations that truncate one or more functional domains in K v 7.2. This analysis revealed that epilepsy mutations are significantly clustered at the voltage sensing S4, the pore loop, and S6 of K v 7.2 (p < 0.005), whereas silent and nonpathogenic mutations did not cluster at any of the functional domains (Fig. 1b, Supplementary Table S2). Importantly, epilepsy mutations of K v 7.2 were significantly associated with the "severe or EE" group (p < 0.001) and not the "mild or BFNE" and "uncertain severity" groups ( Fig. 1a-e, Supplementary Table S3).
Our MHF analysis also revealed that helix B and helix B-C linker have significantly more pathogenic mutations (p < 0.01) than other domains within the K v 7.2 C-terminal tail due to the clustering of "severe or EE" mutations (p < 0.05) (Fig. 2a-d, Supplementary Tables S2-3. Since K v 7.2 binds to CaM through helices A and B (Fig. 2e) 28, 32 , we next tested if the clinical severity of epilepsy mutations is associated with the extent to which K v 7.2 variants bind to CaM. Both "mild or BFNE" and "severe or EE" mutations located at helices A and B decreased the CaM binding energy of K v 7.2 (Fig. 2e,f). Furthermore, EE mutations occur at the positively charged residues in the distal portion of helix B and the helix B-C linker away from the CaM contact site in the modeled highlighted with colored spheres on the C-alpha atoms of one subunit: mild or BFNE (blue), uncertain severity (purple), severe or EE (red). Where more than one mutation occurs at a single position, the residue is colored by the most severe phenotype. (e) Location of amino acids mutated in mild or BFNE (blue), uncertain severity (purple), or severe or EE (red) in S4, the pore loop and S6 of K v 7.2.
Voltage-dependent activation of homomeric K v 7.2 channel is disrupted by selected EE mutations in epilepsy mutation hotspots. To test if EE variants within the mutation hotspots disrupt key functional protein domains of K v 7.2, we selected four EE mutations which have not been previously characterized: L203P at the voltage-sensing S4 23 , L268F at the pore loop 26 , and K552T and R553L at helix B 22,24 (Fig. 3a,b). To determine their effects on voltage-dependent activation of homomeric K v 7.2 channels, we performed whole-cell patch clamp recording in Chinese hamster ovary (CHOhm1) cells, which display very low expression of endogenous K + channels and depolarized resting membrane potential of −10 mV 12,33 . Application of depolarizing voltage steps from −100 to +20 mV in GFP-transfected CHOhm1 cells produces very little voltage-dependent currents that reverse around −26 mV 12,34 . In contrast, the same voltage steps in cells transfected with GFP and K v 7.2 wild-type (WT) generated slowly activating voltage-dependent outward K + currents that reached peak current densities of 17.3 ± 1.1 pA/pF at +20 mV (Fig. 3d,e, Supplementary Fig. S1). The average V 1/2 of WT channels (−26.8 ± 2.1 mV) was similar to the previously published value of −25 ± 1.9 mV 35 . Consistent with increased outward K + current, cell expressing K v 7.2 displayed hyperpolarized resting membrane potential (−35.5 ± 1.1 mV) and reversal potential (−38.8 ± 1.9 mV) (Supplementary Tables S4-5. Cells expressing GFP and K v 7.2-L203P produced K + currents with a large depolarizing shift in their voltage dependence and V 1/2 and increased activation time constants, decreasing peak current densities at voltage steps up to 0 mV. These cells also displayed depolarizing resting membrane potential (−26.8 ± 2.0 mV) and reversal potential (−22.1 ± 1.2 mV) ( Fig. 3d-f, Supplementary Fig. S1, Supplementary Tables S4-5. Surprisingly, their peak current density at +20 mV was larger (27.6 ± 1.88 pA/pF) than that of WT channels despite their slower activation kinetics (Fig. 3c-e, Supplementary Fig. S1). The L268F mutation in the pore loop decreased outward K + currents through K v 7.2 channels but not their protein level ( Fig. 3c-f, Supplementary Fig. S1). While the R553L mutation in distal helix B had no effect on K v 7.2 channels, the K552T mutation reduced both protein and current expression ( Fig. 3c-f, Supplementary Fig. S1). The L268F and K552T mutations did not alter voltage-dependence, activation kinetics, and reversal potential of K v 7.2 currents (Fig. 3g,h, Supplementary Tables S4-5). pip 2 -induced potentiation of K v 7.2 current is blocked by selected EE variants. PIP 2 is a critical cofactor required for the opening of K v 7 channels 14,16,17,36 and is proposed to bind to the intracellular side of S4, the S2-S3 and S4-S5 linkers, and intracellular region from pre-helix A to the helix B-C linker [11][12][13][14]28,[36][37][38] . Therefore, we next tested if selected EE mutations alter gating modulation of K v 7 channels by PIP 2 . Previous studies have shown that the activation of K v 7 channels is far from saturated by the endogenous membrane level of PIP 2 39 and that supplying exogeneous PIP 2 can enhance single-channel open probability and whole-cell current densities of homomeric K v 7.2 channels 12,14,37,40 .
Indeed, inclusion of diC8-PIP 2 (100 μM) in the intracellular pipette solution increased K + currents through K v 7.2 WT channels by 2-fold and caused a modest left shift in voltage-dependence ( Fig. 3d-f, Supplementary Figs. S1-2) as previously shown 12,14,37 . Surprisingly, all selected EE mutations abolished diC8-PIP 2 -induced potentiation of K v 7.2 channels and hyperpolarizing shift in their voltage dependence, resulting in a significant reduction in their current densities compared to WT channels in the presence of diC8-PIP 2 ( Fig. 3d-f, Supplementary Figs. S1-2, Supplementary Table S5).
To increase cellular PIP 2 levels, we transfected phosphatidylinositol-4-phosphate 5-kinase (PIP5K), which catalyzes the formation of PIP 2 via the phosphorylation of phosphatidylinositol-4-phosphate 41 . Consistent with previous reports 11,40,42 , co-expression of PIP5K increased K + currents through K v 7.2 WT channels with a hyperpolarizing shift in their voltage dependence. Consistent with the recording with diC8-PIP 2 inclusion (Fig. 3), this effect was absent in K v 7.2 channels containing L268F, K552T, and R553L variants ( Supplementary Fig. S3), indicating that these mutations abolished current potentiation of K v 7.2 channels upon increasing cellular PIP 2 levels.
Modeled K v 7.2 structure and molecular dynamics simulation suggest that selected EE mutations reside in pip 2 binding regions of K v 7.2. To investigate if selected EE mutations are located in PIP 2 -binding regions, we compared our modeled K v 7.2 structure bound to CaM and the published structure of TRPV1 channel embedded in lipid nanodiscs with phosphatidylinositol bound (PDB: 5irz) (Fig. 4a) 43 . In the modeled K v 7.2 structure, the voltage-sensor (S1-S4) and the pore domain of K v 7.2 form the hydrophobic cavity where L203 and L268 are located (Fig. 4a). Similar to the binding of phosphatidylinositol to TRPV1 channel ( Fig. 4a) 43 , the fatty acid tails of amphiphilic PIP 2 are most likely embedded in this hydrophobic cavity of K v 7.2 where L203P and L268F mutations reside. Furthermore, the bottom of the voltage-sensor (S1-S4) together with proteins. For clarity, cropped gel images are shown. Full-length gels can be found in Supplementary Fig. S7,a. (d) Representative recordings after subtraction of leak currents. Leak current was defined as non-voltagedependent current from GFP-transfected cells. (e) Average peak current densities at all voltage steps. *p < 0.05, ***p < 0.005 based on one-way ANOVA Fisher's test. (f) Average peak current densities at -20 mV (left) and + 20 mV (right). p values are computed from one-way ANOVA Tukey test. (g) Normalized conductance (G/ Gmax) at all voltage steps. (h) Activation time constant (τ) at + 20 mV. The number of GFP-cotransfected cells that were recorded without diC8-PIP 2 : K v 7.2 WT (n = 12), L 2 03P (n = 17), L268F (n = 17), K552T (n = 13), or R553L (n = 13). The number of GFP-cotransfected cells that were recorded with diC8-PIP 2 : K v 7.2 WT (n = 11), L 2 03P (n = 14), L268F (n = 13), K552T (n = 11), or R553L (n = 11). Data shown represent the Ave ± SEM. (2020) 10:4756 | https://doi.org/10.1038/s41598-020-61697-6 www.nature.com/scientificreports www.nature.com/scientificreports/ pre-helix A, helix B and the helix B-C linker of K v 7.2 form a highly basic environment favorable for binding the phosphate headgroup of PIP 2 (Fig. 4a), consistent with previous studies in K v 7.1 28 .
To test if PIP 2 interacts with K552 and R553 in distal helix B, we performed molecular dynamics (MD) simulation. We constructed a homology model of the CaM-bound closed state conformation of K v 7.2 using the structure of K v 7.1 (PDB: 5VMS) 28 as a template, and employed targeted MD to model the open-state conformation of K v 7.2 in the explicit lipid bilayers containing 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and PIP 2 lipids. To extensively sample the lipid-protein interactions, we constructed two independent simulation systems, each containing seven PIP 2 molecules randomly placed around K v 7.2 without CaM in both outer and inner membrane leaflets (~2.2% PIP 2 ) (Fig. 4b). Within the time frame of the simulations, PIP 2 molecules diffused from the periphery of the K v 7.2 structure towards its central region (Fig. 4c).
To examine PIP 2 binding to helix B, we measured the distance between the center of mass of R553 from each monomer and that of the phosphate groups at position 4 and 5 of PIP 2 (Fig. 4d). Binding of PIP 2 molecules towards R553 was observed in 2 out of 4 monomers in the first simulation and in 3 out of 4 monomers in the second simulation (Fig. 4d). In both simulations, PIP 2 molecules interacted with K552-R553-K554 within 100 ns and remained stably bound throughout the simulations (Fig. 4e, Supplementary Video 1, Supplementary Fig. S4), consistent with previous in vitro biochemical studies and molecular docking simulations that demonstrated PIP 2 binding to the corresponding residues in the C-terminal helix A-B fragments of K v 7.1 37 . These findings suggest that K552T and R553L mutations are located in helix B of K v 7.2 that interacts with the phosphate head group of PIP 2 .
To test if selected EE mutations alter PIP 2 affinity, we examined the K v 7.2 current decay upon PIP 2 depletion induced by activation of voltage-sensitive phosphatase (VSP) 11,42,44 . In CHOhm1 cells coexpressing danio rerio VSP 11 , the 10s-depolarization step at voltages from +40 mV decreased peak currents of K v 7.2 channels, reaching the maximal decay of 53.2 ± 4.0% at +100 mV ( Supplementary Fig. S5). Current decay of K v 7.2-K552T was greater than that of WT at +40 mV but was comparable to that of WT from +60 to +100 mV ( Supplementary  Fig. S5), suggesting that the K552T mutation modestly decreased PIP 2 affinity to K v 7.2. Interestingly, the same depolarization steps delayed current decay of K v 7.2 channels containing L203P, L268F, and R553L mutations ( Supplementary Fig. S5), indicating their reduced sensitivity to PIP 2 depletion.

Selected EE variants decrease current expression of heteromeric K v 7 channels and their current potentiation by diC8-PIP 2 .
Since KCNQ2-associated EE is an autosomal dominant epileptic syndrome, we repeated voltage-clamp recording in CHOhm1 cells transfected with plasmids for K v 7.3, wild-type K v 7.2, and mutant K v 7.2 at a 2:1:1 ratio as described 35 Table S5). The L268F variant also increased their activation kinetics (Fig. 5e). None of the tested mutations affected total protein expression of K v 7.2 and K v 7.3 ( Fig. 5f).
When diC8-PIP 2 was added in the intracellular pipette, all tested EE mutations significantly decreased current densities of heteromeric channels at +20 mV compared to WT without altering their voltage-dependence ( Fig. 5a-d, Supplementary Table S5) and their activation time constant was increased by L203P and K552T mutations (Fig. 5e). Importantly, all tested EE variants abolished PIP 2 -induced current potentiation of heteromeric channels (

Selected EE variants variably decrease axonal surface expression of heteromeric K v 7 channels.
The physiologically relevant current through K v 7 channels is controlled by both their function and expression at the neuronal plasma membrane. Given that K v 7.2 interaction with CaM and K v 7.3 are critical for axonal surface expression of K v 7 channels 9,45 , we next tested if selected EE variants of K v 7.2 affect interaction with CaM and K v 7.3 and axonal targeting of K v 7 channels (Figs. 6-7, Supplementary Figs. S9-11). Coimmunoprecipitation assay in HEK293T cell lysate 12,45 revealed that the K552T and R553L mutations in helix B decreased K v 7.2 binding to YFP-tagged CaM whereas the mutations including L203P in S4 and L268F in the pore loop had no effect ( Fig. 6a,b). None of the tested mutations affected K v 7.2 interaction with K v 7.3 ( Fig. 6c,d). Total K v 7.2 expression was also reduced by the L203P and K552T variants in cells co-expressing CaM but not K v 7.3 ( Fig. 6).
To test if selected EE mutations of K v 7.2 affect surface density of K v 7 channels, we transfected rat dissociated hippocampal cultured neurons with K v 7.3 containing an extracellular HA epitope (HA-K v 7.3) and performed surface immunostaining of HA-K v 7.3 as described 9,12,45 (Fig. 7, Supplementary Fig. S11). In cultured neurons, transfection of HA-K v 7.3 alone yields negligible surface expression of HA-K v 7.3 9 . However, co-transfection of K v 7.2 WT results in robust HA-K v 7.3 expression on the plasma membrane of the AIS and distal axons compared to the soma and dendrites (Fig. 7a,b) 9,12,45 , resulting in a surface fluorescence "Axon/Dendrite" ratio of 3.9 ± 0.49 (Fig. 7c).
Although the L203P mutation in S4 did not affect surface and total expression of HA-K v 7.3/ K v 7.2 channels ( Fig. 7a-d), the L268F mutation in the pore loop abolished their preferential enrichment at the axonal surface by severely decreasing their axonal surface density (surface Axon/Dendrite ratio = 0.85 ± 0.10, Fig. 7a-c) and also reduced intracellular K v 7.2 expression in the axon (Fig. 7d). The K552T and R553L mutations in helix B significantly reduced surface expression of heteromeric channels in both distal axon and dendrites (Fig. 7f-i), resulting in similar surface Axon/Dendrite ratios as the WT channels (Fig. 7h).
Disruption of CaM binding to K v 7.2 has been shown to impair axonal enrichment of K v 7 channels by inhibiting their trafficking from the endoplasmic reticulum (ER) 45 . The ability of the L268F mutation to impair axonal K v 7 surface expression without affecting K v 7.2 binding to CaM or K v 7.3 (Figs. 6, 7a-d) suggests a different mechanism. To test if the L268F mutation reduces axonal enrichment of K v 7 channels by increasing their endocytosis, we used dynamin inhibitory peptide (DIP, 50 μM) which blocks dynamin-dependent endocytosis in cultured hippocampal neurons 46 . The DIP treatment for 45 min induced a small increase in surface HA-K v 7.3/ K v 7.2 WT and L268F mutant channels in the soma and dendrites but not axons (Fig. 7e,j), indicating their basal endocytosis in somatodendritic membrane. Although the DIP treatment had no effect on K552T mutant channels, the same treatment modestly increased axonal surface expression of L268F mutant channels (Fig. 7e,j) www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ increase did not reach the axonal level of WT channels (Fig. 7e), suggesting that increased endocytosis is not the main cause for reduced axonal surface expression of L268F mutant channels.

Discussion
In this study, we investigated the pathogenetic mechanisms underlying de novo EE mutations of K v 7.2. Visual inspection in K v 7.2 primary sequence has suggested the enrichment of EE variants at S4, the pore domain from S5 to S6, and helices A and B 12,25,47 . Clustering of epilepsy mutations in the ion transport domain of K v 7.2 has also been detected by identifying its variation-intolerant genic sub-regions 48 . Our novel MHF statistical algorithm interpreted in the context of modeled K v 7.2 atomic structure (Figs. 1-2) supports these earlier observations. We discovered that "severe or EE" missense variants cluster at S4, the pore loop that contains the selectivity filter, S6, helix B, and the helix B-C linker of K v 7.2 ( Fig. 1). A recent study by Goto et al., reported that the EE missense variants cluster at the pore domain, S6, and pre-helix A of K v 7.2 49 . The regional differences in mutation clusters between our study and Goto et al., could be attributed to the use of different algorithms and databases (ExAC and GnomAD) as sources for non-pathogenic mutations. Nonetheless, both studies identified the pore domain and S6 as hotspots of EE variants, supporting the functional importance of these regions.
However, sequence variant interpretation from the prediction algorithms should be used carefully. The presence of both gain-of-function and loss-of-function EE variants in S4 47,50-52 suggest that it is not straight forward to predict the genotype-phenotype correlation of EE. In addition, both EE and BFNE variants exist in each of the identified hotspots and even at the same codon 18,49 , suggesting that different amino acid substitutions at the same   49 , the in vivo impact of a mutation is difficult to predict in patients due to their variable exposures to genetic and environmental factors. Thus, the use of multiple in-silico tools and comprehensive experimental analyses of epilepsy variants are needed to understand their effects on K v 7 channels ex vivo and in vivo.
Our functional characterization of de novo EE variants selected from the mutation hotspots revealed that each mutation impaired the function of its associated protein domain within K v 7.2. The L203P mutation in the main voltage sensor S4 induced a large depolarizing shift in voltage-dependence and slowed activation kinetics of homomeric K v 7.2 channels (Fig. 3) but had no effect on heteromeric channels (Fig. 5). In contrast, the L268F mutation in the pore loop decreased current densities of both homomeric and heteromeric channels without affecting their voltage dependence (Figs. 3, 5). K552T and R553L mutations in CaM-binding helix B decreased the interaction between CaM and K v 7.2 (Fig. 6), which is shown to play critical roles in M-current expression and inhibition of hippocampal neuronal excitability 53 . Current suppression of homomeric channels is a common feature of EE variants of KCNQ2 54 . Given the overlapping distribution of K v 7.2 and K v 7.3 throughout the hippocampus and cortex 4 , the dominant negative functional effects of L268F and K552T variants on heteromeric channels (Figs. 3, 5) are expected to induce neuronal hyperexcitability and may underlie severe symptomatic EE with drug-resistant seizures, psychomotor delay, and profound intellectual disability 22,26 .
Interestingly, our modeled K v 7.2 structure revealed that the distal helix B and helix B-C linker come together with pre-helix A to form a positively charged surface close to the voltage sensor S1-S4 and the base of S6 (Fig. 2). Mutations of basic amino acid residues including H328C, R325G, and R333W at pre-helix A and R560W at the helix B-C linker of K v 7.2 have been shown to impair regulation of K v 7.2 currents by PIP 2 11,12,14 , which couples voltage sensor activation to the opening of the gate 28,36 . Our MD simulations revealed that K552 and R553 in distal helix B bind to the negatively charged head group of PIP 2 (Fig. 4). Importantly, K552T and R553L mutations impaired current enhancement of both homomeric and heteromeric channels upon acute or tonic increase in PIP 2 (Figs. 3, 5, Supplementary Fig. S3), suggesting that these mutant channels cannot respond to the changes in cellular PIP 2 . Since stable binding of CaM to K v 7.2 is crucial for PIP 2 modulation of neuronal K v 7 channels 33 , a decrease in CaM binding (Fig. 6) may also contribute to the loss of PIP 2 -induced current enhancement of K552T and R553L mutant channels (Figs. 3, 5).
The impairment of PIP 2 -induced current enhancement of L268F mutant channels was unexpected (Figs. 3, 5-7) because it is unlikely for the hydrophobic L268 to bind a negatively charged head group of PIP 2 . A comparison between modeled K v 7.2 structure and TRPV1 structure (Fig. 4a) suggests that the amphiphilic chains of PIP 2 may extend to the hydrophobic cavity created by the voltage-sensors S1-S4 and the pore domain of K v 7.2. We speculate that the L268F mutation at this hydrophobic interface impair K v 7.2 interaction with PIP 2 . Furthermore, analogous residue for L268 in the bacterial KcsA structure can secure the proper opening size of the pore 55 . Therefore, it is also possible that the L268F mutation may disrupt PIP 2 -dependent coupling to the pore opening 36 .
Several studies including our own have investigated PIP 2 affinity of K v 7 function by inclusion of diC8-PIP 2 in the intracellular pipette solution 12,14,33,37 . However, caution must be exercised in interpreting their results. Potentiation of K v 7.2-L203P current by tonic elevation of cellular PIP 2 upon PIP5K expression but not acute application of diC8-PIP 2 (Fig. 3, Supplementary Fig. S3), suggest that diC8-PIP 2 inclusion may not readily potentiate the mutant channels that displayed very slow activation kinetics (Fig. 3). Furthermore, the loss of diC8-PIP 2 -induced current potentiation can be either due to decreased PIP 2 affinity or saturated level of interaction with PIP 2 at low PIP 2 concentration. We found that the K552T mutation modestly weakens PIP 2 affinity, whereas other mutant channels were resistant to PIP 2 depletion (Supplementary Fig. S5). Considering multiple proposed PIP 2 binding sites in K v 7.2 including S2-S3 and S4-S5 linkers, pre-helix A, helix B, and helix A-B and helix B-C linkers [11][12][13][14][15] , selected EE mutations may cause conformational change that weakens or enhances PIP 2 affinity to other regions within K v 7.2. As Suh and Hille (2008) pointed out 56 , it is not straight forward to determine PIP 2 affinity of mutant channels by assessing their currents after manipulation of PIP 2 level. Nonetheless, the lack of current potentiation upon increasing cellular or exogenous PIP 2 (Fig. 3, Supplementary Fig. S3) and the location of the EE mutated residues in a region of K v 7.2 that binds to fatty acid tails or polar headgroups of PIP 2 (Fig. 4) strongly suggest that there are multiple ways by which the selected EE variants may influence PIP 2 interaction with K v 7 channels and reduce their currents.
Our investigation of selected EE variants on neuronal expression of K v 7 channels revealed that K552T and R553L mutations in helix B reduced enrichment of K v 7 channels at the axonal surface (Fig. 7), supporting previous observations that the degree of CaM interaction with K v 7.2 correlates with the overall amount of K v 7 channels at the axonal surface 12,45 . Unexpectedly, the L268F mutation at the pore loop severely decreased both surface and intracellular expression of heteromeric channels in axons without affecting K v 7.2 binding to CaM or K v 7.3 (Figs. 6-7), demonstrating the importance of studying K v 7 expression in neurons. Decreased axonal expression of K v 7.2-L268F and minor effects of endocytosis inhibition (Fig. 7) suggest that a severe reduction of L268F mutant channels at the axonal surface is caused by a CaM-and endocytosis-independent mechanism. Given that of transfected neurons that were analyzed in Fig. 7i: K v 7.2 WT (n = 17), K552T (n = 23), R553L (n = 14), UT (n = 13). (e,j) Background-subtracted mean intensities of surface HA fluorescence from the transfected neurons treated with vehicle (CTL) control or dynamin inhibitory peptide (DIP). The number of transfected neurons that were analyzed in Fig. 7e: WT + CTL (n = 14), WT + DIP (n = 13), L268F + CTL (n = 8), L268F + DIP (n = 8). The number of transfected neurons that were analyzed in Fig. 7j: WT + CTL (n = 8), WT + DIP (n = 13), L268F + CTL (n = 6), L268F + DIP (n = 7). Data represent the Ave ± SEM (*p < 0.05, **p < 0.01, ***p < 0.005 against WT channels). (2020) 10:4756 | https://doi.org/10.1038/s41598-020-61697-6 www.nature.com/scientificreports www.nature.com/scientificreports/ misfolded membrane proteins are retained in the ER for chaperone-assisted refolding 57 , the L268F mutation may cause a folding defect that facilitates ER retention and disrupts forward trafficking of heteromeric channels to the axon.
Taken together, we identified EE mutation hotspots in K v 7.2 and discovered that each variant selected from these hotspots impairs the function of its associated protein domain and displays a combination of defects in voltage-and PIP 2 -dependent activation and axonal expression of K v 7 channels (Fig. 8). Such combinations of defects may decrease K v 7 current and its ability to inhibit neuronal excitability in neonatal brain 5,12 , as conditional deletion of K v 7.2 during embryonic development results in hippocampal and cortical hyperexcitability and spontaneous seizures in mice 7 . Continued optimization of prediction algorithms and experimental interrogations to understand pathophysiology of K v 7-associated EE will aid the development of better therapeutic strategies for this disease.

Materials and Methods
The resampling statistical algorithm. A resampling algorithm titled Mutation Hotspot Finder (MHF) was developed to search for mutation clusters that localize to the functional domains in human K v 7.2 (GenBank: NP_742105.1). The complete MHF algorithm can be found in GitHub repository (https://github.com/jerrycchen/ MutationHotspotFinder). The functional domains were annotated based on multiple published sources 3,28,30,32 and the RIKEE database (www.rikee.org). Briefly, the MHF algorithm compares the observed numbers of mutations against the expected numbers of mutations, and computes corresponding statistical significance through bootstrapping within each pre-specified protein functional domains. The following sections explain the MHF algorithm in detail.
The S denote the set of all observed single amino acid mutations for the whole sequence of protein X (e.g. the whole sequence of K v 7.2), or the subset of the sequence of protein X (e.g. the intracellular C-terminal tail of K v 7.  Table S2), D 1 is the number of mutations in the S1, and D 2 is the number of mutations in the S2, and etc. Among 9 functional domains in intracellular K v 7.2 C-terminal tail (Supplementary Table S2), D 1 is the number of mutations in the pre-Helix A, and D 2 is the number of mutations in the Helix A, and etc.
The MHF algorithm assumes that mutations are equally observed at each amino acid position within functional domains when there is no further association between the mutations and the domains. Due to this null hypothesis 58 , the application of MHF algorithm is restricted to single amino acid mutations, and is not suitable for mutations outside of the coding sequence as well as nonsense or frameshift mutations that delete one or more function domains. Under such assumption, we can randomly draw samples (i.e. bootstrap), with size = | | S , from the sequence X to construct the bootstrapped "mutation sets" in order to simulate the distribution of mutations. www.nature.com/scientificreports www.nature.com/scientificreports/ For iteration k of bootstrapping, the  S k ( ) denote the bootstrapped mutation set where k K {1, 2, , } ∈ … . In the context of this paper, we ran 10,000 iterations of bootstrapping (i.e. = K 10, 000). The ∼ D k ( ) j is defined as the number of mutations in  S k ( ) that fall into the functional domain j of protein sequence X. The empirical expected number of mutations within each domain j was constructed from The empirical P-values (P j ) were computed from a right-tailed test to measure the level of statistical significance on the proportion of the bootstrapped mutation sets that had more mutations than the observed mutation set S at each individual protein functional domain.
is the Indicator function. The computed P-values were adjusted for multiple comparisons (J times) using Bonferroni's correction. Mutations were visualized and mapped to K v 7.2 and K v 7.3 primary structures with MutationMapper (http://www. cbioportal.org/mutation_mapper.jsp Mapper). Fisher's Exact Test was implemented using the standard fisher. test() function in R (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/fisher.test.html).
Structure modeling and visualization. The S1-S6 sequence of K v 7.2 (R75-Q326) was threaded to the cryo-EM structure of Xenopus laevis K v 7.1 bound to CaM (PDB: 5VMS) 28 . The loops of K v 7.2 (E86-W91 and K255-T263) were rebuilt in FoldIt (https://fold.it/portal). The structure was relaxed in Rosetta software (https:// www.rosettacommons.org/software) using two rounds of rotamer sampling followed by side chain and backbone minimization, ending with minimization of all degrees of freedom while maintaining C4 symmetry. The lowest scoring decoy with Root mean square deviation (RMSD) < 2.0 Å was chosen as the final model. The amino acid residues mutated in BFNE and EE are indicated in the Rosetta-based model.
To model the interaction between CaM and K v 7.2 helices A and B, the helix A sequence of K v 7.2 (E322-V367) was threaded to the crystal structure of chimeric K v 7.3 helix A -K v 7.2 helix B in complex with Ca 2+ -bound CaM (PDB: 5J03) 32 . The structure was relaxed with Rosetta using two rounds of sequential rotamer, side chain and backbone minimization, followed by rigid body minimization. Mutations were made to the model in Rosetta followed by sequential rotamer, side chain, backbone, and rigid body minimization. The binding energy was calculated from 20 simulations. Structures were visualized using PyMOL 2.0 (Schrödinger, LLC). DNA Constructs and mutagenesis. EYFP-hCaM was a gift from Dr. Emanuel Strehler (Addgene plasmid # 47603). The plasmid pIRES-dsRed-PIPKIγ90 was a gift from Dr. Anastasios. Tzingounis (University of Connecticut) and was previously described 42 . Plasmids pcDNA3 with KCNQ3 cDNA (GenBank: NM004519) encoding K v 7.3 (GenBank: NP_004510.1), HA-K v 7.3, and KCNQ2 cDNA (GenBank: Y15065.1) encoding K v 7.2 (GenBank: CAA 75348.1) have been previously described 9,12,45 . Compared to the reference sequence of K v 7.2 (GenBank: NP_742105.1), this shorter isoform lacks 2 exons which do not harbor pathogenic variants to date. However, the amino acid numbering in the manuscript conforms to the reference sequence of K v 7.2 for clarity. Epileptic encephalopathy mutations (L203P, L268F, K552T, R553L) were generated using the Quik Change II XL Site-Directed Mutagenesis Kit (Agilent).
Electrophysiology. Whole cell patch clamp recordings in Chinese hamster ovary (CHO hm1) was performed as described 12 . To express homomeric K v 7.2 channels, cells were transfected with pEGFPN1 (0.2 μg) and pcDNA3-K v 7. . Leak-subtracted current densities (pA/pF), normalized conductance (G/Gmax), and channel biophysical properties were computed as described 12,35 with the exception that V 1/2 and the slope factor k were calculated as described 35,59 by fitting the plotted points of G/Gmax with a Boltzmann equation To examine the decline of K v 7.2 current upon activation of Dr-VSP, CHO hm1 cells were transfected with pDrVSP-IRES2-EGFP (0.5 μg) and pcDNA3-K v 7.2 WT or mutant (0.5 μg). The pDrVSP-IRES2-EGFP plasmid was a gift from Yasushi Okamura (Addgene plasmid # 80333). Voltage-clamp recording of K v 7.2 current upon depolarization-induced Dr-VSP activation was performed as described 60 with an external solution containing 144 mM NaCl, 5 mM KCl, 2 mM CaCl2, 0.5 mM MgCl2, 10 mM glucose and 10 mM HEPES (pH 7.4). Patch pipettes (3)(4) were filled with intracellular solution containing 135 mM potassium aspartate, 2 mM MgCl2, 1 mM EGTA, 0.1 mM CaCl2, 4 mM ATP, 0.1 mM GTP and 10 mM HEPES (pH 7.2). Cells were held at -70 mV and 10 s step depolarizations were applied in 20 mV steps from -20 to +100 mV with 2 min inter-step intervals to allow PIP 2 regeneration. The extent of K v 7.2 current decay upon Dr-VSP activation during 10 s depolarization was measured as the ratio of current at 10 s over peak current at each voltage step. (2020) 10:4756 | https://doi.org/10.1038/s41598-020-61697-6 www.nature.com/scientificreports www.nature.com/scientificreports/ Molecular dynamics simulation. For modeling of open and closed states of K v 7.2, the closed-state conformation of KCNQ2 in calmodulin-bound form was modeled based on the recent cryo-EM structure of K v 7.1 (PDB code 5VMS) 28 . Multiple sequence alignment of the template and KCNQ2 sequence was performed by using TCoffee web server (https://www.ebi.ac.uk/Tools/msa/tcoffee/). After the alignment, the homology model of closed-state conformation was built with MODELLER 61 . The stability of the closed-state conformation of K v 7.2 was tested by performing all-atom molecular dynamics (MD) simulations in explicit lipid bilayer.
In order to model the open-state conformation of K v 7.2, we performed non-equilibrium MD simulations. Using our stable closed-state conformation of K v 7.2, we performed 20-ns of Targeted MD (TMD) 62 simulations in an explicit lipid bilayer. TMD has been shown to drive the conformational changes by gradually minimizing the RMSD of S4-S5 and S6 helices of the closed-state conformation and the target structure which is K v 1.2/K v 2.1 in open conformation (PDB: 2R9R) 63 . As major structural changes occur in the pore region of the channel, we applied a restraint (force constant = 250 kcal/mol/Å) on the S4-S5 and S6 helices of each monomer to drive it towards the target state which was defined by a highly homologous K v 1.2/K v 2.1 channel in open-state conformation (PDB code 2R9R) 63 . The success of TMD was gauged by measuring the backbone RMSD of S4-S5 and S6 helices with respect to the target (Supplementary Fig. S6). Upon completion of TMD, all the structural restraints were released and the stability of the obtained open-state conformation of K v 7.2 was tested by performing MD simulations in explicit lipid bilayer (Supplementary Fig. S6).
All the MD simulations were performed with NAMD2.12 65 using CHARMM36m force field for lipid/protein 66 and a timestep of 2 fs. Long range electrostatic interactions were evaluated with particle mesh Ewald (PME) 67 and periodic boundary conditions were used throughout the simulations. Non-bonded forces were calculated with a cutoff of 12 Å and switching distance of 10 Å. During the simulation, temperature (T = 310 K) and pressure (P = 1 atm) (NPT ensemble) was maintained by Nose-Hoover Langevin piston method 68 . During pressure control, the simulation box was allowed to fluctuate in all the dimensions with constant ratio in the x-y (lipid bilayer) plane.
Immunoprecipitation. HEK293T cells were plated on 100 mm cell culture dishes (BD Biosciences, 2 × 10 6 cells per dish) and maintained in Minimal Essential Medium containing 10% Fetal Bovine Serum, 2 mM glutamine, 100 U/mL penicillin and 100 U/mL streptomycin at 37 °C and 5% CO 2 . At 24 hr post plating, the cells were transfected with plasmids (total 1.6 μg) containing K v 7.2 and EYFP-hCaM (1:1 ratio), using FuGENE6 transfection reagent (Promega). For coimmunoprecipitation studies of K v 7.2 and K v 7.3, the cells were transfected with K v 7.2 and K v 7.3 containing an extracellular hemagglutinin epitope (HA-K v 7.3) (1:1 ratio). At 48 h post transfection, the cells were washed with ice-cold PBS and lysed in ice-cold immunoprecipitation (IP) buffer containing (in mM): 20 Tris-HCl, 100 NaCl, 2 EDTA, 5 EGTA, 1% Triton X-100 (pH 7.4) supplemented with Halt protease inhibitors (Thermo Fisher Scientific). The lysate containing equal amount of proteins were first precleared with Protein A/G agarose beads (100 μL, Santa Cruz) for 1 hr at 4 °C, and then incubated overnight at 4 °C with Protein A/G-agarose beads (100 μL) and rabbit anti-K v 7.2 antibody (2 μg). This amount of anti-K v 7.2 antibody allowed us to immunoprecipitate the equal amount of K v 7.2 proteins and analyze the effects of mutations on the amount of co-immunoprecipitated EYFP-hCaM and HA-K v 7.3. After washing with IP buffer, the immunoprecipitates were eluted with SDS sample buffer by incubating at 75 °C for 10-15 min, and analyzed by western blot analysis with mouse anti-GFP (1:500 dilution), mouse anti-K v 7.2 (1:200 dilution), mouse anti-HA antibodies (1:500 dilution), and rabbit anti-GAPDH antibodies (1:1000 dilution). Antibodies used in coimmunoprecipitation and immunoblotting include anti-K v 7.2 (Neuromab, N26A/23), rabbit anti-K v 7.2 (Alomone, APC-050), anti-GFP, anti-HA, Immunocytochemistry. All procedures involving animals were reviewed and approved by the Institutional Animal Care and Use Committee at the University of Illinois Urbana-Champaign and conducted in accordance with the guidelines of the U.S National Institute of Health (NIH). Primary rat dissociated hippocampal cultured neurons prepared from 18-day old embryonic rats were plated on 12 mm glass coverslips (Warner Instruments, 10 5 cells per coverslip) coated with poly L-lysine (0.1 mg/mL). These neurons were maintained in neurobasal medium supplemented with B27 extract, 200 mM L-glutamine, and 100 U/mL penicillin and streptomycin in a cell culture incubator (37 °C, 5% CO 2 ). At 5 days in vitro (DIV), neurons were transfected with plasmids (total 0.8 μg) containing K v 7.3 with an extracellular hemagglutinin epitope (HA-K v 7.3) and wild-type or mutant K v 7.2 using lipofectamine LTX as described 12,45 .
Fluorescence and phase contrast images of transfected neurons were viewed using a Zeiss Axio Observer inverted microscope High-resolution gray scale images of healthy transfected neurons were acquired using a 20X objective with a Zeiss AxioCam 702 mono Camera and ZEN Blue 2.6 software and saved as 16-bit CZI and TIFF files. To compare the fluorescence intensity of the neurons transfected with different constructs, the images were acquired using the same exposure time within one experiment.
The image analyses were performed from the healthy transfected neurons using ImageJ Software as described 12,45 and excluded the transfected neurons with broken neurites or soma as well as regions where fasciculation or overlapping processes occurred. The axon was identified as a process that were labeled for the AIS marker 14D4, whereas the dendrites were identified as the processes that were absent for 14D4 in the transfected neurons. ImageJ software was used to trace the all major primary dendrites, the AIS (defined as the first 0-30 μm segment of the axon), and distal axon (defined as the segment between 50 and 80 μm from the beginning of the axon) as 1 pixel-wide line segments, and obtain their mean fluorescent intensities. The perimeter of the neuronal soma was also traced to obtain background-subtracted mean fluorescent intensities of the soma. Statistical analyses. All analyses are reported as mean ± SEM. Using Origin 9.1 (Origin Lab), the Student t test and one-way ANOVA with post-ANOVA Tukey and Fisher's multiple comparison tests were performed to identify the statistically significant difference with a priori value (p) < 0.05. The number of separate transfected cells for immunostaining and electrophysiology was reported as the sample size n.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.