Altered expression of Alzheimer's disease-related genes in the cerebellum of autistic patients: a model for disrupted brain connectome and therapy.

Autism and Alzheimer's disease (AD) are, respectively, neurodevelopmental and degenerative diseases with an increasing epidemiological burden. The AD-associated amyloid-β precursor protein-α has been shown to be elevated in severe autism, leading to the 'anabolic hypothesis' of its etiology. Here we performed a focused microarray analysis of genes belonging to NOTCH and WNT signaling cascades, as well as genes related to AD and apoptosis pathways in cerebellar samples from autistic individuals, to provide further evidence for pathological relevance of these cascades for autism. By using the limma package from R and false discovery rate, we demonstrated that 31% (116 out of 374) of the genes belonging to these pathways displayed significant changes in expression (corrected P-values <0.05), with mitochondria-related genes being the most downregulated. We also found upregulation of GRIN1, the channel-forming subunit of NMDA glutamate receptors, and MAP3K1, known activator of the JNK and ERK pathways with anti-apoptotic effect. Expression of PSEN2 (presinilin 2) and APBB1 (or F65) were significantly lower when compared with control samples. Based on these results, we propose a model of NMDA glutamate receptor-mediated ERK activation of α-secretase activity and mitochondrial adaptation to apoptosis that may explain the early brain overgrowth and disruption of synaptic plasticity and connectome in autism. Finally, systems pharmacology analyses of the model that integrates all these genes together (NOWADA) highlighted magnesium (Mg(2+)) and rapamycin as most efficient drugs to target this network model in silico. Their potential therapeutic application, in the context of autism, is therefore discussed.

Autism is a neurodevelopmental disorder characterized by specific activity patterns and aberrant social interaction and communication.[1] Genetic and environmental factors have been proposed to account for apparently increasing rates of autism diagnostics.[2, 3] Its neurobiological basis is, however, poorly understood and the behavioral phenotypes display a high grade of heterogeneity, leading to the definition of ‘autism spectrum disorders' (ASDs). It has been reported that clinical onset of ASD seems to be preceded by a phase of brain overgrowth in the first years of life followed by decreased growth rates between 5 and 12 years of age, when compared with controls.[4] These abnormalities are paralleled by higher levels of brain-derived neurotrophic factor, insulin growth factors 1 and 2 (IGF-1 and IGF-2), IGF-binding protein 3, growth hormone-binding protein, neurotrophin-4, and a neuroinflammatory response in the ASD brains.[5, 6, 7, 8] Changes in trophic support may account for abnormal brain growth and an altered connectome, thus defining multiple ASD symptoms. Recent pieces of evidence also indicate the role for deregulated processing of the Alzheimer's disease (AD)-associated amyloid-β (Aβ) precursor protein (APP) in ASD-associated pathology. Aβ is considered to be central for the pathophysiology of AD, which in its genetic and sporadic forms appears to be the most common cause of age-associated dementia. AD is characterized by specific pathological hallmarks represented by extracellular deposits of fibrillar Aβ, intraneuronal accumulation of neurofibrillary tangles due to hyperphosphorylation of cytoskeletal Tau filaments, as well as impairment of neurogenesis. The end result is massive neuronal death and consequent decline in cognitive functions.[9, 10, 11] In ASD, the AD-associated APP, in particular its neuroprotective processing product or secreted APPα (sAPPα), was found to be elevated.[12] Current research indicates that abnormalities of the cerebellum, now believed to have a role in cognitive function, are associated with autism; being one of the brain regions where cellular and growth abnormalities are most pronounced in the disorder.[13, 14, 15]

Here we address the question of whether changes of expression in the same set of genes observed in AD neuropathology could also be found in autism; and if this is the case, whether drugs tested not only in the context of autism but also in AD may be considered for therapeutics in these patients. Our study shows altered expression of genes belonging to NOTCH, WNT, AD, and apoptosis pathways in cerebellar samples from autistic individuals. In addition, systems pharmacology analyses of such a pathway that integrates NOTCH, WNT, AD, and apoptosis-related genes suggest magnesium (Mg2+) and rapamycin for further therapeutic exploration in the context of autism; two compounds/drugs that therapeutic use have already been discussed for different AD models.[16, 17]

Results

Subnetwork analyses of AD-related pathways reveal transcriptional changes in the cerebellum of autistic individuals

For analyzing the potential transcriptional changes in genes belonging to pathways typically associated to AD pathology in autism, we constructed four subnetworks models of NOTCH (Supplementary Figure S1A), WNT (Supplementary Figure S2A), AD (Figure 1a), and apoptosis (Supplementary Figure S3A) pathways. These subnetworks contained 41, 136, 171, and 70 respective genes/proteins (Supplementary Tables S1–S4). Next, we performed a focused microarray analysis of the genes belonging to each subnetwork by using samples from the cerebellum of autistic patients versus controls. In 31% (116 out of 374) of the genes belonging to these in silico models were differentially expressed in ASD samples (corrected P-value <0.05) (Table 1). More specifically, ∼30% for NOTCH (12 genes), ∼23% for WNT (31 genes), >40% for AD (69 genes), and ∼26% for apoptosis (18 genes). UQCRC1, NDUFS3, NDUFA9, JUN, CIDEB, CDC42, NDUFA6, NDUFB7, NDUFA1, and COX4I1 (Table 1) were the top 10 genes that displayed most significant changes in their expression. The downregulated genes were cytochrome c reductase and NADH dehydrogenase related. Among the AD-related genes, upregulation of GRIN1, the channel-forming subunit of NMDA glutamate receptors, and MAP3K1, known activator of JNK and ERK pathways with anti-apoptotic effect, were also found together with significantly lower expression of presinilin 2 (PSEN2) and APBB1 (Table 1).

Figure 1

AD subnetwork analysis. (a) General landscape of interactions between genes/proteins belonging to the AD pathway (map05010; KEGG Pathway; http://www.genome.jp/kegg/pathway.html). The in silico network model was developed by using the STRING 9.05 database resource search tool, under a confidence score of 0.600 and using ‘Databases' and ‘Experiments' as input options. (b) Focused microarray analyses of AD-related genes in the cerebellum of autistic patients versus healthy controls over the in silico model. The z axis for representing the relative gene expression was constructed by using the ViaComplex software

Table 1

Differentially expressed genes from the NOTCH, WNT, Alzheimer's disease, and apoptosis subnetworks in the cerebellum of autistic patients

Gene symbol	Alias and/or description	Ensembl ID (ENSP)	Differential gene expression (Up/down)	Corrected P-value (FDR)	Subnetwork contribution/s
UQCRC1	Ubiquinol–cytochrome c reductase core protein I	ENSP00000203407	Down	0.001062657	AD
NDUFS3	NADH dehydrogenase (ubiquinone) Fe-S protein 3, 30 kDa (NADH-coenzyme Q reductase)	ENSP00000263774	Down	0.001722868	AD
NDUFA9	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 9, 39 kDa	ENSP00000266544	Down	0.00179604	AD
JUN	Jun oncogene	ENSP00000360266	Up	0.001921098	WNT
CIDEB	Cell death-inducing DFFA-like effector b	ENSP00000258807	Down	0.001949191	APO
CDC42	Cell division cycle 42 (GTP-binding protein, 25 kDa)	ENSP00000314458	Down	0.002000135	WNT
NDUFA6	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 6, 14 kDa	ENSP00000330937	Down	0.002033295	AD
NDUFB7	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 7, 18 kDa	ENSP00000215565	Down	0.002205663	AD
NDUFA1	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 1, 7.5 kDa	ENSP00000360492	Down	0.002205663	AD
COX4I1	Cytochrome c oxidase subunit IV isoform 1	ENSP00000253452	Down	0.002246413	AD
NCSTN	Nicastrin	ENSP00000294785	Up	0.002354256	AD NOT
UQCRC2	Ubiquinol–cytochrome c reductase core protein II	ENSP00000268379	Down	0.002354256	AD
NDUFS4	NADH dehydrogenase (ubiquinone) Fe-S protein 4, 18 kDa (NADH-coenzyme Q reductase)	ENSP00000296684	Down	0.002459308	AD
CYC1	Cytochrome c-1	ENSP00000317159	Down	0.002656014	AD
COX8A	Cytochrome c oxidase subunit 8A (ubiquitous)	ENSP00000321260	Down	0.002656014	AD
NDUFB4	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 4, 15 kDa	ENSP00000184266	Down	0.002688877	AD
EIF2AK2	Eukaryotic translation initiation factor 2-alpha kinase 2	ENSP00000233057	Up	0.002850857	AD
BID	BH3-interacting domain death agonist	ENSP00000318822	Down	0.002892499	AD APO
NDUFA4	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 4, 9 kDa	ENSP00000339720	Down	0.003092192	AD
WNT3	Wingless-type MMTV integration site family, member 3	ENSP00000225512	Up	0.003122302	WNT
NDUFAB1	NADH dehydrogenase (ubiquinone) 1, alpha/beta subcomplex, 1, 8 kDa	ENSP00000007516	Down	0.003262718	AD
NDUFA11	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 11, 14.7 kDa	ENSP00000311740	Down	0.003262718	AD
COX5A	Cytochrome c oxidase subunit Va	ENSP00000317780	Down	0.003530677	AD
PRKCA	Protein kinase C, alpha	ENSP00000284384	Up	0.003532751	APO WNT
NDUFB9	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 9, 22 kDa	ENSP00000276689	Down	0.003532751	AD
NDUFA2	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 2, 8 kDa	ENSP00000252102	Down	0.003622328	AD
COX5B	Cytochrome c oxidase subunit Vb	ENSP00000258424	Down	0.003627494	AD
BCL2	B-cell CLL/lymphoma 2	ENSP00000329623	Up	0.003686548	APO
PRICKLE2	Prickle homolog 2 (Drosophila)	ENSP00000295902	Up	0.003858642	WNT
SDHB	Succinate dehydrogenase complex, subunit B, iron–sulfur (Ip)	ENSP00000364649	Down	0.003949261	AD
NDUFA10	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 10, 42 kDa	ENSP00000252711	Down	0.004023111	AD
NDUFB2	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 2, 8 kDa	ENSP00000247866	Down	0.00417954	AD
NDUFB6	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 6, 17 kDa	ENSP00000369176	Down	0.00420016	AD
NDUFA7	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 7, 14.5 kDa	ENSP00000301457	Down	0.004449939	AD
NDUFB3	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 3, 12 kDa	ENSP00000237889	Down	0.004449939	AD
UQCRB	Ubiquinol–cytochrome c reductase-binding protein	ENSP00000287022	Down	0.004918915	AD
NCOR2	Nuclear receptor co-repressor 2	ENSP00000348551	Up	0.005223194	NOT
LRP5	Low density lipoprotein receptor-related protein 5	ENSP00000294304	Up	0.005717875	WNT
BAD	BCL2-associated agonist of cell death	ENSP00000309103	Down	0.005763042	AD APO
UQCRH	Ubiquinol–cytochrome c reductase hinge protein	ENSP00000309565	Down	0.005906155	AD
PSEN2	Presenilin 2 (Alzheimer disease 4)	ENSP00000355747	Down	0.005968353	AD NOT
NDUFA13	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 13	ENSP00000380364	Down	0.00609761	AD
NDUFB10	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 10, 22 kDa	ENSP00000268668	Down	0.006262982	AD
COX6C	Cytochrome c oxidase subunit VIc	ENSP00000297564	Down	0.006681564	AD
UQCRQ	Ubiquinol–cytochrome c reductase, complex III subunit VII, 9.5 kDa	ENSP00000367934	Down	0.00745614	AD
NDUFA3	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 3, 9 kDa	ENSP00000398290	Down	0.007600014	AD
HEY1	Hairy/enhancer-of-split related with YRPW motif 1	ENSP00000338272	Down	0.00784934	NOT
DAAM1	Disheveled associated activator of morphogenesis 1	ENSP00000247170	Up	0.007909188	WNT
PPP3CC	Protein phosphatase 3 (formerly 2B), catalytic subunit, gamma isoform	ENSP00000240139	Down	0.008118636	AD APO WNT
COX7B	Cytochrome c oxidase subunit VIIb	ENSP00000417656	Down	0.008265168	AD
NDUFB11	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 11, 17.3 kDa	ENSP00000276062	Down	0.008503556	AD
NFKB1	Nuclear factor of kappa light polypeptide gene enhancer in B cells 1	ENSP00000226574	Up	0.008570685	APO
KAT2B	K(lysine) acetyltransferase 2B	ENSP00000263754	Up	0.009111103	NOT
NDUFA12	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 12	ENSP00000330737	Down	0.009479959	AD
UQCRFS1	Ubiquinol–cytochrome c reductase, Rieske iron–sulfur polypeptide 1	ENSP00000306397	Down	0.009524703	AD
NDUFS5	NADH dehydrogenase (ubiquinone) Fe-S protein 5, 15 kDa (NADH-coenzyme Q reductase)	ENSP00000362058	Down	0.009998423	AD
BAX	BCL2-associated X protein	ENSP00000293288	Down	0.010349538	APO
CYCS	Cytochrome c, somatic	ENSP00000307786	Down	0.010843806	AD
NDUFV1	NADH dehydrogenase (ubiquinone) flavoprotein 1, 51 kDa	ENSP00000322450	Down	0.011017543	AD
RBX1	Ring-box 1	ENSP00000216225	Down	0.01185897	WNT
NFKBIA	Nuclear factor of kappa light polypeptide gene enhancer in B cells inhibitor, alpha	ENSP00000216797	Up	0.012804007	AD APO
AXIN2	Axin 2	ENSP00000302625	Up	0.01310984	WNT
CSNK1G1	Casein kinase 1, gamma 1	ENSP00000305777	Up	0.013381084	WNT
NDUFS8	NADH dehydrogenase (ubiquinone) Fe-S protein 8, 23 kDa (NADH-coenzyme Q reductase)	ENSP00000315774	Down	0.014197263	AD
HSD17B10	Hydroxysteroid (17-beta) dehydrogenase 10	ENSP00000168216	Down	0.014254764	AD
RAC3	Ras-related C3 botulinum toxin substrate 3 (rho family, small GTP-binding protein Rac3)	ENSP00000304283	Down	0.014704894	WNT
DTX3	Deltex homolog 3 (Drosophila)	ENSP00000338050	Down	0.015144277	NOT
CDK5	Cyclin-dependent kinase 5	ENSP00000297518	Down	0.015567567	AD
CABIN1	Calcineurin-binding protein 1	ENSP00000263119	Up	0.015953582	AD WNT
COX7C	Cytochrome c oxidase subunit VIIc	ENSP00000247655	Down	0.016686515	AD
RUVBL1	RuvB-like 1 (Escherichia coli)	ENSP00000318297	Down	0.016765972	WNT
DFFA	DNA fragmentation factor, 45 kDa, alpha polypeptide	ENSP00000366237	Up	0.017000844	APO
NDUFB8	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 8, 19 kDa	ENSP00000299166	Down	0.01782766	AD
COX6B1	Cytochrome c oxidase subunit VIb polypeptide 1 (ubiquitous)	ENSP00000246554	Down	0.017845204	AD
TLE1	Transducin-like enhancer-of-split 1 (E(sp1) homolog, Drosophila)	ENSP00000365682	Up	0.018030294	NOT
NOTCH1	Notch homolog 1	ENSP00000277541	Up	0.018222797	AD NOT
RYR2	Ryanodine receptor 2 (cardiac)	ENSP00000355533	Up	0.019085894	AD
CACYBP	Calcyclin-binding protein	ENSP00000356652	Down	0.02050623	WNT
CSNK2B	Casein kinase 2, beta polypeptide	ENSP00000415615	Down	0.02050623	WNT
SENP2	SUMO1/sentrin/SMT3 specific peptidase 2	ENSP00000296257	Up	0.021603075	WNT
CHD8	Chromodomain helicase DNA-binding protein 8	ENSP00000406288	Up	0.021898019	WNT
DTX2	Deltex homolog 2 (Drosophila)	ENSP00000322885	Up	0.022277856	NOT
SDHC	succinate dehydrogenase complex, subunit C, integral membrane protein, 15 kDa	ENSP00000356953	Down	0.023024554	AD
LRP1	Prolow-density lipoprotein receptor-related protein 1 Precursor (LRP)	ENSP00000243077	Up	0.023159082	AD
DLL1	Delta-like 1 (Drosophila)	ENSP00000355718	Up	0.023184387	NOT
APBB1	Amyloid beta (A4) precursor protein binding, family B, member 1 (Fe65)	ENSP00000299402	Down	0.024237411	AD
MAP3K7	Mitogen-activated protein kinase kinase kinase 7	ENSP00000358335	Down	0.02453243	WNT
NDUFC2	NADH dehydrogenase (ubiquinone) 1, subcomplex unknown, 2, 14.5 kDa	ENSP00000281031	Down	0.025489098	AD
APLP1	Amyloid beta (A4) precursor-like protein 1	ENSP00000221891	Up	0.026382742	AD
BIRC2	Baculoviral IAP repeat-containing 2	ENSP00000227758	Down	0.026521867	APO
NDUFB5	NADH dehydrogenase (ubiquinone) 1 beta subcomplex, 5, 16 kDa	ENSP00000259037	Down	0.027921676	AD
PPARD	Peroxisome proliferator-activated receptor delta	ENSP00000310928	Up	0.029747209	WNT
SKP1	S-phase kinase-associated protein 1	ENSP00000231487	Down	0.03098415	WNT
PRKCB	Protein kinase C, beta	ENSP00000305355	Down	0.031488508	WNT
CTNNBIP1	Catenin, beta-interacting protein 1	ENSP00000366466	Down	0.031553064	WNT
RCAN1	Regulator of calcineurin 1	ENSP00000320768	Down	0.031896183	AD
AIFM1	Apoptosis-inducing factor, mitochondrion-associated, 1	ENSP00000287295	Down	0.031972517	APO
CALM2	Calmodulin 2 (phosphorylase kinase, delta)	ENSP00000272298	Down	0.032263251	AD
CAMK2D	Calcium/calmodulin-dependent protein kinase II delta	ENSP00000339740	Up	0.034576294	WNT
SMG7	Smg7 homolog, nonsense-mediated mRNA decay factor (Caenorhabditis elegans)	ENSP00000340766	Up	0.034586995	APO
TBL1X	Transducin (beta)-like 1X-linked	ENSP00000217964	Up	0.036511455	WNT
CASP3	Caspase 3, apoptosis-related cysteine peptidase	ENSP00000311032	Up	0.037829198	AD APO
FRAT2	Frequently rearranged in advanced T-cell lymphomas 2	ENSP00000360058	Up	0.037835094	WNT
TCF7L2	Transcription factor 7-like 2 (T-cell specific, HMG-box)	ENSP00000358404	Up	0.03801246	WNT
PRICKLE1	Prickle homolog 1 (Drosophila)	ENSP00000345064	Up	0.038401928	WNT
MAP3K1	Mitogen-activated protein kinase kinase kinase 1	ENSP00000382423	Up	0.038406592	AD
MAML1	Mastermind-like 1 (Drosophila)	ENSP00000292599	Up	0.038709571	NOT
GRIN1	Glutamate receptor, ionotropic, N-methyl D-aspartate 1	ENSP00000360616	Up	0.03916101	AD
SNCA	Synuclein, alpha (non A4 component of amyloid precursor)	ENSP00000338345	Down	0.040205801	AD
IKBKG	Inhibitor of kappa light polypeptide gene enhancer in B cells, kinase gamma	ENSP00000358622	Down	0.04087009	APO
PRKCZ	Protein kinase C, zeta	ENSP00000367830	Down	0.041055303	WNT
EP300	MicroRNA 1281	ENSP00000263253	Up	0.041373019	NOT WNT
ENDOG	Endonuclease G	ENSP00000361725	Down	0.042737537	APO
PPP3CB	Protein phosphatase 3 (formerly 2B), catalytic subunit, beta isoform	ENSP00000378306	Down	0.046120071	AD APO WNT
BCL2L1	BCL2-like 1	ENSP00000302564	Up	0.046646273	APO
NDUFV2	NADH dehydrogenase (ubiquinone) flavoprotein 2, 24 kDa	ENSP00000327268	Down	0.046908433	AD

Corrected P-values <0.05 were considered significant. ‘NOT', ‘WNT', ‘AD', and ‘APO' represent NOTCH, WNT, Alzheimer's disease, and apoptosis subnetworks, respectively

The in silico models were subsequently subjected to further analysis with the ViaComplex software that plots the expression mean values derived from each diseased sample (autism) over the expression of healthy control samples (z axis) in the network model (NOTCH, WNT, AD, and apoptosis), presenting them as two-dimensional (2D) ViaComplex-generated landscapes (Supplementary Figures S1B–S3B; and Figure 1b). The generated plot visualizes the general landscape of gene expression of diseased samples (autism) versus healthy controls for each subnetwork. Expression levels are color coded with warm colors (yellow to red) being indicative of increased expression, whereas cold (green to blue) colors indicate decreased gene expression when compared with control samples. The 2D ViaComplex plots revealed a general upregulation of NOTCH pathway (Supplementary Figure S1B) and specific decrease in the expression of mitochondria-related genes (network cluster on the upper left; Figure 1b).

NOWADA is an integrative network model for NOTCH, WNT, AD, and apoptosis subnetworks

As results presented above showed significant changes in expression of genes belonging to NOTCH, WNT, AD, and apoptosis pathways in cerebellar brain biopsies from ASD patients, our next aim was to test the feasibility of an in silico model that is able to integrate these four subnetworks into one single network and therefore to characterize the relevant molecular route. The Venn diagram, which visualizes the level of subnetwork information integration (Figure 2), shows that the four subnetworks contain common nodes (genes/proteins) allowing communication with one another. Then, by using STRING 9.05 (‘Experiments' and ‘Databases' confidence score of 0.600) and plotting with Cytoscape software, we developed a NOTCH-WNT-Alzheimer's Disease-Apoptosis (NOWADA) gene/protein interaction network model able to characterize in silico these molecular interactions (Figure 3). The model is composed of 374 genes/proteins interconnecting through 3665 interactions, with only five nodes (PPP3CA, PPP3CB, PPP3CC, PPP3R1, and PSEN1) contributing to maximum three groups (NOTCH, WNT, AD, or apoptosis subnetworks) at the same time.

Figure 2

Venn diagram showing information integration from the studied subnetworks. NOTCH, WNT, AD, and apoptosis models contain common nodes (genes/proteins) with, at least, one subnetwork

Figure 3

Network model of interactions between genes/protein belonging to NOTCH, WNT, AD, and apoptosis subnetworks (NOWADA network). The present model was developed by using the STRING 9.05 database resource search tool, under a confidence score of 0.600, and by using ‘Databases' and ‘Experiments' as input options and visualized by plotting it with Cytoscape software. The number of subnetwork contributions of each gene/protein within the network is represented in the figure as indicated in the inset

Identification of hey hub genes/proteins within NOWADA network in the context of autism

In order to describe the global characteristics of the network, we measured centralities[18] or topological network properties as follows: (i) ‘stress', representing how much a node (gene/protein) is traversed by a high number of ideal routes or short paths in a network. Then, a gene/protein traversed by higher number of short paths will be, by definition, more stressed. (ii) ‘connectivity or degree' that quantifies the local topology of each node by summing up the number of its adjacent nodes. (iii) ‘Betweenness', which is similar to stress as a topological network property, provides a more informative and elaborated centrality index as it measures how frequently the shortest path, connecting every pair of nodes, is going through a third given node. In other words, both ‘stress' and ‘betweenness' give information about the influence of a node over the spread of information throughout the network. (iv) Finally ‘closeness', defined by the inverse of the average length of the shortest paths to access all other proteins in the network, gives an idea about the level of proximity of a node to other nodes or how long it will take the information to spread from a given node to other nodes in the network. Logically, the larger the value, the faster the information spreads. Our focus in this part of the study is to target vulnerable (central) components in NOWADA network with 31% of its genes displaying significant changes in expression in the cerebellum of autistic patients (corrected P-values <0.05). Upon calculation of each centrality value for each gene/protein (Supplementary Table S5), we identified the hub-bottlenecks (HBs) of the network that also displayed higher values of ‘stress' and closeness centralities. These network nodes were (in alphabetical order) APP, CALM1, CASP3, CTNNB1, DVL2, DVL3, EP300, FZD5, GSK3B, MYC, NOTCH1, PLCB2, PRKCA, and TP53 (Supplementary Table S6 and Figure 4). As central members, these genes/proteins will control the flow of biological information within the network and its disruption (e.g., by drug interaction) could destroy the entire network into small components. Remarkably, CTNNB1, gene encoding for β-catenin protein (Supplementary Table S2), was the node with top values for all the measured centralities (Supplementary Table S6 and Figure 4).

Figure 4

Analysis of the topological properties (a–c) of genes/proteins belonging to the network (NOWADA). Dashed lines are indicating the threshold value for each property. Upregulated and downregulated genes in the subnetwork are marked with red and blue triangles, respectively, as indicated in the inset. Note that NH-Bs and HBs distinguish non-HBs from HBs, respectively

Systems pharmacology analysis of NOWADA network suggests Mg2+ and rapamycin as most efficient drugs to target the model in silico

Previous analysis allowed us to identify central nodes (bottlenecks) with high topological network values. At this stage, as NOWADA network can be considered as a common pathway model for AD and autism, these 14 nodes (genes/proteins) could represent potential targets for both pathologies in silico. With the aim of elucidating the most suitable and efficient drugs targeting the network in the context of autism, differential gene expression values were taken into account. Our selection criteria was then as follows: (i) to select drugs targeting higher number of genes/proteins within the network; (ii) from those, to choose the drugs affecting/associated to higher numbers of bottlenecks (preferably, HBs with high values of ‘stress' and ‘closeness' centralities); and (iii) to select drugs that targets, within the network (NOWADA), include higher number of genes with significant changes in expression. Upon systematic review of the literature (see Materials and Methods section), 47 used/proposed drugs for autism therapy were divided in 15 groups and tested for potential interaction/s with genes/proteins of the network model. The groups are acetylcholinesterase inhibitors, allosteric modulators of metabotropic glutamate receptors, antidepressants (SSRI and tricyclic), antiepileptics/anticonvulsants, antiopioids, antipsychotics, complementary and alternative medicines, inhibitors of mTOR, neurohormones, NMDAR (agonists and antagonists), mood stabilizers, psycoestimulants, and sympatholytic medications. Incidentally, 36 out of the 47 autism-related drugs have already been studied in different models of AD (in humans or in vivo) for potential therapeutic use[19, 20, 21, 22] (Supplementary Table S7). Upon using STITCH 3.1, we confirmed that 16 out of the 47 tested drugs had positive hits with at least one gene/protein of NOWADA network (with ‘Experiments' and ‘Databases' as input options; confidence score of 0.600). These therapeutic agents were (in alphabetical order) as follows: chlorpromazine, clozapine, D-cycloserine, fluoxetine, galantamine, haloperidol, imipramine, lithium, luteolin, Mg2+, melatonin, memantine, olanzapine, rapamycin, valproic acid, and ziprasidone. Most targeted nodes were GSK3B (glycogen synthase kinase 3β, 6 hits), AKT1 (5 hits), CALM1 (3 hits), GRIN1 (3 hits), GRIN2B (3 hits), CHRNA7 (2 hits), GRIN2A (2 hits), and JUN (2 hits) (Supplementary Table S9). Except for CHRNA7, our systematic review of the literature could confirm the expression of the indicated gene products in neuronal and/or glial cells (astroglia, microglia, and oligodendroglia; Supplementary Table S9), thus representing both neuronal and glial targets. When looking into the characteristics of the tested compounds/drugs–genes/protein interactions, Mg2+ and rapamycin are highlighted as the most efficient drugs in silico to target this disrupted pathway, which is here represented in silico by NOWADA network (Supplementary Table S10). More specifically, Mg2+ interacts with 32 genes/proteins; among those, eight differentially expressed genes in the cerebellum of autistic patients: the upregulated EIF2AK2 (corrected P-value=0.002850857) and CAMK2D (corrected P-value=0.034576294), as well as the downregulated CDC42 (corrected P-value=0.002000135), COX5A (corrected P-value=0.003530677), NDUFV1 (corrected P-value=0.011017543), NDUFS8 (corrected P-value=0.014197263), RAC3 (corrected P-value=0.014704894), and CALM2 (corrected P-value=0.032263251) (Figure 5). In order to elucidate the biological processes that Mg2+ could be influencing in this context, a functional enrichment analysis of the genes belonging to the Mg2+ subnetwork was performed. The top 10 affected processes are as follows: phosphorylation, phosphorus metabolic process or phosphate-containing compound metabolic process, neuromuscular process, regulation of neurological system process, protein amino-acid autophosphorylation, small GTPase-mediated signal transduction, mitochondrial electron transport (NADH to ubiquinone), positive regulation of cell projection organization, and intracellular signaling cascade (corrected P-value <0.05; Supplementary Table S11).

Figure 5

Mg2+ subnetwork analysis. (a) In silico network model of the interactions between Mg2+ and genes/proteins belonging to the network (NOWADA), developed by using STITCH 3.1 with ‘Experiments' and ‘Databases' as input options and a confidence score of 0.600. (b) Actions view representation between the network nodes by ‘catalysis' (yellow). (c) Actions view representation between the network nodes by ‘binding' (white color). (d) Actions view representation between the network nodes by ‘inhibition' (red color). Upregulated and downregulated genes in the subnetwork are marked with red and blue triangles, respectively, as indicated in the inset

Discussion

Here we studied the possible neurobiological mechanism involved in the excessive rates of early brain overgrowth observed in autism[4] based on the recent evidence, indicating a possible role for the AD-associated APP processing pathway in ASD.[23] Proliferation in the brain results from balancing proliferative rates and developmental apoptosis (i.e., more proliferation and less death will result in a larger brain).[24] APP is proteolytically processed by two competing pathways such as the amyloidogenic (β-secretase-mediated) and the non-amyloidogenic (α-secretase-mediated) pathways.[25] sAPPβ (secreted APPβ) and the neurotoxic Aβ peptide are generated from APP and released to the extracellular site by proteolytic cleavage of both β- and γ-secretases.[26] Aβ has been suggested to have a critical role in neurotoxicity and apoptotic events seen in AD.[27] In contrast, the non-amyloidogenic product sAPPα and the p3 peptide are formed by proteolytic cleavage of both α- and γ-secretases[26] sAPPα exerts proliferative effects on neural progenitor cells isolated from embryonic brains and is able to enhance synaptogenesis, neurite outgrowth, cell survival, and cell adhesion;[28, 29, 30, 31] and APP seems to be able to modulate the neuronal precursors migration.[32] As γ-secretase cleaves several proteins that include APP and Notch,[33] the potential role of NOTCH signaling pathway in the genesis of AD is also considered.[34] Focused microarray analysis used in this study demonstrated that expression of 40% of the AD-related genes was altered in the cerebellum of autistic patients (corrected P-value <0.05), these changes mainly represented by downregulation of NADH dehydrogenases and cytochrome c oxidases (Table 1). These results reflect deficient mitochondrial respiratory chain that is consistent with recent reports describing mitochondrial dysfunction in ASD;[35] or may indirectly confirm tight Notch–mitochondria interconnections.[36] An upregulation of GRIN1 (corrected P-value=0.03916101), the channel-forming subunit NR1 of NMDA glutamate receptors and MAP3K1 (corrected P-value=0.038406592), also known as MEKK1 or MAPK/ERK Kinase Kinase 1, which activates JNK and ERK pathways with anti-apoptotic effect,[37] were also found among the AD-related genes in ASD samples (Table 1). In vivo studies have shown that stressful events during gestation can increase the expression of NMDA receptors in the brain.[38] The altered expression of NMDA receptors may be related to epigenetic regulation (e.g., DNA hypomethylation) because histone lysine methylation at gene promoters is involved in developmental regulation and maintenance of region-specific expression patterns of ionotropic and metabotropic glutamate receptors in the brain.[39] Furthermore, glutamate receptors expression and/or their activity can modulate the APP processing mechanisms.[40] High glutamate receptors activation and consequent increase in intracellular calcium activate the ERK signaling cascade, enhancing α-secretase cleavage of APP, and thus reducing APP processing into Aβ. An increase in α-secretase activity may lead to higher production of sAPPα.[26] In this study, we hypothesize that deregulated glutamatergic synaptic transmission/plasticity, caused by increased expression of GRIN1 and hence increased density of NMDA receptors, may affect neuronal growth in autism by shifting APP processing to the production of sAPPα through ERK-mediated α-secretase activity. SAPPα, in turn, may facilitate proliferation by activating the PI3K/Akt/mTOR pathway.[41, 42] Based on the present data and existing literature, we proposed this model (Figure 6) as a mechanism that could explain the focal modification of neurogenesis, migration, and alterations of the cytoarchitecture of brain cortex, subcortical structures and cerebellum observed in individuals with autism.[43] Consequently, the decrease in expression of mitochondrial enzymes observed in our array (Table 1) may represent a compensatory mechanism to contain Ca2+ overload-induced apoptosis in these cells. It has already been reported that compromised mitochondria accumulate (observed by increased mitochondrial mass) in pyramidal neurons in temporal cortex of younger ASD subjects.[44] Moreover, PSEN1 and PSEN2 are believed to have almost overlapping cellular functions and contribute to diverse physiological processes, such as regulation of Ca2+ homeostasis, protein transport and turnover, autophagy, cell adhesion, neurotransmitter release, and axon guidance.[45] PSEN2, however, seems to be responsible for endoplasmic reticulum (ER) and mitochondria interactions.[46] Closer ER–mitochondria juxtaposition (regulated by PSEN2) may expose mitochondria to excessive Ca2+ stimulation, triggering the apoptotic cascade by mitochondria Ca2+ overload. In our array analysis, we found that PSEN2 is downregulated (corrected P-value=0.005968353) in the cerebellum of autistic patients (Table 1), which could represent another compensatory mechanism to reduce Ca2+ overload-induced apoptosis. Our hypothesis is further supported by the following observations: (i) mutant mice in which glutamate receptors are overstimulated by knocking out (KO) glutamate transporters GLAST and GLT1, and thus leading to excessive glutamatergic signaling in the prenatal stage, compromising early brain development via overstimulation of NMDARs. NR1 deletion in double KO mice almost completely rescued multiple brain defects including cortical, hippocampal, and olfactory bulb disorganization and defective corticothalamic and thalamocortical axonal projections;[47] (ii) opposite to what it is seen in AD patients, higher levels of sAPPα have been described in plasma from autistic children when compared with controls, favoring an increased α-secretase pathway;[12, 48] (iii) evidence of autoimmunity and persistent systemic immune activation have been reported in ASD,[49] and mice overexpressing human sAPPα indeed exhibit higher T-cell cytokine secretion together with reduced CD4+ and higher CD8+ T-cell populations in splenocytes when compared with wild-type animals. These animals also display hypoactivity, impaired social interaction, elevated levels of brain glial fibrillary acidic protein (GFAP) expression, and altered Notch1 and interleukin-6 levels.[50] In the brain of autistic individuals, increased expression of GFAP in the areas with disturbed neuronal architecture has been detected together with interleukin-6, interleukin-8, tumor necrosis factor-α, interferon-γ, and granulocyte–macrophage colony-stimulating factor, suggesting astrogliotic response and potential alterations in neurogenesis and neuronal migration.[51, 52, 53, 54] (iv) Finally, we observed a significant decrease (corrected P-value=0.024237411) in expression of APBB1 (Aβ A4 precursor protein binding, family B, member 1; or F65) (Table 1). In vitro F65 overexpression induces a marked increase in Aβ secretion, whereas Aβ secretion is decreased in Fe65 knockdown cells and also in hippocampal neurons of Fe65/Fe65L1 KO mice.[55, 56, 57]

Figure 6

Hypothetical model of the synaptic mechanism regulating Aβ production and the link to aberrant brain growth and connectivity in autism. The development of living organisms parallels to a so-called physiological ‘intrinsic developmental stress', which is associated to massive internal changes during morphogenesis, and thus contributes to the time window of vulnerability to the environment. Development and severity of ASD are likely to arise from complex interactions between pre-existing genetic vulnerability/ies in these children, the exposure to noxious environmental factors, and the timing of the stressful event(s); given that prenatal life, infancy, childhood, and adolescence are critical periods characterized by increased vulnerability to stressors.[2, 3] Glial activation and neuroinflammation seem to persist during adulthood in autism, and early brain overgrowth in autism could be a consequence of an over-activation of neural proliferation in a trial to compensate the cellular loss induced by environmental stress/ors, triggering precipitous/uncontrolled migration and misdistribution of neurons among different brain areas. Our model proposed that aberrant epigenetic regulation may lead to increased density of NMDA receptors and, therefore, to increased Ca+2 entry and stimulation of ERK-dependent α-secretase activity. The decrease in mitochondrial enzymatic activity together with the downregulation of PSEN2 at the ER may represent a compensatory mechanism to reduce the Ca2+ overload-induced apoptosis. PSEN2 downregulation may diminish the interaction between ER and mitochondria, reducing its Ca+2 uptake. In turn, mitochondrial mass would increase most likely to maintain the cellular ATP levels. Higher levels of sAPPα would activate the PI3K/Akt/mTOR pathway resulting in proliferation, aberrant brain growth, and disruption of synaptic plasticity and connectome

Systems pharmacology analyses of NOWADA network revealed that Mg2+ and rapamycin could efficiently target this pathway in silico, when compared to other drugs (Supplementary Table S10). RhoGAPs, regulatory molecules of RHO GTPases (e.g., CDC42, RHOA, and RAC1), which altered expression were shown in autism,[3] use Mg2+ as a cofactor to reach catalytic efficiency and specificity in GTP hydrolysis (Supplementary Table S11).[58] Furthermore, GSK3B, a HB in our network model (Supplementary Table S6), is a ubiquitous serine/threonine kinase that catalyzes the transfer of γ-phosphate of ATP to the hydroxyl oxygen of the Ser and Thr residues of a kinase-specific protein substrate is a Mg2+-dependent kinase.[59] Its aberrant function has been linked to diverse neurological disorders such as AD.[60] Lower concentrations of Mg2+ were described in children with autism.[61] Similarly, in AD patients, Mg2+ concentrations in the brain seem to be significantly lower when compared with age-matched normal individuals.[62] Recently, it was suggested that elevation of brain Mg2+ concentration in the brain exerts substantial synaptoprotective effects in an in vivo model of AD, and hence Mg2+ may have therapeutic potential for treating AD in humans.[17] In autism, positive behavioral effects of combined Mg2+ and pyridoxine treatment were described;[63] combination that was also proposed to be beneficial in ASD patients that experience seizures.[64] The therapeutic usage of rapamycin (mTOR inhibitor) in ASD is nowadays discussed. TSC1 and TSC2 are upstream inhibitory regulators of mTOR activity, and studies performed in TSC2+/− mice showed that treatment with rapamycin can rescue synaptic plasticity and reverse the abnormal late phase of long-term potentiation.[20] These results encourage future clinical trials with TOR inhibitors as pharmacological treatment of ASD.

We believe that the computational analyses performed in this study gave valuable clues to (i) develop a model that may explain the aberrant brain overgrowth and impaired neuronal connectivity observed in these children; and (ii) refer to the need of further therapeutic exploration of Mg2+- and rapamycin-based treatments in animal models and clinical trials with autistic patients.

Materials and Methods

Development of gene/protein network models

NOTCH, WNT, AD, and apoptosis subnetwork models were developed through (i) extracting the information provided by the KEGG PATHWAY database (http://www.genome.jp/kegg/pathway.html; NOTCH pathway: map04330, WNT pathway: map 04310; AD pathway: map05010, and apoptosis pathway: map04210); and (ii) by establishing the interactions between genes/proteins with the database resource Search Tool for the Retrieval of Interacting Genes/Proteins STRING version 9.05 (http://string-db.org/); with ‘Experiments' and ‘Databases' as input options and a confidence score of 0.600. STRING is a well-known public database with information about direct and indirect functional protein–protein interactions. The genes/proteins belonging to NOTCH, WNT, AD, and apoptosis subnetwork models were identified by the HUGO Gene Symbol[65] and Ensemble protein ID (Supplementary Tables S1–S4). Once they were selected, the links between two different nodes (genes/proteins), provided by the STRING database, are saved in data files to be processed in the Medusa interface[66] and visualized as a subnetwork (NOTCH, WNT, AD, or apoptosis). The Venn diagram was constructed by using the freely available software system R (http://www.r-project.org)[67] in order to visualize the grade of molecular relation (common genes/proteins) between the generated subnetworks. The in silico network model integrating NOTCH, WNT, AD, and apoptosis subnetworks (NOWADA) was developed by using again STRING version 9.05 with ‘Experiments' and ‘Databases' as input options and a confidence score of 0.600 (0.400 is considered ‘medium confidence') and plotted by using Cytoscape, an open source platform for complex network analysis and visualization (http://www.cytoscape.org/).[68]

Data acquisition and processing, determination of variably expressed genes, statistics, and gene expression network visualization

Microarray data (raw data;.cel files) were obtained from the Geo DataSets database (http://www.ncbi.nlm.nih.gov/geo/). The data set (GSE38322), originally contributed by Ginsberg et al.,[69] is publicly available; it contains data from cerebellar brain tissue of autistic patients and control subjects. Experimental assays are described in full in the original publication and details include the selection criteria of patients (both inclusion and exclusion criteria) and diagnostic profile.[69] For elucidating the differential expression of members from the NOTCH, WNT, AD, and apoptosis subnetworks, expression data were filtered from probes with <0.05 signal detection P-values and normalized by using the lumi package from R and robust spline normalization. For differential gene expression analysis, normalized data of cerebellar samples from patients versus controls were analyzed by using the limma package from R[67] and false discovery rate (FDR) control for statistical assessment of the microarray data (corrected P-values <0.05 were considered significant). Specific corrected P-values for each differentially expressed gene are provided (Table 1). In order to visualize the general landscape of gene expression in silico, gene expression in samples from autism were plotted versus the expression found in healthy controls to be visualized with the ViaComplex software, a tool that was previously developed and validated in our laboratory.[70] This software plots gene expression over the Medusa network (in a color scale) by overlapping functional input data (microarray expression data) with interaction information (NOTCH, WNT, AD, and apoptosis subnetworks) and distributes the microarray signal according to the coordinates of the nodes and its interactions in the network of interest.

Elucidation of network centralities to predict the relevance of genes in the overall network architecture

Topological network properties of the genes/proteins integrating the network (NOWADA) were calculated by using Cytoscape.[68] Numeric values corresponding to the properties of each node are additionally provided (Supplementary Table S5). Thresholds for selection criteria of the most relevant nodes according to their centrality values were established considering one s.d. of the mean for each topological network properties: ‘stress': 48 919.44 (mean=15251.9); ‘connectivity or degree': 35.84815 (mean=19.59893); ‘betweenness': 0.022270086 (mean=0.005765439); and ‘closeness': 0.372488254 (mean=0.325097663). Centrality values were plotted in graphs by using OriginLab (http://www.originlab.com/index.aspx?go=Products/OriginPro). Nodes (genes/proteins) with high connectivity (above the threshold) and low betweenness centrality (below the threshold) are considered hub-non-bottlenecks and those ones with a value above the thresholds for both connectivity and betweenness centrality are referred to as HBs. Similarly, non-HBs (NH-Bs) are nodes with low connectivity and high betweenness centrality, and non-hub-non-bottlenecks are nodes with both connectivity and betweenness centrality values below the thresholds.[18]

In silico development of compound (drug)-gene/protein network model

To construct the compound (drug)-gene/protein network models, we first performed a systematic review of relevant literature associated with autism or AD therapy in the PubMed database (http://www.ncbi.nlm.nih.gov/pubmed/). These articles were obtained by using the two terms ‘autism' and ‘ASD' and combining them with the following terms: ‘Alzheimer's disease', ‘therapy', ‘therapeutics', ‘treatment', and/or ‘drug', respectively. Moreover, compound identifier, formula, description, and selected autism- and AD-related reports (with PMID) for these drugs (Supplementary Table S6) are additionally provided with KEGG DRUG (http://www.genome.jp/kegg/drug/), the Search Tool for Interactions of Chemicals STITCH version 3.1 (http://stitch.embl.de/), and PubMed database as sources of information. Second, these drugs were submitted to in silico analysis by searching the number of genes/proteins within NOWADA network that are targeted by the selected drugs. To achieve this goal, STITCH version 3.1 (http://stitch.embl.de/) was used with ‘Experiments' and ‘Databases' as input options and a confidence score of 0.600 (0.400 is considered ‘medium confidence'). A list with the most drug-targeted network nodes (genes/proteins) is also provided (Supplementary Table S8). Again, the links between two different nodes (drug-genes/proteins), provided by the STITCH database, is saved in data files to be handled in the Medusa interface and visualized as a subnetwork. For screening potential cellular targets of the most drug-targeted network nodes (NOWADA), another systematic review of both original research articles and reviews was performed in PubMed database. These articles were obtained by using the genes/proteins name or symbol and combining them with the following terms: ‘Neuronal', ‘Neuron', ‘Astroglia', ‘Astrocyte', ‘Microglia', ‘Oligodendroglia', and/or ‘Oligodendrocyte' (Supplementary Table S9).The selected references confirm the expression of the indicated gene products in brain cells. Note that the selected references do not cover therapeutic modalities, age or brain region specificity.

Functional enrichment analysis (biological processes)

The biological processes of the genes products belonging to the subnetwork model of interactions was determined by using the database for annotation, visualization, and integrated discovery DAVID v6.7 (http://david.abcc.ncifcrf.gov/), giving a number of functional annotation tools to researchers for better comprehension of the biological meaning behind any large list of genes. Only those biological processes with corrected P-value <0.05 (FDR) were selected. Specific corrected P-values for each biological process are provided (Supplementary Table S9).

A graphical abstract, summarizing the contents and the methodological approaches used for the present study, is also provided (Figure 7).

Figure 7

Abstract workflow summarizing the different approaches used for the present study. The information from the different pathways (NOTCH, WNT, AD, and apoptosis) was gathered from KEGG pathway (http://www.genome.jp/kegg/pathway.html) to develop the subnetworks. The differential expression of the genes belonging to each subnetwork was analyzed in the cerebellum of autistic patients by using R. A Venn diagram showed common genes/proteins between NOTCH, WNT, AD, and apoptosis pathways. Thereafter, a network model was developed for integrating these subnetworks into one single gene/protein interaction network model (NOWADA network) by using STRING 9.05 and Cytoscape tools. Further elucidation of the topological network properties revealed a number of potential molecular targets or ‘vulnerable' points of the in silico model and therapeutic agents used in autism targeting central nodes of the network were studied (STITCH 3.1). Finally, a subnetwork representing the interaction between Mg2+ and genes/proteins from NOWADA network was constructed and a functional enrichment analysis was performed (DAVID v6.7) to elucidate the biological processes potentially affected by this compound