Single-Cell RNA-Seq Dissects the Heterogeneity of Human Microglial Related to Alzheimer’s Disease Based on Gene Regulatory Networks

doi:10.21203/rs.3.rs-1411661/v1

PDF

Research Article

Single-Cell RNA-Seq Dissects the Heterogeneity of Human Microglial Related to Alzheimer’s Disease Based on Gene Regulatory Networks

https://doi.org/10.21203/rs.3.rs-1411661/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted 11 Mar, 2022

You are reading this latest preprint version

Microglia, well recognized as immune effector cells inherent in the central nervous system, plays a fatal role in brain development and neurodegenerative diseases; however, the extent of human microglia heterogeneity and complex molecular functions are still poorly understood. Currently, most approaches to cell heterogeneity determine the function of cell subsets solely through high/low gene expression levels, indirectly ignoring the linkages between genes. Here, from the point of view of inter-gene association, we separately constructed gene regulatory networks of each microglia subtype based on single-cell RNA sequencing (scRNA-seq) and then identified significant genes in each network using tools from graph theory. Surprisingly, disease enrichment analysis of these genes revealed that two microglia subsets were closely related to Alzheimer's disease (AD), in which one subset was ulteriorly discovered involved in the antibiotic synthesis, and the other associated with Aβ plaque phagocytosis not only promoted inflammatory responses by NF-kappa B signaling pathway but also participated in tau protein binding and beta-amyloid binding, suggesting that this microglia subtype might act as tau protein phagocytosis as well. In a word, the approach that gene regulatory network combined with graph theory to identify the function of microglia subtypes makes up for the deficiency of differential expression and provides new insight into the cell type therapies related to AD.

Alzheimer's disease

microglia

Single-cell RNA-seq

gene regulatory network

graph theory

Parenchymal immune cells of the central nervous system not only make an essential contribution in regulating brain development and repairing damaged neurons but also participate in many neurological diseases [1]. Recently, microglia cells have been found to be activated or maladjusted in disease states, resulting in disease severity [2, 3]. Activated microglia could release numerous different cytokines and chemokines after the anti-inflammatory and pro-inflammatory signals imbalances, which makes continued in a state of immune response of the brain, leading to brain neuron loss and the emergence of diverse neurological disorders [4–7]. Therefore, the study of microglia subtypes and their internal molecular functions is of great significance to discover the pathogenesis and delay of these neurological diseases. Currently, most microglia subtypes have been obtained by conventional differential expression analysis [8]. However, gene expression is strictly regulated by transcription factors and cofactors and significantly up-regulated as well as down-regulated genes cannot fully reveal key factors in the phenotype. Single-cell network biology as a feasible scheme could further understand the potential regulator of cellular heterogeneity in human disease [9–10].

Recognizing that single-cell data requires new approaches to GRN modeling, a number of methods have recently been developed, such as SCENIC, PIDC, ACTION, SINCERA and so on. But SCENIC is limited to TF-based relationships [11], PIDC is difficult to scale to sparse single-cell data [12–13], ACTION requires TFs and their targets and only constructs TF-driven networks [14], and SINCERA also considers only TF and its targets and requires gene /TFs expression in most cells at the same time [15]. Only by comprehensively considering the regulatory and co-expression relationships of various genes, TF, etc., can we better identify microglia subtypes and their functions. Therefore, Iacono et al. 's approach, which considers gene expression is strictly regulated by TF, cofactors and signal molecules, is the most suitable model for our microglia subgroup network [16]. In addition, this network model adopts self-developed Z-scores to convert counts and Pearson correlation coefficient to capture the correlation between genes, which greatly improved the size and number of correlations. And the use of noisy probabilistic models also largely eliminates the effects of drop-out events and other technical artifacts. Nevertheless, the gene regulatory network obtained in this way is very dense. In order to identify significant nodes in the network, we introduce four kinds of centralities in graph theory, namely Degree, Betweenness, Closeness, and Pagerank. Each centrality could be regarded as a measure of the importance of a gene in the network and is obtained by a different algorithm. By screening out all the genes with high centralities, we could ulteriorly find the function and role of microglia subsets through enrichment analysis of these genes. Here, the preponderance of graph theory is not only to mine important elements with a regulatory network but also to expand other centralities, which greatly improves the reliability of the network model method.

In this study, we analyzed single-cell transcriptomes of 71,070 human brain parenchymal immune cells across diverse neurological disorders, and these cells were all purified by a validated experimental pipeline [17]. Here, the Harmony algorithm was selected as the batch processing method for single-cell data, and the latest KNN-Louvain iterative algorithm was used to cluster the processed samples, which facilitated successful identifications of functionally diversified rare microglial populations associated with neurodegeneration, especially AD [18–19]. Inspired by the interconnectedness of genes, we combined knowledge of gene regulatory networks and graph theory to analyze the potential function of our microglia subtypes. Inherently, single-cell transcriptome data are highly sparse and extremely noisy; thus, we took the Z-scores obtained by iterative differential expression to evaluate the correlation between genes so as to construct the gene regulatory network of subpopulations. After that, we assessed the importance of genes according to the size of the different centralities of each node in the network and performed enrichment analysis based on the genes with high centrality ranking to explore the function of microglia subtypes. Finally, two microglia subtypes (clusters 5 and 16) with different inflammatory functions were found to be associated with AD. Cluster 16 was closely associated with pro-inflammatory responses, while cluster 5 might play an anti-inflammatory part complying with the brain immune environment. And so, this whole set of ideas gets rid of the constraints of differential expression and provides a new direction for the identification of rare microglia subtypes.

In this section, we described in detail the whole experimental design workflow, which included batch processing, cell clustering, construction of gene regulatory networks, and means of seeking out momentous genes in the networks.

Overall experimental flow design

The main goal of this research is to excavate the models of a different microglia population structure and explore the factors that play a key role in the unknown microglia cell types and ultimately infer the hypothesis of their association with diseases. Figure 1 describes the entire experimental procedure and designing scheme in detail, which included datasets of single-cell RNA sequencing (scRNA-seq) preprocessing (Figure. 1a), subpopulation-based single-cell regulatory network construction, and identification of momentous genes through gene centralities (Figure. 1b).

The human brain parenchymal immune cells across various neurological disorders were all derived from the dorsolateral prefrontal cortex and were contributed from 16 individuals of both sexes (Fig. 1a). The gene expression data from these cells were then executed to Harmony algorithm for eliminating the impact of batches on clustering after quality control, and from results in the graph, the cells from various subjects perfectly blend with others hence no bias produced by batch effects. After that, we made use of the iterative KNN-Louvain approach to cluster the samples, and the population structure consisted of 18 cell types was subsequently obtained by adjusting the proper resolution, of which nine cell subpopulations were subsets of microglia (Fig. 1b). To further explore the molecular heterogeneity of microglia clusters, we constructed gene regulatory networks of microglia subpopulations respectively and calculated the centrality of each node in the network using the concept of graph theory. Meanwhile, we evaluated the importance of each gene according to the centrality size ranking as well. Adds: AD&LBD Alzheimer disease and Lewy body disease, AD Alzheimer disease, GLI gliosis, PD&MDD Parkinson disease, and major depressive disorder, PD&OB&DM Parkinson disease and obesity and dementia, BD&MDD Bipolar Disorder and major depressive disorder, SP&HP schizophrenia and hypoxia, TBI traumatic brain injury, VD vascular dementia.

Batch processing with Harmony algorithm

Considering a variety of experimental and potential biological factors, we finally chose the Harmony algorithm with superior performance for batch processing of vast samples obtained [18]. Specifically, Harmony first embedded the single-cell transcriptome profile into low-dimensional space through principal component analysis (PCA) and then applied an iterative process to remove the effects specific to the dataset. The iterative process mainly consisted of two steps: clustering and calibration of each cell. One step was to calculate the linear correction factor specific to cell type or cell state by self-developed soft K-means clustering algorithm, and the model was

$$\underset{R,Y}{\text{min}}\sum _{i,k}{R}_{ki}2\left(1-{Y}_{k}^{T}{Z}_{i}\right)+\sigma {R}_{ki}\text{log}{R}_{ki}+\sigma \theta {R}_{ki}\text{log}\left(\frac{{O}_{ki}}{{E}_{ki}}\right){\varnothing }_{i}$$

$$s.t.{\forall }_{i}{\forall }_{k}{R}_{ki}>0,{\forall }_{i}\sum _{k=1}^{K}{R}_{ki}=1$$

The formula showed the distance of the cells to their assigned centroids concretely, in which Z represented the eigenspace of some data and Y represented centroids obtained by clustering. Besides, ${R}_{ki}$ denoted the relationship between cell i and cluster k, and $l\text{o}\text{g}\left(\frac{{O}_{ki}}{{E}_{ki}}\right){\varnothing }_{i}$ meant adding an entropy regularization term to R, weighted by a hyperparameter σ. For each batch variable, this algorithm introduced a new parameter θ to measure the dependency between batch membership and cluster allocation. The other step was to eliminate batch-specific variations from each cell cluster according to the linear mixture model, and the model was

$${W}_{k}={{(\varnothing }^{*}diag\left({R}_{k}\right){\varnothing }^{{*}^{T}}+\lambda I)}^{-1}{\varnothing }^{*}diag\left({R}_{k}\right){Z}^{T}$$

The equation above was mainly used to solve the collinearity issue of the model and ${\varnothing }^{*}$ was a design matrix containing intercept items that capture the independent batch variations in each cluster. Z represented the matrix of the original PCA embedding and $diag\left({R}_{k}\right)$ represented the diagonal matrix of cluster member terms in cell cluster k. Additionally, λ as ridge penalty hyperparameter would shrink ${W}_{k}$ gradually approaching 0.

Clustering based on an iterative KNN-Louvain approach

In order to cluster the pre-processed samples accurately, we combined batch processing and adopted a graph-based clustering method, which primarily constructed a KNN graph using the Euclidean distance in PCA space and refined the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). After that, the Louvain algorithm, as one of the most popular algorithms to optimize modularity, was applied to uncover the community structure by means of iterative grouping [19], and the algorithm model was

$$Q= \frac{1}{2m}*\sum _{ij}\left[{A}_{i,j}-\frac{{k}_{i}*{k}_{j}}{2m}\right]*\delta ({C}_{i},{C}_{j})$$

In this equation, m was the total number of edges involved in clustering and ${A}_{i,j}$ referred to the edge weights between nodes i and j. If nodes i and j were connected, then ${A}_{ij}=1$; otherwise, ${A}_{ij}=0$. In addition, ${k}_{i}$ represented the sum of all edge weights pointing to node i, and the same went for k. ${C}_{i}$ was the community where node i belonged and ${C}_{j}$ was the community where node j belonged. When i and j located in the same community, $\delta ({C}_{i}, {C}_{j}) = 1$; else, it was 0.

Construction of gene regulatory network

After clustering, for constructing the gene regulatory network of microglia subsets, we then utilized function subset () in Seurat to extract the gene expression matrix of microglia subsets, severally. On account of the noise and sparsity of single cell data, we determined to employ Z-scores to measure correlations between genes since these correlations might be obscured by drop-out events and other technical artifacts, and the model was

$$Z-score=\sqrt{{Z}_{num.model}^{2}+{Z}_{Wilcoxon}^{2}}$$

This algorithm compensated for the weaknesses of numerical model and Wilcoxon, and the synergy of these two tests broke through the limitation of sample size. ${Z}_{num.model}$ was calculated by the numerical model, homoplastically ${Z}_{wilcoxon}$ was calculated by Wilcoxon rank-sum test.

Identification of significant genes in networks using graph theory

For the sake of accurately pinpointing the key players in microglia subsets, we chose four types of centralities (Degree, Betweenness, Closeness and Pagerank) in graph theory and network analysis to evaluate the status and significance of genes in gene regulatory networks and these four metrics are all obtained by different calculation methods. And if the four centrality values of a certain gene in the network are all large, it suggests that this gene may be associated with many other genes and play a vital role in this subtype of cell. Additionally, Degree refers to the number of undirected edges a node has to other nodes is the most basic means to evaluate the significance in networks, which indicates the possibility of acquiring the information spreading in the network. Betweenness implies the ability to control the information is a centrality based on shortest paths as the interaction between two non-adjacent members depends on the other members of the network, especially those on the path between the two members. It is calculated by enumerating all the shortest paths of the network and quantifying the number of times each node falls into the shortest path and its model is

$${C}_{B}\left(v\right)=\sum _{s\ne v\ne t\in V}\frac{{\sigma }_{st}\left(v\right)}{{\sigma }_{st}}$$

In this equation, ${\sigma }_{st}\left(v\right)$ represents the number of shortest paths of $s\to$v through node v, and ${\sigma }_{st}$ refers to the number of shortest paths, simply Betweenness responses the importance of node v as ‘bridge’. Closeness is from the reciprocal of the sum of the minimal distances of a vertex to all other vertices, which indicates the velocity of getting the information from network. Pagerank centrality created by google results from a random walk of the network and the nodes with high Pagerank shows the prevalence of the gene in the network. Besides, the specific algorithm is as follows:

$${PR}_{\left({p}_{i}\right)}=\frac{1-d}{n}+d\sum _{{P}_{j}\in {M}_{\left(i\right)}}\frac{{PR}_{\left({p}_{j}\right)}}{{L}_{\left(j\right)}}$$

In this formular, ${{p}_{i,}p}_{j}$ point different nodes in regulatory network, ${M}_{\left(i\right)}$ refers the collection of all the nodes link to this gene and ${L}_{\left(j\right)}$ is the number of external links which means other nodes contact to ${p}_{j}$. Moreover, n represents the quantity of nodes and d (0 < d ≤ 1) is damping factor which represents the probability that the link will continue after reaching a node, and the default value is 0.85. All these centralities were specifically calculated with the package igraph 1.2.2 in the end.

Data Sources and preprocessing

The dataset, which mainly consisted of 33 donors with different psychiatric disorders, stages of age and genders in this study, were all from the Synapse database [20]: syn24168322. However, many individuals were not labeled with specific diagnoses, so 16 individuals with explicit diseases, including Alzheimer Disease, Bipolar Disorder, Gliosis, Parkinson Disease and so on were ultimately selected for further analysis, which included six females and ten males. Then, through 10X method in the IlluminaNovaseq6000 platform, the dataset composed of 71,070 parenchymal immune cells from the dorsolateral prefrontal cortex in the brain of selected individuals were obtained after quality control. Furthermore, this microarray data also eliminated the effects of brain region heterogeneity in this way. Afterward, the dataset was loaded into a purified preprocessing pipeline to eliminate dimensional effects, and this process was then accomplished by functions NormalizeData () and ScaleData () via the Seurat package. Since immune cells were obtained under different experimental conditions in the presence of a batch effect, imposing a negative impact on clustering, we considered the combination of PCA reduction and Harmony algorithm for batch processing of standardized data to facilitate subsequent cluster analysis [21–23]. Eventually, we implemented these vital steps through the functions RunPCA () and RunHarmony () in Seurat and Harmony packages and chose the first 20 pcs as input.

Identification of human brain parenchymal immune cells

After the dataset was partially optimized by pretreatment, we used the iterative KNN-Louvain clustering algorithm for unsupervised clustering of cells, and these critical processes were then performed with functions FindNeighbsors () and FindClusters () in Seurat v3, respectively. Meanwhile, the resolution parameter of clustering was set as 0.6. To identify human brain myeloid cell subsets at the genetic level, non-parametric Wilcoxon rank sum test, a reliable and mainstream differential expression test, was ultimately identified to look for features of differential expression between clusters of cells and numerous known signature genes which have been demonstrated these genes could differentiate various cell types in broad categories were higher expression levels (Fig. 2). In addition, to visualize the clustered cells in a two-dimensional plane, we performed the tSNE (t-Distributed Stochastic Neighbor Embedding) nonlinear dimensionality reduction method to display the clustering results, and this step was achieved by function RunTSNE ().

Of the 18 cell clusters obtained by the iterative KNN-Louvain clustering algorithm, 10 clusters (clusters 1,2,4,5,6,9,13,14,16,18) were preliminarily inferred as myeloid cells because of a high level of myeloid marker (AIF1) expression and higher cells proportion expressing the markers through the Wilcoxon rank sum test [24], yet the expression of microglial cell markers (C1QA, C1QB, and C1QC) in cluster 6 was not very high; thus this cluster probably represents monocytes. Whereafter, the surplus clusters were analogously characterized by high expression of unique marker genes from the CellMarker database, such as MT-ND3 (cluster 3), CD3E (cluster 7,15), NKG7 (cluster 8), SPARCL1 (clusters 10,12), CD79A (cluster 11), and TPSAB1 (cluster 17), likely representing monocytes, T cells, natural killer cells, astrocytes, B cells, and mast cells, respectively (Fig. 2a). TSNE plot composed of tens of thousands of heterochromatic dots delineates the population distribution of different microglia subsets and other subgroups, in which each dot is considered as a cell; simultaneously all cells are put colors on according to their classification relations (Fig. 2b). By projecting the expression of selected genes into a model of population structure (Fig. 2c), marker genes were found enriched so observably in certain clusters gene, just as C1QA was significantly aggregated and highly expressed in a microglial cell and gene CD3E was mainly expressed in cluster 7, which implying the feasibility of genetically annotating cell types. Adds: MG microglia, Mono monocyte, TC T cell, NK nature kill cell, AST astrocyte, BC B cell, MC mast cell.

Assessment of microglia heterogeneity within and among individuals

To explore correlations between microglia subsets and individuals, we calculated the proportion of each cell subpopulation in the same individual with the diseases and compared the variations of the percentage of the same microglia subsets across different individuals from the view of the distribution of clusters.

From a global perspective, the bar plot depicts exactly how different clusters of cells are distributed in a donor and the line chart shows how the percentage of the same microglia changes in different individuals with the diseases. Then, the first five cell clusters distinctly show a higher percentage of cells in all individuals with different diagnoses (Fig. 3a), which indicates these cells might be homeostasis microglia associated with the cellular state (fresh or frozen). Moreover, the proportion of three microglia subsets (clusters 9,16 and 18) was dramatically upregulated in some diagnosed individuals, cluster 16 occurs only in Alzheimer's disease and gliosis as well as cluster 18 in bipolar disorder and major depressive disorder, which implies that these cell subtypes could serve as a potential target for disease treatment (Fig. 3b, c, d).

Annotation of microglia clusters using reported features

To further understand the function and role of microglia subtypes, we then turned to annotate our microglia clusters obtained by the iterative KNN-Louvain clustering algorithm (clusters 1, 2, 4, 5, 6, 9, 13, 14, 16, and 18) using microglia types and corresponding genes that have appeared in the literature.

On account of these highly expressed genes in cluster 4 (Fig. 4a), we deduced that our human microglia subpopulation (cluster 4) represented “IEG+” microglia which was produced by the stimulus of environmental change since the genes had recently been described to be induced by tissue dissociation or FACS by van den Brink and colleagues, and most of the genes in the plot were immediate early genes (IEG) [25–26]. After that, in the same way, cluster 9 (Fig. 4b) was found to be enriched in amounts of interferon-related genes, which indicated that this subtype probably represented late-response state microglia [24, 27]. Furthermore, homeostasis microglia markers were highly expressed in the first four cell subsets (Fig. 4c), consistent with our previous assumptions, suggesting that these cells (clusters 1,2, and 5) mainly perform daily nervous system tasks [28]. Finally, cluster 16 was annotated as a subtype of microglia associated with Aβ plaque phagocytosis in regards to a high expression of HIF1A, which might regulate synaptosome phagocytosis in vitro [29]. For the remaining microglia subpopulations (clusters 13, 14, and 18), we subsequently used functional annotation to explore its existence significance.

Single-cell regulatory networks identify essential genes in microglia subsets

Since genes interact with each other, we turned to a network model to explore the underlying genes in microglia subsets (Fig. 5). However, the abnormal noise and sparsity of single-cell data, along with large amounts of dropout events, made it very difficult to construct single-cell regulatory networks. Therefore, we chose the self-developed ${Z}_{-score}$ as the standard to measure the correlation between genes to increase the number of inter-gene associations. In addition, gene regulation networks were extremely complex and node–dense; thus, we filtered nodes in the network according to the four centralities in graph theory. Concretely, this process was completed by the function compute.network () in bigSCale and the adaptive threshold was ultimately set to 0.85 (ρthresh = 0.85) to retain significant correlations. Here, cluster 1 was taken as the example, and the visualization results were as follows.

In this part, we mainly described the construction of the gene regulatory network and the screening criteria of characteristic genes. From Fig. 5a, it could be seen that many genes (TIMM17A, MRTO4, SLC11A1, SLC7A5, NME1, and SORL1) with higher degree centrality tended to accompany by higher betweenness centrality. Therefore, gene screening should not rely on a single indicator but a combination of four indicators to screen genes in the network. Taken cluster 1 as an example (Fig. 5b), we sorted the genes in the gene regulatory network of cluster 1 according to the size of their centralities and then intersected the top genes to select approximately 200 genes for subsequent enrichment analysis. For unknown microglia clusters (Cluster 13, 14, and 18), we did the same and selected about 200 genes in the end.

Exploration of associations between human microglia clusters and diseases

Given too few individuals involved and lack of evidence of an association between microglia subpopulation abundance and Alzheimer's disease, we determined to use the DOSE Bioconductor package to perform disease ontology (DO) enrichment analysis to excavate possible relationship between microglia subsets and pathology of AD based on all genes in the gene regulatory network of microglia subpopulations (clusters 1, 2, 4, 5, 6, 9, 13, 14, 16 and 18) [30]. Living up to expectations, the enrichment results of clusters 5 and 16 both appeared AD pathology, as the following picture showed. Nevertheless, the remaining microglia subtypes (clusters 1, 2, 4, 6, 9, 13, 14, and 18) were not found to be significantly associated with AD. Therefore, we would focus only on clusters 5 and 16 here.

From the DO enrichment results (Fig. 6a, b), we could see that clusters 5 and 16 were obviously associated with AD, especially cluster 16 with a higher significance. And interestingly, we discovered that both Alzheimer's disease and tauopathy occurred in two clusters simultaneously, implying that our microglia subpopulations might affect AD by regulating Tau proteins. Due to the high significance of AD and tauopathy in cluster 16 and a large number of enriched genes, we decided to further discuss the association between AD and Tauopathy through the similarity of the enriched gene sets. Then, by mapping the network and Venn diagram between genes and diseases in cluster 16 (Fig. 6c, d), the genes enriched in Alzheimer's disease and tauopathy showed striking similarities, which indicated that this subtype of microglia associated with Aβ plaque phagocytosis might also influence synaptic growth through tau protein. In addition, the enriched gene sets in cluster 5 also showed great overlap in AD and tauopathy, as shown in Supplementary Table 1.

Functional enrichment analysis based on DAVID database

Eventually, we used the DAVID database to perform functional enrichment analysis of the selected genes at the genetic level [31]. Specifically, for unknown microglia subpopulations (clusters 13, 14, and 18), we screened about 200 genes with high centralities (Degree, Betweenness, Closeness, and Pagerank) ranking in the gene regulatory network of each subgroup for KEGG enrichment analysis. As for cluster 5 and cluster 16, we decided to annotate the function according to the genes enriched to AD in the DO analysis results (134 genes in cluster16 and 22 genes in cluster5). Moreover, since cluster 16 was closely associated with AD, we further performed Gene Ontology (GO) enrichment analysis based on those 134 genes, and the results are shown in the figure below.

From the metabolic pathway enrichment analysis (Fig. 7a and Supplement Table 2), it could be seen that cluster 13 was involved in antigen processing and presentation, cluster 14 was related to the chemokine signaling pathway, and cluster 18 was significantly enriched in cancer transcriptional dysregulation. Indirectly, it could be inferred that cluster 14 had a pro-inflammatory effect, and cluster 18 was closely associated with the formation of cancer because much of the literature had shown an inverse relationship between neurological diseases and cancer [32–35].

Homoplastically, from KEGG enrichment results of cluster 5 and cluster 16 (Fig. 7b and Supplement Table 3), cluster 5 was related to the biosynthesis of antibiotics, and cluster 16 was related to NF-kappa B signaling pathway, which implied that cluster 5 was an anti-inflammatory microglia subtype as well as cluster 16 a proinflammatory microglia subtype in AD pathology.

Due to the dramatic correlation between cluster 16 and AD (Fig. 6b), we then turned to perform the GO enrichment analysis of genes associated with AD in cluster 16 to explore BP, CC, and MF in DAVID (Fig. 7c and Supplement table 4). And as it turns out, cluster 16, a subtype of microglia associated with Aβ plaque phagocytosis, was not only involved in the inflammatory and innate immune response but also associated with tau and beta proteins, which indicated that this microglia subpopulation was involved in various pathogenesis of AD and could be a pivotal target for the treatment of Alzheimer's disease [36].

Additionally, cluster 16 was present only in AD and gliosis (Fig. 3b), which implicitly illustrated that this microglia subtype might also play a similar role in gliosis.

As network modeling approaches are increasingly recognized in biology, in this study, we provide novel thinking, single-cell sequencing combines with gene regulatory networks to uncover the microglia subpopulations related to AD. Specifically, we first clustered the preprocessed samples through the KNN-Louvain iterative clustering algorithm and then constructed gene regulation of subtypes and identified significant genes in the network using relevant knowledge in graph theory. At the same time, the rank-sum test was used to screen for genes significantly expressed between cell clusters to facilitate the identification of known microglia subsets. Since genes didn't act independently of each other, our approach based on the network model largely compensated for the defects caused by differential expression analysis and facilitated the mining of transcription factors that played a dominant role in subtypes.

Fortunately, we eventually identified two microglia subtypes (clusters 5 and 16) closely associated with AD through DO enrichment analysis, and these microglia cells were all related to inflammation but had completely different functions. Cluster 5, an anti-inflammatory microglial subtype, potentially modulated the immune microenvironment of the brain through biosynthesis of antibiotics while cluster 16, a pro-inflammatory microglial subtype, not only closely associated with Aβ plaque phagocytosis but also involved in a variety of pathogenesis of AD, including NF-kappa B signaling pathway, tau protein binding, and beta-amyloid binding, which indicated that cluster 16 possibly participated in tau protein phagocytosis.

Inflammation is a central mechanism in Alzheimer's disease, and our methods further reveal the connection of both at the cellular level. More critically, our approach opened up the possibility of cell therapy for AD and unearthed a new target that played a pivotal role in Alzheimer's disease.

Availability of Data and Materials

The gene expression data of patients across diverse neurological disorders were downloaded from the Synapse database (https://www.synapse.org/#!Synapse:syn24168322)

Data Availability Statements

1 The datasets generated during and/or analyzed during the current study are available in the Synapse database.

2 The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

3 All data generated or analyzed during this study are included in this published article.

Ethics Approval and Consent to Participate

This article does not contain any studies with human participants or animals performed by any of the authors; therefore, the ethical approval and consent to participate are not applicable.

Consent for Publication

The authors have consented to publication of this article.

Competing Interests

The authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China (No. 61803257) and Natural Science Foundation of Shanghai (No. 18ZR1417200).

Authors' contributions

Conception and design of the research: Haidong Wei and WeiKong. Acquisition, analysis, and interpretation of data: Haidong Wei, Wei KongandShuaiqun Wang. Statistical analysis: Haidong Wei and WeiKong. Drafting the manuscript: Haidong Wei and WeiKong. Manuscript revision for important intellectual content: Haidong Wei, Wei Kong andShuaiqun Wang. All authors have read and approved the manuscript.

Acknowledgements

The authors would like to thankthe Synapse Database for open access to the data andthe support from the National Natural Science Foundation of China andNatural Science Foundation of Shanghai.

1 Prinz M, et al. Heterogeneity of CNS myeloid cells and their roles in neurodegeneration. Nat Neurosci. 2011 Sep 27;14(10):1227-35. doi: 10.1038,nn.2923.

2 Saijo K, Glass CK. Microglial cell origin and phenotypes in health and disease. Nat Rev Immunol. 2011 Oct 25;11(11):775-87. doi: 10.1038,nri3086.

3 Paolicelli RC, Bolasco G, Pagani F, Maggi L, Scianni M, Panzanelli P, Giustetto M, Ferreira TA, Guiducci E,et al. Synaptic pruning by microglia is necessary for normal brain development. Science. 2011 Sep 9;333(6048):1456-8. doi: 10.1126,science.1202529. Epub 2011 Jul 21.

4 Grammas P. Neurovascular dysfunction, inflammation and endothelial activation: implications for the pathogenesis of Alzheimer's disease. J Neuroinflammation. 2011 Mar 25;8:26. doi: 10.1186,1742-2094-8-26.

5 Meraz-Ríos MA, et al. Inflammatory process in Alzheimer's Disease. Front IntegrNeurosci. 2013 Aug 13;7:59. doi: 10.3389,fnint.2013.00059.

6 Heneka MT, et al. Neuroinflammation in Alzheimer's disease. Lancet Neurol. 2015 Apr;14(4):388-405. doi: 10.1016,S1474-4422(15)70016-5.

7 Heneka MT, et al. Neuroinflammation in Alzheimer's disease. Lancet Neurol. 2015 Apr;14(4):388-405. doi: 10.1016,S1474-4422(15)70016-5.

8 Stratoulias V, et al. Microglial subtypes: diversity within the microglial community. EMBO J. 2019 Sep 2;38(17):e101997. doi: 10.15252/embj.2019101997. Epub 2019 Aug 2. PMID: 31373067; PMCID: PMC6717890.

9 Blencowe M, et al. Network modeling of single-cell omics data: challenges, opportunities, and progresses. Emerg Top Life Sci. 2019 Aug;3(4):379-398. doi: 10.1042,etls20180176. Epub 2019 Jul 8.

10 Cha J, Lee I. Single-cell network biology for resolving cellular heterogeneity in human diseases. Exp Mol Med. 2020 Nov;52(11):1798-1808. doi: 10.1038,s12276-020-00528-0. Epub 2020 Nov 26.

11 Aibar S, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017 Nov;14(11):1083-1086. doi: 10.1038/nmeth.4463. Epub 2017 Oct 9. PMID: 28991892; PMCID: PMC5937676.

12 Filippi S. and Holmes C.C. (2017) A Bayesian nonparametric approach to testing for dependence between random variables. Bayesian Anal. 12, 919–938 10.1214/16-BA1027

13 Chan TE, et al. Gene Regulatory Network Inference from Single-Cell Data Using Multivariate Information Measures. Cell Syst. 2017 Sep 27;5(3):251-267.e3. doi: 10.1016/j.cels.2017.08.014. PMID: 28957658; PMCID: PMC5624513.

14 Mohammadi S., et al. (2018) A geometric approach to characterize the functional identity of single cells. Nat. Commun. 9, 1516 10.1038/s41467-018-03933-2

15 Guo M., et al. (2015) SINCERA: a pipeline for single-cell RNA-seq profiling analysis. PLoSComput. Biol. 11, e1004575 10.1371/journal.pcbi.1004575

16 Iacono G., Massoni-Badosa R. and Heyn H. (2019) Single-cell transcriptomics unveils gene regulatory network plasticity. Genome Biol. 20, 110 10.1186/s13059-019-1713-4

17 Olah M, et al. A transcriptomic atlas of aged human microglia. Nat Commun. 2018 Feb 7;9(1):539. doi: 10.1038/s41467-018-02926-5.

18 Korsunsky I, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019 Dec;16(12):1289-1296. doi: 10.1038/s41592-019-0619-0. Epub 2019 Nov 18.

19 Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep. 2019 Mar 26;9(1):5233. doi: 10.1038/s41598-019-41695-z.

20 Zhang W, et al. SynDB: a Synapse protein DataBase based on synapse ontology. Nucleic Acids Res. 2007 Jan;35(Database issue): D737-41. doi: 10.1093/nar/gkl876. Epub 2006 Nov 10.

21 Luecken MD, Theis FJ. Current best practices in single-cell RNA-seq analysis: a tutorial. Mol Syst Biol. 2019 Jun 19;15(6):e8746. doi: 10.15252/msb.20188746.

22 Tran HTN, et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biol. 2020 Jan 16;21(1):12. doi: 10.1186/s13059-019-1850-9.

23 Macosko EZ, et al. Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets. Cell. 2015 May 21;161(5):1202-1214. doi: 10.1016/j.cell.2015.05.002.

24 Olah M, et al. Single cell RNA sequencing of human microglia uncovers a subset associated with Alzheimer's disease. Nat Commun. 2020 Nov 30;11(1):6129. doi: 10.1038/s41467-020-19737-2.

25 van den Brink SC, et al. Single-cell sequencing reveals dissociation-induced gene expression in tissue subpopulations. Nat Methods. 2017 Sep 29;14(10):935-936. doi: 10.1038/nmeth.4437.

26 Li Q, et al. Developmental Heterogeneity of Microglia and Brain Myeloid Cells Revealed by Deep Single-Cell RNA Sequencing. Neuron. 2019 Jan 16;101(2):207-223.e10. doi: 10.1016/j.neuron.2018.12.006. Epub 2018 Dec 31.

27 Mathys H, et al. Temporal Tracking of Microglia Activation in Neurodegeneration at Single-Cell Resolution. Cell Rep. 2017 Oct 10;21(2):366-380. doi: 10.1016/j.celrep.2017.09.039.

28 Bennett ML, et al. New tools for studying microglia in the mouse and human CNS.ProcNatlAcadSci U S A. 2016 Mar 22;113(12): E1738-46. doi: 10.1073/pnas.1525528113. Epub 2016 Feb 16.

29 Grubman A, et al. Transcriptional signature in microglia associated with Aβ plaque phagocytosis. Nat Commun. 2021 May 21;12(1):3015. doi: 10.1038/s41467-021-23111-1.

30 Yu G, et al. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. 2015 Feb 15;31(4):608-9. doi: 10.1093/bioinformatics/btu684. Epub 2014 Oct 17.

31 Huang da W, et al. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44-57. doi: 10.1038/nprot.2008.211. PMID: 19131956.

32 Roe CM, et al. Alzheimer disease and cancer. Neurology. 2005;64:895–898. doi: 10.1212/01.WNL.0000152889.94785.51.

33 Driver JA, et al. Inverse association between cancer and Alzheimer’s disease: results from the Framingham Heart Study. BMJ. 2012;344:e1442. doi: 10.1136/bmj.e1442.

34 Realmuto S, et al. Tumor diagnosis preceding Alzheimer’s disease onset: is there a link between cancer and Alzheimer’s disease? J Alzheimers Dis. 2012;31:177–182.

35 Musicco M, et al. Inverse occurrence of cancer and Alzheimer disease: a population-based incidence study. Neurology. 2013;81:322–328. doi: 10.1212/WNL.0b013e31829c5ec1.

36 Cummings J, et al. Drug development in Alzheimer's disease: the path to 2025. Alzheimers Res Ther. 2016 Sep 20; 8:39. doi: 10.1186/s13195-016-0207-9.

PDF

Version 1

posted 11 Mar, 2022

You are reading this latest preprint version

Single-Cell RNA-Seq Dissects the Heterogeneity of Human Microglial Related to Alzheimer’s Disease Based on Gene Regulatory Networks

Status:

Version 1

Abstract

Figures

Introduction

Methods

Overall experimental flow design

Batch processing with Harmony algorithm

Clustering based on an iterative KNN-Louvain approach

Construction of gene regulatory network

Identification of significant genes in networks using graph theory

Results

Data Sources and preprocessing

Identification of human brain parenchymal immune cells

Assessment of microglia heterogeneity within and among individuals

Annotation of microglia clusters using reported features

Single-cell regulatory networks identify essential genes in microglia subsets

Discussion

Exploration of associations between human microglia clusters and diseases

Functional enrichment analysis based on DAVID database

Conclusion

Declarations

References

Supplementary Files

Status:

Version 1