Inherited and multiple de novo mutations in autism/developmental delay risk genes suggest a multifactorial model

Background We previously performed targeted sequencing of autism risk genes in probands from the Autism Clinical and Genetic Resources in China (ACGC) (phase I). Here, we expand this analysis to a larger cohort of patients (ACGC phase II) to better understand the prevalence, inheritance, and genotype–phenotype correlations of likely gene-disrupting (LGD) mutations for autism candidate genes originally identified in cohorts of European descent. Methods We sequenced 187 autism candidate genes in an additional 784 probands and 85 genes in 599 probands using single-molecule molecular inversion probes. We tested the inheritance of potentially pathogenic mutations, performed a meta-analysis of phase I and phase II data and combined our results with existing exome sequence data to investigate the phenotypes of carrier parents and patients with multiple hits in different autism risk genes. Results We validated recurrent, LGD, de novo mutations (DNMs) in 13 genes. We identified a potential novel risk gene (ZNF292), one novel gene with recurrent LGD DNMs (RALGAPB), as well as genes associated with macrocephaly (GIGYF2 and WDFY3). We identified the transmission of private LGD mutations in genes predominantly associated with DNMs and showed that parental carriers tended to share milder autism-related phenotypes. Patients that carried DNMs in two or more candidate genes show more severe phenotypes. Conclusions We identify new risk genes and transmission of deleterious mutations in genes primarily associated with DNMs. The fact that parental carriers show milder phenotypes and patients with multiple hits are more severe supports a multifactorial model of risk. Electronic supplementary material The online version of this article (10.1186/s13229-018-0247-z) contains supplementary material, which is available to authorized users.


Background
In 1943, Leo Kanner first described 11 children with "early infantile autism" as those who manifest "a powerful desire for aloneness" and "an obsessive insistence on persistent sameness" [1]. Since then, the term autism has evolved. Autism spectrum disorder (ASD) is now defined as a group of clinical heterogeneous disorders with deficits in social interaction and repetitive or restrictive behaviors as main clinical characterizations often accompanied by other impairments, such as intelligence and language deficits. Studies conducted on several continents (Asia, Europe, and North America) indicate a prevalence of approximately 0.12-2.64% [2][3][4][5][6].
Changes in diagnostic criteria have been accompanied by large-scale, genome-wide, and targeted sequencing analyses dramatically accelerating the discovery of candidate genes associated with ASD [7][8][9][10][11][12]. Detailed phenotype-genotype correlations have been established for several high-risk genes in the last several years, such as CHD8 [13], ADNP [14], DYRK1A [15], and POGZ [16], emphasizing both the clinical and genetic heterogeneity of ASD. Despite these advances, several limitations remain. Only a small fraction of the genetic risk has been defined, the penetrance of most mutated genes is unknown, and genotype-phenotype correlations in most cases are unclear.
Previously, we presented results from targeted sequencing of 1543 ASD probands from the Autism Clinical and Genetic Resources in China (ACGC) [17]. This study extends the cohort to 2926 probands (~2000 from trios or quads) with the aim to identify novel autism risk mutations, their associated phenotypes, and the transmission characteristics of potentially pathogenic mutations within families. In addition to discovering potential novel ASD risk loci in the Chinese population, we report the discovery of patients with de novo mutation (DNM) in more than one autism risk gene, suggesting that multiple genes may contribute to ASD pathology. Consistent with this model, we identify families with inherited likely gene-disrupting (LGD) mutations in genes primarily associated with DNM emphasizing the importance of clinical evaluation and follow-up genetic counseling.

Study samples
Study samples were selected from the ACGC [17]. This collection includes seven clinical referring centers (Additional file 1: Figure S1) and consists of~10,000 ASD familial DNA samples. Patients were diagnosed primarily according to DSM-IV/V criteria documenting additional comorbid conditions where possible. Of the 3120 ASD probands in the ACGC, 2276 represent complete parent-child trios or quads and the majority are simplex autism with no family history of ASD. Peripheral-blood DNA of all individuals with ASD and their parents, if available, was collected with informed consent by seven coordinating centers. Genomic DNA was extracted from the whole blood. All study procedures were in accordance with the ethical standards of the Institutional Review Board of the School of Life Sciences at Central South University, Changsha, Hunan, China.
Targeted sequencing of the ACGC Previously, we targeted 213 candidate genes for sequencing in 1647 probands (1543 probands after QC) (phase I), as previously described [17] (Table 1). We define recurrent mutations based on the presence of independent mutations in the same gene but not necessarily at the same site. The candidate gene set consisted of recurrent DNM calls (including LGD and missense) identified from exome sequencing of autism families-primarily the Simons Simplex Collection (SSC) and the Autism Sequencing Consortium (ASC). In addition, we also included genes with recurrent de novo (DN) events among intellectual disability (ID) and developmental delay (DD) probands and genes disrupted by copy number variants (CNVs) but for which no DN single-nucleotide variant has yet been identified.
In this study, we sequenced an additional 1473 probands (Table 1) from the ACGC using single-molecule molecular inversion probes (smMIPs) [7,18]. The samples were sequenced and analyzed using a staged approach. During the first stage, we combined the results from phase I sequence data from an additional 851 probands for 211 genes from the original study. During the second stage, we reduced the gene set and focused on targeted resequencing of 85 genes selected from the 211 genes according to the latest progress and our in-house data for the remaining 622 probands. We combined the smMIP data and quality-controlled data from both phase I and phase II for 3120 total patients. After controlling for coverage and uniformity, data from 2926 probands (phase II, 1383 probands) with 2154 from trios or quads (phase II, 1109 probands) passed quality control (QC) ( Table 1, Additional file 2: Data 1, Additional files 3 and 4: Figures S2 and S3). All subsequent analyses are based only on those genes and samples that passed the following QC metrics.

Quality control and variant calling
We applied a similar pipeline as in phase I for QC and variant calling. For QC, only the data from the individuals with more than 75% of the target (with coverage over 8×) and genes with more than 30% of the individuals (and a minimum of eightfold sequence coverage) were used for the following analysis. Sequencing was performed using the Illumina HiSeq2000 platform. Reads were aligned against hg19 with BWA-MEM (v0.7.13) [19] after removing incorrect read pairs and low-quality reads; single-nucleotide variants and indels were called with FreeBayes (v0.9.14). We consider variants exceeding eightfold sequence coverage and read quality over 20 (QUAL > 20) for annotation with Seattle-Seq, as previously described [20]. Variants with allele count (AC) ≤ 3 (0.1%) and allele frequency < 0.1% in the Exome Aggregation Consortium (ExAC) database were considered for subsequent analysis and validation.

Variant validation and microsatellite analysis
We validated variants by PCR amplification (500 bp amplicons) followed by Sanger sequencing. We tested transmission for all validated variants wherever parental DNA was available. To eliminate nonpaternity, we performed an independent microsatellite analysis using 8-13 autosomal microsatellite markers for each family. Microsatellite loci were amplified by PCR using fluorescently labeled primers, and labeled products were analyzed by capillary electrophoresis using GeneMarker software and the ABI 3730XL DNA Analyzer.

Statistical analyses
We applied two statistical models to assess the excess of DNMs for the 187 genes that passed QC. The first, a chimpanzee-human divergence (CH) model [7], considers the length of the gene and divergence between chimpanzees and humans to predict the expected number of DNMs, while the model denovolyzeR [21] estimates mutation rates based on trinucleotide context and accommodates known mutational biases such as CpG hotspots. Default parameters were used for both models with an expected rate of 1.5 DNMs per exome for the CH model [7]. P values in the ACGC-only analysis were corrected for 187 genes. P values in the ACGC, SSC, and ASC combined meta-analysis were corrected for 18,945 genes for which mutation rates have been estimated. A comparison of Combined Annotation Dependent Depletion (CADD) score distribution between missense DNMs within SCN2A and CHD8 from ASD patients and rare missense mutations within SCN2A and CHD8 from ExAC controls was performed using the Wilcoxon signed-rank test. The relationship between affected status and the number of DNMs in the SSC exome data was tested using a logistic regression model with affected status as a response variable and DNM numbers, father's birth age, and gender as predictor variables. This model took the form: logit[P(af-fectedStatus)]~DNM Number + FatherAgeAtBirth + Gender. Affected status represents a binary classification (proband (n = 2508) and sibling (n = 1911)) for each individual from the SSC whole-exome sequencing study cohort [9]. The DNM number represents the number of DNMs in each individual that result in an autosomal amino-acid alteration (LGD, missense, and in-frame indels). FatherAgeAtBirth is defined as the father's age in months at the time of the child's birth. The model was run with the glm function in the R statistical package for all individuals. Subsequently, males and females were separately tested using the same model without gender as a co-variable. The relationship between Social Responsiveness Scale (SRS) score and DNM number was tested using linear regression models with SRS score as the response variable, DNM number as the predictor variable, and gender as a co-variable. The relationship between nonverbal IQ (NVIQ) and DNM number was tested using the same model as for the SRS score. The relationship between seizures and DNM number was tested using the logistic regression model with seizure status (yes/no) as the response variable, DNM number as the predictor variable, and gender as a covariant. P values of these three phenotype analyses were corrected with false discovery rate (FDR) approach.

Mutation discovery
We discovered 2496 rare variants that are predicted to alter the amino acid sequence or disrupt a gene in phase II samples. We selected LGD mutations (n = 99) and missense mutations with a CADD score of equal to or greater than 30 (termed MIS30+; n = 133), because of their higher DN and pathogenic probability [17,22]
We applied two statistical models (see the "Methods" section) to assess the probability of excess of DNMs for the 187 genes. We identified 17 genes that reached significance for an excess of DNMs by the CH model and 13 genes by denovolyzeR (q < 0.05) in the ACGC cohort ( Table 2). Combining the ACGC analysis and ACGC-SSC-ASC analysis, 21 total genes reached significance by the CH model and 18 genes reached significance by denovolyzeR. ZNF292, which has not been reported as significant in previous studies, was implicated as a novel autism risk gene in the ACGC cohort by both models (q = 0.014, CH model; q = 0.016, denovolyzeR model) ( Table 2). We note that one LGD DNM was recently reported in the DDD study [24] and another LGD DNM was reported in an ID patient [25] (Fig. 1) clearly implicating this gene in ID as well as autism. In addition, we established recurrent LGD DNMs for RALGAPB by combining SSC and ASC exome data, although it is still not significant (q = 0.13, CH model; q = 0.33, denovolyzeR model) (Fig. 1, Table 2). Interestingly, an LGD DNM was reported in an epilepsy patient from the EPI4K study [26] (Fig. 1). CTNNBP2 was previously implicated for autism risk by the TADA (Transmission And De novo Association) test [11], and we now report DN significance based on the discovery of one missense and two LGD DNMs (q = 0.005, CH model; q = 0.02, denovolyzeR model) (Fig. 1, Table 2). We did not observe other potential pathogenic mutations in other ASD risk genes sequenced in this study in the patients with ZNF292, RALGAPB, and CTNNBP2 DNMs.

Clinical evaluation of ASD-relevant mutations
For patients carrying DNMs in autism risk genes, we reviewed the clinical details and made an attempt wherever possible to recontact families in order to assess phenotype, perform a physical examination, and assess co-occurring conditions. We observed significant ID, DD, and other comorbidities, such as behavior problems, in the well-defined or syndromic ASD genes-SCN2A, MECP2, FOXP1, ADNP, and ASXL3-which is consistent with the previous genotype-phenotype correlation analysis (Additional file 11: Data 5). Since we detected a relatively large number of probands with DNMs in SCN2A, which raises the possibility of dominant-negative or gain-of-function effects of the missense DNMs, we compared the phenotypes between patients with LGD DNMs (n = 6) and patients with missense DNMs (n = 11). However, we did not observe a significant difference between the two groups (Additional file 11: Data 5). While the number of patients with recurrent mutations in the candidate genes is too few to make definitive genotype-phenotype correlations with specific genes, several interesting Previously well-defined genes in neurodevelopmental disorders trends were observed. First, the majority of patients with severe DNMs and a cognitive assessment showed evidence of some form of intellectual impairment. Only TNRC6B, NCKAP1, and one of the two ZNF292 LGD DNMs occur in autism patients with an IQ in the normal range. Since microcephaly and macrocephaly have long been recognized as a co-occurring condition of ASD, we also assessed patients for abnormalities in head circumference (HC). In addition to CHD8, patients with LGD mutations within WDFY3, KMT5B, and GIGYF2 have notably larger HC. Three patients with LGD mutations in WDFY3 (1 DN, 1 inherited, and 1 undetermined inheritance) were identified with HC Z-scores of 2.5, 3.0, and 2.8, respectively. Two inherited LGD mutations were identified in KMT5B in two ACGC patients and both showed evidence of macrocephaly. Similarly, two GIGYF2 LGD mutations (1 DN and 1 undetermined) were identified in two patients and both were macrocephalic. In contrast, patients with DYRK1A, CDKL5, and MED13L LGD mutations have smaller HC, consistent with previous reports [15,27].

Multiple DNMs in ASD patients
During our analysis of the ACGC cohort, we identified two patients with two LGD DNMs in genes where each individually had reached significance for an excess for DNMs (Fig. 2a). Most notably, patient M01813 carried LGD DNMs in autism risk genes SCN2A and CDKL5, albeit the latter occurs near the terminal portion of the protein. The other patient, GX0477.p1, carried missense DNMs from MECP2 and RALGAPB. We assessed the frequency of such "double-hit" DNMs for our initial target set of 187 genes in both the SSC and ASC exome sequence datasets. We identified four additional SSC probands and four ASC probands with double-hit DNMs of which 6/8 pairs of genes were also classified as autism risk genes in the Simons Foundation Autism Research Initiative (SFARI) Gene database (Fig. 2a). Of those, only two missense mutations were presented in the ExAC database, although variants identified in ExAC do not indicate they are benign [28]. No such double-hits were identified in the unaffected SSC siblings. Although multiple-hit events are expected to occur more frequently in probands than siblings, we further investigated their effect on phenotype and gender differences among probands.
Using the SSC exome datasets, we examined the relationship between affected status and the number of DNMs and correcting for father's age at birth and gender (see the "Methods" section). Individuals with increasing DNM numbers are more likely to be affected Fig. 1 Distribution of DNMs in two potential novel autism risk genes (ZNF292, RALGAPB) and CTTNBP2. Dagger symbol indicates the DNMs reported in this study. ZNF292: p.R89* was reported in an ASC patient, p.L943Qfs*5 was reported in an ID patient, p. N1741Lfs*25 was reported in DDD study. RALGAPB: p.M289Vfs*3 was reported in an ASC patient, p.S1287* was reported in an EPI4K patient. CTTNBP2: p.V706Efs*14 was reported in an SSC patient (p = 7.24 × 10 −6 , OR = 1.19) (Additional file 12: Figure  S7). We further analyzed male and female samples separately by the same analysis. Both male (p = 0.006, OR = 1.15) and female (p = 0.0005, OR = 1.29) probands demonstrate a significant relationship between affected status and the number of DNMs (Fig. 2b), even after correcting for father's age at birth. Interestingly, we observed that, compared to male samples, females demonstrated an increased odds ratio for additional DNMs. This result is consistent with the "female protective model," where female probands may require additional mutational burden to reach a clinical diagnosis of autism [29]. We repeated the same analyses removing cases with no DNM under the same regression models. We observed that individuals with increased DNMs are more likely to be affected (p = 0.028, OR = 1.15) (Additional file 13: Figure S8a). When samples are stratified by genetic sex, we observe a slight increase but no significant effect among males (p = 0.6, OR = 1.09; p values corrected for two tests) (Additional file 13: Figure   S8b), while females demonstrate a stronger (p = 0.037, OR = 1.3) effect than the grouped analysis (Additional file 13: Figure S8c). These associations should be regarded as suggestive until replicated with cohorts of large sample size.
Finally, to test whether patients with more DNMs demonstrate more severe phenotypes, we explored the relationship between three autism-related phenotypes (autism symptom impairment, seizure, and NVIQ) and DNM numbers. Autistic severity, per parent report on the SRS [30], increased with increasing DNM numbers with marginal significance (p = 0.07, q = 0.11, OR = 1.49; Fig. 2c). Similarly, there is also a significant trend for increasing frequency of seizures (p = 0.01, q = 0.03, OR = 1.18; Fig. 2c) with the increase of DNMs. While patients with increased DNMs appear to have decreased IQ, this trend is not significant (p = 0.27; Fig. 2c).

Inheritance of potential high-risk mutations
Although this study focused primarily on DNMs, we also identified 40 LGD mutations within known autism genes a b c Fig. 2 Multiple-hit model for ASD. a Ten families and corresponding double-hit DNMs in the ACGC, SSC and ASC cohorts are shown. Dagger symbol indicates the genes listed as autism risk genes in SFARI; number sign indicates the variants presented in ExAC database. b The logistic histogram plot shows that both males and females have a higher probability of being affected with an increase in the number of DNMs even after correcting for paternal age effect. Females show a higher odds ratio (OR) than males for this additional DNM effect. c The plot shows the distribution of OR and the corresponding 95% CI of regression models, which predict affected status and different phenotype severity by DNM numbers where transmission was observed from supposed unaffected parents to ASD offspring (22 maternal, 18 paternal) (Fig. 3a). Specifically, we identified 12 inherited LGD mutations in genes where a burden of excess DNMs had been previously described, including CHD8 (3), KMT5B  (2), DSCAM (2), FOXP1 (2), SCN2A (1), ADNP (1), and  WDFY3 (1) (Fig. 3b). Similarly, we also discovered a CNV disrupting CHD8 through our clinical work. The 140 kbp deletion was transmitted from a father to both affected siblings and is absent from the Database of Genomic Variants. The deletion was further validated by array comparative genomic hybridization (Fig. 3c). Combined with a previously reported inherited LGD [13], we report five ASD families with inherited CHD8 LGD mutations (Fig. 3b).
To evaluate the ASD phenotypes of parents carrying potential pathogenic LGD mutations, we attempted to recontact all ACGC families carrying inherited CHD8 and KMT5B LGD mutations for clinical reevaluation. Wherever possible, we assessed IQ using the age-appropriate Wechsler battery, HC, and autism-related traits using the Broad Autism Phenotype Questionnaire (BAPQ). Three families with inherited LGD mutation of CHD8 (as well as one family with inherited LGD mutation of KMT5B) were successfully recontacted. All three parents (two mothers and one father) carrying CHD8 LGD mutations show lower NVIQ scores (< 80), which fall in the borderline range (Table 3). Their scores, however, are significantly higher than their affected children, who are generally severely impaired (IQ < 40), and in one case, cognitive deficits were so severe, an estimate of IQ could not be determined. Two of the children show increased HC consistent with the CHD8 phenotype [13], although only one carrier mother could be clinically classified as macrocephalic. Macrocephaly was also observed for the father carrying the 140 kbp deletion of CHD8 (Z-score = a b c Fig. 3 Inheritance of high-risk ASD genes. a The number of inherited and DNMs by autism risk gene within the ACGC. b Families with inherited LGD mutations or gene-disrupting CNVs in ASD high-risk genes, including CHD8 (5), KMT5B (2), DSCAM (2), FOXP1 (2), SCN2A (1), ADNP (1), and WDFY3 (1). c Genomic location of inherited CNV-disrupting CHD8 within an ASD family from the ACGC cohort. LGD, likely gene-disrupting 5.5). BAPQ data for all three CHD8 carriers suggests that parents carrying CHD8 LGD mutations (Additional file 14: Figure S9) show autistic traits with high scores across the domains of autism: behavior (rigid personality) and social communication (pragmatic language deficits). Similar to CHD8, the single parent carrying KMT5B LGD mutations also shows a lower IQ within normal range and features consistent with a broader autism phenotype. The data clearly argue that even among genes where there is strong evidence of increased DNM burden, an LGD mutation may not be necessary or sufficient to develop autism suggesting reduced penetrance or variability of expressivity. Overall, these data support the idea that instead of non-penetrance, the mutations are resulting in variable phenotypes consistent with a range of ASD manifestations.

Discussion
We have sequenced the protein-encoding region of 187 ASD candidate genes in 784 autism patients and 85 genes in 599 additional autism patients and performed a meta-analysis from 2926 ACGC patients to identify novel risk genes, mutations, and genotypephenotype relationships. Severe DNMs in SCN2A and CHD8 (including missense burden) account for 1.5% of the ACGC cohort with an additional 3.33% of the patients showing DNMs in an additional 36 genes, most of which reach DN significance. Patients with recurrent WDFY3, GIGYF2, and KMT5B LGD mutations show evidence of increased HC size, implicating novel macrocephaly-associated ASD genes. Consistent with this observation, it has been reported that loss of Wdfy3 in mouse leads to regional enlargements of the cerebral cortex [31].
ZNF292 was implicated as a novel ASD risk gene in this study. ZNF292 encodes a KRAB C2H2 zinc finger protein thought to function as a growth hormone-dependent transcription factor. Unfortunately, the biological function of ZNF292 is still unclear. Both patients with ZNF292 LGD DNMs meet diagnostic criteria for ASD (DSM-IV). Besides autism-related phenotypes, both showed delayed language development and abnormal EEG patterns. Patient M02463 (p.S832Ifs*28) showed mild ID and attention deficit hyperactivity disorder; however, patient M32023 (p.S1473Ffs*5) presented with normal IQ (108). Besides the three LGD DNMs reported in ASD patients, there are two LGD DNMs reported in DD and ID patients [24,25]. Despite this excess of DNMs, it should be noted that LGD mutations have been reported in ExAC; unfortunately, the phenotype of these individuals cannot be further assessed.
Although not yet significant, RALGAPB is also a promising risk gene for follow-up as recurrent LGD DNMs were identified in ASD. In addition, an LGD DNM was also identified in a patient with epilepsy. RAL-GAPB encodes a Ras-like GTPase-activating protein.
Several genes encoding the GTPase-activating protein have been associated with autism risk, such as SYN-GAP1, TSC2, ARHGAP32, and ARHGAP33 [17, 32, 33]. Of note, dysregulation of the Ras signaling pathway is a well-known etiologic factor from both genetic and functional studies associated with autism [34,35].
CHD8 is a well-described high-impact ASD gene with no LGD mutation identified in well-defined controls [13]. LGD mutations are estimated to be extremely rare in the general population, such as ExAC (minor allele frequency = 5 × 10 −5 ). Here, we describe five families with inherited LGD or CNV events. Parent carriers possess mild neurodevelopmental phenotypes, including borderline IQ and broader autism phenotypes suggesting variable expressivity as opposed to a non-penetrant mutation with no phenotypic consequence. One possibility for this variability may be that a heterozygous loss-of-function mutation can cause a mild phenotype but is, by itself, not necessary and sufficient to result in an autism diagnosis unless it occurs in conjunction with other risk mutation(s). Alternatively, carrier parents who were not previously diagnosed with a neurodevelopmental disorder may harbor protective genetic variants that dampen a more severe clinical presentation of ASD. From a clinical perspective, the two are difficult to discern, but in either scenario, early diagnosis and family counseling are particularly important. Consistent with this observation, our data also indicate that multiple DNMs in different autism risk genes within the same patient play an important role in both ASD etiology as well as disease severity. Although previous studies have shown that DNMs affect a continuum of functional outcomes [36], we investigated broader phenotypes, including occurrence of seizures. The observed association between severity of autism symptomatology and number of DNMs may provide some mechanistic insight into the heterogeneity of impairments in ASD. Such oligogenic effects have been observed previously for large CNVs associated with DD [37] and have been noted in several recent studies of ASD [38][39][40]. The model is distinct from a polygenic one because it puts forward that a relatively small number of rare or DNMs of large effect are primarily responsible for disease etiology and phenotypic severity, although the outcome may still be influenced by other factors, such as common variants, environment, or stochasticity during development. The analysis of SSC whole-exome sequencing data also reveals that, compared to male samples, females demonstrated an increased odds ratio for additional DNMs, although it should be noted that this study was limited to a relatively small number individuals where only exonic mutations were detected. Nevertheless, this result is consistent with the "female protective model," which has been proposed with both genetic and epidemiological evidence in ASD [11,29]. If this multifactorial model and female protective effect are more broadly applicable, the increased sensitivity afforded by whole-genome sequence may become more important than targeted approaches, such as exome or molecular inversion probe (MIP) sequencing, for diagnosis, discovery, and understanding of the genetic architecture and sex bias of ASD.

Conclusions
Targeted sequencing of candidate genes in the ACGC has identified novel ASD risk genes, mutations, and genotype-phenotype relationships. Among well-established autism risk genes primarily associated with DNMs, we identify ASD families where deleterious mutations are transmitted and find that parental carriers most often show a subset of milder phenotypes. We also identify families where patients carried DNMs in two or more autism risk genes and such individuals appear to be more severely affected. Both observations provide further support for a multifactorial model of ASD risk and suggest that a monogenic model of disease will be too simplistic even for the most penetrant causes of ASD.