Association of rare missense variants in the second intracellular loop of NaV1.7 sodium channels with familial autism.

Autism spectrum disorder (ASD) is a complex neurodevelopmental disorder often accompanied by intellectual disability, language impairment and medical co-morbidities. The heritability of autism is high and multiple genes have been implicated as causal. However, most of these genes have been identified in de novo cases. To further the understanding of familial autism, we performed whole-exome sequencing on five families in which second- and third-degree relatives were affected. By focusing on novel and protein-altering variants, we identified a small set of candidate genes. Among these, a novel private missense C1143F variant in the second intracellular loop of the voltage-gated sodium channel NaV1.7, encoded by the SCN9A gene, was identified in one family. Through electrophysiological analysis, we show that NaV1.7C1143F exhibits partial loss-of-function effects, resulting in slower recovery from inactivation and decreased excitability in cultured cortical neurons. Furthermore, for the same intracellular loop of NaV1.7, we found an excess of rare variants in a case-control variant-burden study. Functional analysis of one of these variants, M932L/V991L, also demonstrated reduced firing in cortical neurons. However, although this variant is rare in Caucasians, it is frequent in Latino population, suggesting that genetic background can alter its effects on phenotype. Although the involvement of the SCN1A and SCN2A genes encoding NaV1.1 and NaV1.2 channels in de novo ASD has previously been demonstrated, our study indicates the involvement of inherited SCN9A variants and partial loss-of-function of NaV1.7 channels in the etiology of rare familial ASD.

Introduction

Autism spectrum disorder (ASD) is a behaviorally defined, childhood-onset complex neurodevelopmental disorder with a broad spectrum of symptoms and a wide range of severity. The core symptoms of ASD include deficits in social communication and social interactions, along with restrictive, repetitive patterns of behaviors, interests or activities (DSM5).[1] A recent report by the Center for Disease Control and Prevention (CDC) puts the incidence of autism as 1 in 68 children, making it a growing public health concern. Autism has a strong genetic component that is further complicated by overlap with other neurodevelopmental disorders such as intellectual disability, epilepsy and schizophrenia.[2, 3] Advances in array and sequencing technologies facilitated elucidation of the contribution of de novo copy number variations and single-nucleotide mutations to the molecular basis of sporadic ASD.[4, 5] From recent copy number variations and sequencing studies, it is estimated that there are hundreds to thousands of ASD risk loci[3, 4] and de novo events account for 10–40% of cases.[6] The analysis of genes with de novo variants using systems biology approaches such as pathway enrichment and gene co-expression profiling implicates neuronal signaling, synaptic transmission, chromatin biology and transcription regulation[4, 5, 7] in the etiology of de novo ASD. However, in familial ASD much of heritability remains unexplained.[8] Traditional linkage studies have implicated regions cumulatively covering more than 10% of the genome and only a few loci have been replicated at genome-wide significance level.[9] To identify genes involved in familial autism, new sequencing methods were initially applied in consanguineous affected individuals in order to detect recessive genes.[10, 11] More recently whole-exome sequencing has been used in multiplex families to identify heterozygous-shared variants that are associated with ASD risk.[12, 13] Such studies identified additional candidate genes for autism but lacked statistical power and did not include replication.[14] Although evidence is lacking for the involvement of specific genes, some prior studies indicated an association between voltage-gated sodium channels and familial ASD.[15]

The alpha subunits of voltage-gated sodium (NaV) channels are evolutionarily conserved and responsible for initiation and propagation of action potentials in nerve and muscle.[16] Nine pore-forming alpha subunits of NaV channels have been identified in mammals, with different channel subtypes showing specific expression profiles in various excitable tissues. Mutations in the primary central nervous system NaV channels (encoded by SCN1A, SCN2A, SCN3A and SCN8A) are associated with neurological, psychiatric and neurodevelopmental disorders including epilepsy, autism and cognitive impairment.[17, 18, 19] In the peripheral nervous system, NaV1.7 channels encoded by SCN9A are important for the excitability of sensory neurons and pain perception. Gain-of-function mutations in NaV1.7 lead to primary erythermalgia,[20] paroxysmal extreme pain disorder,[21] small nerve fiber neuropathy[22] and increased pain perception,[23] whereas homozygous complete loss-of-function (LOF) mutations cause congenital insensitivity to pain (CIP).[24, 25, 26, 27] The absence of other neurological symptoms in these sensory disorders suggested that NaV1.7 channels are crucial for relaying sensory information in humans, but their roles outside the nociceptive system are not well understood.

The traditional family-based design for gene discovery uses large extended families with multiple-affected family members and combines linkage and sequencing analyses to limit the search space for variant identification and to obtain statistical evidence for causality. For study of genetically heterogeneous traits, use of extended families in this context is particularly useful.[28] As large multigenerational families with autism are rare, we applied a modified approach. We focused on families with affected distant relatives such as second and third cousins and used stringent frequency filtering of the exome data to identify putatively causal variants under the assumption that for each family the causal variant is private and sufficient to cause the phenotype. Our focus on protein-altering variants further narrowed the number of candidate genes in each family. This approach resulted in the identification of a small number of candidate genes in each family including SCN9A. Functional analysis of the SCN9A variant, NaV1.7C1143F, revealed partial LOF effects on excitability of cultured neurons. Further case–control analysis provided additional evidence for association of NaV1.7 mutations in the second intracellular loop with ASD, and functional analysis of one of these variants, NaV1.7M932L/V991L, also demonstrated reduced firing of cortical neurons. Together, our results suggest that SCN9A is important for normal brain function and that rare variants in this gene, encoding partial LOF mutations in the second intracellular loop of NaV1.7, are involved in familial autism.

Materials and methods

Study sample

For the family-based exome sequencing study we identified three second cousin and two third cousin families in the NIMH repository (https://www.nimhgenetics.org/) and selected two affected cousins from each family for exome sequencing (N=10). All families were of Caucasian background.

For the association study, we selected 855 familial Caucasian cases and 960 Caucasian controls from the NIMH repository. An additional unrelated 288 familial Caucasian cases were included from a University of Washington (UW) multiplex autism collection and 208 Caucasian subjects from a UW study on dyslexia served as control samples. For both NIMH and UW ASD subjects, the diagnostic status was determined by gold-standard ADOS and ADI-R assessment followed by expert clinical judgment using all available information. The NIMH controls sample consists of adults who completed an online short self-report assessment to exclude severe psychiatric disorders, but were not specifically screened for autism. The UW controls were assessed for dyslexia by comprehensive testing, and screened for absence of ASD with a questionnaire. Samples from NIMH participating institutions and UW were utilized in concordance with Institutional Review Boards approvals for the participation in the NIMH repository, consenting of subjects and data sharing.

Whole-exome and targeted capture sequencing, variant identification and annotation

Exome capture using the Nimblegen SeqCap EZ Human Exome Library v2.0 (Roche, Basel, Switzerland), paired-end 50 bp sequencing on an Illumina HiSeq2000 sequencing platform (San Diego, CA, USA), and variant calling and annotation were performed in the UW Genome Sciences Center for Mendelian Genomics as described before.[29] For targeted sequencing we designed 117 single-molecule molecular inversion probes for the protein-coding region of four candidate genes and performed sequencing and analysis as described before.[30] Sequence reads were aligned to NCBI human reference genome GRCh37 (hg19) using Burrows–Wheeler Aligner[31] (BWA v0.7.10) and Genome analysis toolkit[32] (GATK v2.4) was used for variant calling. Single-nucleotide variants were annotated as nonsynonymous, splice, stop gain, stop lost or synonymous using ANNOVAR[33] and further classified as rare or private based on their frequencies in variant databases. We defined variants with frequency <0.01 in 1000 Genomes (1KG) 'European' samples as rare and variants not present in any of the data sets (dbSNP132, 1KG (April 2012) and Exome Sequencing Project (ESP) European American data sets) as private. Detailed methods are provided in the Supplementary Information.

Primary cultures of mouse cortical neurons, transient expression of NaV1.7 channels and electrophysiology

Primary cultures of cortical neurons were prepared from C57BL/6 mice on postnatal days 0–3 as described.[34] The cerebral cortices were dissected, digested, triturated and plated on coverslips. Cortical neurons were transfected with NaV1.7WT, NaV1.7C1143F and NaV1.7M932L/V991L using electroporation (Neon, Life Technologies, Carlsbad, CA, USA) or Lipofectamine 2000 (Life Technologies). Co-electroporation of cDNA encoding GFP allowed identification of the transfected cells and all the recordings were performed on GFP-labeled cells only. Recordings were made 3–5 days after transfections. Detailed methods are described in the Supplementary Information.

Quantitative real-time PCR

Cerebral cortices and hippocampi were isolated from P21-P24 mice and used for quantitative real-time PCR (comparative CT method) using TaqMan primers and probes (Applied Biosystems, Foster City, CA, USA). See Supplementary Information for full description.

Genetic association analysis

For exome-sequenced cases, the level of relatedness was confirmed using KING.[35] We used PLINK[36] to identify potentially duplicated samples in case–control analysis as some UW cases were also present in the NIMH sample. Pairs with PI_HAT (Proportion IBD) more than 0.9 were considered identical/duplicates. For such pairs only one sample was considered for further analysis. For unrelated subjects we performed gene-based variant-burden and single-variant association analyses with a X2-test (1 df, Yates correction) and calculated a P-value for each gene/region. For the association study we included protein-altering variations (nonsynonymous, splice, stop gain and stop loss) with a frequency <0.01 in 1KG EA. For the SCN9A LII–III region that was associated with autism in our case–control study we performed transmission disequilibrium test (TDT) analysis to control for spurious association that might be due to population substructure. TDT uses X2 goodness of fit statistics to determine if there is a preferential transmission of a risk allele to the affected cases. In the TDT analysis we used affected siblings of cases with rare protein-altering variants in the SCN9A LII–III region, under the assumption that if SCN9A LII–III variants contribute to the phenotype we should see overtransmission of risk variants to affected siblings. This analysis incorporated a rare-variant extension of the TDT as described.[37] In a manner equivalent to our case–control study we used a TDT-burden of rare variants that counts the number of minor-allele-transmission events to affected siblings.

Variant validation and evaluation in family members

Selected private variants identified by exome sequencing that were shared between affected cousins were independently validated using Sanger sequencing and tested for co-segregation in additional available family members. Rare variants identified in case–control-targeted sequencing that code for the second intracellular loop of SCN9A were also confirmed with Sanger sequencing. For each proband, available family members were genotyped to determine if they carry the variant allele. Primers for PCR amplification and sequencing were designed using Primer3 (v 4.0.0) (http://bioinfo.ut.ee/primer3-0.4.0/primer3/) and bidirectional sequencing on an ABI 3730 DNA analyzer (Applied Biosystems) was done as previously described.[13]

Bioinformatics analysis of gene expression, co-expression and biological function enrichment

Expression data based on RNA sequencing of 53 different tissues, including 13 brain regions, were downloaded from the GTEx portal.[38] To identify genes with a brain expression pattern that correlates with that of SCN9A, we calculated the Pearson’s correlation coefficient for each gene and SCN9A. Genes with strong correlation (−0.5⩾r⩾0.5) with P<0.05 were retained. Pearson correlations and P-values were calculated using the rcorr function of library Hmisc in R. To identify biological implications for this set of genes we performed Gene Ontology Biological Processes enrichment analysis using Database for Annotation, Visualizations and Integration Discovery (DAVID).[39] For DAVID analysis Benjamini–Hochberg adjusted P<0.01 and fold enrichment >2 were considered to be significant.

Genotype–phenotype correlation analysis

To investigate whether SCN9A LII–III variants associate with a phenotypic subtype of ASD we compared available phenotypes in 29 subjects who carried LII–III variants to a comparison sample of 2334 individuals from the full NIMH sample. We examined demographics, cognitive abilities (nonverbal intelligence quotient and receptive vocabulary), parental report of autism symptoms domains as assessed with the Autism Diagnostic Interview,[40] milestones (age of walking, use of first words and phrases), as well as specific clinical characteristics assessed with parental questionnaire (regression, aggression, self injury and history of seizures). To detect differences between the groups, analysis of variance was used for the quantitative variables age, cognition and parental report of autism symptom domains, whereas X2 nonparametric tests were used for gender and rates of specific phenotypic characteristics (for example, regression). All analyses were performed using IBM SPSS version 23.0 (Armonk, NY, USA).

Results

Family exome analysis

We sequenced 10 exomes from five families with cousins affected with autism, which are deposited in NDAR (https://ndar.nih.gov/edit_collection.html?id=1919). On average, we generated 3.8 Gb of mapped sequence data per individual and >92% of bases had >8x coverage across all samples. After quality filtering, on average 23 685 +/− (429) variants were identified in each individual. After frequency and function filtering was performed, each subject on average had 104 +/− (26) gene disrupting (splice, frameshift and stop) and nonsynonymous private variants (Supplementary Table 1). The affected cases in each family shared between one to four private variants that were all nonsynonymous. We evaluated the impact of missense variants on the structure and function of the proteins using physical and comparative considerations as implemented in SIFT, Polyhen2 and GERP to identify variants more likely to have effects on the protein function (Table 1). Variants were confirmed and evaluated if they were present in additional affected family members. In this way we identified four candidate genes SCN9A, PLEC, KANK3 and CCDC73, which we further evaluated for association with autism. Given the strong evidence for involvement of voltage-gated sodium channels in ASD, and the high impact of the SCN9AC1143F variant based on bioinformatics analysis, we performed functional studies on NaV1.7C1143F.

Table 1

Private variants shared by the affected cousins of each exome-sequenced family

Family id	Position	Gene	Nucleotide/ AA change	SIFT	Polyphen-2	GERP
74-0668	chr14:104029378	APOPT1	c.C79A:p.P27T	0.01	0	−1.81
	chr19:35850022	FFAR3	c.A230G:p.N77S	0.88	0.003	−6.53
152MM0304	chr1:3328298	PRDM16	c.G1537A:p.G513S	0.75	0.852	2.22
	chr2:167108286	SCN9A	c.G3428T:p.C1143F	0	1	5.82
156-3860	chr6:167594182	TCP10L2	c.A831T:p.E277D	1	0	−1.58
	chr8:144991583	PLEC	c.G12364C:p.V4122L	0.99	1	5.08
	chr11:32636075	CCDC73	c.G1789C:p.E597Q	0.41	0.047	1.63
	chr19:55998317	NAT14	c.C615G:p.D205E	0.1	0.009	0.636
156-3897	chr19:8399390	KANK3	c.C1321G:p.P441A	0.23	0.131	3.36
	chr19:58117083	ZNF530	c.A190T:p.T64S	0.79	0.001	−3.39
152MM0122	chr2:233712266	GIGYF2	c.G3651C:p.Q1217H	0.06	0.22	−2.51

Abbreviations: AA, amino acid; GERP, genomic evolutionary rate profiling; PolyPhen-2, polymorphism phenotyping v2; SIFT, sorting intolerant from tolerant.

SIFT predicts impact of amino acid substitutions based on the degree of conservation in sequence alignments derived from closely related sequences. Scores <0.05 are considered deleterious. PolyPhen-2 predicts impact of a variant on the structure and function of a human protein using eight sequence-based and three structure-based predictive features. Scores >0.95 are considered probably damaging. GERP identifies functional constraint of a sequence variant by quantifying substitution deficits in multiple alignments. Substitution deficits represent a natural measure of constraint that reflects the strength of past purifying selection. Higher GERP scores are more deleterious.

Genes selected for case–control study are noted in bold.

Functional analysis of the NaV1.7C1143F variant

To investigate the biophysical ramifications of this mutation, we first expressed NaV1.7WT and NaV1.7C1143F in human embryonic kidney tsA-201 cells and examined the biophysical properties of the sodium currents. When NaV1.7 channels were expressed without auxiliary β subunits, the voltage dependence and kinetics of activation and fast inactivation were unchanged by the mutation (Supplementary Figures 1a–d), but we observed a trend toward more negative voltage dependence of slow inactivation (Supplementary Figure 1e). Activation, fast inactivation and slow inactivation were unaffected by the mutation when NaV1.7 channels were co-expressed with β1 and β2 subunits (Supplementary Figures 1f–j).

In neurons, NaV1.7 channels have many additional interacting partners including multiple combinations of β subunits,[41] fibroblast growth factor homology factors,[42] protein kinases[43] and potentially other regulatory proteins. Therefore, further characterization of the NaV1.7C1143F variant was carried out in cultured cortical neurons from mice. We confirmed that SCN9A is expressed in frontal cortex and hippocampus of mice using quantitative PCR. SCN9A expression is clearly detectable although 52.6±5.5-fold lower compared with the expression of SCN1A, a gene with a well-demonstrated role and expression in central nervous system (Supplementary Figure 2). Our results of SCN9A expression are in full agreement with recently published results of mouse cortex and hippocampus single-cell RNA sequencing studies,[44] as well as with previous results from monkeys.[45] In addition evidence for SCN9A brain expression in mice has been documented, although not published (Mouse Genome Database (MGD), MGI: 3631994, (URL: http://www.informatics.jax.org/) (April 2016)).

NaV1.7 currents recorded in voltage clamp mode with endogenous sodium channels blocked revealed that the NaV1.7C1143F mutation had no effect on peak current amplitude, voltage dependence of activation or voltage dependence of fast inactivation (Figures 1a–c). However, the rate of recovery from fast inactivation was two-fold slower in NaV1.7C1143F (Figure 1d). Sodium channel availability, during repetitive depolarization at 20 Hz and 50 Hz, declined ~20% in amplitude at both stimulation frequencies for NaV1.7WT. Notably, for NaV1.7C1143F, the decline was significantly increased, with a 35% reduction in current amplitude at 20 Hz and a 67% reduction at 50 Hz (Figures 1e and f).

Figure 1

Properties of NaV1.7WT and NaV1.7C1143F expressed in cultured cortical neurons. (a–f) Biophysical characterization of NaV1.7WT and NaV1.7C1143F in cultured cortical neurons. Because cortical neurons express endogenous NaV1.7 channels, as well as many other types of sodium channels, we used a tetrodotoxin (TTX)-resistant form of NaV1.7 (Supplementary Information, NaV1.7Y362S) and blocked endogenous channels with 500 nm TTX. (a) Representative set of sodium current traces from cortical neurons expressing NaV1.7WT and NaV1.7C1143F. (b) Mean current–voltage (I–V) relationships of peak currents. (c) Voltage dependence of activation (right curves, the V1/2 of channel activation was 19.7±4.4 mV for NaV1.7WT and 19.3±9.3 mV for NaV1.7C1143F, P>0.05) or the voltage dependence for steady-state fast inactivation (left curves, V1/2 for NaV1.7WT was −62.6±1.6 mV and for NaV1.7C1143F −69.7±3.62 mV, P>0.05). (d) Recovery from fast inactivation (τ=3.72±0.6 ms in NaV1.7WT and 7.63±1.35 ms in NaV1.7C1143F, P=0.017). (e and f) Mean normalized currents during 100 depolarizations to 0 mV at 20 Hz (e) and 50 Hz (f). NaV1.7WT, n=7, NaV1.7C1143F, n=5. Firing of cortical neurons expressing NaV1.7WT and NaV1.7C1143F. (g) Representative firing in response to 160 pA depolarizing current injection. (h) Membrane input resistance. (i) Rheobase, the minimal current needed to evoke action potentials (AP) (j) average number of AP in response to 1 s depolarizing current injection at the indicated intensity. The resting membrane potential was unaffected by the C1143F mutation (−65.8±2.7 mV in NaV1.7WT and −60.7±3.6 mV in NaV1.7C1143F, P>0.05). Similarly, the threshold for action potential (−32.7±1.7 mV for NaV1.7WT and −27.68±5.4 mV for NaV1.7C1143F, P>0.05) and cell capacitance (95.6±22 pF in NaV1.7WT and 85.3±37 pF in NaV1.7C1143F, P>0.05) were unchanged. *P<0.05. NaV1.7WT, n=7, NaV1.7C1143F, n=5.

The measures of cortical neuron membrane capacitance, resting membrane potential and threshold for action potential generation were comparable between neurons expressing NaV1.7WT and NaV1.7C1143F (Figure 1g). For neurons expressing NaV1.7C1143F, the input resistance was lower and rheobase increased, requiring larger current injection to evoke action potentials (Figures 1h and i; P<0.05). In neurons expressing NaV1.7C1143F the number of action potentials was decreased during trains evoked by a broad range of injected current intensities (Figure 1j).

Genetic association analysis and functional analysis of the NaV1.7M932L/V991L variant

To further test the significance of mutations in SCN9A in autism, we conducted extensive genetic association analyses. After removal of duplicated and poorly captured samples, 1004 unrelated familial cases (NIMH cases=800, UW cases=204) and 1127 controls (NIMH controls=924, UW controls=203) remained. Gene-based variant-burden analysis did not identify significant results for any of the four evaluated genes (Table 2). However, single-variant analysis was significant for the SCN9AV991L rs4369876 variant (P=0.00004) that is located in the same cytoplasmic loop (LII–III) as the C1143F variant that we observed in a single ASD exome-sequenced family. The V991L variant was previously reported to be in complete linkage disequilibrium with an M932L variant.[46] In our analysis the M932L variant was initially filtered out due to poor quality of reads for that position. Capillary sequencing confirmed the presence of M932L in all subjects with the V991L variant. Eight additional rare protein-altering variants were observed in LII–III and none in controls (P=0.0083). For LII–III, we identified a total of 29 rare variants in cases and 2 in controls for a combined variant burden LII–III, P=5.1 × 10−7 (Table 3). Variants identified in all four genes are reported in Supplementary Table 2.

Table 2

Gene-based variant-burden association study for SCN9A, PLEC, KANK3 and CCDC73 in 1004 familial cases and 1127 unscreened controls

Gene	Case	Control	X2-test (Yates’ P-value)
SCN9A	67/936	60/1067	0.18
SCN9A (LII–III)	29/975	2/1119	5.1 × 10−7
KANK3	5/999	8/1119	0.72
CCDC73	1/1003	5/1122	0.27
PLEC	7/997	5/1122	0.62

Abbreviation: ASD, autism spectrum disorder.

For PLEC only 1 kb coding region surrounding original variant identified in affected family was sequenced. Association for the second intracellular loop (LII–III) of SCN9A gene is shown in bold. Case: number of ASD case samples with variant/number of ASD case samples without variant. Control: number of control samples with variant/number of control samples without variant.

Table 3

All variants identified in second intracellular loop (LII–III) of SCN9A

Position	Nucleotide change	AA change	ESP6500 EU	1000Genome EUR	dbSNP132	Case	Control	SIFT	Polyphen-2	GERP
chr2:167099083	c.T3523A	p.Y1175N	0	0	NA	1/1003	0/1127	0	0.969	5.01
chr2:167099157	c.G3449A	p.W1150X	0	0	NA	1/1003	0/1126	0	0	1.74
chr2:167108345	c.G3369T	p.L1123F	0.000364	0.0013	NA	3/1001	0/1126	0.05	0.13	−2.18
chr2:167108386	c.C3328T	p.R1110W	0.000854	0.0026	NA	1/901	0/1091	0.01	0.6	3.92
chr2:167129091	c.G3136A	p.D1046N	0.000121	0	NA	1/1003	0/1126	1	0.003	5.42
chr2:167129135	c.C3092T	p.T1031I	0	0	NA	1/1001	0/1110	0.23	0.054	5.42
chr2:167129256	c.G2971T	p.V991L	0.004177	0	rs4369876	21/976	2/1119	0.08	0.004	−0.949
chr2:167133540a	c.A2794C	p.M932L	0.004070	0	rs12478318	21/976	2/1119	0.53	0.043	5.65

Abbreviations: AA, amino acid; ASD, autism spectrum disorder; dbSNP, database of short genetic variation; ESP, Exome Sequencing Project; GERP, genomic evolutionary rate profiling; NA, variant is not present; PolyPhen-2, polymorphism phenotyping v2; SIFT, sorting intolerant from tolerant.

ESP and 1000Genome frequencies are only for 'European American' and 'European' population, respectively. Case: number of ASD case samples with variant/number of ASD case samples without variant. Control: number of control samples with variant/number of control samples without variant. SIFT, Polyphen-2 and GERP are described in Table 1.

This variant was initially filtered out due to poor quality of reads for that position. Capillary sequencing confirmed the presence of M932L in all subjects with the V991L variant.

The TDT was performed on 32 affected siblings of cases that had rare gene disrupting variants in the SCN9A LII–III region. For 25 affected siblings, the proband had the SCN9AM932L/V991L variant, whereas for 7 siblings, the proband had the other rare LII–III variants. We found overtransmission of rare alleles to affected siblings for SCN9AM932L/V991L variant (17 minor alleles, 8 major alleles, X2=3.24, P=0.071) and for other LII–III variants (5 minor alleles, 2 major alleles, X2=1.28, P=0.256). Combining SCN9AM932L/V991L and other rare LII–III variants the overtransmission was statistically significant (22 minor alleles, 10 major alleles, X2=4.5, P=0.033).

As NaV1.7M932L/V991L variant, which was present in 21 cases and 2 controls, carried the bulk of genetic association, we performed functional analysis of the NaV1.7M932L/V991L variant in cultured cortical neurons (Figure 2). NaV1.7M932L/V991L variant did not affect peak current amplitude, voltage dependence of activation, voltage dependence of fast inactivation or channel availability during repetitive depolarization at 50 Hz (Figures 2a–d). However, neurons expressing NaV1.7M932L/V991L fired fewer action potentials during 1 s long depolarizing current injection, indicating reduced excitability (Figures 2e and f).

Figure 2

Properties of NaV1.7WT and NaV1.7M932L/V991L expressed in cultured cortical neurons. (a) Representative set of sodium current traces from cortical neurons expressing NaV1.7WT and NaV1.7M932L/V991L. (b) Mean current–voltage (I–V) relationships of peak currents. (c) Voltage dependence of activation (right curves) and the voltage dependence of steady-state fast inactivation. (d) Mean normalized currents during 100 depolarizations to 0 mV at 50 Hz. NaV1.7WT, n=4, NaV1.7M932L/V991L, n=5. (e and f) Firing of cortical neurons expressing NaV1.7WT and NaV1.7M932L/V991L. (e) Representative firing in response to 100 pA depolarizing current injection. (f) Average number of action potentials (AP) in response to 1 s depolarizing current injection at the indicated intensity. NaV1.7WT, n=7, NaV1.7M932L/V991L, n=7.

Phenotype analysis

Available phenotypic information for the exome-sequenced families is presented in Supplementary Information. Genotype–phenotype correlation analysis for 29 subjects that had variants in SCN9A LII–III did not reveal any specific phenotype characteristics for available measures. Phenotype information is shown in Supplementary Table 3.

Bioinformatics analysis of SCN9A expression, gene co-expression and biological function enrichment

In humans, the GTeX data[38] also confirms the expression of SCN9A in the brain (Supplementary Figure 3). Brain co-expression analysis of GTeX data identified 538 genes positively co-expressed with SCN9A. DAVID analysis identified seven enriched biological processes for SCN9A co-expressed genes: forebrain development and neuronal differentiation, regulation of secretion and hormone levels, cell–cell signaling, and feeding and adult behavior (Supplementary Table 4). The most significant enrichment was observed for feeding behaviors (P=8.8 × 10−6, 7.91-fold enrichment, 15 genes) and forebrain differentiation (P=0.002, 4.06-fold enrichment, 17 genes).

Discussion

Although recent studies have identified multiple genetic factors involved in sporadic, de novo autism, genetic factors contributing to familial autism are less well characterized. Here, to identify genes that harbor rare variants that contribute to familial autism, we adopted a three-step design that combines family-based whole-exome sequencing for candidate gene identification, followed by electrophysiology studies to show functional relevance of identified variants, and additional genetic analysis to provide statistical evidence for the association with ASD. For the exome family analysis, we narrowed the number of identified genes by focusing on variants that are both private and protein disrupting. In this way we identified 11 candidate variants in 5 families. Bioinformatics analysis suggested that four of these variants were more likely to have functional effects. A variant-burden analysis of four genes was significant for rare protein disrupting variants in the second intracellular loop of SCN9A. This is the same LII–III loop that harbored a SCN9AC1143F variant in one family that was exome-sequenced. In further support of the association with autism, for LII–III variants we show preferential transmission to affected siblings as well. SCN9A is in tight linkage disequilibrium with a cluster of sodium channel genes on chromosome 2q24.3 (SCN1A, SCN2A, SCN3A and SCN7A) and SCN1A is <200 kb away. This raises the possibility that the association we have identified with SCN9A might be due to mutations in SCN1A and SCN2A genes that are implicated in autism as well. Although we were not able to directly sequence SCN1A and SCN2A genes, we feel that this scenario is very unlikely as the main phenotype associated with SCN1A (OMIM* 182389) and SCN2A (OMIM *182390) mutations are seizures, which are not reported in our subjects carrying SCN9A variants. However we cannot exclude the possibility of SCN1A involvement, as mutations in this gene were associated with autism without seizures as well.[47, 48] In addition to genetic association, both NaV1.7C1143F variant identified in a single exome-sequenced family and NaV1.7M932L/V991L variant identified in association analysis demonstrated a partial LOF effect when expressed in cortical neurons.

We did not find genetic evidence for association in the four other families studied by exome sequencing. This is not surprising as ASD is highly heterogeneous. In these families other genetic mechanisms, that we have not evaluated, might be contributory. Such mechanisms include disruptions of noncoding regulatory DNA,[49] genetic heterogeneity within the family[50] and polygenic inheritance.[51]

Functional consequences of the NaV1.7C1143F and NaV1.7M932L/V991L variants

Electrophysiological analysis of NaV1.7C1143F expressed in tsA-201 human embryonic kidney cell lines with and without auxiliary β1 and β2 subunits suggests, at most, a small effect on basic sodium channel functions. In neurons, NaV1.7 channels have many additional interacting partners, and indeed, in cortical neurons the rate of recovery from fast inactivation was two-fold slower for NaV1.7C1143F channels. Such slowed recovery is expected to result in longer periods of reduced sodium channel availability during and following trains of stimuli. In agreement with this expectation, sodium current amplitude conducted by NaV1.7C1143F expressing neurons declined significantly more compared with NaV1.7WT during repetitive depolarization at 20 and 50 Hz. These biophysical alterations are both expected to result in reduced firing and excitability of cortical neurons. The experiments in current-clamp mode revealed that neurons expressing NaV1.7C1143F had lower input resistance and increased rheobase, so larger current injection was needed to evoke action potentials, resulting in decreased action potential firing (Figure 1).

NaV1.7M932L/V991L variant has been implicated before in neuropathic pain syndromes[22] and is associated with partial deletion of pain perception,[52] suggesting both increased and decreased neuronal firing. Here, electrophysiological analysis of NaV1.7M932L/V991L did not detect changes in the biophysical properties of the currents. However, similarly to NaV1.7C1143F, cortical neurons expressing NaV1.7M932L/V991L fired fewer action potentials, indicating reduced excitability (Figure 2). Although we do not completely understand the causes for reduced firing in the NaV1.7M932L/V991L variant, it might be related to changes in cellular localization or neuronal maturation. LII–III contains binding sites for ankyrin-G that are crucial for the localization of sodium channels at the axon initial segment.[43] Neuronal excitability was therefore examined in mature neurons (after the formation of axon and dendrites). In contrast, biophysical characterization was performed on immature cultures (3–4 days in vitro) in order to avoid space clamp problems, but potentially before the manifestation of functional alterations. In addition, unlike the C1143F variant, which was found to be highly conserved and predicted to be damaging with three bioinformatics tools (Table 1), strong evidence for conservation and potential impact on protein function is lacking for the M932L and V991L variants (Table 3). Furthermore, functional implication of M932L/V991L variant might strongly depend on genetic background, as indicated by its frequency in different populations. On the basis of ExAC database that aggregates multiple studies (www.ExAC.broadinstitute) M932L and V991L variants are rare in European (MAF=0.00236) and African American (MAF=0.0029) populations with frequencies that are comparable to our Caucasian controls (MAF=0.0016). However, M932L and V991L variants have higher frequency in East (MAF=0.061) and South Asians (MAF=0.012), and in Latino population variant frequency (MAF=0.23) far exceeds the frequency we have observed in our Caucasian cases (MAF=0.022). However, for sodium channel genes, genetic background is known to play a role in seizure phenotype expression[53] and in mouse model SCN1A mutations resulting in severe epilepsy and impaired cognitive abilities in the C57BL/6 strain and in extremely rare spontaneous seizures with normal cognitive function on the 129/SvJ background.[54]

Nevertheless, although the evidence for the disease-causing role of NaV1.7M932L/V991L variant is tempered due to high population frequency of the variant in Asian and Latino populations and lack of strong evolutionary conservation, we have identified six additional rare LII–III variants that support association with ASD. Overall our findings indicate that ASD-associated SCN9A variants result in partial LOF of NaV1.7 channel, causing decreased action potential firing and altering the input–output relationships in circuits whose neurons express SCN9A.

How might NaV1.7 channels contribute to autism?

Although traditionally considered as a gene that is important for the peripheral nervous system, there is increased understanding that SCN9A has a role in the central nervous system as well. Missense mutations in SCN9A are associated with familial febrile seizures,[55] a benign form of epilepsy. Similarly, mutations in SCN9A can modify the severity of Dravet syndrome, a rare genetic intractable epilepsy syndrome caused by complete LOF mutations in NaV1.1 channels, which are prominent in interneurons of the brain.[55, 56] One emerging concept from studies of mouse models of ASDs, including idiopathic autism,[57] Rett Syndrome,[58] Fragile X Syndrome[59] and Dravet Syndrome,[60] is that core autistic-like behaviors are associated with an increased ratio of excitatory to inhibitory neurotransmission in the brain. Most of the experimental evidence suggests increased excitation. In contrast, in Dravet syndrome mice, autistic-like behaviors are caused by LOF mutations in NaV1.1 channels and consequent selective impairment of action potential firing in GABAergic inhibitory neurons.[60, 61] These autistic-like behaviors are rescued by enhancement of inhibitory neurotransmission by treatment with low doses of clonazepam, a positive allosteric modulator of GABAA receptors.[60] Data presented here (Supplementary Figure 2) confirm the expression of NaV1.7 in the cortex. The single-cell neuronal gene transcription analysis by Zeisel et al.[44] demonstrated that the NaV1.7 channel is expressed in both GABAergic and glutamatergic neurons, but with higher levels of expression in interneurons.[44] Thus, by analogy to LOF mutations in NaV1.1, a hypothesis can be made that rare variants in NaV1.7 decrease the firing of a specific set of GABAergic neurons that are important in control of social behaviors. Furthermore, using GTeX data we show that SCN9A is co-expressed with genes that are enriched for forebrain neuronal differentiation. These genes include transcription factors DLX1 and DLX2 that are essential for the production of forebrain GABAergic interneurons during embryonic development[62] and have been implicated in autism in a genetic association study.[63] These observations support a hypothesis that functional variants in SCN9A channels contribute to autism phenotype by affecting a specific set of GABAergic neurons that are important in control of social behaviors.

Our analyses of GTeX data indicate that SCN9A is highly expressed in the hypothalamus, a brain region that has been implicated in autism.[64, 65] Among other roles, the hypothalamus is responsible for synthesizing the behaviorally important hormone oxytocin that has a role in social-bonding behaviors and has been implicated in autism as well.[66] In addition, the hypothalamus is a crucial regulator of feeding behavior.[67] We have found that genes co-expressed with SCN9A are enriched for functions related to feeding behaviors and the function of NaV1.7 channels in hypothalamic neurons was shown to be important for body weight regulation.[68] This indicates that, in addition to the role in GABAergic neurotransmission, the role of SCN9A in autism might be mediated through changes in hypothalamic functions, which in turn can affect multiple hormonally regulated processes that are frequently disrupted in autism such as oxytocin-mediated social interactions[69] and feeding behaviors.[70]

Further directions

Until now, truncating mutations in NaV1.7 with full LOF effects that completely abolish NaV1.7 currents were found to be responsible for rare autosomal-recessive CIP.[24, 25, 26, 27] Fewer than 100 individuals with CIP were reported in the literature and phenotype description in these reports indicate that beyond CIP other sensory modalities are unaffected. Interestingly, extensive neuropsychological testing was reported for a single 9-year-old girl with partial CIP, caused by compound heterozygous partial LOF SCN9A variant. On the basis of parental report this individual had impaired social skills and below average levels of empathy.[71] Although in this report presence of ASD was not evaluated, presence of deficits in social skills and empathy indicates the need for more detailed evaluation for ASD in individuals with CIP. In addition, our study has shown the association of heterozygous partial LOF variants with ASD. This indicates that carriers of SCN9A CIP variants might be at increased risk for autism as well, and formal ASD screening and evaluation in patients with CIP will allow to establish the diagnosis, which in turn, will improve the care and treatments such individuals are receiving.

In our sample, we were unable to detect any specific autism subtype based on demographic, ADOS variables and clinical characteristics. However, our analysis of SCN9A brain co-expressed genes indicates that feeding behaviors might be a prominent component of the SCN9A-associated phenotype. We did not have information about feeding behaviors for our cohort. Our results indicate that ASD subjects with deregulated feeding behaviors might be enriched for SCN9A variants. In addition dysregulated-feeding behaviors could be a phenotype that is suitable for genetic studies and it might be valuable to add such phenotype to genetic studies of autism.

The physiological role of LII–III region in voltage-gated sodium channels is not well understood. Earlier human genetic studies implicates variants in LII–III in a growing number of hereditary disorders including cardiac arrhythmias[72, 73, 74] and epilepsy;[56] and we now add the association with ASD. Reported functional analysis of these variants demonstrated effects on channel inactivation,[72, 73, 74, 75] and we report alterations in the rate of recovery from fast inactivation and neuronal excitability (Figures 1–2). Further functional studies on SCN9A and other brain expressed NaV family genes that examine electrophysiological properties and subcellular localization in specific cell subtypes like GABAergic neurons are warranted to further elucidate the mechanism by which NaV channels LII–III variants increase the risk for ASD. In addition, pharmacological interventions in SCN1A-null mice that reverse ASD-like symptoms[60] indicate that for individuals with ASD and NaV mutations ASD, symptoms might be responsive to pharmacological interventions as well.