Microglia Mediate the Occurrence and Development of Alzheimer’s Disease Through Ligand-Receptor Axis Communication

Alzheimer’s disease (AD) is a common neurodegenerative disease. Its onset is insidious and its progression is slow, making diagnosis difficult. In addition, its underlying molecular and cellular mechanisms remain unclear. In this study, clustering analysis was performed on single-cell RNA sequencing (scRNA-seq) data from the prefrontal cortex of 48 AD patients. Each sample module was identified to be a specific AD cell type, eight main brain cell types were identified, and the dysfunctional evolution of each cell type was further explored by pseudo-time analysis. Correlation analysis was then used to explore the relationship between AD cell types and pathological characteristics. In particular, intercellular communication between neurons and glial cells in AD patients was investigated by cell communication analysis. In patients, neuronal cells and glial cells significantly correlated with pathological features, and glial cells appear to play a key role in the development of AD through ligand-receptor axis communication. Marker genes involved in communication between these two cell types were identified using five types of modeling: logistic regression, multivariate logistic regression, least absolute shrinkage and selection operator (LASSO) and support vector machine (SVM). LASSO modeling identified CXCR4, EGFR, MAP4K4, and IGF1R as key genes in this communication. Our results support the idea that microglia play a role in the occurrence and development of AD through ligand-receptor axis communication. In particular, our analyses identify CXCR4, EGFR, MAP4K4, and IGF1R as potential biomarkers and therapeutic targets in AD.


INTRODUCTION
Alzheimer's disease (AD) is the most common cause of dementia and the most frequent type of progressive neurodegenerative disease. The major clinical characteristics of AD involve cognitive and behavioral dysfunction, and its major pathological feature is unregulated accumulation of amyloid β (Aβ) and intraneuronal neurofibrillary tangles (NFTs) in neurons (Stancu et al., 2014;Hogh, 2017;Zhu et al., 2017). At present, about 47 million people worldwide suffer from dementia, and the number is projected to increase to 131 million by 2,050 (Tiwari et al., 2019).
Due to its insidious onset and slow progression, diagnosis of AD can be difficult (Hogh, 2017). So far, five drugs have been approved by the US Food and Drug Administration (FDA) for the treatment of AD (Cummings et al., 2014). However, these drugs only alleviate disease progression or symptoms, without curing the disease. Thus, there is a critical need to explore the underlying pathogenesis of AD in order to develop therapeutic and even preventive interventions.
AD involves alterations in the interactions between neurons and glial cells (De Strooper and Karran, 2016). Brain damage caused by neuronal death may lead to cognitive decline, memory loss, learning disability and emotional changes (Vasic et al., 2019). Improving the viability of endogenous neuronal progenitor cells (NPCs) can prevent neural injury and degenerative disease (Zou et al., 2016). Activation of the MAPK signaling pathway in neurons may contribute to degenerative neuronal death (Kheiri et al., 2018).
Oligodendrocyte progenitor cells (OPCs), which express the proteoglycan NG2, are the main proliferating cells in the adult central nervous system. OPCs differentiate during postnatal development into myelin-forming oligodendrocytes (de Faria et al., 2019), which wrap around axons in the central nervous system, forming an insulating myelin sheath that assists in signal transmission as well as maintains and protects neuronal function (Nasrabady et al., 2018). During AD development, oligodendrocytes exhibit morphological alterations, deterioration of myelin integrity and axonal destruction (Cai and Xiao, 2016). Apoptosis among oligodendrocytes renders neurons vulnerable to injury and loss.
Microglia, the resident immune cells in the central nervous system, account for approximately 10% of all cells in that system, and they play a subtle and complex role in the development of AD (Lawson et al., 1990;Cameron and Landreth, 2010). Increasing evidence indicates that microglia act as a "double-edged sword, " both helping and harming neurons (Hanisch and Kettenmann, 2007;Sierra et al., 2013). On the one hand, microglia act as scavengers in the brain, clearing substances such as Aβ and thereby preventing the formation of amyloid plaques in the brain. On the other hand, prolonged microglial activation leads to the release of pro-inflammatory cytokines, which trigger a pro-inflammatory cascade that damages and destroys neurons (Perry et al., 2010;Jiang et al., 2012). This cascade involves extensive, complex signaling via chemokines, cytokines and Tolllike receptors (Zhou et al., 2021).
AD development may involve dysregulation of the normally complex yet balanced communication among different cell types in the central nervous system (Vainchtein and Molofsky, 2020). For example, astrocytes and microglia express many cytokines, chemokines and signaling molecules that can activate or destroy neighboring cells, and dysregulation of their production can promote the development of AD (Orre et al., 2014;Ceyzeriat et al., 2018).
Further elucidation of which cell types communicate with one another and of how they communicate may help us understand how AD occurs and progresses. This may then guide the development of novel diagnostic and therapeutic strategies. Toward this end, the present study analyzed singlecell RNA sequencing (scRNA-seq) data from AD brains. Different from previous studies, we provide novel answers and molecular evidences for the cellular communication at the single-cell transcriptome level to mediate the occurrence and development of AD. In particular, microglia play a leading role in the intercellular communication network of AD to regulate the biological signal cascade of target cells through ligand-receptor axis contributing to AD development.

Data Sources
Postmortem human brain samples were obtained from 48 participants in the Religious Order Study (ROS) or the Rush Memory and Aging Project (MAP), a longitudinal cohort of aging and dementia that included clinical data, detailed post-mortem evaluations, and "omics" tissue profiling (Bennett et al., 2018). Half of these individuals showed substantial Aβ burden and other pathological features of AD, while the other half showed no or extremely low Aβ burden or other pathological features of AD. The prefrontal cortex tissue from these 48 people underwent single-cell RNA sequencing (scRNA-seq), and the data came from Hansruedi Mathys' group of professors (Mathys et al., 2019).
To identify genes in neuron-glia communication, the following AD datasets from the Gene Expression Omnibus (GEO) database (Barrett et al., 2013) were used as a training set: GSE16759, GSE18309, GSE28146, GSE4757, GSE48350, GSE5281, GSE84422, and GSE9770 on the GPL570 platform. The datasets GSE33000 and GSE44772 on the GPL4372 platform served as a validation set.

Cell Clustering and Differentially Expressed Genes
The Seurat Package 1 was used to cluster cells based on scRNA-seq data, and the clusters were visualized using the UMAP package (Becht et al., 2018). In order to define marker genes of clusters, genes differentially expressed between AD patients showing clear pathology and individuals showing no or minimal pathology were determined in each cell type using the "FindAllMarkers" function with the settings min.pct = 0.25 and logfc.threshold = 0.25. Differences associated with P < 0.05 were considered significant.

Pseudo-Time Analysis
We employed the Monocle 2 using variable genes selected by Seurat as the input to determine the evolutionary state of different cell types (Qiu et al., 2017). The gene-cell matrix in the scale of unique molecular identifier (UMI) counts was provided as input to Monocle, and then the "newCellDataSet" function was used to create an object with the parameter expressionFamily = negbinomial.size. After dimension reduction and cell ordering, cell trajectory was inferred using default parameters of Monocle.

Correlation Between Pathological Features and Differentially Expressed Genes, and Functional Enrichment Analysis
Correlations between pathological features and DEGs were calculated, and those DEGs that were associated with pathological features were screened for functional enrichment based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment using the clusterprofiler package . Results associated with P < 0.05 were considered significant.

Intercellular Communication Analysis
Intercellular communication was performed using the iTALK routine within the R package. The receptor ligand interaction was determined using the STRING database of protein-protein interactions (Szklarczyk et al., 2017).

Construction of a Diagnostic Model and Analysis of Receiver Operating Characteristic Curves
A logistic regression model was constructed based on the optimal features of DEGs that correlated with AD pathology. Subsequent multivariate logistic regression was performed using genes that were significant by logistic regression. In addition, least absolute shrinkage and selection operator (LASSO) logistic regression was performed using the "glmnet" routine in R (project.org/package=glmnet). The LASSO model that was able to group samples based on the smallest number of DEGs was selected, and the output was analyzed using the "plot" function. The "tune" function in R (e1071) was used to optimize gamma values and generational values via 10-fold cross validation from an initial range from 10 −20 to 10 20 (Jiang et al., 2020). Finally, a support vector machine (SVM) model was also constructed.
To assess the diagnostic efficacy of the three models, we performed ROC curve analysis using the pROC package (Robin et al., 2011), allowing calculation of the area under the ROC curve (AUC).

Molecular Docking of Ligands Into Receptors
PDB files of key candidate ligands and receptors were downloaded from the Protein Database 2 (Burley et al., 2017), 2 www.rcsb.org and molecular docking studies were performed using Hex8.0.0 software (Macindoe et al., 2010). Docking models were visualized using Pymol software (Mooers, 2020).

Single-Cell Transcriptomic Landscape in Brain Tissue From Alzheimer's Disease Patients and Controls
In this study, scRNA-seq data were analyzed from the prefrontal cortex of individuals with clear AD pathology or with no or minimal pathology. The experimental design flow chart is shown in Figure 1A. These cells were divided into eight types ( Figure 1B): excitatory neurons (Ex), oligodendrocytes (Oli), inhibitory neurons (In), microglia (Mic), OPCs, astrocytes (Ast), endothelial cells (End), and pericytes (Per). The compositions of different cell types in each sample are shown in Figure 1C. We found that in individuals showing no or little AD pathology, IL1RAPL1, PCDH9, CNTNAP2, LSMP, and RORA could serve as cellular markers in most cell types ( Figure 1D). Figure 1E shows the marker gene expression patterns in each cell type that were shared between individuals with AD or with no or little AD pathology.
To further explore cell specificity between individuals with or without strong AD pathology, potential marker genes of AD pathology were extracted, genes overlapping with common marker genes were removed, and those specifically expressed in individuals with strong AD pathology were identified ( Figure 1F).

Evolution of Expression Dysregulation During Alzheimer's Disease Progression
We identified nine clusters in excitatory neurons, oligodendrocytes and microglia, and we also observed the abundance of these cell subsets in AD and control groups (Figure 2A). Interestingly, cluster 5 abundance was highest in microglia. Then a trajectory was constructed for each cell type using Monocle. We found that the various cell types followed different trajectories. Cell changes over pseudo-time were shown in Figure 2B, and different clusters were enriched along different paths (Figure 2C). We explored DEGs across pseudotime in different cell types (Figure 2D and Supplementary  Figure 1).

Pathological Process of Alzheimer's Disease at the Single-Cell Level
In order to understand the pathological changes during the process of AD, we analyzed the clinical data, pathological features and cell types in all 48 samples. We found that compared with individuals showing weak or no AD pathology, those showing strong pathology showed significantly higher global AD pathology burden (gpath), overall amyloid level (amyloid) and neuritic plaque burden (plaq_n) (Figure 3A). In contrast, the two groups of individuals did not differ significantly in cognitive state at death (cogdx). Then we calculated the correlation between the abundance of each cell type in each sample and the four pathological characteristics (amyloid, cogdx, gpath, plaq_n). We found that microglia significantly correlated positively with amyloid and the gpath. OPCs correlated negatively with amyloid and positively with cogdx. Oligodendrocytes correlated positively with cogdx, while inhibitory neurons correlated negatively with gpath and plaq_n ( Figure 3B).
In order to achieve a broad description of the AD ecosystem, we identified DEGs in each cell type between the two groups of individuals (Figure 3C), as well as their correlations Variables included age_death, defined as age at death; education, years of education; msex, self-reported sexual orientation; "1," male sex; "0," female sex; gpath, global AD pathology burden; plaq_n, neuritic plaque burden; amyloid, overall amyloid level; and cogdx, cognitive state at death. (B) Bubble heatmaps demonstrating cell abundances in individuals with strong AD pathology (gpath, amyloid, plaq_n, Cogdx). (C) Manhattan map showing genes within each cell type that were differentially expressed between the AD-pathology and no-pathology groups. (D) Quadrant diagram showing the correlation between pathological characteristics and differentially expressed genes. The ordinate represents correlation coefficient, the abscissa represents log [fold change (FC)], and color represents cell type.
with pathological features (Figure 3D). Those DEGs showing a correlation with such features were analyzed further as described below.

Microglia Mediate Alzheimer's Disease Occurrence and Development
Among the KEGG pathways enriched with DEGs in each cell type (Figure 4A), the pathways MAPK signaling, AD, necroptosis, ferroptosis and Notch signaling attracted our attention, because they have been related to AD progression. We analyzed intercellular communication in the two groups of individuals in order to explore differentially expressed ligands and receptors in different cell types. This allowed us to construct a comprehensive regulatory network that may contribute to AD (Figure 4B). The iTALK analysis identified significant differences in the abundance of ligand-receptor pairs. Microglia targeted excitatory neurons through the GPC3-IGF1R axis, Microglia are postulated to communicate with neurons through the GPC3-IGF1R axis, activating the MAPK signaling pathway, promoting neuronal apoptosis and thus mediating the occurrence and development of AD. Microglia may target oligodendrocytes through the APP-NGFR and CXCL12-CXCR4 axes, activating AD signaling pathways and thereby promoting disease occurrence and development. Microglia may target oligodendrocyte progenitor cells through the FGL1-EGFR axis to activate the necrotic apoptotic pathway, triggering the death of oligodendrocytes and rendering neurons vulnerable to injury. This, in turn, contributes to AD occurrence and development.
FIGURE 5 | Hub genes as biomarkers for the diagnosis of Alzheimer's disease. (A) Forest map showing single-factor logistic diagnostic efficacy of potential mechanistic genes. CXCR4, EGFR, MAP4K4, and IGF1R emerged as disease risk factors. (B) AUCs for multivariate logistic, LASSO, random forest, and SVM modeling. All AUCs were more than 0.5 and therefore considered significant. (C) AUCs for the four models based on GSE33000 and GSE44772 datasets on the GPL4372 platform. AUCs were more than 0.5. (D) Box diagram showing the expression of core genes (model genes) in individuals with strong AD pathology or no/weak pathology.
oligodendrocytes through the APP-NGFR and CXCL12-CXCR4 axes and OPCs through the FGL1-EGFR axis. Molecular docking studies predicted negative docking energies for GPC3 and IGF1R, APP and NGFR, CXCL12 and CXCR4, as well as FGL1 and EGFR, which suggests that such ligand-receptor complexes can form in vivo (Figure 4C).
Based on these results, we postulate that microglia communicate directly with neurons through the GPC3-IGF1R axis, which then activates the MAPK signaling pathway to promote neuronal apoptosis and thereby mediate the development and progression of AD. In contrast, microglia target oligodendrocytes through the APP-NGFR and CXCL12-CXCR4 axes to activate the AD pathway and thereby promote the development and progression of the disease. In addition, microglia target OPCs through the FGL1-EGFR axis to activate necrotic apoptosis, which reduces the abundance and activity of oligodendrocytes and thereby renders neurons vulnerable to injury, contributing to AD development ( Figure 4D).

Hub Genes as Biomarkers for Alzheimer's Disease Diagnosis
The hub genes for intercellular communication involved in the progression of AD development were extracted for univariate logistic analysis, and 16 genes with a significant P-value were obtained ( Figure 5A). These 16 genes were then used to conduct multi-factor logistics modeling, LASSO modeling, random forest modeling, and SVM modeling. These models gave respective AUCs of 0.726, 0.719, 1.000, and 0.970 for discriminating individuals with severe or minimal AD pathology ( Figure 5B).
Then we used GSE33000 and GSE44772 datasets of the GPL4372 platform to verify the four models ( Figure 5C). AUCs were 0.699 for logistic modeling, 0.705 for LASSO modeling, 0.61 for random forest modeling, and 0.548 for SVM modeling. Since the LASSO model showed AUC > 0.7 for both the training and validation sets, it was selected as the optimal model. Nine characteristic genes obtained from the LASSO modeling were extracted, and CXCR4, EGFR, MAP4K4, and IGF1R were found to be expressed at higher levels in individuals with AD than in the control group ( Figure 5D). Thus, these four genes may be useful in the diagnosis of AD.

DISCUSSION
AD is one of the most common and devastating agingrelated neurodegenerative disorders. Previous studies on AD have focused mainly on data from single cells, but few studies have examined intercellular communication affecting biological processes within cells (Hansen et al., 2018;Wang and Colonna, 2019). For example, studies on AD pathophysiology have focused on gene transcription (Mantzavinos and Alexiou, 2017;Lashley et al., 2018), showing that transcription factors regulate disease processes in AD (Zou C. et al., 2019). Other studies have shown that the microRNA miR-34a plays an important role in AD pathology, mainly by inhibiting the amyloidogenic processing of APP (Jian et al., 2017). Here we report single-cell transcriptomics for AD-related pathology based on a cohort of 48 individuals showing strong disease features or no or minimal pathology. Of eight major cell types in the brain, we identified microglia, OPCs, oligodendrocytes, and excitatory neurons as playing crucial roles in the progression of AD. Neuronal and glial cells correlated with several pathognomonic features of AD: global AD pathology burden, neuritic plaque burden, overall amyloid level, and cognitive state at death. Most hallmark genes were expressed specifically in individuals with strong AD pathology. Our results are consistent with the idea that AD onset and progression involve alterations not in a single gene but in intercellular connections.
The ligand-receptor axis is a fundamental means of communication between cells. Ligand binding to specific cell surface receptors can initiate intracellular signaling cascades that lead to various cellular responses (Fafian-Labora and O'Loghlen, 2020). The receptor-ligand axis can mediate the development of diseases (Allard et al., 2012). For instance, many ligand-receptor signaling patterns have been described between hepatocytes, illustrating the crucial role of intrahepatic crosstalk in tissue homeostasis and injury response (Krenkel and Tacke, 2017). Our data indicate extensive association of microglia with neurons, oligodendrocytes, and OPCs in AD. This association may involve the p38 MAPK pathway, since such signaling has been shown to mediate neuroinflammation triggered by microglia and astrocytes, and it has also been associated with autophagy-type AD . Studies have shown that the GLP-1 ligand-receptor axis is involved not only in AD but also in Parkinson's and Huntington's diseases (Janssens et al., 2014). The present study extends this literature by clarifying additional ligand-receptor axes involved in AD.
In the present work, the marker genes CXCR4, EGFR, MAP4K4, and IGF1R were highly expressed in the prefrontal cortex of individuals with strong AD pathology. CXCR4 has previously been shown to regulate neuronal firing and neuron/glia communication (Bezzi et al., 2001). In recent studies, EGFR has been validated as a possible target in the treatment of AD (Tsuji et al., 2021). Development of MAP4 kinase inhibitors can act as motor neuron protectors (Bos et al., 2019). Long-term inhibition of IGF1R signaling can attenuate AD progression and promote neuroprotection (George et al., 2017). These studies and our results argue strongly that these four genes may serve as biomarkers for AD diagnosis.
Our results should be interpreted with caution in light of some limitations. Overfitting is possible, given the strong performance of the model on the training set (AUC = 0.719), yet the model worked well with the validation set (AUC = 0.705), suggesting robustness and repeatability. Nevertheless, our results were obtained entirely through bioinformatic analysis, so they should be verified and extended in functional biochemical studies.

CONCLUSION
Our study suggests that microglia may mediate the occurrence and development of AD through ligand-receptor axis communication. We also identified CXCR4, EGFR, MAP4K4, and IGF1R as potential biomarkers of the disease and as candidate therapeutic targets.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
CJ, LW, RM, YL, and DZ conceived and designed the study. LW, RL, LC, LL, CZ, and YM collected the data. All authors analyzed the data, prepared figures and tables, wrote the manuscript, reviewed the manuscript, and approved its submission.