Rare coding variation provides insight into the genetic architecture and phenotypic context of autism

Abstract

Some individuals with autism spectrum disorder (ASD) carry functional mutations rarely observed in the general population. We explored the genes disrupted by these variants from joint analysis of protein-truncating variants (PTVs), missense variants and copy number variants (CNVs) in a cohort of 63,237 individuals. We discovered 72 genes associated with ASD at false discovery rate (FDR) ≤ 0.001 (185 at FDR ≤ 0.05). De novo PTVs, damaging missense variants and CNVs represented 57.5%, 21.1% and 8.44% of association evidence, while CNVs conferred greatest relative risk. Meta-analysis with cohorts ascertained for developmental delay (DD) (n = 91,605) yielded 373 genes associated with ASD/DD at FDR ≤ 0.001 (664 at FDR ≤ 0.05), some of which differed in relative frequency of mutation between ASD and DD cohorts. The DD-associated genes were enriched in transcriptomes of progenitor and immature neuronal cells, whereas genes showing stronger evidence in ASD were more enriched in maturing neurons and overlapped with schizophrenia-associated genes, emphasizing that these neuropsychiatric disorders may share common pathways to risk.

Fig. 1 |. Overview of SNV/indel and CNV rates in ASD by mode of inheritance and constraint.

a, The ASD cohort consisted of 49,049 family-based samples (15,036 cases) and 14,188 case-control samples (5,591 cases). One sample was a proband in one trio and a mother in another. b, The relative difference in PTV frequency between cases and unaffected controls (top) and average per sample variant count in unaffected controls (bottom) across inheritance classes (color) and LOEUF deciles (5,446 genes in top three deciles of LOEUF). Using a binomial test, cases were enriched for PTVs among the most constrained genes (lower LOEUF deciles), which weakened as negative selection against PTVs was relaxed (higher LOEUF deciles). c, Equivalent analyses were performed for missense variants annotated by MPC score and synonymous variants. Synonymous variants were not enriched in cases or controls, as evaluated via binomial tests. d, Benchmarking of the GATK-gCNV exome CNV discovery pipeline compared against WGS on overlapping samples achieved a sensitivity of 86% and PPV of 90% for rare CNVs (<1% site frequency) that spanned more than two captured exons (red line). e, The relative difference in variant frequency between cases and controls for deletions. Using binomial tests, we found that the enrichment of deletions (Del.) overlapping genes in the lowest LOEUF decile were stronger than PTVs in the same LOEUF deciles. f, Equivalent analysis for duplications (Dup.) demonstrated a similar pattern of enrichment compared with deletions but with more subtle relative differences. Statistical tests in b, c, e and f were two-sided binomial tests with 95% confidence interval error bars shown, P values (not corrected for multiple tests) and sample sizes are given in Supplementary Table 22.

Fig. 2 |. Contribution of CNVs to ASD by mechanism and genomic location.

a, CNVs included deletions (DEL) or duplications (DUP) of genomic segments and involved a subset of recurrent sites known as GD loci. GDs mediated by NAHR harbored recurrent breakpoints localized to flanking segmental duplications (Seg. Dup.), whereas non-NAHR GDs did not. b, De novo CNVs were highly enriched in affected cases compared with unaffected offspring (Fisher’s exact test) and the effect size was greater than that observed in de novo PTVs or de novo missense variants (logistic regression). c, ORs for de novo GD CNVs in probands compared with unaffected siblings, a subset of which have no observed de novo CNVs in unaffected individuals in this cohort (for example, 16p11.2 deletions, 15q11.2-q13 duplications). d, Analysis of all GDs (de novo and inherited) in ASD cases compared with GDs in a population-based cohort (UKBB) discovered using GATK-gCNV with identical parameters, with LOESS-smoothed bands of the 95% confidence interval of the OR in gray. e, Parent-of-origin analysis of de novo CNVs using binomial tests showed maternal bias for NAHR-mediated CNVs at GD regions, which was most pronounced for the 16p11.2 GD as previously described. Statistical tests in b and c were Fisher’s exact test with 95% confidence interval plotted as error bars, P values (not corrected for multiple tests) and sample sizes are located in Supplementary Table 22; statistical tests in d were Fisher’s exact test of carrier status in 13,786 unique ASD cases and 143,532 unique UK biobank controls, P values (not corrected for multiple tests) are located in Supplementary Table 10; statistical tests in e were binomial test with 95% confidence interval plotted as error bars, sample sizes and P values (not corrected for multiple tests) are located in Supplementary Table 22.

Fig. 3 |. Integrating variant types and inheritance classes improves association power and reveals mutational biases within candidate genes.

a, Our new implementation of the TADA model included de novo, case-control and rare inherited modules for each variant type: PTV, MisB, MisA, deletion and duplication. We leveraged information from ASD probands as well as unaffected siblings in evaluating the effect of de novo variants. b, The evidence of ASD association contributed by each variant type for each of the 72 ASD genes with FDR ≤ 0.001. Some genes were associated predominantly with missense variants and duplications (for example, PTEN and SLC6A1), suggesting mechanisms other than haploinsufficiency. c, Applying TADA to our aggregated ASD dataset yielded 72 genes at FDR ≤ 0.001, compared with 32 and 19 genes at the same threshold in previous studies on a subset of the samples^,. Our expanded TADA model improved the integration of available evidence of association and increased gene discovery at equivalent statistical thresholds on the same datasets. d, We quantified the relative contribution of variant class and mode of inheritance to these 72 ASD-associated genes, demonstrating that de novo PTVs and MisB variants represented the strongest contributions to the association signals. e, Association evidence (BF) was predominantly driven by de novo variants. The statistical test used in b was the extended TADA model.

Fig. 4 |. Relative contribution of evidence types in ASD risk genes.

a, The relative evidence of ASD association in the extended TADA model (log₁₀BF) for the 72 ASD risk genes (FDR ≤ 0.001) shown for likely loss-of-function mechanism (PTVs and deletions) on the x axis versus variants that may act via alternative mechanisms (missense variants and duplications) on the y axis. b, Plot of the relative association evidence from de novo (y axis) versus inherited (x axis) variation for the 72 ASD risk genes. c, Evidence for ASD association for the gene PLXNA1 was derived from de novo and inherited missense variants localized to the Plexin domain at the C-terminus of the Plexin-A1 protein. Statistical test in c was the transmission disequilibrium test. Obs/exp, observed/expected.

Fig. 5 |. Integration of ASD and DD datasets.

We performed meta-analysis of the ASD cohort with the 31,058 DD trios reported in Kaplanis et al. (n = 46,094 combined NDD cases). a, Relative difference of de novo PTVs in ASD and DD cohorts across deciles of constraint as measured by LOEUF. b, Relative difference of de novo missense variants in DD and ASD cohorts across categories of MPC scores. c, To explore overlap in association evidence across ASD and DD risk genes, we considered the 477 TADA-DD genes with FDR ≤ 0.05. We evaluated their P value distributions converted from TADA-ASD FDRs and observed nonuniformity suggesting that, in aggregate, 70.1% of these genes also had evidence of association with ASD. d, In the complementary analysis of the 185 TADA-ASD genes with FDR ≤ 0.05, we looked at their P value distributions converted from TADA-DD FDRs and again observed high nonuniformity, suggesting that, in aggregate, 86.6% of these genes had evidence of association with DD. e, Using PTV and MisB variant data, we devised a chi-squared statistic, denoted the C statistic, to measure if a gene has more observed variants in one cohort relative to the other. A mixture model was used to deconvolve the commingled distributions. f, We transformed the fitted mixture distribution into posterior probability for ASD enrichment. Using a cutoff of <0.01, we found 82 DD-predominant genes, while using a cutoff of >0.99 we found 36 ASD-predominant genes. Statistical tests: a,b, two-sided binomial test, with 95% confidence interval error bars shown; P values (not corrected for multiple tests) and sample sizes located in Supplementary Table 22; c,d, R v.3.5.3 package limma_3.38.3; e,f, mixture model.

Fig. 6 |. Single-cell data reveals differential neuronal layers impacted by ASD and DD genes.

a, A UMAP plot visualization after integrating two studies^, that provided single-cell gene expression of human fetal brains consisting of 37,000 cortical cells at 6–27 weeks postconception. Similar cell types from the two batches were grouped together while preserving cells unique to either study. See Supplementary Table 17 for unabbreviated cell type labels and classifications to ‘progenitor’ and ‘neuron’ types. b, Both ASD- and DD-predominant genes (right and left, respectively) were found to be enriched in interneurons and excitatory neurons compared to glial cells. Compared with DD-predominant genes, ASD-predominant genes were relatively more neuron-enriched than progenitor-enriched. The developmental trajectory of excitatory neurons was approximately recapitulated in the UMAP, starting with OPC and other progenitor cells and ending with maturing upper-layer enriched and deep layer excitatory neurons. For interneurons, InCGE and InMGE were precursors to InN. Statistical test used in b was Fisher’s exact test.