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 sum of included 58,902 astrocytes, 14,265 microglia, 77,466 oligodendrocytes, 3,520 endothelial, 25,252 excitatory neurons, 31,268 inhibitory neurons and 9,422 oligodendrocyte progenitor cells (Figure 1C). Average expressions of the known cell type markers for each glia were projected on the UMAP plots to identify the cell populations (Figure 1D).

The sex-specific DEGs in astrocytes
58,902 astrocytic nuclei were analyzed (Figure 2A), with 27,504 from male AD (46.69%) and 31,398 from female AD (53.31%). Sex-specific DEGs revealed upregulation of 138 genes, including DST, CACNA2D3, and AC016831.7, downregulation of 105 genes, including CNTN5, RORA, RASSF8, and CADM2, and unchanged 4995 genes (Figure 2B). Further, GO and KEGG analyses identified that those DEGs primarily concentrated in the pathways of neurons, synapses, and hormones, with neuronal pathways including regulation of neuron project development, neuron spine, etc., synaptic pathways including synaptic organization and glutamatergic/cholinergic synapse, and hormone pathways including thyroid hormone synthesis and insulin secretion (Figure 2C,D). Thirty genes with top frequency were obtained (Figure 2E), with PLCB1 ranking the first. PPI networks were constructed to explore the relationship between these genes (Figure 2F) and identified DLG2, CAMK2D, CALM2, and PRKACB as core genes. UMAP plots displayed differences of selected DEGs: KCND2, CAMK2D, MT3, LINC00278, and XIST were higher in female AD and lower in male AD, while NLGN4Y, DLG2, PRKACB, CALM2, UTY, and TTTY14 were opposite (Figure 2G). MT3, CALM2, DLG2, KCND2, PRKACB, CAMK2D, and NLGN4Y were finally determined as sex-specific DEGs of astrocytes as shown in Table 1.

The sex-specific DEGs in microglia
14,265 microglial nuclei (Figure 3A) were analyzed, with 5,327 (37.34%) from male AD and 8,938 (62.66%) from female AD. Sex-specific DEGs revealed upregulation of 224 genes, including KCNIP4 and LRRTM4, and downregulation of 930 genes, including APOE, MT-CO3, and FTL, and unchanged 13,111 genes (Figure 3B). Further, GO and KEGG analyses identified that those DEGs mainly concentrated in the pathways of neurons, phagosome, hormones, and others, with neuronal pathways including neuron-to-neuron synapse, phagocytic pathways including phagosome and regulation of phagocytosis, hormone pathways including estrogen signaling pathway, oxytocin signaling pathway, etc., and others including regulation of inflammatory response, learning or memory, amyloid-beta clearance, etc. (Figure 3C,D). Thirty genes with top frequency were obtained, with TLR2 and TREM2 tied for second place (Figure 3E). PPI networks were constructed to explore the relationship between these genes (Figure 3F) and identified ACTB, APP, and FYN as core genes. UMAP plots displayed differences of selected DEGs: APP, FOS, XIST, and CTSD were higher in female AD and lower in male AD, while NLGN4Y, TREM2, LINC0028, APOE, UTY and TTTY14 were opposite (Figure 3G). TREM2, FOS, APOE, APP, and NLGN4Y were finally determined as sex-specific DEGs of microglia, as shown in Table 2.

The sex-specific DEGs in oligodendrocyte
77,466 oligodendrocytic nuclei were analyzed (Figure 4A), with 42,469 from male AD (54.82%) and 34,997 from female AD (45.18%). Sex-specific DEGs revealed upregulation of 384 genes, including PCDH9, MT-CO1, NEAT1, and NPAS3, downregulation of 188 genes, including FRMD4A, PLP1 and LSAMP, and the rest 76,894 genes unchanged (Figure 4B). Further, GO and KEGG analyses identified that those DEGs primarily concentrated in the pathways of neurons, synapses, and hormones, with neuronal pathways including neuron spine, synaptic pathways including neuron-to-neuron synapse and glutamatergic/dopaminergic synapse, and hormone pathways including neurotrophin signaling pathway, aldosterone synthesis and secretion, calcium signaling pathway, etc (Figure 4C,D). Thirty genes with top frequency were obtained (Figure 4E), with GRIN2A and PSEN1 tied for second place. PPI networks were constructed to explore the relationship between these genes (Figure 4F) and identified GRIN2A and GRIA2 as core genes. UMAP plots displayed differences in selected DEGs: GRIN2A, ITPR2, GNAS, and XIST were higher in female AD and lower in male AD, while NLGN4Y, UTY, and TTTY14 were opposite (Figure 4G). GRIN2A, ITPR2, GNAS, and NLGN4Y were finally determined as sex-specific DEGs of oligodendrocytes, as shown in Table 3.

Interaction between key genes and the top 30 Genes of glial cells, and NLGN4Y as the common shared gene
The Venn diagrams and the PPI networks provided an overview of the close interaction between key genes (Figure 5A,B) and the top 30 genes (Figure 5C, D) of each glial cell. Results showed that ACYB, APP, JUN, PRKACB, and DLG2 were at the core of the PPI network, and NLGN4Y was a common shared gene for sex-specific DEGs in all glial cells.

Gene-TF-miRNA network construction
The NLGN4Y-TF-miRNA network contained 13 nodes and 12 edges (Figure 6A). NLGN4Y was regulated by 1 TF, namely CTCF, and 11 miRNAs, including hsa-miR-185, hsa-miR-137 and hsa-miR-9.

Displaying target drug and TCMs of NLGN4Y with network
A total of 1 NLGN4Y target drug and 64 indirect target TCMs were retrieved in Coremine Medical. TCMs with statistical significance in the results were marked with blue. Drug Antithrombin III and five TCMs with the network, namely Heikunbu, Wulingzhi, Xiazhicao, Shuizhi, and Mahuang, were visualized, which were considered statistically significant (Figure 6B).

Effects of target TCMs and corresponding active ingredients on AD
For Kunbu, 10 ingredients in TCMIP and 48 ingredients in TCMSP were retrieved. Through searching the PubMed database, 5 ingredients were retrieved related to AD: fucosterol, saringosterol, thiamine, stearidonic acid, and phlorofucofuroeckol-A, of which the oral bioavailability (OB) of the first two was ≥30% and similar to the drug sex (DL) was ≥0.18. Details are displayed in Table 4. As for Mahuang, 28 ingredients in TCMIP and 363 ingredients in TCMSP were retrieved. With the same method as Kunbu, 25 ingredients related to AD were retrieved. Among them, 6 ingredients with OB ≥30% and DL ≥0.18: quercetin, eriodictyol, naringenin, taxifolin, stigmasterol, and luteolin, are listed at the top of Table 5.

Hurb-gene-disease network and enrichment analysis
Hurb-gene-disease network is shown in Figure 6C. The shared target gene of Kunbu and Mahuang was ACHE, whose associated diseases were AD and cognitive deficits. The target genes of Kunbu alone were ALKBH3 and ELOVL4, diseases related to which were aging and brain atrophy, respectively.

The results of GO enrichments (BP, MF, and CC) for Kunbu (Figure 7A, Figure 8A, Figure 9A) and Mahuang (Figure 7B, Figure 8B, Figure 9B) are displayed with bar charts. Reactome pathways for Kunbu (Figure 10A) and Mahuang (Figure 10B) are also shown. The enriched terms marked with arrows were related to the "hormone-synapse-neuron" axis, such as steroid hormone receptor activity, neuron projection, and steroid hormone-mediated signaling pathway in Kunbu and chemical synaptic transmission, estrogen-dependent gene expression and nervous system process in Mahuang.

Figure 1: Gene expression data acquiring, single-nuclei RNA-seq profiling, and cell type characterization. (A) Samples obtained from GEO DataSets for analysis preparation. (B) 2-dimensional UMAP plot of sum nuclei (N = 116,101 for male; N = 103,994 for female). (C) Proportions for each type of cell split by gender. (D) Average expression of 5 well-established cell type markers projected on the UMAP plot. Please click here to view a larger version of this figure.

Figure 2: Astrocytes are heterogeneous and have sex-specific transcriptomic changes in Alzheimer's disease. (A) UMAP plot of astrocytic nuclei (N = 59,010). (B) Sex-specific DEGs. (C,D) Circle plots illustrated the significant functionally enriched terms of astrocytic DEGs obtained from GO and KEGG databases (GO: C, KEGG: D). (E) The top 30 most frequent astrocyte-sex DEGs in GO and KEGG pathways. (F) PPI network of the top 30 most frequent astrocyte-sex DEGs in GO and KEGG pathways. (G) Average expression of remarkable sex-specific DEGs projected on the UMAP plots. Please click here to view a larger version of this figure.

Figure 3: Microglia are heterogeneous and have sex-specific transcriptomic changes in Alzheimer's disease. (A) UMAP plot of microglial nuclei (N = 14,265). (B) Sex-specific DEGs. (C,D) Circle plots illustrated the significant functionally enriched terms of microglial DEGs obtained from GO and KEGG databases (GO: C, KEGG: D). (E) The top 30 most frequent microglia-sex DEGs inGO and KEGGpathways. (F) PPI network of the top 30 most frequent microglia-sex DEGs in GO and KEGG pathways. (G) Average expression of remarkable sex-specific DEGs projected on the UMAP plots. Please click here to view a larger version of this figure.

Figure 4: Oligodendrocytes are heterogeneous and have sex-specific transcriptomic changes in Alzheimer's disease. (A) UMAP plot of oligodendrocytic nuclei (N = 77,466). (B) Sex-specific DEGs. (C,D) Circle plots illustrated the significant functionally enriched terms of oligodendrocytic DEGs obtained from GO and KEGG databases (GO: C, KEGG: D). (E) The top 30 most frequent oligodendrocyte-sex DEGs in GO and KEGG pathways. (F) PPI network of the top 30 most frequent oligodendrocyte-sex DEGs in GO and KEGG pathways. (G) Average expression of remarkable sex-specific DEGs projected on the UMAP plots. Please click here to view a larger version of this figure.

Figure 5: Overlaps of DEGs and PPI networks of top frequency DEGs. (A,B) Venn diagram illustrating the key genes for each glia and the corresponding PPI network. (C,D) Venn diagram illustrating the top 30 frequent sex-specific DEGs for each glia enriched in GO and KEGG pathways and the corresponding PPI network. Please click here to view a larger version of this figure.

Figure 6: Network diagram with NLGN4Y and its targeting TCMs as core. (A) Gene-TF-miRNA coregulatory network. (B) Gene-drug-TCM network for NLGN4Y. (C) TCM-gene-disease network of Kunbu and Mahuang. Please click here to view a larger version of this figure.

Figure 7: GO (BP) Enrichment analysis for Kunbu and Mahuang. (A) Top 20 enriched biological processes of Kunbu. (B) Top 20 enriched biological processes of Mahuang. Please click here to view a larger version of this figure.

Figure 8: GO (MF) Enrichment analysis for Kunbu and Mahuang. (A) Top 20 enriched molecular functions of Kunbu. (B) Top 20 enriched molecular functions of Mahuang. Please click here to view a larger version of this figure.

Figure 9: GO (CC) Enrichment analysis for Kunbu and Mahuang. (A) Top 20 enriched cellular components of Kunbu. (B) Top 20 enriched cellular components of Mahuang. Please click here to view a larger version of this figure.

Figure 10: Reactome pathways for Kunbu and Mahuang. (A) Reactome pathways of Kunbu. (B) Reactome pathways of Mahuang. Please click here to view a larger version of this figure.

Table 1: Key genes of the sex-specific DEGs in astrocytes. Please click here to download this Table.

Table 2: Key genes of the sex-specific DEGs in microglia. Please click here to download this Table.

Table 3: Key genes of the sex-specific DEGs in oligodendrocytes. Please click here to download this Table.

Table 4: Kunbu' active ingredients on AD. Please click here to download this Table.

Table 5: Mahuang' active ingredients on AD. Please click here to download this Table.

Supplementary Figure 1: Screenshot for the R software usage. Please click here to download this File.

Supplementary Table 1: Sample information. Please click here to download this File.

Supplementary Table 2: Marker genes. Please click here to download this File.

Supplementary File 1: The R CODE. Please click here to download this File.

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. Furthermore, the target TCMs Kunbu and Mahuang were considered to potentially affect the "hormone-synapse-neuron" axis and were identified to act as therapeutic candidates for AD by regulation of gender specificity.

The Harmony algorithm provided distinct advantages over other integration algorithms by integrating data in rare cells, optimizing memory usage and computational speed for large samples, and accommodating complex experimental designs with diverse cellular sources and technology platforms¹⁸. QC and data integration sections were critical components in enhancing the accuracy and reliability of the analysis, enabling researchers to glean more profound insights into the biology of AD and other neurodegenerative illnesses. The meticulous execution of these steps is pivotal in achieving the overall quality and validity of the results derived from single-nuclei RNA sequencing experiments.

One of the modifications made in this protocol was the use of the SCTransform function to normalize and standardize the data²⁰. This step ensured that the data was comparable and standardized for further analysis. To mitigate the batch effect, an additional step of data integration was taken, which helped to correct the batch effect and improve the accuracy of the results. Troubleshooting was conducted when one of the samples (GSM5106107) showed unusual performance compared to the other samples, and it was excluded from the analysis.

One limitation of this study pertained to the relatively small sample size employed, which may restrict the representativeness of the findings to the broader population. Furthermore, the utilization of data from multiple laboratories and research groups introduced inherent heterogeneity, potentially limiting the reproducibility and generalizability of the results. Although the identification of key regulatory factors for NLGN4Y was achieved through a gene-TF-miRNA coregulatory network, further experimental verification was necessary to ascertain the functional roles and underlying mechanisms attributed to these regulatory factors. While enrichment analysis suggested Mahuang and Kunbu as potential candidates for gender-specific AD treatment, relying solely on bioinformatics predictions was insufficient to establish their efficacy and safety. Additional experiments and clinical studies were imperative to validate the therapeutic effects of these promising drugs.

The approach delineated in this study exhibited considerable potential for diverse research domains, notably in the discernment of sex-specific DEGs within glial cells in AD. Such investigations hold promise for elucidating the nuanced pathogenic disparities between male and female individuals, consequently facilitating the formulation of sex-specific therapeutic interventions. Moreover, the method's applicability extended to other neurodegenerative conditions, offering opportunities to identify sex-specific DEGs and augment comprehension of the diseases' underlying mechanisms.

In conclusion, the study presented a comprehensive and robust protocol for analyzing single-nuclei RNA sequencing data from AD patients, potentially paving the way for the development of tailored therapeutic interventions. Furthermore, the methodology employed in this study may be extended to other neurodegenerative conditions, offering opportunities to uncover sex-specific gene expression patterns and deepen the understanding of these complex diseases.

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