Blood Transcriptomic Mega-Analysis of Autism Spectrum Disorder – Supplemental Material – Page 21

Supplementary Material

Table of Contents

Supplementary Results

Sex-Specific Analyses and Cross-Sex Comparisons Pages 2-3

Supplementary Discussion Pages 4-9

References Pages 10-17

Supplementary Figures Pages 23-35

Supplementary Figure 1 Page 18-20

Supplementary Figure 2 Page 21-23

Supplementary Figure 3 Page 24-28

Supplementary Figure 4 Page 29

Supplementary Figure 5 Page 30-31

Supplementary Results

Sex-Specific Analyses and Cross-Sex Comparisons

When data were analyzed separately for females (nASD = 90, ncomparison = 117), 19,508 gene-level mixed-effect models produced valid results and 1,609 genes showed at least nominal evidence for dysregulation (uncorrected p < 0.05), with 802 over-expressed and 807 under-expressed in ASD, and none reached a corrected level of significance (FDR q < 0.05). Among males (nASD = 536, ncomparison = 330), 17,268 gene-level models produced valid statistical results and 2939 showed at least nominal evidence for dysregulation (uncorrected p < 0.05), with 1338 over-expressed and 1601 under-expressed in ASD, and 594 reached a corrected level of significance (FDR q < 0.05). When examining the intersection of nominally dysregulated gene lists (Supplementary Figure 4; center of Venn diagram), and accounting for the common background of genes analyzed, we observed significant cross-sex overlap of 489 dysregulated (hypergeometric p < 5.3x10-61), 233 over-expressed (p < 1.4x10-76) and 240 under-expressed genes (p < 4.7x10-64). Interestingly, we observed a few genes showing apparently sex-discordant ASD-related effects (Supplementary Figure 4). Also notably, when we compared the absolute values of the differences of least squared means among the top 1% of dysregulated genes for each sex, we observed larger magnitude ASD-related mean differences among females (mean |differencef| = 0.48 + 0.06, mean |differencem| = 0.30 + 0.03; p < 6.4x10-73) and larger F-values among males (mean Ff = 10.9 + 2.4, mean Fm = 18.4 + 4.6; p < 4.2x10-60); these differences could potentially be explained by differences and sample size and relative statistical power. Full lists of sex-specific summary statistics will be provided upon request to the corresponding author.

Gene-set analysis was performed using the single-gene statistics derived from each sex (Supplementary Figure 4; side panels). Among the 7,350 sets assessed for females, 1,397 reached a liberal threshold of statistical significance (family-wise BH correction) and 368 reached a conservative threshold (family-wise Bonferroni correction). Among males, 7,358 gene-sets were assessed and 817 and 243 met liberal and conservative thresholds, respectively. We compared the results across sexes with respect to the common background of interrogated gene-sets, in order to identify transcriptomic effects that appear to be conserved across sexes; expectedly, many effects observed in the main analysis showed evidence for preservation in both females and males (Supplementary Figure 4; center of Venn diagram). We also highlight gene-sets that reached a conservative threshold for directional dysregulation in one sex, but that showed no evidence for dysregulation (either conservative or liberally-defined) in the opposite sex (Supplementary Figure 1; side panels); we further verified that no biologically similar (or semantically interchangeable) terms show evidence for dysregulation in the opposite sex. Potentially female-specific findings included over-expression of genes involved in DNA repair and under-expression of genes involved in heme metabolism, IL-2 signaling via STAT5, IL-15 signaling, transcriptional targets of MEK and miRNA species, and several chromosomal cytobands (Supplementary Figure 4, left panel). Potentially male-specific findings included over-expression of genes involved in O-linked glycosylation and transcriptional targets of MEK (Supplementary Figure 3, right panel). Full lists of sex-specific gene-set summary statistics will be provided upon request to the corresponding author.

We repeated network-level analyses separately for females (nASD = 90, ncomparison = 117) and males (nASD= 536, ncomparison = 330). Among females, modules were bi-directionally highly preserved between ASD cases and comparison subjects (two modules showed moderate preservation; Supplementary Figure 1B). When female cases and comparison subjects were analyzed together, we identified forty-one network modules; linear mixed models identified 8 as nominally associated with ASD, though none survived BH correction for multiple testing (Supplementary Figure 5, left side). Among males, modules are also bi-directionally highly preserved (Supplementary Figure 1C). When cases and comparison subjects were analyzed together, we identified thirty-two modules, 8 of which were nominally associated with ASD, and 4 survived BH correction for multiple testing (Supplementary Figure 5, right side). On the whole, ASD-associated modules among females (composed of 1236 genes) shared a significant number of genes (299, hypergeometric p < 6.0x10-75) with those implicated among males (composed of 1208 genes). Two ASD-associated modules among females (dark turquoise and midnight blue) did not show greater-than-chance over-representation of genes that were implicated in either the single-gene differential expression analysis or network analysis among males (Supplementary Figure 5, pink bracket). The dark turquoise module was diminished in affected females and was enriched with erythroblast markers and genes involved in heme metabolism, among other functions. The midnight blue module was enhanced in affected females and was enriched with CD4+ and CD+ T-cell-related functions, estrogen-responsive genes, cell proliferation and WNT-beta-catenin signaling. Among males, four modules that were diminished in ASD cases (cyan, dark turquoise, midnight blue, dark turquoise, and yellow) did not show greater-than-chance over-representation of genes that were implicated among female samples (Supplementary Figure 5, blue bracket). Notably, the midnight blue module was enriched with genes involved in neurotrophic signaling, cell adhesion, and B-cell-related functions, and showed a trend toward enrichment with GWAS signal (permuted p < .09). The dark turquoise module was also enriched with B-cell related functions and transcriptional targets of NKX2. The cyan module was enriched with genes involved in cell cycle, estrogen responses, MHCII antigen presentation, and several transcriptional targets, as well as GWAS signal (permuted p < 0.04). The large yellow module was enriched with genes involved in innate immune signaling functions, neurodevelopmental processes, neurotrophic signaling, glycerophospholipid metabolism, and numerous transcriptional targets. Full lists of sex-specific summary statistics and module enrichments will be provided upon request to the corresponding author.


Supplementary Discussion

We sought to understand our findings (focusing on functional gene-sets and cell types) in the context of relevant literature. As alluded to in the main text discussion, we found several instances where the literature describes higher levels of circulating signaling protein (e.g., cytokine, trophic factor) in ASD, yet we observed reduced expression of leukocyte mRNAs that encode intracellular proteins that serve in the protein’s intracellular second messenger or regulatory cascade and/or a group of mRNAs reflect that are transcriptionally stimulated by the circulating protein’s signaling cascade. We also found instances where ASD-associated transcriptomic differences in blood appear to be opposite those previously described in post-mortem brain tissue. We have attempted to depict these discrepancies, focusing on studies of human tissues, in the main text Figure 3. The studies that support this depiction, as well as other pertinent information obtained from model systems are discussed in subsequent sections. These discrepancies are curious and beg further investigation. Future mechanistic studies that examining ASD cases and comparison subject blood samples could simultaneously measure circulating signal protein, intracellular second messenger protein and signaling activation, and intracellular expression of key transcriptional targets in order to help confirm the present observation and shed light on their mechanisms. Similarly, within subject studies examining both blood and brain tissues could help confirm whether opposing effects exist across tissue at the level of the transcriptome, proteome and phosphorylome. We conclude by highlighting literature supporting the rationale for explicitly considering sex in biological studies of ASD.

RAS-MEK-MAPK and ERK1 & 2 Signaling

As discussed in the main text, we observed reduced expression of the RAS-MEK-MAPK and ERK1 and 2 signaling cascades and their transcriptional targets in blood samples; these effects were apparent in both sexes. Relatively few studies have reported on the transcriptional evidence for this cascade’s activation. For context, it is worth noting that EGFR-family and platelet-derived growth factors receptors (PDGFRs) exert effects through this second messenger cascade (Roberts and Der, 2007; Navas et al., 2012; Klinghoffer et al., 1996; Gilbertson et al., 2006; Satoh et al., 1993), and these cascades appear to play critical roles in survival for some neurons (Wagner et al., 2006; Botella et al., 2003; Read et al., 2009; Gilbertson et al., 2006), and influence behavior (Van Buskirk and Sternberg, 2007). Mutations resulting in hyper-activity of the RAS-MAPK signaling axis are well-established causes of combined cranio-facial-neurodevelopmental syndromes (e.g. neurofibromatosis type 1, Noonan Syndrome, LEOPARD Syndrome, Costello Syndrome, and cardiofaciocutaneous syndrome) known to have a high incidence of ASD (Sanchez-Ortiz et al., 2014; Alfieri et al., 2014; Stornetta and Zhu, 2011). There’s reason to suspect that over- or under-active RAS signaling could contribute to neurodevelopmental abnormalities (Stornetta and Zhu, 2011). In neuronal cells, these cascades play roles in cell proliferation, migration, differentiation, and death. A brief review of the neurbiological roles of these signaling axes are available in Yang et al.’s introduction section (Yang et al., 2013). Microdeletions and microduplications that included ERK1 have been identified in individuals with ASDs (Weiss et al., 2008; Kumar et al., 2008). RAS/RAF/ERK1/2 signaling activities were reported to be significantly up-regulated in the frontal cortex of autistic individuals (though theses paper were recently retracted; Yang et al., 2011; Zou et al., 2011) and in the BTBR murine model of autism (Yang et al., 2013; Faridar et al., 2014; Yin et al., 2014), with one study further suggesting coordinate over-expression of p-ERK in both prefrontal cortex and lymphocytes of an animal model (Faridar et al., 2014). Studies of cultured neurons indicate that RAF/ERK up-regulation stimulates the migration of cortical and cerebellar granular cells, and impairs the formation of excitatory synapses in the former cells. RAF/ERK up-regulation also suppressed the development and maturation of cortical neurons(Yang et al., 2013). In animal models of Fragile X Syndrome, heightened RAS signaling appears to undermine glutamate receptor 1-mediated long term potentiation; this effect could be restored with up-regulation of RAS-PI3k-AKT signaling (Hu et al., 2008). Animal models involving a deletion of region syntenic to the human 16p11.2 region (harboring ERK1) show a host of behavioral abnormalities that parallel aspects of ASD and attention deficit-hyperactivity disorder (Horev et al., 2011; Pucilowska et al., 2015; Portmann et al., 2014; Yang et al., 2015), as well as abnormalities in the development of a variety of brain structures and cortical cytoarchitecture. Remarkably, one model organism for this microdeletion showed reduce the amount of ERK1 protein (with ERK2 protein unchanged from control levels), but also hyper-phosphorylation (activation) of both ERK1 and ERK2 relative to controls, suggesting that a complex regulatory mechanism might result in paradoxical hyperactivity of this signaling pathway (Pucilowska et al., 2015). Another recent study indicated that pharmacological blockage of MEK activity at specific developmental time points can also produce ASD-like behavioral abnormalities (Yufune et al., 2015). Notably, a conditional knock-out model of ERK2 also featured abnormalities in social behavior (Satoh et al., 2011). Distal 22q11.2 microdeletions affecting the gene encoding ERK2 (among others) have been reported in association with cognitive impairment and a DiGeoge-syndrome like phenotype (Shaikh et al., 2007; Samuels et al., 2009). Thus, there is substantial evidence that a controlled range of signaling activity through this axis may be important for typical brain development.

PI3K-AKT-MTOR Signaling

We observed under-expressions of gene-sets related to PI3K-AKT-MTOR signaling in females and males. Another recent microarray-based study found nominally significant enrichment of MTOR-related genes among those showing significant evidence of alternative transcript splicing in ASD blood samples (Stamova et al., 2013). It’s worth noting that EGF (through interactions with EGFR and the MET receptor) are capable of signaling via AKT. One recent study found lower levels of total p-AKT in leukocytes from children with ASD; p-AKT levels were negatively correlated with leukocyte EGFR and Hepatocyte Growth Factor levels (Russo, 2015). A study of lymphocytes from individuals with Fragile X Syndrome (another genetic syndrome with high rates of ASD) identified no differences in RNA levels of translational control genes like MTOR, AKT, and EIF4E, but they observed higher levels of p-AKT and higher levels of MTOR-specific phosphorylation of p70 S6K1, suggesting increased activity of the signaling axis (Hoeffer et al., 2012). For broader context, there is considerable literature supporting the interactions of the RAS-MAPK and PI3K-AKT signaling systems in controlling cell survival (Castellano and Downward, 2011; Aksamitiene et al., 2012; Lee et al., 2006; Mendoza et al., 2011) and these systems also appear to interact to influence neuronal morphology (Kumar et al., 2005). With respect to the AKT-MTOR signaling axis, they are thought to play a role in complex signaling integration, cell survival, and autophagy, as well as the coupling of neurotransmission and protein synthesis for plasticity (Hoeffer and Klann, 2010). There is strong evidence that mutations leading to the over-activation of this signaling axis are associated with neuroanatomical phenotypes and ASDs in humans (Crino, 2011). However, studies of post-mortem brain from idiopathic ASD cases (and studies of brain from valproic acid-exposed rats) suggests reduced expression and activation of AKT-MTOR signaling proteins (including reduced expression of p706k and eIF4B proteins) may be present (Nicolini et al., 2015). Chen et al. discuss evidence implicating dysregulated insulin-like growth factors (IGF-I) upstream of AKT-MTOR in ASD (Chen et al., 2014). Another bioinformatic analysis identified TNF-α and β-estradiol as potential regulators of gene networks related to AKT-MTOR signaling (Lee et al., 2012). A study of ncRNAs differentially expressed in post-mortem brain of ASD cases suggested that AKT-MTOR signaling mediators were enriched among ncRNA regulatory targets (Ander et al., 2015).

Epithelial Growth Factor Receptor (EGFR) and ERBB2 Signaling

We observed a main effect consistent with reduced EGFR and ERBB2 (a.k.a. HER2/neu, human epidermal growth factor receptor 2) signaling is ASD blood samples. Our search did not identify studies reporting transcriptional differences explicitly related to this signaling axis in ASD. Circulating EGF concentrations have been reported as higher in ASD (Tobiasova et al., 2011; İşeri et al., 2010), while others reported reduced EGF in ASD samples (Manzardo et al., 2012; Russo, 2013; Suzuki et al., 2007). Another recent study found lower levels of plasma MIG6 protein, a negative regulator of EGF signaling, in ASD samples (Russo, 2014). One microarray study of cortical tissue identified EGF signaling as an enriched function among differentially expressed genes, though they did not shed light on the direction of this effect (Garbett et al., 2008). Genetic studies have shown a haplotypic association involving the locus of the EGF ligand (Toyoda et al., 2007). Notably, the neuregulin (NRG) gene-family members are ligands for members within the EGFR family; neuregulin 3 (NRG3, ligand for ERBB4) is implicated as a schizophrenia risk locus by common variant association studies (Kao et al., 2010). An animal model of NRG3 over-expression during development resulted in adult phenotype characterized by increased anxiety and diminished social behavior (Paterson and Law, 2014). More generally, EGFR signaling is involved in cell growth, proliferation, and differentiation; it also plays roles in the developing and adult CNS (Page, 2003; Libermann et al., 1984; Mazzoni and Kenigsberg, 1992). EGFRs exert signaling effects through MAPK, PI3K/AKT, and STAT signaling cascades.