Convergence and Divergence of Common and Rare Variants of Autism Spectrum Disorders in Tissue-specific Pathways and Gene Networks

TL;DR

This study explores how both common and rare genetic variants contribute to autism by analyzing gene networks in brain and peripheral tissues.

Contribution

The study integrates multi-tissue, multiomics data to reveal shared and distinct roles of common and rare variants in autism.

Findings

Brain regions involved in synaptic signaling and neurodevelopment are enriched for both common and rare ASD variants.

Peripheral tissues like the digestive and immune systems are mainly influenced by common variants.

Key regulators like SYT1 and ADD2 may coordinate the effects of both variant types in ASD.

Abstract

The genetic heterogeneity of autism spectrum disorder (ASD) presents significant challenges in understanding its pathogenic mechanisms, as the genetic risk involves numerous common variants and rare de novo or inherited variants. Prior research has mainly focused on identifying rare variants and their impact on neurodevelopment and neuronal functions in cortical brain regions. By contrast, common variants, which contribute substantially to ASD heritability, remain understudied, suggesting a need to consider both variant types to understand ASD’s genetic mechanisms. Previous studies have also implicated subcortical brain regions and peripheral digestive and immune systems, but tissue-specific mechanisms remain unclear. We address these knowledge gaps by identifying gene networks, pathways, and key regulators informed by ASD common variants in brain and peripheral tissues, further examining whether these networks also capture genes informed by rare variants. Our approach integrates genome wide association study (GWAS) summary statistics, tissue-level genetics of gene expression, and gene coexpression and transcriptional regulatory networks across ~50 tissues. Our multitissue, multiomics analysis reveals that key brain regions and networks crucial for synaptic signaling and neurodevelopment are enriched for both rare and common variants, whereas peripheral tissues, such as the digestive and immune systems, are primarily informed by common variants. This partitioning of key tissues and biological pathways into core (targeted by both variant types) and modifying components provide insight into ASD heterogeneity. We also identified central gene network regulators, such as SYT1 and ADD2, which may orchestrate the effects of both common and rare ASD genetic risk factors on ASD pathogenesis.

Figures

C](#F2); [Figure 3](#F3)). The abundance of both brain and peripheral coexpression modules enriched for ASD GWAS signals ([Figure 2B](#F2)) supports the growing body of research that suggests ASD relevance within not only the brain but the digestive, endocrine, and immune systems [[15](#R15), [16](#R16), [68](#R68), [71](#R71), [74](#R74), [85](#R85)], and further provides potential causal inference for their involvement due to the incorporation of genetic variants in our analysis. The notable convergence of known rare and common variants in the brain and the identified key drivers regulating

C](#F4) shows the percentage of key drivers across both brain and peripheral tissues that are also either known genes harboring rare variants from the SFARI database, or contain common variants from the GWAS. We observed a higher percentage of known rare ASD variants among brain tissue-derived KDs compared to peripheral tissue-derived ones. As the rare variants mostly affect brain development and neuronal functions and have larger effect sizes, our results align with the central role of brain tissues in the pathogenesis of ASD particularly for ASD driven by rare variants. By contrast, there wa

A](#F5)–[B](#F5)). As expected, a greater number of brain tissue key driver subnetworks contained a significant number of both rare and common variants, which suggest that these key drivers are crucial regulators within ASD pathogenesis given the convergence of these two types of variants ([Figure 5C](#F5)). The rankings of key driver subnetworks in terms of their rare or common variant enrichment highlight the central role of the brain in ASD, yet also point towards the digestive system as the leading peripheral region for common variant abundance ([Figure 5A](#F5)–[B](#F5)). Our analysis als

B](#F6)), *SCN8A, AMPH, ATP9A* (Supplementary Figure 1A–C), demonstrated enrichment of both variant types, thereby emphasizing the convergence between these variants orchestrated by these brain KDs. For instance, the ADD2 subnetwork had both rare (e.g., *SHANK2, SCN8A*) and common variants (e.g., *ACTN2, DMN1*) in its neighborhood ([Figure 6B](#F6)). Adducin genes encode cytoskeleton proteins that are critical for osmotic rigidity and cell shape by regulating the formation of the spectrin-actin membrane skeleton [[50](#R50)]. *ADD2* is expressed in the brain, and its knock-out results in the l

Full text

INTRODUCTION

Autism spectrum disorder (ASD) is a complex neurodevelopmental condition that manifests through social and communication deficits and various behavioral abnormalities [1, 2]. The prevalence of ASD has grown in recent years, likely due to improved diagnosis, with the latest epidemiological estimates suggesting more than 50 million people worldwide with some form of autism [3]. The current autism diagnosis criteria from the American Psychiatric Association’s Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition include social-communication challenges (e.g., difficulties in social-emotional reciprocity, nonverbal communication, and relationships) and behavioral characteristics (e.g., repetitive movements, strict routines, and sensory sensitivities), and it categorizes ASD into three severity levels [4].

The wide range of social and communication deficits, behavioral abnormalities, and varying levels of severity indicate a diversity in ASD pathogenesis. This is likely rooted in substantial genetic heterogeneity as numerous rare mutations and common genetic variants with varying effect sizes have been identified in ASD [5]. Despite possessing a greater risk of causing significant autism syndromes, rare variants comprise of only ~1% of patients with autism [6]. Common variants, in contrast, are more prevalent and collectively contribute to >50% of ASD heritability despite smaller effect sizes associated with individual variants [5, 7], and can be summarized into an individual’s polygenic risk score to quantify the genetic predisposition to ASD [8].

Previous studies have identified key brain regions and their associated molecular and cellular pathways involved in ASD pathogenesis. The frontal and temporal cortical regions have shown abnormal gene expression patterns in autistic patients compared to typically developing individuals [9]. Additionally, two regions vital for complex cognitive processes, the anterior cingulate cortex and the amygdala, have exhibited decreased neuronal activity and abnormal growth, respectively, in autistic patients [10, 11]. There has also been evidence for non-symmetric development in the lateral ventricles and hippocampus [12]. Other implicated brain regions include the prefrontal, parietal, and visual cortices, the cerebellum, and the caudate nucleus [13, 14]. Interestingly, peripheral systems have also been implicated in ASD. For example, immune dysregulation has been associated with altered neurodevelopment and behavior [15]. The microbiota-gut-brain axis has also been explored given the interactions between commensal bacteria, immune cells, enteric nerves, and neurotransmitters, and the observation that ASD patients frequently present with gastrointestinal complications [16].

Despite these findings, the specific the tissues, genes, and pathways underlying ASD are incompletely understood, and no effective treatments currently exist [17]. To address this, examining available omics data in a tissue-specific context is crucial as it may reveal relevant and potentially causal disease-associated mechanisms informed by common or rare variants. With this in mind, we hypothesize that ASD follows an omnigenic model of pathogenesis, where hub genes with large effect sizes interact with peripheral genes with smaller effect sizes via interconnected networks [18]. It is plausible that rare variants are enriched among the hub genes across essential tissues whereas common variants are enriched among peripheral genes within a broader range of tissues. Therefore, elucidating how common and rare variants of ASD converge and diverge in tissue-specific gene networks will identify key tissues and gene drivers within gene regulatory networks, providing further insights for future mechanistic and therapeutic studies.

In the current study, we integrate common variants from a genome-wide association studie (GWAS) of ASD with tissue-specific gene regulation in the form of expression and splicing quantitative trait loci (eQTLs, sQTLs) and gene coexpression and regulatory networks to identify networks, pathways, and key regulators affected by common variants. We further compare these findings with ASD rare variants to uncover commonalities and differences between the two types of variants in terms of tissue-, pathway-, and network specificity to distinguish between core and modifying mechanisms.

METHODS

Analysis overview

We utilize a multiomics integration tool, Mergeomics, for our analysis of ASD [19, 20], as its focus on integrating genetic variants with functional genomics to retrieve potential causal mechanisms is well suited for the goal of the study. The robustness of this tool has been substantiated by experimental validation of its computational predictions, and it has been successfully applied to the analysis of multiple complex diseases [21–24]. Briefly, we integrated full summary statistics of an ASD genome-wide association study (GWAS) with tissue-specific expression and splicing quantitative trait loci (eQTLs/sQTLs) and tissue-specific gene coexpression networks to allow for the ranked identification of pathways and gene subnetworks most associated with ASD based on common variants examined in the GWAS (Figure 1). The pipeline then performs a key driver analysis to determine network hub genes, termed “key drivers”, whose neighboring networks are enriched for disease-associated genes within interconnected gene regulatory networks. The outputs of Mergeomics include a ranking of biological pathways and subnetworks informed by ASD GWAS common variants in a tissue-specific manner, and a visualization of key drivers within tissue-specific disease subnetworks. Advantages of Mergeomics are: i) it utilizes the full GWAS summary statistics representing the spectrum of disease association strengths of genetic variants, and ii) it contains a unique test statistic that summarizes disease association enrichment at multiple quantile thresholds to derive stable statistics that are less dependent on any given GWAS significance cutoff to account for variability in sample size and statistical power. Figure 1 depicts our overall pipeline, datasets utilized, and the three steps of the Mergeomics analysis that are discussed in further detail below: marker dependency filtering (MDF), marker set enrichment analysis (MSEA), and key driver analysis (KDA).

ASD common variant GWAS summary statistics

We used ASD GWAS summary statistics with association p-values for over 9 million SNPs, derived from a study that included >18,000 ASD cases and ~28,000 controls from a Central European (CEU) population [25]. Quality control and GWAS analyses were performed using Ricopili, PLINK, and METAL pipelines, adjusting for batch effects and principal components to account for population stratification [26–28]. We chose to use this dataset because it is the largest ASD GWAS to date. However, we note that the population is primarily CEU and may have limited generalizability to other populations.

Marker dependency filtering (MDF)

SNP-to-gene mapping was performed using ASD GWAS SNPs and tissue-specific expression and splicing quantitative trait loci (eQTLs and sQTLs) from the Genotype Tissue Expression (GTEx) project database [29], and distance-based mapping (±20 kb) was used as an alternative mapping method to capture cis-regulatory relations. The complete list of tissue-specific eQTL/sQTL data for 49 tissues used as input for the Mergeomics analysis are in Supplementary Table 1. We corrected for linkage disequilibrium (LD) to filter known dependencies between SNPs based on an LD cutoff threshold of r^2^ > 0.5 (the default threshold for Mergeomics based on performance benchmarking) from the CEU population, as the ASD GWAS population is mainly CEU. The output of MDF contained tissue-specific SNP-to-gene mappings based on eQTL, sQTL, and distance-based methods, with their ASD association strengths represented by −log10 p-values after filtering out SNPs in LD by keeping the SNP with the strongest GWAS association within each LD block. The LD clumping step was done to remove redundant SNPs entering downstream analysis, and the mapping of SNPs to genes via eQTLs and sQTLs was performed to offer tissue and gene annotations for GWAS SNPs based on data-driven potential tissue-specific gene regulation. No GWAS p-value cutoffs were applied, ensuring a comprehensive range of association signals for downstream null distribution estimation in MSEA.

Tissue-specific Weighted Gene Coexpression Network Analysis (WGCNA) to define data-driven functional gene sets in each tissue

To group genes with functional relevance in individual tissues in a data-driven manner, we used the transcriptome data from the GTEx database to construct tissue-specific WGCNA gene coexpression modules [30]. Typically, these modules contain genes that are coexpressed and functionally related. This provided a means of placing our ASD-associated GWAS genes mapped from MDF into categories that have biological relevance in individual tissue contexts. These modules were functionally annotated through pathway enrichment analysis using KEGG [31], Reactome [32], and BioCarta [33] databases. We constructed WGCNA networks for all tissues in GTEx. To facilitate the interpretation of tissues that generated results after marker set enrichment analysis, we categorized the tissues into seven systems: Adipose/Immune (subcutaneous adipose, visceral adipose, EBV transformed lymphocytes, spleen, whole blood), Brain (amygdala, anterior cingulate cortex, caudate, cerebellar hemisphere, cortex, frontal cortex, hippocampus, hypothalamus, nucleus accumbens, putamen, pituitary, spinal cord, substantia nigra), Cardiovascular (aortic artery, atrial appendage, coronary artery, left ventricle, tibial artery), Digestive (esophagus mucosa, esophagus muscularis, gastroesophageal junction, liver, minor salivary gland, pancreas, sigmoid colon, small intestine terminal ileum, stomach, transverse colon), Endocrine (adrenal, ovary, pituitary, testis, thyroid), Female Reproductive (breast mammary tissue, endocervix, ectocervix, fallopian tube, ovary, uterus, vagina), and Male Reproductive (prostate, testis).

Marker set enrichment analysis (MSEA)

We next used MSEA to assess whether tissue-specific WGCNA coexpression modules from GTEx tissues are enriched for SNPs showing stronger ASD GWAS associations compared to a null distribution of random gene sets. GWAS-mapped genes from MDF and tissue-specific WGCNA coexpression modules were used as input, with a chi-like statistic used as the enrichment test statistic to determine coexpression modules enriched in the ASD GWAS:

$$[eqn]$$

Briefly, n represents the number of quantile points used to divide the GWAS SNPs into significant and non-significant groups. Rank-based quantiles were used instead of specific p-value cutoffs to mitigate the sample size influence on p-values. We used 10 quantiles ranging from 0.5 to an upper limit (qupper) adjusted by the median module size (μ). This upper limit was calculated as qupper = 1.0 – 1.0/μ, allowing the quantile distribution to account for module size.

Oi and Ei represent observed and expected positive associations at each quantile point, with the difference divided by the square root of the expected count plus a stability parameter, κ = 1, to account for extremely low expected counts. Thus, the calculated χ value is a sum of the output at each quantile point from a given test gene set or module. A null distribution was created by generating random gene sets matching the gene number of each coexpression module and calculating its χ value, with the following null hypothesis: Given all distinct markers from a set of N genes, these markers contain an equal proportion of positive associations when compared to all the distinct markers from a set of N random genes [19]. Z-scores were then calculated to measure the deviation of each module’s χ value from the null distribution’s mean, and enrichment scores were determined from these Z-scores to rank tissue-specific modules by their enrichment for ASD GWAS signals. P-values were calculated from the Z scores using a Gaussian distribution, and a false discovery rate (FDR) was estimated using the Benjamini-Hochberg method.

Tissue-specific Bayesian gene regulatory networks to define gene regulatory relations

WGCNA coexpression modules focus on correlation but not regulatory information between genes. Therefore, to elucidate directional gene regulatory relations, we also constructed Bayesian networks from tissue-specific GTEx databases using the RIMBANet package [34] for all tissues, which incorporates priors from eQTLs, transcription factors to target relations, and mutual information to model gene-gene regulation. As Bayesian networks from individual datasets are typically sparse, networks from similar tissues were subsequently merged to derive composite networks for brain, digestive, cardiovascular, endocrine, immune, adipose, and reproductive tissues to reduce sparsity and ensure each network contains at least 10,000 genes. We further merged the immune and adipose tissue networks given the role of adipose tissue in immune regulation [35]. The merged networks and their corresponding constituent tissues and network sizes are shown in Supplementary Table 2.

Key driver analysis (KDA)

We used KDA to identify key drivers (i.e., hub network genes) of ASD-associated gene coexpression modules from MSEA and their associated neighbors within tissue-specific Bayesian gene regulatory networks constructed from GTEx tissues. Key drivers were predicted by identifying hub nodes in the top 25% of edge connections. For each hub node and its subnetwork, KDA assessed enrichment in ASD-associated gene sets using a calculation similar to MSEA. The proportion of subnetwork nodes overlapping with ASD-associated modules was compared to a null distribution generated from reshuffled subnetworks of the same size to calculate an enrichment statistic. The outputted networks were ranked by their enrichment for ASD GWAS-informed gene sets from MSEA. Using EnrichR, a gene set enrichment annotation tool [36–38], we performed pathway enrichment analyses to understand the functions of predicted key driver subnetworks.

Assessment of rare variants in tissue-specific key driver (KD) networks of ASD common variants

ASD rare variants were obtained from the Simons Foundation Autism Research Initiative (SFARI), which stratifies genes into four confidence levels: Level Syndromic (high confidence in both ASD and a specific syndrome beyond the characteristics of ASD), Level 1 (high confidence in their implication in ASD), Level 2 (strong candidate for ASD association), and Level 3 (moderate evidence based on previous research) [39].

Within each tissue-specific KD subnetwork identified for ASD common variants, genes in the subnetworks that contain known rare variants were annotated based on information from the SFARI database. We considered the rare variants both as a collective across each of the ASD levels and at each individual stratification for the analysis. Enrichment for rare variant genes within each common variant-informed gene subnetworks was assessed using Fisher’s exact test.

KD subnetwork visualization

Cytoscape (version 3.10.1) [40] was used to visualize each of the top ranked KDs and their subnetworks. Genes mapped to common and rare ASD variants were annotated for each KD subnetwork.

RESULTS

Marker set enrichment analysis (MSEA) reveals tissue-specific coexpression modules enriched for ASD GWAS signals

The MSEA analysis identified 48 tissues (11 brain regions, 37 peripheral tissues) out of the 53 tissues analyzed from which at least one coexpression module was significantly enriched for ASD associations in GWAS at FDR < 5% (Figure 2A). Across these tissues, there were 196 ASD-enriched WGCNA modules. Supplementary Table 3 shows the top 10 modules based on the statistical significance of ASD GWAS enrichment, where coexpression modules with diverse functional annotations from 5 brain regions and 5 peripheral tissues were observed. The anterior cingulate cortex, which is involved in emotional regulation and cognitive control, and the amygdala, another region crucial for emotional response, contained two coexpression modules with the highest ASD enrichment. The top functional annotations for these two modules included protein activity, immune system, and neuroplasticity. The other top brain tissue coexpression modules were from the frontal cortex (Brodmann Area 9, [BA9]), cerebellum, and cortex, and were enriched for pathways involved in neuronal signal transduction and immune regulation. Together these findings support how impaired neuronal and immune functions play a significant role in ASD pathogenesis. Interestingly, peripheral tissues from the digestive system (liver, minor salivary gland), reproductive system (uterus, testis), and immune system (lymphocytes) contained highly significant coexpression modules relevant to mRNA splicing, immune pathways, cell cycle, and mammalian target of rapamycin (mTOR) signaling, which have been implicated in ASD [15, 41, 42]. When further examining the genes in these pathways with the strongest GWAS association p-values, we note an abundance of non-coding RNAs (ncRNAs) in addition to genes important for neuronal and immune functions. These results align with recent studies that implicate the roles of ncRNAs in ASD and other related disorders [43–47].

We next examined all significant modules across brain and peripheral tissues, expecting significant modules from brain tissues to collectively show higher levels of significance for ASD GWAS enrichment compared to modules from peripheral tissues given that ASD is a neurodevelopmental disorder. Surprisingly, we did not find significant difference in the average FDRs between brain and peripheral tissue modules (Figure 2B), suggesting that peripheral tissues are not less informative than brain tissues in capturing ASD common variants in coherent functional modules. Across all significant modules, pathway annotation revealed a broad range of consistent pathways including cell cycle, gene regulation (particularly splicing), neuronal signaling, oxidative phosphorylation, immune system, and mTOR signaling (Figure 2C).

To further explore the pathways within the enriched coexpression modules of tissues with the most relevance to ASD, we categorized tissues into seven systems: Adipose/Immune, Brain, Cardiovascular, Digestive, Endocrine, Female Reproductive, and Male Reproductive. There were nearly 800 pathways across all tissue systems, among which 80 were shared across all seven systems (Figure 3), mainly encompassing immune pathways, cell cycle regulation, cellular signaling, growth and proliferation, and protein interactions. Among the system-specific pathways, we observed crucial neuronal and signaling processes (postsynaptic activation via glutamate binding to NMDA receptor, Ras activated CREB phosphorylation, olfactory signaling and transduction) and neurodegenerative conditions (Alzheimer’s disease, Amyotrophic Lateral Sclerosis) for the brain; immune regulation involving the CD40 pathway and circadian clock for the Adipose/Immune system; glycosphingolipid metabolism and neuroactive ligand-receptor interaction for the Endocrine system; and hyaluronan uptake/degradation and MAPK signaling for integrins within the digestive system. The tissue-specific modules and processes suggest diverse and unique contributions of various systems to ASD pathogenesis.

Key driver analysis (KDA) identifies distinctions between brain and peripheral tissue associations in ASD

Using the significant tissue-specific coexpression modules identified from MSEA, we performed KDA to identify tissue-specific key drivers. We also assessed whether the key drivers and their subnetworks capture both common and rare variants of ASD by intersecting the subnetworks with genes containing these variants. First, as seen in Figure 4A, we observed higher numbers of rare and common variants in the key drivers from the brain Bayesian networks compared to the peripheral networks, supporting the importance of the brain in ASD, as expected. Besides brain KDs, digestive system KDs contain the second highest numbers of rare and common variants. However, when normalizing the gene counts of rare and common variants against the total number of genes in the peripheral and brain Bayesian networks, the normalized count was higher for peripheral tissue networks (Figure 4B), which is likely due to the much smaller sizes of the peripheral tissue networks compared to more complex brain networks (Supplementary Table 2). This finding supports that peripheral tissue genes possess sizable contributions to the overall genetic burden of ASD.

We next focused on the identified KDs and sought to assess the differences in the degrees of statistical enrichment of common and rare variants among KDs between the brain and peripheral tissues. Figure 4C shows the percentage of key drivers across both brain and peripheral tissues that are also either known genes harboring rare variants from the SFARI database, or contain common variants from the GWAS. We observed a higher percentage of known rare ASD variants among brain tissue-derived KDs compared to peripheral tissue-derived ones. As the rare variants mostly affect brain development and neuronal functions and have larger effect sizes, our results align with the central role of brain tissues in the pathogenesis of ASD particularly for ASD driven by rare variants. By contrast, there was a similar percentage of common variants among KDs between brain and peripheral tissues, suggesting that common variants of ASD are less discriminative between the brain and the peripheral tissues.

We then statistically assessed whether the KDs were enriched for common and rare variants of ASD using the hypergeometric test. As shown in Figure 4D, brain tissue KDs showed much stronger enrichment for both the rare and common ASD variants than peripheral KDs, further supporting the importance of the brain in ASD. Moreover, the brain KDs showed a significantly higher enrichment for rare variants than for common variants. As rare variants have larger effect sizes than common variants, the observation of stronger enrichment of rare variants among brain KDs also supports the stronger influence of brain networks in ASD.

Prioritization of KDs based on rare/common variant enrichment in KD subnetworks

We prioritized KDs based on the significance of their subnetworks for common and rare variant enrichment. Among the top 15 KD subnetworks in terms of their rare (Figure 5A) and common (Figure 5B) variant enrichment, respectively, there was a higher abundance of brain tissue KD subnetworks (12 out of 15) that were enriched for rare variants, and an increasing representation of digestive tissue KD subnetworks (5 out of 15) that had enrichment for common variants. There was also a higher proportion of brain tissue KD subnetworks enriched for both rare and common variants (Figure 5C; Supplementary Table 4), supporting convergence between rare and common variants in these KD subnetworks.

Two brain tissue KD subnetworks stood out due to their high enrichment of either common or rare variants. Synaptotagmin 1 (SYT1, Figure 6A) held the highest rare variant enrichment in its subnetwork across all KDs. Our key driver analysis corroborates its role in ASD as a known syndromic rare variant, and its subnetwork also contains both rare (e.g., RELN, NCKAP1) and common (e.g., CACNB3, SV2B) variants, suggesting a regulatory role in both variant types. SYT1 is a synaptic vesicle membrane protein involved in calcium-mediated neurotransmitter release, thus playing a vital role in neuronal communication [48, 49]. It is a known syndromic ASD rare variant causing severe neurodevelopmental abnormalities. The top annotation terms from the SYT1 subnetwork include synaptic transmission, cation channel activity regulation, and nervous system development, agreeing with this gene’s importance in neurodevelopment (Figure 6A).

The subnetworks of several other top brain KDs, Adducin 2 (Figure 6B), SCN8A, AMPH, ATP9A (Supplementary Figure 1A–C), demonstrated enrichment of both variant types, thereby emphasizing the convergence between these variants orchestrated by these brain KDs. For instance, the ADD2 subnetwork had both rare (e.g., SHANK2, SCN8A) and common variants (e.g., ACTN2, DMN1) in its neighborhood (Figure 6B). Adducin genes encode cytoskeleton proteins that are critical for osmotic rigidity and cell shape by regulating the formation of the spectrin-actin membrane skeleton [50]. ADD2 is expressed in the brain, and its knock-out results in the loss of activity-dependent connection formation between neurons [51]. The top functional annotations of the ADD2 subnetwork include synaptic vesicle processes and dendrite extension regulation, which are key components of neurodevelopment (Figure 6B). The function of ADD2 itself and the functional annotations of its subnetwork support its relevance to learning and development and its potential contribution to ASD. SCN8A, a sodium voltage-gated channel gene essential for neuronal excitability, has been associated with neurodevelopmental disorders including ASD [52]. AMPH, a gene involved in synaptic vesicle cycling, is crucial for neurotransmitter release and thus may influence neuronal communication [53]. ATP9A, a phospholipid transporter involved in membrane trafficking and neuronal homeostasis [54, 55], is also a key candidate for further investigation of ASD pathogenesis.

In addition to the top brain KDs discussed above, our analysis also revealed a top KD within the endocrine Bayesian network, PRR36, which showed strong enrichment for common variants (e.g., DNM1, MYT1, MAPK15) but not for rare variants (Supplementary Figure 2A). While its precise role remains understudied, PRR36, or proline rich 36, has been associated with ASD-related genes such as SHANK3, which is implicated in synaptic organization and neurodevelopment [56]. Another top KD, PCDH7 (encoding protocadherin; Supplementary Figure 2B), was identified from digestive tissue networks and its subnetwork is enriched for both rare (SFARI Level 1 gene SETBP1; Level 2 genes CHRM3, KCMNA1, GPC4, STK39, PLN; Level 3 gene PDE3B etc.) and common variants (e.g., MITF, PDE3A, TSPAN2, THRB etc.). PCDH7 is important for cell adhesion, cell-cell interaction and signaling, and brain plasticity and synapse maturation [57, 58]. Its identification as a key driver in the ASD-enriched sigmoid colon module, along with its association with rare ASD variants and its enrichment for cell signaling pathways, suggests shared functions between brain and digestive tissues, with potential contributions to ASD via either or both systems.

DISCUSSION

ASD is a neurodevelopmental disorder that possesses complex heterogeneity in its genetic architecture [17]. Previous research has uncovered both common and rare variants in the pathogenesis of ASD, and have implicated the central nervous system as well as various peripheral tissues and biological pathways that interact to affect ASD [5, 9–16, 59–61]. However, a leading challenge in our understanding of ASD is whether the common and rare variants inform on similar or different mechanisms that can be targeted for therapy. If the former, generalized therapeutic approaches may be developed to target ASD; if the latter, precision approaches are needed to target individual subtypes of ASD.

In our study, we utilize summary statistics from an ASD GWAS [25], tissue-specific eQTL/sQTL data from the GTEx database [29], and tissue-specific gene coexpression and regulatory networks to investigate ASD. Considering both common and rare variants, we sought to understand the molecular interactions in a broad range of brain and peripheral tissues to uncover tissues, genes, biological pathways, and gene networks that are enriched for either common or rare variants of ASD in a completely data-driven manner without using prior knowledge as input. Furthermore, we aimed to predict key drivers within gene regulatory networks that may be of interest for further mechanistic and experimental studies.

From our marker set enrichment analysis (MSEA) of ASD GWAS, we found a diverse range of tissues and gene coexpression modules enriched for ASD common variants. While top enriched modules are mainly from brain tissues such as the anterior cingulate cortex and the amygdala, there were also peripheral tissues whose modules showed relevance to ASD GWAS signals (Supplementary Table 3). The coexpression module enrichments of the anterior cingulate cortex and amygdala are consistent with previous findings that observed structural abnormalities, dysregulations and alterations in both ASD mouse models and adult ASD individuals [62–64]. Other brain regions (frontal cortex [BA9], cerebellum, broader cortical tissue) which had an enrichment of modules for ASD, were previously found to exhibit abnormal development in the disorder [65–67]. Thus, our findings on these brain regions harboring ASD common variants provide additional evidence supporting their importance in cognitive function, motor coordination, and neuronal activity in the context of ASD.

In addition to these confirmatory outcomes of various brain regions, we found that tissues in the digestive, reproductive, endocrine, and immune systems also displayed ASD GWAS enrichment within particular coexpression modules. The digestive system has recently become heavily implicated in ASD pathogenesis as research has continued to explore the interactions within the microbiota-gut-brain axis, and how there are distinctive gastrointestinal complications present in those with ASD [16, 68]. Our findings on the reproductive system support research on how maternal factors and their influence on fetal development contribute to ASD risk [69–74], which can occur via hormones or immune pathways. Our genetics-driven findings on the multi-system involvement highlight the complexity of the genetic architecture across both brain and peripheral regions of the body in the pathogenesis of ASD.

Our exploration of pathway annotations from these coexpression modules enriched for ASD GWAS signals further emphasize how the complexity of ASD spans across diverse biological pathways. We found that immune, cellular signaling, and cell growth and regulation pathways were represented among significant modules from various tissues (Figure 2C; Figure 3). Pathways related to Alzheimer’s disease and Amyotrophic Lateral Sclerosis were consistently enriched among modules from various brain regions, which interestingly aligns with their neurological dysfunction and potential associations with ASD [75–79].

We also noticed that tissues and pathways involved in the endocrine and immune systems were highly represented among the ASD-associated modules. The male and female reproductive, digestive, and adipose systems contained various pathways involved in immune regulation (e.g., Interleukin-6, CD40 pathways), metabolic functions (e.g., cysteine and methionine metabolism, fatty acid synthesis, bile acid metabolism), structural integrity and development (hyaluronan uptake and degradation), or cellular signaling, processes, and regulation (e.g., MAPK pathways, EDG1 pathways, ENOS activation). mTOR signaling, critical for the regulation of metabolism, immune function, and cellular homeostasis [42], is also among the consistent pathways across systems. Our results support the role of immune dysregulation and metabolic dysfunction in ASD pathogenesis [80–84] and implicates the interconnected nature between pathways across systems and their collective contribution to the diversity and complexity of ASD pathogenesis. As a whole, our MSEA findings recapitulate previous studies that highlight the impact of altered neural functioning in various brain regions, and also provide causal inference through genetic variants supporting the roles that digestive, reproductive, endocrine, and immune system have in ASD development.

In our key driver analysis, we discovered that the brain outperforms peripheral tissues in terms of both the abundance of significant rare and common variants among the key drivers identified and the enrichment in their subnetworks for genes with either rare or common variants (Figure 4; Figure 5A–B). As expected, a greater number of brain tissue key driver subnetworks contained a significant number of both rare and common variants, which suggest that these key drivers are crucial regulators within ASD pathogenesis given the convergence of these two types of variants (Figure 5C). The rankings of key driver subnetworks in terms of their rare or common variant enrichment highlight the central role of the brain in ASD, yet also point towards the digestive system as the leading peripheral region for common variant abundance (Figure 5A–B). Our analysis also highlighted top brain key driver genes such as SYT1 (Figure 6A) and ADD2 (Figure 6B). Our detailed analysis of these key driver subnetworks strongly support their potential role in orchestrating genes with either rare or common ASD variants to perform key neuronal functions. While SYT1 has previously been found to be a gene with known rare variant of ASD [39], and a common variant of ADD2 has been implicated in the disorder as well [25], our network analysis highlights the unexplored tissue-specific network interactions between these key drivers and other ASD genes with key functions in neuronal excitability, homeostasis, and neurotransmitter release (e.g., SCN8A, AMPH, ATP9A, shown in Supplementary Figures 1A–C) carrying both types of variants. By contrast, PRR36 (Supplementary Figure 2A), a gene associated with other ASD-related genes such as SHANK3 (implicated in synaptic organization and neurodevelopment [56]) highlights potential ASD mechanisms driven by common variants given the absence of rare variants in its network. Furthermore, the identification of shared KDs between the brain and digestive tissues as seen in PCDH7 (Supplementary Figure 2B), implicates converging mechanisms across tissues that may individually or collectively influence ASD. Overall, these networks provide the opportunity for future experimentation and validation to test how multiple ASD genes interact to affect ASD pathogenesis.

Altogether, our systems genetics approach to studying both ASD common and rare variants across the brain and periphery not only confirmed numerous previous findings but revealed new avenues for future exploration. Each of our top 10 coexpression modules (Supplementary Table 3) captured tissue origins and pathways with known relevance in ASD pathogenesis, with cellular processes and immune pathways being the most replicated across brain and peripheral tissues (Figure 2C; Figure 3). The abundance of both brain and peripheral coexpression modules enriched for ASD GWAS signals (Figure 2B) supports the growing body of research that suggests ASD relevance within not only the brain but the digestive, endocrine, and immune systems [15, 16, 68, 71, 74, 85], and further provides potential causal inference for their involvement due to the incorporation of genetic variants in our analysis. The notable convergence of known rare and common variants in the brain and the identified key drivers regulating ASD genes offer prioritized targets for future mechanistic and therapeutic studies.

Limitations of the study

We note the following limitations in our study. First, while we attempt to use the largest ASD GWAS along with functional genomics datasets such as GTEx with the most comprehensive coverage of tissue-specific gene regulation, the demographics of both the GWAS and the GTEx datasets is primarily of central European origin. ASD research is constantly updating and new discoveries continue to be made, so future studies with expanded population coverage are warranted when such data becomes available. Resources that are actively updated such as the Simons Foundation for Autism Research Initiative [39], the Australia Autism Biobank [86], and the MSSNG database [87] will undoubtedly provide additional useful data for future studies. Secondly, our studies mainly integrated genetic and transcriptomic data to model genetically perturbed gene networks in ASD. Further incorporation of transcriptome, metabolome, and epigenome datasets of ASD [88–90] is necessary to expand our understanding. Third, we also acknowledge that we have only performed a tissue-level analysis. Cell-level analysis is a fast-growing area of research across numerous fields, and ASD is a disease that demands this level of cellular specificity. Lastly, there is a male bias in the ASD prevalence [85], which emphasizes that a sex-specific analysis is likely to reveal additional key insights. Expanding data coverage and performing cell-level, sex-specific analyses are pivotal directions for future research to provide a higher-resolution understanding of ASD’s genetic complexities. Overall, our findings will assist in the generation of new hypotheses and experimental validations to uncover new mechanisms and guide therapeutic development for neurodevelopmental disorders such as ASD.

Supplementary Material

1

This is a list of supplementary files associated with this preprint. Click to download.

• SuppFigure2.tif
• SuppFigure1.tif