
==== Front
Genes Dis
Genes Dis
Genes & Diseases
2352-4820
2352-3042
Chongqing Medical University

S2352-3042(24)00134-X
10.1016/j.gendis.2024.101337
101337
Full Length Article
Gene interactions analysis of brain spatial transcriptome for Alzheimer's disease
Wang Shengran abc
Greenbaum Jonathan c
Qiu Chuan c
Swerdlow Russell H. d
Haeri Mohammad d
Gong Yun c
Shen Hui c
Xiao Hongmei hmxiao@csu.edu.cn
ef⁎⁎
Deng Hongwen hdeng2@tulane.edu
c⁎
a Reproductive Medicine Center, The Sun Yat-sen Memorial Hospital, Sun Yat-sen University, Guangzhou, Guangdong 510120, China
b Center for System Biology, Data Sciences and Reproductive Health, School of Basic Medical Science, Central South University, Changsha, Hunan 410008, China
c Tulane Center of Biomedical Informatics and Genomics, Deming Department of Medicine, Tulane University School of Medicine, Tulane University, New Orleans, LA 70112, USA
d Department of Pathology and KU Alzheimer's Disease Research Center, University of Kansas Medical Center, Kansas City, KS 66160, USA
e Institute of Reproductive & Stem Cell Engineering, School of Basic Medical Science, Central South University, Changsha, Hunan 410008, China
f Center of Reproductive Health, School of Basic Medical Science, Central South University, Changsha, Hunan 410008, China
⁎ Corresponding author. hdeng2@tulane.edu
⁎⁎ Corresponding author. Institute of Reproductive & Stem Cell Engineering, School of Basic Medical Science, Central South University, Changsha, Hunan 410008, China. hmxiao@csu.edu.cn
22 5 2024
11 2024
22 5 2024
11 6 10133720 11 2022
25 12 2023
21 2 2024
© 2024 The Authors. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co., Ltd.
2024
Chongqing Medical University
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Recent studies have explored the spatial transcriptomics patterns of Alzheimer's disease (AD) brain by spatial sequencing in mouse models, enabling the identification of unique genome-wide transcriptomic features associated with different spatial regions and pathological status. However, the dynamics of gene interactions that occur during amyloid-β accumulation remain largely unknown. In this study, we performed analyses on ligand-receptor communication, transcription factor regulatory network, and spot-specific network to reveal the dependence and the dynamics of gene associations/interactions on spatial regions and pathological status with mouse and human brains. We first used a spatial transcriptomics dataset of the AppNL-G-F knock-in AD and wild-type mouse model. We revealed 17 ligand-receptor pairs with opposite tendencies throughout the amyloid-β accumulation process and showed the specific ligand-receptor interactions across the hippocampus layers at different extents of pathological changes. We then identified nerve function related transcription factors in the hippocampus and entorhinal cortex, as well as genes with different transcriptomic association degrees in AD versus wild-type mice. Finally, another independent spatial transcriptomics dataset from different AD mouse models and human single-nuclei RNA-seq data/AlzData database were used for validation. This is the first study to identify various gene associations throughout amyloid-β accumulation based on spatial transcriptomics, establishing the foundations to reveal advanced and in-depth AD etiology from a novel perspective based on the comprehensive analyses of gene interactions that are spatio-temporal dependent.

Keywords

Alzheimer's disease
Amyloid-β
Ligand-receptor
Spatial transcriptomics
Transcription factor
==== Body
pmcIntroduction

The relationship of amyloid plaques to the neurodegenerative process and related gene expression changes is a central question in Alzheimer's Disease (AD) research. Spatial transcriptomics (ST) is an attractive approach that provides unbiased transcriptome profiling of sequenced molecules with spatial localization and focused image signals like amyloid plaques.1,2 Previous studies used mouse brains from AppNL-G-F (AD model) and C57BL/6 (wild-type/WT) mice at 3, 6, 12, and 18 months of age to explore ST with the corresponding quantitative amyloid plaque index information, and demonstrated alterations in the plaque-induced and oligodendrocyte gene networks in different phases of AD. However, the dynamic changes of gene interactions between ligand-receptor (L-R) pairs, mediated by transcription factor (TF) regulation, and various degrees of gene associations/interactions at the transcriptome level throughout the disease process in different brain regions remain largely unexplored.

Defining a more detailed landscape of how different types of gene interactions are altered in the neuronal structures of the brain may contribute to a better understanding of AD etiology.2, 3, 4 To accomplish this goal, multiple approaches such as CellPhoneDB5 and single-cell regulatory network inference and clustering (SCENIC)6 have been used to resolve spatio-temporal differences in gene interactions in the same ST dataset. In addition, ST allows gene expression profile selection to be anchored to some particular brain regions of interest, which in this study were the hippocampus (HP) and entorhinal cortex (ENTI). The HP and ENTI play central roles in AD pathology. Hippocampal deficits in amyloid-β (Aβ)-related rodent models of AD have revealed synaptic, behavioral, and circuit-level defects.7,8 ENTI is one of the earliest affected brain regions.9 Aβ proteins are seen primarily in the ENTI in mild AD and “spread” to the HP and other cortical areas as the disease progresses.10 Selective overexpression of mutant amyloid precursor protein (APP) in ENTI has been reported to cause an aberrant excitatory cortico-hippocampal network activity leading to behavioral abnormalities.10,11

Here we used the ST dataset of AppNL-G-F and C57BL/6 mice with annotations for brain region,12 AD phenotype, and Aβ index to explore the genome-wide transcriptomic features associated with different spatial regions and pathological status. In general, our work revealed 17 L-R pairs with opposite tendencies throughout the Aβ accumulation process. We also showed the specific L-R interactions across the HP layers at different extents of pathological changes which may be related to their unique functions cross these region pairs. We searched the convergent functional genomic (CFG) scores of these unique L-R pairs in the AlzData database13 to validate their AD-related functions. Next, we identified nerve function related TFs in the HP and ENTI in order to explore their changes with disease progression. We then calculated genes with different degrees (interactive connection with other genes) based on the network degree matrix (NDM) from spot-specific network (SSN) analysis. The purpose was to reveal unique gene connections in different brain regions and age groups of different phenotypes (WT/AD). These results identified novel and known L-R pairs/genes that may contribute to disease pathology and revealed the strong and coordinated cellular responses in the amyloid plaque niche based on a comprehensive analysis of gene interactions that are spatio-temporal dependent.

Materials and methods

Data sources and workflow

The whole brain ST dataset used in the present study was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds) with accession number GSE152506. This dataset12 integrated 10,327 transcriptomics profiles (10,327 spots/tissue domains) from AppNL-G-F and C57BL/6 male mice at 3, 6, 12, and 18 months of age as one gene expression matrix annotated with Brodmann area, pathological status (Aβ index), and age information. Twelve mouse samples were used in this dataset: mouse number (n) = 2 per phenotype (WT or AD) for 3 and 18 months of age; n = 1 per phenotype (WT or AD) for 6 and 12 months of age; n = 12 for total. Two ST slides were included per brain sample from each of the 3- and 18-month mice. Only one ST slide was included per brain sample from each of the 6- and 12-month mice. The total ST slide number was 20, including 4 in WT_03 (3 months), 1 in WT_06 (6 months), 1 in WT_12 (12 months), 4 in WT_18 (18 months), 4 in AD_03, 1 in AD_06, 1 in AD_12, and 4 in AD_18. Spatial transcriptomics experiments were performed following the Library Preparation Manual (Spatial Transcriptomics, Stockholm, Sweden).14 The diameter of a tissue domain is 100 μm. Data integration and normalization (by the EdgeR “cpm” function) for the multiple slides for the same age and phenotype groups were done in the original paper12 as well as the log-cpm matrix downloaded from GEO.

The HP ST dataset used in the present study was available at Mendeley https://doi.org/10.17632/6s959w2zyr.1y. This dataset came from a study2 which included two types of AD models. The triple transgenic AD model (3xAD) expresses three human gene variants: APP K670N/M671L, MAPT P301L, and PS1 M146V. 3xPB mouse is a modified version of the 3xAD strain formed by adding a deficiency in DNA polymerase beta (Polb). The 3xPB AD model displays more severe deficits in memory, learning, and long-term potentiation as well as higher levels of DNA damage and cell death. All mice are male (11–13 months of age) on a C57BL/6J background. For HP ST, three mice per model and two sections per mouse were included.

The ENTI single-nuclei RNA sequencing (snRNA-seq) dataset used in the present study was obtained from the GEO database with accession number GSE138852. This dataset came from a study15 in which post-mortem AD and non-diseased age-matched human ENTI tissues were collected from the Victorian Brain Bank. Disease or control group allocation was based on overall amyloid and tau pathology. In total, this dataset included six AD and six control gene expression matrices (6 AD patients and 6 control individuals) from snRNA-seq with annotations for cell type and phenotype.

The functions of genes in significant L-R pairs or TF regulatory networks were validated in an AD database, AlzData (http://www.alzdata.org/).13 AlzData provides various types of information for target genes, including the expression of the target gene regulated by AD genetic variants; significant physical interactions of proteins of the target genes with AD core genes such as APP, PSEN1, PSEN2, APOE, or MAPT13; and differentially expressed in AD mouse models before AD pathology emergence and the correlation of target gene expression with AD pathology in amyloid beta or tau line AD mouse models. One CFG score is added if any of the above evidence is significant. The CFG score ranged from 0 to 5.

In brief, we used the whole brain ST dataset (GSE152506) to perform L-R interaction analysis within spots and across brain regions (HP layers). Next, TF network and SSN analysis were performed in HP and ENTI areas. Then AlzData database was used to validate the target genes' AD-related functions by the CFG score. We also used the independent HP ST dataset to validate the expression of unique L-R interactions that are only differentially expressed in one HP region pair in the whole brain ST dataset. The ENTI snRNA-seq dataset was used to validate the different expression of important TFs in human AD patients compared with the control (Fig. 1A).Figure 1 L-R pair clustering and enrichment analysis. (A) Workflow of the research. (B) Dynamic changes of three L-R clusters in the 3-month spots through Aβ accumulation. The X-axis represents the expression changes and the Y-axis represents different Aβ accumulation groups. Each line represents one L-R pair. The line color represents the membership score for each L-R pair from low (blue) to high (red). (C–F) Neuron (C), axon (D), synapse (E), and amyloid (F) related GO analysis results of L-R genes of three L-R clusters in 3 months. The outline color represents different clusters. The fill color represents adjusted P-values from low (blue) to high (yellow). (G) Dynamic changes of four L-R clusters in the 18-month spots through Aβ accumulation. The X-axis represents the expression changes and the Y-axis represents different Aβ accumulation groups. Each line represents one L-R pair. The line color represents the membership score for each L-R pair from low (blue) to high (red). (H–K) Neuron (H), axon (I), synapse (J), and amyloid (K) related GO analysis results of L-R genes of four L-R clusters in 18 months. The outline color represents different clusters. The fill color represents adjusted P-values from low (blue) to high (yellow).

Figure 1

Analysis of L-R communication

We used CellPhoneDB,5 a public knowledge base of ligands-receptors that considers interacting partners as binary interactions annotated by IUPHAR (http://www.guidetopharmacology.org/) and cytokines, hormones, and growth factors interacting with receptors annotated by the iMEX consortium (https://www.imexconsortium.org/) to systematically analyze the intra-spots/cross brain regions in their communication network, and calculate the average log gene expression level and communication significance of each L-R pair. To perform statistical inference of L-R specificity, a null distribution for each L-R pair mean was generated by random permutation of spot labels (1000 times by default). By calculating the proportion of the means that are “as or more extreme” than the actual mean, CellPhoneDB obtained a P-value for the likelihood of spot-type specificity of a given L-R pair. Predicted interaction pairs with P-values <0.05 were considered as significant as defined in the original analyses.5

The intra-spot molecular interactions analysis only used CellPhoneDB to get the expression value of each L-R pair within each spot. Specifically, CellPhoneDB only provided the L-R database and calculated the average log gene expression level of ligand and receptor genes in the same spot in this part.

In the cross-region L-R analysis, CellPhoneDB compared the average log gene expression level of each L-R pair among different region pairs. P-values <0.05 were considered statistically significant (Fig. S1A). The output results (Fig. 4A, B) were the total significant L-R pair numbers (Fig. S1B) in each pair of regions.

Our L-R communication analysis needed more detailed subgroup classification (regarding certain brain regions in the same age, Aβ index, and phenotype group) information than TF regulatory network and SSN analyses as input. However, the expression data missing rate in some brain regions at certain Aβ accumulation levels was higher in 6- and 12-month-age groups than the two other age groups. This is because the total numbers of spots in these groups were limited in the whole brain ST dataset compared with the 3- and 18-month-age groups. So, the 6- and 12-month-age groups were not involved in our L-R communication analysis.

Soft clustering analysis

Soft clustering methods such as Mfuzz can assign genes to several clusters using the fuzzy c-means algorithm with time/process-course data on the gene expression and may be more accurate and robust than hard clustering.16 To analyze the expression trend of intra-spot L-R pairs with Aβ accumulation, the Mfuzz package in R was used to perform soft clustering analysis for visualizing the Aβ level-dependent expression patterns of L-R pairs and assign L-R pairs to several cluster modules with different expression tendencies through Aβ accumulation. Parameter m was calculated by the “mestimate” function.

Gene enrichment analysis

To identify the significantly enriched pathways of significant L-R gene pairs, we used the clusterProfiler R package to perform Gene Ontology (GO) enrichment analysis. Only terms showing adjusted P-values less than 0.05 (adjusted for multiple testing by using the BH method) were considered as significantly enriched.

SCENIC analysis

First, putative TF motifs were identified by RcisTarget (hg19-tss-centered-10kb-7species.mc9nr.feather, https://resources.aertslab.org/cistarget) software in Bioconductor. The SCENIC algorithm was then used to construct gene networks and scored regulon (gene regulatory modules) activities.6 We used GENIE3 in Bioconductor to identify co-expressed gene modules and infer potential TF targets for each module with significantly enriched motifs. RcisTarget identifies enriched transcription factor binding motifs and candidate TFs for a gene list. The motifs that are annotated to the corresponding TF and obtain a normalized enrichment score >3.0 are retained. Regulons were identified from co-expression and DNA motif analyses. Then, the AUCell function calculated an area under the curve across the ranking of all genes in a particular spot. By calculating whether the input gene subset (TFs and their putative targets in a particular regulon) was enriched at the top of the ranking, AUCell scored the regulon activity in each spot as a whole, instead of only the TF or individual genes. Regulatory network of high confident annotation interactions (the direct TF-target gene interactions annotated by the cisTarget databases, highConfAnnot = “True”) were visualized in Cytoscape (v3.9.1).

Protein-protein interaction network

The protein-protein interaction network was constructed by the Search Tool for the Retrieval of Interacting Genes (STRING) database. The String database provides information about known and predicted protein interactions based on correlation, regulation relationships, or protein binding validated in Co-IP or other experiments. The Cytoscape software v.3.9.1 was used to construct the protein-protein interaction network based on STRING analysis results.

SSN analysis

We performed SSN analysis in Matlab software as previously described.17 Briefly, the SSN method calculates a statistic ρˆxy(k) to assess the inter-relationships among gene x and gene y in spot k. “n” was the total spot number. We set n x(k) = n y(k) = 0.1 n as default in the scatter plot.17 The coefficient 0.1 denotes the box size. After hypothesis testing based on ρˆxy(k) to determine the significance of each gene interaction, the total number of significant interactions including gene x is returned as NDM value for gene x. We then used the Wilcoxon rank-sum test to identify genes with significantly different NDM degrees (false discovery rate <0.05) in WT phenotype compared with AD, in different brain regions and age groups.

Results

CellPhoneDB analysis reveals neural related intra-spot L-R interactions

To explore the intra-spot molecular interactions, we clustered L-R pairs by different tendencies through the Aβ accumulation process. All ST tissue domains were separated into WT and four quantiles of AD according to the Aβ index (increasing Aβ accumulation level from L1 to L4; Table S1). All complex information involved in L-R pairs is provided in Table S2. There were three clusters (Fig. 1B and Table S3) in the 3-month group. Cluster 1 showed a downward tendency from WT to L4. Cluster 2 showed an obvious up-regulation from L3 to L4. Compared with any other spot with Aβ accumulation (L1–L4), L-R pairs in cluster 3 were lower expressed in WT spots. GO analysis revealed the different enrichment results of neuron, axon, synapse, and amyloid related terms in each cluster of the 3-month group (Fig. 1C–F). There are four clusters (Fig. 1G and Table S3) in the 18-month group. Cluster 1 showed an upward tendency from WT to L4. Cluster 2 showed an obvious down-regulation in L1, while cluster 3 was mainly up-regulated in L1. Genes in cluster 4 L-R pairs were expressed at lower levels in WT spots than any other spots with Aβ accumulation (L1–L4). GO analysis revealed the different enrichment results of neuron, axon, synapse, and amyloid related terms in each cluster of the 18-month group (Fig. 1H–K). We explored L-R pairs with both genes involved in the same kind of term related to neuron, axon, synapse, or amyloid and found that cluster 3 in the 3-month group and cluster 1 in the 18-month group have the most L-R pairs related to amyloid (Fig. S2). These results showed the dynamic changes of neural related intra-spot L-R interactions through the Aβ accumulation process.

17 L-R pairs showed opposite tendencies through Aβ accumulation in the 3-month and 18-month groups

Cluster 1 in the 3-month group and cluster 1 in the 18-month group had 17 common L-R pairs (Table 1). However, these pairs showed opposite tendencies with a general inhibition from WT, L1/L2 (L1 and L2 Aβ index group), to L3/L4 (L3 and L4 Aβ index group) in 3 months and a general activation in 18 months among different brain regions (Fig. 1B, G; Fig. 2A). The heatmap showed a similar clustering result of these L-R pairs (except for EPHB2_EFNB2; Fig. 2B) in different brain regions. These results suggested that the different active patterns of 17 common L-R pairs were mainly related to Aβ accumulation rather than to different brain regions. EPHA3_EFNA3, EPHA1_EFNA3, and EPHA8_EFNA3 have high correlation scores with each other in the 3-month and 18-month groups, respectively, which revealed the existence of potential functional or regulatory interactions between them (Fig. 2C–E).Table 1 List of 17 L-R pairs that showed opposite tendencies through the Aβ accumulation process in the 3-month and 18-month groups.

Table 1L-R pair	Complex	
C3_aMb2	aMb2: ITGB2/ITGAM	
EFNA1_EPHA8		
FZD9_WNT1		
SPP1_a4b1	a4b1: ITGB1/ITGA4	
INS_INSR		
EPHB4_EFNB1		
TNF_PTPRS		
EPHA4_EFNB1		
EFNB2_EPHA4		
ADGRL1_NRG1		
NTF3_NTRK3		
EPHB6_EFNB2		
NRP1_VEGFA		
EPHB2_EFNB2		
EPHA3_EFNA3		
EPHA1_EFNA3		
EPHA8_EFNA3		

Figure 2 L-R pairs showed opposite tendencies through Aβ accumulation at 3 and 18 months in intra-spot CellPhoneDB analysis. (A) Average expression of 17 L-R pairs with opposite tendencies. Labels showed sample group information. (B) Expression heatmap of L-R pairs with opposite tendencies in different brain regions at 3 (left) and 18 (right) months. (C) Expression correlation between L-R pairs with opposite tendencies at 3 (left) and 18 (right) months. (D) Correlation between EPHA3_EFNA3, EPHA1_EFNA3, and EPHA8_EFNA3 at 3 months. (E) Correlation between EPHA3_EFNA3, EPHA1_EFNA3, and EPHA8_EFNA3 at 18 months. (F) The protein-protein interaction network of 17 L-R pairs with opposite tendencies. Arc color represents four kinds of GO terms. The directed edge is from ligand to receptor. The red double-line marked edge represents complex. (G) Average expression of EPH-related L-R pairs with opposite tendencies. Labels showed sample group information. TH, thalamus; HY, hypothalamus; FB, fiber tract; HPd, dendritic layer of the hippocampus; HPs, somatic layer of the hippocampus; CNU, cerebral nucleus; CTXsp, cortical subplate; OLF, olfactory area; ENTI, entorhinal area; TEP, temporal association area, ectorhinal area, and perirhinal area; AUD, auditory area; SSp, primary somatosensory area; PLT, posterior parietal association area; RSP, retrosplenial area.

Figure 2

Next, we performed protein-protein interaction analysis to further explore the interactions in these 17 L-R pairs (Fig. 2F). In this network, 11 genes from EFNA/EFNB and erythropoietin-producing hepatocellular (EPH) gene families as L-R pairs showed strong connections between each other. Their function was mainly related to neuron, axon, and synapse (Fig. S2). Average expression levels of 9 EFNA/EPFB/EPH-related L-R pairs were concordant with opposite tendencies that the 3-month group showed an inhibition and the 18-month group showed an activation when progressing from WT, L1/L2, to L3/L4 (Fig. 2G).

C3_ITGAM/ITGB2 complex was the only L-R pair that was involved in amyloid related terms in the protein-protein interaction network. ITGAM was a highly expressed marker gene of microglia (Fig. 3A; single-cell RNA-seq result of different cell types from the human brain in AlzData13 (http://www.alzdata.org/, GSE67835)). CD11b encoded by ITGAM and CD18 encoded by ITGB2 formed leukocyte-specific complement receptor 3. iC3b is the activation product of complement C3. Complement receptor 3's space conformation will change when interacting with iC3b to initiate Aβ clearance, phagocytosis, adhesion, or other biological processes18 (Fig. 3B). C3_ITGAM/ITGB2 complex pair was mainly down-regulated at L3/L4 in the 3-month group, and conversely up-regulated at L3/L4 in the 18-month group (Fig. 3C). These results showed that the C3_ITGAM/ITGB2 complex L-R pair has different correlations with Aβ index at different ages.Figure 3 Active pattern of C3-ITGAM/ITGB2 complex. (A)ITGAM expression in different nerve cell types, from the AlzData database. (B) Pattern graph of ITGAM/ITGB2 complex's conformation changes when interacting with C3. (C) Average expression of C3-ITGAM/ITGB2 complex pairs at 3 (top) and 18 (bottom) months. Labels showed sample group information.

Figure 3

Dynamic changes of L-R interactions across HP

L-R interaction networks of HP showed unique features in different age groups or Aβ levels (Fig. 4A, B). The L-R interactions in 3 months and 18 months of both WT and AD groups were both most concentrated on the granule cell layer of the dentate gyrus in Hps (HPs_DG_sp) related region pairs. WT spots in both the 3-month and 18-month groups have more significant interactions across HPs_DG_sp and the pyramidal layer of CA3 in Hps (HPs_CA3_sp) compared with L1/L2 and L3/L4 spots in the same age group. Interactions within the DG polymorph layer in HPd (HPd_DG_po) were more specific in the 18-month spots of the WT group compared with the 3-month spots of the same phenotype group. L3/L4 spots in 3 months showed a general increase in differentially expressed L-R pair numbers compared with WT and L1/L2 spots at this age, especially in HPs_DG_sp itself and HPs_DG_sp to other HP regions. This tendency does not fit in the 18-month group as more region pairs showed less differentially expressed L-R numbers in L3/L4 sports compared with WT and L1/L2 spots. Dot plots showed the top 10 unique L-R pairs (by P-values) that were only significantly differentially expressed in one region pair (Fig. 4C). To explore these gene member's functions, we used the AlzData database to search their scores in AD related terms. FLT1 in the 18-month group demonstrated the highest CFG score (5). CXCL13, C3, FGFR2, and CNR2 in the 3-month group and WNT4, VCAM1, and VEGFA in the 18-month group also showed 4 in the CFG grading (Table 2). These results revealed the potential of unique L-R pairs' AD related functions in different HP layers.Figure 4 L-R interactions in cross the hippocampal regions in CellPhoneDB analysis. (A) L-R interaction networks in different age groups or Aβ levels. The width of the edge represents the number of significantly differentially expressed L-R pairs. The edge direction is from ligand to receptor. Labels showed sample group information. (B) Heatmap of significantly differentially expressed L-R pair numbers between two hippocampal layers. (C) Unique top 10 L-R pairs (by P-values) that are only significantly differentially expressed in one region pair, marked by dots with black outlines. Labels showed sample group information.

Figure 4

Table 2 CFG scores of top 10 unique L-R pair genes in hippocampus analysis.

Table 2Gene	eQTL	GWAS	PPI	Early_DEG	Abeta	Tau	CFG	Note	
FLT1	1	2	PSEN1, PSEN2, APOE	Yes	0.466∗∗	−0.119, ns	5	HP18	
WNT4	4	2	-	Yes	−0.634∗∗∗	−0.326, ns	4	HP18	
VCAM1	2	0	APOE	Yes	0.575∗∗∗	0.437, ns	4	HP18	
VEGFA	2	0	APP, PSEN2, APOE	Yes	−0.215, ns	−0.758∗∗	4	HP18	
CXCL13	1	0	APP, PSEN1, APOE	Yes	0.602∗∗∗	0.727∗∗	4	HP03	
C3	2	0	APP, PSEN1, PSEN2, MAPT	Yes	0.850∗∗∗	0.761∗∗∗	4	HP03	
FGFR2	3	1	PSEN2, APOE	Yes	0.179, ns	−0.387, ns	4	HP03	
CNR2	1	0	APP, PSEN1, APOE	Yes	0.854∗∗∗	0.750∗∗	4	HP03	
ITGA3	1	0	APP, APOE	Yes	−0.284, ns	−0.511, ns	3	HP03, HP18	
AREG	NA	0	PSEN2	Yes	0.321∗	−0.287, ns	3	HP18	
ANTXR1	2	6	-	No	0.535∗∗∗	−0.003, ns	3	HP18	
FN1	1	2	APP, MAPT, APOE	NA	NA	NA	3	HP03, HP18	
ITGAV	2	0	APP, PSEN1, MAPT, APOE	Yes	−0.091, ns	0.419, ns	3	HP03, HP18	
GP1BA	2	32	APP	NA	NA	NA	3	HP18	
IL13RA1	1	NA	-	Yes	0.727∗∗∗	0.588∗	3	HP18	
TGFB1	1	0	APP, MAPT, APOE	NA	0.871∗∗∗	0.681∗∗	3	HP18	
EGFR	2	2	PSEN1, PSEN2, APOE	NA	0.119, ns	0.012, ns	3	HP18	
F11R	1	0	NA	Yes	0.801∗∗∗	0.639∗	3	HP18	
VWF	2	0	APP, MAPT, APOE	NA	0.640∗∗∗	0.493, ns	3	HP18	
FZD2	9	0	APP, PSEN1, MAPT, APOE	NA	−0.419∗∗	−0.702∗∗	3	HP18	
COL4A5	4	NA	APOE	No	0.396∗∗	0.359, ns	3	HP03	
COL16A1	1	0	-	Yes	0.470∗∗	0.176, ns	3	HP03	
FGF1	0	17	PSEN2, APOE	Yes	0.145, ns	−0.415, ns	3	HP03	
FGF18	10	1	PSEN2, MAPT, APOE	NA	NA	NA	3	HP03	
CCL4	2	0	-	Yes	0.873∗∗∗	0.887∗∗∗	3	HP03	
CD70	1	0	-	Yes	0.714∗∗∗	0.691∗∗	3	HP03	
WNT7A	1	0	PSEN2	Yes	−0.055, ns	0.241, ns	3	HP03	
TNFRSF17	2	0	PSEN2	NA	0.775∗∗∗	0.574∗	3	HP03	
Notes: eQTL, expression of the target gene is regulated by AD genetic variants (genetic variants: IGAP GWAS P-value < 1E-3; eQTL: P-value < 1E-3). GWAS, IGAP P-value < 1E-3. PPI, the target gene has significant physical interaction with APP, PSEN1, PSEN2, APOE, or MAPT (P-value <0.05). Early_DEG, the target gene is differentially expressed in AD mouse models before AD pathology emergence. Pathology cor (abeta), correlation of target gene expression with AD pathology in abeta line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). Pathology cor (tau), correlation of target gene expression with AD pathology in tau line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). CFG is the total convergent functional genomic score of the target gene; 1 CFG point is added if any of the above evidence is significant, and the CFG point ranges from 0 to 5.

Next, we further validated our HP L-R interaction results using the additional spatial transcriptome datasets of two other types of AD model mice. L-R interaction networks of HP showed unique features in different phenotypes (WT, 3xAD, and 3xPB). Compared with 3xAD and 3xPB mice, WT spots had more significant interactions within DG and CA1 subregions and cross CA1/DG (Fig. 5A, B). On the other hand, WT spots had less significant interactions within CA2-3 and cross CA1/CA2-3 compared with spots of 3xAD and 3xPB mice (Fig. 5A, B). Dot plots showed the top 10 unique L-R pairs (by P-values) that were only significantly differentially expressed in one region pair within a certain phenotype group (Fig. S3). Next, all unique L-R pairs in 3xAD and 3xPB (11–13 months) HP spots were compared with the first dataset (AppNL-G-F knock-in mice) in L1/L2 and L3/L4 HP spots at 18 months respectively. Although 3xAD and 3xPB mice were younger than the AppNL-G-F knock-in mice, they were all in the cognitive impairment disease stage with deficits in memory and learning. “COL1A1_a10b1 complex”, “COL8A2_a10b1 complex”, “NRG3_ERBB4”, “MADCAM1_a4b7 complex”, and “EPHA3_EFNA3” were significantly differentially expressed in the same region pairs in 3xAD and 3xPB spots. Meanwhile, these L-R pairs were significantly differentially expressed in L1/L2 subregion pairs at 18 months which belong to the same region that the unique L-R pairs in 3xAD and 3xPB spots showed different expression (Fig. 5C). “NAMPT_P2RY6”, “COL13A1_a1b1 complex”, “MERTK_GAS6", “NPR1_NPPC”, and “NRG3_ERBB4” were significantly differentially expressed in the same region pairs in 3xAD and 3xPB spots. Meanwhile, these LR pairs were significantly differentially expressed in L3/L4 subregion pairs at 18 months which belong to the same regions that the unique L-R pairs in 3xAD and 3xPB spots showed different expression (Fig. 5D). These results showed the potential functional interactions in the related region pairs and further validated the spatial conservation of these L-R pairs in different AD models.Figure 5 L-R interactions cross the hippocampus in 3xAD and 3xPB mice. (A) L-R interaction networks in WT, 3xAD, and 3xPB mice. The width of the edge represents the number of significantly differentially expressed L-R pairs. The edge direction is from ligand to receptor. (B) Heatmap of significantly differentially expressed L-R pair numbers between two hippocampal layers. Different colors represented the numbers of differentially expressed L-R pairs from less (blue) to more (red). (C) Significantly differentially expressed L-R interactions in the same region pairs in 3xAD, 3xPB slides, and corresponding subregion pairs in the L1/L2 group at 18 months. (D) Significantly differentially expressed L-R interactions in the same region pairs in 3xAD, 3xPB slides, and corresponding subregion pairs in the L3/L4 group at 18 months.

Figure 5

TF regulatory analysis identified a functional TF network in the HP

To further explore other kinds of interactions in HP, we used SCENIC to perform TF regulatory analysis. According to the TF regulatory network activity heatmap, different phenotype groups at 3 months and 6 months were clustered together, whereas the transcriptomics regulatory profile showed a clear separation between WT and AD groups at 12 months and 18 months (Fig. 6A). This separation was concordant with the progression of AD pathology over these periods. According to the TF regulatory activity heatmap of different HP regions, HPs (somatic layer of the hippocampus), and HPd (dendritic layer of the hippocampus) spots were clustered as two groups as expected (Fig. 6B). Regulons were clustered into three groups according to their activities. Cluster 1 regulons were mainly active in the somatic layer, whereas cluster 2 regulons were mainly active in the DG_po subregion of the dendritic layer. Genes in cluster 1 regulons were involved in neurotransmitter (regulation of neurotransmitter levels, neurotransmitter transport, neurotransmitter secretion, etc.), synaptic (synaptic vesicle cycle, synapse organization, synaptic vesicle transport, etc.), and cognition/memory related functions (Fig. 6C). Genes in cluster 2 regulons were involved in stimulus (especially visual stimulus) response related functions like visual perception, detection of light stimulus, etc. (Fig. 6D). Further, we used Cytoscape for network visualization and marked the categories of gene function, i.e., cognition/memory, neurotransmitter, and synaptic in cluster 1 regulon network (Fig. 6E) and stimulus response in cluster 2 regulon network (Fig. 6F).Figure 6 Transcription factor regulatory analysis in the hippocampus. (A) Regulons active heatmap in different ages and phenotypes. (B) Regulon active heatmap in different hippocampal layers. (C) GO analysis results in cluster 1 regulons in Figure 7B. (D) GO analysis results in cluster 2 regulons in Figure 7B. (E) Transcription factor regulatory network in cluster 1 regulons. (F) Transcription factor regulatory network in cluster 2 regulons.

Figure 6

In the TF regulatory heatmap of different Aβ index groups, we found a cluster (cluster A) that included 6 up-regulated regulons from WT-to-L4 Aβ levels (Fig. 7A). GO analysis showed that their functions were mainly involved in “negative regulation of humoral immune response” and “DNA-binding transcription activator activity, RNA polymerase II-specific” (Fig. 7B). Different gene functions were also marked in the regulatory network of cluster A (Fig. 7C). PRKCB and NPTX2 in HP network and NR2F2 in cluster A demonstrated the highest CFG score, suggesting these genes may be crucial in AD related pathological process (Table 3).Figure 7 Transcription factor regulatory analysis in Aβ index group. (A) Regulon active heatmap in different Aβ index groups in the hippocampus. (B) GO analysis results in cluster A regulons in Figure 8A. (C) Transcription factor regulatory network in cluster A regulons. (D) Regulon active heatmap in different Aβ index groups in the entorhinal cortex. (E) GO analysis results in cluster 3 regulons in Figure 7D. (F) Transcription factor regulatory network in cluster 3 regulons.

Figure 7

Table 3 CFG scores of target genes in hippocampal transcription factor analysis.

Table 3Gene	eQTL	GWAS	PPI	Early_DEG	abeta	tau	CFG	Cluster	
PRKCB	1	5	PSEN1, MAPT, APOE	Yes	−0.231, ns	−0.199, ns	4	HPs	
NR2F2	0	1	PSEN2, APOE	Yes	−0.487∗∗∗	0.075, ns	4	A	
NPTX2	1	1	-	Yes	−0.688∗∗∗	−0.783∗∗∗	4	HPs	
SYT1	0	7	MAPT	Yes	−0.209, ns	−0.471, ns	3	HPs	
SV2A	2	0	NA	Yes	−0.199, ns	−0.746∗∗	3	HPs	
STAT5A	1	0	PSEN2	NA	0.477∗∗	0.606∗	3	A, HPd	
SREBF2	1	2	APOE	NA	NA	NA	3	HPs	
HSF1	1	0	-	Yes	0.462∗∗	0.120, ns	3	HPs	
HOXA5	1	NA	-	Yes	0.372∗	0.052, ns	3	HPs	
HOMER1	1	0	-	Yes	−0.338∗	−0.341, ns	3	HPs	
GRIA1	3	30	MAPT	No	0.036, ns	0.145, ns	3	HPs	
GAD1	3	0	-	Yes	−0.207, ns	−0.562∗	3	HPs	
ESRRA	2	0	PSEN2, MAPT, APOE	Yes	−0.043, ns	−0.309, ns	3	A, HPs	
EGR1	0	2	PSEN2, MAPT	NA	−0.403∗∗	−0.559∗	3	HPs	
ARC	2	1	-	Yes	−0.194, ns	−0.259, ns	3	HPs	
Notes: eQTL, expression of the target gene is regulated by AD genetic variants (genetic variants: IGAP GWAS P-value < 1E-3; eQTL: P-value < 1E-3). GWAS, IGAP P-value < 1E-3. PPI, the target gene has significant physical interaction with APP, PSEN1, PSEN2, APOE, or MAPT (P-value <0.05). Early_DEG, the target gene is differentially expressed in AD mouse models before AD pathology emergence. Pathology cor (abeta), correlation of target gene expression with AD pathology in abeta line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). Pathology cor (tau), correlation of target gene expression with AD pathology in tau line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). CFG is the total convergent functional genomic score of the target gene; 1 CFG point is added if any of the above evidence is significant, and the CFG point ranges from 0 to 5.

TF regulatory analysis identified functional TF network in ENTI

Since ENTI is an important affected region in the early AD process,10,11 we further performed SCENIC analysis in this brain area. In the TF regulatory heatmap of different Aβ index groups, we found cluster 3 with 32 regulons was mainly active in L1 and L2 groups (Fig. 7D). This result was concordant with the progression of AD pathology in this region. GO analysis showed functions of regulons in cluster 3 were mainly involved in neurotransmitter (regulation of neurotransmitter levels, neurotransmitter transport, etc.), synaptic (synapse organization, synaptic vesicle cycle, synaptic vesicle localization, etc.), and cognition/memory related terms (Fig. 7E). Regulatory network of cluster 3 was marked with the corresponding gene functions (Fig. 7F). MEF2C and PRKCB ranked first in the CFG score list of cluster 3 TFs (Table 4). These findings suggest that MEF2C and PRKCB may potentially regulate the neuropathology process in ENTI.Table 4 CFG scores of target genes in entorhinal cortex transcription factor analysis.

Table 4Gene	eQTL	GWAS	PPI	Early_DEG	abeta	tau	CFG	
MEF2C	1	2	PSEN2, MAPT	Yes	−0.159, ns	−0.198, ns	4	
PRKCB	1	5	PSEN1, MAPT, APOE	Yes	−0.231, ns	−0.199, ns	4	
SV2A	2	0	NA	Yes	−0.199, ns	−0.746∗∗	3	
VGF	2	0	-	Yes	−0.792∗∗∗	−0.760∗∗	3	
HSF1	1	0	-	Yes	0.462∗∗	0.120, ns	3	
SREBF2	1	2	APOE	NA	NA	NA	3	
EGR1	0	2	PSEN2, MAPT	NA	−0.403∗∗	−0.559∗	3	
STAT5A	1	0	PSEN2	NA	0.477∗∗	0.606∗	3	
CUX2	2	2	PSEN2	NA	NA	NA	3	
GRIA1	3	30	MAPT	No	0.036, ns	0.145, ns	3	
Notes: eQTL, expression of the target gene is regulated by AD genetic variants (genetic variants: IGAP GWAS P-value < 1E-3; eQTL: P-value < 1E-3). GWAS, IGAP P-value < 1E-3. PPI, the target gene has significant physical interaction with APP, PSEN1, PSEN2, APOE, or MAPT (P < 0.05). Early_DEG, the target gene is differentially expressed in AD mouse models before AD pathology emergence. Pathology cor (abeta), correlation of target gene expression with AD pathology in abeta line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). Pathology cor (tau), correlation of target gene expression with AD pathology in tau line AD mouse models (r, P-value; ns: P > 0.05; ∗P < 0.05; ∗∗P < 0.01; ∗∗∗P < 0.001). CFG is the total convergent functional genomic score of the target gene; 1 CFG point is added if any of the above evidence is significant, and the CFG point ranges from 0 to 5.

To determine if any of the high CFG score genes related to ENTI were differentially expressed in human AD datasets, we obtained a human snRNA-seq dataset composed of AD and healthy ENTI cortex samples and probed for our genes.15 Six cell types were identified in this dataset (Fig. 8A, left). The UMAP plot of MEF2C expression is shown in Fig. 8A (right). Compared with other cell types, MEF2C was highly expressed in microglia and neuron (false discovery rate <0.05; Fig. 8B). A subset of AD neuron (Fig. 8C) and endothelial cell (Fig. 8D) also appeared to have up-regulated MEF2C compared with control. Presumably, genes related to ENTI with higher CFG scores may impact AD cell biology, especially in neuron to a greater extent.Figure 8 snRNA-seq in human ENTI and SSN analysis results. (A) UMAP manifold of snRNA-seq colored by cell type (left) and MEF2C expression (right). (B)MEF2C expression in different cell types. (C) UMAP of neuron subset colored by phenotype (left) and MEF2C expression (right). (D) UMAP of endothelial cell subset colored by phenotype (left) and MEF2C expression (right). (E) Genes with significantly different NDM degrees (WT vs. AD) in the hippocampus. (F) Genes with significantly different NDM degrees (WT vs. AD) in ENTI. snRNA-seq, single-nuclei RNA sequencing; ENTI, entorhinal cortex; SSN, spot-specific network; WT, wild-type; AD, Alzheimer's disease; NDM, network degree matrix.

Figure 8

SSN analysis revealed specific high NDM degree genes in different phenotypes

To explore gene connections at the whole transcriptome level, we constructed SSN for each single spot (i.e., one network for one spot) in HP and ENTI regions. Gene expression changes dynamically or even drastically with time and conditions, however, the functions or linear/nonlinear associations among genes generally remain unchanged or change gradually with small perturbations.17 So, we used SSN to transform the “unstable” gene expression matrix into a “stable” gene association based NDM. We found genes with significant differences in NDM degrees between WT and AD in HP (Fig. 8E) and ENTI (Fig. 8F) in different age groups. For example, Ttr showed a higher NDM degree in WT spots than AD spots in the 3-month group of HP. Conversely, Ttr showed a higher NDM degree in AD spots than WT spots in the 12-month group of HP. Trem2 showed a higher NDM degree in WT spots than AD spots in the 18-month group of HP. Conversely, Trem2 showed a higher NDM degree in AD spots than WT in the 18-month group of ENTI. These results showed the specific high NDM degree genes that may have wide interactions across the whole transcriptomes in certain phenotypes in different brain regions and age groups.

Discussion

In the current study, we performed L-R communication, TF regulatory network, and SSN analyses to reveal the dependence and dynamics of gene associations from multiple perspectives based on samples from mouse and human brains. Using an ST dataset of App knock-in mouse model,12 we revealed the specific L-R interaction status across the whole brain areas or HP layers within and across spots in different spatial regions and pathological status. We also further revealed important TF regulatory and SSN gene associations in HP and ENTI, respectively. Next, we used another independent ST dataset from different AD mouse models and human ENTI snRNA-seq data and AlzData database (that includes both human and mouse data) to partially validate our results. This study provides a framework for the comprehensive analyses of gene-gene associations or transcriptional networks that are spatio-temporal dependent within or among spots in ST, thereby laying the foundation for exploring characteristic gene functions and in-depth AD etiology from the perspective of spatio-temporal dependent gene interactions.

In our study, Mfuzz assigned neuron, axon, synapse, and amyloid related intra-spot L-R pairs into different clusters. Interestingly, we found that 17 L-R pairs showed opposite tendencies in the 3-month and 18-month groups through the Aβ accumulation process. The opposite tendencies of these 17 L-R pairs may be related to their different functions in different pathological statuses. For example, interactions in the protein-protein interaction network were most concentrated among EFNA/EFNB and EPH gene families. EPH receptors belong to the tyrosine kinase receptor gene family that plays a pivotal role in brain development.19 EPHA4 activation can lead to dendritic spine retraction and synaptic dysfunction.19 Several reports also support the protective role of EPHA4 in AD pathogenesis.20,21 EPHA4 promotes the proliferation of neural progenitor cells and alleviates cognitive impairment in a PDGFRβ dependent manner.20 These studies showed the contradictory functions of EPHA4 in AD pathogenesis. Since our results also showed opposite tendencies in EPHA4 related L-R pairs through Aβ accumulation in the 3-month and 18-month groups, further investigation may focus on the potentially different functions of EPHA4 in the different stages or processes of the disease during different age groups.

C3_ITGAM/ITGB2 complex was the only Aβ related L-R pair between spots in the GO enrichment analysis results that has the opposite tendencies between 3-month-old and 18-month-old mice. Interaction of complement component iC3b with its receptor complement receptor 3 (CD11b/CD18) on the surface of microglia involved in the uptake and clearance of Aβ.22 C3 deficiency in APP/PS1 mice increased plaque burden in the mouse brain with decreased age- and plaque-related synapse and neuron loss, glial reactivity, and spared cognitive decline.18 Our results showed the potential of the C3_ITGAM/ITGB2 complex's different functions in plaque-related neurodegeneration and cognitive health at different age groups.

L-R associations across brain regions showed different states between different ages, region pairs, and phenotypes in both expression levels and different numbers of significantly differentially expressed L-R pairs. For example, VCAM1_a4b7 complex in HPd_CA1_so to HPs_CA1_sp, L3/L4 level in 18 months. The low levels of VCAM1 on brain endothelial cells allowed microglia to remain ramified, calm, and proceed to neurogenesis. Increased VCAM1 expression in the old brain tethered more leukocytes to the luminal side and activated microglia to suppress neurogenesis which could be reversed by systemic anti-VCAM1 antibody.23 In addition, some L-R pairs showed spatial conservative differential expression (differentially expressed in the same region pair when compared with other region pairs) in different AD models (AppNL-G-F, 3xAD, and 3xPB) including collagen related genes like COL1A1, COL8A2, and COL13A1. Post-mortem evidences were reported in transgenic animal models and brains of individuals with AD, that increased collagen content was related to microvascular morphology changes along with the cognitive status worsening in AD mouse model and AD human patients.1,24 Our results supplied the spatio-temporal information for the interactions that contained these AD related genes. These results may furnish a foundation for future functional research based on L-R interactions instead of only individual genes, with the novel high spatial resolution empowered by ST.

According to the TF regulatory network analysis results of HP in different age groups, transcriptomics regulatory profiles of WT mice at 12 and 18 months of age were clustered together, whereas the AD mice transcriptomics regulatory profiles still changed between these two time points. These tendencies were in concordance with the progression of pathology over these periods, as well as the original research results of this dataset.12 GO enrichment analysis results showed that regulons active in the somatic layer were involved in neurotransmitter, synaptic, and cognition/memory functions, whereas regulons active in the DG_po subregion of the dendritic layer were involved in stimulus (e.g., visual stimulus) response related functions. These results showed region-specific active patterns of such regulons in HP with different potential functions.

Next, we explored the regulon active patterns through Aβ deposition in HP and found that a gradually up-regulated regulon cluster (cluster A) was involved in “negative regulation of humoral immune response” and “DNA-binding transcription activator activity, RNA polymerase II-specific” GO terms. These dynamic changes showed cluster A regulons may be related to the humoral immune response and the catastrophic failure of the transport of transcription-related molecules (especially RNA polymerase II) between the cytoplasm and the nucleus in the late AD neurotoxic process.25 GO analysis results of 32 active regulons in L1 and L2 groups of ENTI, an early affected area in AD, revealed their wide influences at the early AD stage including neurotransmitter, synaptic, and cognition/memory related functions. MEF2C was one of the high CFG score TFs in these active regulons. MEF2C was highly expressed in AD neuron and endothelial cell clusters compared with control neuron and endothelial cell clusters in the snRNA-seq data of human ENTI (Fig. 8C, D). The MEF2 TF family modulates the structural and synaptic plasticity underlying memory formation, and MEF2C gene is related to synaptic function alteration in AD.26, 27, 28, 29 Overexpression of MEF2 has been reported to impede the increase of learning-induced spine and impair memory formation via an Arc-mediated reduction of GluA2-AMPA-type glutamate receptor expression.30,31

Gene connection analysis at the whole transcriptome level showed genes with high NDM degrees were specific in certain phenotypes of different regions and age groups. Ttr is an AD neuroprotection gene that is known to bind to Aβ and facilitate its clearance from the brain,32 SSN analysis showed the transcriptome-based influences of Ttr were unique in certain phenotypes in different age groups, which underlying the Aβ clearance function of Ttr may also depend on age conditions. Trem2 is an AD risk factor due to its impacts on central/peripheral lipid metabolism and the integrity of the blood-brain barrier.33 SSN analysis showed the transcriptome-based influences of Trem2 were unique in certain phenotypes in different brain areas. These results showed these genes' unique transcriptome-based influences in different phenotypes, regions, and age groups that further suggested their known impact on AD may be unique under such different conditions. Here, this added information was explored at the novel ST resolution for the first time.

CellPhoneDB and SCENIC were developed for the single-cell RNA-seq data. Both methods have been used in the bulk-like ST data before.34, 35, 36 However, the related results were limited to the spot resolution, without specified cell type information.

In the case of the validation in the HP 3xAD/3xPB datasets, we tried to identify the L-R interactions that were significantly differentially expressed in the same region pairs in 3xAD, 3xPB slides, and the same corresponding subregion pairs in L1/L2 or L3/L4 (Fig. 4C) group at 18 months of the AppNL-G-F dataset, respectively. These results primarily validated that the significantly differentially expressed L-R interactions were consistent in different AD mouse models. In the case of the validation in the human ENTI dataset, we primarily validated that some important TFs in the AppNL-G-F dataset were also differentially expressed in certain cell types of human AD patients compared with control (Fig. 8A–D). However, a limitation of the present study is the low resolution of the ST dataset. The research object of this study is the Aβ niche (based on spots). Both cell proportion and regulon activity of certain cell types can affect the regulon activity level of spots. Under the current resolution of this dataset, it is more suitable for conducting preliminary research in the Aβ niche (100 μm in diameter), rather than excessive distinguishing between the two reasons mentioned above. A higher-resolution ST dataset and related molecular experiments will provide more thoroughness and more solid validation of the current results.

Conclusions

The current study revealed gene associations in multi-aspects based on an ST dataset of App knock-in mouse model with the validation on the Alzdata database and supplementary spatial sequencing/snRNA-seq datasets from mice and human brains. The unique genome-wide spatio-temporal dependent association features in different spatial, temporal, and pathological statuses establish the foundations to reveal the advanced and in-depth AD etiology and mechanisms in different stages of the neurodegeneration processes on the amyloid plaque niche in AD.

Author contributions

SW, HS, and HD conceived and designed the study. SW analyzed data and wrote the manuscript. JG, CQ, RS, MP, YG, HS, HX, and HD revised the manuscript. HX contributed to funding acquisition. All authors read and approved the manuscript.

Conflict of interests

The authors declare no conflict of interests.

Funding

HMX was partially supported by grants from the 10.13039/501100001809 National Natural Science Foundation of China (No. 81902277 ) and the National Key R&D Program of China (No. 2017YFC1001100 ).

Data availability

The whole brain ST dataset can be accessed with accession number GSE152506 in GEO database. The HP ST dataset was available at Mendeley https://doi.org/10.17632/6s959w2zyr.1y. The ENTI snRNA-seq dataset can be accessed with accession number GSE138852 in GEO database.

Appendix A Supplementary data

The following are the Supplementary data to this article.Table S1

The cutoff values of Aβ accumulation and the number of tissue domains per quantile.

Table S1

Table S2

Complex information of all L-R pairs involved in the analysis.

Table S2

Table S3

L-R cluster information in the 3-month and 18-month groups.

Table S3

Figure S1 The details of cross-region L-R analysis. (A) Expression level of each L-R pair. P-values <0.05 were considered statistically significant. (B) The output results of total significant L-R pair numbers in each region pair corresponded to Fig. 4A, B.

Figure S1

Figure S2 Enrichment analysis results of different L-R pair clusters. (A) L-R pairs enriched in neuron/axon/synapse/amyloid related terms in the 3-month group. These were three L-R pair clusters. The dot color represents four kinds of GO terms. Dot size is proportional to the membership score. (B) L-R pairs enriched in neuron/axon/synapse/amyloid related terms in the 18-month group. There were four L-R pair clusters. The dot color represents four kinds of GO terms. Dot size is proportional to the membership score.

Figure S2

Figure S3 Unique top 10 L-R pairs (by P-values) that only significantly different expressed in one hippocampal region in WT, 3xAD, and 3xPB mice, marked by dots with black outlines.

Figure S3

Peer review under responsibility of Chongqing Medical University.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.gendis.2024.101337.
==== Refs
References

1 Alvarez-López M.J. Castro-Freire M. Cosín-Tomás M. Long-term exercise modulates hippocampal gene expression in senescent female mice J Alzheimers Dis 33 4 2013 1177 1190 23168450
2 Navarro J.F. Croteau D.L. Jurek A. Spatial transcriptomics reveals genes associated with dysregulated mitochondrial functions and stress signaling in alzheimer disease iScience 23 10 2020 101556
3 Hohman T.J. Bush W.S. Jiang L. Discovery of gene-gene interactions across multiple independent data sets of late onset Alzheimer disease from the Alzheimer Disease Genetics Consortium Neurobiol Aging 38 2016 141 150 26827652
4 Quan X. Liang H. Chen Y. Qin Q. Wei Y. Liang Z. Related network and differential expression analyses identify nuclear genes and pathways in the hippocampus of Alzheimer disease Med Sci Monit 26 2020 e919311
5 Efremova M. Vento-Tormo M. Teichmann S.A. Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes Nat Protoc 15 4 2020 1484 1506 32103204
6 Aibar S. González-Blas C.B. Moerman T. SCENIC: single-cell regulatory network inference and clustering Nat Methods 14 11 2017 1083 1086 28991892
7 Vyas Y. Montgomery J.M. Cheyne J.E. Hippocampal deficits in amyloid-β-related rodent models of Alzheimer's disease Front Neurosci 14 2020 266 32317913
8 Mu Y. Gage F.H. Adult hippocampal neurogenesis and its role in Alzheimer's disease Mol Neurodegener 6 2011 85 22192775
9 Arendt T. Brückner M.K. Morawski M. Jäger C. Gertz H.J. Early neurone loss in Alzheimer's disease: cortical or subcortical? Acta Neuropathol Commun 3 2015 10 25853173
10 Criscuolo C. Fontebasso V. Middei S. Entorhinal Cortex dysfunction can be rescued by inhibition of microglial RAGE in an Alzheimer's disease mouse model Sci Rep 7 2017 42370
11 van Hoesen G.W. Hyman B.T. Damasio A.R. Entorhinal cortex pathology in Alzheimer's disease Hippocampus 1 1 1991 1 8 1669339
12 Chen W.T. Lu A. Craessaerts K. Spatial transcriptomics and in situ sequencing to study Alzheimer's disease Cell 182 4 2020 976 991.e19 32702314
13 Xu M. Zhang D.F. Luo R. A systematic integrated analysis of brain expression profiles reveals YAP1 and other prioritized hub genes as important upstream regulators in Alzheimer's disease Alzheimers Dement 14 2 2018 215 229 28923553
14 Ståhl P.L. Salmén F. Vickovic S. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics Science 353 6294 2016 78 82 27365449
15 Grubman A. Chew G. Ouyang J.F. A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation Nat Neurosci 22 12 2019 2087 2097 31768052
16 Kumar L. Futschik M.E. Mfuzz: a software package for soft clustering of microarray data Bioinformation 2 1 2007 5 7 18084642
17 Dai H. Li L. Zeng T. Chen L. Cell-specific network constructed by single-cell RNA sequencing data Nucleic Acids Res 47 11 2019 e62
18 Shi Q. Chowdhury S. Ma R. Complement C3 deficiency protects against neurodegeneration in aged plaque-rich APP/PS1 mice Sci Transl Med 9 392 2017 eaaf6295
19 Ganguly D. Thomas J.A. Ali A. Kumar R. Mechanistic and therapeutic implications of EphA-4 receptor tyrosine kinase in the pathogenesis of Alzheimer's disease Eur J Neurosci 56 9 2022 5532 5546 34989046
20 Chen Q. Song H. Liu C. The interaction of EphA4 with PDGFRβ regulates proliferation and neuronal differentiation of neural progenitor cells in vitro and promotes neurogenesis in vivo Front Aging Neurosci 12 2020 7 32116646
21 Tamura K. Chiu Y.W. Shiohara A. Hori Y. Tomita T. EphA4 regulates Aβ production via BACE1 expression in neurons Faseb J 34 12 2020 16383 16396 33090569
22 Fu H. Liu B. Frost J.L. Complement component C3 and complement receptor type 3 contribute to the phagocytosis and clearance of fibrillar Aβ by microglia Glia 60 6 2012 993 1003 22438044
23 Yousef H. Czupalla C.J. Lee D. Aged blood impairs hippocampal neural precursor activity and activates microglia via brain endothelial cell VCAM1 Nat Med 25 6 2019 988 1000 31086348
24 Bourasset F. Ouellet M. Tremblay C. Reduction of the cerebrovascular volume in a transgenic mouse model of Alzheimer's disease Neuropharmacology 56 4 2009 808 813 19705573
25 Mastroeni D. Chouliaras L. Grover A. Reduced RAN expression and disrupted transport between cytoplasm and nucleus; a key event in Alzheimer's disease pathophysiology PLoS One 8 1 2013 e53349
26 Nettiksimmons J. Tranah G. Evans D.S. Yokoyama J.S. Yaffe K. Gene-based aggregate SNP associations between candidate AD genes and cognitive decline Age 38 2 2016 41 27005436
27 Allen M. Kachadoorian M. Carrasquillo M.M. Late-onset Alzheimer disease risk variants mark brain regulatory loci Neurol Genet 1 2 2015 e15
28 Ren J. Zhang S. Wang X. MEF2C ameliorates learning, memory, and molecular pathological changes in Alzheimer's disease in vivo and in vitro Acta Biochim Biophys Sin 54 1 2022 77 90 35130621
29 Tang S.S. Wang H.F. Zhang W. MEF2C rs190982 polymorphism with late-onset Alzheimer's disease in Han Chinese: a replication study and meta-analyses Oncotarget 7 26 2016 39136 39142 27276684
30 Cole C.J. Mercaldo V. Restivo L. MEF2 negatively regulates learning-induced structural plasticity and memory formation Nat Neurosci 15 9 2012 1255 1264 22885849
31 Carmichael R.E. Wilkinson K.A. Craig T.J. Ashby M.C. Henley J.M. MEF2A regulates mGluR-dependent AMPA receptor trafficking independently of Arc/Arg3.1 Sci Rep 8 1 2018 5263 29588465
32 Gião T. Saavedra J. Cotrina E. Undiscovered roles for transthyretin: from a transporter protein to a new therapeutic target for Alzheimer's disease Int J Mol Sci 21 6 2020 2075 32197355
33 Li R.Y. Qin Q. Yang H.C. TREM2 in the pathogenesis of AD: a lipid metabolism regulator and potential metabolic therapeutic target Mol Neurodegener 17 1 2022 40 35658903
34 Wu R. Guo W. Qiu X. Comprehensive analysis of spatial architecture in primary liver cancer Sci Adv 7 51 2021 eabg3750
35 Wang M. Hu Q. Lv T. High-resolution 3D spatiotemporal transcriptomic maps of developing Drosophila embryos and larvae Dev Cell 57 10 2022 1271 1283.e4 35512700
36 Chen A. Liao S. Cheng M. Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays Cell 185 10 2022 1777 1792.e21 35512705
