Transcriptomic analysis of autistic brain reveals convergent molecular pathology.

Autism spectrum disorder (ASD) is a common, highly heritable neurodevelopmental condition characterized by marked genetic heterogeneity. Thus, a fundamental question is whether autism represents an aetiologically heterogeneous disorder in which the myriad genetic or environmental risk factors perturb common underlying molecular pathways in the brain. Here, we demonstrate consistent differences in transcriptome organization between autistic and normal brain by gene co-expression network analysis. Remarkably, regional patterns of gene expression that typically distinguish frontal and temporal cortex are significantly attenuated in the ASD brain, suggesting abnormalities in cortical patterning. We further identify discrete modules of co-expressed genes associated with autism: a neuronal module enriched for known autism susceptibility genes, including the neuronal specific splicing factor A2BP1 (also known as FOX1), and a module enriched for immune genes and glial markers. Using high-throughput RNA sequencing we demonstrate dysregulated splicing of A2BP1-dependent alternative exons in the ASD brain. Moreover, using a published autism genome-wide association study (GWAS) data set, we show that the neuronal module is enriched for genetically associated variants, providing independent support for the causal involvement of these genes in autism. In contrast, the immune-glial module showed no enrichment for autism GWAS signals, indicating a non-genetic aetiology for this process. Collectively, our results provide strong evidence for convergent molecular abnormalities in ASD, and implicate transcriptional and splicing dysregulation as underlying mechanisms of neuronal dysfunction in this disorder.

We analyzed post-mortem brain tissue samples from 19 autism cases and 17 controls from the Autism Tissue Project and the Harvard brain bank (Supplementary Table1) using Illumina microarrays. For each individual, we profiled three regions previously implicated in autism[5]: superior temporal gyrus (STG /Brodmann’s area [BA] 41/42), prefrontal cortex (BA9) and cerebellar vermis. After filtering for high-quality array data (Methods), we retained 58 cortex samples (29 autism, 29 controls) and 21 cerebellum samples (11 autism, 10 controls) for further analysis (see Methods for detailed sample description). We identified 444 genes showing significant expression changes in autism cortex samples (DS1, Figure 1b), and only 2 genes differentially expressed (DE) between the autism and control groups in cerebellum (Methods), indicating that gene expression changes associated with autism were more pronounced in cerebral cortex, which became the focus of further analysis (Supplementary Table 2). There was no significant difference in age, post mortem interval (PMI), or RNA integrity numbers (RIN) between autism and control cortex samples (Supplementary Figure 1, Methods).

Figure 1

Gene expression changes in autism cerebral cortex

A. Heatmap of top 200 genes differentially expressed between autism and control cortex samples. Scaled expression values are color-coded according to the legend on the left. The dendrogram depicts hierarchical clustering based on the top 200 DE genes. The top bar (A/C) indicates the disease status: red-autism, black-control. The bottom bars show additional variables for each sample: sex (grey-male, black-female), brain area (black-temporal, grey-frontal), co-morbidity of seizures (green-autism case with seizure disorder, red-autism case without seizure disorder, black-control), age, RNA integrity number (RIN) and post mortem interval (PMI). The corresponding scale for quantitative variables is shown on the left.

B. Top: Venn diagram depicting the overlap between genes DE in frontal and temporal cortex. Bottom: Venn diagram describing the overlap between genes DE in the initial cohort (DS1) and the replication cohort (DS2). Differential expression in the initial cohort was assessed at an FDR<0.05 and fold changes>1.3. The statistical criteria were relaxed to p<0.05 for the replication dataset since it involved fewer samples.

C. Expression fold changes for all genes DE in the initial cohort are plotted on the x-axis against the fold changes for the same genes in the replication cohort on the y-axis. Green-genes downregulated in the autism group in both datasets, red-genes upregulated in the autism group in both datasets, grey-genes with opposite direction of variation in the two datasets. Horizontal lines show fold change threshold for significance.

D. (left) Diagram depicting the number of genes showing significant expression differences between frontal and temporal cortex in control samples (top) and autism samples (bottom) at FDR<0.05. (right) Top 20 genes differentially expressed between frontal and temporal cortex in control samples. All of the genes shown are also differentially expressed between frontal and temporal cortex in fetal midgestation brain[10], but show no significant expression differences between frontal and temporal cortex in autism. The horizontal bars depict p values for differential expression between frontal and temporal cortex in the autism and control groups.

Supervised hierarchical clustering based on the top 200 DE genes showed distinct clustering of the majority of autism cortex samples (Figure 1a), including one case that was simultaneously found to have a 15q duplication (Methods, Supplementary Table 1), which is known to cause 1% of ASD[6]. Cortex samples from ten of the cases coalesced in a single tight-clustering branch of the dendrogram. Clustering was independent of age, sex, RIN, PMI, co-morbidity of seizures, or medication (Figure 1a; Supplementary Figure 2c). It is interesting to note that the 2 ASD cases that cluster with controls (Figure 1a) are the least severe cases, as assessed by global functioning (Supplementary Table 12). We observed a highly significant overlap between DE genes in frontal and temporal cortex (p = 10E-44; Figure 1b), supporting the robustness of the data and indicating that the autism-specific expression changes are consistent across these cortical areas. We also validated a cross section of the DE genes by RT-PCR and confirmed microarray-predicted changes in 83% of the genes tested (Methods, Supplementary Figure 2b). Gene ontology (GO) enrichment analysis (Methods) revealed that the 209 genes downregulated in autistic cortex were enriched for GO categories related to synaptic function, whereas the upregulated genes (N=235) showed enrichment for GO categories implicated in immune and inflammatory response (Supplementary Table 3).

To test whether these findings were replicable, and to further validate the results in an independent dataset, we obtained tissue from an additional frontal cortex region (BA44/45) from nine ASD cases and five controls (DS2; Supplementary Table 4). Three of the cases and all of the controls used for validation were independent from our initial cohort. Ninety-seven genes were differentially expressed in BA44/45 in DS2, and 81 of these were also DE in our initial cohort (p=1.2E-93, hypergeometric test; Figure 1b, c). Remarkably, the direction of expression differences between autism and controls was the same as in the initial cohort for all but 2 of the 81 overlapping DE probes. Hierarchical clustering of DS2 samples based on either the top 200 genes DE in the initial cohort or the 81 overlapping genes showed distinct separation of cases from controls (Supplementary Figure 6). In addition, comparison of these DE results with another, smaller study of the STG in ASD[7], revealed significant consistency at the level of DE genes, including down-regulation of DLX1 and AHI1 (Supplementary Table 5). Thus, DE analysis produced robust and highly reproducible results, warranting further refined analysis.

We next applied weighted-gene co-expression network analysis (WGCNA)[8; 9] to integrate the expression differences observed between autistic and control cerebral cortex into a higher order, systems level context. We first asked whether there are global differences in the organization of the brain transcriptome between autistic and control brain by constructing separate co-expression networks for the autism and control groups (Methods). The control brain network showed high similarity with the previously described human brain co-expression networks (Supplementary Table 7), consistent with the existence of robust modules of co-expressed genes related to specific cell types and biological functions[8]. Similarly, the majority of the autism modules (87%) showed significant overlap with the previously described human brain modules (Supplementary Table 6), indicating that many features reflecting the general organization of the autism brain transcriptome are consistent with that of the normal human brain.

The expression levels of each module were summarized by the first principal component (the module eigengene), and were used to assess whether modules are related to clinical phenotypes or other experimental variables, such as brain region. Two of the control module eigengenes (cM6, cM13) showed significant differences (p<0.05) between the two cortical regions as expected, whereas none of the ASD modules showed any differences between frontal and temporal cortex. This led us to explore the hypothesis that the normal molecular distinctions between the two cortical regions tested were altered in ASD compared with controls. Remarkably, while 174 genes were DE between control BA9 and BA41 (FDR<1%), none of the genes were DE in the same regional comparison among the ASD cases. This was not simply an issue of statistical thresholds, as relaxing the statistical criteria for DE to an FDR of 5%, identified over 500 DE genes in controls, and only 8 in ASD brains, confirming the large difference observed in regional cortical differential gene expression between ASD cases and controls (Figure 1d, Methods). Analysis of DE from a data set[10] of gene expression in developing fetal human brain showed a highly significant (p = 5.8E-9) overlap of DE genes with those found in controls in this study, independently confirming that these genes differentiate normal temporal and frontal lobes. We evaluated the homogeneity of gene expression variance across the autism and control groups using Bartlett’s test (Methods) which indicated that increased variance was not the major factor responsible for the striking difference in regional gene expression between ASD and controls (Supplementary Figure 7 and Supplementary Data).

These data suggest that typical regional differences, many of which are observed during fetal development[10], are attenuated in frontal and temporal lobe in autism brain, pointing to abnormal developmental patterning as a potential pathophsyiological driver in ASD. This is especially interesting in light of a recent anatomical study of 5 cases with adult autism which demonstrated a reduction in typical ultrastructural differences between 3 frontal cortical regions in autism[11]. Together, these independent studies provide both molecular and structural evidence suggesting a relative diminution of cortical regional identity in autism.

To identify discrete groups of co-expressed genes showing transcriptional differences between autism and controls, we constructed a co-expression network using the entire dataset, composed of both autism and control samples (Methods). As previously shown for complex diseases[12; 13] co-expression networks allow analysis of gene expression variation related to multiple disease-related and genetic traits. We assessed module eigengene relationship to autism disease status, age, gender, cause of death, co-morbidity of seizures, family history of psychiatric disease, and medication, providing a complementary assessment of these potential confounders to that performed in the standard DE analysis (Supplementary Table 9).

The comparison between autism and control groups revealed two network modules whose eigengenes were highly correlated with disease status, and not any of the potential confounding variables (Supplementary Table 9). We found that the top module (M12) showed highly significant enrichment for neuronal markers (Supplementary Table 9), and high overlap with two neuronal modules previously identified as part of the human brain transcriptional network[8]: a PVALB+ interneuron module and a module of genes involved in synaptic function. The M12 eigengene was under-expressed in autism cases, indicating that genes in this module were down-regulated in the autistic brain (Figure 2). Consistent with the pathways identified to be downregulated in autism by DE analysis (Supplementary Table 3), the functional enrichment of M12 included the GO categories involved in synaptic function, vesicular transport and neuronal projection.

Figure 2

Gene co-expression modules associated with autism

A, D Top-heatmap of genes belonging to the co-expression module. Bottom-corresponding module eigengene values (y-axis) across samples (x-axis) . Red-autism, grey-controls.

B, E Visualization of the M12 and M16 modules respectively. The top 150 connections are shown for the each module. Genes with the highest correlation with the module eigengene value (i.e intramodular hubs) are shown in larger size.

C, F Relevant gene ontology categories enriched in the M12 and M16 modules

Remarkably, unlike DE genes, M12 showed significant overrepresentation of known autism susceptibility genes[2] (Supplementary Table 10; p = 6.1E-4), including CADPS2, AHI1, CNTNAP2, and SLC25A12, supporting the increased power of the network-based approach to identify disease-relevant transcriptional changes. A further advantage of network analysis over standard analysis of DE is that it allows one to infer the functional relevance of genes based on their network position[9]. The hubs of M12, i.e. the genes with the highest rank of M12 membership[8], were A2BP1, APBA2, SCAMP5, CNTNAP1, KLC2, CHRM1 (Supplementary Data). The first three of these genes have previously been implicated in autism[14-16], while the fourth is a homologue of the autism susceptibility gene CNTNAP2[17]. We highlight the group of genes most strongly connected to the known ASD genes (Supplementary Figure 5) and emphasize the down-regulation of several interneuron markers, such as DLX1 and PVALB, as candidates for future genetic and pathologic investigations.

The second module of co-expressed genes highly related to autism disease status, M16, was enriched for astrocyte markers and markers of activated microglia (Supplementary Table 9), as well as for genes belonging to immune and inflammatory GO categories (Figure 2). This module, which was upregulated in ASD brain, showed significant similarity to two modules identified in previous studies of normal human brain[8]: an astrocyte module and a microglial module. Consistent with this functional annotation, two of the hubs of the M16 module were known astrocyte markers (ADFP and IFITM2).

One of the hubs of the M12 module was A2BP1/FOX1, a neural- and muscle-specific alternative splicing regulator[18] and the only splicing factor previously implicated in ASD[16]. Since A2BP1 was down-regulated in several ASD cases (Supplementary Figure 8), this observation provided a unique opportunity to identify potential disease relevant A2BP1 targets. While A2BP1-regulated alternative exons have been predicted genome-wide[19], few genes have been experimentally validated as A2BP1 targets[20]. To identify potential A2BP1-dependent differential splicing (DS) events in ASD brain, we performed high-throughput RNA sequencing (RNA-Seq) on 3 autism samples with significant down-regulation of A2BP1 (average fold change by qRT-PCR: 5.9) and 3 control samples with average A2BP1 levels. We identified 212 significant alternative splicing events (Supplementary Data). Among these, 36 had been defined by Zhang et al.[19] as predicted targets of A2BP1/2, which represents a highly significant overlap (p=2.2E-16). In addition, five previously validated A2BP1 targets showed evidence of alternative splicing, four of which (ATP5C1, ATP2B1, GRIN1 and MEF2C) were confirmed as having differential splicing between ASD samples with low A2BP1 expression and control samples, indicating that we were able to identify a high proportion of the expected A2BP1-dependent DS events. We also observe that alternative exons with increased skipping in ASD relative to control cases are significantly enriched for A2BP1 motifs in adjacent, downstream intronic sequences (p = 1.09E-7, Fisher exact test), consistent with previous data[19].

The top GO categories enriched among ASD DS genes highly overlapped with the GO categories found to be enriched in the M12 module (Figure 3b). In addition, A2BP1 target genes showed enrichment for actin-binding proteins and genes involved in cytoskeleton re-organization (Figure 3b). Among top predicted A2BP1-dependent DS events (Figure 3a) are CAMK2G, which also belongs to the M12 module, as well as NRCAM and GRIN1. The latter are proteins involved in synaptogenesis, in which allelic variants have been associated with autism and schizophrenia respectively[21; 22].

Figure 3

A2BP1-dependent differential splicing events

A. Top A2BP1-specific differential splicing events. DS events showing the most significant differences in alternative splicing between low-A2BP1 autism cases and controls as well as DS differences consistent with the A2BP1 binding site position. The horizontal axis depicts the percentage of transcripts including the alternative exon. Red-autism samples, black-control samples.

B. Relevant gene ontology categories enriched in the set of genes containing exons differentially spliced between low-A2BP1 autism cases and controls.

RT-PCR assays confirmed a high proportion (85%) of the tested DS changes involving predicted A2BP1 targets (Supplementary Figure 8). We further tested the DS events validated by RT-PCR in 3 independent ASD cases with decreased A2BP1 levels and confirmed the predicted changes in alternative splicing (Supplementary Figure 8), indicating that the observed DS events are indeed associated with reduced A2BP1 levels, rather than due to inter-individual variability. The RNA-Seq data thus provides validation of the functional groups of genes identified by co-expression analysis, and evidence for a convergence of transcriptional and alternative-splicing abnormalities in the synaptic and signaling pathogenesis of ASD.

To test whether our findings are more generalizable, and determine whether the autism-associated transcriptional differences observed are likely to be causal, versus collateral effects or environmentally-induced changes, we tested whether our co-expression modules or the DE genes show enrichment for autism genetic association signals. M12 showed highly significant enrichment for association signals (p = 5E-4), but neither M16 nor the list of DE genes showed such enrichment (Figure 4). As a negative control, we performed the same set-enrichment analysis using two GWAS studies for non-psychiatric disease performed on the same genotyping platform: a genome-wide association for hair color[23], and a GWAS study of Warfarin maintenance dose[24] finding no significant enrichment of the association signal (Figure 4b, Supplementary Figure 4). These results indicate that (a) M12 consists of a set of genes that are supported by independent lines of evidence to be causally involved in ASD pathophysiology, and (b) the upregulation of immune response genes in the autistic brain observed by us and others[25] has no evidence of a common genetic component.

Figure 4

GWAS set enrichment analysis

A. GWAS set enrichment analysis using the discovery AGRE cohort from Wang et al.[27] For each gene set (DE genes, M12 and M16) the null distribution of the enrichment score generated by 10000 random permutations is shown (x-axis) and the enrichment score for the gene set is depicted by a red vertical line. A p value < 0.01 was considered significant to correct for multiple comparisons.

B. GWAS signal enrichment of DE genes and the autism-associated co-expression modules M12 and M16. Enrichment p values are shown for an autism GWAS dataset (Wang et al.[27], AGRE discovery cohort) as well as two control datasets consisting of GWAS studies of non-psychiatric traits: Han et al.[23] (Negative control 1) and Cooper et al.[24] (Negative control 2). The red line marks the p value threshold for significance.

Our system-level analysis of the ASD brain transcriptome demonstrates the existence of convergent molecular abnormalities in ASD for the first time, providing a molecular neuropathological basis for the disease, whose genetic, epigenetic, or environmental etiologies can now be directly explored. The genome-wide analysis performed here significantly extends previous findings implicating synaptic dysfunction, as well as microglial and immune dysregulation in ASD[6] by providing an unbiased systematic assessment of transcriptional alterations and their genetic basis. We show that the transcriptome changes observed in ASD brain converge with GWAS data in supporting the genetic basis of synaptic and neuronal signaling dysfunction in ASD, while immune changes have a less pronounced genetic component and thus are most likely either secondary phenomena or caused by environmental factors. Since immune molecules and cells such as microglia play a role in synaptic development and function[26], we speculate that the observed immune up-regulation may be related to abnormal ongoing plasticity in ASD brain. The striking attenuation of gene expression differences observed here between frontal and temporal cortex in ASD likely represents a defect of developmental patterning and provides a strong rationale for further studies to assess the pervasiveness of transcriptional patterning abnormalities across the ASD brain. We also demonstrate for the first time alterations in differential splicing associated with A2BP1/FOX1 levels in ASD brain, and show that many of the affected exons belong to genes involved in synaptic function. Finally, given current evidence of genetic overlap between ASD and other neurodevelopmental disorders including schizophrenia and ADHD, our data provide a new pathway-based framework from which to assess the enrichment of genetic association signals in other psychiatric disorders.

METHODS SUMMARY

Brain Tissue

Postmortem brain tissue was obtained from the Autism Tissue Project and the Harvard Brain Bank as well as the MRC London Brain bank for Neurodegenerative Disease. Detailed information on the autism cases included in this study is available in Methods.

Microarrays and RNAseq

Total RNA was extracted from 100 mg of tissue using a Qiagen miRNA kit according to the manufacturer’s protocol. Expression profiles were obtained using Illumina Ref 8 v3 microarrays. RNAseq was performed on the Illumina GAIIx, as per the manufacturer’s instructions. Furthere detailed information on data analysis is available in Methods. All microarray and RNAseq data is deposited in GEO under accession numbers which are pending.

Full detailed Methods accompany this paper as supplemental material.