|Home | About | Journals | Submit | Contact Us | Français|
The identification of rare inherited and de novo copy number variations (CNVs) in human subjects has proven a productive approach to highlight risk genes for autism spectrum disorder (ASD). A variety of microarrays are available to detect CNVs, including single-nucleotide polymorphism (SNP) arrays and comparative genomic hybridization (CGH) arrays. Here, we examine a cohort of 696 unrelated ASD cases using a high-resolution one-million feature CGH microarray, the majority of which were previously genotyped with SNP arrays. Our objective was to discover new CNVs in ASD cases that were not detected by SNP microarray analysis and to delineate novel ASD risk loci via combined analysis of CGH and SNP array data sets on the ASD cohort and CGH data on an additional 1000 control samples. Of the 615 ASD cases analyzed on both SNP and CGH arrays, we found that 13,572 of 21,346 (64%) of the CNVs were exclusively detected by the CGH array. Several of the CGH-specific CNVs are rare in population frequency and impact previously reported ASD genes (e.g., NRXN1, GRM8, DPYD), as well as novel ASD candidate genes (e.g., CIB2, DAPP1, SAE1), and all were inherited except for a de novo CNV in the GPHN gene. A functional enrichment test of gene-sets in ASD cases over controls revealed nucleotide metabolism as a potential novel pathway involved in ASD, which includes several candidate genes for follow-up (e.g., DPYD, UPB1, UPP1, TYMP). Finally, this extensively phenotyped and genotyped ASD clinical cohort serves as an invaluable resource for the next step of genome sequencing for complete genetic variation detection.
Rare inherited and de novo copy number variations (CNVs) contribute to the genetic vulnerability in autism spectrum disorder (ASD) in as many as 5–10% of idiopathic cases examined (Devlin and Scherer 2012). Smaller intragenic CNVs (often called indels) can be involved, or CNVs can encompass an entire gene, and in some cases they can affect several genes as part of a genomic disorder (Lee and Scherer 2010). Screening for CNVs using microarrays has proven to be a rapid method to identify both large and small genomic imbalances associated with ASD susceptibility (Jacquemont et al. 2006; Sebat et al. 2007; Autism Genome Project Consortium 2007; Christian et al. 2008; Marshall et al. 2008; Pinto et al. 2010; Shen et al. 2010; Sanders et al. 2011).
A trend emerging from recent investigations in autism genetics is that significant heterogeneity and complexity exists. However, some highly penetrant risk genes (e.g., the SHANK, NRXN, and NLGN family members) and CNV loci (e.g., 1q21.1, 15q13.3, and 16p11.2) are now known (Devlin and Scherer 2012). Moreover, there has been some progress in identifying multiple mutations in single individuals (Schaaf et al. 2011), suggesting possible multigenic threshold models for ASD (Cook and Scherer 2008). Scientific designs enabling the discovery of a comprehensive set of genetic variants will be required to fully assess the genetic burden in ASD.
CNV detection in case and control cohorts have been performed using several different microarray platforms, including those utilizing single-nucleotide polymorphisms (SNPs) and others using comparative genomic hybridization (CGH) arrays (Redon et al. 2006; Graubert et al. 2007; Christian et al. 2008; Conrad et al. 2010; Pinto et al. 2011). There are advantages to each microarray platform. SNP-based microarray platforms are typically cost-effective and have the potential to analyze the genomic DNA sample at both the SNP genotype and copy number level, whereas CGH arrays are often used in a clinical diagnostic setting because of the better signal-to-noise-ratio achieved in comparison to SNP arrays (Shinawi and Cheung 2008). Multiple algorithms, some specific to certain array types, are available to detect CNVs, and they can vary considerably in the number of CNV calls made even for an identical array experiment. In recent comprehensive studies, authors have provided performance assessments of CNV detection platforms and methods by examining control DNA samples (Lai et al. 2005; Curtis et al. 2009; Hester et al. 2009; Winchester et al. 2009; Pinto et al. 2011). The consensus is that using multiple microarray technologies and prediction algorithms increases CNV detection rates.
In this study, we used a high-resolution Agilent CGH array comprising 1 million (1M) probes to assay for CNVs in a Canadian cohort of 696 unrelated ASD cases, 615 of which were genotyped previously on SNP arrays, including Illumina 1M single/duo (Pinto et al. 2010), Affymetrix 500K (Marshall et al. 2008), Affymetrix 6.0 (Lionel et al. 2011; and A. C. Lionel, unpublished data), and Illumina Omni 2.5M (A. C. Lionel, unpublished data) arrays. Our high-quality data, interpreted for the first time with CNV control data generated on 1000 individuals using the same 1M CGH array, enabled detection of numerous rare CNVs in each ASD sample examined and identified new potential ASD risk genes. These data allowed us to significantly expand the profile of genetic variants that are potentially causative of ASD and to identify novel molecular pathways affecting ASD vulnerability in this cohort.
The array CGH component of this study included 696 unrelated ASD cases, 39 affected siblings, and 42 parents (17 complete trios) that passed array quality control. All the ASD cases met the criteria for autism on one or both diagnostic measures – Autism Diagnostic Interview-Revised training and Autism Diagnostic Observation Schedule training. The samples came from three Canadian sites: Hospital for Sick Children in Toronto, Ontario; McMaster University, Hamilton, Ontario; and Memorial University of Newfoundland, St. John’s, Newfoundland. For the 696 unrelated ASD cases (571 males and 125 females), CGH experiments were performed on genomic DNA derived from blood for 354 cases, lymphoblastoid cell line DNA was used for 340 cases, and in two instances saliva was the DNA source used.
The SNP microarray component of the study used for comparison included 615 unrelated ASD cases, and the data from 433 of these experiments have been published (Marshall et al. 2008; Pinto et al. 2010). The CNVs detected in these samples were found using previously described stringent calling criteria defined for two different computer algorithms, which yields reliable CNV calls that can be experimentally validated >95% of the time (Redon et al. 2006; Marshall et al. 2008; Conrad et al. 2010; Pinto et al. 2010).
A control cohort of 1000 DNA samples from reportedly healthy donors (502 males and 498 females) were acquired from BioServe (Beltsville, MD) and are part of a collection of >12,000 control samples originally banked by Genomics Collaborative, Inc. (acquired by BioServe in 2007) for the purpose of large-scale genomic studies (Ardlie et al. 2002). Donors were consented and deidentified via a protocol approved by the institutional review board. The control DNA samples were derived from apparently healthy white donors older than 45 years of age. Health history information, documented at the time of consent, was used to select the samples based on the following attributes: body mass index between 15 and 35, blood glucose level <125 mg/dL, total cholesterol level between 100 and 300, systolic blood pressure between 100 and 150, and no major diseases (e.g., cancer and neurodegenerative diseases) or psychiatric disorders (e.g., alcoholism, mental illness, and depression). The control cohort DNA samples (also called PDx controls) used in this study are proprietary and currently biobanked at Population Diagnostics, Inc. (Melville, NY) for future microarray, sequencing and genotyping validation experiments.
The ASD test and reference samples were labeled with Cy5 and Cy3, respectively using Invitrogen BioPrime CGH labeling kit (Invitrogen, Carlsbad, CA). A pool of 50 sex-matched Caucasian control samples was used as a reference. The PDx control test and reference samples were labeled with Cy3 and Cy5, respectively and one sex-matched sample was used as a reference. All samples have been run on Agilent 1M CGH array according to manufacturer’s protocol (Agilent Technologies, Santa Clara, CA). The Agilent 1M CGH array consists of a total of 974,016 probes providing relatively uniform whole genome coverage. The arrays were scanned using the Agilent microarray scanner, and data were extracted using Feature extraction software version 10.5.1.1. The array experiments for ASD cases and PDx controls were run at The Centre for Applied Genomics (Toronto, Canada) and in the service laboratory of Oxford Gene Technology (Oxford, UK), respectively.
All array CGH data (both ASD cases and PDx controls) were analyzed in precisely the same manner using two programs—DNA Analytics v.4.0.85 (Agilent Technologies, Santa Clara, CA) and DNAcopy (Venkatraman and Olshen 2007)—to obtain high-confidence calls. At least five consecutive probe sets were used to call a CNV. For DNA Analytics, the aberration algorithm of Aberration Detection Method-2 was used with a threshold of 6.0, a minimum absolute average log2 ratio per region of 0.25, and maximum number of aberrant regions of 10,000 was used to identify all aberrant intervals. A nested filter of two was applied with subsequent removal of nested child calls using a custom script and retaining only the parent aberrations.
The other algorithm used to call CNVs was DNAcopy, from R Bioconductor package, which is a circular binary segmentation algorithm. The default settings and a log2 ratio cutoff of −0.41 and 0.32 for loss and gain, respectively, was used to call CNVs. Any segment with absolute median log2 ratio to median absolute deviation value less than 2 was filtered out. The CNVs detected by DNA Analytics and DNAcopy for each individual were merged using outer probe boundaries. As in our previous SNP microarray studies (Marshall et al. 2008; Pinto et al. 2010; Lionel et al. 2011), we defined a CNV as being ‘stringent’ if it was detected by both algorithms at the sample level. The stringent dataset was utilized for novel rare CNV discovery.
Experiments with poor derivative log ratio spread (DLRs > 0.3) were discarded. Any sample that had an excessive number of CNVs detected using either algorithm, measured by mean plus 3 SDs, was identified as an outlier and removed from further analysis. Any experiment that was a gender-mismatch was removed, and we excluded CNVs that were within centromere proximal cytobands. Twenty ASD cases that had CNVs larger than 5 Mb in size were removed from further downstream analyses of overlap with SNP microarrays, global rare CNV burden analysis and gene-set association tests.
Stringent CNV calls were classified as rare based on three separate stepwise comparisons with control datasets. First, a subset of stringent CNVs found at a frequency <0.5% for the total sample set (including 676 ASD cases and 1000 PDx controls) was compiled. Second, CNV data from an additional 4139 in-house controls were used to filter CNVs found at ≥0.1% frequency and where 50% by length overlapped with CNVs in the in-house controls. The additional controls consisted of 1782 subjects from the Study on Addiction: Genetics and Environment [SAGE (Bierut et al. 2010)], a total of 1234 unrelated controls from the Ottawa Heart Institute study (Stewart et al. 2009), and 1123 European controls from the PopGen study (Krawczak et al. 2006). The SAGE controls were genotyped with Illumina Human 1M-single BeadChip arrays and a subset of stringent CNVs detected by both iPattern (Pinto et al. 2010) and QuantiSNP (Colella et al. 2007) were used. The Ottawa Heart Institute and PopGen controls were genotyped with Affymetrix Genome-Wide Human SNP 6.0 arrays and the stringent subset consisted of regions that were detected by at least two of the three different CNV calling algorithms, Birdsuite (Korn et al. 2008), iPattern, and Affymetrix Genotyping Console. Third, the list of CNVs that overlapped 50% or more by length with known polymorphic regions in the genome (Conrad et al. 2010; McCarroll et al. 2008) were excluded. The final list of rare CNVs consisted of 1,884 CNVs in ASD cases and 2,299 CNVs in PDx controls.
For each sample for which CNV calls from SNP microarrays were available, stringent CNVs detected using the Agilent 1M array were overlapped with the stringent CNVs detected by corresponding SNP microarray experiments (Marshall et al. 2008; Pinto et al. 2010; A.C. Lionel, unpublished data). The CNVs from the SNP arrays were filtered to include only the regions with five probes or more. The CNVs were considered to be novel when 50% or more by length of the detected call was unique to a platform.
The experimental validation of novel CNV calls was carried out using SYBR-Green-based (Stratagene) real-time quantitative polymerase chain reaction (qPCR) method by at least two independent PCR assays. Each assay was conducted in triplicate, with one set of primers corresponding to the region of interest and the other mapping to a control region on FOXP2 at 7q31.1 (serving as a negative diploid control). The ratio between the test and control regions were then determined using standard curve method, and a fold-change less than 0.7 was confirmed as deletion and greater than 1.3 was confirmed as duplication. The parents and siblings were also tested for inheritance and segregation of CNVs, respectively.
Of the 676 unrelated ASD cases that passed QC, 615 of them had SNP microarray data available [Affymetrix Genome-Wide Human SNP Array 6.0, Affymetrix GeneChip Human Mapping 500K Array Set (Commercial and Early Access), Illumina Human1M and Human1M-Duo, and Illumina HumanOmni2.5-4 BeadChip]. We combined this subset of samples with the HapMap3 data set (International HapMap Consortium 2010), which includes individuals of different ethnicity (e.g., 183 Utah residents with Northern and Western European ancestry, 91 Toscani from Italy, 89 Han Chinese from Beijing, 181 Yorubans from Nigeria). For the ancestry analysis, only a subset of ~30,000 autosomal, non-MHC SNPs that were common to all the platforms were used. We analyzed this data set using multidimensional scaling (MDS) as implemented in PLINK (Purcell et al. 2007) to identify a putative group of European ASD cases. To avoid confounding ancestry issues, we only used European ASD cases for global rare CNV burden analyses and gene-set association tests.
CNVs mapping to genes were considered for burden analysis whenever at least one exon was encompassed by or impacted by the CNV. For each subject, we calculated the number of rare CNVs, the log10 of the sum of rare CNV sizes, and the total number of genes harboring a rare CNV. CNV burden differences were assessed by comparing the distribution of these statistics between case and control subjects. Specifically, differences in distribution over case and control subjects were tested using the nonparametric Wilcoxon exact test (as implemented in the R package exactRankTests), whereas the magnitude of the difference was assessed looking at the ratio of case and control means. CNV burden differences were assessed for all rare variants together, as well as for deletions and duplications only.
Gene-sets representing pathways, functional annotations and protein domains were tested if they were more frequently affected by rare CNVs in ASD cases compared with controls using the Fisher exact test [FET (Pinto et al. 2010)]. Gene-sets were compiled from Gene Ontology annotations (downloaded from NCBI in April 2011), pathway databases (KEGG, Reactome, bioCarta, NCI in March 2011), and protein domains (PFAM, March 2011). We only tested gene-sets with member genes numbering between 25 and 750: 2456 in total, 1939 from Gene Ontology, 414 from pathways, and 103 from PFAM domains. We found that small gene-sets decrease the statistical power of the analysis, whereas larger gene-sets tend to have a very broad biological scope without much useful biological meaning. For each gene-set, we built a contingency matrix with subjects as sampling units. Subjects were categorized as (1) cases or controls and (2) having at least one gene-set harboring a rare CNV or not. On the basis of this contingency table, we tested higher prevalence of rare CNVs in autism probands vs. controls using a one-tailed FET. Any subject in which a rare CNV affected more than 20 genes was not considered for the gene-set analysis because these individuals may have a broader set of gene functions perturbed by rare variants, which may consequently reduce the quality of gene-set analysis results. The FET nominal p value was corrected for multiple tests using a case/control class permutation procedure to estimate an empirical false discovery rate (FDR). We performed 5000 permutations of case-control labels and for each permutation we tested gene-sets following exactly the same procedure. The p values were ranked from lowest (most significant) to highest (least significant) and for each p value we computed the empirical FDR as the average number of gene-sets with equal or smaller p value over permutations. We have selected 25% as the empirical FDR significance threshold. Significant gene-sets were visualized using the Enrichment Map Cytoscape plugin (Merico et al. 2010).
Of the 696 unrelated ASD cases examined by array CGH, 20 were found to carry CNVs larger than 5 Mb (Supporting Information, Table S1) and excluded from further analysis. Four of these cases were known to have Down syndrome as well as an ASD diagnosis, 14 carried large chromosomal abnormalities previously detected through karyotyping and genotyping with SNP microarrays, and two cases harbored cell line artifacts. The rare stringent CNVs from the remaining 676 unrelated ASD cases were used for gene burden analysis, gene-set association tests and comparison with CNV data from other SNP microarrays to identify novel CNVs (Figure 1, Table 1).
We also examined the CNV content of the 1000 PDx controls and observed a significantly higher average number of CNVs per sample (p value < 2.2e-16) compared with the ASD cases (Table 1). This finding is likely attributable to the use of a different reference DNA, a single sex-matched individual, in the PDx CGH experiments rather than a pool of 50 sex-matched individuals used as a reference in the ASD CGH experiments. However, when we focused on rare variants (as defined in the materials and methods section), we found a smaller yet significant difference in the opposite direction (p value 0.002287; Table 2). The presence of a significant change in the relation between class (ASD, control) and CNV number after restricting to rare variants was further confirmed using a quasi-Poisson and linear regression model (class variable and rare variant filter variable interaction p value < 2.2e-16). Of the 676 unrelated ASD cases and 1000 PDx controls, 630 ASD cases (93%) and 896 PDx controls (90%) had at least one rare variant detected on the CGH 1M array (Table 2).
Of the 676 unrelated ASD cases, 615 were genotyped previously with SNP microarrays including 26 cases on Illumina Human Omni 2.5M-quad (2.5 million probes), 234 cases on Affymetrix Genome-Wide Human SNP Array 6.0 (1.8 million probes), 262 cases on Illumina Human 1M single infinium chip (1 million probes), 11 cases on Illumina 1M duo array (1 million probes), and 82 ASD cases were genotyped on lower resolution Affymetrix Mapping 500K chip set (500,000 probes). The other 61 unrelated ASD cases were only run on the Agilent 1M array and thus not included in this CNV platform comparison analysis. For the 615 samples run on both SNP and CGH array platforms, we performed a sample-level 50% one-way overlap of stringent Agilent 1M CNVs with the stringent CNVs from the SNP array platforms. We found that 64% of the Agilent 1M CNVs are novel with respect to the CNVs detected on the SNP microarrays (Figure 2). A more detailed comparison of CNVs detected by the 1M CGH array vs. the various SNP array platforms is shown in Table 3. For example, a comparison between CNVs detected using similar resolution arrays Agilent 1M and Illumina 1M showed that, on average, 24 novel CNVs/sample were detected by the Agilent 1M CGH array whereas only eight novel CNVs/sample were detected by the Illumina 1M SNP array. This platform comparison analysis suggests that use of multiple microarray platforms provides complementary data as every microarray platform detects a unique set of novel CNVs.
Rare ASD CNVs that were not detected by SNP microarrays but were detected by the Agilent 1M CGH array were defined as novel rare CNVs (946 of 1884 rare CNVs detected on the CGH 1M array). We examined the size distribution of these 946 novel rare CNVs detected in the ASD cohort and found that approximately 75% of them are <30 kb in size (Figure S1). These smaller CNVs that went undetected by SNP microarrays were missed due to insufficient probes coverage at these loci and because of the lower signal-to-noise often observed for SNP array platforms (e.g., SNP-based arrays are optimized for robust genotyping call rates, which may minimize quantitative probe response to copy number variation). Only 25% of the novel, rare CNVs were >30 kb in size and were mostly missed by previous studies due to insufficient probe coverage, CNV-calling parameters and analysis algorithms chosen to define a CNV as stringent.
Rare ASD-specific CNVs were defined using a total of 5139 controls (1000 PDx controls summarized in Table 2 and 4139 in-house controls previously reported in Krawczak et al. 2006, Stewart et al. 2009, and Bierut et al. 2010). A total of 1884 rare CNVs were detected in ASD cases, of which 946 of them were novel as compared with previous SNP microarray studies (Table S2). Using qPCR, we validated 117 of 132 (88.6%) novel and rare stringent CNVs that were tested. Of the 946 novel rare ASD CNVs, 57 CNVs are reported in Table 4 that correspond to overlapping CNVs in two or more unrelated ASD cases (32 cases at 14 loci), recurrent CNVs (i.e., same breakpoints) in two or more unrelated ASD cases (24 cases at 11 loci), or are a de novo event [1 case (Table 4)]. Some of the overlapping/recurrent CNVs impacted previously identified ASD genes such as DPYD, RGS7, NRXN1, CNTNAP5, ERBB4, GRM8, NRXN3, YWHAE, and DMD (Serajee et al. 2003; Wu et al. 2005; Autism Genome Project Consortium 2007; Marshall et al. 2008; Bruno et al. 2010; Pagnamenta et al. 2010; Pinto et al. 2010; Carter et al. 2011; Vaags et al. 2012), whereas others were novel, including RERE, NCKAP5, ROBO2, DAPP1, POT1, LEP, PLXNA4, CHRNB3, ZNF517, MIR3910-1/MIR3910-2, CIB2, MMP25/IL32, MYH4, RAB3A/MPV17L2, SAE1, and SYAP1 (Figures 3 and and44).
One of the ASD-specific CNVs was a maternally inherited duplication at 15q25.1 in three unrelated ASD cases (117395L, 94478, 132199L; Figure 4A) disrupting the exon of CIB2 (Calcium and integrin binding family member 2; 3/696 cases vs. 0/5139 controls; FET two-tailed p = 0.001691). The transcript and protein of CIB2 gene is found to be present mainly in the hippocampus and cortex of the brain (Blazejczyk et al. 2009). The encoded protein of this gene is shown to be involved in Ca2+ signaling, which controls a variety of processes in many cell types. In neurons, Ca2+ signaling maintains synaptic transmission, neuronal development and plasticity (Blazejczyk et al. 2009).
In another two unrelated ASD cases (115813L, 117463L), we identified a recurrent CNV of size 45.7 kb (2/696 cases vs. 0/5,139 controls; FET two-tailed p = 0.01421). The duplication disrupts four exons of the DAPP1 (dual adaptor of phosphotyrosine and 3-phosphoinositides) gene at the 4q23 (Figure 4B) and the gene is suggested to be involved in signal transduction processes.
Another interesting recurrent CNV of size 24.8 kb was detected in two unrelated ASD cases (68672, 50800L), a duplication at 17p13.3 disrupting two exons of YWHAE (tyrosine 3/tryptophan 5-monooxygenase gene) gene, which presumably act via haploinsufficiency (Figure 4C). It was maternally inherited in both ASD cases, and an equivalent event based on our definition of overlap was not present in controls (2/696 cases vs. 0/5139 controls; FET two-tailed p = 0.01421). YWHAE belongs to the 14-3-3 family of proteins, which mediate signal transduction, and is highly conserved in both plants and mammals. Only microduplications in YWHAE gene have been reported in ASD. It has been shown that the phenotype of patients with a 17p13.3 microduplication involving YWHAE gene show autistic manifestation, behavioral symptoms, speech and motor delay, subtle dysmorphic facial features, and subtle hand-foot malformations (Bruno et al. 2010). It is also noteworthy that there were larger CNVs found in two PDx controls to intersect YWHAE. The two ASD cases were both of Asian descent, and we also have found other cases and controls of Asian descent bearing YWHAE-CNVs (M. Gazzellone, unpublished data), suggesting that the role of this gene in ASD need to be further explored.
In three unrelated male ASD probands (44644, 124475, 45554), we observed a recurrent novel CNV, a 24.7-kb duplication encompassing two exons of the SAE1 (SUMO1 activating enzyme subunit 1) gene at the 19q13.32 locus (Figure 4D). The same CNV also was found in one control (3/696 cases vs. 1/5139 controls; FET two-tailed p = 0.00616). Interestingly, another duplication of size 50.8 kb disrupting six exons of SAE1 was observed in a fourth unrelated ASD case (90278) in the present study and was also detected by previous SNP microarray study (Lionel et al. 2011). The SAE1 gene is involved in protein sumoylation process and is shown to interact with the ARX gene, which is involved in Autistic disorder (Sherr 2003; Rual et al. 2005; Ewing et al. 2007; Gareau and Lima 2010; Wilkinson et al. 2010; Szklarczyk et al. 2011).
In two unrelated cases (69180, 59144), we identified overlapping CNVs impacting the PLXNA4 (plexin A4) gene at the 7q32.3 locus (Figure 4E). One CNV is a 14.2-kb loss encompassing an exon of the gene and based on our overlap criteria is not observed in 5,139 controls, while the second CNV is a 15.5-kb gain encompassing untranslated regions of the gene and is observed in 2 of 5139 controls. PLXNA4 is involved in axon guidance as well as nervous system development (Suto et al. 2003; Miyashita et al. 2004).
In the present study, only one de novo CNV was found (all other qPCR-validated CNVs, 116 of 117, were inherited from either parent) a 36-kb loss encompassing the intron of the GPHN (Gephyrin) gene at the 14q23.3 locus. This de novo CNV (Figure 4F) was found in a male ASD proband (103018L) and was not picked up on the previous SNP array (Pinto et al. 2010), and it was not found in any of the 5139 controls. Gephyrin is suggested to play a central organizer role in assembling and stabilizing inhibitory postsynaptic membranes in human brain (Waldvogel et al. 2003). In our other unpublished data, we have also identified another deletion encompassing several exons of GPHN in an unrelated ASD case and a de novo deletion encompassing exons of the gene in a schizophrenia case suggesting that GPHN gene could be a novel susceptibility gene playing a more general role in neurodevelopmental disorders. We believe the lack of novel, rare de novo CNVs captured in the present study is simply due to our study design because nearly all the de novo CNVs in this ASD cohort are relatively larger in size and therefore were already detected using SNP microarrays (Marshall et al. 2008; Pinto et al. 2010, Lionel et al. 2011, A. C. Lionel, unpublished data, i.e., not reported in the present study).
We also detected other novel, rare CNVs present in only one unrelated ASD case (Table 5) in previously identified genes associated with ASD such as CTNND2, CDH18, PARK2, NXPH1, MTHFD1, and NF1 [Table 5 (Williams and Hersh 1998; Marui et al. 2004; Glessner et al. 2009; Pinto et al. 2010; Griswold et al. 2011; Salyakina et al. 2011)]. Novel, rare CNVs occurring in genes known to be associated with other neurodevelopmental disorders (e.g., KIRREL3) or potentially playing a role in neurodevelopment also were found (e.g., CTNNA2, NDST1, SLC24A2, NFIB, APLP2, ATP2C2, CECR2, DAGLA, and UPB1).
A paternally inherited 7.9-kb deletion disrupting an exon of the SLC24A2 [solute carrier family 24 (sodium/potassium/calcium exchanger), member 2] gene was observed at 9p22.1 in a male ASD case (61180-L; Figure 5A) but in none of the controls. The SLC24A2 gene may play a role in neuronal plasticity (Li et al. 2006).
In a female ASD case (118909L), we observed a 7.2-kb loss disrupting an exon of the CECR2 (cat eye syndrome chromosome region, candidate 2) gene at the 22q11.21 locus (Figure 5B), which was not observed in controls. We also detected a CNV in the same gene in ASD case 124632L reported in another study (Pinto et al. 2010). CECR2 is a chromatin remodeling factor that has been proposed to play a role in embryonic nervous system development (Banting et al. 2005).
In another unrelated male ASD proband (72816L), we have identified a 15-kb paternally inherited deletion (Figure 5C) disrupting seven exons of the DAGLA (diacylglycerol lipase, alpha) gene at the 11q12.2 locus that is not found in controls. DAGLA is known to synthesize an endocannabinoid that has been associated with retrograde synaptic signaling and plasticity (Gao et al. 2010).
Another interesting gene is UPB1 (ureidopropionase, beta) located at the 22q11.23 locus, in which a 6.7-kb deletion disrupting two exons has been found in a male ASD proband (154266L; Figure 5D). The deletion was not observed in any of 5139 controls. This gene is involved in the last step of the pyrimidine degradation pathway and deficiencies in UPB1 have been associated with developmental delay (van Kuilenburg et al. 2004).
Finally, our study includes analysis of two ASD cases not previously run on SNP microarrays. In both cases, large-sized CNVs previously associated with ASD were found, a 546-kb maternally inherited 16p11.2 duplication (ASD case 100564) and a 656-kb maternally inherited 22q11.22-q11.23 duplication (ASD case MM0177-3) that overlaps with the 22q distal deletion region (Figure S2).
To perform more robust downstream analyses (global rare CNV burden and gene-set association) on the ASD CNV findings in the present study, we used the SNP genotype data available on the 615 ASD cases to identify those of European ancestry based on MDS. We visualized the HapMap and ASD samples in a two-dimensional space using the first and second MDS components. The plot displayed three major clusters (Europeans, Eastern Asians, and Africans), with Mexicans and Indians distributed between the European and Eastern Asian cluster (Figure S3A). According to MDS components 1 and 2, ASD samples mostly clustered with Europeans, although some of them clearly clustered with other ethnicities. We defined a set of boundaries on the first and second MDS component, which identified all European samples, and only two Mexican samples (Figure S3B). Any ASD sample falling within these boundaries was considered to be European (505 of 615 total genotyped unrelated samples). We did not use the third MDS component, which is better at resolving ethnic subgroups but does not account faithfully for differences among major groups (e.g., Europeans and Eastern Asians were very close to each other; Figure S3C). The rare variants from 505 European ASD cases along with 1000 European PDx controls were used further for rare CNV burden and gene-set association analyses to avoid any confounding issues due to different ancestries of samples.
To investigate global differences in CNV burden, we assessed the distribution of three CNV statistics in ancestry-matched European ASD case subjects compared with control subjects: (1) subject CNV number, (2) subject total CNV length, and (3) subject total number of CNV genes.
We observed significant differences in CNV number and total gene number for deletions but not for duplications (significance threshold: Wilcoxon test p value < 0.01; Table 6). In terms of magnitude, the total number of deletion genes was the largest difference found between the ASD and control subjects (ratio of means: 1.77). Considering that subjects were matched by platform type and other essential parameters, and also considering that previous authors found a stronger difference in deletion burden rather than duplication, the differences observed very likely translate to real biological differences. In addition, the relative difference in total gene number was larger than the relative difference in CNV number, an effect that is harder to explain by technical or experimental factors; because of the relatively large sample size, it is important to consider both the significance and magnitude of burden differences.
We also assessed whether rare CNVs in genes that are causally implicated in ASD were enriched in cases over controls. The ASD gene list used comprised 110 genes compiled from the peer-reviewed literature (Betancur 2011). We observed significant enrichment of deletions impacting genes implicated in ASD [p = 8.896e-05, odds ratio = 20.59 with 95% confidence interval = 2.95-888.96] in ASD cases than controls (Table 7). There were 11 ASD cases with deletions in ASD candidate genes (Table S3) and all were experimentally validated by qPCR except three false-positive CNVs overlapping the SYNGAP1 gene. The one SYNGAP1 CNV that did validate is a de novo 112 kb deletion, which was described previously in the Pinto et al. 2010 study that disrupts SYNGAP1 and encompasses four other genes. After removing the three false-positive CNVs in the SYNGAP1 gene and testing for enrichment again, we still observed significant enrichment of deletions in genes implicated in ASD (p = 0.001585, odds ratio = 14.74 with 95% CI = 1.95−656.52) in ASD cases as compared with controls. The other ASD cases where we observed rare exonic loss in ASD candidate genes were PTCHD1 (Marshall et al. 2008; Pinto et al. 2010), VPS13B (Pinto et al. 2010), DMD (Pinto et al. 2010), DPYD (novel to this study and as described in Carter et al. 2011), SHANK2 (Pinto et al. 2010), NF1 (novel to this study) and NRXN1 (ILMN Omni 2.5 M array; A. C. Lionel, unpublished data).
We tested functional gene-sets for enrichment in ASD cases over controls to identify biological processes potentially involved in ASD. We found significant results only for deletions, which is in agreement with global burden results (Table 6). There were 23 gene-sets that had a permutation-based FDR < 25% in deletions and these were used to construct a functional map of ASD (Figure 6, Table S4). We identified several gene-set clusters, some of which were associated with ASD by earlier studies. For example, gene-sets involved in GTPase/Ras signaling pathways were previously reported in ASD (Meechan et al. 2009; Pinto et al. 2010). The novel pathway that has been discovered from this study is the gene-set enriched for nucleotide metabolism. The list of genes in the nucleotide metabolism pathway is shown in Table S5. To strengthen the conclusion that ASD variants impact genes involved in the nucleotide metabolism pathway, we performed experimental qPCR validation on this set of CNVs and found all were valid except one deletion arising due to cell line artifact in the FHIT gene. The confirmation of 14 of 15 genes (Table S5) based on qPCR validation of the microarray findings suggests that this pathway is not the result of false positives in the dataset.
Our high-resolution CGH data have identified multiple novel, rare CNVs in ASD cases that went undetected by SNP microarrays run on the same individuals; thus providing an additional valuable resource for ASD risk gene discovery and validation. Although hundreds of ASD susceptibility genes/loci (Betancur 2011) are now known, no single gene or locus accounts for more than 0.8% frequency of cases in a given cohort, with most contributing to less than 0.1% of ASD cases (Devlin and Scherer 2012). Moreover, other CNV studies (Pinto et al. 2010; Sanders et al. 2011) and recent exome sequencing studies (Neale et al. 2012; O’Roak et al. 2012; Sanders et al. 2012) suggest that perhaps several hundred autism risk genes may exist. These data indicate that to delineate all ASD risk genes and alleles, different experimental strategies will likely be required when assessing even larger sets of ASD patient collections.
We find that the Agilent CGH array is sensitive for detection of many smaller (<30 kb) CNVs, which are often missed by SNP microarrays even when probe coverage is sufficient at these loci. The results of this study suggest that currently available microarray platforms are complementary (i.e., not all CNVs are captured by one platform/array design) and that the number and type of CNVs detected varies depending on microarray probe distribution, sample labeling and hybridization chemistries, and CNV detection algorithms used. The Agilent 1M array is designed exclusively for CNV detection and utilizes a more uniform probe distribution across the genome as compared with SNP arrays such as the Illumina 1M and Affymetrix 6.0 arrays. The median probe spacing of the Agilent 1M CGH array is 2.1 kb and a small proportion of probes correspond to regions that are non-unique in nature; that is, they map to multiple locations in the genome, with most targeted to segmentally duplicated regions. Therefore, the Agilent 1M CGH array is better at detecting CNVs in segmental duplicated regions compared to SNP arrays. Also, a better signal-to-noise ratio has been observed in CGH arrays compared to SNP arrays and is attributed to the use of longer probes [60-mers (Pinto et al. 2011)]. The SNP arrays are also biased toward known CNVs. However, in addition to copy number analysis, one of the advantages of using SNP arrays is the ability to genotype SNPs which is not possible using the CGH array from this study. Until a single technology (e.g., whole-genome sequencing) is sufficiently robust to capture most genetic variants, including structural variants such as CNVs, use of multiple platforms will be advantageous. In some instances, the identification of previously undetected CNV variants will reveal pathogenic events relevant for clinical diagnosis of ASD (Miller et al. 2010; Scherer and Dawson 2011).
The discovery of additional rare variants in this present study led to the identification of a novel disease pathway of triphosphate nucleotide metabolism; specifically, purine and pyrimidine metabolism as a potential biological process involved in ASD. Several genes found to harbor a rare exonic loss in cases but not in controls have been previously associated to a neuropsychiatric or neurodevelopmental phenotype in human [DPYD, UPB1 (van Kuilenburg et al. 2004, 2009)] and mouse [UPP1, TYMP (Haraguchi et al. 2002; López et al. 2009)] studies. Another interesting functional category revealed by this study is small GTPase signaling pathways, with genes such as ARHGAP15, CDH13, NF1, RALGDS, SYNGAP1, and VAV3 harboring rare losses in cases but not in controls (Table S2). Rare mutations in the ARHGAP15 gene have been reported to be associated with ASD (O’Roak et al. 2011); NF1 and SYNGAP1 genes are also implicated in ASD (Williams and Hersh 1998; Marui et al. 2004; Pinto et al. 2010). The protein encoded by the CDH13 gene is one of numerous cadherins expressed in the brain, which have been shown to regulate many neural processes (Redies et al. 2012) and CDH13 has been implicated in ADHD through genome-wide association and extended pedigree linkage studies (Lesch et al. 2008; Rivero et al. 2012). The RALGDS and VAV3 genes are involved in multiple signaling pathways, including nerve growth factor receptor signaling pathway (Ferro and Trabalzini 2010; Keilhoff et al. 2012). Therefore, the genes CDH13, RALGDS, and VAV3 could represent interesting examples of ASD candidate genes for further follow-up. One of the genes, DLC1, contributed to both the “cytoskeleton and contractile fiber” and “small GTPase signaling” pathways. Rare exonic loss in this gene was observed in an ASD case but not in any of the controls. The DLC1 gene is involved in multiple phenotypes in mouse, including nervous system development (Sabbir et al. 2010) and thus is an interesting ASD candidate gene for further study. Examples of some of the genes contributing to the “cytoskeleton and contractile fiber” pathway are DMD, EPB41L2, MYH7, PGM5, and TRIM32, where we observed rare exonic losses in cases but not in controls. The gene DMD has been shown to be associated with ASD in several studies (Wu et al. 2005; Hinton et al. 2009; Erturk et al. 2010). The TRIM32 gene is suggested to be involved in muscle and nervous system development mouse phenotypes (Kudryashova et al. 2009) and deletions in this gene are also reported in cases with attention deficit hyperactivity disorder (Lionel et al. 2011). The gene EPB41L2 binds glutamate (Shen et al. 2000) and dopamine receptors (Binda et al. 2002) and is also suggested to be involved in mouse reproductive and nervous system development phenotypes (Ivanovic et al. 2012). The PGM5 gene interacts with the DMD gene (Wakayama et al. 2000; Szklarczyk et al. 2011). Thus, we believe this set of ‘Cytoskeleton and Contractile fiber’ pathway genes, which are impacted by ASD-specific variants in the present study, warrant further follow-up.
Whereas much of the focus in identifying ASD genes has been in studying de novo events, assessing the role of both autosomal and X-linked rare inherited CNVs also has been fruitful, yielding new susceptibility loci such as NRXN3 (Vaags et al. 2012), SHANK1 (Sato et al. 2012), and the Xp11.2 PTCHD1-PTCHD1AS region (Noor et al. 2010). However, in taking de novo and rare inherited variants together, no one gene accounts for more than 1% of the etiology in ASD highlighting the complex genetic heterogeneity of the disorder. It is, therefore, essential to capture the entire spectrum of genetic variation contributing to ASD risk to account for unexplained heritability. Moreover, given the rarity of some of the putative risk variants identified in this study, it is likely that high resolution genome-wide scans of tens of thousands of ASD cases will be needed to validate and contextualize these findings. The results from the present study will complement ASD genetic risk factors being identified through current whole exome and genome sequencing efforts. The public availability of the new rich resource of rare CNVs uncovered in this study will serve as an important resource for further prioritization of putative ASD risk genes for subsequent genetic and functional characterization.
We thank Sylvia Lamoureux, Quen Tran, and The Centre for Applied Genomics (TCAG) for technical assistance. This work was supported by grants from NeuroDevNet, the University of Toronto McLaughlin Centre, Genome Canada, and the Ontario Genomics Institute, the Canadian Institutes for Health Research (CIHR), the Canadian Institute for Advanced Research, the Canada Foundation for Innovation, the government of Ontario, and The Hospital for Sick Children Foundation. A.P. is supported by a fellowship from the CIHR Autism Training Program. S.W.S. holds the GlaxoSmithKline-CIHR Chair in Genome Sciences at the University of Toronto and The Hospital for Sick Children. Two provisional patent applications originating from this data have been filed with the U.S. Patent and Trademark Office.
Communicating editor: R. Cantor