Immune gene network of neurological diseases: Multiple sclerosis (MS), Alzheimer's disease (AD), Parkinson's disease (PD) and Huntington's disease (HD).

Neurological diseases, such as MS, AD, PD and HD, are a major health concern of the elderly population, but still therapeutic options are limited. Recent advances in genomic sequencing and bioinformatics, present an opportunity to understand mechanisms of these diseases for identification of therapeutic targets. Several studies have shown association of immune dysfunction with immune system mediated neurological disease, MS, as well as neurodegenerative diseases (AD, PD and HD). However, similarities and differences in role of the immune system, immune pathways and immune cell types in these diseases remains unknown. In this study, immune cell type signature genes in gene networks associated with neurological diseases, MS, AD, PD and HD was investigated using meta-analysis and bioinformatics methods. Application of Weighted Gene Co-expression Network Analysis (WGCNA) on publicly available gene expression datasets (microarray and RNA-seq) revealed a ModArray_04 module (microarray) or ModRNAseq_06 module (RNA-seq), significantly associated with MS, AD, PD and HD. Hypergeometric enrichment test revealed significant enrichment of immune cell type genes in these neurological disease modules. This study demonstrates that immune system mediated neurological disease, MS and neurodegenerative diseases (AD, PD and HD), share a common gene network characterized by immune cell type signature genes (microglia, monocytes and macrophages) and are probable targets for therapeutic intervention. In summary, this work shows a connection between MS, a disease where the role of the immune system and inflammation is established, and neurodegenerative diseases (AD, PD and HD) where the role of inflammation is still a hypothesis.

Introduction

It is estimated that number of adults above 65 years of age will outnumber children in U.S. by 2035. This presents a significant challenge to healthcare because age increases susceptibility to life changing neurodegenerative diseases and cancers [1]. Alzheimer's disease (AD) is characterized by appearance of extracellular amyloid-beta plaques, intracellular neurofibrillary tangles and neuronal cell death particularly in cortex, and hippocampus [2]. AD is most common form of dementia occurring in elderly people and overall prevalence of AD is expected to double in next 20 years [3]. Aging is risk factor for Parkinson's disease (PD) and its hallmarks are neurodegeneration in substantia nigra and sever motor dysfunction [4]. PD effects 1% of U.S. population above age 60 and increases to 5% in individuals above age 85 [4]. Huntington's disease (HD) is an inherited neurodegenerative disease caused by unstable CAG trinucleotide repeat expansion in huntingtin (HTT) gene and manifests as cognitive, psychiatric and/or movement dysfunction [5, 6]. Adult onset HD is most common form of this disease and can occur between age 2 to 85 [7]. Therefore, aging is a major risk factor for neurodegenerative diseases, AD, PD and HD.

Several studies have sought to identify mechanistic links between aging and neurodegenerative diseases. The goal of this introduction is to provide readers with an overview of all known mechanisms of AD, PD, HD and MS. Mitochondrial dysfunction, immune system and inflammation are emerging as crucial mechanistic links between aging and neurodegeneration. Mitochondria are cellular organelles that produce cellular energy in form of ATP by oxidative phosphorylation of fats, and carbohydrates in mitochondria [8]. Reduced size and number of mitochondria, altered mitochondrial oxidative phosphorylation, and altered mitochondria morphology occur in AD, PD and HD [9, 10, 11]. ‘Neuroinflammation hypothesis’ suggests dysfunction of immune system response and inflammation in central nervous system as a driver of neurodegenerative diseases [12]. Microglia, the resident immune cells of brain, are widely studied in context of neurodegeneration and aging. Several studies have shown that microglia signature genes are altered in three neurodegenerative diseases, AD, PD and HD, in context of aging [13, 14, 15, 16, 17, 18]. Due to increased leakiness of blood brain barrier and systemic alterations in immune system with aging, immune cells from peripheral blood enter the brain and influence brain inflammation [19, 20]. However, the role played by different immune cell types in neurodegenerative diseases needs further investigation. Especially, it will be interesting to tease out these regulatory immune mechanisms in context of neurodegeneration independent of aging.

Multiple Sclerosis (MS) is most common chronic inflammatory brain disease mediated by immune system. MS is characterized by demyelination, and axonal degeneration [21]. MS typically occurs in adults between 20 to 40 years of age and is diagnosed based on demyelination lesions [22]. Unlike neurodegenerative diseases, such as AD, PD and HD, old age does not necessarily increase risk of MS. Though no driver gene(s) is known to cause MS, studies have shown that genetic predisposition and environmental factors increase risk of MS [23]. Consistent with immune origin of MS, early life infection by Epstein-Barr virus (EBV) or human herpesvirus 4 and genetic variants of immune genes increase risk of MS [24, 25]. Chronic activation of immune pathway in neurodegenerative diseases (AD, PD and HD) is well known. However, these studies were done on one disease at a time, which makes it challenging to compare the role of immune system across all neurodegenerative diseases (AD, PD and HD) [26, 27, 28, 29]. Two studies Durrenberger et al and Mukherjee et al, compared multiple neurodegenerative diseases and found a neurodegenerative disease associated immune pathway, and immune cell type (microglia), respectively [13, 30]. However, the role played by other immune cells, such as innate immune cell types, monocytes and macrophages, in neurodegenerative diseases (AD, PD and HD) and how it compares to classic immune system mediated disease, MS, remains unknown.

In this work, a side-by-side comparison was made between neurodegenerative diseases (AD, PD and HD) and immune system disorder, MS. Bioinformatics analysis revealed that three immune cell types, microglia, monocytes and macrophages, underlie neurodegenerative diseases (AD, PD and HD) and immune dysfunction mediated neurological disease MS, independent of age.

Methods

Code availability

Computational pipeline used in this work is available on github https://github.com/smukher2/GithubSubmittedADPDHDMSimmuneSept2021.

Source of RNA-seq and microarray samples

Publicly available datasets on NIH GEO https://www.ncbi.nlm.nih.gov/geo/ were used in this study. Human microarray datasets were obtained from cortex (control and AD), substantia nigra (control and PD) and whole brain (control and MS). Human RNA-seq datasets were obtained from prefrontal cortex (control, AD, HD, PD), optic chiasm (control and MS) and hippocampus (control and AD). Mouse RNA-seq datasets were used to identify brain cell type genes (astrocytes, neurons, oligodendrocytes, microglia) and blood immune cell type genes (monocytes and macrophages). The series numbers for these datasets are provided at the end of this manuscripts under acknowledgements.

RNA-seq raw reads processing

RNA-seq analysis pipeline was adapted and developed from previously published standard methodologies [13, 31]. Briefly, RNA-seq raw reads were aligned to GRCh38 and GRCm38 genomes with TopHat using gene annotation information from Homo_sapiens.GRCh38.92.gtf and Mus_musculus.GRCm38.92.gtf files, for human datasets (AD, PD, HD and MS), and mouse datasets (brain cell types and immune cell types), respectively, on cyverse high-performance computational cluster [32]. Mapped output bam files were used as input for HTSeq-Count to count reads overlapping genes for estimation of gene abundance [33]. All further analyses, including graphical visualizations were done in R environment with R packages and custom R functions. Raw gene counts were converted to TPM values (filtered to keep values > 1 in at least 2 samples) to normalized for variations in library size and mRNA abundance. TPM values were further scaled to log2TPM+1 scale and used as input for rest of the bioinformatics pipeline. The log2TPM+1 values were visualized with volcano plots, barplots and density plots using ggplot2_3.1.1 R package [34]. Built-in R functions, prcomp stats4_3.5.0 and corrplot_0.84, were used to perform PCA analyses and correlation analysis, respectively, to visualize variability of replicate samples before and after TPM normalization [34, 35].

Microarray data processing

Previously published microarray analysis pipeline was used here [13]. GEOquery R package was used to obtain microarray datasets and GPL annotation files [36]. Expression values from all datasets were quantile normalized using limma R package and converted to log2+1 scale [37]. Probe IDs were converted to gene symbols and in case of duplicates, the first probe ID gene expression value was kept for all datasets. The log2+1 expression values were used as input for estimation of DEGs (Differentially Expressed Genes) and construction of gene co-expression networks (WGCNA network). The log2+1 expression values were visualized with volcano plots, barplots and density plots using ggplot2_3.1.1 R package [34].

Consensus differential gene expression detection from limma, edgeR and simple comparison of means

Differentially expressed genes (DEGs) between a given cell type and all other cell types, were identified using limma_3.38.3, edgeR_3.24.3 and simple comparison expression means [37, 38, 39]. In limma and edgeR, model design included condition, study or batch, gender, age and tissue. For calculation of DEGs in limma and edgeR, contrasts were made for the two conditions being compared holding all other variables in design model constant. Limma and edgeR methods included correction for false discovery and multiple testing using Benjamini and Hochberf (BH) to correct for the presence of >30K transcripts (genes, pseudogenes and antisense) in the analysis [40]. Though DEG analysis outputs both upregulated and downregulated genes, here only upregulated genes were selected as cell type specific signature genes. Signature or marker genes were defined by significant presence of genes (upregulation) and not significant absence of genes (downregulation) because generally cells are identified by presence of marker genes for techniques such as, Fluorescence-activated cell sorting and immunohistochemistry. BH corrected adjP-values <0.05 and a fold changes >5 were cut-offs used for significantly upregulated genes. Gene expression of DEGs were visualized with volcano plots, barplots and density plots using ggplot2_3.1.1 R package [34]. Consensus DEGs were those genes that were present in DEG lists of at least two of the three methods, limma, edgeR and simple comparison of means. Significance of overlap between DEG lists were calculated using GeneOverlap_1.18.0 R package and number of genes overlapping were visualized as Venn-diagrams using VennDiagram_1.6.20 R package [41, 42]. UserListEnrichment R function was used to uncover any significant overlap of cell type genes between cell types [29].

Consensus gene co-expression network detection from RNA-seq and microarray datasets with WGCNA

To build gene co-expression networks or modules of highly co-expressed genes significantly associated with traits of interest WGCNA_1.66 R package was used [43, 44]. To model and remove gene expression effects of covariants, such as study or batch, gender, age, tissue and surrogate variables (hidden variables), and retain only disease condition effects, surrogate variable adjustment and linear model (SVA + LM) methods were applied with sva_3.30.1 and limma_3.38.3 R packages [13, 37, 45, 46, 47]. SVA + LM adjustment and WGCNA analysis were done in parallel on both RNA-seq and microarray datasets. SVA + LM adjusted values were converted to log2+1 normalized gene expression and used as input for WGCNA to construct gene networks or modules based on gene co-expression. In WGCNA, gene expression correlations were determined to calculate a topology overlap matrix (TOM) and hierarchical clustering method was applied to build weighted gene co-expression networks or modules with default parameters and minimum module size 100 [43,44]. WGCNA moduleEigengenes R function was used to represent each module with its own eigengene [43, 44]. While Pearson correlation is good for linear relations between continuous variables, Spearman correlation also works for monotonic relations and ordinal (categorical) variables. As biological interactions are complex, often monotonic and the disease state trait is categorical, here Spearman correlation was used to determine module disease association. To determine module trait association, Spearman correlation was calculated for pairs of module eigengenes and traits. Module eigengenes are computed by the WGCNA package and represent the gene expression pattern for all genes in a given module. Disease covariant or trait was split into pairs of AD-control, PD-control, HD-control and MS-control, to identify modules associated with a given disease independent of other diseases. Consensus modules between RNA-seq and microarray datasets were identified using module preservation and useListEnrichment R functions from WGCNA [29, 43, 44]. Preservation analysis was used to determine if connectivity of genes in microarray modules held true in RNA-seq and significance was reported as a Z-summary score. Z-summary scores higher than 5 are recommended to be significant as per creators of WGCNA package. Only disease associated modules that were preserved between the two platforms (microarray and RNA-seq) (Figure 3E) and had significantly overlapping genes (Figure 3F) were selected. We used this method of selection as irrespective of platform the modules showed disease association. The gene overlap also helps identify the name of the equivalent module in RNA-seq platform given that preservation analysis only takes the input of the reference microarray network module names.

Figure 3

Identification of AD, PD, HD and MS associated WGCNA modules in microarray and RNA-seq datasets. A: Bar-plot showing number of genes in AD, PD, HD and MS associated microarray WGCNA modules. B: Spearman correlation coefficients and p-values for significant AD, PD, HD, and/or MS associated microarray WGCNA modules. C: Bar-plot showing number of genes in AD, PD, HD and MS associated RNA-seq WGCNA modules. D: Spearman correlation coefficients and p-values for significant AD, PD, HD and MS associated RNA-seq WGCNA modules. E: Preservation analysis between all microarray and RNA-seq WGCNA modules using microarray WGCNA module labels. F: Table showing microarray WGCNA modules associated with AD, PD, HD and MS that significantly overlapped with RNA-seq WGCNA modules. For this figure the following R packages were used: WGCNA_1.66 and ggplot2_3.1.1 [34,43,44].

Characterization of gene co-expression disease associated modules

UserListEnrichment R function was used to determine enrichment of cell type signature genes in WGCNA modules [29]. EnrichR_1.0 a R package was used to determine biological process gene ontology (GO) characterization of the modules [48]. Cell type signature genes underlying only significantly preserved disease associated modules, as determined by preservation analysis in the previous step, were reported in Table 2.

Table 2

Significant cell types enriched in microarray and RNA-seq WGCNA modules associated with AD, PD, HD and MS.

WGCNA_Array_ module_names	DiseaseAD	DiseaseHD	DiseasePD	DiseaseMS	p-value_ Cell_Type_ Enrichment_ in_Array_ WGCNA_ Modules	Cell_Type_ Marker_Genes	WGCNA_RNA-seq_ module_names	p-value_ Cell_Type_ Enrichment_ in_RNA-seq_ WGCNA_ Modules
	p-value	p-value	p-value	p-value
ModArray_03	9.47E-09	2.89E-13	1.86E-102	3.16E-08	3.26E-16	neuronBYothers	ModRNAseq_14	3.78E-80
ModArray_04	1.92E-36	6.27E-47	0.0005408	4.04E-103	7.77E-08	macrophagesBYothers	ModRNAseq_06	0.00389653
ModArray_04	1.92E-36	6.27E-47	0.0005408	4.04E-103	2.91E-21	microgliaBYothers	ModRNAseq_06	8.32E-39
ModArray_04	1.92E-36	6.27E-47	0.0005408	4.04E-103	7.38E-06	monocytesBYothers	ModRNAseq_06	5.30E-02
ModArray_06	5.23E-05	5.95E-08	9.88E-52	3.21E-81	3.30E-09	myeoligoBYothers	ModRNAseq_08	1.92E-26
ModArray_10	2.17E-42	4.70E-61	0.0216352	9.59E-53	3.85E-17	astroBYothers	ModRNAseq_07	1.50E-142
ModArray_14	1.66E-10	6.97E-10	3.05E-70	2.09E-43	5.43E-30	endoBYothers	ModRNAseq_18	8.54E-85
ModArray_15	1.09E-64	1.50E-84	2.85E-17	1.36E-25	0.00020016	neuronBYothers	ModRNAseq_02	NA
ModArray_20	1.61E-09	9.87E-11	2.28E-12	1.24E-11	1.49E-19	astroBYothers	ModRNAseq_07	1.50E-142

Results

SVA + LM approach reduces batch effects across AD, PD, HD and MS microarray and RNA-seq studies

To identify neurodegenerative disease specific transcriptome, publicly available control and neurological disease (AD, PD, HD and MS) brain gene expression datasets (microarray and RNA-seq) were normalized in preparation for meta-analysis. SVA + LM model's ability to effectively reduce batch effects or cross study variability in microarray and RNA-seq, was visualized by comparing pre-SVA + LM normalized (before) and SVA + LM normalized (after) expression data plots. Density plot comparison of expression data colored by study, shows a greater overlap, indicating lower variability between studies after SVA + LM normalization (Figure 1A for microarray, Figure 2A for RNA-seq). In Box-Whiskers plot, slightly skewed mean expression in GSE33000 and GSE108000 microarray studies, and GSE67333 RNA-seq study, were fixed by SVA + LM normalization (Figure 1B for microarray, Figure 2B for RNA-seq). Wide dispersion of samples from same study (same color) in PCA plot was reduced by SVA + LM normalization (Figure 1C for microarray, Figure 2C for RNA-seq). Correlation plots used to evaluate correlation of gene expression values between studies, showed increased positive correlation after SVA + LM normalization (Figure 1D for microarray, Figure 2D for RNA-seq). Thus, SVA + LM normalization reduced variability unrelated to disease in AD, PD, HD, MS and control gene expression datasets (microarray and RNA-seq), making the data suitable for meta-analysis.

Figure 1

Effect of SVA + LM adjustment on AD, PD, HD, MS and control brain microarray datasets (GSE33000, GSE20333, GSE26927, GSE20164, GSE20292, GSE108000 and GSE43490). A: Density plot representation of gene expression per study before and after SVA + LM adjustment B: Box and whisker plot representation of gene expression per study before and after SVA + LM adjustment. C: Principal component (PC) projection of datasets before and after SVA + LM adjustment. D: Pearson correlation plot of studies before and after SVA + LM adjustment. For this figure the following R packages were used: prcomp built-in R function in stats4_3.5.0, corrplot_0.84, ggplot2_3.1.1, sva_3.30.1 and limma_3.38.3 [13,34,35,37,45, 46, 47].

Figure 2

Effect of SVA + LM adjustment on AD, PD, HD, MS and control brain RNA-seq datasets (GSE53697, GSE64810, GSE68719, GSE100297 and GSE67333). A: Density plot representation of gene expression per study before and after SVA + LM adjustment B: Box and whisker plot representation of gene expression per study before and after SVA + LM adjustment. C: Principal component (PC) projection of datasets before and after SVA + LM adjustment. D: Pearson correlation plot of studies before and after SVA + LM adjustment. For this figure the following R packages were used: prcomp built-in R function in stats4_3.5.0, corrplot_0.84, ggplot2_3.1.1, sva_3.30.1 and limma_3.38.3 [13,34,35,37,45, 46, 47].

WGCNA network analysis reveals a common neurodegenerative disease (AD, PD and HD) and MS specific gene co-expression module

Gene co-expression networks were constructed with WGCNA to uncover molecular mechanisms in neurodegenerative diseases (AD, PD and HD) and MS. WGCNA network analysis of microarray data revealed 23 modules, where ModArray_01 was largest module of size 2141 genes and ModArray_22 was smallest module of size 29 genes (Table1). WGCNA network analysis of RNA-seq data revealed 29 modules, where ModRNAseq_01 was largest module of size 3240 genes and ModRNAseq_28 was smallest module of size 105 genes (Table1). Spearman correlations between module eigengenes and disease trait (AD, PD, HD and MS) in microarray dataset, showed 22 WGCNA modules significantly associated with disease trait (Figure 3A,B). Spearman correlation between module eigengenes and disease trait (AD, PD, HD and MS) in RNA-seq dataset, showed 27 WGCNA modules significantly associated with disease trait (Figure 3C,D). Few WGCNA modules were distinctly associated with only one of four diseases, such as PD only module ModArray_17 (microarray), MS only modules ModRNAseq_04/13/16/23/25/28 (RNA-seq), HD only module ModRNAseq_19 (RNA-seq), AD only module ModRNAseq_20 (RNA-seq) (Figure 3B,D). Widespread 100% preservation of microarray and RNA-seq WGCNA modules was observed (Figure 3E). Hypergeometric enrichment test showed significant overlap of disease specific (AD, PD, HD and MS) WGCNA modules from microarray and RNA-seq datasets (Figure 3F). Thus, microarray and RNA-seq WGCNA gene network analyses complemented each other in identifying AD, PD, HD and MS disease specific modules, that were largely preserved.

Table 1

Number of genes in microarray and RNA-seq WGCNA modules.

	WGCNA_microarray _module_names	Number_of_genes		WGCNA_RNA-seq _module_names	Number_of_genes
1	ModArray_01	2141	1	ModRNAseq_01	3240
2	ModArray_02	1486	2	ModRNAseq_02	1850
3	ModArray_03	783	3	ModRNAseq_03	1690
4	ModArray_04	335	4	ModRNAseq_04	1668
5	ModArray_05	297	5	ModRNAseq_05	1638
6	ModArray_06	245	6	ModRNAseq_06	1586
7	ModArray_07	194	7	ModRNAseq_07	1177
8	ModArray_08	186	8	ModRNAseq_08	1115
9	ModArray_09	172	9	ModRNAseq_09	1112
10	ModArray_10	139	10	ModRNAseq_10	956
11	ModArray_11	134	11	ModRNAseq_11	875
12	ModArray_12	108	12	ModRNAseq_12	812
13	ModArray_13	107	13	ModRNAseq_13	714
14	ModArray_14	99	14	ModRNAseq_14	669
15	ModArray_15	99	15	ModRNAseq_15	581
16	ModArray_16	80	16	ModRNAseq_16	424
17	ModArray_17	80	17	ModRNAseq_17	384
18	ModArray_18	66	18	ModRNAseq_18	314
19	ModArray_19	53	19	ModRNAseq_19	298
20	ModArray_20	52	20	ModRNAseq_20	294
21	ModArray_21	37	21	ModRNAseq_21	276
22	ModArray_22	29	22	ModRNAseq_22	247
23	ModArray_00	2	23	ModRNAseq_23	229
			24	ModRNAseq_24	224
			25	ModRNAseq_25	224
			26	ModRNAseq_26	200
			27	ModRNAseq_27	186
			28	ModRNAseq_28	105
			29	ModRNAseq_00	2

Differential gene expression analysis reveals brain cell type and immune cell type specific marker genes

To identify brain cell type (neuron, astrocyte, oligodendrocyte and endothelial) and immune cell type (microglia, monocyte and macrophage) specific signature genes, differential gene expression meta-analysis was performed with three DEG analysis methods (edgeR, limma and simple comparison of means). Comparison of astrocytes with other cell types revealed 1509 DEGs (fold change >1.25, p-value <0.05) in at least two of the three DEG methods (Figure 4A). Endothelial cells compared to other cell types showed 477 DEGs (fold change >1.25, p-value <0.05) in at least two of the three DEG methods (Figure 4B). Each immune cell type, when compared with other cell types showed 592 DEGs for macrophages, 344 DEGs for microglia and 863 DEGs for monocytes, with fold change >1.25 and p-value <0.05, in at least two of the three DEG methods (Figure 4C, D, 4E). Oligodendrocytes, when compared with other cell types showed 124 DEGs for mature oligodendrocytes, 128 DEGs for precursor oligodendrocytes and 157 DEGs for new oligodendrocytes, with fold change >1.25 and p-value <0.05, in at least two of the three DEG methods (Figure 4F-H). Comparison of neurons with other cell types revealed 889 DEGs (fold change >1.25, p-value <0.05) in at least two of the three DEG methods (Figure 4I). UserListEnrichment R function showed that only a) new oligodendrocyte and precursor oligodendrocyte b) myelinated oligodendrocytes and precursor oligodendrocyte c) macrophages and monocytes, shared significant number of common signature genes.

Figure 4

A–I: Venn Diagram gene overlap and heat-map of overlap significance for differential gene expression calculation methods (limma, edgeR and simple comparison of means). Upregulated DEGs for contrasts Astrocytes by others (A), Endothelial cells by others (B), Macrophages by others (C), Microglia by others (D), Monocytes by others (E), Myelinated oligodendrocytes by others (F), Precursor oligodendrocytes by others (G), new oligodendrocytes by others (H) and Neuron by others (I). (J) Only significant overlap of cell type genes identified (A–I) with another cell type are shown. For this figure the following R packages were used: GeneOverlap_1.18.0, VennDiagram_1.6.20, limma_3.38.3, edgeR_3.24.3 and simple comparison expression means [37, 38, 39, 41, 42].

Immune cell type (microglia, macrophage and monocyte) marker genes were enriched in AD, PD, HD and MS associated WGCNA module

Brain cell type (neuron, astrocyte, oligodendrocyte and endothelial) and immune cell type (microglia, monocyte and macrophage) specific marker genes identified by differential gene expression analysis were used to estimate relative enrichment of these cells in AD, PD, HD and MS WGCNA modules. Hypergeometric enrichment analysis revealed cell type enrichment in AD, PD, HD and MS associated WGCNA modules (Figure 5A,B, Table2). ModArray_03 microarray WGCNA module and ModRNAseq_14 RNA-seq WGCNA module, were most significantly enriched with neuron cell type markers (Table2). Interestingly, ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module, were the only modules that were significantly enriched with all three-immune cell type (macrophage, microglia and monocyte) markers (Table2, Figure 5A,B). ModArray_06 microarray WGCNA module and ModRNAseq_08 RNA-seq WGCNA module, were significantly enriched with myelinated oligodendrocyte cell type markers (Table2). ModArray_10 and ModArray_20 microarray WGCNA modules, and ModRNAseq_07 WGCNA module, were significantly enriched with astrocyte cell type markers (Table2). ModArray_14 microarray WGCNA module and ModRNAseq_18 RNA-seq WGCNA module, were significantly enriched with endothelial cell type markers (Table2). Though microarray WGCNA module ModArray_15 was enriched with neuron cell type signature genes, its equivalent ModRNAseq_02 RNA-seq WGCNA module was not significantly enriched with any cell type genes (Table2). Thus, each brain cell type occupies separate WGCNA modules, while all three immune cell types (macrophage, microglia and monocyte) occupy the same WGCNA module. This suggests that mechanistically each brain cell type utilizes a separate gene network or module, while all immune cell types utilize the same gene network or module.

Figure 5

Immune cell type (macrophage, microglia and monocyte) marker genes were enriched in WGCNA modules associated with AD, PD, HD and MS. A: Gene overlap between AD, PD, HD and MS associated WGCNA modules from microarray (ModArray_04) and RNA-seq (ModRNAseq_06 RNA-seq). B: Table showing most significant overlap of cell type signature genes in ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module. C, D: Gene Ontology (GO) Biological Process analysis of ModArray_04 microarray WGCNA module (C) and ModRNAseq_06 RNA-seq WGCNA module (D). For this figure the following R packages were used: EnrichR_1.0, WGCNA_1.66 and ggplot2_3.1.1 [34,43,44,48].

Gene ontology (GO) annotation of immune cell type enriched AD, PD, HD and MS associated WGCNA modules

Immune system was the most significant biological process gene ontology (GO) in both ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module (Figure 5C,D). Immune process, GO “neutrophil activation involved in immune response” (GO:0002283) was largest in both ModArray_04 (overlap 71/484 p-value 1.43E-42) microarray WGCNA module and ModRNAseq_06 (overlap 83/484 p-value 7.66E-08) RNA-seq WGCNA module (Figure 5C,D). Other large GOs in ModArray_04 microarray WGCNA module were related to cell signaling, such as “cytokine-mediated signaling pathway (GO:0019221) overlap 70/634 p-value 3.46E-34” and “positive regulation of intracellular signal transduction (GO:1902533) overlap 36/480 p-value 5.58E-14” (Figure 5C). Large GOs related to cell signaling were also present in ModRNAseq_06 RNA-seq WGCNA module, such as “protein phosphorylation (GO:0006468) overlap 57/471 p-value 0.000911972” and “toll-like receptor signaling pathway (GO:0002224) overlap 25/87 p-value 8.78E-09” (Figure 5D). These results are consistent with immune cell type enrichment in ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module.

Discussion

In MS patients, immune system attacks myelin surrounding neurons and is the most studied central nervous system immune dysfunction [49]. Several recent studies now indicate that a dysfunctional immune system also underlies pathology of other central nervous system diseases, such as neurodegenerative diseases AD, PD and HD [13, 16, 17, 29, 50, 51]. Recently, Mukherjee et al showed that resident brain immune cells called microglia, are a major cell type enriched in gene networks associated with aging and neurodegenerative diseases, AD, PD and HD [52]. Growing evidence suggests that with aging there is an increase in systemic inflammation and collapse of blood brain barrier [19, 20]. However, how immune dysfunction in neurodegenerative diseases compares to MS, a disease known to be caused by immune dysfunction, is not well understood. To address this question, regulatory molecular networks underlying neurodegenerative diseases (AD, PD and HD) and MS were constructed, and immune system characteristics of these networks were evaluated.

To enable comparison of disease effects on gene expression, transcriptome gene expression data were normalized with SVA + LM adjustment to remove effects of age and other covariants on gene expression. A caveat of this work is that association of RNA-seq WGCNA modules with PD was not always significant compared to microarray WGCNA modules. This is probably because RNA-seq PD samples were from cortex while microarray samples were form substantia nigra. However, in spite of this limitation cell type enrichment results were same for significant PD associated microarray WGCNA modules and RNA-seq WGCNA modules (Table2).

Different disease (AD, PD, HD and MS) associated WGCNA modules were enriched by each of the major brain cell type signature genes (neurons, astrocytes, oligodendrocytes and endothelial cells) (Table2). This suggests that all major brain cell type occupy separate AD, PD, HD and MS disease modules and gene networks. Interestingly, all immune cells type signature genes (microglia, monocyte and macrophage) were enriched in the same ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module, associated with AD, PD, HD and MS diseases (Table2, Figure 5A,B). Though macrophages and monocytes shared significant number of common signature genes, microglia signature genes did not significantly overlap with either of them (Figure 4J). This rules out that the three innate immune cell types occur in the same module simply because they have significant number of common signature genes. Taken together, these results suggest that all immune cells underlie a common shared network, which is distinctly different from networks occupied by other brain cell types in AD, PD, HD and MS. One caveat here is that the cell type signature genes and disease modules were obtained from different species, mice and human, respectively.

The AD, PD, HD and MS associated ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq WGCNA module, enriched by all three immune cells (macrophages, microglia and monocytes) was characterized by several immune system related GO terms, “neutrophil activation involved in immune response”, “cytokine-mediated signaling pathway”, “positive regulation of intracellular signal transduction”, “protein phosphorylation” and “toll-like receptor signaling pathway”. Taken together, these results suggest a significant role of immune system in the context of neurodegenerative diseases (AD, PD, HD) and MS. In this study, the module or network where immune system is significantly involved in these diseases has been identified. Though more research is required to understand the interaction between immune system and neurodegenerative diseases (AD, PD and HD), one study suggests that microglia are activated by monocytes and macrophages through systemic chronic inflammation in neurodegenerative diseases [53].

Drug targeting strategies involve systemic targeting (specific targeting of organs, tissue and cells) and intracellular targeting (specific targeting of genes and pathways). The basis of network medicine is that diseases occur due to perturbation of entire gene networks and not due to a few individual genes [54]. Several clinical trials targeting the immune system have been undertaken in AD, PD, HD and MS patients. In these trials, anti-inflammatory drugs were found to be ineffective in AD patients, while inhibitors targeting specific inflammatory signaling proteins were found beneficial [55, 56]. A recent study showed that immune-suppressor drugs, inosine monophosphate dehydrogenase inhibitors and corticosteroids, lowered the risk of PD [57]. In this work, several immune pathways were identified in immune cell type enriched in AD, PD, HD and MS associated ModArray_04 microarray WGCNA module and ModRNAseq_06 RNA-seq module. These pathways and WGCNA modules can potentially serve as excellent candidate for intracellular drug targeting in AD, PD, HD and MS treatment. Additionally, identification of a common underlying immune response in AD, PD, HD and MS diseases in this study, presents a great opportunity to repurpose drugs developed for one of the diseases for another disease. An example of successful drug repurposing, is the approval of a MS treatment drug for immunomodulator targeting macrophages called fumaric acid ester dimethylfumarate (DMF) for HD treatment [58].

In 2015, Durrenberger et al found common immune pathways associated with multiple neurological diseases (AD, PD, HD, ALS, MS and schizophrenia) using microarray datasets [30]. Though this work requires validation with RNA-seq datasets, it raises the exciting possibility that immune cell types (microglia, monocytes and macrophages) associated with neurodegenerative diseases (AD, PD and HD), may also be associated with non-neurodegenerative psychiatric diseases, such as schizophrenia. More research is required to understand the role of innate and adaptive immune system in neurodegenerative diseases, and non-neurodegenerative diseases to identify a safe and effective therapeutic window for immune system drug targeting.

Declarations

Author contribution statement

Shradha Mukherjee: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

For datasets: AD, PD, HD, MS and control microarray (GSE33000, GSE20333, GSE26927, GSE20164, GSE20292, GSE108000 and GSE43490), and RNA-seq (GSE53697, GSE64810, GSE68719, GSE100297 and GSE67333) datasets were obtained from NCBI GEO.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.