Revealing the Molecular Mechanisms of Alzheimer’s Disease Based on Network Analysis

The complex pathology of Alzheimer’s disease (AD) emphasises the need for comprehensive modelling of the disease, which may lead to the development of efficient treatment strategies. To address this challenge, we analysed transcriptome data of post-mortem human brain samples of healthy elders and individuals with late-onset AD from the Religious Orders Study and Rush Memory and Aging Project (ROSMAP) and Mayo Clinic (MayoRNAseq) studies in the AMP-AD consortium. In this context, we conducted several bioinformatics and systems medicine analyses including the construction of AD-specific co-expression networks and genome-scale metabolic modelling of the brain in AD patients to identify key genes, metabolites and pathways involved in the progression of AD. We identified AMIGO1 and GRPRASP2 as examples of commonly altered marker genes in AD patients. Moreover, we found alterations in energy metabolism, represented by reduced oxidative phosphorylation and ATPase activity, as well as the depletion of hexanoyl-CoA, pentanoyl-CoA, (2E)-hexenoyl-CoA and numerous other unsaturated fatty acids in the brain. We also observed that neuroprotective metabolites (e.g., vitamins, retinoids and unsaturated fatty acids) tend to be depleted in the AD brain, while neurotoxic metabolites (e.g., β-alanine, bilirubin) were more abundant. In summary, we systematically revealed the key genes and pathways related to the progression of AD, gained insight into the crucial mechanisms of AD and identified some possible targets that could be used in the treatment of AD.


Introduction
Alzheimer's disease (AD) is the most frequently diagnosed neurodegenerative disorder worldwide and the sixth leading cause of death in the United States and is on the rise [1]. The disorder is characterised by amyloid plaques and neurofibrillary tangles; cell loss, vascular damage and dementia follow as a direct result of the vicious cycle of their deposition (Figure 1) [2]. Along with age and family history, inheritance also plays an essential role in the development of AD [3]. Among several genetic risk factors, the APP, PSEN1, PSEN2 had been identified as a causative factor for early-onset AD, the APOE-ε4 allele (encodes a protein that transports cholesterol in the bloodstream) was shown to have a strong impact on late-onset AD. However, studies based on APOE status among different racial and ethnic groups have shown inconsistent results [1,4,5]. Despite a number of treatments being approved by the U.S. Food and Drug Administration (FDA), Systems biology-based methodologies have proven effective in providing biological insight into genetic elements and in identifying the most promising biomarkers associated with complex diseases [8][9][10]. For instance, weighted gene co-expression network analysis (WGCNA) was successfully used to construct gene networks in various diseases and identify centrally connected hub genes as promising biomarkers or therapeutic targets [11,12]. Additionally, the use of WGCNA can provide a deeper understanding of which functional regulators may be driving transcriptional signatures in the development of disease, such as transcription factors (TFs) [13]. Moreover, GEnome-scale metabolic Models (GEMs) provide quantitative information on how the different metabolites are linked to each gene and each reaction in the network [8,14,15]. GEMs condense all information about the known associations of protein-encoding genes and how these genes/proteins interact with other bioactive compounds and associated reactions.
In the present study, we first conducted differential expression analysis (DEA) and gene set enrichment analyses (GSEA) to collect baseline biological information. Second, we constructed gene co-expression networks for each brain region and searched for biologically relevant modules to identify key genes. We evaluated the importance of these key genes with respect to literature, baseline biological information and protein-protein interactions (PPIs). Finally, we reconstructed AD-specific GEMs and performed analysis to depict reporter metabolites and reveal the changes in metabolic pathways. A workflow for this study is presented in Figure 2.
In the present study, we first conducted differential expression analysis (DEA) and gene set enrichment analyses (GSEA) to collect baseline biological information. Second, we constructed gene co-expression networks for each brain region and searched for biologically relevant modules to identify key genes. We evaluated the importance of these key genes with respect to literature, baseline biological information and protein-protein interactions (PPIs). Finally, we reconstructed AD-specific GEMs and performed analysis to depict reporter metabolites and reveal the changes in metabolic pathways. A workflow for this study is presented in Figure 2.

Distinct Gene Expression and Functional Profiles in Different Brain Regions of AD
We performed DEA to characterise regional gene expression changes in AD compared to control. We found 2885 differentially expressed genes (DEGs) (1472 up, 1413 down) in dorsolateral later prefrontal cortex (DLPFC), 477 DEGs (205 up, 272 down) in temporal cortex (TCX) and 1515 DEGs (944 up, 571 down) in cerebellum (CBE) regions (Supplementary Dataset S2). The DEGs between each pair of regions were significantly overlapped as showed by the significance of shared DEGs (Table 1), and most of these genes were changed in the same direction (Supplementary Dataset S2). Furthermore, 34 genes were found to be significantly differentially expressed in all three datasets. Interestingly, 33 of these 34 genes were found to be expressed in the same direction and involved in many pathways ranging from secretion to signalling and to cancer-associated pathways ( Table 2). We then performed Kyoto Encyclopedia of Genes and Genomes (KEGG) [16] pathways enrichment analysis for the 34 genes as well as overlapped genes between each pair of DEGs from these three datasets to investigate their functional relevance. As a re-

Distinct Gene Expression and Functional Profiles in Different Brain Regions of AD
We performed DEA to characterise regional gene expression changes in AD compared to control. We found 2885 differentially expressed genes (DEGs) (1472 up, 1413 down) in dorsolateral later prefrontal cortex (DLPFC), 477 DEGs (205 up, 272 down) in temporal cortex (TCX) and 1515 DEGs (944 up, 571 down) in cerebellum (CBE) regions (Supplementary Dataset S2). The DEGs between each pair of regions were significantly overlapped as showed by the significance of shared DEGs (Table 1), and most of these genes were changed in the same direction (Supplementary Dataset S2). Furthermore, 34 genes were found to be significantly differentially expressed in all three datasets. Interestingly, 33 of these 34 genes were found to be expressed in the same direction and involved in many pathways ranging from secretion to signalling and to cancer-associated pathways ( Table 2). We then performed Kyoto Encyclopedia of Genes and Genomes (KEGG) [16] pathways enrichment analysis for the 34 genes as well as overlapped genes between each pair of DEGs from these three datasets to investigate their functional relevance. As a result, we found that these overlapped DEGs are over-represented in many cancer-associated pathways, signalling pathways, cell differentiation and apoptosis (Table 3), suggesting a potential link to cell senescence and cell death.   We also investigated the functional profiles of DEGs from each brain region using the KEGG pathway and Gene Ontology (GO) [17] term enrichment analysis ( Figure 3). Most generally, we found that the cancer pathways were enriched in DEGs of all three regions, while anabolic reactions, e.g., oxidative phosphorylation and amino acid biosynthesis, and synaptic activity tended to diminish. As we mentioned above, all shared DEGs were associated with many pathways enriched with DEGs from brain tissues ( Figure S5a,b). The reduction in GABAergic and glutamatergic synapse pathways in CBE is consistent with our knowledge of the GABAergic and glutamatergic neuron loss in AD. We observed reduced oxidative phosphorylation, small molecule metabolism (pyruvate, butanoate, propanoate, beta-alanine, fatty acid) and amino acid degradation (valine, leucine, isoleucine, lysine) that reflects crucial AD-associated changes in brain metabolism. Supported by the literature, the enrichment of HIV-1 infection [18,19] and hepatitis B infection [20] associated pathways, NF-κβ [21,22], VEGF [23] and toll-like receptor signalling [24] pathways, diminishment of peroxisome [25,26] and morphine addiction [27] associated pathways in CBE was found among interesting AD-associated abnormalities. We also observed a reduction in the retrograde endocannabinoid signalling pathway [28] and nicotine addiction [29]  In parallel, DLPFC and CBE showed abundant GO term enrichment. Cytoskeleton organisation (e.g., actin cytoskeleton organisation), GTPase activity, cell membrane-linked signalling pathways and associated processes (e.g., G protein-coupled receptor signalling pathway, neurotransmitter transport, glutamate receptor signalling pathway) and angiogenesis (e.g., blood vessel morphogenesis, epithelial cell differentiation) were enriched significantly in DLPFC data. We also found the enrichment for cytoskeleton, neuronal regions and respective molecular activity (e.g., microtubule-binding, cadherin binding, Rho GTPase binding) ( Figure S6a-c). These enriched pathways indicate changes supporting cholinergic signalling and vascularisation. The enriched pathways in the CBE dataset were similar but scarce compared to DLPFC results ( Figure S6d-f). GO enrichment analysis did not yield any significant results for TCX data.

Co-Expression Network Analysis in Different Brain Regions of AD
We also constructed region-specific co-expression networks to investigate functional gene modules (See methods). We identified one large module for DLPFC (node size = 2034), one large module for TCX (node size = 2422) and 11 modules for CBE (total node size = 574) complying given size and connectivity criteria. Even though the centrality measurement was not as dramatic as in DLPFC and TCX modules (Supplementary Dataset S4), the 11 modules for CBE were merged as one module gene group for analyses hereafter as they have high connectivity (Supplementary Dataset S5). We also observed significant and consistent overlap between modules from different datasets (Supplementary Dataset S5). Then we investigated the association between the DEGs that we identified in different brain regions and the genes involved in each module. We found that the identified modules were overlapped with the down-regulated DEGs from their respective region, suggesting that down-regulated genes in these regions follow a similar regulation pattern. We also generated a merged set of down-regulated genes (n = 1975) which were significant in at least one dataset based on DEA for studying AD changes in general. We found that each module contains nearly one-third of all identified down-regulated DEGs (609 in DLPFC module, 798 in TCX module, 184 in merged CBE module) and a modest amount of these genes were hub genes in the modules (86 in DLPFC module, 89 in TCX module, 17 in merged CBE module). Modules that we determined were then annotated for a list of numerous biological processes and molecular functions, ranging from amino acid biosynthesis to vesicular transport. We classified these annotations into energy metabolism and synaptic activity as shown in the figures (Figure 4, Figure S7a-i).   M8 (n = 2034), the single module from DLPFC, and M11 (n = 2422), the single module from TCX, were enriched for nearly all metabolic pathways partly due to their size. For instance, genes involved in vitamin, glycan and leukotriene metabolisms were abundant specifically for M8 and M11. Nevertheless, oxidative phosphorylation was the only significant enrichment for both modules (hypergeometric p-value 0.00034 and 0.00054). Aminoacyl t-RNA metabolism was enriched significantly (hypergeometric p-value < 0.0476) for M8 genes. Whereas some CBE modules (M57, M90, M195 and M244) were not enriched for any given annotation, others shared genes associated with synaptic activity and energy metabolism.
We then investigated the status of genes related to amyloid and tau hypotheses in modules and differential expression analysis. These are comprised of genes responsible for amyloid precursor protein (APP) synthesis and catabolism down to Aβ (i.e., APP, BACE1, PSEN1, PSEN2, NCSTN, ADH1 and ADAM10) and tau synthesis (i.e., MAPT). Only two genes, APP and PSEN2, were found in our modules, DLPFC and TCX, and they were down-regulated in these datasets. Despite not being in the modules, significantly changed expression of MAPT (down) and ADAM10 (up) in DLPFC were noteworthy.
We detected two hub genes shared in DLPFC, TCX and CBE modules: GPRASP2 and AMIGO1. GPRASP2 is a G-protein coupled receptor (GPCR) regulator. It facilitates endocytosis of GPCRs from the plasma membrane to lysosomes for degradation [30]. AMIGO1 is important for the growth of neurites and may contribute to the myelination of neural axons [31]. AMIGO1 and GPRASP1, which come from the same protein family of GPRASP2, are significantly down-regulated in DLPFC. Notably, GPRASP2 is also down-regulated with bolder line significance (adjusted p-value in DLPFC = 0.0568). According to Brain Atlas data on Human Protein Atlas (HPA) [32], GPRASP2 has relatively high expression levels in brain tissues overall compared to other tissues (https://www.proteinatlas.org/ENSG000001583 01-GPRASP2/tissue. Accessed 5 May 2021). There is a 10-15% decrease in the expression of GPRASP2 in all datasets. We also observed 5% down-regulation in DLPFC, a 10% decrease in CBE and nearly no change in TCX for AMIGO1. Still, these small changes may be critical considering that AMIGO1's highest expression is in the cerebral cortex and cerebellum among all tissues (https://www.proteinatlas.org/ENSG00000181754-AMIGO1/tissue. Accessed 5 May 2021).
Moreover, other shared hub genes among two datasets include PAK1 in CBE and DLPFC, ATP6V1E1, CDC42, DCTN2, ERLEC1 and MPP1 in CBE and TCX, and 86 other genes between DLPFC and TCX. These genes are mostly associated with glucose-dependent energy metabolism, ubiquitin-proteasome system (UPS), synaptic activity and plasticity, cytoskeleton organisation as well as intra-golgi and retrograde golgi to endoplasmic reticulum traffic.

Alterations in Energy Metabolism, Chaperones and Synaptic Activity
Based on the network analysis and DEA results, we found a wide change in energy metabolism, chaperones and synaptic activity. Regarding the energy metabolism, we found a batch of ATPase genes that are responsible for ATP catalysis and solute carrier family genes which are responsible for the transport of ATP and its substrates through membranes were affected. ATP1B3 encodes non-catalytic β-3 subunit of ATPase Na + /K + transporting enzyme, which is one of the shared DEGs and expressed 15-40% higher in all AD tissues. ATPAF1, which encodes ATP synthase mitochondrial F1 complex assembly factor, and 15 ATPase cation transporting subunit encoding genes (mostly H + , one Na +/ K + and one plasma membrane Ca +2 transporting), were significantly down-regulated in DLPFC. ATP1B2, ATP1A3, ATP13A4 and ATP6V0C were significantly up-regulated in CBE. There were no significant expression changes in ATPases in TCX other than ATP1B3. The expression of the hub genes ATP6V1A and ATP6V1B2 shared by DLPFC and TCX, and ATP6V1E1 shared by DLPFC and CBE were not significantly changed. Solute carrier family member genes SLC9A6, SLC25A4 and SLC25A14 were shared by the DLPFC and TCX, but their expression was not significantly changed (Table 4). In addition, we found some solute carrier genes showed different expression patterns associated with AD. SLC25A30, SLC7A2 and SLCO4A1 were up-regulated in DLPFC and TCX; SLC6A12, SLC38A2 and SLC9A3R2 were up-regulated and SLC16A6 was down-regulated in DLPFC and CBE; SLC9A3 was down-regulated in TCX and CBE (Supplementary Dataset S2). These results suggest that the expression of many genes associated with energy production were changed in AD and these genes were mostly up-regulated. Chaperone proteins, also known as heat shock proteins (HSPs), are crucial for proteostasis. Chaperones can be divided into two large families in eukaryotes: the Hsp70 family (gene symbols: HSPs) and the Hsp40 family (gene symbols: DNAJs). According to our findings, HSPA12A and DNAJC5 showed higher centrality in DLPFC and TCX than other HSPs (Table 4). Seven members of the DNAJ C subfamily were significantly down-regulated in DLPFC but not in TCX. The other three DNAJ genes were significantly up-regulated in CBE. Seven members of HSPs had changes in DLPFC, whereas three other HSP genes were significantly up-regulated in TCX. There are few interactions between HSPs and other proteins. DNAJC30 encodes a mitochondrial protein that is enriched in neurons and regulates mitochondrial respiration. This protein interacts with a mitochondrial rRNA methyltransferase encoded by MRM1, which has lower expression in DLPFC AD samples. It was reported that MRM1 catalysed methyltransferase which is important for the mitoribosome assembly and synthesis of electron transport chain submits [33]. In addition, our findings show us an interaction between DNAJC30 and RNF170 proteins. RNF170 was down-regulated in all tissues, especially in TCX. A strong association of RNF170 and active inositol 1,4,5-trisphosphate (IP3) receptors-the former mediates the ubiquitination of the latter-was shown in an animal proteomics study [34]. IP3 receptors are located on the endoplasmic reticulum as calcium channels and are degraded by UPS upon activation. HSP90B1 interacts with SGTB on the protein level, which is a co-chaperone and regulates Hsp70 ATPase activity. Both HSP90B1 and SGTB were down-regulated in DLPFC.
We also observed a wide dysregulation of synapse-associated genes in AD. The results showed that more than a third of synaptic genes were differentially altered (mostly down-regulated) in DLPFC and CBE tissues for each synapse type. Of 27 genes shared in all synapse types, two genes, GNG12 and GNGT2, were up-regulated in DLPFC, while seven genes, namely GNG2, GNG3, GNG4, GNB5, CACNA1D, PRKACA and PRKCG, were down-regulated. GNB4 and CACNA1D were down-regulated in CBE. We found eight synaptic activity-associated genes (AP2M1, ARHGEF9, CALM3, DLG3, GABRB3, L1CAM, NPTN and NSF) which interact with other DEGs on protein level and these genes were shared in DLPFC and TCX modules. Most of these genes were down-regulated in DLPFC except for L1CAM and NSF. AP2M1 encodes a subunit of assembly protein complex 2, which is crucial for the acidification of endosomes and lysosomes through proton pumping. CALM3 encodes calmodulin, which is critical for synaptic activity since it mediates voltage-dependent calcium channels [35,36]. Calmodulin also takes part in a protein kinase complex, CaMKIIα, and is critical for learning and memory [37]. The protein encodes by DLG3 interacts with NMDA receptor subunits at the synapse; thus, it is required for synaptic plasticity and learning [38]. General down-regulation of synaptic activity associated genes both in DLPFC and CBE implies that the cholinergic hypothesis is more plausible to explain AD in our case. CREB3, a member of cholinergic and dopaminergic synapses is significantly up-regulated in TCX.

Metabolic Alterations in Different Brain Regions of AD
To investigate the metabolic alterations in different regions in AD, we performed the GEM analyses. First, we reconstructed the region-specific GEM based on the transcriptomic profiles of each region in the brain. Numbers of reactions, metabolites and genes for each GEM are given in Table 5. Region-specific GEM construction allows us to understand the changes in the activity of metabolites, reactions, subsystems/pathways and metabolic functions/tasks. In terms of reactions, all groups showed a great similarity (≥95%); however, groups from the same brain region compared to other groups had a higher similarity. Group correlations based on gene expressions supported this finding ( Figure 5). DLPFC groups were more distinct from others, as we observed in previous analyses. We analysed the metabolic tasks ( Figure 6) and the differentially active pathways (Figure 7) by comparing GEMs to enhance our understanding of functional changes in AD from plausible metabolite perturbations. Metabolic tasks were mostly different in DLPFC. Lactosylceramide de novo biosynthesis (a class of glycosphingolipid) and bilirubin conjugation failed in the AD DLPFC group, while they were modelled to be performed successfully in the control DLPFC group. NAD, NADP, adrenic acid and CMP-N-acetylneuraminate (CMP-Neu5Ac) de novo biosynthesis failed in both DLPFC groups. Adrenic acid biosynthesis also failed in other groups apart from the control TCX group; hence, its failure is more likely linked to aging rather than AD [39].   Regional differences in subsystem activities were remarkable. A drastic decrease in lacto-glycosphingolipid biosynthesis, lipoic acid metabolism, O-glycan metabolism, pentose and glucuronate interconversions, blood group biosynthesis and protein assembly were exclusively founded in the DLPFC group. Interestingly, O-glycan metabolism and pentose and glucuronate interconversions were increased in TCX but decreased in DLPFC. In addition, we observed a weak increase in propanoate metabolism, histidine metabolism and protein modification and a decrease in β-alanine metabolism and serotonin and melatonin biosynthesis in TCX. Different from other tissues, CBE samples showed higher activity of linoleate, retinol, xenobiotics and tryptophan metabolisms. In addition, serotonin and melatonin biosynthesis and propanoate metabolism were decreased in CBE but increased in TCX.

Figure 7.
Heatmap for metabolic pathways expressed at least 10% differently in one group of samples. Each colour tone closer to blood-red refers to 10% increase compared to other metabolic pathways, while each colour tone closer to deep blue 10% decrease compared to other metabolic pathways.
Then we performed the report metabolites analysis based on the result of differentially expressed analysis (Datasets S6-8). With the use of only down-regulated genes, we observed the significant changes in the three medium-chain fatty acids (pentanoyl-CoA, hexanoyl-CoA, (2E)-hexenoyl-CoA), acetate and H 2 O for all tissues (Figure 8). We observed that fatty acids were abundant in control groups for all three tissues, suggesting a general decrease in fatty acids including adrenic acid. The lower abundance of (R)-3hydroxybutanoate in DLPFC was noteworthy, since its importance for feeding neurons from astrocytes. Histidine was reported based on the up-regulated genes in TCX, which aligns with its lower metabolism. Retinoates were reported based on the down-regulated genes in DLPFC. Cholesterol, choline and maltose derivatives were reported based on the down-regulated genes in TCX. Vitamin D derivatives and metabolites responsible for the oxidation of leukotriene B4 (LTB4) were reported based on the down-regulated genes in CBE. LTB associated metabolites were reported based on the up-and down-regulated genes in TCX, implying an interplay. Similarly, proteoglycans were reported based on the up and down-regulated genes in DLPFC. Proteoglycan enrichments were opposing (enriched hsa05205_Proteoglycans_in_cancer and decreased O-glycan metabolism), suggesting that they may cover a different set of genes. CMP and CMP-Neu5Ac were reported based on the down-regulated genes in DLPFC, while CMP-Neu5Ac synthesis showed to be failed in the tissue.

Discussion
We performed network analysis to reveal the key genes and metabolites involved in the progression of AD. We found that shared DEGs (e.g., ATP1B3, RAF1 and STAT5B), hub genes (GPRASP2 and AMIGO1) and reporter metabolites (pentanoyl-CoA, hexanoyl-CoA, (2E)-hexenoyl-CoA and acetate) associated with AD progression are involved in signalling and energy metabolism. For instance, lower levels of unsaturated fatty acids were observed together with reduced activity of ATPases. Fatty acids enter into beta-oxidation, converted smaller acyl-CoA molecules, such as pentanoyl-CoA, in hepatocyte (liver) and astrocyte (brain). This continues down to the acetyl-CoA and 3-hydroxybutyrate which can diffuse into neurons. In neurons, they are not only used in energy production but also induce BDNF indirectly that promotes synaptic plasticity and stress resistance [40]. We also found that genes encoding ATPases, synaptic activity and plasticity, cytoskeleton organisation, vesicular transport (including endocytosis) and HSPs were tended to be down-regulated. Taken together with the depletion of unsaturated fatty acids and repression in energy metabolism may be affected by or cause synaptic loss and demyelination.
Myelination is becoming more central in our understanding of AD and aging, although studies have converged different opinions. One opinion is that AD only comes after the completion of myelination. Another claims that myelin damage precedes amyloid plaque deposition, becomes the first neuropathological abnormality in AD; even more, Aβ and tau proteins may be produced during myelin repair [41]. Brain sphingolipids are being considered important in AD due to the large conformational changes of Aβ in the presence of glycosphingolipids, synthesised from lactosylceramides, possibly ending in Aβ aggregation. In addition to them, membrane microdomains rich in cholesterol, ceramides and glycosphingolipids may control APP processing [42,43]. Sphingolipid metabolism and sphingolipid signalling pathway were found to be enriched in AD DLPFC. Ceramide, ceramide-1-phosphate, sphingosine, sphingosine-1-phosphate and D-gal-Nacylsphingosine were reported for this region for up-regulated genes (Figure 9). We can conclude an overall increase in these metabolites due to up-regulation of CERS2, KDSR and SPHK. The increase in sphingosine-1-phosphate and its receptors on the cellular membrane (S1PR1 and S1PR5) may induce the up-regulation of MAPK and PI3K/Akt signalling pathways. On the other hand, the down-regulation of B4GALT6 and B3GALNT1 may create a bottleneck for glycosphingolipid biosynthesis from lactosylceramides and nucleotide sugars, e.g., CMP-Neu5Ac (Figures 3a and 7 for corresponding enrichments). In the case of myelination, the down-regulation of the hub gene AMIGO1 is also critical. AMIGO1 expression is specific to the central nervous system and reduced in a pro-inflammatory environment, whereas its isomers are expressed more widely. It is promoted by neuronal growth-promoting factor amphoterin and has shown to be effective for dendritic growth. AMIGO1 was shown to be down-regulated in the post-mortem multiple sclerosis brain as well [44]. The unsolved question is thus how amyloid-β (Aβ) is related to brain degeneration. The modern perspective indicates the interplay between protein aggregation and other events that we provided evidence of. Aβ oligomers can trigger mitochondrial dysfunction and synaptic impairment, in three ways. First, they bind to the cellular form of prion protein on lipid rafts on post-synaptic neuron; this activates N-methyl-D-aspartate receptors (NMDARs), causing Ca 2+ overload. Ca 2+ overload triggers neurofibrillary tangle formation, mitochondrial dysfunction, leading to synaptic impairment and then apoptosis [45]. Second, oligomers interact with several NMDAR-inducing membrane receptors, such as EphB2 (Ephrin type-B receptor 2) and α7nAChr (acetylcholine receptor), leading to proteasomal degradation and NMDAR-induced synaptic impairment. Ca 2+ overload may occur in a third way by binding a specific form of Aβ to Na + /K + -ATPase α3 and its activation of voltage-gated Ca 2+ channels (VGCC) [46]. Supporting the expected low levels of Aβ due to the diminishment of APP processing genes, NMDAR genes (i.e., GRINA and GRIN2A in DLPFC), VGCC genes (i.e., CACNA2D1, CACNA1D, CACNA1E, CACNA1G), many ATPases, also DLG3 and CALM3, which mediate NMDARs and VGCCs respectively, were down-regulated. Possibly, HSPs were down-regulated in response to a decrease in protein aggregation. Considering these, the mechanisms of Ca 2+ overload in neurons should be lower.
Interestingly, we found no significant expressional changes for Aβ processing genes despite the tendency of lower expression of MAPT, APP and downstream APP processing genes in three brain regions other than underexpression of PSEN2 and overexpression of ADAM10 in DLPFC. While changes in DLPFC were more striking, there were counterintuitive points beginning with APP processing. The shared hub gene GPRASP2 plays an important role in APP processing as well. GPRASP2 encodes GPCR-associated sorting proteins (GPRASPs). These proteins interact with a broad spectrum of GPCRs and facilitate their degradation. Taken together, explaining the association between lower GPRASP2 expression and amyloid processing depends on further studies about the functionality and interactions of GPRASP2.
Another noteworthy aspect of our findings was related to cancer pathways. The antagonist relationship between cancer and degeneration diseases in aging including cardiovascular and neurodegenerative disorders is intriguing. The accumulating evidence from transcriptomics and genomics put forth an aging-driven trade-off between them; the loss of proteostasis, epigenetic changes and immune-response associated responses mostly occur in a proliferative way in tumour cells, whereas they occur in a way that shortens the lifespan of aged cells. Genomic instability, the secretion of proinflammatory cytokines and deregulated nutrient sensing, which interferes with insulin, mTOR and AMPK signalling pathways, are common to both paths [47][48][49]. In the case of AD, the tradeoff was observed at the population level. Retrospective and long-term prospective cohort studies have demonstrated clearly the lower probability of AD among cancer survivors and the lower probability of cancer in AD patients [49,50]. The up-regulation of many cancer-associated genes contributes to the overall enrichment, which plays a decisive role in cell fate. As given in Table 3, RAF1, SMAD4, RALBLP1 and STAT5B were up-regulated in both cortical regions and cerebellum. Firstly, RAF1 is a master regulator for cell fate decisions, such as proliferation, differentiation and apoptosis. It can activate the MAPK signalling pathway and may promote cell survival [51]. SMAD4 induced inhibition of TGF-β activity suppresses tumour activity [52]. In addition, TGF-β signalling pathway activation is important for synaptic transmission; therefore, high SMAD4 activity may be linked to the synaptic deficit as well [53]. RALBP1 is challenging to interpret; cause high RALPB1 activity is considered important for metastasis [54], while its depletion may lead to mitochondrial dysfunction and synaptic deficits [55]. Finally, STAT5B responds to cytokines and affects the transcription of transcription factors, such as FOXP3, CD25 and Bcl-2. These factors may induce cell development and prevent apoptosis, depending on the activated tissues [56]. In addition to cancer-associated genes, we find evidence for the switch towards degeneration in metabolic findings. Cancer cells are well known for their high energy and monomer consumption rates [57]. The diminishment in oxidative phosphorylation, which was reflected by lower ATPase expressions, and amino acid and fatty acid biosynthesis lowers the possibility of cancer co-occurrence, whereas the enrichment in signalling pathways explains the high enrichment in cancer pathways.
We assessed some metabolites from reporter metabolite analysis because of their neuroinflammatory or neuroprotective features. Adrenic acid is one of the most abundant fatty acids in the early human brain and is critical for inflammatory response together with arachidonic acid. Arachidonic acid is catalysed to hydroxy-peroxy-eicosatetraenoic acids (HPETEs) and hydroxy-eicosatetraenoic acids (HETEs), which are metabolised into leukotrienes. Leukotrienes have possible inflammatory effects; LTB4 may also enhance Aβ production by modulating γ-secretase activity [58]. On the other hand, it was shown adrenic acid inhibits LTB4 production by neutrophils and also enhances macrophage phagocytosis [59]. Here, we observed that both AD up-regulated and down-regulated genes contribute to the production of leukotrienes, HETEs and HPETEs in TCX and CBE. This may emphasise brain damage and neuroinflammation in the aging brain in general [60,61]. On the other hand, the poor levels of adrenic acid may show the ineffective response of the brain against AD-caused inflammation. Retinoids resemble adrenic acid in this regard. They are vitamin A derivatives that are present abundantly in the hippocampus, cortex and other brain regions. They inhibit LTB production by oxidation, thus, have pro-inflammatory effects. Additionally, they enhance choline acetyltransferase activity; therefore, they are beneficial for acetylcholine production and neurogenesis [60]. Reporter metabolite analysis (RMA) showed a decrease in retinoates in DLPFC, while retinol metabolism increased in CBE. Therefore, it is safer to assume that brain regions have different preferences for the use or production of these metabolites. Other vitamins and vitamin derivatives were mentioned for their neuroprotective effects in previous studies and interplay with fatty acids [62][63][64][65][66].
Propanoate and its products are another group of metabolites that have neuroprotective value. Excessive production of propionyl-CoA, due to high levels of propanoate in the gut, was shown to have neurotoxic effects such as mitochondrial dysfunction, neuroinflammation, glutamate excitotoxicity and decrease in serotonin levels. Propanoate metabolism is also linked to β-alanine metabolism and TCA cycle. Propionyl-CoA is metabolised to malonate semialdehyde, which is converted to acetyl-CoA and is interconverted to β-alanine by transaminase activity and involvement of ketone bodies and amino acids [67]. β-alanine is known to be neurotoxic [68]. It is also reported by both DLPFC and TCX up-regulated genes for this pathway. Propanoate metabolism was found to be more active in the AD TCX group compared to its control, which was represented by the conversion of 3-oxopropionyl-CoA to 3-hydroxypropionyl-CoA in TCX. The lower activity of serotonin biosynthesis in AD is also consistent. In contrast, the propanoate pathway was repressed and similarly, serotonin biosynthesis increased in the AD group of CBE. In addition, more metabolites for this pathway were reported for CBE, e.g., malonyl-CoA. Nevertheless, no DEGs were recorded for this pathway.
Another group of metabolites is bile acids. We observed AD selective failure of bilirubin conjugation in DLPFC. Bilirubin and biliverdin were also found in TCX upregulated genes. Bilirubin conjugation is the detoxification of degraded bilirubin with glucuronic acid, which is followed by its secretion. Bile acids facilitate dietary fat digestion and were linked to AD through their importance in cholesterol homeostasis and lipid absorption [68,69].
It seemed that many metabolites were severely affected by AD. Among them, unsaturated fatty acids (especially hexanoyl-CoA and pentanoyl-CoA), retinoids and other vitamins are known to improve cognitive abilities [68][69][70][71]. Their neuroprotective effects were studied extensively and put forth clinically [64][65][66]. Thus, reintroducing them may ameliorate cognitive decline. For instance, oral administration of a mixture including L-serine and L-carnitine tartrate has recently been shown to be effective at improving cognitive functions in AD patients [72]. In this regard, further studies are required to test the effects of metabolites that we found against cognitive decline.
In conclusion, we performed co-expression network analysis and metabolic modelling, and applied reporter metabolite analysis. Robustness of methods, analysis of a large number of samples and the use of different brain regions enhanced the reliability of our findings. In summary, we found a decline in energy production, myelination and synaptic activity in all brain regions. We also proposed some of the key genes as biomarkers. Among them, we found that AMIGO1 and GPRASP2 are highly critical for AD progression, and these genes should be investigated for their association to other pathological markers in large cohorts.

Data Sources and Sample Selection
Transcriptomic data (i.e., transcript abundances, transcript per million (TPM) values and metadata) of the Religious Orders Study and the Memory and Aging Project (ROSMAP) and Mayo Clinic RNAseq Study (MayoRNAseq) were obtained from the Accelerating Medicines Partnership-Alzheimer's Disease (AMP-AD) database (https://www.nia.nih. gov/research/amp-ad. Accessed 24 March 2021). Specimens in the ROSMAP study originated from DLPFC, while specimens in the MayoRNAseq study were from TCX and CBE, mostly from the same postmortem individuals. In the ROSMAP dataset, we excluded samples of MCI and other conditions contributing to cognitive impairment in order to remove potential confounders arising from other dementias; therefore, in total 423 samples (AD = 222, control = 201) were included (Supplementary Dataset S1). In MayoRNAseq dataset, TCX and CBE samples were obtained as matching samples from the same individuals. Similarly, we excluded samples of individuals with pathological aging or progressive supranuclear palsy. In total 150 samples (AD = 80, control = 70) for TCX and 144 samples (AD = 79, control = 65) for CBE were selected. For the sake of clarity in writing, ROSMAP dataset is abbreviated as DLPFC, and MayoRNAseq datasets are abbreviated as either TCX or CBE throughout the article.

Differential Expression Analysis and Gene Set Enrichment Analysis
The biomaRt package was used to access gene symbol and accession IDs from EN-SEMBL. Only transcripts of protein-coding genes were evaluated. Genes with low expressions (sum of TPM of a gene lower than 1) both in AD and control groups were discarded. The limma [73] removeBatchEffect function was used to remove batch effect by using flow cell number for DLPFC, source for TCX and both for CBE as covariates by inspecting top 1000 genes by variance on UMAP plots [74]. The effect of removing these covariates was visible on the plots (Figures S1a-f, S2a-I and S3-f). Differential expression analyses (DEAs) were performed using DESeq2 [75] "parametric" model, which was showed to be fit for mean of normalized gene counts ( Figure S4), and adding respective batch covariates into design matrix. Hypergeometric tests were performed to evaluate the concordance between each pair of gene lists [76]. GSEAs were performed by using GO and KEGG datasets. GO enrichments from UniProt [77] were checked using clusterProfiler package [78], enrichGO function. KEGG pathway enrichments were performed using fgsea package [79].

Co-Expression Network Analysis
Co-expression networks are suited well to reveal genes with similar regulation patterns across modules, which are highly connected gene sets, and ideal for making biological interpretations. Here, co-expression networks were built using combined AD and control data normalised for batch effects and based on Spearman correlations between genes (Dar S1). The genes involved in the top 1% of correlations by correlation coefficient were accepted for downstream analysis. Spearman correlation coefficients were converted to adjacency matrices and used to build undirected networks using the igraph package [80].
Using the random walk algorithm, densely connected subgraphs, henceforth called modules, were determined. Modules with a node size ≥ 20 and transitivity (clustering coefficient) > 0.5 were chosen for downstream analysis. Modules complying with these criteria from each network were matched to modules from other tissues and tested for significant overlaps by hypergeometric tests, which we used to indicate consistency between them and the presence of key genes within them. Then, differentially expressed gene sets were matched to these overlapping modules to understand whether co-expressed genes in these modules follow differential expression patterns, which was also tested by hypergeometric test. Next, we inspected the overlap of modules overlapping to gene set to see the coherence between tissues. We merged modules detected from each tissue. We narrowed genes of interest to hub genes (top 10% in terms of degree) assuming the higher probability of them being key genes or associated with them. GSEA was performed on them to obtain enriched biological annotations. Then, we looked at consistently shared hub genes among modules and searched for enriched metabolic pathways from HMR3 database.
We consulted the Interactome Atlas from the Human Reference Protein Interactome Mapping Project [81] to identify PPI partners of the key genes from co-expression analysis and their family members. For this purpose, we downloaded HuRI, HI-union, Test space screens-19 and Lit-BM tables, which have in total 96,370 interactions between 20,209 genes. We converted ENSEMBL IDs to gene symbols in reference to GENCODE v27. It must be noted that interactomes from these tables only show binary protein interactions, being blind to additional interaction partners and indirect associations in protein complexes. Even more, database is estimated to cover only a low proportion of human binary PPIs, hypothetically millions, and cellular function of them are mostly unknown; therefore, we evaluated PPIs as mostly complementary and informative rather than explanatory for results. Genes of interest were investigated in HPA to find out their normal protein and gene expression levels in brain tissues compared to other organs and their biological annotations.

Genome-Scale Metabolic Modelling and Reporter Metabolite Analysis
For each of DLPFC, TCX, and CBE datasets, we calculated AD and control consensus TPMs for each gene by taking arithmetic means across samples, without filtering out lowly expressed genes. Context-specific GEMs were constructed for each diagnosis group and for each brain region by overlaying group consensus TPMs onto HMR3 reference human model (8160 reactions and 3819 genes) and providing a list of 57 metabolic tasks that are compulsory for cellular growth. Modelling was performed on MATLAB v2020b (MathWorks Inc, USA), using tINIT algorithm in RAVEN 2.3.1 [82] toolbox for model construction and GUROBI (Gurobi Inc, TX) and MOSEK (Mosek ApS, Denmark) as the optimization solvers.
Constructed metabolic models were compared based on three aspects. First, the general reaction content similarities were compared, namely using Hamming distance. This shows the similarities of diagnosis groups on metabolite level. Second, the differences in subsystem (i.e., metabolic pathway) coverage were shown for subsystems with at least 10% difference in one or more GEMs to show pathways significantly changed between conditions or tissues. Third, the functionality of GEMs was compared to perceive the differences in their metabolic activity.
RMA shows the metabolites around which transcriptional changes occur significantly between different conditions by the use of DEA results (Datasets S6-S8). The method is sensitive to small perturbations and easily reveals patterns in metabolic network [83]. We applied the RAVEN implementation of RMA on DESeq2 results for each brain region, i.e., log2 fold changes in AD group compared to control group and p-values for these changes. We focused on significant metabolites with p-value ≤ 0.05 and evaluated metabolites reported from down-regulated and up-regulated genes separately.