Erosion of Gene Co-expression Networks Reveal Deregulation of Immune System Processes in Late-Onset Alzheimer’s Disease

We have applied a novel and integrative analysis framework for next-generation sequencing (NGS) data to 503 human subjects provided by the Religious Orders Study and Memory and Aging Project (ROSMAP) to examine changes in transcriptomic organization and common variants in association with late-onset Alzheimer’s disease (LOAD). Our framework identified seven reproducible, co-regulated modules after quality control (QC), clinical segregation, preservation filtering, and functional ontology analysis. These modules were specifically enriched in several innate and adaptive immune system processes, the synaptic vesicle cycle, and Hippo signaling. Topological and functional erosion of these modules due to shedding of genes and loss of in-module connectivity was diagnostic of disease progression. Perturbation analysis revealed that only 1% of eQTLs overlapped genes participating in these co-regulated modules. Common variants nevertheless identified components of the immune systems like human leukocyte antigen (HLA) complex and microtubule-associated protein tau (MAPT) regions in association with LOAD. Our results implicate microglial function, adaptive immune response, and the structural degeneration of neurons as contributors to the transcriptional deregulation observed along with common genetic variants in the progression of LOAD.

Introduction

Late-onset Alzheimer’s disease (LOAD) is a complex condition involving tau protein aggregates or tauopathy, amyloid and lipid processing, aging, immune system response, metabolism, lysosomal processing, and cerebrovascular health (Rogers et al., 1988; Braak et al., 2011; Jevtic et al., 2017; Wang et al., 2017). Progress in understanding and describing this large and diverse set of biological systems is in part determined by our ability to fully integrate clinical neuropathological data with comprehensive models that combine several modes of next-generation sequencing (NGS) data. To this end, we have applied our novel and integrative analysis framework to 503 subjects (305 cases/198 controls) provided by the Religious Orders Study and Memory and Aging Project (ROSMAP) study (Bennett et al., 2012) to develop a detailed landscape of the genetic and regulatory systems involved in LOAD, specifically with respect to clinical scores. Our study accomplished the following four objectives: (1) identified changes in transcriptional organization in association with clinical phenotypes; (2) characterized the systematic transcriptomic and functional changes accompanying LOAD through clinical segregation co-expression analysis; (3) identified common genomic loci involved in LOAD; and (4) tested the relationship between predicted expression quantitative trait loci (eQTLs) and systematic changes in gene expression.

One of the core elements of our approach is a weighted gene co-expression analysis (WGCNA), enabling the classification and identification of highly correlated and connected modules of genes grouped by co-expression (Zhang and Horvath, 2005; Langfelder and Horvath, 2008). Network modules can be described as series of interrelated nodes and edges. Here, nodes are messenger RNA (mRNA) transcripts. Edges represent the correlation coefficients between two or more given nodes, where degrees are the number of edges shared by nodes. Genes usually have many regulators, so we chose a hierarchical co-expression model. Our combined approach increases specificity by reducing large co-expression networks to only functionally significant and highly reproducible modules. Functional significance is defined as the Gene Ontology biological process p-value and reproducibility is defined as the module preservation Z-score (Langfelder et al., 2011). Co-expression analysis has been successfully applied in Alzheimer’s disease (AD), incorporating clinical scores and differential expression to identify co-regulated modules changing with disease in an “all-in-one” analytical design (Miller et al., 2008; Liang et al., 2018; Meng and Mei, 2019). However, the common approach to co-expression modeling does not include clinical segregation analysis. Here we provide clinical segregation for three groups: no cognitive impairment (NCI), mild cognitive impairment (MCI), and AD subjects. Additionally, strong genetic associations have been observed in LOAD (Naj et al., 2011; Lambert et al., 2013; Sims et al., 2017) along with systematic changes in gene expression profiles and transcriptional organization (Miller et al., 2008, 2010; Zhang et al., 2013; Ramasamy et al., 2014). Therefore, we hypothesize that genetic variation should account for changes in gene expression observed in transcriptomic analyses. We systematically tested the relationship between predicted eQTLs and transcriptomic organization to show underlying perturbations in gene networks that can partially account for changes observed in co-expression analyses.

Methods

The Accelerating Medicines Partnership (AMP) provides a variety of multi-platform next-generation sequence (NGS), clinical, and other –omics data. We selected all subjects from the ROSMAP study with overlapping clinical, RNA-seq, and DNA-seq data from the prefrontal cortex from a total of 503 elderly individuals varying from cognitively healthy to diagnosed AD (Bennett et al., 2012; De Jager et al., 2018). All subjects reported race as Caucasian. According to study details, RNA was extracted from the gray matter of the dorsolateral prefrontal cortex and quantified using the NanoDrop spectrophotometer.;101-bp paired-end, Illumina HiSeq reads were aligned to the human reference genome 19 (hg19). Genotype data were generated using the Affymetrix GeneChip 6.0 platform and filtered based on the following quality control (QC) criteria: genotype call rate less than 99%, minor allele frequency (MAF) less than 2%, and a Hardy–Weinberg equilibrium threshold below 1%. A total of 619,377 single nucleotide polymorphisms (SNPs) passed QC and were used in this analysis.

Analysis Framework

Our analytical framework was previously introduced (Malamon and Kriete, 2018) and extended here to include additional features, such as clinical segregation, module preservation, gene set, and functional enrichment analyses (see Supplementary Figure 1 for workflow diagram and full description of methods). This workflow consists of four main components: QC, co-expression modeling, functional enrichment, and eQTL analysis. We perform a comprehensive, three-tiered QC process to normalize and reduce the RNA-seq dataset to the 20,000 most informationally dense and connected transcripts. Co-expression networks are constructed using the WGCNA toolkit (Langfelder and Horvath, 2008). Next, we apply WGCNA’s module preservation testing procedure to measure statistical reproducibility in all modules. We exclude all modules with preservation Z-scores below 10 standard deviations. Higher Z-scores signify modules that reoccur despite changing input conditions. These become candidate modules. Functional term and enrichment analyses are performed on all candidate modules. Gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was used to examine larger functional network trends and reproduce candidate modules and genes. We provide a novel approach leveraging clinical segregation co-expression analysis to examine and compare alterations in network and module structure and organization with disease progression. For segregation analysis, QC, co-expression modeling, and functional enrichment were repeated for all clinical subgroups. Finally, genome-wide association (GWA) and perturbation analyses were performed. GWA provides all genomic loci (SNPs) predicted in association with disease status. eQTL analysis provides the predicted effects of SNPs on gene expression. Perturbation analysis was performed by overlapping co-expression module genes with eQTLs.

Clinical Segregation Analysis

Figure 1 outlines our clinical segregation protocol, which was designed to assess how transcriptomic differences are presented in clinical subgroups. We segregated samples by extracting sample data based on COGDX and CERAD scores and processing each group independently in WGCNA (see Supplementary Table 1 for clinical definitions and Supplementary Figure 2 for data plot). COGDX collapses 19 different neuropsychological tests into a single “Global Cognitive Score” (De Jager et al., 2018). The CERAD protocol provides neuropathological classifications for disease based on a wide variety of life-style, neuropsychological, and cognitive tests (Mirra et al., 1991). For COGDX, we segregated samples leaving 167 subjects with NCI, 131 subjects with MCI, and 205 subjects with an Alzheimer’s diagnosis (AD). For CERAD, we segregated samples leaving 130 subjects with no AD (CERAD_1), 226 subjects with possible or probable AD (CERAD_2), and 147 subjects with confirmed AD (CERAD_3). We independently processed and analyzed all six clinical subgroups using WGCNA with the same network parameters for all experiments.

FIGURE 1

Figure 1. Overview of clinical segregation co-expression analysis. An outline of our novel approach for independently analyzing and comparing co-expression networks and module characteristics in regard of clinical disease progression scores. The three vertical lanes represent COGDX segregation for three different cognitive scores (NCI, MCI, and AD) as defined in Supplementary Table 1. All subgroups were processed independently. First, quality control (QC) was applied to each set to retain only the 20,000 most informationally dense and variable transcripts. Next, networks were constructed with identical modeling parameters for all three subgroups. Module preservation (MP) testing was used to filter modules to only those that were highly reproducible (Z-score > 10), leaving seven modules. Within these modules, we observed a significant loss in the total transcripts classified, within-module connectivity, and functional term enrichment in association with disease progression. Heatmap tiles in the bottom lane refer to functionally significant GO biological process terms.

Results

Network Construction

We calculated the transcriptomic network’s total connectivity using the median-based bi-weight mid-correlation, which is more accurate than Pearson’s method for gene co-expression modeling (Zheng et al., 2014). Raising the soft-threshold to a power of 6 produced an overall R² value of 0.895, as seen in Supplementary Figure 3. Note that the R² value rises sharply and quickly flattens out with a slope of −1.080 at just six iterations. WGCNA was used to construct the initial co-expression network (Supplementary Figure 4) using all 503 subjects. Careful consideration was used in selecting the criteria for module identification, also known as branch trimming. We identified 26 distinct modules, totaling 4,429 transcripts with an average of 201 genes per module. Modules contain directionally signed groups of genes. In other words, genes in the same module are co-expressed in the same direction and well correlated with one another. We selected to lean on the side of specificity by not partitioning around medoids, leaving a total of 15,571 (77.85%) transcripts out of modules (unclassified), as indicated in gray (Supplementary Figure 4). Overall, the dendrogram shows clean, distinct clustering with sufficient levels of local dissimilarity. WGCNA arbitrarily assigns module names by color, i.e., gray and magenta. Supplementary Figure 5 shows all module-to-module and module-to-eigentrait (eigenvector of clinical metric) correlations for each of the four clinical NP traits. Supplementary Spreadsheet S1 contains all co-expressed genes grouped by module.

Clinical Segregation and Module Preservation Analysis

To investigate network and module characteristics with respect to disease progression, we segregated samples according to COGDX and CERAD scores and analyzed each of the six clinical subgroups independently in WGCNA. Clinical subgroups were assigned according to Supplementary Table 1. For example, subjects with COGDX scores of 0 or 1 were assigned to the NCI group. For COGDX, we uniquely classified 4,542, 3,998, and 2,992 genes for the NCI, MCI, and AD groups, respectively. For CERAD, we uniquely classified 3,426, 3,957, and 3,991 genes for the CERAD_1, CERAD_2, and CERAD_3 groups, respectively. WGCNA’s module preservation function allowed us to accurately measure module reproducibility through permutation testing. We calculated module preservation Z-statistics using 200 permutations for all six subgroups. See Supplementary Figure 6 for preservation statistics. Modules with Z-scores above 10 are not obtained by random chance and can be reliably reproduced (Langfelder et al., 2011; Li et al., 2015). A total of seven candidate modules (Table 1) survived preservation testing. Segregation by COGDX showed increased reproducibility and stability in module preservation over segregation based on CERAD assessment scores; therefore, we selected COGDX modules for further analysis.

TABLE 1

Table 1. Statistically significant functional terms for seven well-preserved modules sorted by adjusted p-value.

Functional Enrichment of Co-expression Modules

Biologically relevant, functional pathways should be reproducible and overlap known LOAD pathologies. To this end, we queried the GO database to examine the functional ontologies of the seven candidate modules. Table 1 provides statistically significant biological process terms involving known LOAD pathologies. The “magenta” module, which showed the strongest functional association, was well-correlated with COGDX and highly enriched with many immune-related genes including ABI3, APBB1IP, CD33, CD86, DOCK2, human leukocyte antigen (HLA)-DRA, HLA-DMB, MS4A4A, MS4A6A, MS4A7A, TREM2, and TYROBP. Other “magenta” GO terms include “complement pathway,” “cytokine signaling,” “neutrophil degranulation,” and “Toll-Like receptor activation.” Additional modules involving the immune system included the “turquoise” and “brown” modules, which were both enriched for the “regulation of complement activation.” A GO “cellular components” query revealed the “dendrite membrane” as significant for the “yellow” module (p-value = 3.73E-06). This observation is consistent with the “biological process” query results, which provided several synaptic processes including “neuronal projection,” “vesicle cycle,” and “synaptic maintenance.” The “blue” module was functionally enriched for genes in the Hippo signaling pathway, including AMOT, FAT4, LAT2, TJP1, TJP2, STK3, and YAP1. “Fatty acid oxidation” was also significant for the “blue” module. Additionally, cell-specific enrichment was performed on all seven modules (Uhlen et al., 2015; Kuleshov et al., 2016; Lachmann et al., 2018). See Supplementary Table 2 for results.

Organizational Changes in Immune Module

Figure 2 shows the erosion of the “magenta” module by comparing the network characteristics of the three co-expression networks segregated by COGDX. “Magenta” contained 191, 145, and 99 genes for NCI, MCI, and AD, respectively (Figure 2A) and shared 86 genes across all clinical groups. The mean intramodular degrees for the NCI, MCI, and AD subgroups provided in Figure 2B were 47.13, 39.27, and 31.27, respectively. The “blue” module shared 85 genes in all three subgroups with 859, 693, and 222 genes, respectively. The mean intramodular degrees for each subgroup were 409.32, 341.86, and 127.42 for the “blue” module. The “yellow” module shared 27 genes with 859, 693, and 222 genes, respectively. The mean intramodular degrees for “yellow” were 90.48, 74.69, and 12.66, respectively (p-value = 2.2E-16). Similar trends were observed in the other four modules. ANOVA and Bartlett’s test for heteroscedasticity were performed for all transcripts by COGDX subgroup revealing a significant (p-value < 0.05) increase in the expression of 22 “magenta,” 31 “yellow,” and 70 “blue” genes. Heteroscedasticity was significant (p-value < 0.001) for 36, 48, and 74 genes, respectively. These data provide supporting evidence for the deregulation of gene networks in these three modules. Supplementary Spreadsheet S2 contains all ANOVA and Bartlett’s testing results.

FIGURE 2

Figure 2. Erosion of nodes and edges for top three functional modules by CODGX segregation. (A) Venn diagram of genes in the immune-enriched module (magenta) for the three COGDX subgroups, NCI, MCI, and AD. (B) Boxplots with the number of intramodular connections (degrees) grouped by COGDX. (C) Venn diagram for “blue” module (Hippo Signaling). (D) Boxplots of degrees grouped by COGDX for “blue” module. (E) Venn diagram for “yellow” module (synaptic vesicle cycle). (F) Boxplots of degrees grouped by COGDX for the “yellow” module.

Gene Set Enrichment Analysis

The Broad Institute’s GSEA toolkit was used to identify disease-associated pathways via the KEGG biological pathway database (Kanehisa and Goto, 2000; Kanehisa et al., 2016, 2017). We performed a pre-ranked analysis with 10,000 permutations to discover differences in functional gene networks with regard to disease status. Supplementary Figure 7 provides the top and bottom five KEGG pathways sorted by p-value. The top five most significant (p-value < 1.5E-03) pathways positively enriched or under-represented in cases contained several immune-related genes also observed in co-expression modeling, including HLA-DRA, HLA-DMB, and CD86. Interestingly, cases exhibited deregulation in many immune system-related genes, which is consistent with the shedding of co-expressed genes revealed in the previous section. Negative enrichment scores denote an overrepresentation of pathway gene expression in cases. “Alzheimer’s, Parkinson’s, and Huntington’s disease” pathways showed high overrepresentation in cases.

Transcription Factor Analysis

Finally, we asked whether transcription factors may be influential for the observed changes in modules. Transcription factor binding site interrogation was performed using human single-site analysis (oPOSSUM) (Ho Sui et al., 2005) carried out at 8-bit minimum specificity, 40% conservation cutoff, 5,000 bp upstream/downstream the transcription start site, 85% matrix threshold, against a background of 24,752 genes. “Magenta” genes were highly enriched (p-value < 0.001) for the SPI1 and Interferon Regulatory Factor 8 (IRF8) transcription factor binding motif. Genes with SP1 binding site were also enriched in genes lost from “magenta” and included CD4, CYBA, HAMP, HCST, HLA-DMA, IL18, TLR10, and TREM2. Genes with PPARG:RXRA binding site included CD4, CYBA, HAMP, HCST, HLA-DMA, IL18, TLR10, and TREM2.

Expression Quantitative Trait Locus Analysis

We used MatrixEQTL (Shabalin, 2012) to test the linear associations between changes in gene expression and genotype for the same 503 individuals used in co-expression modeling. Interestingly, 90% of the top 100 eQTLs (sorted by adjusted p-value) occurred in the microtubule-associated protein tau (MAPT) region. Several HLA loci were statistically significant, including HLA-A, HLA-C, HLA-DOB, HLA-DP1, HLA-DRB1, and HLA-DRB5. Allele-specific changes in expression were observed not only on MAPT but also on MAPT-AS1, CRHR1, KANSL1-AS1, LRRC37A2, MAPK8IP1P1, and MAPK8IP1P2. Regional association plots for the MAPT and HLA-DPB2 regions were generated using LocusZoom (Pruim et al., 2010), provided in Figure 3. Linkage data were provided by the International HapMap Project (The International HapMap Consortium, 2003). Supplementary Figures 8–10 provide genome-wide association and box-plots of gene expression by genotype for four MAPT and four HLA-region SNPs identified in eQTL analysis. Supplementary Spreadsheet S3 contains all significant eQTLs with SNP (rsID), location, and p-value.

FIGURE 3

Figure 3. Genome-wide association plots for eQTL analysis. Regional association plots for MAPT and HLA-DPB2 regions. (A) Recombination rates (cM/Mb) (right vertical axis) and –log₁₀ (p-value) (left vertical axis) for SNPs with linkage peaks in blue for the MAPT and (B) HLA-DPB2 regions. SNPs are colored by the linkage disequilibrium correlation coefficient (R²). Genetic linkage data were provided by the International HapMap Project.

Perturbation Testing

To determine specific sources of genetic variation and their effects on the transcriptome, we overlapped all predicted eQTLs with all genes classified in co-expression modeling. Less than 1% of eQTLs (N = 5,392 gene/SNP pairs) across the 522 genes overlapped genes identified in co-expression network modules. We observed no discernable pattern in eQTLs and classified genes. Although observed differences in co-expression based on segregation are largely unexplained by individual eQTLs, functional ontology and transcription factor enrichment analysis provided polymorphisms in multiple genes sharing the transcription factor IRF motifs M08887 and M00972, which regulate many HLA genes.

Discussion

Our analysis revealed several key functional domains and pathways through which systematic deregulation occurs in LOAD. Co-expressed transcripts, transcription factors, and genomic loci were statistically significant contributors to LOAD progression via deregulation along several immune system pathways. In segregation co-expression analysis, we observed a substantial reduction in the organizational structure of several well-preserved, functional modules, as indicated by fewer classified genes and lower intramodular connectivity in MCI and AD subjects as compared to controls. In the course of this experiment, we improved specificity in detecting functionally relevant co-expression modules in a complex disease through rigorous QC protocols and data reduction schemes, namely, module preservation testing. This is significant because co-expression analyses produce very large networks with dozens of modules. This much data can be cumbersome and difficult to interpret. Network module erosion and gene shedding were observed in the microglia (“magenta”), synaptic (“yellow”), and Hippo signaling (“blue”) modules.

We chose to make the “magenta” module the focus of this discussion based on two statistical facts: (1) it was the most statistically significant module in functional gene enrichment analysis (Table 1) and (2) this module has been observed before in a similar study. Zhang et al. (2013) identified a module (“light cyan”) containing 537 genes in the human prefrontal cortex which were highly associated (p-value = 2.1e-87) with the same immune-related GO terms. Remarkably, 98 “light cyan” genes overlapped our “magenta” module. Assuming a hypergeometric distribution, the probability of identifying the same 98 genes from a total gene pool of 20,000 is 3.54e-129. We initially hypothesized that co-expression analysis would reveal cell-specific expression modules. Co-expression segregation analysis allowed us to compare specific changes in network and module organization.

The “magenta” module contained genes such as ABI3, CD33, MS4A46, MS4A6S, TREM2, and TYROBP, which have been previously linked to AD through protein-coding mutations (Naj et al., 2011; Sims et al., 2017) and are all critical to microglial activation and response (Satoh et al., 2017). Microglia are the principal innate immune cells of the brain and ingest and degrade amyloid plaques (Koenigsknecht-Talboo and Landreth, 2005). Segregation analysis based on COGDX showed that CD33 and TREM2 were co-expressed in NCI and MCI subjects but not in AD subjects, and CD4 was only co-expressed in the MCI module. CD33 and CD4 are associated with reactive microglia and have been linked to AD (Griciuc et al., 2013). TREM2 is activated by ligand binding and increases Aß clearance through the apoptosis-related phosphatidylinositol-3 kinase (PI3K) signaling pathway, while activating CD33 attenuates the innate immune response and Aß clearance. CD33 and TREM2 showed high heteroskedasticity in AD subjects and have been suggested as cross-talking Alzheimer’s genes (Chan et al., 2015). Taken together, our data suggest that CD33 and TREM2 co-regulation are important to maintaining healthy brain aging. ANOVA analysis of the “magenta” module showed increased expression in many genes including IL10RA and CD37. CD37 is activated by Aβ and mediates both humoral and cellular immune responses (van Spriel et al., 2004, 2009). This module also included HLA-DMA, HLA-DMB, and HLA-DRA.

Since the purpose of this study was to compare normal brain aging with AD, underlying aging pathways were not directly assessed. However, we noticed an interesting overlap with previous findings in a WGCNA study on the aging of the prefrontal cortex (Hu et al., 2018). Hu et al. reported a module enriched in the synaptic vesicle cycle function associated with brain aging progression. This module (“blue”) overlaps with the enrichment of our “yellow” module defined here. Within the GO term “synaptic vesicle cycle,” seven genes (AP2M1, ATP6V0D1, DNM1, RAB3A, STX1A, UNC13A, and VAMP2), involved in vesicle transport, endocytosis, and exocytosis, are shared between both studies. The difference in platforms and sample sizes makes this similarity remarkable, suggesting a further manifestation of synaptic dysfunction and impaired cognition in LOAD. The “blue” module was significantly enriched for Hippo signaling, which not only has implications on cell growth and autophagy but also the immune system (Zhang et al., 2018). Aß has been shown to initiate nuclear pro-apoptotic transcription factors in the Hippo signaling pathway, resulting in neuronal death (Sanphui and Biswas, 2013).

Our second motivation for the study was to examine the common genetic variants associated with LOAD. LOAD is likely influenced by the interaction of many polygenic, low- and moderate-effect variants. In our study, less than 1% of eQTLs overlapped genes classified in co-expression modeling. Of course, this did not directly explain changes in coexpression; however, eQTL analysis provided perturbations in multiple, functionally related genes (HLA-A, HLA-C, HLA-DOB, HLA-DP1, HLA-DRB1, and HLA-DRB5), all sharing transcription factor motif IRF. Interferon-regulatory factors modulate the interferon system in innate and adaptive immunity, and INF-γ induces differential expression of MHC class II HLA-DR and HLA-DP genes (Helbig et al., 1991). IFN-γ is expressed by infiltrating Th1 cells, resident microglia, and neurons and has been implicated in the development of AD and systemic autoimmunity. IFN-γ signaling is known to adversely affect AD pathologies and cognitive function (Mastrangelo et al., 2009; Monteiro et al., 2016). Activation of microglia by INF-γ inhibits Aβ clearance (Bate et al., 2006; Browne et al., 2013). HLA region eQTLs and changes in IFN-γ signaling can partially explain transcriptomic immune deregulation observed in cases.

The high concentration of eQTLs in the MAPT region highlights the impact of genetic variation on disease risk not only through MAPT haplotypes but also in several neighboring genes. MAPT pathologies provide a mechanistic link between the immune system and neurodegeneration involving microglia activation (Bhaskar et al., 2010). Splice-variants of MAPT-AS1 actively suppress MAPT translation (Coupland et al., 2016) and could prove to be a useful therapeutic target by reducing hyperphosphorylated tau levels. KANSL1 is critical to brain development (Koolen et al., 1993) and has been linked to AD (Jun et al., 2016). Corticotropin Releasing Hormone Receptor 1 (CRHR1) agonists have been shown to increase Aβ production (Futch et al., 2017). Determining the precise nature of the relationship between genetic variation and the expression of MAPT region genes will undoubtedly provide additional insights into tauopathy and thus LOAD risk.

Immune network architectures account for desirable immune system properties such as inducibility, adaptability, and robustness (Schrom et al., 2017). Data segregation combined with co-expression analysis sheds light onto these processes in LOAD, revealing adaptations during disease onset and erosion of networks in the later stages. Observed increases in transcriptional heterogeneity resemble observations in Parkinson’s disease (Mar et al., 2011), but can only partially account for module erosion since many highly variable genes are still present in the AD modules. Taken together, this study provides insights into a complex and dynamic landscape of genetic and regulatory processes centered around innate and adaptive immune system function. Systematic reductions in co-regulated genes and intramodular connectivity were diagnostic of increasing variability in several critical LOAD pathologies, including neuroinflammation, adaptive immunity, synaptic loss, and apoptosis. We propose that a reduction in regulatory and compensatory systems could also account for decreased robustness during disease progression, but the underlying mechanisms and the combined role of genetic variants are far from clear. This study highlights the adequacy of combining multi-omics NGS data types with longitudinal clinical and other developing, deep-phenotype data to decipher the complex molecular dynamics underlying complex diseases like LOAD.

Data Availability Statement

All data analyzed were obtained from the Accelerating Medicines Partnership for Alzheimer’s Disease (AMP-AD) Data Portal and can be accessed at https://www.synapse.org/#!Synapse:syn2580853/wiki/409844.

Author Contributions

JM performed co-expression, functional, GSEA, and eQTL analyses. AK performed TF analysis and supervised the study. JM and AK analyzed the functional outcomes and wrote the manuscript.

Funding

The study was supported in part by a Drexel Areas of Research Excellence (DARE) initiative on Cell2Society Aging Research (AK).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the many people who have contributed in myriad ways to the ROSMAP and Accelerating Medicines Partnership. We would also like to sincerely thank Drs. Ahmet Sacan and Catherine von Reyn (Drexel University’s School of Biomedical Engineering, Science and Health Systems), Dr. Daniel Marenda (Drexel University’s Department of Biology), and Dr. Peter Clark (University of Pennsylvania’s School of Medicine) for their advisement on this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2020.00228/full#supplementary-material

References

Bate, C., Kempster, S., Last, V., and Williams, A. (2006). Interferon-gamma increases neuronal death in response to amyloid-beta1-42. J. Neuroinflammation 3:7.

PubMed Abstract | Google Scholar

Bennett, D. A., Schneider, J. A., Arvanitakis, Z., and Wilson, R. S. (2012). Overview and findings from the religious orders study. Curr. Alzheimer Res. 9, 628–645. doi: 10.2174/156720512801322573

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhaskar, K., Konerth, M., Kokiko-Cochran, O. N., Cardona, A., Ransohoff, R. M., and Lamb, B. T. (2010). Regulation of tau pathology by the microglial fractalkine receptor. Neuron 68, 19–31. doi: 10.1016/j.neuron.2010.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Braak, H., Thal, D. R., Ghebremedhin, E., and Del Tredici, K. (2011). Stages of the pathologic process in Alzheimer disease: age categories from 1 to 100 years. J. Neuropathol. Exp. Neurol. 70, 960–969. doi: 10.1097/NEN.0b013e318232a379

PubMed Abstract | CrossRef Full Text | Google Scholar

Browne, T. C., McQuillan, K., McManus, R. M., O’Reilly, J. A., Mills, K. H., and Lynch, M. A. (2013). IFN-gamma Production by amyloid beta-specific Th1 cells promotes microglial activation and increases plaque burden in a mouse model of Alzheimer’s disease. J. Immunol. 190, 2241–2251. doi: 10.4049/jimmunol.1200947

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, G., White, C. C., Winn, P. A., Cimpean, M., Replogle, J. M., Glick, L. R., et al. (2015). CD33 modulates TREM2: convergence of Alzheimer loci. Nat. Neurosci. 18, 1556–1558. doi: 10.1038/nn.4126

PubMed Abstract | CrossRef Full Text | Google Scholar

Coupland, K. G., Kim, W. S., Halliday, G. M., Hallupp, M., Dobson-Stone, C., and Kwok, J. B. (2016). Role of the long non-coding RNA MAPT-AS1 in regulation of microtubule associated protein Tau (MAPT) expression in Parkinson’s disease. PLoS One 11:e0157924. doi: 10.1371/journal.pone.0157924

PubMed Abstract | CrossRef Full Text | Google Scholar

De Jager, P. L., Ma, Y., McCabe, C., Xu, J., Vardarajan, B. N., Felsky, D., et al. (2018). A multi-omic atlas of the human frontal cortex for aging and Alzheimer’s disease research. Sci. Data 5:180142. doi: 10.1038/sdata.2018.142

PubMed Abstract | CrossRef Full Text | Google Scholar

Futch, H. S., Croft, C. L., Truong, V. Q., Krause, E. G., and Golde, T. E. (2017). Targeting psychologic stress signaling pathways in Alzheimer’s disease. Mol. Neurodegener. 12:49. doi: 10.1186/s13024-017-0190-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Griciuc, A., Serrano-Pozo, A., Parrado, A. R., Lesinski, A. N., Asselin, C. N., Mullin, K., et al. (2013). Alzheimer’s disease risk gene CD33 inhibits microglial uptake of amyloid beta. Neuron 78, 631–643. doi: 10.1016/j.neuron.2013.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Helbig, H., Kittredge, K. L., Palestine, A. G., Coca-Prados, M., and Nussenblatt, R. B. (1991). Gamma-interferon induces differential expression of HLA-DR, -DP and -DQ in human ciliary epithelial cells. Graefes Arch. Clin. Exp. Ophthalmol. 229, 191–194. doi: 10.1007/bf00170556

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho Sui, S. J., Mortimer, J. R., Arenillas, D. J., Brumm, J., Walsh, C. J., Kennedy, B. P., et al. (2005). oPOSSUM: identification of over-represented transcription factor binding sites in co-expressed genes. Nucleic Acids Res. 33, 3154–3164. doi: 10.1093/nar/gki624

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Y., Zhao, T., Zang, T., Zhang, Y., and Cheng, L. (2018). Identification of Alzheimer’s disease-related genes based on data integration method. Front. Genet. 9:703. doi: 10.3389/fgene.2018.00703

PubMed Abstract | CrossRef Full Text | Google Scholar

Jevtic, S., Sengar, A. S., Salter, M. W., and McLaurin, J. (2017). The role of the immune system in Alzheimer disease: etiology and treatment. Ageing Res. Rev. 40, 84–94. doi: 10.1016/j.arr.2017.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jun, G., Ibrahim-Verbaas, C. A., Vronskaya, M., Lambert, J. C., Chung, J., Naj, A. C., et al. (2016). A novel Alzheimer disease locus located near the gene encoding Tau protein. Mol. Psychiatry 21, 108–117. doi: 10.1038/mp.2015.23

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361. doi: 10.1093/nar/gkw1092

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., and Goto, S. (2000). KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M., and Tanabe, M. (2016). KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462. doi: 10.1093/nar/gkv1070

PubMed Abstract | CrossRef Full Text | Google Scholar

Koenigsknecht-Talboo, J., and Landreth, G. E. (2005). Microglial phagocytosis induced by fibrillar beta-amyloid and IgGs are differentially regulated by proinflammatory cytokines. J. Neurosci. 25, 8240–8249. doi: 10.1523/JNEUROSCI.1808-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Koolen, D. A., Morgan, A., and de Vries, B. B. A. (1993). “Koolen-de vries syndrome,” in GeneReviews((R)), eds M. P. Adam, H. H. Ardinger, R. A. Pagon, S. E. Wallace, L. J. H. Bean, K. Stephens, et al. (Seattle, WA: OMIM).

Google Scholar

Kuleshov, M. V., Jones, M. R., Rouillard, A. D., Fernandez, N. F., Duan, Q., Wang, Z., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachmann, A., Torre, D., Keenan, A. B., Jagodnik, K. M., Lee, H. J., Wang, L., et al. (2018). Massive mining of publicly available RNA-seq data from human and mouse. Nat. Commun. 9:1366. doi: 10.1038/s41467-018-03751-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, J. C., Ibrahim-Verbaas, C. A., Harold, D., Naj, A. C., Sims, R., Bellenguez, C., et al. (2013). Meta-analysis of 74,046 individuals identifies 11 new susceptibility loci for Alzheimer’s disease. Nat. Genet. 45, 1452–1458. doi: 10.1038/ng.2802

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is my network module preserved and reproducible? PLoS Comput. Biol. 7:e1001057. doi: 10.1371/journal.pcbi.1001057

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., Zhang, Y., Yu, Y., Wang, P., Wang, Y., Wang, Z., et al. (2015). Quantitative assessment of gene expression network module-validation methods. Sci. Rep. 5:15258. doi: 10.1038/srep15258

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, J. W., Fang, Z. Y., Huang, Y., Liuyang, Z. Y., Zhang, X. L., Wang, J. L., et al. (2018). Application of weighted gene co-expression network analysis to explore the key genes in Alzheimer’s disease. J. Alzheimers Dis. 65, 1353–1364. doi: 10.3233/JAD-180400

PubMed Abstract | CrossRef Full Text | Google Scholar

Malamon, J. S., and Kriete, A. (2018). Integrated systems approach reveals sphingolipid metabolism pathway dysregulation in association with late-onset Alzheimer’s Disease. Biology 7:16. doi: 10.3390/biology7010016

PubMed Abstract | CrossRef Full Text | Google Scholar

Mar, J. C., Matigian, N. A., Mackay-Sim, A., Mellick, G. D., Sue, C. M., Silburn, P. A., et al. (2011). Variance of gene expression identifies altered network constraints in neurological disease. PLoS Genet. 7:e1002207. doi: 10.1371/journal.pgen.1002207

PubMed Abstract | CrossRef Full Text | Google Scholar

Mastrangelo, M. A., Sudol, K. L., Narrow, W. C., and Bowers, W. J. (2009). Interferon-{gamma} differentially affects Alzheimer’s disease pathologies and induces neurogenesis in triple transgenic-AD mice. Am. J. Pathol. 175, 2076–2088. doi: 10.2353/ajpath.2009.090059

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, G., and Mei, H. (2019). Transcriptional dysregulation study reveals a core network involving the progression of Alzheimer’s disease. Front. Aging Neurosci. 11:101. doi: 10.3389/fnagi.2019.00101

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, J. A., Horvath, S., and Geschwind, D. H. (2010). Divergence of human and mouse brain transcriptome highlights Alzheimer disease pathways. Proc. Natl. Acad. Sci. U.S.A. 107, 12698–12703. doi: 10.1073/pnas.0914257107

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, J. A., Oldham, M. C., and Geschwind, D. H. (2008). A systems level analysis of transcriptional changes in Alzheimer’s disease and normal aging. J. Neurosci. 28, 1410–1420. doi: 10.1523/JNEUROSCI.4098-07.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirra, S. S., Heyman, A., McKeel, D., Sumi, S. M., Crain, B. J., Brownlee, L. M., et al. (1991). The consortium to establish a registry for Alzheimer’s disease (CERAD). Part II. Standardization of the neuropathologic assessment of Alzheimer’s disease. Neurology 41, 479–486. doi: 10.1212/wnl.41.4.479

PubMed Abstract | CrossRef Full Text | Google Scholar

Monteiro, S., Ferreira, F. M., Pinto, V., Roque, S., Morais, M., de Sa-Calcada, D., et al. (2016). Absence of IFN-gamma promotes hippocampal plasticity and enhances cognitive performance. Transl. Psychiatry 6:e707. doi: 10.1038/tp.2015.194

PubMed Abstract | CrossRef Full Text | Google Scholar

Naj, A. C., Jun, G., Beecham, G. W., Wang, L. S., Vardarajan, B. N., Buros, J., et al. (2011). Common variants at MS4A4/MS4A6E, CD2AP, CD33 and EPHA1 are associated with late-onset Alzheimer’s disease. Nat. Genet. 43, 436–441. doi: 10.1038/ng.801

PubMed Abstract | CrossRef Full Text | Google Scholar

Pruim, R. J., Welch, R. P., Sanna, S., Teslovich, T. M., Chines, P. S., Gliedt, T. P., et al. (2010). LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26, 2336–2337. doi: 10.1093/bioinformatics/btq419

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramasamy, A., Trabzuni, D., Guelfi, S., Varghese, V., Smith, C., Walker, R., et al. (2014). Genetic variability in the regulation of gene expression in ten regions of the human brain. Nat. Neurosci. 17, 1418–1428. doi: 10.1038/nn.3801

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers, J., Luber-Narod, J., Styren, S. D., and Civin, W. H. (1988). Expression of immune system-associated antigens by cells of the human central nervous system: relationship to the pathology of Alzheimer’s disease. Neurobiol. Aging 9, 339–349. doi: 10.1016/s0197-4580(88)80079-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanphui, P., and Biswas, S. C. (2013). FoxO3a is activated and executes neuron death via Bim in response to beta-amyloid. Cell Death Dis. 4:e625. doi: 10.1038/cddis.2013.148

PubMed Abstract | CrossRef Full Text | Google Scholar

Satoh, J. I., Kino, Y., Yanaizu, M., Tosaki, Y., Sakai, K., Ishida, T., et al. (2017). Microglia express ABI3 in the brains of Alzheimer’s disease and Nasu-Hakola disease. Intractable Rare Dis. Res. 6, 262–268. doi: 10.5582/irdr.2017.01073

PubMed Abstract | CrossRef Full Text | Google Scholar

Schrom, E., Prada, J. M., and Graham, A. L. (2017). Immune signaling networks: sources of robustness and constrained evolvability during coevolution. Mol. Biol. Evol. 35, 676–687. doi: 10.1093/molbev/msx321

PubMed Abstract | CrossRef Full Text | Google Scholar

Shabalin, A. A. (2012). Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28, 1353–1358. doi: 10.1093/bioinformatics/bts163

PubMed Abstract | CrossRef Full Text | Google Scholar

Sims, R., van der Lee, S. J., Naj, A. C., Bellenguez, C., Badarinarayan, N., Jakobsdottir, J., et al. (2017). Rare coding variants in PLCG2, ABI3, and TREM2 implicate microglial-mediated innate immunity in Alzheimer’s disease. Nat. Genet. 49, 1373–1384. doi: 10.1038/ng.3916

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

The International HapMap Consortium, (2003). The International HapMap Project. Nature 426, 789–796. doi: 10.1038/nature02168

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlen, M., Fagerberg, L., Hallstrom, B. M., Lindskog, C., Oksvold, P., Mardinoglu, A., et al. (2015). Proteomics. Tissue-based map of the human proteome. Science 347:1260419. doi: 10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

van Spriel, A. B., Puls, K. L., Sofi, M., Pouniotis, D., Hochrein, H., Orinska, Z., et al. (2004). A regulatory role for CD37 in T cell proliferation. J. Immunol. 172, 2953–2961. doi: 10.4049/jimmunol.172.5.2953

PubMed Abstract | CrossRef Full Text | Google Scholar

van Spriel, A. B., Sofi, M., Gartlan, K. H., van der Schaaf, A., Verschueren, I., Torensma, R., et al. (2009). The tetraspanin protein CD37 regulates IgA responses and anti-fungal immunity. PLoS Pathog. 5:e1000338. doi: 10.1371/journal.ppat.1000338

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Gu, B. J., Masters, C. L., and Wang, Y. J. (2017). A systemic view of Alzheimer disease - insights from amyloid-beta metabolism beyond the brain. Nat. Rev. Neurol. 13:703. doi: 10.1038/nrneurol.2017.147

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Gaiteri, C., Bodea, L. G., Wang, Z., McElwee, J., Podtelezhnikov, A. A., et al. (2013). Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer’s disease. Cell 153, 707–720. doi: 10.1016/j.cell.2013.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., and Horvath, S. (2005). A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 4:Article17. doi: 10.2202/1544-6115.1128

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zhang, H., and Zhao, B. (2018). Hippo signaling in the immune system. Trends Biochem. Sci. 43, 77–80. doi: 10.1016/j.tibs.2017.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, C. H., Yuan, L., Sha, W., and Sun, Z. L. (2014). Gene differential coexpression analysis based on biweight correlation and maximum clique. BMC Bioinformatics 15(Suppl. 15):S3. doi: 10.1186/1471-2105-15-S15-S3

PubMed Abstract | CrossRef Full Text | Google Scholar