A novel computational biostatistics approach implies impaired dephosphorylation of growth factor receptors as associated with severity of autism

Abstract

The prevalence of autism spectrum disorders (ASDs) has increased 20-fold over the past 50 years to >1% of US children. Although twin studies attest to a high degree of heritability, the genetic risk factors are still poorly understood. We analyzed data from two independent populations using u-statistics for genetically structured wide-locus data and added data from unrelated controls to explore epistasis. To account for systematic, but disease-unrelated differences in (non-randomized) genome-wide association studies (GWAS), a correlation between P-values and minor allele frequency with low granularity data and for conducting multiple tests in overlapping genetic regions, we present a novel study-specific criterion for 'genome-wide significance'. From recent results in a comorbid disease, childhood absence epilepsy, we had hypothesized that axonal guidance and calcium signaling are involved in autism as well. Enrichment of the results in both studies with related genes confirms this hypothesis. Additional ASD-specific variations identified in this study suggest protracted growth factor signaling as causing more severe forms of ASD. Another cluster of related genes suggests chloride and potassium ion channels as additional ASD-specific drug targets. The involvement of growth factors suggests the time of accelerated neuronal growth and pruning at 9-24 months of age as the period during which treatment with ion channel modulators would be most effective in preventing progression to more severe forms of autism. By extension, the same computational biostatistics approach could yield profound insights into the etiology of many common diseases from the genetic data collected over the last decade.

Figure 1

μGWAS QR plot (curved) vs traditional ssGWAS QQ plot (straight), left: AGP I, right: AGP II. Each point represents the most significant result among all diplotypes centered at the same SNP ranked by significance (low to high). Dashed blue curve: projection. Solid blue curve: loess estimation (see Materials and Methods). Vertical lines connect the most significant s-values (−log₁₀ P) of a gene (dot) with its expected value (solid blue line). Light red and gray vertical lines indicate genes with unknown function and results with low reliability (either low μIC or reliance on a single SNP), respectively. Top and bottom gene list (by significance, right to left, excluding genes with unknown function): μGWAS and ssGWAS results, respectively. Shaded genes are among the genes highlighted in Figure 2. Full and open triangles mark genes with an identical match or family member of SFARI genes, respectively (see Supplementary Table 1 for details). The dotted horizontal lines represent the projected WG apex (6.272 and 6.064) and an exploratory 100 gene cutoff (4.835 and 4.480) for AGP I and AGP II, respectively (Supplementary Table 1). The horizontal solid blue line indicates the proposed study-specific GWS.

Figure 2

Ras/Ca²⁺ signaling in ASD and CAE. (a) AGP I, (b) AGP II, (c) CAE. Pathway-related genes among the top 20, 50 and 100 are circled in bold, double and thin blue lines, respectively. Genes included in SFARI Gene (ASD) and CarpeDB (CAE), respectively, are shaded in red (see Supplementary Table 1 for details); the five genes identified in previous GWAS (see Introduction) are indicated in turquoise and underlined. Upon GF binding to cell-surface receptors (for example, IGFR, MET, PDGFR, ERBBn), formation of receptor complexes initiates proliferation, cytoskeletal organization and survival along Ras downstream effectors. GFs are immediately deactivated by PTPRs. Downstream activities are modulated by agonists binding to G-protein-coupled receptors (GPCR) activating phospholipase C (PLC) to form membrane diacylglycerol (DAG) and inositol trisphosphate (IP₃). While DAG activates Ras directly, IP₃ stimulates (‘winged' arrows) the release of Ca²⁺ from the endoplasmic reticulum (ER), starting a process of Ca²⁺-dependent activation of Ras involving several feedback loops. The fall of Ca²⁺ concentration in internal stores (dotted areas) leads, via STIM1, to the opening of store-operated Ca²⁺ channels (SOCC) in the plasma membrane. ITPKB phosphorylates IP₃ into IP₄, which opens voltage-operated Ca²⁺ channels (VOCC). CaCCs can either directly activated by Ca²⁺ elevation or through Ca²⁺/calmodulin kinase II (CaMKII)-mediated phosphorylation. Other plasma membrane ion channels involved are Ca²⁺ channels operated by NMDA and kainate ligands, voltage-operated potassium channels (VOPC). GABA-operated Cl⁻ channels reverse from excitatory efflux to inhibitory influx during maturation. Overall Ca²⁺ levels are limited by plasma membrane Ca²⁺ ATPase (PMCA). Known drug interactions are indicated in green.

Figure 3

Extended Manhattan plot of μGWAS results for PTPRT (top), PPFIBP1 (bottom/left), and PTPRB (bottom/right) by AGP stage. The X-axis shows base pairs within chromosome. Black dots indicate significance in ssGWAS, lines indicate significance in diplotypes of width 2 (dotted) .. 6 (solid), red color indicates low μ-scores for reliability, suggesting a potential artifact, unless supported in both populations. Green dots and s-values indicate univariate results for SNPs within the most significant region. s-Values in gray indicate nearby results. Below the panels are gene annotations, LD blocks, and recombination rate from HapMap. The PTPRT region comprises rs6102794, rs6072693, rs6072694, rs6102795, rs6016759 and rs6102798. The ‘x' and box at the bottom indicate a somatic mutation at rs146825584 and a deletion at 41,036,259–41,300,521 (Supplementary Table S2, AU018704), respectively. The PTPRB: region comprises rs3782377, rs2567137, rs2567133, rs2278342, rs2116209 and rs2278341. KCNMB4 results driven by a single SNP in one population only (rs787931, red ‘x') are indicated as a potential artifact, but the related KCNMA1 was the most significant gene in AGP II (Table 1).

Figure 4

Comparison of HFA cases against all parental controls. Genes shaded in dark and light blue are members of or associated with the Ras/Ca²⁺ pathway (Figure 2), respectively (see Figure 1 for legend). PFN2 inhibits the formation of IP₃ and DAG by binding to PIP2, LIMK1 is phosporylated by ROCK1 and PAK1, downstream of RHOA and RAC1, respectively. PTENP1 acts as a decoy for PTEN-targeting miRNAs. SLC25A21 may be involved in 2-oxodipate acidemia, which is accompanied by mental retardation and learning disabilities. Cytogenic bands: ANO2: 12p13.3, ANO4: 12q23.3.

Figure 5

Comparison of PTPRT allelotype profiles between SDA cases, melanoma controls, and HFA cases. Rows indicating individual subjects' SNP profile (orange/green: homozygous; yellow: heterozygous) are sorted within each population by diplotype μ-score (dark green to dark read) computed from the three consensus SNPs (rs6102794, rs6072694 and rs6102795 out of the six-SNP PTPRT region of Figure 3), which are highlighted as more saturated. Dotted lines are added for visual guidance.

Figure 6

Hypothesized interventions to prevent regression in children with ASD. During the critical period of developing cortical structures for social interactions the risk of stress-induced regression might be reduced through a combination of strategies including a protective environment where exposure to unfamiliar faces is limited and pharmacological interventions to reduce hyperexcitability related to Ca²⁺ signaling by targeting ion channels determined through genetic testing of genes known to be involved in Ca²⁺ signaling among children with a risk phenotype.