Introduction

Autism and autism spectrum disorders (ASDs) are characterized by deficits in social interactions, communication difficulties and the presence of restricted, repetitive and stereotyped patterns of behavior.1, 2 Autism is co-occurrent in disorders with known cytogenetic etiologies such as Fragile X, suggesting a role for genome structural aberrations; however, these account for <10% of cases to date.3 The remainder, often referred to as ‘idiopathic’ autism, are considered highly heritable with a 5–10% recurrence rate in siblings and a 60–90% concordance rate in monozygotic twins.4, 5 To date, genetic linkage and association studies have had limited success in identifying genetic risk loci for ASDs.6, 7, 8, 9, 10, 11 Most recently, genome copy-number variation (CNV) analyses have identified several candidate loci, including 16p11.2 and NRXN1 (neurexin 1).10, 12 The lack of common predisposition genes has led to the suggestion that rare variants constitute the majority of ASD risk, and that these variants collectively affect a limited number of common functional pathways. In particular, genetic and functional studies of candidate genes suggest the involvement of synaptic transmission in ASD etiology.13, 14 This evidence served as the basis for our hypothesis that inherited rare structural variants in genes from a number of closely related neurological pathways collectively contribute to a significant portion of disease risk in idiopathic ASDs.

Materials and methods

Study cohorts

The Autism Genetic Resource Exchange (AGRE) collection includes a total of 5431 affected and parental samples from 1000 families. These samples are grouped into four different sets based on the time when they were recruited: set 1 (660), set 2 (1373), set 3 (598) and set 4 (2800). Pedigree and affected status information on 4411 subjects are provided in the AGRE pedigree file at the AGRE website. Approximately 95% of these families are multiplex. The autism discovery and replication cohorts were derived from AGRE set 4 and sets 1–3, respectively. Cases designated as ‘autism’ and corresponding parents were included; cases designated as ‘spectrum’ and ‘not quite autism’ were not considered. Cases not of European descent and cases with known genetic syndromes or other non-idiopathic autistic causes were excluded. After applying these exclusion criteria and our genotyping and CNV detection quality control criteria (see below), the final discovery cohort included 1793 subjects (631 cases and 1162 parental samples), and the replication cohort included 1702 subjects (593 cases and 1109 parents). A number of parents (90 parents in the discovery cohort and 77 parents in the replication cohort) were excluded for failing to meet genotype quality measures. For the discovery cohort, 262 of the 335 families have multiple siblings included, for a total of 558 siblings. For the replication cohort, 218 of the 354 families have multiple siblings included, for a total of 457 siblings.

Healthy control samples were recruited from well-child visits conducted within the healthcare network of the CHOP (Children's Hospital of Philadelphia). Criteria for discovery control cohort inclusion were: age within age 3–18 years; screened for having no chronic illness as well as autism; and an extensive and continuous electronic health record with no indications of chronic health issues by diagnosis, prescribed medications, diagnostic tests or subspecialty encounters. All controls were unrelated as determined by health record and genotype assessment. After quality metrics used for the autism cohort were applied, the discovery control cohort comprised 1005 individuals of European descent, 723 of African descent and 47 of Asian descent. The replication control cohort consisted of 2026 healthy individuals comprising the CHOP CNV resource, as previously described.15 Gene-specific CNV analyses included healthy control individuals confined to European descent, whereas functional analyses included all ethnicities in order to identify and remove a larger constellation of nonpathogenic variants. The congenital heart defect cohort used as a control comparison consisted of 610 subjects with conotruncal defects and 896 parents of European descent who were uniformly assessed, genotyped and analyzed at CHOP. The congenital heart defect cohort was processed using the same protocols, quality measures and CNV size thresholds described for the ASD cohort, including restricting consideration of CNVs to those with unequivocal inheritance from a parent.

Genotyping

AGRE and control samples were assayed on the Illumina Infinium II Human-Hap550 BeadChip (Illumina, San Diego, CA, USA) as previously described.16 The standard Illumina cluster file was used for the analysis, which is generated at Illumina by genotyping 120 HapMap samples, running the BeadStudio clustering algorithm and reviewing single-nucleotide polymorphisms (SNPs) with poor performance statistics, including call frequency, cluster separation and Hardy–Weinberg equilibrium. For both control sets, ancestry informative markers available on the Human-Hap550 BeadChip17 were used to evaluate eligible subjects to determine ethnicity. Genotypes that failed quality measures for accurate CNV detection (call rate >98%) were excluded.

CNV detection and analysis

The Illumina BeadStudio 3.0. software package was used for initial CNV detection analysis. Log R ratios and B allele frequencies were first exported from BeadStudio. Log R ratio values were used as an additional sample-wide genotype quality control measure, and log R ratios with the s.d. above 0.35 were excluded from the study. CNV predictions were generated using CNV Workshop as described previously.15, 18, 19 To reduce the possibility of type I error, only deletions spanning 10 consecutive SNPs and duplications spanning 20 consecutive SNPs were included. These CNV size thresholds were selected based upon previous studies from our group15, 18, 19 and experience with samples undergoing array-based clinical diagnostics at our institution.20 Specifically, in these previous studies, our CNV detection method was shown to have a negligible false-positive rate using these thresholds. An increased threshold for duplications was used because of the difference in the degree of expected signal intensity variation between copy states for each probe: a twofold decrease for heterozygous deletion versus only a 1.5-fold increase for duplication. For comparative analysis purposes, two CNVs were considered equivalent if their genomic regions overlapped at least 80% reciprocally. All CNVs meeting quality and experimental design metrics were independently assessed by visual inspection with BeadStudio. For pathway analyses, all olfactory receptor genes were removed before consideration, as this gene family is rapidly evolving21, 22, 23 and tends to be associated with gene clusters that can over-represent such genes in functional enrichment studies. Careful analysis of enriched gene sets indicated no genomic region-specific bias from other functional gene clusters. For single gene analyses, apparently identical CNVs (>80% reciprocal overlap in genomic regions) from related siblings were removed before analysis. Significance testing for individual gene analyses used Fisher's exact test.

DNA for all AGRE and congenital heart defect samples was derived from immortalized cell lines, whereas all healthy control DNAs were derived from blood. Although acknowledging that the difference in DNA origin between cases and controls could theoretically have an effect on the results, we consider this to be highly unlikely, for the following reasons: (1) with the exception of the increased frequency of autism subject CNVs detected in the IgH-juxtaposed RNA gene LOC10013346, we saw no significant difference in genome-wide CNV rates in autism cases, and no other single gene was significantly enriched in cases; (2) real and presumptive de novo CNVs are excluded from our functional analyses; (3) our functional analyses would likely not be affected by increased cell line artifacts, as each gene is counted in a frequency-independent manner in these analyses.

Gene ontology analysis

Gene ontology analysis was carried out using the Explain Analysis System 2.4.2 (BIOBASE, Wolfenbüttel, Germany). RefSeq genes overlapping rare and inherited CNVs observed in autism cases were extracted from the UCSC (University of California, Santa Cruz) Genome Browser and input into the requisite pathway analysis models using default settings. Fisher's exact test was used when testing for enrichment of any gene ontology terms in the input gene lists. Mouse gene annotations were used for analysis because of a richer annotation set than for humans. All analyses were repeated with human annotations and yielded comparable results, although significance values were lessened because of the paucity of annotations for some genes. A Benjamini–Hochberg correction for multiple testing was applied for all functional enrichment analyses using the p.adjust function with parameter BH in the statistical package R.

Mouse phenotype ontology analysis

Mouse phenotype analysis was performed using Mouse Phenotype Ontology (MP) term annotations from the MGI (Mouse Genome Informatics) resource (Jackson Laboratory; Bar Harbor, ME, USA). Assigned MP terms for mouse models of genes overlapping rare and inherited autism CNVs were extracted from MGI. If a MP term was associated with a gene, all parental MP terms were assigned to the gene as well. Fisher's exact test was employed to measure the enrichment of autism CNV genes annotated with a given MP term compared with all genes annotated with a MP term. The P-values were used to rank the significance of enrichments.

CNV validation by quantitative PCR

Universal Probe Library (UPL; Roche, Indianapolis, IN, USA) probes were selected using the ProbeFinder v2.41 software (Roche). Probes and primer sequences selected for each CNV region are listed in Supplementary Table 12. Quantitative PCR was performed on an ABI 7500 Real Time PCR Instrument or on an ABI Prism 7900HT Sequence Detection System (Applied Biosystems, Foster City, CA, USA). Each sample was analyzed in quadruplicate in a 10 μl reaction mixture (100 nM probe, 200 nM each primer, 1 × Platinum Quantitative PCR SuperMix-Uracil-DNA-Glycosylase with ROX (Invitrogen, Carlsbad, CA, USA) and 25 ng genomic DNA). Values were evaluated using Sequence Detection Software v2.2.1 (Applied Biosystems). Data analysis was further performed using either the ΔΔCT method or qBase.24 Reference genes, chosen from COBL, GUSB, PPIA and SNCA, were included based on the minimal coefficient of variation. Data were subsequently normalized by setting a normal control to a value of 1.

Results

We conducted a modified case–control CNV study, restricting inclusion to probands and siblings of European descent with a diagnosis of idiopathic autism. We focused our analysis on rare inherited structural variants in the autistic subjects and potential enrichment of any functional categories assigned to genes overlapping these CNVs. Genotypes of autistic subjects were obtained from the AGRE25 (see Materials and methods) and were divided into discovery and replication cohorts based upon the time of recruitment.

In the discovery stage of the study, 631 autism subjects and 1162 parents of these subjects were compared with a healthy control set of 1775 children. We observed no significant differences in the frequency of CNV rate (P=0.109), CNV size (P=0.136), deletions (P=0.987) or duplications (P=0.144) between the subject and control groups. A total of 395 inherited CNVs (135 duplications, 256 heterozygous deletions and 4 hemizygous deletions of the X chromosome) were identified in 286 autism subjects that were not present in healthy controls (Supplementary Table 1). CNV transmission to probands was independent of parent gender (185/395 with paternal inheritance; P=1). For the purpose of this study, we considered and refer to these CNVs as ‘rare’, as the estimated prevalence of each is <0.05% in an ethnically matched control population. These CNVs overlapped a total of 382 genes excluding olfactory receptor genes. In all, 14 genes had CNVs in 5 unrelated individuals, and the remaining 368 genes had CNVs in 4 individuals. This autism-enriched, CNV-associated gene set (autism gene set) was analyzed to determine if it was enriched for particular biological or disease processes.

Gene ontology analysis revealed significant enrichment of the autism gene set for a number of biological processes (BPs) after correction for multiple testing (Table 1 and Supplementary Table 2). Many of these processes were associated with synaptic function and neurotransmission, including cell–cell signaling (P=1.97 × 10−7), transmission of nerve impulse (P=7.06 × 10−5), synaptic transmission (P=8.05 × 10−5), neuron adhesion (P=2.06 × 10−3), central nervous system development (P=4.30 × 10−3), regulation of neurotransmitter levels (P=4.88 × 10−3), neurotransmitter secretion (P=1.30 × 10−2) and neurotransmitter transport (P=2.95 × 10−2). Several cellular component (CC) terms were similarly enriched in the autism gene set, including presynaptic membrane (P=5.65 × 10−4) and synapse (P=1.11 × 10−3) (Table 2 and Supplementary Table 3). Both the BP and CC terms of significance were represented by sizable and only modestly overlapping gene sets (BP: n=114 unique genes, mean=3.2 terms/gene; CC: n=18 unique genes, mean=1.4 terms/gene), indicating considerable independence in implication of individual terms (Supplementary Figure 1). Autism CNV genes frequently associated with terms related to synapse function and neurotransmission included previous candidates CNTN4, NRXN1 and PARK2, as well as genes CNTN6, CTNND2, GRIN2A, NCAM2, GRM7, GRM8 and RCAN1.

Table 1 GO biological process terms significantly enriched in genes overlapping inherited rare autism CNVs
Table 2 GO cellular component terms significantly enriched in genes overlapping inherited rare autism CNVs

To demonstrate the specificity of the approach, the same CNV inclusion criteria and analysis strategies were re-applied using a large set of parent–child trios for congenital heart disease. The enriched pathway spectrum was almost entirely dissimilar in the cardiovascular cohort, and no function classes consistent with neuropsychiatric function were found significant, even without application of a multiple test correction (HM Xie et al., in preparation).

To replicate the functional analyses, we applied the same CNV inclusion criteria and methods as used for the discovery cohort to 593 additional idiopathic autism subjects of European descent, and 1109 corresponding parental samples, from AGRE. These were compared with CNVs from 2026 healthy individuals comprising the CHOP healthy control cohort.15 A total of 392 inherited CNVs (131 duplications, 256 heterozygous deletions, 1 homozygous deletion and 4 hemizygous X chromosome deletions) were exclusive to 270 of the autistic subjects (Supplementary Table 1). CNV transmission was again independent of parent gender (205/392 with paternal inheritance; P=1). Excluding olfactory receptors, these CNVs overlapped 387 unique genes, of which only 43 (11.1%) were present in the discovery autism CNV gene set.

A total of 12 BP and 3 CC terms reaching corrected significance in the discovery cohort were also enriched in the replication cohort after multiple test correction. As with the discovery cohort, the genes with CNVs in these enriched terms showed only moderate overlap between terms (Supplementary Figure 1). Furthermore, the autism CNV gene sets representing the replicated enriched terms demonstrated only a small overlap between the two cohorts, with only 23.7% (19/99) and 29.4% (5/17) of replication genes also present in the discovery cohort for BP and CC, respectively. For example, of 21 and 18 genes assigned to the BP term synaptic transmission with autism-specific CNVs in the discovery and replication cohorts, respectively, only 6 genes had CNVs in both cohorts. Among the genes with autism-specific CNVs found in both cohorts were known autism candidate genes CTNND2, CNTN4, NRXN1 and PARK2, along with glutamate receptors GRIN2A, GRIN2C and GRM7 (Table 5).

The discovery and replication cohorts were then combined into a single set of 1224 samples. Analysis of this combined cohort yielded results that were consistent with the initial findings. Autism candidate genes NRXN1 (P=1.46 × 10−2) and CNTNAP2 (P=3.51 × 10−3) demonstrated individual gene significance but failed to retain significance after multiple test correction (Table 4). In functional enrichment analysis, all BP and CC terms enriched in the discovery cohort retained multiple test-corrected significance in the combined cohort. The most significant functions were BP terms cell–cell signaling (P=6.14 × 10−7), transmission of nerve impulse (P=4.98 × 10−6) and synaptic transmission (P=4.98 × 10−6) (Table 1 and Supplementary Table 2); and CC terms synapse (P=1.49 × 10−5), presynaptic membrane (P=1.10 × 10−4) and dystrophin-associated glycoprotein complex (P=2.36 × 10−3) (Table 2 and Supplementary Table 3). There was again broad representation of genes among significant terms (Supplementary Figure 1).

As a complementary approach, we determined whether our autism CNV gene sets were preferentially associated with particular phenotypes assigned to mouse transgenic and null mice previously developed for these genes. For this analysis, we used MP term assignments from the MGI group26 that have been reported for 113 of the 382 discovery cohort autism CNV genes, and 98 of 387 such genes in the replicate cohort. Of 6406 functional terms, 12 were significantly enriched in the autism gene set both for the discovery and replication cohorts, and subsequently in the combined cohort (Table 3 and Supplementary Table 4). Many of the terms were highly consistent with an autistic phenotype, including abnormal (±CNS) synaptic transmission, abnormal motor capabilities/coordination/movement, abnormal motor coordination/balance, reduced NMDA-mediated synaptic currents, abnormal impulse conducting system conduction, decreased startle reflex, abnormal miniature excitatory postsynaptic currents and abnormal spatial learning. After correction for multiple testing, two closely associated terms retained significance in the replication and combined cohorts (abnormal synaptic transmission, combined P=2.52 × 10−2; abnormal CNS synaptic transmission, combined P=2.67 × 10−2). The generally lower significance values than observed for the BP and CC analyses likely reflect the relatively small number of genes with well-annotated phenotypes; however, the general trend of enrichment for synaptic function and neurotransmission was consistent with our previous analyses. Only 5 of 33 (15.2%) autism CNV genes within the 12 significant MP assignments were represented in both the discovery and replication cohorts, again suggesting widespread genetic heterogeneity.

Table 3 Mouse phenotype ontology terms significantly enriched in genes overlapping inherited rare autism CNVs

To determine the relative contribution of the unreported autism candidate genes to the enriched functional groups, each functional analysis was repeated after excluding a set of 205 genes either previously reported as autism risk candidates, or that reside in a genomic rearrangement associated with autism susceptibility (Supplementary Table 5). Despite the reduced power of these analyses, most BP and all CC and MP terms retained multiple test corrected significance (BP and CC) or uncorrected significance (MP) in both the discovery and replication cohorts (Supplementary Tables 6–11). Several additional terms consistent with autism pathophysiology were newly significant, including cell–cell adhesion (BP), response to mechanical stimulus (BP), cell junction (CC) and cytoplasmic vesicle membrane (CC). These results support representation of previously unattributed autism risk in the enriched autism gene sets.

To search for common variants, we determined the extent to which any single genes were represented by CNVs in autism cases relative to matched healthy individuals. We first compared our CNV set with findings from previous autism genomic studies, including a recent study by Glessner and colleagues11, 27, 28 that reported CNV findings for a set of autism spectrum disorder subjects partially overlapping our cohort. Combining the two cohorts, 25 inherited CNVs overlapped one or more previously implicated candidate genes (NRXN1, CNTN4, SUMF1, DISC1, PARK2, AUTS2 and CNTNAP2) or genomic regions (Table 4 and Supplementary Table 5). No gene had evidence of over-representation in the ASD cases, although formal significance testing was not possible because of the unavailability of inheritance status of CNVs in the control groups. The majority of CNVs detected in these genes did not substantially overlap CNVs detected in healthy controls.

Table 4 Analysis of inherited CNVs identified in genes previously implicated in autism

When expanding this approach genome-wide and considering all ASD subject CNVs regardless of inheritance status and presence in controls, only a single RNA gene on 14q32, LOC10013346, was significantly enriched for CNVs in ASD cases (27 overall CNVs versus none in the discovery control cohort, and 29 versus 4 in the replication control cohort). However, CNVs disrupting this gene in autism cases also spanned the adjoining immunoglobulin heavy chain locus, suggesting that these CNVs were likely cell culture artifacts. Moreover, only one of these CNVs in the discovery cohort and two in the replication cohort were present in parental samples. No other individual genes were significantly over-represented in either autism cohort.

CNVs overlapping five genomic regions (CTNND2, GRIN2A, NCAM1, KCNE2/FAM165B/KCNE1/RCAN1 and ICA1/NXPH1) were selected for experimental validation of CNVs based upon the following criteria: higher frequency seen in autism cases relative to controls; evidence for involvement in neurotransmission processes; and little or no previous implication for autism risk. As a validation control, we included two autism probands with known CNTN4 deletions and one proband with a known CNTN4 duplication.28 All 17 CNVs identified in probands and parental samples in these loci were experimentally validated (Supplementary Figure 2).

Discussion

Our functional enrichment analysis of the AGRE cohort identified a number of genes enriched for inherited CNVs in ASD subjects, several of which are known to have functions consistent with ASD etiology. Of particular interest were the ionotropic glutamate receptor GRIN2A and CNVs spanning both the islet cell autoantigen ICA1 (islet cell autoantigen 1) and the adjacent α-neurexin ligand NXPH1 (neurexophilin-1), as these genes were disrupted by CNVs in subjects from both autism cohorts and had no evidence of structural variation in healthy controls (Supplementary Figure 3). ICA1 is involved in glutamate receptor-mediated transmission,29 whereas the neurexophilin NXPH1 is a ligand for the protein product of the autism candidate gene NRXN1.30 Barnby et al.31 implicated GRIN2A (glutamate receptor, ionotropic, N-methyl D -aspartate 2A) as an autism candidate by gene-specific association, and further evidence suggests involvement in susceptibility to schizophrenia and attention deficit hyperactivity disorder.32, 33, 34 GRIN2A-null mice exhibit impairments in contextual memory and synaptic plasticity.35

Also of special interest were 12 autism CNV genes associated with significantly enriched functions each in the BP, CC and MP analyses (Table 5). Three of these genes (DOC2A, NRXN1 and PARK2) have previous support as autism susceptibility loci and validate our approach. The receptor tyrosine kinase ERBB4 controls maturation and plasticity of glutamatergic synapses, myelination and dopaminergic function, and has been suggested as having a role in schizophrenia risk.36 GRIN2A (described above), GRIN2C and GRM7 function as glutamate receptors, whereas UNC13A-null mice demonstrate impairment of glutamatergic synaptic vesicle maturation.37 The neurexin homolog NRXN3 is not well characterized but has roles in synaptic plasticity and synaptic vesicle exocytosis.38 SYN3 functions in synaptogenesis and the modulation of neurotransmitter release,39, 40 and SYN3-null mice exhibit deficits in social learning, contextual and cued fear conditioning and fear-potentiated startle.41 KCNMA1 encodes a calcium-activated potassium channel linked to glutamatergic synapses,42 whereas mice null for this gene exhibit cerebellar learning deficiencies.43 Disruption of KCNMA1 was identified in an autistic subject with a balanced translocation t(9;10)(q23;q22).44

Table 5 Autism CNV genes implicated in significant functional processes in both the discovery and replication cohorts

During review of this manuscript, Pinto et al.45 reported an analysis of rare CNVs in ASD cases relative to healthy controls. Both studies are consistent with the presence of many novel and previously reported loci, each with only a small number of distinct CNVs observed in the subjects considered. In comparing our results with the study of Pinto et al,45 there is little overlap between the sets of individual genes suggested as ASD candidates. These differences may be because of our exclusive focus on inherited events or might further support the presence of a substantial pool of potential ASD risk loci. Encouragingly, pathway enrichment findings between the two studies show greater consistency, with enrichment found in cell adhesion, central nervous system development, synaptic function and glutamate-mediated neuronal signaling processes. A follow-up study comparing inherited and de novo CNV functional attributes would allow determination of the degree to which each CNV class is functionally distinct.

Our results are consistent with the hypothesis that inherited autism risk is genetically highly heterogeneous, both from our observations of a lack of autism-specific CNVs overlapping any single gene in our analysis, and the lack of overlap between gene sets represented by the same autism-enriched functional terms in our two cohorts. Better estimates of the number of genes that can potentially confer risk, and the number of risk alleles typically required to confer sufficient risk in an individual, await higher-resolution genomic studies in larger cohorts, as well as platforms that more fully consider common CNVs. Our method, which has yielded results similarly consistent with disease-specific pathophysiology for attention deficit hyperactivity disorder,18 may be a viable approach for identifying risk signatures across functionally related processes. Caveats of this approach include insufficient specificity to identify individual candidate genes within moderate cohort sizes, and reliance upon incomplete and potentially biased gene knowledge. In addition, our approach only considered ASD variants not found in controls, which may oversimplify the biological role that certain variants might contribute to ASD etiology and is also dependent upon control cohort size. However, our results do consistently and independently implicate possible pathways, functions and gene sets that can be readily investigated by follow-up genomic and experimental investigation. Recent studies have implicated individual genes associated with synaptic function, including several neuroligands, as well as indications of glutamatergic pathway involvement.8, 13 Our current results extend these findings and, more importantly, for the first time provide statistical evidence that synaptic function and glutamate-mediated neurotransmission are contributing factors in autism etiology.