Identification of key transcriptome biomarkers based on a vital gene module associated with pathological changes in Alzheimer’s disease

WGCNA and functional enrichment analysis of gene modules

To determine if all samples in GSE29378 were appropriate for network analysis, we investigated the sample dendrogram. The results indicated that all 63 samples mainly yielded two clusters, where GSM726085, GSM726092, GSM726098, GSM726100, GSM726101, GSM726102, GSM726106, GSM726108, GSM726113, GSM726114, GSM726115, GSM726120, GSM726121, GSM726122, GSM726125, GSM726126, GSM726127, GSM726128, GSM726129, GSM726130, GSM726133, GSM726138, GSM726140, and GSM726141 became one group while the remaining samples formed the other. And there was no potential outlier in all 63 samples in GSE29378 (Supplementary Figure 4A).

We then performed a network topology analysis to explore the appropriate soft-thresholding power for WGCNA. The results suggested that when the power value was set as 7 (scale-free R² = 0.85), the network possessed scale-free topology with integral modular features (Supplementary Figure 3A). By setting the cut height as 0.3, 12 modules were eventually identified (non-clustering genes were shown in gray) (Supplementary Figure 3B). Details of the 12 modules were provided in Supplementary Table 1.

To confirm the independence of each module, we conducted an interaction relation analysis. 1000 genes were randomly extracted from the 12 modules to design a network heatmap. As shown in the network heatmap (Supplementary Figure 4B), each block’s color represents the overlap degree of two corresponding genes on the horizontal axis and vertical axis. Dark color represents low overlap, and progressively lighter yellow color indicates higher overlap. The results indicated that genes from different modules had low overlap degrees, while blocks along the diagonal showed lighter colors, which indicated that each module was independent of the others. Hierarchical clustering dendrogram of the 12 module eigengenes (MEs) revealed two main clusters: one contained 3 modules (tan, blue, and yellow modules), while the other contained 9 modules (black, magenta, brown, green, purple, green-yellow, turquoise, pink, and red modules; Supplementary Figure 4C).

KEGG pathway enrichment analysis was conducted to investigate each module’s potential functions (Supplementary Figure 5). The results indicated that the “protein digestion and absorption” pathway was enriched in the black module. Genes in the blue module and brown module were enriched in the “neurodegeneration” and “valine, leucine and isoleucine degradation” pathways, respectively. “Glycine, serine and threonine metabolism” and “focal adhesion” were severally enriched in the green module and green-yellow module, while genes in the magenta module and pink module were respectively associated with the “tuberculosis” and “ribosome” pathways. The “Huntington’s disease”, “protein processing in endoplasmic reticulum”, and “neuroactive ligand-receptor interaction” pathways were respectively enriched in the purple, red, and tan modules. The turquoise module and yellow module were singly enriched in the “herpes simplex virus 1 infection” and “axon guidance” pathways. In summary, we identified 12 gene modules by performing WGCNA on dataset GSE29378 and parts of them were enriched in functional pathways associated with AD.

Identification of the target gene module

To identify the target module from the 12 resulting gene modules, we performed a module-trait analysis, and the analytical results were visualized in a heatmap, where the green module was detected positively correlated with both plaque score (Aβ deposition) (r = 0.43, p = 4×10⁻⁴) and BRAAK stage (p-tau accumulation) (r = 0.39, p = 0.002) (Figure 1A). In addition, the correlation between the MM value and GS value of each gene in the green module was analyzed. The results indicated that genes in the green module showed strong relationships with both plaque score (r = 0.57, p =1.7e−134) and BRAAK stage (r = 0.47, p = 3.7e−86) (Figure 1B, 1C). Thus, the green module was chosen as the target module for further study.

Figure 1. Identification of the target module. (A) Heatmap of the correlations between MEs and clinical traits of AD. Each cell contains the correlation coefficient and p-value (within brackets). (B, C) Scatter plots of GS value of plaque score (B) and BRAAK stage (C) vs. the MM value of each gene in the target module. (D) Heatmap of genes in the target module of entorhinal cortex samples in GSE118553. (E) Heatmap of genes in the target module of hippocampus samples in GSE48350. AsymAD, asymptomatic AD.

We further displayed expression heatmaps of genes in the target module based on two independent AD-related datasets (GSE118553 and GSE48350). Each heatmap’s vertical axis represents “gene symbol”, and the horizontal axis means “sample identity”. The results showed that genes in the target module were capable of discriminating AD and control individuals (Figure 1D, 1E).

Visualization of expression pattern and chromosome location of genes from the target module

We selected 100 genes with the highest MM value to explore the chromosome distribution and overall expression pattern of the target module (Figure 2). Among the 100 genes, the top 10 genes possessing the highest MM value, including GNA13, CYP2U1, PLSCR4, CAT, GJA1, PRDX1, CCDC109B, RNASE4, CYBRD1, and RAB23, were located on chromosomes 17, 4, 3, 11, 6, 1, 14, and 2.

Figure 2. Circular visualization of connectivity, expression patterns, and chromosomal positions of the 100 genes with the highest MM value in the target module. The expression profiles of control (outer ring) and AD (inner ring) individuals of GSE29378 were presented in the circular heatmap. “Red” indicates upregulation, “blue” represents downregulation, and “white” denotes genes that are not present in a given dataset. The outer circle represents chromosomes; lines coming from each gene point to their specific chromosomal locations. The ten genes with the highest MM value were shown in green font and they were connected with green lines in the center of the circles.

Identification of hub genes

To screen out hub genes in the target module, the criteria for selection were set as follows: MM value > 0.9, GS value to plaque score > 0.3, and GS value to BRAAK stage > 0.3. Finally, 6 hub genes were picked out, including GNA13, CYP2U1, PLSCR4, CAT, GJA1, and PRDX1. The specific information of the 6 genes was provided in Table 1.

Refinement of hub genes

To refine hub genes, we firstly explored the expression levels of those 6 hub genes in GSE118553 and GSE48350 (Figure 3A, 3B). Results showed that all hub genes were differently expressed in AD and control groups except for CAT. Further refinement was conducted through the receiver operating characteristic (ROC) curve analysis based on GSE118553 and GSE48350 (Figure 3C). The top 3 genes with the highest area under curve (AUC) value in each ROC analysis were picked out, and the intersecting genes were selected. As a result, GNA13 and GJA1 were identified as refined hub genes of the target module (Figure 3D).

Figure 3. Identification of the refined hub genes. (A, B) Expression of the hub genes in GSE118553 (A) and GSE48350 (B). (C) ROC analysis of the 6 hub genes in GSE118553 and GSE48350. (D) The intersection of the top 3 genes with the highest AUC value in ROC analysis based on GSE118553 and GSE48350. Data were presented as the mean ± SD in each group. ns, p > 0.05, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Functional sub-module analysis of the target module

We employed ClusterONE, a plug-in of Cytoscape, to identify functional sub-modules in the target module. Firstly, genes in the target module were employed to construct a protein-protein interaction (PPI) network using the Search Tool for the Retrieval of Interacting Genes (STRING). In the PPI network, 6693 pairs and 1147 genes were represented by edges and nodes, and the node size reflected the MM value of each gene. Six functional sub-modules containing 375 related genes were finally excavated, which were displayed in different colors in the PPI network (Figure 4A).

Figure 4. Identification of the functional sub-modules and construction of a ceRNA network. (A) The PPI network containing 6 functional sub-modules in the target module. Red nodes indicate genes in functional sub-module 1; “orange” indicates genes in functional sub- module 2; “yellow” indicates genes in functional sub-module 3; “green” indicates genes in functional sub-module 4; “bule” indicates genes in functional sub-module 5; “indigo” indicates genes in functional sub-module 6; “purple” indicates genes shared by over 2 functional sub-modules; and “gray” indicates genes not belonging to any functional sub-module. The size of certain node reflects the corresponding gene’s MM value. (B) A ceRNA network describing lncRNA–miRNA–functional sub-module interaction. “Indigo” indicates functional sub-module; “orange” indicates lncRNA; and “gray” indicates miRNA. The size of certain node reflects its degree value in the network.

Then gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses were conducted to annotate genes in functional sub-modules. Significant results with the most gene-ratio value were shown in Supplementary Figure 6. When the quantity of KEGG pathways or that of terms of each GO namespace, including biological processes (BP), molecular functions (MF), and cellular components (CC), was greater than 5, only the top 5 enriched results were displayed. GO enrichment analysis showed that genes in functional sub-module 1 (M1) were significantly enriched in “nuclear division”, “spindle”, and “protein serine/threonine phosphatase activity” terms. Genes in functional sub-module 2 (M2) were mainly enriched in “fatty acid metabolic process”, “mitochondrial matrix”, and “coenzyme binding”. The “neutrophil degranulation”, “secretory granule lumen”, and “coenzyme binding” terms were enriched in functional sub-module 3 (M3). Genes in functional sub-module 4 (M4) were enriched in “histone modification”, “transcription regulator complex”, and “transcription coactivator activity”. Genes in functional sub-module 5 (M5) were mainly enriched in “RNA splicing”, “spliceosomal complex”, and “structural constituent of the nuclear pore”. For each GO namespace, genes in functional sub-module 6 (M6) were enriched in “NF-κB signaling”, “external side of plasma membrane”, and “cysteine-type peptidase activity”. KEGG pathway enrichment analysis showed that genes in the “cell cycle” pathway were enriched in M1. Genes in M2 and M3 were enriched in the “PPAR signaling pathway” and “carbon metabolism”, respectively. The “lysine degradation” and “spliceosomes” pathways were severally enriched in M4 and M5, while genes in M6 were associated with the “NOD-like receptor signaling pathway”. Notably, these processes and pathways have been previously reported to be associated with various aspects of AD, including neuronal loss, accumulation of Aβ and p-tau, microglia dysfunction, and neurodegeneration [26–31].

Construction of a ceRNA network and identification of hub ncRNAs

To select hub ncRNAs, we constructed a lncRNA-miRNA-functional sub-module ceRNA network. The interaction analysis of lncRNA-functional sub-module and miRNA-functional sub-module predicted 84 lncRNAs and 254 miRNAs strictly participating in 96 and 316 pairs of interplays within all 6 functional sub-modules. After a comprehensive analysis of interactions among functional sub-modules and screened ncRNAs, 5 functional sub-modules, 7 lncRNAs, and 148 miRNAs were selected for the ceRNA network construction (Figure 4B). In the ceRNA network, we found that one lncRNA, MEG3, and 22 miRNAs, which were listed in Table 2, could regulate more functional sub-modules than other ncRNAs, and they were determined as hub ncRNAs in our study.

Refinement of hub ncRNAs

To refine hub ncRNAs, we performed ROC analyses in three datasets, where GSE118553 and GSE29378 were employed for lncRNA refinement, and GSE120584 was used for miRNA refinement. The results showed that the AUC value of MEG3 was greater than 0.85 based on GSE118553, which indicated its potential role in AD diagnosis (Figure 5C). As for miRNA, the top 3 miRNAs with the highest AUC value, including miR-106a-3p, miR-1-3p, and miR-24-3p, were filtered from 22 hub miRNAs for further identification (Table 2 and Figure 5D). We then displayed the expression of MEG3, miR-106a-3p, miR-1-3p, and miR-24-3p in AD and control groups based on the three datasets. The results suggested that the expression of MEG3 was significantly decreased in AD patients (Figure 5A). And miR-106a-3p and miR-24-3p were differently expressed between AD and control groups, while miR-1-3p had no expression difference in the two groups (Figure 5B). Finally, MEG3, miR-106a-3p, and miR-24-3p were selected as refined hub ncRNAs in our study.

Figure 5. Identification of the refined hub ncRNAs. (A) Expression of MEG3 in GSE118553 and GSE29378. (B) Expression of hub miRNAs in GSE120584. (C) ROC analysis of MEG3 in GSE118553 and GSE29378. (D) ROC analysis of the hub miRNAs in GSE120584. Data were presented as the mean ± SD in each group. ns, p > 0.05, ** p < 0.01 and **** p < 0.0001.

Pathological changes in 5XFAD mice

To confirm the pathological characterization of 5XFAD and LM mice employed in this study, we detected the expression of Aβ₄₂-related β-C-terminal fragment (β-CTF) and p-tau (Ser396) in hippocampus homogenates using western blotting assay. The results indicated that 4 mo 5XFAD mice had higher levels of β-CTF and p-tau (Ser396) than age-matched LM mice (p < 0.05) (Supplementary Figure 2A–2C). Besides, 5XFAD mice exhibited additively increased production of β-CTF and p-tau (Ser396) with aging (p < 0.05) (Supplementary Figure 2D–2F).

Western blotting analysis of GNA13, ROCK2, and GJA1

The western blotting analysis was applied to confirm the expression difference of GNA13, ROCK2, and GJA1 between AD and control groups. Consequently, compared with LM mice, GNA13, ROCK2, and GJA1 were all upregulated in the hippocampus of 5XFAD mice (p < 0.05) (Figure 6A, 6B).

Figure 6. Expression and functional analyses of GNA13, ROCK2 and GJA1. (A, B) Representative immunoblots (A) and densitometry (B) analysis of GNA13, ROCK2, and GJA1 expression in the hippocampus of 5XFAD mice and LM mice. (C, D) Representative immunoblots (C) and densitometry (D) analysis of GNA13, ROCK2, and GJA1 expression in the hippocampus of 5XFAD mice of various ages. (E) Top 3 gene sets (according to NES) enriched in the high-expression group of GNA13. (F) Top 3 gene sets (according to NES) enriched in the high-expression group of GJA1. Data were presented as the mean ± SD of five mice in each group. * p < 0.05, ** p < 0.01 and **** p < 0.0001.

To further verify the association between GNA13, ROCK2, GJA1, and aging in 5XFAD mice. We examined the expression of GNA13, ROCK2, and GJA1 in 1 mo, 4 mo, and 8 mo 5XFAD mice. As expected, the results indicated that compared with 1 mo 5XFAD mice, the expression of GNA13, ROCK2, and GJA1 was higher in the 4 mo and 8 mo groups, while 8 mo 5XFAD mice had the highest expression of all three genes (p < 0.05), which suggested that GNA13, ROCK2, and GJA1 were significantly upregulated in 5XFAD mice with the increase of age (Figure 6C, 6D).

Quantitative real-time PCR (qRT-PCR) analysis of MEG3, miR-106a-3p, and miR-24-3p

QRT-PCR assay was applied to confirm the expression difference of MEG3, miR-106a-3p, and miR-24-3p between AD and control groups. The results indicated that compared with LM mice, MEG3, miR-106a-3p, and miR-24-3p were all downregulated in the hippocampus of 5XFAD mice (p < 0.05) (Figure 7A–7C).

Figure 7. Expression and functional analyses of MEG3, miR-106a-3p and miR-24-3p. (A–C) The expression of MEG3, miR-106a-3p, and miR-24-3p in the hippocampus of 5XFAD mice and LM mice determined by qRT-PCR assay. (D–F) The expression of MEG3, miR-106a-3p, and miR-24-3p in the hippocampus of 5XFAD mice of various ages determined by qRT-PCR assay. (G–I) KEGG pathway analysis of target genes of MEG3, miR-106a-3p, and miR-24-3p. Data were presented as the mean ± SD of five mice in each group. * p < 0.05, ** p < 0.01 and **** p < 0.0001.

We further explored the association between MEG3, miR-106a-3p, miR-24-3p, and aging of 5XFAD mice. The results indicated that MEG3 was decreased in 4 mo and 8 mo groups (p < 0.05), compared to 1 mo group, but there was no expression difference between 4 mo group and 8 mo group (p > 0.05) (Figure 7D). MiR-106a-3p and miR-24-3p were both downregulated in 8 mo group, compared to 1 mo group (p < 0.05), however, there was no difference between 4 mo group and the other two groups (p > 0.05) (Figure 7E, 7F). In summary, the above results indicated that the expression of MEG3, miR-106a-3p, and miR-24-3p was partly associated with the aging of 5XFAD mice.

Gene set enrichment analysis (GSEA) of refined hub genes

To further investigate the potential functions of GNA13 and GJA1 in AD, GSEA was performed based on GSE29378. There were 12 and 3 significant gene sets enriched in high-expression groups of GNA13 and GJA1, respectively (false discovery rate (FDR) q value < 0.25) (Supplementary Table 3). The top 3 gene sets with the highest normalized enrichment score (NES) were thought to reflect the potential functions of GNA13 and GJA1. The results indicated that the “proteasome”, “spliceosome”, and “ribosome” gene sets were enriched in the high-expression group of GNA13, and the “basal transcription factors”, “spliceosome”, and “ribosome” gene sets were enriched in the high-expression group of GJA1 (Figure 6E, 6F). Details of results were provided in Supplementary Table 3.

KEGG pathway enrichment analysis of target genes of refined hub ncRNAs

We performed KEGG pathway enrichment analysis on downstream targets of the three refined hub ncRNAs to explore their functions. 4338 target genes of MEG3 were recognized based on the RAID (version 2.0) and LncRNA2Target databases (version 2.0). The target genes of miR-106a-3p and miR-24-3p were predicted based on the mirCode and mirTarBase databases, and there were 8473 and 11555 predicted targets of miR-106a-3p and miR-24-3p, respectively (data not shown). Enrichment analysis indicated that the most target genes of MEG3 were enriched in the “MAPK signaling pathway”, “proteoglycans in cancer”, and “cAMP signaling pathway”, in which the gene-ratio of the “MAPK signaling pathway” reached up to 0.059 (Figure 7G). The most targets of both miR-106a-3p and miR-24-3p were enriched in the “PI3K-Akt signaling pathway”, “human papillomavirus infection”, and “focal adhesion”. The miR-106a-3p related gene-ratio of the “PI3K-Akt signaling pathway” was 0.059, while the miR-24-3p related gene-ratio was 0.054 (Figure 7H, 7I).

Data collection and preprocessing

Microarray datasets, including GSE29378, GSE118553, GSE48350, and GSE120584, were collected from the GEO database. Details of each microarray dataset were provided in Table 3. GSE29378 was employed for the identification of target module, validation of MEG3, and GSEA. GSE118553 was used for validation of hub genes and MEG3. GSE48350 and GSE120584 were employed for validation of hub genes and hub miRNAs, respectively.

All expression profiles were normalized by the normalizeBetweenArrays function in the limma package in R for batch effects management [47]. In the process of mapping genes to symbols, if numerous probes were mapped to a similar symbol, their mean value was regarded as the gene manifestation value. Probes with more than one gene and empty probes were removed. The workflow used in this study was shown in Figure 8.

Figure 8. Workflow used in this study.

WGCNA

Firstly, the goodSamplesGenes function in the WGCNA package were employed to detect samples with missing values. Legitimate samples were then delivered to the cluster analysis by the hclust function in the stats package to identify and eliminate the outliers. The WGCNA package was used to construct the co-expression network [25]. The pickSoftThreshold function was utilized to find the soft-threshold power β in accordance with standard scale-free networks. The constructed adjacency matrix was then transformed into the topological overlap matrix (TOM) with the soft-threshold power β set as 7 (scale-free R² = 0.85). In addition, the dynamic tree cut method was employed, with the cut height set as 0.3 and minimal module size set as 200, to identify gene modules.

Identification of the target module and hub genes

MEs, defined as the primary critical constituents of certain gene modules, were used to evaluate the possible relationship of gene modules with clinical traits, including plaque score and BRAAK stage. The gene module with the highest correlation index was identified as the target module. The expression heatmaps of the target module in the two independent datasets, GSE118553 and GSE48350, were displayed using the ggplot2 package in R.

The hub genes were selected based on the evaluation of MM value and GS value. MM value represents the distance from the expression profile of a gene to that of the ME, which computes how adjacent a gene is to a gene module. GS value represents the log10 transformation of the p-value (GS = lg p) in the linear regression between gene expression and AD pathological traits. In this study, hub genes were defined as genes in the target module with MM value greater than 0.9 and GS values to both plaque score and BRAAK stage greater than 0.3 [48].

Visualization of the gene expression pattern and chromosome location

The OmicCircos package in R was utilized to visualize the expression pattern and chromosomal location of the 100 genes with the highest MM value in the target module.

PPI network construction and functional sub-module recognition

STRING (version 11.0) online software (https://string-db.org/) was used to search the interrelationships of genes from the target module for the construction of the interaction network [49]. Cytoscape was applied for PPI network visualization [50], in which node represented protein and edge represented interaction between proteins.

Furthermore, the network was analyzed with ClusterONE, which could contribute to the classification of proteins (genes) in the network, to identify the overlapped proteins based on a cohesion algorithm. In the ClusterONE algorithm, the higher the cohesion score between two proteins is, the more likely they could produce interaction and form a protein complex [51].

Construction of the ceRNA network

The interactions between certain lncRNA and its target genes were downloaded from RAID (version 2.0, http://rna-society.org/raid/) and LncRNA2Target database (version 2.0, http://123.59.132.21/lncrna2target/) [52, 53], the interactions between certain miRNA and its target genes were downloaded from mirCode (http://www.mircode.org) [54], and mirTarBase database (https://bio.tools/mirtarbase) [55]. And the interactions of lncRNAs and miRNAs were downloaded from mirCode database. We defined a lncRNA or miRNA as a regulator of certain functional sub-module if the genes in the functional sub-module significantly overrepresented the target genes of the lncRNA or miRNA (p < 0.01, hypergeometric test). Finally, the ceRNA network was constructed and visualized in Cytoscape.

Validation of the hub genes and ncRNAs

Independent gene expression profiles containing AD hippocampus (GSE48350) and entorhinal cortex (GSE118553) samples were employed for diagnostic validation of hub genes and hub ncRNAs. ROC analyses were carried out to evaluate a certain hub gene’s sensitivity and specificity using the pROC package [56]. GSE29378 and GSE118553 were employed for the ROC analysis of hub lncRNAs. GSE120584 was employed for ROC analysis of hub miRNAs. The result of a certain ROC analysis was reflected by the corresponding AUC value.

Animals

All progeny mice employed in this study were obtained by breeding parental male 5XFAD mice with parental female C57BL/6xSJL mice. The parental male 5XFAD mice, as described by Oakley et al., originated from Jackson Laboratory (Bar Harbor, ME, USA, Stock Number: 006554), were obtained from the agent JOINN Laboratories (Suzhou, China) [21]. The parental female C57BL/6xSJL mice were obtained from Sippr-BK laboratory animal Co. Ltd (Shanghai, China).

Genotyping of progeny mice was performed by PCR analysis of tail sample DNA, according to the supplier’s protocol. The primer details were presented in Supplementary Table 2. The results of genotyping were provided in Supplementary Figure 1, where the progeny mice with the expression of APP (377 bp) and PS1 (608 bp) were identified as 5XFAD mice, while ones with the expression of only reference DNA (324 bp) were employed as LM mice in this study.

All mice were housed in individually ventilated cages with specific pathogen-free conditions and had access to water and food pellets ad libitum. All animal experiments were in compliance with the relevant animal ethics regulations of the Animal Use Committee of Shanghai General Hospital, Shanghai Jiao Tong University School of Medicine.

Western blotting assay

Mice were sacrificed under sevoflurane anesthesia, and the hippocampus tissues were collected and homogenized in RIPA lysis buffer (Cat# WB3100, NCM, China) with protease and phosphatase inhibitors (Cat# P002, NCM, China) for 30 min. The supernatants were harvested after centrifugation at 15000 rpm for 15 min at 4° C. The concentrations of the supernatants were measured using the bicinchoninic acid protein assay kit (Cat# NCI3225CH, Thermo, USA). After quantification, equal amounts of protein were loaded in each well and separated with 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The separated proteins were electrophoretically transferred onto a polyvinylidene difluoride membrane (PVDF). The membranes were then blocked with 5% skim nonfat milk (Cat# A600669, Sangon, China) or 3% bovine serum albumin (Cat# A500023, Sangon, China) for 2h at room temperature (RT). After blocking, the membranes were incubated with primary antibody, including rabbit anti-GNA13 (44 KD, 1:1000, Cat# ab128900, Abcam, UK), rabbit anti-ROCK2 (160 KD, 1:1000, Cat# 9029S, CST, USA), rabbit anti-GJA1 (43 KD, 1:1000, Cat# ab235585, Abcam, UK), rabbit anti-Aβ₄₂ (raised against a peptide corresponding to amino acids 707-713 of P05067, 17 KD, 1:1000, Cat# 700254, Invitrogen, USA), rabbit anti-p-tau (Ser396) (79 KD, 1:1000, Cat# ab109390, Abcam, UK), rabbit anti-GAPDH (36 KD, 1:2000, Cat# AF7021, Affinity, China), and rabbit anti-Tubulin-α (50 KD, 1:2000, Cat# AF7010, Affinity, China), at 4° C overnight. The membranes were then incubated with anti-rabbit secondary antibodies (1:4000, Cat# 111-035-003, Jackson, USA) at RT for 1 h. The protein bands were detected with enhanced chemiluminescence (ECL) western blot kit (Cat# P10100, NCM, China) and visualized using a ChemiDoc XRS System with Image Lab software (Bio-Rad, USA).

QRT-PCR assay

Total RNA was extracted from hippocampus tissues by TRIzol (Cat# 15596026, Invitrogen, USA) following the manual. The RNA quality was evaluated via the A260/A280 ratio. For quantification of MEG3, cDNA was synthesized from total RNA by reverse transcriptase using the random primer (Cat# RR037A, Takara, Japan). For miR-106a-3p and miR-24-3p, RNA was reversely transcribed to cDNA using the corresponding miRNA stem-loop RT primer (Supplementary Table 4). Each sample was prepared in triplicate for a total reaction volume of 20 μl, with 250 nM forward and reverse primers, 10 μl SYBR Green (Cat# Q204, NovaBio, China), and 20 ng cDNA.

Details of the primer sequence were provided in Supplementary Table 5. All reactions were carried out in a QuantStudio 3 Real-Time PCR system. The expression of MEG3 was normalized to GAPDH, while the expression of miR-106a-3p and miR-24-3p was normalized to U6.

GSEA

GSEA software (version 4.0.3) was used to perform GSEA [57]. Based on the median expression of certain hub genes, 31 AD samples were divided into high-expression and low-expression groups. FDR q value < 0.25 was regarded as statistically significant. The reference gene set, “c2.cp.kegg.v7.1.symbols.gmt”, was downloaded from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp).

Functional enrichment analysis

GO and KEGG pathway enrichment analyses were conducted with a criterion of adjust-p-value < 0.05, using the clusterProfiler package in R [58].

Data analysis

All data were expressed as the mean with standard deviation (mean ± SD). The statistical analysis was performed using Prism (Version 8.0) Software. The results of western blotting and qRT-PCR analyses were analyzed based on the student 2-tailed unpaired t-test and one way analysis of variance (ANOVA), followed by Tukey post hoc test. Statistical significance was considered to occur at p < 0.05.