|Home | About | Journals | Submit | Contact Us | Français|
Autism spectrum disorder is a severe early onset neurodevelopmental disorder with high heritability but significant heterogeneity. Traditional genome-wide approaches to test for an association of common variants with autism susceptibility risk have met with limited success. However, novel methods to identify moderate risk alleles in attainable sample sizes are now gaining momentum.
In this study, we utilized publically available genome-wide association study data from the Autism Genome Project and annotated the results (P <0.001) for expression quantitative trait loci present in the parietal lobe (GSE35977), cerebellum (GSE35974) and lymphoblastoid cell lines (GSE7761). We then performed a test of enrichment by comparing these results to simulated data conditioned on minor allele frequency to generate an empirical P-value indicating statistically significant enrichment of expression quantitative trait loci in top results from the autism genome-wide association study.
Our findings show a global enrichment of brain expression quantitative trait loci, but not lymphoblastoid cell line expression quantitative trait loci, among top single nucleotide polymorphisms from an autism genome-wide association study. Additionally, the data implicates individual genes SLC25A12, PANX1 and PANX2 as well as pathways previously implicated in autism.
These findings provide supportive rationale for the use of annotation-based approaches to genome-wide association studies.
Autism spectrum disorder (ASD) is characterized by impairments in three core domains: development of social interaction; language and communication development; and restricted or circumscribed interests and repetitive behaviors . By definition, children with autism show impairments in all three areas, though the degree to which each is affected varies. The study of autism genetics began with twin studies conducted as early as the 1970s , indicating a strong heritable component. Since then, the autism genetics field has progressed to include linkage, association, copy number variation and candidate gene studies targeting the broad autism diagnosis as well as specific endophenotypes of autism. The discovery of genetic variants involved in autism has been confounded by phenotypic, genetic and allelic heterogeneity, multiple modes of transmission, inconsistent application of diagnostic criteria and unknown environmental risks. Progress has been made in our understanding of the genetic architecture of autism despite these complexities. However, few genes harboring common risk factors have been convincingly implicated in the etiology of autism.
Multiple traditional genome-wide association studies (GWAS) of ASD have been published with limited power to detect common variants contributing to risk. Genome-wide significantly associated risk alleles have been identified on chromosome 5p14 [3,4], 5p15  and 20p12.1 . The latter study was the largest to date and included over 1,500 parent–child-trio families in the analysis. The single genome-wide significant result found on chromosome 20p12.1 fell within the MACRO domain containing protein 2 gene. However, at least one subsequent study has failed to replicate this association .
Recent studies have shown that novel analytic approaches using GWAS data can increase the power to identify statistically significant relationships. For example, Lu and Cantor describe a method that allows for sex-specific associations, resulting in highly significant findings at the ryanodine receptor 2 locus (P <3.9×10-11) in a multiplex autism dataset . Additionally, several pathway-based analyses of autism GWAS have recently been published, showing enrichment of sub-genome-wide significant associations with pathways of prior interest in autism. For example, Yaspan and colleagues applied their method (pathway analysis by randomization incorporating structure) to the available Autism Genetic Resource Exchange GWAS dataset, subsequently identifying over-representation of pathways involved in ubiquitination (supporting a pathway analysis in Glessner et al. ) as well as synthesis and degradation of ketone bodies . Anney and colleagues also published a pathway-based analysis, which successfully applied the SNP ratio test methodology to genome-wide data from Autism Genome Project (AGP) families, resulting in a number of convincingly implicated pathways . Mining of existing autism GWAS data is proving to be a valuable undertaking with the potential to provide significant insight into the biology and common risk assessment of ASD.
Recent studies providing thorough investigation of top signals from GWAS of a number of complex disorders have shown enrichment for expression quantitative trait loci (eQTL) in relevant tissues [11-13]. These findings suggest that SNPs with significant functional impact on the expression of genes are more likely to be associated with disease phenotypes . Based on this finding, our group sought to determine if brain-specific eQTL were over-represented in the most significant association from a recently published autism GWAS and whether the specific genes implicated as enriched might highlight new, or help cement existing, candidate genes and pathways. We used publically available data from the AGP GWAS , the SNP and CNV Annotation Database  and genome-wide expression datasets in brain  to test for enrichment of eQTL (those previously identified in transformed lymphocytes and two regions of the brain, parietal lobe and cerebellum), among top nominally associated SNPs from reported autism GWAS.
This study was designed to detect enrichment of brain eQTL amongst top signals from the recent AGP GWAS by Anney et al. . All human subjects research in the study by Anney et al. was conducted in accordance with the Helsinki Declaration and with oversight by the appropriate site-specific Internal Review Board . A schematic of the method utilized in this study is presented in FigureFigure1.1. Briefly, four total lists of SNPs were collected from the additional files of the paper by Anney et al. These include one SNP list from each primary analysis (Spectrum diagnosis [Spec]|All Ancestry [All]; Strict diagnosis [Strict]|All; Spec|Western European Ancestry [WestEur]; Strict|WestEur) meeting a GWAS P-value threshold of 10-3. Three of these analyses reflect strata of the broadest diagnostic (spectrum) and ancestry (all) categories and therefore these analyses are overlapping and non-independent. The top SNPs from each of these GWAS were analyzed for eQTL enrichment in two different brain tissues (parietal and cerebellum) as well as transformed lymphoblastoid cell lines (LCLs), as described below. Comparison of each SNP list to a randomized null distribution of SNPs (described below) within each eQTL database yielded an empirically determined enrichment P-value. Furthermore, the targets for those eQTL present within the four AGP GWAS lists and showing enrichment in brain tissues were identified. Additional analyses were conducted on these eQTL targets to determine if particular genes, gene families or pathways were strongly implicated by eQTL present in any of the primary AGP GWAS analyses.
The methods and samples used for the AGP GWAS are described in depth in the publication by Anney et al. . Additional files published with that analysis contained all SNPs significant at P <0.001 under an additive model from which the primary analyses were extracted. The lists of significant SNPs came from family-based primary analyses including those stratified by diagnosis subtype and population. In total, four lists of SNPs, significant at a GWAS P <10-3 threshold, were taken from the publically available GWAS results published by Anney et al.  (Table(Table11).
eQTL data from LCLs were generated from 176 HapMap European (CEU) and African (YRI) cell lines as described previously (GSE7761) . Cerebellar (GSE35974) and parietal cortex (GSE35977) cis and trans eQTL were generated from 153 individuals of European ancestry obtained from the Stanley Medical Research Institute . ComBat and Surrogate Variable Analysis were used to adjust for batch and both known and unknown covariate effects [16,17]. SNP genotyping was performed on the Affymetrix Genome-Wide Human SNP Array 5.0 (Affymetrix, Santa Clara, CA, USA) and imputed, using MaCH to include all HapMap SNPs [18,19]. The Affymetrix Human Gene 1.0 ST Exon Array was used for gene expression profiling. Imputed genotype dosage data were analyzed for association with expression phenotypes using PLINK . For cis eQTL, defined as SNPs within 4Mb of the probe site, the corrected significance threshold for determining eQTL status (eQTL P-value) was P <0.0001. For trans eQTL, the eQTL P-value was corrected as 0.05/n, where n=25,834 probes for the expression data.
All ASD GWAS SNPs with a GWAS P ≤10-3 in the four primary analyses were subsequently annotated with eQTL information, including the strength of the evidence for the impact of the polymorphism on expression. To test for enrichment of eQTL among these top SNP associations, 1,000 randomized SNP sets were generated. Each of these 1,000 lists, comprising the simulated set of SNPs, included the same number of SNPs as the original list of ASD GWAS associations and contained only SNPs matching the minor allele frequency distribution of the original list of ASD GWAS SNPs. The random lists were sampled without replacement from the set of typed SNPs on the Illumina HumanHap610 array. Minor allele frequency matching was conducted by classifying all Illumina Hap610 platform SNPs into discrete minor allele frequency bins at 5% intervals (0 to 5%, 5 to 10%, …,45 to 50%), followed by random selection of SNPs from the same allele-frequency bins as those in the top signals. The number of eQTL in each simulated set yields an empirical distribution and an enrichment P-value, calculated as the proportion of randomized sets in which the eQTL count matches or exceeds the actual observed count in the list of top SNP associations.
The program PANTHER (http://www.pantherdb.org) was used to identify networks or gene families represented by top autism GWAS eQTL. PANTHER’s classification of genes and proteins provides functionality to determine if gene annotations (that is, membership in pathway or Gene Ontology) are over-represented in a test gene list compared with expectation, given all genes in the human genome. P-values generated by the program are Bonferonni-corrected for the number of comparisons conducted. Of the 140 unique genes implicated as eQTL targets from the four primary autism GWAS, 130 were mapped to identifiers in the PANTHER database and analyzed for pathway, molecular function and biological process terms that were over-represented compared to expectation.
The focus of this study was to determine if enrichment of brain eQTL (that is, SNPs that regulate brain expressed genes) existed in the top signals from published autism GWASs. To this end, we analyzed SNPs at a GWAS P <10-3 for four primary analyses from the published AGP GWAS .
The primary GWAS enrichment analyses consist of two stratifications including diagnosis and ancestry. We observed significant enrichment of parietal (P=0.004) and cerebellar (P <0.003) eQTL, but not LCL eQTL (P=0.502) among the top signals from the most broadly inclusive dataset of spectrum diagnosis including all ancestries (Figure(Figure2).2). Similarly, enrichment of both parietal (P=0.005) and cerebellar (P=0.012) eQTL persisted in the spectrum diagnosis when limited to individuals of Western European descent. Significant enrichment of parietal eQTL (P <0.001), but not cerebellar (P=0.872) or LCL eQTL (P=0.110) was identified in the Strict|All dataset (Figure(Figure3).3). We again observed significant enrichment of parietal (P <0.001) eQTL but not cerebellar (P=0.45) or LCL (P=0.49) eQTL in the Strict diagnosis when limited to individuals of Western European descent. Results from all analyses are summarized in TableTable11.
SNPs from all four primary non-independent SNP lists that were eQTL in either cerebellum or parietal tissues contributed to a total number of 539 brain eQTL found in the primary AGP GWAS top signals. In brief, of the total 539 eQTL SNPs, 256 eQTL SNPs were uniquely identified, meaning that the remaining 283 eQTL SNPs were redundancies of the 256 unique eQTL SNPs. Of the unique 256 eQTL SNPs, 214 (84%) were found to be eQTL for only one gene. Forty-two (16%) eQTL SNPs acted as eQTL for more than one gene. Of the 214 SNPs that target only one gene, 124 (58%) act in cis and 88 (42%) act in trans. Of the 42 SNPs that target more than one gene, 21 (50%) SNPs regulate all target genes in cis, 5 (12%) SNPs regulate all target genes in trans, and 16 (38%) SNPs are both trans and cis regulators of different target genes. All eQTL from enrichment analysis are provided in Additional file 1: Table S1.
From this total set of SNPs, 140 genes were uniquely implicated as eQTL targets. Of these 140 genes, 62 (44%) were cis implicated and 78 (56%) implicated in trans; the proportion of cis implicated genes (in contrast to trans implicated genes) is greater than expected in either cerebellum or parietal tissues based on the observed number of such genes (31% and 28%, respectively). Of the 140 uniquely implicated genes, 71 (51%) genes were implicated in only one analysis; 69 (49%) genes were detected in two or more analyses; 20 (14%) genes were identified in three or more analyses; and 14 (10%) genes were found in all four analyses. There was relatively little overlap in the eQTL target genes across tissues. In an analysis of the broadest diagnostic and ancestry category, only 18 (13%) of the 140 brain-enriched genes were also targeted by eQTL found in LCLs. Of these 18, only three (2%) were identified in parietal, cerebellum and LCL tissues. Of the 140 brain-enriched genes, only 10 genes (7%) were found in both cerebellum and parietal tissues. In all pairwise comparisons of tissues, the overlap in the eQTL target genes was not statistically significant (P >0.05, Fisher’s exact test). All genes implicated in enrichment analyses are given in Additional file 2: Table S2.
We reviewed the function of those genes implicated in the enrichment analysis to identify common pathways, gene families or functions of said genes. The genes solute carrier family 25 member 12 (SLC25A12), pannexin 1 (PANX1) and pannexin 2 (PANX2) were strongly implicated in our analyses. SLC25A12 was multiply-implicated by a unique set of 31 SNPs, which was by far the largest number of SNPs implicating a single transcript and was the only transcript implicated in such a way (Figure(Figure4).4). These SNPs all show a high degree of linkage disequilibrium and appear to implicate an extended haplotype (Figure(Figure5).5). Twenty-eight of the 31 SNPs were found in more than one analysis and two of those 28 SNPs were found in three analyses (Figure(Figure6).6). All SNPs implicating SLC25A12 were cis regulators. All eQTL SNPs implicating SLC25A12, and additional information such as SNPs in linkage disequilibrium (r2>0.8), are provided in Additional file 3: Table S3.The genes PANX1 and PANX2 were also multiply-implicated. PANX1 was identified as an eQTL target in all four of the primary analyses. No genes were implicated by both cis and trans eQTL; however, PANX1 was targeted by nine cis eQTL and PANX2 was targeted by three trans eQTL.
Additionally, we sought to identify any overlap between the eQTL implicated in this study and those recently published by Voineagu and colleagues , who identified a number of genes differentially expressed in autism and matched control brains. We found that four of the 140 brain eQTL target genes from our list overlapped with the 1,153 differentially expressed genes identified by Voineagu et al. . These included SLC25A12PANX2calponin 3 acidic and cleavage and polyadenylation specificity factor subunit 1.
Based on PANTHER analysis, we found modest evidence of significant over-representation of 14 processes and pathways of interest in autism, including structural components of ribosomes (P <0.009), methyltransferase activity (P<0.05), the nicotinic acetylcholine receptor signaling pathway (P<0.02), sensory perception (P<0.01), and receptor activity (P<0.05) among others (Table(Table22).
This study sought to identify biologically meaningful associations buried within top signals that did not meet genome-wide significance thresholds from autism GWAS results. Traditional GWAS in autism have been burdened with problems of heterogeneity and limited sample sizes. To date, multiple traditional autism GWAS have been published although it is important to note that some of these studies use overlapping datasets [3-6]. We believe that, while these studies might be underpowered to identify signals that survive genome-wide correction, their reported results do not reflect the full value of the data.
There are excellent examples to date of such novel methodologies providing additional insight into the role of common genetic variants in complex human disorders, including autism. Anney and colleagues in the AGP also conducted a novel analysis of the present genome-wide data utilizing the SNP ratio test in an effort to detect over-transmitted SNPs in predefined pathways of interest, which subsequently identified multiple pathways and/or Gene Ontology terms enriched for associated signals . Interestingly, two of these pathways, ribosomal components and methyltransferase activity, were also identified in the current study as pathways regulated by eQTL-associated SNPs. Our study complements these pathway-based approaches and demonstrates a method of utilizing expression based SNP annotations to mine already existing analyses for disease associations with functional variants in relevant tissues.
Our study identified significant enrichment of parietal and cerebellum eQTL, but not LCL eQTL, among the top SNP associations in all four of the primary AGP GWAS identified by Anney et al. This pattern of results showing significant enrichment in brain (the affected and relevant tissue) and not in tissues peripheral to the main pathology was also seen in a study of cis-regulatory SNPs in bipolar disorder . Similarly, Below et al.  report enrichment of top signals from a type 2 diabetes GWAS in tissues involved in pathogenesis (muscle and adipose) but not LCLs. Additionally, we found no significant difference in proximity to genes between our GWAS implicated SNPs and the SNPs forming the null distribution. Taken together, the minor allele matching, the lack of a significant difference in distance to nearest gene, and the fact that enrichment is found in some tissues and not in others provide strong evidence against possible sources of systematic bias.
There was very little overlap in the eQTL targets found in parietal, cerebellum and LCLs. Only three target genes, protein NIP-SNAP homologue 3A (NIPSNAP3A), TM2 domain containing protein 2 (TM2D2) and copine 1 (CPNE1), were present in all three tissues in the analysis of the broadest diagnostic and ancestry categories. These findings also strongly suggest that the most appropriate tissue for further functional work is the tissue of pathology. This poses clear challenges for difficult to access tissues, such as neurons. However, there is significant evidence that human induced pluripotent stem cells, once transformed into neurons, provide an invaluable source of neuronal tissue from live patients in which to study cellular phenotypes [24-26]. Our study highlights the importance of furthering the development of methods that allow access to tissues involved in pathology.
While we expected, and observed, overlap in top signals between each primary GWAS, our findings showed an abundance of unique eQTL and implicated a total of 140 target genes. Additionally, three genes, SLC25A12PANX1 and PANX2, were strongly implicated by our results. Two of these genes, SLC25A12 and PANX2, were also identified in a recent study of genes differentially expressed in the brains of individuals with autism compared to control brains . PANX2 is a protein implicated by trans eQTL that modulates the timing of neuroprogenitor commitment to a neuronal lineage in the hippocampus . It is located 500kb proximal to SH3 and multiple ankyrin repeat domains protein 3 in the Phelan-McDermid Syndrome 22q13.33 terminal deletion region that has been implicated in ASD . PANX1 was implicated in cis by multiple SNPs in multiple analyses and across multiple tissues and has recently been shown to play a role in N-methyl-D-aspartic acid-mediated epileptiform activity . In murine models, Panx1 and Panx2 are strongly co-expressed early in the developing brain, specifically in the hippocampus .
SLC25A12 is a 109kb gene found on chromosome 2q31.1 and functions as a calcium-binding mitochondrial protein, integral in the exchange of aspartate for glutamate across the mitochondrial membrane. Several studies have investigated the role of SLC25A12 in the development of autism and SLC25A12 SNPs have been implicated [31-34]. Due to the platform used in the original AGP GWAS and quality-based filtering of imputed eQTL data, only one of the previously published SLC25A12 risk variants (rs908670) was included in our study . Rs908670 was identified as an SLC25A12 eQTL in parietal tissue. Previous work has found SLC25A12 to be differentially expressed in autism brains compared to controls [22,35]. One such recent study found SLC25A12 expression decreased in autism brains compared to controls . We found that the risk alleles of the SLC25A12 eQTL SNPs identified in the AGP GWAS are correlated with decreased expression of SLC25A12 in the parietal lobe. We add our findings to this growing body of evidence and suggest that common variation in SLC25A12 auto-regulates expression and may contribute to autism susceptibility.
Our study is limited by the available analyses and the sample sizes used in the original GWAS. It is with some caution that the comparison of results (for example, Strict|All versus Spec|All) must be interpreted. Enrichment analyses, such as those conducted here, include only the top SNP signals and therefore are sensitive to power and sample size differences between the GWAS scans that provide the data to go into an enrichment analysis. It is not possible to make direct comparisons between GWAS enrichment results that differed with respect to sample size, ranging from 720 subjects (Strict|WestEur) to 1,385 subjects (Spec|All). However, within a single sample (for example, Spec|All or Strict|WestEur), comparisons can easily be made between all tissues including cerebellum, parietal and LCLs.
Additional mining of available GWAS data may provide new insight into the biology of autism while allowing the genetics community to leverage data from smaller studies of GWAS. Our study provides evidence for the hypothesis that SNPs below the genome-wide significant threshold (P=10-8) are functionally relevant to the development of autism and may yet contribute to risk. Additionally, our results point specifically to SLC25A12 and PANX1/2 and to pathways (such as, methyltransferase, ribosomal components) previously implicated in autism.
Our results demonstrate that loci moderately associated with autism are more likely to be brain-specific eQTL than randomly sampled loci conditioned on minor allele frequency. Moreover, targets for the eQTL identified among top signals from autism GWAS include compelling candidate genes such as SLC25A12 and the pannexins. These findings suggest that while very few SNPs from a GWAS reach genome-wide significance, it is likely that an abundance of SNPs at sub-significant thresholds are indeed associated with disease risk. This study provides further rationale for annotation-based approaches to GWAS methods.
ASD: autism spectrum disorder; AGP: Autism Genome Project; eQTL: expression quantitative trait loci; GWAS: genome-wide association study; LCL: lymphoblastoid cell line; PANX1: pannexin 1; PANX2: pannexin 2; SNP: single nucleotide polymorphism; SLC25A12: solute carrier family 25 member 12.
The authors declare they have no competing interests.
LKD conducted analyses and drafted the manuscript. EG conducted analyses and assisted in drafting the manuscript. EKG assisted in drafting the manuscript. JAB and CL conducted analyses yielding brain eQTL data and participated in the drafting of the manuscript. EHC directed analyses and assisted in drafting the manuscript. JS directed analyses of the study and assisted in drafting the manuscript. NJC conceived of the study, directed analyses and helped to draft the manuscript. All authors read and approved of the final manuscript.
Table S1. Includes all eQTL from enrichment analysis identified in parietal and cerebellar tissues, from all of the four primary GWAS analyses. Worksheets provided in the database include information on the eQTL target gene as well as the eQTL host gene.
Table S2. Provides additional information on each gene that was targeted in cerebellum and parietal tissues by an autism-associated eQTL SNP.
Table S3. Provides additional supplementary data regarding SLC25A12 eQTL including the original GWAS p-values and z-scores, parietal tissue eQTL p-values and beta scores and tagged and proxy SNPs for significant eQTL.
We wish to acknowledge the efforts made by the families who participated in the original AGP GWAS. We also wish to acknowledge the AGP-GPAC for their collection of the sample, dissemination of results and thoughtful discussion of the methods and results in our paper. The project described was supported in part by Grant Number by Grant Number R01 MH090937, P50MH094267, and by Grant Number P50 HD055751 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development. Finally we would like to thank Suma Jacob, MD, PhD, for her review of the manuscript.