Biomarker Identification for Gender Specificity of Alzheimer's Disease Based on the Glial Transcriptome Profiles

Summary

This study analyzed single-nuclei transcriptomes of thirty-three individuals with Alzheimer's disease (AD), revealing sex-specific DEGs in glial cells. Functional enrichment analysis highlighted synaptic, neural, and hormone-related pathways. Key genes, namely NLGN4Y and its regulators, were identified, and potential therapeutic candidates for gender-specific AD were proposed.

Abstract

Many sex-specific biomarkers have been recently revealed in Alzheimer's disease (AD); however, cerebral glial cells were rarely reported. This study analyzed 220,095 single-nuclei transcriptomes from the frontal cortex of thirty-three AD individuals in the GEO database. Sex-specific Differentially Expressed Genes (DEGs) were identified in glial cells, including 243 in astrocytes, 1,154 in microglia, and 572 in oligodendrocytes. Gene Ontology (GO) functional annotation analyses and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses revealed functional concentration in synaptic, neural, and hormone-related pathways. Protein-protein interaction network (PPI) identified MT3, CALM2, DLG2, KCND2, PAKACB, CAMK2D, and NLGN4Y in astrocytes, TREM2, FOS, APOE, APP, and NLGN4Y in microglia, and GRIN2A, ITPR2, GNAS, and NLGN4Y in oligodendrocytes as key genes. NLGN4Y was the only gene shared by the three glia and was identified as the biomarker for the gender specificity of AD. Gene-transcription factor (TF)-miRNA coregulatory network identified key regulators for NLGN4Y and its target TCMs. Ecklonia kurome Okam (Kunbu) and Herba Ephedrae (Mahuang) were identified, and the effects of the active ingredients on AD were displayed. Finally, enrichment analysis of Kunbu and Mahuang suggested that they might act as therapeutic candidates for gender specificity of AD.

Introduction

Alzheimer's disease (AD) is a global disease with high incidence, and it accounts for 60%-80% of dementia¹. Despite its high incidence, the mechanistic pathogenesis of AD is not clearly delineated, and there have been no effective therapeutics until now². The main pathologies in AD were identified as neuronal atrophy and the accumulation of pathological debris, mainly microtubule-associated protein Tau, and β-amyloid (Aβ)^3,4. The pathogenesis of AD is associated with abnormal autophagy, oxidative stress, mitochondrial dysfunction, inflammation, and energy metabolism disorder⁵. Prevalence surveys proved that two-thirds of AD patients were women⁶. Sex-specific differences in AD exist in the etiology, clinical manifestations, prevention, and treatment. Thus, revealing the biological mechanism that causes sex-specific differences in AD and targeting traditional Chinese medicine (TCM) can potentially provide a more comprehensive theoretical framework to understand the pathogenesis of AD, and to further guide accurate treatment strategy.

Neuroglial cells, especially microglia, astrocytes, and oligodendrocytes, potentially contribute to the pathogenesis of AD. In AD, microglia are activated and genetically altered, which contribute to inflammatory response, phagocytosis, and Aβ clearance^7,8; astrocyte is genetically altered, which affects synaptic activity, ion homeostasis, and energy and lipid metabolism⁹; oligodendrocyte is genetically altered with sex specificity, which contributes to neuronal loss, neurofibrillary tangles, and white matter lesions^10,11.

In this study, we employed single-nuclei RNA sequencing (snRNA-seq) as a superior technique. Compared to single-cell RNA sequencing (scRNA-seq), snRNA-seq offers advantages in terms of sample richness, cell type integrity, and data reliability^12,13. SnRNA-seq has been extensively utilized in studies focusing on AD and exploring the role of glial cells^14,15,16. Its wide adoption in these research areas highlights its effectiveness in providing valuable insights into the transcriptional characteristics of glial cells in AD. By leveraging the advantages of snRNA-seq, researchers have been able to uncover crucial information regarding the involvement of glial cells in AD pathology and identify potential therapeutic targets. In order to explore sex-specific neuroglial transcriptional characteristics in AD and potential TCMs for sex specificity of AD, this study analyzed snRNA-seq data from the frontal cortex of AD patients from the NCBI GEO public database. Sex-specific differentially expressed genes (DEGs), Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Protein-protein interaction (PPI) network, and gene-TF-miRNA network are further analyzed to reveal key biomarkers and potential pathogenesis. Finally, potential TCMs were suggested, and their active ingredients were displayed with tables by searching the Coremine Medical, TCMIP, and TCMSP databases.

Protocol

Steps 2 to 9 of the analysis were implemented using R software (see Supplementary Figure 1 and Supplementary File 1), while the remaining steps were executed on the online platforms. The details of the databases used in this protocol (along with the weblinks) are provided in the Table of Materials.

1. Data acquisition

1. Access the publicly available Gene Expression Omnibus (GEO) database at the National Center of Biotechnology Information.
2. Search the GEO data named Alzheimer's disease in the search box.
3. Select the Top Organisms as Homo sapiens on the right side.
   NOTE: The search results were data on Alzheimer's disease in Homo sapiens.
4. After filtering the searched information, download the GSE167490 and GSE183068 data files, which encompass features.tsv, barcode.tsv, and matrix.mtx for each individual single nuclei sample. The datasets comprised 34 samples of AD originating from the frontal cortex, with an equal distribution of 17 male samples and 17 female samples (Supplementary Table 1).

2. Sample merging

1. Configure the data paths and sample names accordingly on the computer. Import the 34 downloaded samples and assign gender-specific names to the samples using the function names.
2. Generate Seurat objects for all samples in a batch-processing manner using the functions list and Read10X, specifying the parameters as min.cells = 3 and min.features = 200.
3. Use the RenameCells function to add sample IDs as prefixes to the cell barcodes to preserve cell barcodes during the merging process. This ensured that each cell retained its unique identity and could be traced back to its original sample source after merging.

3. Quality control (QC)

1. Employ the PercentageFeatureSet function to calculate the mitochondrial gene ratios, erythrocyte gene ratios, and ribosome gene ratios for each cell.
2. Store these computed ratios in the metadata using the [[ ]] operator to attach this information directly to each cell's metadata.
3. Utilize the subset function to conduct cell filtration, specifying the parameters as nFeature_RNA > 200, nFeature_RNA < 10000, nCount_RNA < 60000, percent.mt < 10, percent.rb < 5, and percent.HB < 75.
4. Exclude GSM5106107 from the analysis.

4. Batch effect checking

1. Perform data processing.
   1. Normalize the data using the NormalizeData function.
   2. Identify the top 2000 variable features in the dataset using the FindVariableFeatures function.
   3. Conduct principal component analysis (PCA)¹⁷ on the data using RunPCA, retaining 50 principal components.
   4. Generate an elbow plot using the ElbowPlot function to determine the optimal number of dimensions for subsequent analysis. Consider the first 50 dimensions.
   5. Scale the data using ScaleData to ensure that all features are on a comparable scale.
   6. Identify nearest neighbors using FindNeighbors based on 30 dimensions.
   7. Apply the UMAP algorithm using RunUMAP to reduce the dimensionality of the data to 30 dimensions.
2. Visualize the processed data using the DimPlot function with the reduction parameter set to umap and the group.by parameter set to orig.ident.
   NOTE: This step could generate a plot visualizing the data in the reduced UMAP space, grouped by the original cell identities. Upon examining the UMAP plots, it became apparent that there was a presence of batch effect. The distinct clustering or separation of cells based on their batch or experimental origin suggested that the experimental batches had influenced the gene expression profiles.

5. Data integration

1. Normalize and standardize the data using the SCTransform function.
2. Apply the harmony algorithm¹⁸ to integrate the remaining 33 single-nuclei data. Use the SCT assay for integration and set the maximum number of harmony iterations to 20.
3. Use FindClusters function with a resolution parameter set to 0.07 to identify distinct clusters within the data.
4. Employ the RunUMAP function with a specified number of dimensions (dims = 30) to further reduce the dimensionality of the data and visualize the clusters in a lower-dimensional space.

6. Cell type annotation

1. Collect the marker genes (Supplementary Table 2) of cells through an extensive review of the existing literature.
2. Following the identification of cellular cluster heterogeneity, classify the type of each cluster cell by the marker genes expressed specifically.
3. Present various cellular types with UMAP visualization using ggplot2 package, where oligodendrocyte was highlighted with the color code #DB7093, excitatory neuron with #FF69B4, astrocyte with #1874CD, microglia with #63B8FF, oligodendrocyte precursor cell with #DB7093, inhibitory neuron with #FFC0CB, and endothelial cell with #FF69B4.
4. Calculate the proportions of each cell type stratified by gender.

7. Glial cell data extraction

1. Extract astrocyte data from the integrated bulk data using subset function.
2. Extract microglia data from the integrated bulk data using subset function.
3. Extract oligodendrocyte data from the integrated bulk data using subset function.

8. Glial sex-specific differentially expressed genes (DEGs) capturing

1. Identify sex-specific DEGs of astrocytes using the FindMarkers function (ident.1 = male, ident.2 = female, group.by = group.sum, assay = RNA) with threshold values: p-value < 0.05 and |avg_log2FC| > 30. Label the up-regulated DEGs as Up, down-regulated DEGs as Down, and the rest as Stable.
   1. Visualize the DEGs using the ggplot function, with the x-axis representing the difference in percentage between two conditions (pct.1 - pct.2), and the y-axis depicting the avg_log2FC. The up-regulated genes were highlighted by using the color PaleVioletRed, the down-regulated genes with Pink, and the stable genes with DodgerBlue3.
2. Identify sex-specific DEGs of microglia using the FindMarkers function (ident.1 = male, ident.2 = female, group.by = group.sum, assay = RNA) with threshold values: p-value < 0.05 and |avg_log2FC| > 1. Label the up-regulated DEGs as Up, down-regulated DEGs as Down, and the rest as Stable.
   1. Visualize the DEGs using the ggplot function, with the x-axis representing the difference in percentage between two conditions (pct.1 - pct.2), and the y-axis depicting the avg_log2FC. The up-regulated genes were highlighted by using the color OrangeRed, the down-regulated genes with LightSalmon, and the stable genes with SteelBlue1.
3. Identify sex-specific DEGs of oligodendrocyte using the FindMarkers function (ident.1 = male, ident.2 = female, group.by = group.sum, assay = RNA) with threshold values: p-value < 0.05 and |avg_log2FC| > 10. Label the up-regulated DEGs as Up, down-regulated DEGs as Down, and the rest as Stable.
   1. Visualize the DEGs using the ggplot function, with the x-axis representing the difference in percentage between two conditions (pct.1 - pct.2), and the y-axis depicting the avg_log2FC. The up-regulated genes were highlighted by using the color DeepPink, the down-regulated genes with HotPink, and the stable genes with DeepSkyBlue3.

9. Functional enrichment analyses of sex-specific DEGs

1. Perform Gene ontology (GO) enrichment analysis on sex-specific DEGs for each glial cell type using the enrichGO function. Set the following parameters: OrgDb = org.Hs.eg.db, keyType = SYMBOL, ont = ALL, pAdjustMethod = BH, pvalueCutoff = 0.01, and qvalueCutoff = 0.05.
2. Convert the gene symbols into corresponding gene IDs using function bitr. Conduct Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on sex-specific DEGs for each glial cell type by using enrichKEGG function. Adjust the settings as follows: organism = hsa, keyType = kegg, pAdjustMethod = BH, pvalueCutoff = 0.01, and qvalueCutoff = 0.05.

10. Frequency statistics of glial DEGs in go and kegg pathways, venn diagrams of each glial sex-specific DEGs, and PPI network construction

1. Calculate the frequency of glial sex-specific DEGs in GO and KEGG pathways using a frequency histogram.
2. Access the STRING database to construct the PPI networks.
3. Choose the Multiple proteins. Search for the List Of Names in the search box. Set "Organisms" as Homo sapiens.
4. Review the list of proteins obtained from the search. Click on Continue to proceed.
5. Export the PPI networks by selecting the download option, preferably in PNG format with higher resolution.
6. Visualize the distribution of co-expression for the top sex-specific genes by using the Venn diagrams.
7. Identify the shared gene(s) as key gene(s) in the study based on the Venn diagram analysis.

11. Multifactor regulatory network construction

1. Access the NetworkAnalyst.
2. Click on Gene List Input and specify organism as H. sapiens (human). Set ID type as Official Gene Symbol. Enter the gene name in the search field, and then click on Upload and Proceed.
3. Select Gene-miRNA Interactions and choose miRTarBase v8.0. Confirm the selection by clicking on OK.
4. Proceed to TF-gene Interactions and select the ENCODE database. Click on OK to confirm the selection.
5. Next, navigate to TF-miRNA Coregulatory Network and click on OK to proceed.
6. Finally, choose Proceed to generate the multifactor regulatory network incorporating gene-miRNA interactions and TF-gene interactions.

12. Gene and target TCM analysis

1. Access the Coremine Medical online database.
2. Enter the specific gene name and select the corresponding gene with the suffix gene/protein, human into the search box under the Explore section.
3. Navigate the Drugs section and identify the TCMs associated with the searched drugs.
   NOTE: Statistically significant drugs were marked in blue.
4. Determine the top five TCMs based on their "significance" value as therapeutic TCMs.

13. Research summary of TCM ingredients in targeting key gene

1. Access the Integrative Pharmacology-based Research Platform of Traditional Chinese Medicine (TCMIP) and Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP). Enter the names of the herbs into the search bar to retrieve their corresponding ingredients.
2. Retrieve the ingredients in the PubMed database with a time limit up until April 10^th, 2023. The search terms used included the Molecule Name in TCMSP and Chemical Components in TCMIP as search terms, and were limited to articles published in English.
3. Summarize and analyze the herbs and their corresponding ingredients acting on AD.

14. Confirming of targeting TCMs' treatment function in the sex specificity of AD

1. Import herbs into the TCMIP and navigate to the corresponding description page.
2. Utilize the Export data function and select the CSV format to download the enrichment terms of GO - Biological Process, GO - Cellular Component, GO - Molecular Function, and Reactome Pathway.
3. Visualize the downloaded enrichment terms for each herb using bar charts.

Results

SnRNA-seq analysis of frontal glial transcriptome profiles and the annotation of cell types
In total, 220,095 nuclei and 32,077 genes in the frontal cortex of 17 male AD and 17 female AD were obtained (Figure 1A). UMAP plot visualized the total single-nuclei frontal transcriptomes displaying distinct types of nuclei after dimension reduction analysis (Figure 1B). Total numbers of annotated nuclei captured by gender were shown, which made the...

Discussion

Gender-specificity has been identified in epidemiology, pathology, and clinical manifestation of AD¹⁹. Here, we confirmed the potential pathological mechanism of the "hormone-synapse-neuron axis" from gender-specific glial genes and related pathways in AD patients. NLGN4Y was the only shared gene in the three glia and was chosen as the biomarker for the gender specificity of AD. TF and miRNAs regulating NLGN4Y were strongly linked to gender differences and the development of the nervous system. Fu...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Database
Coremine Medical database	Jointly developed by Norway, the Chinese Academy of Sciences, the Chinese Academy of Medical Sciences, the National Medical Library of the United States and other institutions		When you explore concepts in CoreMine Medical you access a database that is structured to relate important concepts, ranked by statistical relevance, to your topic. For example, if you type in "Alzheimer disease," in addition to retrieving documents and resources that discuss the disease, you will be able to view networks and lists that show how your query concept is related to other bio-medical concepts. This provides an overview of concepts that relate to your search as well as being an interface for navigating information on these concepts. Weblink: https://coremine.com/medical/
Gene Expression Omnibus (GEO)	National Center for Biotechnology Information in the United States (NCBI)		GEO is a public functional genomics data repository supporting MIAME-compliant data submissions. Array- and sequence-based data are accepted. Tools are provided to help users query and download experiments and curated gene expression profiles. Weblink: https://www.ncbi.nlm.nih.gov/geo/
Integrative Pharmacology-based Research Platform of Traditional Chinese Medicine (TCMIP, version: 2.0)	None		Introduction to the Integrated Pharmacology Based Network Computational Research Platform for Traditional Chinese Medicine [TCMIP v2.0], http://www.tcmip.cn/ ) It is an intelligent data mining platform based on the online database of the Encyclopedia of Traditional Chinese Medicine (ETCM), which integrates medical big data management and pharmacological computing services. It aims to reveal the scientific connotation of traditional Chinese medicine theory and the scientific value of original thinking in traditional Chinese medicine, summarize and pass on the experience of famous doctors, control the quality of traditional Chinese medicine, explain the principles of traditional Chinese medicine action, research and development of new Chinese medicine, especially the discovery and optimization of modern drug combinations, Provide a strong data foundation and analytical tools. Based on TCMIP v1.0, a comprehensive upgrade is implemented, including five major databases and seven functional modules. Through system integration and module integration, a comprehensive analysis of the multi-level correlation of the "disease syndrome prescription" interaction network can be quickly achieved. As an intelligent data mining platform, TCMIP v2.0 will provide a strong data foundation and analysis platform for revealing the scientific connotation of traditional Chinese medicine theory and the scientific value of original thinking in traditional Chinese medicine, summarizing and inheriting the experience of famous doctors, quality control of traditional Chinese medicine, elucidating the principles of traditional Chinese medicine action, research and development of new traditional Chinese medicine drugs, especially modern drug combination discovery and optimization. Weblink: http://www.tcmip.cn/TCMIP
NetworkAnalyst	None		Networkanalyze is an online visualization analysis platform for gene expression analysis and meta-analysis. It can perform comparative, quantitative, differential and enrichment analysis of gene expression, protein-protein interaction analysis, integration analysis of multiple datasets, and can also draw high-value images such as PCA, protein-protein interaction network diagram, heatmap, volcano diagram, Wayne diagram, etc. Weblink: https://www.networkanalyst.ca/NetworkAnalyst/
PubMed database	National Center for Biotechnology Information in the United States (NCBI)		The Pubmed database is a biomedical literature database maintained by the National Library of Medicine (NLM) in the United States, aimed at providing the latest medical research results to scientists, doctors, researchers, and students worldwide. This database collects biomedical literature from around the world, including journal articles, papers, books, etc. As of now, the Pubmed database has collected over 30 million articles and is continuously updated every week. Weblink: https://pubmed.ncbi.nlm.nih.gov/
R software	Ross Ihaka and Robert Gentleman		R is a language and environment for statistical computing and graphics. It is a GNU project which is similar to the S language and environment which was developed at Bell Laboratories (formerly AT&T, now Lucent Technologies) by John Chambers and colleagues. R can be considered as a different implementation of S. There are some important differences, but much code written for S runs unaltered under R. Weblink: https://www.r-project.org/
Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP)	Zhejiang Jiuwei Health Co., Ltd		TCMSP is not only a data repository, but also an analysis platform for users to comprehensively study Traditional Chinese Medicines (TCM): including identification of active components, screening of drug targets and generation of compounds-targets-diseases networks, as well as the detailed drug pharmacokinetic information involving drug-likeness (DL), oral bioavailability (OB), blood-brain barrier (BBB),intestinal epithelial permeability (Caco-2), ALogP,fractional negative surface area (FASA-) and number of  H-bond donor/acceptor  (Hdon/Hacc). So far, TCMSP has attracted broad attentions and several groups have published more than 10 papers by using our TCMSP database within about one year. Weblink: https://tcmsp-e.com