Synaptic, transcriptional, and chromatin genes disrupted in autism

Summary

Features of subjects with autism spectrum disorder (ASD) include compromised social communication and interaction. Because the bulk of risk arises from de novo and inherited genetic variation^1-10, characterizing which genes are involved informs on ASD neurobiology and on what makes us social beings.

Whole-exome sequencing (WES) studies have proved fruitful in uncovering risk-conferring variation, especially by enumerating de novo variation, which is sufficiently rare that recurrent mutations in a gene provide strong causal evidence. De novo loss-of-function (LoF) single-nucleotide variants (SNV) or insertion/deletion (indel) variants^11-15 are found in 6.7% more ASD subjects than in matched controls and implicate nine genes from the first 1000 ASD subjects^11-16. Moreover, because there are hundreds of genes involved in ASD risk, ongoing WES studies should identify additional ASD genes as an almost linear function of increasing sample size¹¹.

Here, we conduct the largest ASD WES study to date, analyzing 16 sample sets comprising 15,480 DNA samples (Supplementary Table 1; Extended Data Fig. 1). Unlike earlier WES studies, we do not rely solely on counting de novo LoF variants, rather we use novel statistical methods to assess association for autosomal genes by integrating de novo, inherited and case-control LoF counts, as well as de novo missense variants predicted to be damaging. For many samples original data from sequencing performed on Illumina HiSeq 2000 systems were used to call SNVs and indels in a single large batch using GATK (v2.6). De novo mutations were called using enhancements of earlier methods¹⁴ (Supplementary Information), with calls validating at extremely high rates.

After evaluation of data quality, high-quality alternate alleles with a frequency of < 0.1% were identified, restricting to LoF (frameshifts, stop gains, donors/acceptor splice site mutations) or probably damaging missense (Mis3) variants (defined by PolyPhen-2¹⁷). Variants were classified by type (de novo, case, control, transmitted, non-transmitted) and severity (LoF, Mis3), and counts tallied for each gene.

Some 13.8% of the 2270 autism trios (two parents and one affected child) carried a de novo LoF mutation – significantly in excess of expectation¹⁸ (8.6%, P<10⁻¹⁴) or what is observed in 510 control trios (7.1%, P=1.6×10⁻⁵) collected here and previously published¹⁵. Eighteen genes (Table 1) were hit by 2 or more de novo LoF mutations. These genes are all known or strong candidate ASD genes, but given the number of trios sequenced, we expect approximately two such genes by chance given gene mutability^14,18. While we expect only 2 de novo Mis3 events in these 18 genes, we observe 16 (P=9.2×10⁻¹¹, Poisson test). Because much of our data exist in cases and controls and because we observed an additional excess of transmitted LoF events in the 18 genes, it is evident that the optimal analysis framework must involve an integration of de novo mutation with variants observed in cases and controls and transmitted or untransmitted from carrier parents. Going beyond de novo LoFs is also critical given that many ASD risk genes and loci have mutations that are not completely penetrant.

Transmission and De novo Association

We adopted TADA (for ‘Transmission and De novo Association’), a weighted, statistical model integrating de novo, transmitted and case-control variation¹⁹. TADA uses a Bayesian gene-based likelihood model including per gene mutation rates, allele frequencies, and relative risks of particular classes of sequence changes. We modeled both LoF and Mis3 sequence variants. Because no aggregate association signal was detected for inherited Mis3 variants, they were not included in the analysis. For each gene, variants of each class were assigned the same effect on relative risk. Using a prior probability distribution of relative risk across genes for each class of variants, the model effectively weighted different classes of variants in this order: de novo LoF > de novo Mis3 > transmitted LoF, and allowed for a distribution of relative risks across genes for each class. The strength of association was assimilated across classes to produce a gene-level Bayes Factor (BF) with a corresponding False Discovery Rate or FDR q-value. This framework increases the power compared to use of de novo LoF alone (Extended Data Fig. 2).

TADA identified 33 autosomal genes with an FDR < 0.1 (Table 1) and 107 genes with an FDR < 0.3 (Supplementary Tables 2 and 3 and Extended Data Fig. 3). Of the 33 genes, 15 (45.5%) are known ASD risk genes⁹; 11 have been reported previously with mutations in ASD patients but were not classed as true risk genes owing to insufficient evidence (SUV420H1^11,15, ADNP¹², BCL11A¹⁵, and CACNA2D3^15,20, CTTNBP2¹⁵, GABRB3²⁰, CDC42BPB¹³, APH1A¹⁴, NR3C2¹⁵, SETD5^14,21, TRIO¹¹); and 7 are completely novel (ASH1L, MLL3, ETFB, NAA15, MYO9B, MIB1, VIL1). ADNP mutations have recently been identified in 10 patients with ASD and other shared clinical features²². Two of the newly discovered genes, ASH1L and MLL3, converge on chromatin remodeling. MYO9B plays a key role in dendritic arborization²³. MIB1 encodes an E3 ubiquitin ligase critical for neurogenesis²⁴ and is regulated by miR-137²⁵, a microRNA that regulates neuronal maturation and is implicated in risk for schizophrenia²⁶.

When the WES data from genes with FDR < 0.3 were evaluated for the presence of deletion copy number variants (such CNVs are functionally equivalent to LoF mutations), 34 CNVs meeting quality and frequency constraints (Supplementary Information) were detected in 5781 samples (Extended Data Fig. 1). Of the 33 genes with FDR < 0.1, three contained deletion CNVs mapping to three ASD subjects and one parent. Of the 74 genes meeting the criterion 0.1 ≤ FDR < 0.3, about a third could be false positives. Deletion CNVs were found in 14 of these genes and the data supported risk status for 10 of them (Extended Data Table 1, Extended Data Fig. 4). Two of the 10, NRXN1 and SHANK3, were previously implicated in ASD^2,3,10. The risk from deletion CNVs, as measured by the odds ratio, is comparable to that from LoF SNV in cases versus controls or transmission of LoF from parents to offspring.

Estimated odds ratios of top genes

Inherent in our conception of the biology of ASD is the notion that there is variation between genes in their impact on risk: for a given class of variants (e.g., LoF), some genes have large impact, others smaller, and still others have no effect at all. Yet mis-annotation of variants, among other confounds, can produce false variant calls in subjects (Supplementary Information). These confounds can often be overcome by examining the data in a manner orthogonal to gene discovery. For example, females have greatly reduced rates of ASD relative to males (a so-called ‘female protective effect’). Consequentially, and regardless of whether this is diagnostic bias or biological protection, females have a higher liability threshold, requiring a larger genetic burden before being diagnosed^21,27,28. A corollary is that if a variant has the same effect on autism liability in males as it does in females, that variant will be at higher frequency in female ASD cases compared to males. Importantly, the magnitude of the difference is proportional to risk as measured by the odds ratio (OR); hence, the effect on risk for a class of variants can be estimated from the difference in frequency between males and females.

Genes with FDR < 0.1 show profound female enrichment for de novo events (P=0.005 for LoF, P=0.004 for Mis3), consistent with de novo events having large impact on liability (OR ≥ 20; Extended Data Fig. 5). Genes with FDR between 0.1 and 0.3, however, show substantially less enrichment for female events, consistent with a modest impact for LoF variants (OR range 2-4, whether transmitted or de novo) and little to no effect from Mis3 variants. The results are consistent with inheritance patterns, LoF mutations in FDR < 0.1 genes are rarely inherited from unaffected parents while those in the 0.1 < FDR < 0.3 group are far more often inherited than de novo.

By analyzing the distribution of relative risk over inferred ASD genes¹⁹, the number of ASD risk genes can be estimated. The estimate relies on the balance of genes with multiple de novo LoF mutations versus those with only one: the larger the number of ASD genes, the greater proportion that will show only one de novo LoF. This approach yields an estimate of 1,150 ASD genes (Supplementary Information). While there are many more genes to be discovered, many will have a modest impact on risk compared to the genes in Table 1.

Enrichment analyses

FDR < 0.3 gene sets are strongly enriched for genes under evolutionary constraint¹⁸ (P=3.0×10⁻¹¹, Fig. 1a, Supplementary Table 4), consistent with the hypothesis that heterozygous LoF mutations in these genes are ASD risk factors. Indeed over 5% of ASD subjects carry de novo LoF mutations in our FDR < 0.3 list. We also observed that genes in the FDR < 0.3 list had a significant excess of de novo LoF events detected by the largest schizophrenia WES study to date²⁹ (P=0.0085, Fig. 1a), providing further evidence for overlapping risk loci between these disorders and independent confirmation of the signal in the gene sets presented here.

Open in a separate window

Figure 1

ASD genes in synaptic network

a. Enrichment of 107 TADA genes in: FMRP targets from two independent datasets and their overlap; RBFOX targets; RBFOX targets with predicted alterations in splicing; RBFOX and H3K4me3 overlapping targets; genes with de novo mutations in schizophrenia; human orthologues of Genes2Cognition mouse synaptosome or PSD genes; constrained genes; and, genes encoding mitochondrial proteins (as a control). Red bars indicate empirical P-values. b. Synaptic proteins encoded by TADA genes. c. De novo Mis3 variants in Na_v1.2 (SCN2A). The four repeats (I-IV) with P-loops, the EF-hand, and the IQ domain are shown, as are the four amino acids (DEKA) forming the inner ring of the ion selectivity filter. d. Relevant variants in Ca_v1.3 (CACNA1D). Part of the channel is shown, including helices one and six (S1 and S6) for the I-IV domains, NSCaTE motif, EF-hand domain, pre-IQ, IQ, PCRD, DCRD, proline-rich region, and PDZ-binding motif.

We found significant enrichment for genes encoding mRNAs targeted by two neuronal RNA-binding proteins: FMRP³⁰ (also known as FMR1), mutated or absent in fragile X syndrome (P=1.20×10⁻¹⁷, 34 targets³⁰, of which 11 are corroborated by an independent data set³¹), and, RBFOX (RBFOX1/2/3) (P=0.0024, 20 targets, of which 12 overlap with FMRP), with RBFOX1 shown to be a splicing factor dysregulated in ASD^32,33 (Fig. 1a). These two pathways expand the complexity of the ASD neurobiology to post-transcriptional events, including splicing and translation, both of which would sculpt the neural proteome.

We found nominal enrichment for human orthologs of mouse genes encoding synaptic (P=0.031) and postsynaptic density (PSD) proteins³⁴ (P=0.046, Fig. 1a, 1b, Supplementary Tables 4, 5 and 6). Enrichment analyses for InterPro, SMART, or Pfam domains (FDR < 0.05 and a minimum of 5 genes per category) reveal an overrepresentation of DNA/histone-related domains: 8 genes encoding proteins with InterPro zinc finger (Znf) FYVE PHD domains (142 such annotated genes in the genome; FDR=7.6×10⁻⁴), and five with Pfam Su(var)3-9, Enhancer-of-zeste (SET) domains (39 annotated in the genome; FDR=8.2×10⁻⁴).

Integrating complementary data

To implicate additional genes in risk for ASD, we use a model called DAWN³⁵. DAWN evokes a hidden Markov random field framework to identify clusters of genes that show strong association signal and highly correlated co-expression in a key tissue and developmental context. Previous research suggests human mid-fetal prefrontal and motor-somatosensory neocortex is such a critical nexus for risk¹⁶, thus we evaluated gene co-expression data from that tissue together with TADA scores for genes with FDR < 0.3. Because this list is enriched for genes under evolutionary constraint, we generalized DAWN to incorporate constraint scores (Supplementary Information). When (a) TADA results, (b) gene co-expression in mid-fetal neocortex, and (c) constraint scores are jointly modeled, DAWN identifies 160 genes that plausibly affect risk (Fig. 2), 91 of which are not in the top 107 TADA genes. Moreover, the model parameter describing evolutionary constraint is an important predictor of clusters of putative risk genes (P=0.018).

Open in a separate window

Figure 2

ASD genes in neuronal networks

Protein-protein interaction network created by seeding TADA and DAWN predicted genes. Only intermediate genes that are known to interact with at least two TADA and/or DAWN genes are included. Four natural clusters (C1-C4) are demarcated with black ellipses. All nodes are sized based on degree of connectivity.

A subnetwork obtained by seeding the 160 DAWN genes within a high-confidence protein-protein interactome¹⁴ confirmed that the putative genes are enriched for neuronal functions. We kept the largest connected component, containing 95 seed DAWN genes, 50 of which were in the FDR < 0.3 gene set. The DAWN gene products form four natural clusters based on network connectivity (Fig. 2). We visualized the enriched pathways and biological functions for each of these clusters on canvases³⁶ (Extended Data Fig. 6). Many of the previously known ASD risk genes fall in cluster C3, including genes involved in synaptic transmission and cell-cell communication. Cluster C4 is enriched for genes related to transcriptional and chromatin regulation. Many TADA and DAWN genes in this cluster interact tightly with other transcription factors, histone modifying enzymes and DNA binding proteins. Five TADA genes in the cluster C2 are bridged to the rest of the network through MAPT, inferred by DAWN. The enrichment results for C2 indicate that genes implicated in neurodegenerative disorders could also play a role in neurodevelopmental disorders.

Emergent results

Amongst critical synaptic components found mutated in our study are voltage-gated ion channels involved in fundamental processes including propagation of action potentials (e.g., Na_v1.2 channel), neuronal pacemaking, and excitability-transcription coupling (e.g., Ca_v1.3 channel) (Fig. 1b). We identified, 4 LoF and 5 Mis3 variants in SCN2A (Na_v1.2), 3 Mis3 in CACNA1D (Ca_v1.3), 2 LoF in CACNA2D3 (α-δauxiliary subunits of L-type voltage-gated Ca²⁺ channels, including Ca_v1.3). Remarkably, three de novo Mis3 variants in SCN2A hit residues mutated in homologous genes in patients with other syndromes, including Brugada syndrome (SCN5A) or epilepsy disorders (SCN1A) (p.R379H and p.R937H). These arginines, as well as the threonine mutated in p.T1420M, cluster to the P-loops forming the ion selectivity filter, in proximity of the inner ring (DEKA motif) (Fig. 1c). Because homologous channels mutated in these arginines do not conduct inwards Na⁺ currents^37,38, p.R379H and p.R937H might have similar effect.

Two de novo CACNA1D variants (p.G407R and p.A749G) hit positions proximal to residues mutated in patients with primary aldosteronism and neurological deficits (Fig. 1d). The reported mutations interfere with channel activation and inactivation³⁹. Amongst variants found in cases, p.A59V maps to the NSCaTE domain, also important for Ca²⁺-dependent inactivation, while p.S1977L and p.R2021H co-cluster in the C-terminal proline-rich domain, the site of interaction with SHANK3, a key PSD scaffolding protein. Mutations in RIMS1 and RIMBP2, which can associate with Ca_v1.3, were found in our cohort (but with an FDR.0.3).

Chromatin remodeling involves histone-modifying enzymes (encoded by histone modifier genes, HMGs) and chromatin remodelers (‘readers’) that recognize specific histone post-translational modifications (PTMs) and orchestrate their effects on chromatin. Our gene set is enriched in HMGs (9 HMGs out of 152 annotated in HIstome⁴⁰, Fisher's exact test, P=2.2×10⁻⁷). Enrichment in the GO term ‘histone-lysine N-methytransferase activity’ (5 genes out of 41 so annotated; FDR=2.2×10⁻²) highlights this as a prominent pathway.

Lysines on histones 3 and 4 can be mono-, di-, or tri-methylated, providing a versatile mechanism for either activation or repression of transcription. Of 107 TADA genes, five are SET lysine methyltransferases, four are Jumonji (JmjC) lysine demethylases, and two are readers (Fig. 3a). RBFOX1 co-isolates with H3K4me3⁴¹, and our dataset is enriched in targets shared by RBFOX1 and H3K4me3 (P=0.0166, Fig. 1a, Supplementary Table 4). Some de novo missense variants targeting these genes map to functional domains (Extended Data Fig. 7).

Open in a separate window

Figure 3

ASD genes in chromatin remodeling

a. TADA genes cluster to chromatin remodeling complexes. Amino terminals of histones H3, H4 and part of H2A, are shown. Lysine methyltransferases add methyl groups, while lysine demethylases remove them. b. De novo Mis3 and LoF variants in CHD8. The box shows the outcome of RT-PCR and Sanger sequencing in lymphoblastoid cells for two newly identified de novo splice-site variants. The first mutation hits an acceptor splice site (red arrow), causing the activation of a cryptic splice site (red box), a four-nucleotide deletion, frame shift and a premature stop. The second mutation hits a donor splice site (red arrow), causing exon skipping, frame shift and a premature stop.

For the H3K4me2 reader CHD8, we extended our analyses in search of additional de novo variation in the cases of the case-control sample. By sequencing complete parent-child trios for many CHD8 variants, five variants were found to be de novo, two of which affect essential splice sites and cause loss of function by exon skipping or activation of cryptic splice sites in lymphoblastoid cells (Fig. 3b).

Given the role of HMGs in transcription, we reasoned that TADA genes might be interconnected through transcription “routes”. We searched for a connected network (seeded by 9 TADA HMGs) in a transcription factor interaction network (ChEA)⁴². We found that 46 TADA genes are directly interconnected in a 55-gene cluster (Extended Data Fig. 8) (P=0.002; 1,000 random draws), for a total of 69 when including all known HMGs (Fig. 4) (P=0.001; 1,000 random draws).

Open in a separate window

Figure 4

Transcription regulation network of TADA genes

Edges indicate transcription regulator (source node) and its gene targets (target node) based on ChEA network; interactions among only HMGs are ignored.

Examining the Human Gene Mutation Database we found that the 107 TADA genes included 21 candidate genes for intellectual disability, 3 for epilepsy, 17 for schizophrenia, 9 for congenital heart disease and 6 for metabolic disorders (Fig. 5).

Open in a separate window

Figure 5

Involvement in disease of ASD genes

Venn diagram to visualize the overlap in disease involvement for the TADA genes.

Conclusions

Complementing earlier reports, ASD subjects show a clear excess of de novo LoF mutations over expectation, with a pile-up of such events in a handful of genes. While this handful has a large effect on risk, most ASD genes have much smaller impact. This gradient emerges most strikingly from the contrast of risk variation in male and female ASD subjects. Unlike some earlier studies, but consistent with expectation, the data also show clear evidence for effect of de novo missense SNV on risk; for risk generated by LoF variants transmitted from unaffected parents; and for the value of case-control design in gene discovery. Indeed, by integrating data on de novo, inherited and case control variation, the yield of ASD gene discoveries was doubled over what would be obtained from a count of de novo LoF alone. Almost uniformly ASD genes show large constraint against variation, a feature we exploit to implicate other genes in risk.

Three critical pathways for typical development are damaged by risk variation: (1) chromatin remodeling, (2) transcription and splicing, and, (3) synaptic function. Chromatin remodeling controls events underlying the formation of neural connections, including neural neurogenesis and neural differentiation⁴³, and relies on epigenetic marks as histone PTMs. Here we provide extensive evidence for HMGs and readers in sporadic ASD, implicating specifically lysine methylation and extending the mutational landscape of the emergent ASD gene CHD8 to missense variants. Splicing is implicated by the enrichment of RBFOX targets in the top ASD candidates. Risk variation also hits multiple classes and components of synaptic networks, from receptors and ion channels to scaffolding proteins. Because a wide set of synaptic genes is disrupted in idiopathic ASD, it seems reasonable to conjecture that altered chromatin dynamics and transcription, induced by disruption of relevant genes, leads to impaired synaptic function as well. De novo mutations in ASD^11-15, intellectual disability⁴⁴ and schizophrenia²⁹ cluster to synaptic genes, and synaptic defects have been reported in models of these disorders⁴⁵. Integrity of synaptic function is essential for neural physiology, and its perturbation could represent the intersection between diverse neuropsychiatric disorders⁴⁶.

Extended Data

Extended Data Figure 1

Open in a separate window

Workflow of the study

The workflow began with 16 sample sets, as listed in Supplementary Table 1. DNA was obtained, and exomes were captured and sequenced. After variant calling QC was performed: duplicate subjects and incomplete families were removed; and subjects with extreme genotyping, de novo, or variant rates were removed. Following cleaning, 3,871 subjects with ASD remained. Analysis proceeded separately for SNVs and indels, and CNVs. De novo and transmission/non-transmission were obtained for trio data (published de novo from 825 trios^11,13-15 were incorporated). This path led to the TADA analysis, which found 33 ASD risk genes with q < 0.1; and 107 with q < 0.3. CNV were called in 2,305 ASD subjects.

Extended Data Figure 2

Open in a separate window

Expected number of ASD genes discovered as a function of sample size

The Multiple LoF test (red) is a restricted version of TADA that uses only the de novo LoF data. TADA (blue) models de novo LoF, de novo Mis3, LoF variants transmitted/not transmitted and LoF variants observed in case/control samples. The sample size (N) indicates either (i) N trios, for which we record de novo and transmitted variation, or (ii) N trios, for which we record only de novo events, plus N cases and N controls.

Extended Data Figure 3

Open in a separate window

Heat map of the numbers of variants used in TADA analysis from each dataset in genes with q < 0.3

Left panel, variants in affected subjects; right panel, unaffected subjects. For the counts, we only focus on de novo LoF and Mis3 variants, transmitted/un-transmitted and case/control LoF variants. These variant counts are normalized by the length of coding regions of each gene and sample size of each dataset (|trio|+|case| for left panel, |trio|+|control| for the right panel).

Extended Data Figure 4

Open in a separate window

Genome browser view of the CNV deletions identified in ASD affected subjects

The deletions are displayed in red if with unknown inheritance, in grey if inherited, and in black in un unaffected subjects. Deletions in parents are not shown. For deletions within a single gene, all splicing isoforms are shown.

Extended Data Figure 5

Open in a separate window

Frequency of variants by gender

Frequency of de novo (DN) and transmitted (TR) variants per sample in males (black) and females (white) for genes with q < 0.1 (upper panel), q < 0.3 (central panel), or all TADA genes (lower panel). The P values were determined by a one-tailed permutation test (*P < 0.5; **P < 0.01; ***P < 0.01).

Extended Data Figure 6

Open in a separate window

Enrichment terms for the four clusters identified by protein-protein interaction network

P-values using Mouse-Genome-Informatics/Mammalian-Phenotype (MGI-MP, blue), Kyoto Encyclopedia of Genes and Genomes pathways (KEGG, red), and Gene Ontology biological processes (GO, yellow) are indicated.

Extended Data Figure 7

Open in a separate window

De novo variants in SET lysine methyltransferases and JmjC lysine demethylases

Mis3 are in black, LoF in red, and variants identified in other disorders in grey (Fig. 5). JmjC, Jumonji C domain; JmjN, Jumonji N domain; JmjC, PHD, plant homeodomain; ARID, AT-rich interacting domain; SET, Su(var)3-9, Enhancer-of-zeste, Trithorax domain; FYR N, FY-rich N-terminal domain; FYR C, FY-rich C-terminal domain; PWWP, Pro-Trp-Trp-Pro domain; HMG, high mobility group box; AWS, associated with SET domain; Bromo, bromodomain; BAH, bromo adjacent homology.

Extended Data Figure 8

Open in a separate window

Transcription regulation network of TADA genes only

Edges indicate transcription regulator (source node) and its gene targets (target node) based on ChEA network.

Supplementary Material

Acknowledgements

Footnotes

References