The microarray dataset was downloaded from Gene Expression Ominibus (GEO)¹. The dataset, GSE15222, measured on GPL2700 using Illumina Sentrix Human Ref-8 Expression Bead Chip which covers 24354 probesets consists of human cerebral cortex samples of 188 controls and 176 patients with diagnosis of late-onset AD. The original expression data were rank-invariant by BeadStudio software available from Illumina and the Illumina custom error model was used. Rank-invariant-normalized expression data were log₁₀ transformed, and missing data were encoded as missing, rather than as a zero level of expression. Chips with average detection scores less than 0.99 (5% of control chips, 8% of late-onset AD chips) were excluded from the analysis. Transcripts that were detected in less than 90% of the case or 90% of the control series were excluded from our study. Finally 8650 high quality transcripts were obtained (Webster et al., 2009). Of the 188 controls, 40 were APOE ε4 carriers and 148 were non-carriers. Of the 176 patients, 121 were APOE ε4 carriers and 55 were non-carriers.

The WGCNA package which provides a robust set of R functions was used to detect sub-dataset-specific modules (Langfelder and Horvath, 2008). GSE15222 was divided into four sub-datasets: patients carrying APOE ε4, patients not carrying APOE ε4, controls carrying APOE ε4 and controls not carrying APOE ε4. Outliers were filtered out using Euclidian distance as the similarity measure and average linkage as an agglomeration method. For each sub-group, subjects were clustered based on their dissimilarity, and any arrays with average inter-subject correlation less than 2 standard deviations below the mean were removed. This process was repeated until no arrays needed to be removed (Oldham et al., 2008). One outlier was detected in AD APOE ε4 non-carriers and was excluded before the construction. A demographic depiction of the subjects was performed to examine the similarity of sub-datasets. To construct the network, Pearson correlation coefficients were calculated for all possible gene pairs, and then the coefficients were powered by an exponent β in each sub-group stratified by APOE ε4 and phenotypes. A high β maintains high adjacencies but pushes lower adjacencies toward zero. An optimal exponent β can reduce the false positive rate within the network to the utmost extent (Zhang and Horvath, 2005). The dynamic tree-cutting algorithm was then used to identify modules after the summarization of similar patterns of connectivity among genes as topological overlap matrix (TOM) (Langfelder et al., 2008). Modules with similar module eigengenes (ME) were merged together. Unsupervised average linkage hierarchical clustering identified 18 modules in AD APOE ε4 carriers, 21 modules in AD APOE ε4 non-carriers, 10 modules in control APOE ε4 carriers and 8 modules in control APOE ε4 non-carriers. Next, consensus modules (AD APOE ε4 carriers and non-carriers, AD and control APOE ε4 carriers, and AD and control APOE ε4 non-carriers) were constructed with a similar procedure to that illustrated above. A detailed description of consensus module construction can be found at the WGCNA tutorial website².

The modules constructed using the sub-dataset of AD APOE ε4 carriers were compared with the consensus modules of AD APOE ε4 carriers and non-carriers, and with the consensus modules of AD and control APOE ε4 carriers to detect specific modules of AD APOE ε4 carriers. The modules constructed using the sub-dataset of AD APOE ε4 non-carriers were compared with the consensus modules of AD APOE ε4 carriers and non-carriers, and with the consensus modules of AD and control APOE ε4 non-carriers to detect specific modules of AD APOE ε4 non-carriers. Then, the module preservation statistic Z_summary was utilized to further examine the sub-dataset-specificity of the modules (Langfelder et al., 2011). Unlike the cross-tabulation test, Z_summary takes into account both overlaps in module membership and the density and connectivity patterns of the modules. The following recommended significant thresholds for Z_summary were adopted: Z_summary <2 implies no evidence of module preservation, 2 < Z_summary < 10 implies weak to moderate evidence of module preservation, and Z_summary > 10 implies strong evidence for module preservation. Greater module preservation corresponds to lesser specificity of the module to the sub-dataset. Thus, a Z_summary less than 2 indicates strong evidence of specificity of the module, a Z_summary between 2 and 10 indicates modest evidence of specificity of the module and a Z_summary more than 10 indicates no evidence of specificity of the module. Module networks were graphically depicted using the program Cytoscape³.

Then random samplings with two different levels (50 and 75%) were performed with each level carried out five times. The criterion of replication is that more than 80% of genes within one of the modules (violet, dark magenta and light cyan) detected among the total samples should be included in one module detected in random samplings. Moreover, modules that were obtained from random samplings and whose genes are overlapped more than 80% with violet, dark magenta or light cyan modules should be AD APOE ε4 carrier or non-carrier-specific.

Hub genes were identified by WGCNA via measures of intramodular connectivity and module membership (Zhang and Horvath, 2005; Dewey et al., 2011; Langfelder et al., 2011; de Jong et al., 2012). Intramodular connectivity measures how a given gene is connected, or co-expressed, with the genes of a particular module. Module membership measures the membership of the i-th gene with respect to a given module. Hub genes tend to own high values of intramodular connectivity and module membership.

Module eigengene values adjusted for phenotype, age and sex were used to test differential expression of conserved modules between AD and control APOE ε4 carriers and between AD and control APOE ε4 non-carriers, and the values adjusted for APOE ε4 status, age and sex were used to test differential expressions of conserved modules between AD APOE ε4 carriers and non-carriers.

R package graph was utilized to detect whether there was/were complete graph(s) comprised of hub genes in sub-dataset-specific modules.

Primary neurons were derived from the hippocampus of Sprague–Dawley (SD) rats at postnatal day 1 as previously described (Kam et al., 2010). Neurons were plated at a density of 2 × 10⁶ cells/well on poly-L-lysine-treated 6-well plate. The neurons were treated with recombinant human APOE ε4 or ε3 at a concentration of 5 μg/ml or without APOE treatment on DIV19 and collected for RNA extraction after 24, 48, and 72 h of treatment. Total RNA was extracted using a Total RNA kit 1 from Sigma (St. Louis, MO, USA). RNA was reverse transcribed to synthesize cDNA by using M-MLV reverse transcriptase (Invitrogen, Carlsbad, CA, USA). Real-time PCR was performed using SYBR Green (Roche Diagnostics, Basel, Switzerland). GAPDH was used as the normalization control. The relative mRNA levels were calculated by a comparative Cp value. The primers used were: ENO3: sense: AGCTGCTACCTAGGCAC TCT, antisense: GGTTCCGTCCAGCTCAATCA; GNB3: sense: TTTCACTGGCCACGAGTCAG, antisense: CTCTCGTGGGAG TAGGCTGT; XPO4: sense: GGAATTCAGCAGACGGGAGA, antisense: TAGTGTTTTGGAGGGAGAAAATTCC; GDF10: sense: AATCATCAAGGCTGCCCGAA, antisense: CTGGACC AGAACTCGTGCTT; ISOC1: sense: GCTGCACTAACAAAAC GCCA, antisense: TCATGGGACGGCAGGATAGA; ACLY: sense: GGTAAGCTGGTGCTTACGGA, antisense: TCTGGA TGGCTGAGGTGGTA; DUSP5: sense: GCCGACATTAGCTCC CACTT, antisense: GCCAAAGTTGGGAGAGACCA; TNFR SF18: sense: CACGTGTCCCCGAGATACC, antisense: GTCCC CCAGACGACACTTTT; DNAJB5: sense: CTCCTCACCGCA GCACC, antisense: TTCTCCTCAGCGTTGGGTTC; HILPDA: sense: GCCTGCACGATCTAGTGTGA, antisense: GCACTCC TCTGGATGGATGG; TBC1D8: sense: TGTATTCTCCCATAGC ATGTGGT, antisense: GTCTCCTCAGCGATCAGAGC; BDNF: sense: ACTGTCCTGCTACCGCAGTTG, antisense: GGGTCGC AGAACCGCTAAA; ST8SIA5: sense: ACTTCGTCTTCCGGTG CAAT, antisense: GGAAGTCGTCCAGCATGTACT; CHST12: sense: GGTCTCCTTCGCCAACTTCA, antisense: AGCATCC TCATCCAGGGTCT; ATF5: sense: GTGCCTAGGGTACAGGA GGA, antisense: GCAGAGGGGAGACCTAGACA; CSNK1D: sense: CACCTCACAGATTCCCGGTC, antisense: GCTCTTG GAGCCTGTCCATT; PLIN2: sense: ATTCGCCAGGAAGAAT GTGC, antisense: TGGCATGTAGTGTGGAGCTG.

Functional annotations of sub-dataset-specific modules were obtained by applying ingenuity pathway analysis (IPA)⁴. IPA is a large manually curated database of published information on mammalian biology and diseases. It is widely used for high quality gene set enrichment analyses. Fisher’s exact tests were utilized to calculate the P-values for each functional annotation by comparing the number of genes from the module of interest that participated in the specified IPA term against the total number of genes from the term in the background set (Khatri et al., 2012).

The GWAS (Myers) was carried out in the same population as that of the microarray dataset we utilized, and the dataset was directly downloaded from the Myers laboratory⁵. For replication, the GWAS dataset from the Multi-Site Collaborative Study for Genotype-Phenotype Associations in Alzheimer’s Disease (GenADA) (Li et al., 2008) was downloaded from dbGaP⁶, the GWAS dataset from the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (Mueller et al., 2005) was downloaded from the ADNI database⁷, and the dataset from the International Genomics of Alzheimer’s Project (IGAP) (Lambert et al., 2013) was downloaded from http://www.pasteur-lille.fr/en/recherche/u744/igap/igap_download.php.

We excluded single nucleotide polymorphisms (SNPs) with a missing genotype rate > 0.1, a minor allele frequency (MAF) < 0.01 and a Hardy–Weinberg equilibrium (HWE) < 0.001. Finally, after quality control, 364048, 431284, and 567199 SNPs from the Myers, GenADA and ADNI datasets, respectively, were included for APOE ε4 carriers, and 374149, 427678, and 566459 SNPs from the Myers, GenADA and ADNI datasets, respectively, were included for APOE ε4 non-carriers. GWAS dataset preparation and quality control procedures were performed using the software package PLINK⁸, release v 1.07 (Purcell et al., 2007). Processed data of IGAP were downloaded from the following link: http://www.pasteur-lille.fr/en/recherche/u744/igap/igap_download.php.

A previously reported procedure of GWAS signal enrichment test, called i-GSEA4GWAS, was employed (Zhang et al., 2010). First, the maximum -log(P-value) of the SNPs located between 20 kb upstream and downstream of a gene boundary was assigned to represent the gene, and the min P-value method, which does not excessively penalize large genes if causative associations are proportionally more common in smaller genes, has been widely used in systems biology (Baranzini et al., 2009; Jia et al., 2010). Instead of the commonly used phenotype label permutation, SNP label permutations were implemented to generate the distribution of the enrichment score (ES), and then the gene set’s significance proportion based enrichment score (SPES) was calculated based on genes’ rank using the following equation: SPES = k/K ×ES, where k is the proportion of significant genes of the gene set and K is the proportion of significant genes of the total number of genes in the GWAS (Zhang et al., 2010). Gene-length bias was eliminated by applying adaptive permutation in PLINK before i-GSEA4GWAS. False discovery rate was applied for multiple comparison correction. Then, a comparative quantile-quantile (QQ) plot was used to demonstrate the differences among all genes and genes within the interested modules of APOE ε4 carriers, and the genomic dispersion factor, λ, was used to assess the strengths of genetic association signals and to quantify the differences when compared with the signals of all genes. The comparative QQ plot was produced by plotting the ranked -log₁₀(P-value) against the expected order statistic, -log₁₀[i/(L+1)], where i is the rank for each SNP P-value (1 = smallest and L = largest), and L is the number of SNPs. Function estlambda() within the R package GenABEL was utilized to estimate the genomic dispersion factor λ. If a co-expression module was genetically associated, significant SNPs within 20 kb upstream and downstream of the gene boundaries within the module in the sub-group in which it belongs (APOE ε4 carrier sub-group or APOE ε4 non-carrier sub-group) were identified, and the significance of these SNPs in its opposite sub-group were also calculated.

The conserved promoter motif(s) was/were identified by PhyloCon. PhyloCon, which stands for Phylogenetic Consensus, is one of the first motif-finding algorithms to combine the power of phylogenetic conservation and gene co-regulation (Wang and Stormo, 2003; Wang, 2007). PhyloCon first aligns conserved regions of orthologous sequences into multiple sequence alignments, or profiles, and then compares profiles representing non-orthologous sequences. Then, motifs represented by matrices or IUPAC strings emerge as common regions in these profiles (Wang, 2007). An online web server named WebLogo was used to generate the sequence logos of the motifs (Crooks et al., 2004). In the present study, orthologous sequences of the species Homo sapiens, Rattus norvegicus, Mus musculus, and Canis lupus familiaris were utilized. For promoter regions, the 4 kb segments centered on the annotated transcription start site (TSS) of each human RefSeq gene were extracted according to a previous study’s method (Xie et al., 2005). If the annotated translation start codon was within 2 kb of the TSS, the shorter region that did not overlap the protein-coding sequence was selected. For genes with alternatively spliced first exons, all promoters were included. Details for promoter region detection are given in a previous study (Xie et al., 2005).

JASPAR CORE database and miRWalk were utilized to identify common human transcription factors (TFs), TF binding sites (TFBS) and common miRNAs of the conserved motif(s) and these highly connected genes. The JASPAR CORE database, which contains a curated, non-redundant set of profiles, is derived from published collections of experimentally defined transcription factor binding sites for eukaryotes, The matrices and IUPAC strings were then input into the JASPAR CORE database to identify common human TFs and the TF binding sites (TFBS) of these co-expressed genes (Mathelier et al., 2014). miRWalk, a comprehensive database that incorporates miRNA-targets interactions information produced by eight established miRNA prediction programs on 3′ UTRs of all known genes, was utilized to evaluate whether some miRNAs exist as the common regulators of the mRNAs of genes within the same complete graph post-transcriptionally (Dweep et al., 2011). In the present study, five default prediction programs, miRanda, miRDB, miRWalk, RNA22, and TargetScan, were used. If a miRNA was predicted by at least two prediction programs, it was regarded to bind to the gene with high probability.

GEO dataset GSE20295⁹ for PD and GSE12654¹⁰ for BD were extracted. The procedures for PD and BD co-expression network construction and module detection were the same as those of AD. Then, each PD- or BD-specific module was compared with AD-specific modules via hypergeometric tests which examined whether the number of genes overlapped between PD or BD and AD was significantly larger than that expected by chance.

The demographic characteristics of the enrolled subjects after excluding outlier samples demonstrated no differences (Supplementary Table S1). By applying WGCNA, we identified two specific co-expression modules (violet and dark magenta) in AD APOE ε4 carriers and one module (light cyan) in AD APOE ε4 non-carriers (Figure 1). The permutation-based preservation statistic, Z_summary, demonstrated that the Z_summary statistics of the identified three modules were all less than 2 when assessing the violet and dark magenta modules of AD APOE ε4 carriers in AD APOE ε4 non-carriers and control APOE ε4 carriers and when assessing the light cyan module of AD APOE ε4 non-carriers in AD APOE ε4 carriers and control APOE ε4 non-carriers, showing no evidence of module preservation in these sub-datasets (Figure 1; Supplementary Table S2). These results strongly indicated that the three co-expression modules were sub-dataset-specific. Network visualizations of the violet and dark magenta modules in AD APOE ε4 carriers and the light cyan module in AD APOE ε4 non-carriers are shown in Figures 2A–C, respectively (Supplementary Tables S3–S5 for full list of genes of violet, dark magenta and light cyan modules respectively). The identified hub genes (most centered genes) are shown in Supplementary Tables S3–S5. Moreover, six of the seven hub genes (ISOC1, ENO3, GDF10, GNB3, XPO4 and ACLY, Supplementary Table S3) of the violet module and all the hub genes of the light cyan module (Supplementary Table S5) were found to be highly connected as demonstrated to form a complete graph respectively, which represents the most highly connected gene cluster in biological network and tends to be very important (Jeong et al., 2001; Horvath and Dong, 2008) (Figure 2D; Supplementary Table S6 for the pairwise Pearson correlation coefficients for the violet module and Figure 2E; Supplementary Table S7 for the pairwise Pearson correlation coefficients for the light cyan module).

Furthermore, we examined whether the non-specific co-expression modules were differentially expressed between AD and control APOE ε4 carriers, between AD and control APOE ε4 non-carriers, and between AD APOE ε4 carriers and non-carriers. In the results, no differences were found (data not shown), demonstrating that these non-specific co-expression modules were completely conserved.

Due to the unavailability of other AD-associated datasets with APOEε4 status, we carried out random samplings for replication. All three modules could be replicated (Supplementary Table S8).

The highly connected hub genes in the complete graph implicate that the expression of these genes is correlated. Thus, we examined the co-expression patterns of these genes in primary cultured neurons treated with APOE ε3 or APOE ε4. As shown in Figure 3A, the mRNA expressions of all six genes (ISOC1, ENO3, GDF10, GNB3, XPO4 and ACLY) in the complete graph of the violet module of AD APOE ε4 carriers demonstrated a similar time-dependent patterns in neurons treated with physiological concentration of APOE ε4, whereas, no such patterns were observed in neurons treated with APOE ε3 or in untreated neurons (Figures 3B,C). The pairwise Pearson correlation coefficient matrix of the expressions showed the high correlations between either two genes expressed in the APOE ε4-treated neurons (Supplementary Table S9) with a mean absolute Pearson correlation coefficient of 0.902, which is markedly higher than the mean absolute Pearson correlation coefficients for neurons treated with APOE ε3 (0.613) or for untreated neurons (0.571) (Supplementary Tables S10 and S11 for the pairwise Pearson correlation coefficient matrices for APOE ε3-treated neurons and for untreated neurons, respectively). These results confirm the co-expression of these six genes detected in the violet module.

The mRNA expression of 12 genes in the complete graph of light cyan module of AD APOE ε4 non-carriers displayed similar or opposite time-dependent patterns in primary cultured neurons treated with APOE ε3 (Figure 3E). The pairwise Pearson correlation coefficient matrix showed positive or negative correlations between most gene pairs (Supplementary Table S12), with a mean absolute Pearson correlation coefficient of 0.834. These results validated the co-expression of genes in the complete graph of the light cyan module. The co-expression patterns were much weaker in neurons treated with APOE ε4 and in untreated neurons (Figure 3D, Supplementary Table S13 for the pairwise Pearson correlation coefficient matrix for APOE ε4-treated neurons and Figure 3F, Supplementary Table S14 for untreated neurons, the mean Pearson correlation coefficients were 0.616 and 0.629, respectively).

Co-expressed genes in a module may directly interact with each other or take part in the same biological processes and signaling pathways. Hub genes within each module may play pivotal roles in module function. Thus, IPA was employed for further analysis.

As shown in Figure 4A, the genes of the violet module identified in AD APOE ε4 carriers were enriched in hereditary disorders, neurological diseases, psychological disorders, and nervous system development and function (Supplementary Table S15). The first two hub genes, ISOC1 and ENO3, were enriched in the first annotation, hereditary disorders. These two hub genes and another hub gene GNB3 were enriched in neurological diseases (Supplementary Table S16). The primarily enriched signaling pathways were acetyl-CoA biosynthesis III (from citrate), NAD biosynthesis III, NAD salvage pathway III and NAD biosynthesis from 2-amino-3-carboxymuconate semialdehyde (Figure 4B; Supplementary Table S17), which are mainly energy metabolism-associated signaling pathways. Two relatively important energy metabolism-associated signaling pathways, gluconeogenesis I and glycolysis I, were enriched by the hub gene ENO3 (Supplementary Table S17). Genes of the dark magenta module were mostly enriched in hereditary disorders, neurological diseases, psychological disorders, and G protein-coupled receptors (GPCRs) and second messenger-associated signaling pathways of GPCRs (Figures 4C,D, see Supplementary Tables S18–S20 for all enriched IPA terms). The protein level of the hub gene, GAD2, which catalyzes glutamate to γ-aminobutyric acid, has been reported to be reduced in some cerebral regions of AD patients (Schwab et al., 2013).

Unlike the modules identified in AD APOE ε4 carriers, genes of the light cyan module identified in AD APOE ε4 non-carriers were enriched in some other types of diseases, such as immunological and cardiovascular diseases, in addition to the primarily enriched neurological diseases (Figure 4E, see also Supplementary Table S21). The top-ranked hub gene, DUSP5, was enriched in neurological diseases (Supplementary Table S22). Moreover, another hub gene, BDNF, has been widely reported to be related to AD (Tapia-Arancibia et al., 2008). Our analysis also demonstrated that approximately 66% of enriched genes found in either immunological or cardiovascular diseases were also found in neurological diseases (Supplementary Figure S1). However, genes within the two modules of AD APOE ε4 carriers were not highly enriched in immunological and cardiovascular diseases (Supplementary Tables S16 and S19). In addition, genes of the light cyan module were mainly enriched in the following signaling pathways: tRNA charging, endoplasmic reticulum (ER) stress pathway, TR/RXR (thyroid hormone receptor/retinoid X receptor) activation and TNFR1 (tumor necrosis factor receptor type 1) signaling (Figure 4F; Supplementary Table S23).

Next, we examined whether the three sub-dataset-specific co-expression modules had genetic bases. The violet module in AD APOE ε4 carriers showed significant enrichment of signals from the GWAS dataset of Myers of the corresponding APOE ε4 carriers in the gene expression data [P < 0.001 and P < 0.001 after false discovery rate (FDR) correction, Figure 5A]. Moreover, the significant enrichment of signals was replicated in the other two GWAS datasets of APOE ε4 carriers (P = 0.014 for the data from GenADA and P = 0.034 for the data from ADNI, Figures 5B,C). The dark magenta module of AD APOE ε4 carriers only showed a significant enrichment of signals from the GenADA GWAS dataset of APOE ε4 carriers (P = 0.015, Supplementary Figure S2). Comparative QQ plots showed that genes within violet and dark magenta modules, especially the violet module, deviated from the expected values even further compared with the black dots, which represent all genes (Figures 5D–F). The larger genomic dispersion factors (λ_violet and λ_darkmagenta) further validated their genetic bases (Figures 5D–F). Furthermore, genes in these two modules showed significant enrichment of signals from the IGAP GWAS dataset, which is the biggest AD-associated GWAS dataset to date (FDR-corrected P-value < 1 × 10^-4 for the violet module and P-value = 0.009 for the dark magenta module, Supplementary Figure S3). The QQ plots also showed a deviation of genes from the expected values (Supplementary Figure S4). Moreover, almost all significant SNPs within the 20 kb gene boundaries of the violet module in APOE ε4 carriers were insignificant in non-carriers (Supplementary Table S24). The light cyan module in AD APOE ε4 non-carriers did not show a significant enrichment of signals.

Complete graphs were detected and validated in the violet module of AD APOE ε4 carriers and in the light cyan module of AD APOE ε4 non-carriers, suggesting common transcription modulation. Thus, we first identified conserved promoter motif(s) by employing promoter homologous sequences of these genes from four species, Homo sapiens, Rattus norvegicus, Mus musculus, and Canis lupus familiaris. Two conserved promoter motifs with extremely low conservative P-values (P = 6.539 × 10^-122 for the 49-base motif and P = 1.003 × 10^-74 for the 14-base motif) were detected for the genes of the complete graph of the violet module of AD APOE ε4 carriers (Figures 6A,B), and an adenine-rich conserved promoter motif was identified with a P-value of 1.935 × 10^-125 among the genes of the complete graph of the light cyan module of AD APOE ε4 non-carriers (Figure 6C).

Then, common human TFs of two motifs in the violet module were identified (Table 1). Interestingly, two TFs, ZNF263 and ESR1, were shared by the two motifs, and SNPs within ESR1 have been widely reported to be associated with AD (Mattila et al., 2000; Corbo et al., 2006). RREB1, a TF of 49-base motif, which can potentiate the transcriptional activity of NeuroD1/β2 (Ray et al., 2003), has been reported to be related with AD and potentially with the NFκB1 gene (Pasluosta et al., 2011). Another neuro-associated TF of 49-base motif, named REST, which can repress neuronal genes in non-neuronal tissues (Yeo et al., 2005), has been demonstrated to both silence and repress neuronal genes (Greenway et al., 2007). SNPs proximal to TP63, a TF of 14-base motif, which is present in the aging human hippocampus (Yang and Wang, 1994), are reported to be associated with brain morphometric measures of AD (Shen et al., 2010). Another TF of 14-base motif, named EGR1, up-regulates the PSEN2 gene in neuronal cells (Renbaum et al., 2003). However, currently no highly scored human TFs could be predicted for the motif identified in the light cyan module of AD APOE ε4 non-carriers.

Furthermore, we investigated whether the mRNAs of these genes in the complete graph could be regulated by some common miRNAs. For the six genes in the violet module of AD APOE ε4 carriers, hsa-miR-194, hsa-miR-199a-5p, hsa-miR-199b-5p, hsa-miR-30a, hsa-miR-30d, and hsa-miR-30e were indicated to bind with at least four of the five hub genes except GNB3 (Supplementary Table S25). hsa-miR-194 has been reported to be down-regulated in white matter of AD patients and to be negatively associated with neurofibrillary tangles in gray matter and in neuritic plaques and neurofibrillary tangles in white matter of AD patients (Wang et al., 2011). miRNAs in the hsa-miR-199 family have been shown to target genes involved in neurodegenerative diseases (Roshan et al., 2009). Five members of the miR-30 family are up-regulated in AD patients (Leidinger et al., 2013). However, currently no highly scored common miRNAs were predicted for the genes in the complete graph of the light cyan module of AD APOE ε4 non-carriers.

Parkinson’s disease is also a neurodegenerative disorder with Lewy bodies deposited in neurons. Approximately 41.3% of PD patients are complicated with dementia according to a large population-based investigation (Mayeux et al., 1992). Pathologically, many PD patients show senile plaques and fibrillary tangles within the cerebral cortex (Boller et al., 1980), and many AD patients display Lewy bodies in cortical and subcortical regions (Hansen et al., 1990). BD is associated with an increased risk of dementia (Wu et al., 2013). Some symptoms and neuropathology of AD and BD overlap, such as brain atrophy, cognitive impairment, and emotional disturbances (Rao et al., 2012), however, no plaques and fibrillary tangles have been reported to be characteristic of BD. Thus, to further validate the identified modules, we applied the WGCNA analysis to these two neurological diseases to investigate whether the three identified co-expression modules were related to these diseases.

Transcriptomic data from three brain regions– Brodmann area 9 (BA9) of the prefrontal cortex, the putamen (PT) and the entire substantia nigra (SN) of PD patients were available for analysis. Among the detected PD-specific co-expression modules, only the purple module detected in BA9 of the prefrontal cortex which is the most PD-specific module (Supplementary Figure S5) was significantly correlated with the light cyan module of AD APOE ε4 non-carriers with an extremely low P-value (Figure 7A, P = 9.522 × 10^-24, Supplementary Table S26 for detailed overlapped genes). Two of the overlapping genes, DUSP5 and TBC1D8, were hub genes of the light cyan module, and DUSP5 was the top hub gene of the light cyan module (Supplementary Table S5). Furthermore, many annotations were the same between the two modules. The top diseases and functions and canonical pathways enriched in the light cyan module were also enriched in the purple module (Figures 7B,C). However, no BD-specific modules were found to be related to the three modules identified in AD (data not shown).