Figure 1 shows the overall workflow of this study. All microarray datasets were downloaded from the GEO database¹. We searched the GEO database for microarray datasets using the keyword “Alzheimer.” Datasets were included if they met the following criteria: (1) were from humans; (2) included expression data from the temporal cortex of both AD and cognitively normal (CN) samples, expression data from the temporal cortex of AD samples with different Braak stages, or blood expression data from AD, MCI, and CN samples; (3) the number of rows in each platform was >30,000; (4) the number of AD samples was ≥10, and the number of CN samples was ≥10; and (5) there were no repeated samples among datasets. Finally, five datasets from the temporal cortex of AD and CN samples; one dataset from the temporal cortex of AD samples with different Braak stages; and two datasets from the blood of AD, MCI, and CN samples were selected. Detailed information for these datasets, including GEO accession ID, dataset country, sample numbers, platform ID, and number of genes in each platform, as well as usage in the current study and references, was recorded and is shown in Table 1. Because GSE132903 had far more samples (98 AD, 97 CN) than the other four datasets from the temporal cortex of AD and CN samples, combining this dataset would have blurred the results of the combined analysis. Thus, GSE132903 (Piras et al., 2019) was used to validate further the differential expression of hub genes. Datasets GSE122063 (28 AD, 22 CN; McKay et al., 2019), GSE36980 (10 AD, 19 CN; Hokama et al., 2014), and GSE5281 (16 AD, 12 CN; Liang et al., 2007, 2008a,b; Readhead et al., 2018) had similar sample sizes and were chosen for the combined analysis in the current study. Dataset GSE118553 (52 AD, 31 CN; Patel et al., 2019) was used to build the coexpression network and to perform GSEA. Dataset GSE106241 (Marttinen et al., 2019) was used to test the expression of hub genes in different Braak stages and to explore the correlation with β-secretase activity and Aβ42 levels. Datasets GSE63060 and GSE63061 (Sood et al., 2015) were used to validate the expression of hub genes in the blood of AD samples. Series matrix files of these datasets and their corresponding platform files were downloaded for the current analysis.

Combined analysis was performed for three datasets, including GSE122063 (McKay et al., 2019), GSE36980 (Hokama et al., 2014), and GSE5281 (Liang et al., 2007, 2008a,b; Readhead et al., 2018). Converting probes to gene symbols for series matrix files of each dataset and merging the gene expression data of these three datasets were conducted using Perl (version 5.18.4) (Wall, 1994, Wall et al., 2000). Batch normalization of the merged file was conducted using the ComBat method from the R package “sva” (Johnson et al., 2007; Leek et al., 2019; R Core Team, 2020). DEG screenings were conducted using the R package “limma” (Ritchie et al., 2015). P-values were adjusted using the false discovery rate (FDR) method. Genes with adjusted P <0.05 and | log₂ fold change (FC)| >0.5 were considered as DEGs in the combined analysis. A heatmap of all DEGs was made by the R package “pheatmap” (Kolde, 2019). The R package “OmicCircos” was used to show the chromosomal locations, as well as expression patterns of the top 100 DEGs from the combined analysis.

GO and KEGG pathway analyses were conducted utilizing the R package “clusterProfiler” (Yu et al., 2012). GO terms and KEGG pathways with adjusted P <0.05 were considered to be significant and were exhibited using the R package “GOplot” (Walter et al., 2015).

We used the DEGs from the combined analysis to perform WGCNA using expression data from GSE118553 (Patel et al., 2019). The R package “WGCNA” was used to conduct this analysis and to find clinical trait–related modules and hub genes among the DEGs (Langfelder and Horvath, 2008, 2012). To transform the adjacency matrix to a topological overlap matrix, a soft-threshold power with a scale-free R² near 0.9 and a slope near 1 was selected. We set the soft-threshold power as 5 (scale-free R² = 0.9, slope = -0.96), cut height as 0.25, and the minimal module size as 10 for network construction and module detection. The module with the highest correlation with AD was considered the key module. The network of the key module from WGCNA was generated by Cytoscape (version 3.7.1)². Genes in the key module were selected to perform GO and KEGG analyses to explore their biological functions. Hub genes in WGCNA were defined as those with a gene significance (GS) >0.4 and a module membership (MM) >0.8.

Expression data of hub genes extracted from Dataset GSE132903 (Piras et al., 2019) and the other four aforementioned datasets were utilized to validate the differential expression of these hub genes in the temporal cortex. Hub gene expression data, β-secretase activity, and Aβ42 levels in AD samples from GSE106241 (Marttinen et al., 2019) were used to validate the expressions of hub genes in different Braak stages and to investigate their correlations with β-secretase activity and Aβ42 levels. Moreover, GSE63060 and GSE63061 (Sood et al., 2015) were used to explore the differential expressions of hub genes in the blood of AD and MCI samples. Using the two datasets, receiver operating characteristic (ROC) curves were made by GraphPad Prism (version 8.0.0)³. Bar plots, box plots, and correlation analysis in this section were all generated using GraphPad Prism (version 8.0.0)³.

The animal study protocol was reviewed and approved by the Ethics Committee of Capital Medical University. 5xFAD mice express five AD-related mutations in human forms of APP and PSEN1, including three in APP (K670N/M671L, I716V, and V717I) and two in PSEN1 (M146L and L286V). Eight-month 5xFAD and wild-type (WT, non-transgenic littermates) mice with mixed genders were used in the experiments.

Total RNA was extracted from cortices of 8-month 5xFAD and WT mice by RNAsimple Total RNA Kit (#DP419, TIANGEN, China). RNA 1 μg was used in the following reverse transcription (#RR047, Takara, Japan). Quantitative real-time polymerase chain reaction (qRT-PCR) was done on Applied Biosystems ViiA^TM 7 Real-Time PCR System using TB Green^® Premix Ex Taq (#RR420, Takara, Japan). β-Actin was used as internal control, and relative expression was determined using 2^–ΔΔCT method. Primers were designed and synthesized by BGI Tech Solutions (Beijing Liu He) Co., Ltd. Sequences of primers were shown in Supplementary Table 1.

GSEA (version 4.0.3, Broad Institute, Inc., Massachusetts Institute of Technology, and Regents of the University of California) was conducted to explore the possible biological functions of the hub genes (Mootha et al., 2003; Subramanian et al., 2005). AD samples in GSE118553 (Patel et al., 2019) were divided into two groups according to the expression level of the hub genes. The database “c2.cp.kegg.v7.0.symbols.gmt” was chosen for enrichment. Terms with P <0.05 and FDR <0.25 were defined as significant.

The normality test and homogeneity of variance test were performed on data extracted from GEO datasets. Data that passed these two tests then underwent t-testing for comparisons between two groups and analysis of variance (ANOVA) testing for comparisons among three or more groups. After ANOVA analysis, a Dunnett multiple-comparisons test was used for post-hoc testing. Data that passed the normality test, but did not pass the homogeneity of variance test, underwent t-testing with Welch correction for comparisons between two groups and the Brown–Forsythe ANOVA test for comparisons among three or more groups. Data that did not pass the normality test underwent non-parametric testing, using the Kruskal–Wallis test and Dunn multiple-comparisons test for comparisons among three or more groups. GraphPad Prism (version 8.0.0) (see text footnote 4) was utilized to perform the above tests.

Datasets GSE122063 (McKay et al., 2019), GSE36980 (Hokama et al., 2014), and GSE5281 (Liang et al., 2007, 2008a,b; Readhead et al., 2018) were included in the combined analysis. After combined analysis, 850 DEGs (223 up-regulated and 627 down-regulated) were identified and are shown in heatmap and volcano plots (Supplementary Figures 1, 2 and Supplementary Table 2). The polarity of genes described as “up-regulated” or “down-regulated” in this article is with respect to AD vs. CN. The top 100 DEGs (according to | log₂FC|, including top 50 up-regulated genes, as well as top 50 down-regulated genes) were chosen to visualize their chromosomal locations and expression patterns across the three datasets used for combined analysis, as well as their logarithmic adjusted P-values shown in the inner layer (Figure 2). The top five up-regulated genes were SERPINA3, FCGBP, MAFF, SCIN, and CD163, whereas the top five down-regulated genes were CARTPT, SST, PCSK1, PPEF1, and NPTX2 (Figure 2).

All DEGs were utilized to perform GO and KEGG analyses, and the top five of these terms based on their adjusted P-values are shown in chord plots (Figure 3). We found enrichments in several biological process terms for GO analysis. The top five terms were neurotransmitter transport, synaptic vesicle cycle, neurotransmitter secretion, signal release from synapse, and vesicle-mediated transport in synapse, which are shown in Figure 3A. The top five cellular component terms for GO analysis were presynapse, glutamatergic synapse, synaptic membrane, axon part, and exocytic vesicle (Figure 3B). The top five molecular function terms for GO analysis were neurotransmitter receptor activity, voltage-gated ion channel activity, voltage-gated channel activity, ion gated channel activity, and gated channel activity (Figure 3C). For KEGG pathway analysis, DEGs were mostly enriched in neuroactive ligand-receptor interaction, nicotine addiction, GABAergic synapse, synaptic vesicle cycle, and amphetamine addiction pathways (Figure 3D). Further, we did enrichment analyses in up-regulated genes and down-regulated genes of AD brains separately and found that up-regulated genes were mostly enriched in regulation of angiogenesis, whereas down-regulated genes were mostly related to synaptic functions (Supplementary Figures 3, 4).

Expression data of 823 DEGs were extracted from GSE118553 (27 DEGs not available; Patel et al., 2019) and used to conduct WGCNA (Figure 4). By setting the soft-threshold power as 5 (scale-free R² = 0.9, slope = -0.96; Figure 4B and Supplementary Figure 5) and cut height as 0.25, we acquired eight modules (Figures 4C,D), with non-clustering DEGs in the gray module. Genes in each module are shown in Supplementary Table 3. From the heatmap of module–trait correlations, we found that the green module was the most highly correlated with AD (correlation coefficient = -0.57, P = 2 × 10^–8; Figures 4E,F). We also found that the blue, yellow, and brown modules had similar correlation coefficients with the green module. However, after comparing the relationships between GS and MM, we considered the green module (correlation coefficient = 0.74, P = 2 × 10^–6; Figure 5A and Supplementary Figure 6) to be the key module associated with AD and therefore analyzed it further in detail. The key module contained 31 genes (Figure 5A and Supplementary Table 3), and the interaction of the genes in the key module is shown in Figure 5B. To explore the potential biological functions of genes in the key module, GO and KEGG analyses were performed, and the enrichments are shown in Figures 5C–F. The results indicated that genes in the key module were mostly related to biological processes like learning or memory and cognition (Figure 5C) and enriched in pathways like nicotine addiction and retrograde endocannabinoid signaling (Figure 5F). Under the threshold of MM >0.8 and GS >0.4, we identified the following 19 down-regulated hub genes in the key module: AP3B2, CAMKK2, CHGB, CLSTN1, CRYM, GABRD, GPR158, KIAA0513, MAL2, NPTX1, NRXN3, PHYHIP, RASGRF1, RPH3A, SCN2A, SCN2B, SLC8A2, STMN2, and TERF2IP.

All 19 hub genes underwent expression validation in GSE132903 (Piras et al., 2019), GSE122063 (McKay et al., 2019), GSE36980 (Hokama et al., 2014), GSE5281 (Liang et al., 2007, 2008a,b; Readhead et al., 2018), and GSE118553 (Patel et al., 2019) datasets. Except for CRYM, GABRD, PHYHIP, SCN2B, and TERF2IP in GSE36980 and CAMKK2, PHYHIP, and GPR158 in GSE132903, all other hub genes were significantly down-regulated in AD samples from the five datasets (P < 0.05, 0.01, 0.001, or 0.0001; Figure 6 and Supplementary Figure 7). Among the 19 hub genes, we selected AP3B2, GABRD, GPR158, KIAA0513, and MAL2, which have been little studied on the associations with AD, in order to validate their expressions in cortices of 5xFAD and WT mice and in different Braak stages of AD samples (GSE106241; Marttinen et al., 2019). All five hub genes were down-regulated in cortices of 5xFAD mice compared to those of WT (Figure 7). We also found that the five hub genes were down-regulated as AD progressed, especially AP3B2, KIAA0513, and MAL2, which were significantly down-regulated in Braak III–IV and Braak V–VI when compared with Braak 0 (P < 0.05, 0.01, respectively, Figure 8A). β-Secretase activity and Aβ42 levels in GSE106241 were identified to be up-regulated during AD progression (Figure 8A). We found that all five hub genes were negatively associated with β-secretase activity (P < 0.0001, P = 0.001, P = 0.0004, P < 0.0001, and P = 0.0002, respectively, Figure 8B) and that AP3B2 and KIAA0513 were negatively associated with Aβ42 levels (P = 0.019, 0.032, respectively, Figure 8C). Next, we tested the differential expressions of the five hub genes in the two AD GEO datasets from the blood (data not shown; Sood et al., 2015) and found that KIAA0513 was significantly up-regulated in MCI and AD samples when compared with CN samples (P < 0.05, 0.01, 0.0001, respectively, Figure 9A). We also found that KIAA0513 had the ability to differentiate MCI and AD from CN in the two datasets (Figure 9B). In GSE63060, the area under the curve (AUC) for differentiating MCI and CN samples is 0.69 [95% confidence interval (CI) = 0.61–0.77], and AUC for AD and CN samples is 0.63 (95% CI = 0.56–0.70). In GSE63061, AUC for MCI and CN is 0.62 (95% CI = 0.55–0.69), and AUC for AD and CN is 0.58 (95% CI = 0.51–0.65). Furthermore, KIAA0513 was found to be enriched in neurons of healthy human brains, as determined using an online database AlzData (Xu et al., 2018; Supplementary Figure 8)⁴.

We conducted GSEA of the five little-studied hub genes in the expression data of GSE118553. AD samples in GSE118553 were divided into a “high-expression group” and a “low-expression group.” We acquired 17, 9, 1, 7, and 33 significant gene sets enriched in the expression groups of AP3B2, GABRD, GPR158, KIAA0513, and MAL2, respectively (Supplementary Table 4). According to normalized enrichment scores, genes in the low-expressed KIAA0513 group were mostly related to ribosome function, antigen processing and presentation, and graft-vs.-host disease (Figure 9C).