Genes implicated in ASD are highly co-expressed during human cortical development

We next asked whether genes associated with risk for ASD converge on common biological processes. We compiled a set of 155 ASD genetic risk candidates from the Simons Foundation Autism Research Initiative (SFARI) AutDB database (Basu et al., 2009), which we refer to as SFARI ASD. The SFARI ASD list is a manually curated set of candidate genes implicated by common variant association, candidate gene studies, genes within ASD-associated CNV, and, to a lesser extent, syndromic forms of ASD (Experimental Procedures). We mapped this gene set to the protein coding genes in the developmental co-expression network and observed that SFARI ASD genes are most over-represented in M16 (p = 0.0024, odds-ratio (OR) = 2.9, 95% confidence interval = [1.4-5.5], false discovery rate (FDR) < 0.05), and less so in M13 and M17 (Figure 2A),

Open in a separate window

Figure 2

Enrichment of SFARI ASD, asdM12 and ID genes in developmental networks

(A) Module-level enrichment for gene sets from a curated set of ASD risk genes (SFARI ASD), a curated set of ID genes (“ID all”), and an unbiased set of ASD risk genes (asdM12). Overlapping (ASD/ID overlap) and non-overlapping sets (“ASD only”, “ID only”) are also shown. All enrichment values for over-represented lists with p < 0.05, OR > 1 are shown to demonstrate enrichment trends (*p < 0.05, **FDR < 0.05).

(B), (C) and (D) show network plots for M13, M16, and M17 respectively. Most hub genes overlapping with SFARI ASD and asdM12 enrichment are not the same, showing that enrichment of these two sets is not driven by a narrow shared subset of genes.

Network plots comprise the top 200 connected genes (based on kME, a measure of intramodular connectivity) and their top 1000 connections in the sub-network. Genes with membership in SFARI ASD, asdM12, or the “ID all” list are labeled and plotted according to multidimensional scaling (MDS) of gene expression correlations, which graphs genes with similar expression patterns closer to each other. See also Table S2.

We also examined a set of ASD genes previously shown to be dysregulated in postmortem ASD temporal and frontal cortex (asdM12; Voineagu et al., 2011), which represents a shared molecular pathology in ASD brain identified in an unbiased, genome-wide manner. TheasdM12gene set was strongly enriched in the same three modules as SFARI ASD genes, M13, M16 and M17 (asdM12-M13, p = 3.0×10^-15, OR 3.6 [2.7-4.8]; asdM12-M16, p = 3.5×10^-15, OR 3.9 [2.8-5.3]; asdM12-M17, p = 1.0×10^-7, OR 2.5 [1.8-3.4]; each at FDR< 0.05). A remarkable 42% of asdM12 and 25% of the SFARI ASD sets are found in one of these three modules. Our analysis, which uses gene sets identified based on different methods (only 15 genes overlap between SFARI ASD and asdM12), converges onto three modules involved in prenatal and early postnatal synaptic development.

We next hypothesized that mapping ID genes to this network would enable us to assess whether ASD susceptibility genes show any specificity in their developmental expression patterns. We compiled an extensive set of high confidence genes implicated in monogenic forms of ID from multiple publications (In low and Restifo, 2004; Lubs et al., 2012; Ropers, 2008; van Bokhoven, 2011), referred to as “ID all” (see Experimental Procedures). Remarkably, this set of364 genes expressed in human neocortexis not enriched in any of the 12co-expression modules. Importantly, this lack of enrichment is at a relaxed threshold that reduces the risk of false negatives (uncorrected p > 0.05). Removing the small set of 37 genes (<10%) that overlap between ASD and ID to establish exclusive sets (“ASD only”, “ID only”) further confirms that ASD genes exhibit enrichment, while ID genes do not (Figure 2A, Table S2B). Thus, it is genes connected with the ASD phenotype that are enriched in 3 specific transcriptional modules related to synaptic function during development, but not those that have been related solely to ID.

RDNVs are highly enriched in two co-expression modules in early fetal development

Additional evidence implicating specific genes in ASD comes from whole-exome sequencing in families (Iossifov et al., 2012; Neale et al., 2012; O'Roak et al., 2012b; Sanders et al., 2012), which has identified many rare protein disrupting variants (nonsense, splice-site, frameshift) over-represented in individuals with ASD compared to their unaffected siblings (OR >2). This evidence is largely distinct from the evidence implicating genes in SFARI ASD and asdM12, as it is from purely non-inherited, rare variation discovered in an unbiased, genome-wide manner. We therefore asked whether RDNV-affected genes found in ASD probands shared biological function. We also tested silent RDNVs since they should not exhibit a similar pattern of functional enrichment, providing a key control for gene size, GC content, and other features affecting mutability (Michaelson et al., 2012).

We first tested for enrichment using RDNVs from three studies sharing similar coverage criteria and variant calling methodology (Neale et al., 2012; O'Roak et al., 2012b; Sanders et al., 2012), representing 622 ASD probands and 222 unaffected siblings. Strikingly, genes expressed during development and affected by protein disrupting RDNVs in probands (60 genes, Table S2A, Discovery Set) are significantly enriched in two modules, M2 and M3, which exhibithighly similar developmental trajectories and functional enrichment, indicative of remarkable biological specificity. Eight genes harboring protein disrupting RDNVs are enriched in M2 (p = 0.006, OR = 3.2 [1.3-6.8]; FDR < 0.05) and 10 are enriched in M3 (p = 0.0011, OR = 3.6 [1.6-7.2]; FDR < 0.05). A trend for enrichment is observed for M16 as well, but this does not pass the FDR threshold. For comparison, genes affected by RDNVs in unaffected siblings or affected by silent mutations are not enriched in any modules (Table S2B, p > 0.05). Since missense RDNVs are only weakly over-represented in ASD (Sanders et al., 2012),we reasoned that overlap with network modules might prioritize specific subsets of this RDNV class. We find that a subset of missense RDNV affected genes is over-represented in the same pathways as the more deleterious protein disrupting RDNVs (M2 and M3, Table S2B). Taken together, out of 385protein disrupting or missense RDNV affected genes expressed in brain, 34 are found in M2 (p = 2.9 × 10^-4, OR = 2.1 [1.4-3.0], FDR < 0.05) and 41 in M3 (p = 2.3 × 10^-5, OR = 2.2 [1.5-3.1], FDR < 0.05). There is no enrichment for this combined set in any other modules. Furthermore, the combined set of protein disrupting and missense RDNVs from unaffected siblings was not found enriched any modules (p > 0.05).

We further validated the observed RDNV enrichment pattern in M2 and M3 in an independent set of patients from a study with more stringent RDNV calling criteria (Iossifov et al., 2012). In this additional set of 343 ASD probands and unaffected siblings, we found that the patterns of RDNV enrichment replicated, with the set of protein disrupting and missense RDNVs from ASD probands enriched specifically in M2 and M3 (p < 0.05), and RDNVs from siblings and silent RDNVS not enriched in any set (Table S2B, Replication Set). Combining across all studies, we find that out of 598protein disrupting or missense RDNV affected genes expressed in brain, 52 are in M2 (p = 9.6 × 10^-6, OR = 2.0 [1.5-2.8]) and 61 are in M3 (p = 8.5 × 10^-7, OR = 2.1 [1.6-2.8]). Importantly, the enrichment pattern across modules is not only replicated in the independent set, but is stronger in the combined set, is robust to perturbations in module composition (Figure S3A), and is not driven by variants from any one study (Table S2C-D). We show the enrichment pattern of this combined set across 965 ASD probands and 565 unaffected siblings in Figure 3A and use this combined set for the remainder of our analyses.

Open in a separate window

Figure 3

Enrichment of genes affected by RDNVs in developmental networks

(A) Module-level enrichment for multiple categories of RDNV in ASD affected probands and unaffected siblings combined across four studies. All enrichment values for over represented lists with p < 0.05, OR > 1 are shown to demonstrate enrichment trends (*p < 0.05, **validated in replication set) marked. M2 and M3 are strongly enriched for protein disrupting and missense RDNV affected genes in probands. Enrichment for genes affected by silent RDNVs in probands and RDNV gene sets affected in siblings represent control gene sets and do not show enrichment.

(B) and C) show network plots for M2 and M3, with all genes plotted and all genes carrying RDNVs displayed.

Network plots show all genes in the module with protein disrupting or missense RDNV-affected genes highlighted. Those with high intramodular connectivity (kME > 0.75) are labeled in black, with the rest in grey. The top 1000 connections are shown, and genes are plotted according to the multidimensional scaling (MDS) of co-expression as in Figure 2. See also Figure S2-3 and Table S2.

We next asked whether M2 and M3 prioritized functional subsets of genes with RDNVs. We confirmed that RDNV-affected genes in M2 and M3 are significantly enriched for interactions at a protein level (Figure S2A-D), and highlight genes that are both PPI hubs and co-expression hubs in Figure 3B-C. Furthermore, M2 and M3 genes harboring RDNVs are also more dosage sensitive, as evidenced by the significant increase in the probability of haploinsufficiency (P(HI), Extended Experimental Methods) among genes affected by these mutation classes(Huang et al., 2010; Luo et al., 2012). This is consistent with the heterozygous state of variants observed in ASD probands. Overall, a remarkable proportion, 113/598 (19%) of genes affected by known RDNVs are co-expressed in two modules reflecting similar temporal trends of highexpression in cortex during the neurodevelopmental period of early neuronal fate determination, migration, and cortical lamination. Of note, as with M13, M16, and M17, which were enriched for asdM12 and SFARI ASD, ID genes showed no enrichment in M2 or M3 (p > 0.05).

We also observed that the SFARI ASD genes and asdM12 genes, which are enriched for inherited common variants in ASD (small average effect size), affect the synaptic modules, M13,M16, and M17. In contrast, the non-inherited (larger average effect size) RDNVspreferentially affect the early transcriptional regulation modules (Extended Experimental Methods). We emphasize that this is not absolute, as M16 includes some genes harboring RDNVs (e.g. in SCN2A, SHANK2, NRXN1; Figure 2A). To formally assess common variant enrichment using independent data, we compared ASD GWA signals across these modules (Extended Experimental Methods). Genes inM13 and M16 were more strongly affected by common variation in at least one of two ASD GWA studies(Anney et al., 2012; Wang et al., 2009) than M2 or M3 (Figure S3E). This is consistent with susceptibility of distinct biological processes for different mutational classes, and that in general more severe biological consequences would result from early transcriptional dysregulation during neuronal proliferation and differentiation, compared with later disruption of synaptic development and neuronal function.

ASD gene enriched modules are linked by translational and transcriptional regulation

Upregulated and downregulated modules are highly anti-correlated throughout development, so we hypothesized that common molecular regulatory relationships could potentially link genes within these modules. We first used a set of FMRP-RNA interactors from a cross-linking and immunoprecipitation (CLIP) experiment (Darnell et al., 2011), since Iossifov et al. (2012) had previously shown that RDNVs identified in their exome sequencing study were enriched in this class of genes. Remarkably FMRP targets are specifically enriched in modules that also contain ASD-related genes M2, M16, and M17 (FMRP-M2 p = 1.6×10^-13, OR = 3.0 [2.3-3.9]; FMRP-M16 p = 2.4×10^-29, OR =5.7 [4.3-7.6]; FMRP-M17 p = 9.3×10^-10, OR =2.4 [1.8-3.1]; all at FDR < 0.05; Figure 4A). This provides a strong, independent line of evidence that translational regulation by FMRP not only affects genes harboring RDNVs, but links different molecular pathways that are co-expressed during early fetal cortical development and susceptible to diverse classes of ASD genetic mutation.

Open in a separate window

Figure 4

Translational and transcriptional co-regulation connects developmentally distinct ASD-affected modules

(A) Co-expression based network plot of FMRP interactions with genes in M2, M16, and M17 that are either affected by RDNVs or in an ASD candidate list. Genes are plotted as in Figure 2 and ​and3,3, but now across modules, with FMRP placed at the center.

(B) Summary of TF binding site (TFBS) enrichment in modules for TFs that have evidence for function in a neurodevelopmental context. Dashed lines indicate enrichment in the module for predicted binding sites.

(C-G) MEF2A, MEF2C, SATB1, FOXO1, and ELF1 are all enriched for their binding motifs in the upstream regions of ASD gene-enriched modules following anti-correlated developmental patterns. Network plots highlight genes with a predicted binding site (light dashed arrow) contributing to this enrichment that are also affected by RDNVs or in an ASD candidate list. Arrows representing a TFBS found in a ChIP experiment are marked in dark blue.

For network plots, the top 1000 positive connections between genes are plotted and node size is proportional to connectivity within the genes' assigned module, therefore larger nodes are more central hubs. See also Table S3.

We next tested whether ASD associated modules are also linked at the transcriptional level (Experimental Procedures). We found17 TFs that are predicted to link at least one upregulated and one downregulated module based on binding site enrichment (Figure 4B, Table S3A-B). Many of these TFs are expressed during fetal development (Table S1A), have been previously implicated in relevant neuronal functions, and have DNA binding targets have been experimentally characterized (Table S3B). For example, MEF2A and MEF2C, both members of a TF family regulating synaptic plasticity and glutamatergic synapse number (Ebert and Greenberg, 2013), are enriched for binding targets in M2 and M17, which are anti-correlated across development (Figure 4C-D). SATB1, which is required for the development of cortical interneurons(Close et al., 2012), ELF1, which is involved in axonal guidance, and FOXO1 which regulates neuronal polarity(de la Torre-Ubieta and Bonni, 2011) also link these two modules (Figure 4E-F). To provide further evidence that these are experimentally plausible binding sites, we overlaid these bioinformatic predictions with chromatin immunoprecipitation (ChIP) data where available, supporting many of these predicted interactions, including 39% of MEF2A, 23% of MEF2C and 87% of ELF1binding sites (Figure 4C, 4D, 4G; Extended Experimental Procedures). These results implicate existing and novel TFs as putative co-regulators of ASD-associated gene networks during neocortical development.

Weighted Gene Co-expression Network Analysis

We used the R package WGCNA (Langfelder et al., 2008) to construct co-expression networks, as previously done (Voineagu et al., 2011) and described in detail in Extended Experimental Methods. The modules were characterized using GO Elite to control the network-wide false discovery rate, with all enriched pathways comprising at least 10 genes at Z > 2 and FDR < 0.01 (Zambon et al., 2012). All network plots were constructed using the igraph package in R (Csárdi and Nepusz, 2006).

Protein-protein interaction enrichment

When assessing PPI enrichment in modules, a degree-matched permutation analysis was applied in order to control for biological and methodological biases in PPI data (see Extended Experimental Procedures for details).

Gene sets

The SFARI ASD set was compiled using the online SFARI gene database, AutDB. We used the “Gene Score,” which classifies evidence levels, to restrict our set to those categorized as S (Syndromic) and evidence levels 1-4 (high confidence - minimal evidence). We obtained asdM12 and adsM16 from a prior, independent gene expression study that profiled expression changes in ASD cortex and applied WGCNA to identify modules of dysregulated genes ASD (Voineagu et al., 2011). We curated ID genes from four reviews cataloging genes causing “ID all” ( In low and Restifo, 2004; Lubs et al., 2012; Ropers, 2008; van Bokhoven, 2011) resulting in 401 genes. For candidate lists, we used the HUGO gene nomenclature to find updated gene symbols. We obtained RDNVs from four publications (Iossifov et al., 2012; Neale et al., 2012; O'Roak et al., 2012b; Sanders et al., 2012), and split them into discovery and validation sets as discussed in the results (see Extended Experimental Procedures for further details about gene sets).

Gene set over-representation

All enrichments of gene sets were performed using a two-sided Fisher exact test with 95% confidence calculated according to the R function fisher.test. The false discovery rate was controlled across candidate ASD gene set enrichments, the discovery RDNV set enrichment, and FMRP target enrichment (Table S2B). For RDNV enrichment, we required an OR > 1 and FDR < 0.05 for enrichment in the discovery set, and OR > 1 with p < 0.05 for validation in the replication set. When showing a lack of enrichment, we require p > 0.05 to reduce false negatives, as future studies that add expression time points for networks and genes for enrichment may find enrichment in pathways not enriched here.

Transcription Factor Binding Site Enrichment

The top 200 genes in each module (ranked by kME) were used for TF motif enrichment analysis. Enrichment for each TF motif in TRANSFAC (Matys et al., 2003) was compared to 3 background datasets to ensure robustness: 1000 bp sequences upstream of all human genes, human CpG islands, and the sequence of human chromosome 20 (Extended Experimental Procedures). Only TFs with p < 0.05 across all backgrounds are considered enriched. ChIP data was obtained from ENCODE (Consortium, 2011) and the ChIP enrichment analysis (Lachmann et al., 2010) resource.

We thank the Simons Foundation Autism Research Initiative (www.sfari.org) and the Autism Genetic Resource Exchange (www.agre.org) and all families involved for making this work possible. We gratefully acknowledge data resources from the BrainSpan consortium (NIMH grant 5RC2 MH089921) and the Allen Brain Institute (www.brainspan.org, www.brain-map.org). This work supported by a NIMH Training Grant and NRSA Fellowship (T32MH073526 and F30MH099886, NNP), an R01 (5R01MH060233:05, DHG), an Autism Center for Excellence network grant (5R01MH100027), and the Medical Scientist Training Program at UCLA. We also acknowledge the support of the NINDS Informatics Center for Neurogenetics and Neurogenomics (P30NS062691) at UCLA. We thank Jason Chen, Michael Gandal, and Jason Stein for critically reading the manuscript, as well as Willsey et al. for sharing their results prior to publication.