Introduction

Neurodevelopmental disorders (NDDs) are a group of disorders primarily associated with neurodevelopmental dysfunction that include autism spectrum disorder (ASD), developmental delay (DD), intellectual disability (ID), and attention-deficit/hyperactivity disorder (ADHD)1. Children with NDDs experience difficulties with motor skills, learning and/or memory, language and/or nonverbal communication, and/or other neuropsychiatric problems. Considerable heterogeneity is common at both the phenotypic and genetic levels. With the advent of next-generation sequencing technologies, such as targeted sequencing2,3,4,5, exome sequencing6,7,8,9, genome sequencing10,11,12, and copy number variation (CNV) studies13,14, hundreds of genes and genomic regions have been implicated in NDDs almost exclusively based on the enrichment of de novo mutations (DNMs). But relatively few genes or loci have enough cases identified to prove statistical significance at the genome-wide level.

Ultra-rare and de novo gene-disruptive variants have been shown to play important roles in NDDs15. While DNMs from over 10,000 NDD families have been identified and cataloged16, the number of sequenced samples is still insufficient to reach the most stringent genome-wide significance levels, and samples from different ancestries and regions around the world are required to capture the whole picture of the genetics. Sample sizes in excess of 20,000 are projected to be necessary to reach significance levels by standard case-control criteria17. The discovery of large numbers of families with a disruptive variant in a specific gene, nevertheless, has facilitated establishing more meaningful genotype–phenotype correlations, such as in CHD818, POGZ19, and ADNP20. However, relatively few ASD or NDD genes have been interrogated at this level, emphasizing the need for conducting more candidate gene studies where patients and their families can be reassessed21.

Using single-molecule molecular inversion probes (smMIPs) is a relatively cheap and efficient approach to target sequence candidate genes in a large number of individuals where exome or genome sequencing is not feasible, or in situations where the amount of DNA is limited2. Here, we present targeted sequencing using smMIPs and analysis of the coding and splicing regions of 125 NDD candidate genes in a cohort with over 16,000 NDD patients from the international Autism Spectrum/Intellectual Disability (ASID) network, which includes 18 clinical groups across the world3. We identify 48 genes (25 newly reported) showing significant mutation burden of ultra-rare (MAF < 0.01%) gene-disruptive mutations (FDR 5%) by comparing to ExAC nonpsychiatric controls. Among these 125 targeted genes, we also identify 90 genes enriched for DNMs (FDR 5%) by reevaluating DNM excess in 17,426 NDD trios, including 6499 new autism trios. With this large-scale targeted sequencing effort, we further double the number of patients for many NDD risk genes and present deep phenotype–genotype correlations for seven NDD risk genes (CTCF, HNRNPU, KCNQ3, ZBTB18, TCF12, SPEN, and LEO1).

Results

Targeted sequencing and variant discovery

We initially selected 127 genes for targeted sequencing based primarily on published cases of recurrent DNM16, dividing the genes into two targeted sequencing panels (Fig. 1, Supplementary Data 1). The first panel (NDD1) consisted of 65 candidate genes selected for the first time in our study for sequencing in 17,832 NDD cases; the second panel (hcNDD) represented 62 genes, generally regarded as higher confidence NDD risk genes that had already been sequenced in a smaller subset (12,000–14,000) of ASID samples3,4,5. We applied this second panel to an additional 6,666 NDD cases in this study. We selected patient samples from the international ASID network of 18 clinical groups where ASD and DD/ID samples existed but neither exome nor genome sequence had been generated (Supplementary Fig. 1, Supplementary Table 1).

Fig. 1: Overview of study design.
figure 1

Targeted sequencing was performed in probands for two gene panels: NDD1 (63 genes) and hcNDD (62 genes). Gene and variant counts are after QC. The same categories of variants were retrieved from three previously published smMIP studies for 62 hcNDD genes. All smMIP variants were combined; redundant samples were eliminated and compared to the same category of variants from ExAC non-psych controls. The number of variants is after the exclusion of false positive variants and variants with insufficient coverage in ExAC. Mutation burden analysis identified 48 FDR significant genes (qmutBurden < 0.05, Benjamini–Hochberg correction for 125 genes), of which six reached FWER significance (pmutBurden < 1.25E−06, Bonferroni correction for 20,000 genes and two tests); DNMs of the 125 genes used in this study were identified from exome sequencing in 10,927 published NDD trios and 6,499 new ASD trios that combined as 17,426 NDD parent–child trios. A separate de novo enrichment analysis, using two statistical methods (CH model and denovolyzeR), identified 90 FDR significant genes (qdnEnrich < 0.05, Benjamini–Hochberg correction for 18,946 genes in CH model and 19,618 genes in denovolyzeR), of which, 61 genes reach FWER significance (pdnEnrich < 3.64E−07, Bonferroni correction for 19,618 genes and seven tests) for excess DNM. There is a significant overlap (40 genes) of the significant genes suggested by the two approaches. Then we performed genotype–phenotype correlation analysis for seven NDD risk genes (CTCF, HNRNPU, KCNQ3, ZBTB18, TCF12, SPEN, and LEO1) and present a clearer clinical picture of each gene.

Full size image

In panel NDD1, we designed 2,400 smMIPs to sequence the coding and splicing regions (exons plus five bases at each end) for 65 NDD candidate genes (Supplementary Data 2) among 17,832 NDD cases (8,738 and 9,094 cases with the primary diagnosis of ASD and DD/ID, respectively) (Supplementary Table 1). There were 1,538 samples (784 ASD and 754 DD/ID) and two genes (KCNQ2 and PAXX) that failed quality control (QC) based on read-depth coverage statistics (Supplementary Figs. 2, 3); these samples and genes were removed from subsequent downstream analyses. In total, we identified 31,659 putative single-nucleotide variants (SNVs) or insertions/deletions (indels) for 63 genes in 16,294 samples after QC. This included 586 ultra-rare (minor allele frequency [MAF] < 0.01%, i.e., allele count [AC] ≤ 3 in this study) severe variants, where 212 were likely gene-disruptive (LGD) variants (either a frameshift, nonsense, or canonical splice donor/acceptor variant) in 241 patients, and 374 were missense variants with a Combined Annotation Dependent Depletion (CADD) score22 greater than or equal to 30 (MIS30) in 465 patients. Using Sanger sequencing, we validated 183 LGD variants in 204 patients and 196 MIS30 variants in 233 patients with an overall validation rate of 96.7% (379/392) (Supplementary Data 3). Transmission was successfully assessed for 110 variants where we identified 40 DNMs with 29 de novo LGD (dnLGD), 11 de novo MIS30 (dnMIS30) variants, and 70 inherited variants in 73 families (three inherited MIS30 variants observed in two unrelated families) with maternally inherited variants in 37 families (30 MIS30 and 7 LGD) and paternally inherited variants in 36 families (23 MIS30 and 13 LGD). The majority (50/70) of the inherited variants were missense mutations. Limited clinical data are available for 28 carrier parents (Supplementary Data 5). Among the families where the parental phenotype data is available, one proband also carries a de novo missense variant (p.Arg1241Gln, CADDv1.3 = 15.4) in SHANK2 in addition to the paternally transmitted stop-gain variant (p.Arg860Ter) in CDK13, although the de novo variant is more likely to contribute to the proband’s autism. Most of the carrier parents (24/28) were classified as unaffected with no cognitive impairment, autism, or other psychiatric problems. The remaining four carrier parents show some clinical features related to the variant. One father, for example, who transmitted a MIS30 variant (p.Ser242Phe) in HNRNPR, had special education needs as he attended a school for individuals with learning disabilities but showed no obvious dysmorphic features. Similarly, a mother who transmitted a MIS30 variant (p.Arg339Gln) in CTCF showed a similar facial phenotype as the child but did not present with a clinical diagnosis of ID or ASD and was known to have attended regular school. A mother who transmitted a severe missense variant (p.Arg330Leu) in KCNQ3 was diagnosed with epilepsy but no cognitive impairment (Supplementary Data 5). Finally, one mother who transmitted a splice acceptor variant (c.1189-2 A > G) in TCF12 was diagnosed with long QT syndrome and glaucoma (like the patient) but this shared feature is unlikely related to DD observed in the child or the variant in question. These findings are consistent with the idea that such transmitted variants are by themselves not necessary and sufficient to develop DD but may rather be predisposing variants with a subset of parents manifesting more subtle phenotypes23.

In panel hcNDD, we resequenced 62 genes selected from our previous smMIP panels (Supplementary Data 1) for targeted sequencing with 3,575 smMIPs in 6,666 newly recruited NDD cases (3,562 ASD and 3,104 DD/ID) (Supplementary Table 1). All genes passed QC, but 455 DNA samples (199 ASD and 256 DD/ID) failed QC based on sequence coverage and were excluded from downstream analyses (Supplementary Figs. 2, 3). In total, we identified 72,811 SNV/indel variants for 62 genes in 6,211 patients after QC, including 213 LGD variants in 242 patients and 345 MIS30 variants in 426 patients. We validated 161 LGD variants in 172 patients and 170 MIS30 variants in 196 patients with a validation rate of 98.2% (331/337) for variants where Sanger sequencing was performed (Supplementary Data 3). Inheritance was assessed for 81 variants identifying 29 DNMs (21 dnLGD and 8 dnMIS30 variants) and 52 inherited (34 maternal and 18 paternal) variants. Ultra-rare severe variants were enriched ~2.5-fold among the hcNDD genes when compared to NDD1 genes for LGD (p = 4.82E−24, OR = 2.56 [2.14–3.08, 95% CI]) and MIS30 (p = 8.35E−39, OR = 2.49 [2.17–2.86, 95% CI]) variants (two-sided Fisher’s exact test), which reconfirms that these high-confidence genes usually have more severe variants in NDD cases.

Genes with an excess burden of ultra-rare severe variants

Since the 62 hcNDD genes were also previously sequenced in a subset (12,000–14,000) of ASID cases3,4,5, where we retrieved the same category of 1,120 ultra-rare severe variants with an overall similar validation rate of 97% (519/535) (Supplementary Data 4). We combined all of the retrieved data with our current sequencing in this study. Surveying the 125 genes across 16,000–19,000 NDD cases, there was a total of 2,113 ultra-rare severe variants (843 LGD and 1,270 MIS30 variants) from 2,621 patients (cases, Supplementary Data 5). In order to assess mutation burden, we extracted the same category of mutations corresponding to the smMIP capture regions for the 125 genes from ExAC (r0.3) controls24 without psychiatric disorders (n = 45,376) (controls, Supplementary Data 6). To quantify the population structure captured by our smMIPs, we conducted a principal component analysis (PCA) using the ultra-rare variants identified from our targeted sequencing, and also all the available single-nucleotide polymorphisms (SNPs) that overlap with our smMIPs from the 1000 Genomes Project (phase III high coverage) samples. We did not observe population-specific PCA clusters, suggesting that our ultra-rare variants are not stratified by different world populations (Methods). We excluded false positive variants and controlled for platform differences by removing variants with insufficient coverage between smMIP cases and ExAC controls (Methods). In total, 755 LGD and 1,177 MIS30 variants from smMIP cases, and 524 LGD and 1,810 MIS30 variants from ExAC controls were applied in the mutation burden analysis. We identified 48 genes with a significant excess of LGD and/or MIS30 (qmutBurden < 0.05, corrected ngenes = 125, variant count > 1) (Table 1, Fig. 2, Supplementary Data 10) in cases. Of these, six genes (ADNP, CHD8, DYRK1A, GRIN2B, POGZ, and SCN2A) also reached a more stringent significance threshold that pass exome-wide Bonferroni correction at the family wise error rate (FWER) for LGD variants (pmutBurden < 1.25E−06, corrected ngenes = 20,000, variant count > 1). Among the 48 significant genes, we identified 25 genes that show evidence of ultra-rare LGD and/or MIS30 (FDR 5%) burden for the first time in this large-scale case-control study, although 21 of these have been shown previously to show enrichment for DNMs (Supplementary Data 10).

Table 1 Genes with a significant burden for ultra-rare severe variants.
Full size table
Fig. 2: Significant genes identified from mutation burden and de novo enrichment analyses.
figure 2

a Mutation burden analysis identified 48 genes significant for LGD and/or MIS30 variants in smMIP sequencing compared with the ExAC (r0.3) non-psych subset controls; each dot indicates a gene and the color indicates the category of variant showing significance for the gene (red for LGD, blue for MIS30, and black for both LGD and MIS30). b The CH model and denovolyzeR show high concordance for genes with significant excess of DNM at both FDR and FWER levels. c A union set of 90 genes showing excess DNM (FDR 5%) in de novo enrichment analysis. Gray dashed box in top panel is shown in bottom panel for a zoom view. See Supplementary Data 10 for underlying data.

Full size image

Reevaluation of genes for excess DNMs

As the parent–child exome sequencing for ASD and DD/ID families has increased since the original selection of candidate genes, we also reassessed each of the 125 genes for excess DNM in a larger NDD combined set. In addition to the 537 dnLGD variants and 420 de novo missense (dnMIS) variants from previously published 10,927 NDD cases25 (Supplementary Data 8), we identified 99 dnLGD and 104 dnMIS (including 31 dnMIS30) variants in 6,499 new ASD patients from 5,911 complete families (4,761 simplex and 1,150 multiplex families) in our recent analysis of 27,270 SPARK exomes (unpublished data, https://sparkforautism.org/) (Supplementary Data 9). In total, there are 636 dnLGD and 524 dnMIS (including 201 dnMIS30) variants in the 125 genes from 17,426 NDD (12,123 ASD and 5,303 DD/ID) cases. We reevaluated the genes for excess DNM (dnLGD, dnMIS, dnMIS30, or de novo protein alteration [dnALT] variants that include dnLGD and dnMIS) using two statistical models (Fig. 1): a modified chimpanzee–human divergence model (CH model)4 and the denovolyzeR26 model as previously described25. Correcting for the total number of genes in each model, 81 genes show excess DNM in NDD patients according to the CH model (qdnEnrich < 0.05, corrected ngenes = 18,946, DNM count > 1) compared to 74 genes predicted to be enriched by denovolyzeR (qdnEnrich < 0.05, corrected ngenes = 19,618, DNM count > 1) (Fig. 2, Supplementary Data 10). The combination of both models identified 90 significant NDD candidate genes (union), and 65 genes were seen by both models (intersection). Applying a more stringent FWER significance (pdnEnrich < 3.46E−07, corrected ngenes = 19,618 in seven tests, DNM count > 1) identifies 61 union genes and 39 intersect genes (Fig. 2, Supplementary Data 10). This includes two genes (UIMC1 and GABRG2) firstly significant at a 5% FDR and seven genes (ANK2, TBR1, PHF12, TCF7L2, SETD2, CASZ1, and NSD2), which were significant at 5% FDR previously, that firstly reach FWER significance in this larger NDD cohort (Table 2, Supplementary Data 10).

Table 2 Genes reaching new de novo enrichment significance.
Full size table

Genotype–phenotype correlations

We successfully collected clinical records for 41 probands that carry ultra-rare severe variants in seven significant genes (CTCF, HNRNPU, KCNQ3, ZBTB18, TCF12, SPEN, and LEO1) from families that were available for recontact (Figs. 3 and 4, Supplementary Data 11). We also obtained clinical information for nine probands with dnMIS variants (2 in CTCF, 4 in KCNQ3, and 3 in ZBTB18) identified from the clinical trio exome sequencing at Baylor Genetics, and one DD patient with a dnLGD variant in CTCF that was identified from trio exome sequencing by the Antwerp group (Supplementary Data 11). We integrated the above clinical records with previously published reports and present a more comprehensive genotype–phenotype correlation assessment within the context of each gene (Table 3, Supplementary Data 12–18).

Fig. 3: Severe variants and the genotype–phenotype correlations in CTCF.
figure 3

a LGD (red) and MIS30 (blue) variants are depicted against a protein model for CTCF. Variants new to this study are shown above the protein while published DNMs from denovo-db (v1.5) are below. Variants are flagged with yellow lightning bolt if de novo. Annotated protein domains are shown (colored blocks) for the largest protein isoforms. b Heatmap depicts the common clinical features for patients carrying CTCF severe variants by using the specific HPO annotation (rows), which were retrieved from published studies and our cohort (columns). Phenotypic enrichment is shown according to the features’ recurrence labeled by the increment of color degree. The items with no data available were labeled with “-” and were excluded in the frequency analysis.

Full size image
Fig. 4: Distribution of severe patient variants in six genes.
figure 4

Protein diagrams are shown for HNRNPU (a) KCNQ3 (b) ZBTB18 (c) TCF12 (d) SPEN (e), and LEO1 (f) with the same display metrics that applied in Fig. 3. Validated LGD (red) and MIS30 (blue) variants are plotted. Variants listed above the protein model are new to this study, while the ones below were published previously. Paternal (black arrow) and maternal (green arrow) inheritance are shown if determined. A yellow lightning bolt denotes a de novo mutation.

Full size image
Table 3 Clinical recontact and detailed genotype–phenotype correlations.
Full size table

Germline deleterious variants in CTCF have recently been implicated in autosomal dominant DD/ID syndromic disorder (OMIM #615502) (Supplementary Data 12) with clustering of dnMIS30 variants occurring near the zinc-finger DNA binding domains associated with this protein27. We assessed 13 additional probands (including six with clustered dnMIS variants) from our study (Fig. 3). They are characterized by craniofacial dysmorphisms (9/10), thin vermillion border and lips (4/7), and feeding difficulties (6/11), and exhibit neonatal hypotonia (7/10). Along with these features, patients with CTCF mutations display a broader spectrum of developmental anomalies, including cardiac congenital malformations (1/8) and skeletal anomalies of toes/fingers (2/10). In addition to DD/ID (11/12), 54.5% (6/11) of the patients have a diagnosis of ASD and/or ADHD. The incidence of each phenotype in our probands (n = 13) is representative of the combined dataset, including published reports (n = 56) (Fig. 3).

HNRNPU mutations are now recognized as causative for early infantile epileptic encephalopathy-54 (EIEE54) syndrome (OMIM #617391), also referred to as HNRNPU-related disorder28. We observed seizures (3/3), DD/ID and ASD comorbidities (3/3), movement disorders such as stereotypies, e.g., hand flapping (1/3), and severe speech impairment (1/3) among our patients (Supplementary Data 13). We observed high ASD comorbidity (5/9) in patients carrying KCNQ3 mutations extending the phenotype which primarily associated with benign familial neonatal epilepsy. In our study, about half of the patients were diagnosed with benign familial infantile epilepsy (4/9) or DD (5/9) with or without seizures and cortical visual impairment (Supplementary Data 14). In contrast to HNRNPU, all mutations associated with KCNQ3 were severe missense mutations with no observation of a potential LGD mutation29. ZBTB18 is responsible for autosomal dominant mental retardation-22 (MRD22) syndrome (OMIM #612337), which is characterized by the features also seen in our patients such as moderate to severe DD/ID (7/7), ASD (2/7), speech delay (2/4), variable facial dysmorphisms (3/3), growth delay (2/4), and poor fine-motor skills (2/4) (Supplementary Data 15). TCF12 has been associated with craniosynostosis-3 syndrome (OMIM #615314). This phenotypic feature was observed in two of our patients, as well as other neurobehavioral phenotypes (DD/ID in 3/8 and ASD in 4/8 patients) (Supplementary Data 16).

We also investigated two additional candidate genes: SPEN and LEO1. To our knowledge, SPEN is newly identified in this study with a significant burden only for LGD variants (Table 1), while LEO1 shows excess DNM at both FDR and FWER levels (Supplementary Data 10). All patients with deleterious variants in SPEN show neurobehavioral impairment (Supplementary Data 17) (e.g., DD/ID in 6/7 and ASD in 5/7 patients in this study). Patients with a deleterious variant in SPEN show a more complicated clinical picture with other features, such as mild facial dysmorphism (4/4), muscular hypotonia, tall stature, poor motor coordination, and ocular abnormalities (3/4). Paternally inherited deletions of the LEO1 promoter were recently detected in three individuals with ASD11. Only two patients with disruptive mutations in LEO1 from our cohort could be recontacted, one showed some dysmorphic features and a minor cardiopathy plus global DD, while the other showed rather non-syndromic neurobehavioral features (Supplementary Data 18).

Discussion

Here, we report the results of large-scale targeted sequencing of 125 genes in over 16,000 pediatric NDD patients, with more than half the genes being screened in over 19,000 patients. We investigate these genes under a case-control mutation burden design and also test for DNM enrichment. Our comparison to ExAC controls identifies 48 genes as significantly enriched for ultra-rare severe variants in NDD patients (LGD and/or MIS30 variants, qmutBurden < 0.05, corrected ngenes = 125, variant count > 1). Additionally, 90 of the genes are enriched for DNMs in combined exomes of 17,426 NDD parent–child trios. There are 40 genes significant in both tests defining a subset of genes particularly relevant for future diagnosis of disease irrespective of inheritance patterns or availability of parental data. Overall, 78.4% (98/125) of the genes show some evidence of mutational burden in patients; notably, 61 genes remain significant at a more stringent level of FWER significance (61 with de novo enrichment, six of which were also detected from the case-control design) (Supplementary Data 10). In our targeted sequencing, 76% (95/125) of these genes have ultra-rare LGD variants identified in both patients with a primary diagnosis of ASD and DD/ID suggesting that these particular genes should be regarded as NDD genes as opposed to solely ASD or DD/ID risk genes.

In addition to the 98 genes significant by mutation burden analysis, or the de novo enrichment analysis, or both, there are additional candidates that trend toward increased mutational burden or de novo enrichment among NDD cases. For example, there are seven additional genes if considering a less stringent threshold (FDR 10%). One gene, NCKAP1, shows evidence of increased mutational burden for LGD variants (qmutBurden = 0.07), while six genes show excess DNM, namely SF3B1 (dnMIS qdnEnrich = 0.068 and dnALT qdnEnrich = 0.074), H2AC6 (dnMIS qdnEnrich = 0.053), and NFIA (dnALT qdnEnrich = 0.086) in the CH model and ARID2 (dnLGD qdnEnrich = 0.094), TNRC6B (dnLGD qdnEnrich = 0.097), and DNM1 (dnLGD qdnEnrich = 0.071) under the denovolyzeR model. Given the reported function of these genes and published case reports, it is likely that with increasing sample size these genes may achieve significance in the near future25. To test this, we expanded the number of parent–child trio exome sequencing cases with those from the SPARK pilot study30 and two recent publications from the ASC study8 and DDD study31 for a total of 48,281 NDD trios (excluding sample overlap and redundancy). Across those samples, four of the seven candidate genes reach some level of significance: ARID2 and DNM1 are significant for excess DNM at FWER significance, and H2AC6 and SF3B1 show excess DNM (FDR 5%). Overall, in this expanded de novo enrichment analysis, we estimate that at least 102 of the 125 genes in this study show a significant excess of DNM after adding the SPARK pilot, ASC, and DDD cohorts. Importantly, as additional genes become significant, our targeted sequencing studies will provide an important resource for future follow-up with clinicians and additional families to further investigate these genes.

We followed up clinically on seven candidates with the aim to develop or extend genotype–phenotype correlations. For example, CTCF, the CCCTC-binding factor, is a highly conserved zinc-finger protein that forms a multifunctional complex functioning in defining topologically associated domains, which are important for genome regulation and gene expression32. DNMs in CTCF have been described in patients with ID27. In this study, we identified three dnMIS30 variants based on smMIP screening (Supplementary Data 12) and characterized three additional DD patients with DNM in CTCF from the clinical trio exome sequencing at Baylor Genetics and the Antwerp group. Phenotypic assessments confirm features of the disorder and the importance of germline mutations in CTCF as causative for an autosomal dominant DD/ID syndromic disorder. The aggregate data highlight a striking clustering of deleterious missense mutations between the 2nd and 5th zinc-finger domain27 (Fig. 3). These functional domains have been described as the most important for making contact between the CTCF complex and DNA molecules and, as such, may represent useful targets for future therapeutic intervention33.

Other genes, such as KCNQ3, show a preponderance of severe missense mutations with half of the mutations mapping to the ion transport domain of the protein (Fig. 4). In our study, 5/9 of our patients with clinical information and a KCNQ3 variant are diagnosed with ASD (Supplementary Data 14), expanding the phenotypic spectrum of this gene as well as the main features of DD/ID and benign familial neonatal epilepsy34. All three of our recontacted patients with HNRNPU variants present with seizures (Supplementary Data 13), consistent with its association with epileptic encephalopathy and DD28. All four of our patients with a putative ZBTB18 (also known as RP58 or ZNF238) LGD variant present with DD/ID (Supplementary Data 15); this particular KRAB C2H2 zinc-finger protein has been described as a transcriptional repressor critical during brain development and neuronal differentiation35. Besides the previously reported large number of patients with TCF12 mutations36, we identified eight patients with a generally similar phenotype showing comorbid conditions of ASD and DD/ID in about half of the cases while craniosynostosis, which was originally primarily associated with this gene, was observed in only one-third of affected individuals (Supplementary Data 16).

Some of the newer candidates that have now reached or are nearing statistical significance for mutational burden still require much more extensive clinical follow-up and additional cases to further establish variant pathogenicity and refine the associated phenotype. Such is the case for RNA polymerase-associated protein LEO1, recently implicated in ASD11, although there are relatively few patients reported to date. We identified two additional individuals with stop-gain variants in LEO1, albeit with limited clinical information. Both of them are male—one patient presented with DD and the other with autistic behavior and ADHD with bilateral cryptorchidism (Supplementary Data 18). LEO1 is particularly intriguing in light of the recent observation that LEO1 interacts with the PAF1C complex in Drosophila to selectively transcribe expanded GGGGCC repeats in C9orf72-associated frontotemporal degeneration37. In addition, paternally inherited deletions of the LEO1 promoter11 and dnLGD variants in LEO1 have been reported in large cohort testing of DD and ASD patients7,9.

SPEN is another interesting candidate for further investigation. Haploinsufficiency of SPEN is considered a candidate for the 1p36 deletion syndrome phenotype38 and complete knockout of the gene in mice results in postnatal growth retardation and hypoplasia of the brain, especially involving the hippocampus and cerebral cortex39. We identified seven individuals in our study with DD and/or ASD with variable degrees of clinical information (Supplementary Data 17). Families with probands with SPEN LGD variants have no family history of DD/ID, learning disabilities, or neurological disease. For two patients where clinical data are more extensive, there is an indication of potential dysmorphology and skeletal abnormalities similar to previous reports. While the data, taken together, support the pathogenicity of SPEN LGD mutations, they also highlight a challenge going forward for the community. Unlike genes such as CHD8, POGZ, and ADNP, where large-scale screening has uncovered dozens of affected individuals for clinical evaluation and proved statistical significance at every level, the next tier of genes with ultra-rare and gene-disruptive DNMs will likely require screening of over 100,000 people. If only a handful of individuals with mutations in such genes are available, either from disparate labs with different standards of clinical reporting, or with incomplete family data, the pathogenicity determination may languish for years. Since we estimate that this next tranche of genes may account for more than half of the de novo gene burden associated with NDDs6, a more systematic effort involving targeted resequencing of large cohorts, database coordination (e.g., GeneMatcher), and dedicated researchers/clinicians willing to adopt such orphan genes and collate the clinical data are key. To help avoid false associations, whole-genome sequencing of such patients, their families, and controls may be particularly important to eliminate other genetic causes as contributing to disease and to understand the penetrance of the mutations under study.

Methods

Candidate genes

We considered two sets of genes: new candidates (NDD1) for investigation and high-confidence genes (hcNDD) that have been previously implicated in NDDs. Different criteria were used in selecting these two groups. In panel NDD1, we ranked and selected candidate genes for which no smMIP sequencing had been performed previously. We initially ranked all genes based on the DNMs from published NDD trios cataloged in denovo-db (v1.5), but excluding the following: genes associated with well-known syndromes based on OMIM, genes with extremely high-GC content, and genes with high counts of LGD and MIS30 variants in the ExAC non-psych controls. In total, 65 genes were selected for screening with: (i) 43 genes showing excess DNM25; (ii) 14 genes with evidence of autism sex bias40; (iii) six genes from a network analysis of high-functioning autism indicated previously3; (iv) and two genes (H2AC6 and H1-4) that were considered within a CNV candidate. In panel hcNDD, we continually reselected 62 top candidate genes from our previous smMIP panels3, mainly ranked by the reported number of DNMs from the published NDD trios in denovo-db (v1.5) and number of ultra-rare severe LGD and MIS30 variants identified in targeted sequencing of >13,000 NDD cases. We sequenced an additional 6666 newly recruited NDD cases that had not been previously sequenced using smMIPs. These served as positive controls of known disease genes in this study allowing for the discovery of additional cases for phenotypic evaluation. During the selection of these 125 genes, we evaluated the success rate of all smMIPs for each gene as part of our optimization experiments. We excluded genes, for example, where >20% of smMIPs failed to provide sufficient coverage even after 50-fold spike-in. We also balanced the total number of smMIPs per gene in each panel needed to achieve sufficient sequence depth. In particular, large genes requiring more than 200 smMIPs were triaged to allow a greater number of more moderate-sized genes to be considered. Supplementary Data 1 lists the genes with detailed selection criteria.

Study samples

Patient samples were obtained from the ASID network with informed consent. Only those not previously exome or genome sequenced were selected for targeted sequencing in this study. ASID is an international consortium3 that has expanded to include 18 centers around the world (Supplementary Fig. 1, Supplementary Table 1). The majority of samples were recruited from four sites (Adelaide, ACGC, Troina, and Leuven), as well as three new recruitment centers: an ASD collection from the University of Iowa (Iowa), an ID cohort from Charles University of Czech Republic (Charles), and an ASD cohort from the Italian Autism Network (ITAN). All targeted sequencing, Sanger variant validation, transmission analysis, and clinical recontact performed on the individuals in this study were approved by the University of Washington Institutional Review Board (IRB), in accordance with the ethical standards of the responsible local institutional and national committees. A PCA was used to quantify the population structure captured by our smMIPs. Samples in NDD1 generated two clusters; however, each cluster was composed of samples with mixed ancestries, and 15,659 samples (i.e., 96.1% of the total) were located under one heterogeneous cluster. In the case of hcNDD, a total of three clusters are observed; however, one of the clusters contains a heterogeneous mixture of 6161 samples (99.2% of the total). Overall, these observations suggest that the ultra-rare variants assayed by our targeted sequencing do not capture underlying population structure. Indeed, when we used all the available SNPs that overlap with our smMIPs from the 1000 Genomes Project (phase III high coverage) samples, we observed one large PCA cluster where 2,484 (99.2%) of the samples were included, once again supporting our previous observation that the genotypes of our ultra-rare variants are not stratified by different world populations. Hence, we expect our downstream case-control mutation burden analyses of the ultra-rare variants to not be confounded by population structure.

Targeted sequencing

All of the smMIP capture experiments, HiSeq 2500 sequencing, and Sanger validation experiments were performed at the University of Washington (Seattle, WA, USA), except for the ACGC cohort where experiments were carried out at the Center of Medical Genetics, School of Life Sciences, Central South University, Changsha, Hunan, China. In NDD1, 2,400 smMIPs were designed using MIPgen41 to cover all annotated RefSeq protein-coding exons and the splicing portions within 5 bp of flanking intronic sequence for all 65 genes. Oligos were ordered from Integrated DNA Technologies (IDT, https://www.idtdna.com/). smMIPs were pooled, rebalanced, and spiked-in at a relative concentration of 10X or 50X to improve sequence coverage for poorer-performing smMIPs where possible (Supplementary Data 2). A total of 17,832 NDD cases were sequenced using the balanced NDD1 panel. For hcNDD, 3,575 smMIPs from 62 genes were re-pooled from previous designs3 and tested for 6,666 newly collected NDD cases. smMIP capture libraries were barcoded and pooled with ~288 (3 × 96) samples and sequenced on a lane using an Illumina HiSeq 2500.

Variant annotation and validation

HiSeq data were processed according to the manufacturer’s instructions for base calling; variants were called using FreeBayes (version 1.0.2-6-g3ce827d) with its simplest operation (freebayes -f ref.fa aln.bam > var.vcf). Variants were filtered (QUAL > 20 and DP > 8) excluding common variants in dbSNP142 and then annotated using Ensembl’s Variant Effect Predictor42 (VEP, Ensembl GRCh37 release 94 - October 2018) with assembly GRCh37.p13 as the reference genome. Variants were annotated for all isoforms by VEP and those with the most severe consequence were selected for follow-up. Sanger validations were performed with ~300 bp PCR amplicons. CADD (v1.3) is a tool for scoring the deleteriousness of SNVs as well as indels in the human genome, and MIS30 variants are among the top 0.1% of the ~8.6 billion SNVs of the GRCh37/hg19 reference genome. LGD and MIS30 variants for the 62 genes in hcNDD were obtained from three previously published smMIP studies with same criteria applied (QUAL > 20, DP > 8, and MAF < 0.01% (AC ≤ 3)) (Supplementary Data 4). Similar variants from the targeted regions of 125 genes were obtained from the ExAC non-psych subset as controls with same filtering, i.e., QUAL > 20, DP > 363,008 (Avg. DP > 8), and MAF < 0.01% (AC ≤ 9) (Supplementary Data 6). All smMIP variants (this study and published) were merged with redundancy removed as variants with AC ≤ 3 retained for all subsequent analyses. dnLGD and dnMIS variants in the de novo enrichment analysis were extracted from SPARK-27K cases with ASD (n = 6,499) from complete families and the denovo-db (v1.5) NDD subset (n = 10,927). The published exome DNMs from SPARK pilot and ASC, together with recently released exome DNMs from DDD, were also included in the extended de novo enrichment analysis with sample overlap and redundancy removed. For cohorts like SSC and SPARK, for which the underlying exome data are available, duplicates were identified by running the KING software43, which uses identical by state (IBS) to estimate pairwise relatedness between samples. Any samples with a kinship value > 0.35 were considered to be identical and counted only once. Identical samples from the same cohort were also checked for reported monozygotic twin status. We identified one pair of SSC samples and eight pairs of SPARK samples as having a kinship value > 0.35. Note, samples in SPARK that overlapped with SSC samples were already removed in the final release by the SPARK Consortium. For other published cohorts, for which the underlying exome data are unavailable, the potential sample overlap identification, if applied, was described in each corresponding study. Like in the current DDD study, a total of eight duplicate samples were identified by collecting genotypes at 47 common exonic SNPs for every sample with a DNM found in another individual in the joint set; only one individual from each duplicate pair was kept with a final set of 31,058 samples analyzed. We also excluded sample overlaps reported in the literature. We excluded DD/ID samples in denovo-db (v1.5), which are also included as part of the current DDD study, and also excluded all 2,384 SSC samples in the ASC paper for potential redundancy with denovo-db (v1.5). These measures yielded a total of 48,281 NDD trios in the extended de novo enrichment analysis. To ensure uniformity, the same version of CADD score and VEP annotation were applied, and the analysis was restricted to the canonical transcript with the most deleterious annotation.

Statistical analyses

All statistical tests were performed using the R programming language (version 3.6.1). Benjamini–Hochberg FDR or Bonferroni FWER was applied when appropriate for multiple testing correction as described in the relevant sections. For mutation burden analysis, Fisher’s exact test (one-tailed) was used to compare the number of LGD and MIS30 variants from smMIP sequencing (cases) with those from the ExAC non-psych subset (controls), false positive variants by Sanger validation and variants with insufficient coverage (<90% samples with at least 10X coverage) in ExAC were excluded. The FDR significance threshold was set as qmutBurden < 0.05 where the q-value was corrected by Benjamini–Hochberg method for the total number of genes in this study (ngenes = 125); the FWER significance threshold was set as pmutBurden < 1.25E−06, which was calculated by 0.05/(20,000*2) and corrected by Bonferroni method for 20,000 genes in human genome and two tests performed (LGD and MIS30 variants). For de novo enrichment analysis, we applied both the CH model2 and denovolyzeR26 methods to assess the enrichment for four classes of DNM: dnLGD, dnMIS, dnMIS30, and dnALT. We applied denovolyzeR (v0.2.0) using default settings where dnMIS30 variants are not assessed; a modified CH model4 was applied to include the evaluation of dnMIS30 variants. Both methods apply their own underlying mutation rate estimates to generate the prior probabilities for observing a specific number and class of mutations for a given gene. Briefly, the CH model estimates the number of expected DNMs by incorporating chimpanzee–human coding sequence divergence and the length of the gene; denovolyzeR estimates mutation rates based on trinucleotide context, mutational biases such as CpG hotspots, and macaque–human gene comparisons. Default parameters were used for both methods, and the expected mutation rate of 1.8 DNMs per exome was set to the CH model as an upper bound baseline. The FDR significance threshold was set as qdnEnrich < 0.05 and corrected by the Benjamini–Hochberg method for the number of genes in each model (18,946 for CH model and 19,618 for denovolyzeR). The FWER significance threshold was set as pdnEnrich < 3.64E−07, which was calculated by 0.05/(19,618*7) and corrected by the Bonferroni method for 19,618 genes (the larger number of genes in two models) in seven tests performed (dnLGD, dnMIS, dnMIS30, and dnALT variants in CH model, and dnLGD, dnMIS, and dnALT variants in denovolyzeR).

Phenotypic assessment

Additional de-identified clinical records were obtained with informed consent for probands with ultra-rare severe mutations where the families were available for recontact (Supplementary Data 11). Clinical data were reviewed in consultation with the corresponding clinicians regarding the patient phenotypes and by analyzing existing or published clinical reports (Supplementary Data 12–18). For CTCF, we clustered and translated proband phenotype data into the corresponding Human Phenotype Ontology (HPO) annotation by using the Charité Browser; phenotypic enrichment analysis was performed based on the recurrence of the specific phenotype out of the total available clinical reports according to the HPO code (Fig. 3).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.