The GSE181279 (Xu and Jia, 2021) and GSE63063 (Sood et al., 2015) datasets were obtained from the Gene Expression Omnibus (GEO) database.¹ GSE181279 includes scRNA-seq data of peripheral blood mononuclear cells (PBMCs) from three patients with AD and two cognitively normal controls (NC), based on the GPL24676 platform. GSE63063 is a superseries composed of GSE63060 and GSE63061 datasets. GSE63060 includes gene-chip datasets of blood samples from 49 patients with AD, 39 patients with mild cognitive impairment (MCI), and 67 NC based on the GPL6947 platform. GSE63061 included gene-chip datasets of blood samples from 40 patients with AD, 30 patients with MCI, and 72 NC, based on the GPL10558 platform.

During preprocessing, we used the Seurat R package (v3.1.2) to filter out empty droplets (those containing only ambient RNA) based on the criteria of expressing fewer than 200 genes. We also applied stringent thresholds for unique molecular identifier (UMI) count and mitochondrial gene content to exclude low-quality cells. Cells expressing more than 30% mitochondrial genes were filtered out to avoid cells undergoing apoptosis or stress. Additionally, we applied a lower UMI threshold of 200 and an upper threshold of 8,000 genes. After normalizing the data using the Seurat R package (v3.1.2), the FindClusters function was used to identify the major clusters. Subsequently, t-distributed stochastic neighbor embedding (tSNE) (Pont et al., 2019) was used to visualize major clusters. Cell types were defined based on the reported marker genes (Supplementary Table 1) and differential expression analysis of unique markers. Intercellular communication was analyzed using the CellChat R package (Jin et al., 2021). Pseudotime trajectory analyses were performed using monocle2 R package (Qiu et al., 2017). The contribution of cells to AD was calculated using the Seurat R package.

Differences between AD and NC in cell types were analyzed using the Seurat R package. Differentially expressed genes (DEGs) were identified at P < 0.05. Differences between cell types were analyzed using the limma R package (Ritchie et al., 2015). Enrichment analysis of biological processes (BP) in the Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) for DEGs was performed using the ClusterProfiler R package (Yu et al., 2012).

The co-expression networks of the 2000 most variable genes were constructed using scRNA-seq with weighted gene correlation network analysis (WGCNA) (Langfelder and Horvath, 2008) and high-dimensional WGCNA. We constructed meta-cells for each cell and normalized their expression matrix. The soft threshold power was calculated, and the optimal β was selected to obtain a scale-free network. Then the adjacency matrix was transformed into a topological overlap matrix (TOM) and the co-expression network was constructed. Module eigengenes (MEs) were calculated to represent module expression levels. The correlation between genes and module eigengenes was calculated to obtain eigengene-based connectivity (kME) and to identify highly connected genes (hub genes) in each module. The top 25 hub genes in each module were scored using the moduleexprscor function. Correlations between the modules and cell types were calculated using Pearson’s correlation.

DEGs between the AD and NC groups were analyzed using the limma R package with a permutation-based P-value of < 0.05. Common DEGs were obtained by intersectional analysis of DEGs that were up - or downregulated in both the GSE63060 and GSE63061 datasets.

Enrichment analyses of the Gene Ontology and KEGG databases for common DEGs were performed using Metascape.² Gene set enrichment analysis (GSEA) was performed using the clusterprofiler R package to detect which gene sets were significantly enriched in AD. Adjusted P < 0.05 were considered statistically significant. The area under the receiver operating characteristic curve (AUC) values were calculated using the pROC package (Robin et al., 2011) for both the GSE63060 and GSE63061 datasets. Common DEGs with the top ten AUC values were selected as candidate genes.

We performed differential analyses between AD and MCI or between MCI and NC samples. The generated genes were provided as inputs to STEM (Ernst and Bar-Joseph, 2006) to observe variations in gene expression. Hierarchical clustering was used to observe activated or inhibited variation trends in the KEGG pathways.

Graphical Abstract highlights the experimental design, including the datasets used (GSE181279, GSE63060, and GSE63061) and their respective contributions to identifying immune cell types and gene expression profiles in AD and normal controls. The cell clustering based on scRNA-seq data is depicted, illustrating the differences in immune cell populations between AD patients and controls. The identification of key DEGs and the co-expression network analysis are also represented, emphasizing the role of immune cells in AD progression.

The raw GSE181279 dataset was read using the Seurat R package, and 22,776 individual cells in AD and 14,074 individual cells in normal tissue were obtained (Supplementary Figure 1). This was followed by quality control leading to 21,791 high-quality individual cells and 13,877 individual cells from normal individuals identified (Supplementary Figure 2). Using graph-based tSNE, we identified 13 clusters of major immune cells containing 35,668 total cells. Subsequently, we annotated 13 cell clusters with marker genes for major cell types and found that they were annotated as B cells, T cells, monocytes, and natural killer (NK) cells (Figures 1A, C). The distribution of T cells was clear and loose. Thus, we subdivided the T cell population and identified five T cell subtypes: gamma-delta T cells, CD8+ T cells, Tem/Temra cytotoxic T cells, Tem/Trm cytotoxic T cells, and mucosal-associated invariant T (MAIT) cells (Figures 1B–D). Different cell types were enriched in patients with AD and NC. Monocytes, CD8+ T cells, Tem/Temra cytotoxic T cells, and Tem/Trm cytotoxic T cells were enriched in patients with AD, whereas MAIT and NK cells were enriched in NCs (Figure 2E). Cell type proportions were significantly different between AD and NC, with monocytes, Tem/Temra cytotoxic T cells, and Tem/Trm cytotoxic T cells predominantly present in the AD samples, and MAIT cells, NK cells, and B cells predominantly present in the NC samples (Figure 1F).

We sought to explore the communication networks between immune cells. The interactions between Tem/Temra cytotoxic T cells, Tem/Trm cytotoxic T cells, and CD8 cells were stronger, and CD8 cells were signal recipients (Figure 2A). Three cell branches were identified in the pseudotime analysis: gamma-delta T cells and CD8+ T cells concentrated at the end of branch 1, monocytes concentrated at the end of branch 2, and Tem/Temra cytotoxic T cells and Tem/Trm cytotoxic T cells bifurcating into branch 3 (Figure 2B). The gene expression of immune cells in these branches is shown in Figure 2C. These results revealed the potential transcriptional heterogeneity between cell types.

To explore the co-expression networks of genes in immune cell types, we performed WGCNA. The power parameter range of 1–30 was filtered the power of β = 8 (scale-free R2 = 9) was used as the optimal screening soft threshold to construct a scale-free network and obtain five co-expression modules (Figures 3A, B). The modules spanned multiple cell populations by scoring the top 25 hub genes of each module, then mapping them to single cells (Figure 3C). We identified the top 10 hub genes in each module (Figure 3D). Among these modules, yellow (M3) and green (M5) modules showed significant positive correlations with CD8+ T cells and significant negative correlations with gamma-delta T cells in AD, while turquoise module (M1) showed significant positive correlations with monocytes and MAIT cells and significant negative correlations with CD8+ T cells, gamma-delta T cells, and B cells in AD (Figure 3E). Moreover, all immune cells were significantly different in the M1 group (Figure 3F).

We investigated the contribution of immune cells in AD. The results indicated that MAIT cells had the highest contribution to AD, followed by B cells and NK cells (Figure 4A). By comparing gene expression changes between AD samples and normal controls, we found that MAIT cells and monocytes were primarily upregulated in AD (Figure 4B). Enrichment analysis of DEGs in various cell types revealed that the DEGs were mainly enriched in growth hormone synthesis, secretion, and action (Figure 4C). In the quantitative analysis, we found that Butanoate metabolism was significantly activated in MAIT cells (Figure 4D).

In addition, to explore the disease mechanism in patients with AD from a molecular perspective, we analyzed the genes in immune cells identified by scRNA-seq for differential expression. We identified 1683 DEGs in B cells, 951 DEGs in CD8+ T cells, 1407 DEGs in gamma-delta T cells, 250 DEGs in MAIT cells, 1615 DEGs in monocytes, 644 DEGs in NK cells, 1595 DEGs in Tem/Temra cytotoxic T cells, and 1582 DEGs in Tem/Trm cytotoxic T cells (Figure 5A). Enrichment analysis revealed that these DEGs were enriched for regulation of T cell proliferation, B cell activation, and inflammatory responses (Figure 5B). In KEGG pathway analysis, we found that the DEGs were mainly enriched for ribosomes, Parkinson’s disease, and AD (Figure 5C). Of these, MAIT cells are involved in the fewest signaling pathways. However, all are implicated in the nervous system.

Notably, we identified 110 DEGs from the GSE63060 dataset (Figure 6A) and 429 DEGs from the GSE63061 dataset (Figure 6B). Analyzing the upregulated and downregulated DEGs separately, we found eight DEGs that were upregulated and 89 DEGs that were downregulated in both GSE63060 and GSE63061. Both were considered common DEGs (Figure 6C). Enrichment analysis showed that common DEGs were involved in SRP-dependent co-translational proteins targeting the membrane, oxidative phosphorylation, and ribonucleoprotein complex biogenesis (Figure 6D).

We further extracted candidate genes (CETN2, COX17, MRPL51, NDUFA1, NDUFS5, RPA3, RPL36AL, RPS25, SHFM1, and UFC1) with high AUC values in both GSE63060 and GSE63061 (Figure 7A). RPS25 and RPL36AL were highly expressed in the eight immune cell types identified by scRNA-seq (Figures 7B, C). Interestingly, all candidate genes exhibited decreased expression in AD patients compared to normal controls in both the GSE63060 and GSE63061 datasets (Figure 7D).

In addition, we analyzed the KEGG pathways in AD of GSE63060 (Figure 8A) and GSE63061 (Figure 8B) using GSEA. The results showed that in AD patients, osteoclast differentiation, lysosomes, TNF signaling pathway, JAK-STAT signaling pathway, growth hormone synthesis, secretion, and action were activated; while ribosomes, oxidative phosphorylation, coronavirus disease-COVID-19, Parkinson’s disease, and thermogenesis were inhibited. We found that ribosomes (hsa03010) and coronavirus disease-COVID-19 (hsa05171) were mainly enriched in the eight immune cell types identified by scRNA-seq (Figure 8C).

Temporal expression analysis was performed using STEM software of the differentially expressed genes between AD and MCI and between MCI and controls to further explore the expression patterns of genes in the progression of AD. We found that the expression trends of the 13 genes in the course of the normal control to MCI and AD were consistent in the GSE63060 (Figure 9A) and GSE63061 (Figure 9B) datasets. Among these 13 genes, UFC1 was identified as one of the candidate genes.

Cluster heatmap analysis of common DEGs in GSE63060 (Figure 10A) and GSE63061 (Figure 10B) revealed that the genes upregulated in AD progression were mainly involved in the regulation of neuronal apoptotic processes and the actin cytoskeleton. In contrast, the genes downregulated in AD progression were mainly involved in ribosome and oxidative phosphorylation.