Bioinformatic analysis and experimental validation of hub autophagy-related genes as novel biomarkers for type 2 diabetes mellitus and Alzheimer’s disease

Data collection and ARGs datasets

Datasets relevant to AD and T2DM were systematically curated from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) to align with the primary objective of this study, which focuses on these two distinct diseases. The inclusion criteria for dataset selection were as follows: (1) a minimum of 15 samples in total, (2) at least five individuals in both the patient and control groups, and (3) the availability of either raw data or processed gene expression profiles. Consequently, the following GEO accession datasets were utilized: (1) GSE109887 for AD (training set), consisting of mRNA expression profiles from brain tissue of 46 AD patients and 32 healthy controls; (2) GSE104674 for T2DM (training set): consisting of mRNA expression profiles of subcutaneous adipose tissue from 24 T2DM patients and 24 healthy controls; (3) GSE122063 for AD (validation set), consisting of mRNA expression profiles from brain tissue of 44 AD patients and 56 healthy controls; (4) GSE64998 for T2DM (validation set): mRNA expression profiles from liver tissue of 14 T2DM patients and seven healthy controls.

ARGs were obtained on October 28, 2021, from multiple databases, including HADb (http://www.autophagy.lu/), AUTOPHAGY DATABASE (http://autophagy.info), HAMdb (http://hamdb.scbdd.com/), and MsigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). For MsigDB, 19 gene sets closely associated with autophagy were selected using “autophagy” as the keyword, encompassing a total of 557 genes. Additionally, GeneCards (http://www.genecards.org/) was queried using the same keyword “autophagy”, and genes with a relevance score of ≥1.0 were retained. This search yielded a total of 2,181 ARGs (Table S1).

Differential gene expression analysis

Following data preparation for each disease, differential expression analysis of autophagy-related genes (DEARGs) was performed using the “limma” R package for the AD and control groups in the GSE109887 dataset, and the “edgeR” package for the T2DM and control groups in the GSE104674 dataset. In both analyses, genes were considered significantly differentially expressed if they met the thresholds of |log₂(fold change)| > 0.3 and a false discovery rate (FDR)-adjusted P-value < 0.05. The FDR threshold (<0.05) was applied to control the expected proportion of false positives among the identified DEARGs to ≤5%, a standard and stringent correction for multiple testing in transcriptomic studies. The |log₂FC| > 0.3 threshold (corresponding to approximately 1.23-fold linear changes) was selected to: (1) filter out technical noise and focus on biologically relevant changes; (2) adopt a relatively inclusive approach to detect potentially important genes with subtle expression differences, which is particularly relevant in complex comorbid conditions; and (3) account for the fact that modest expression changes (e.g., 1.2–1.5 fold) in key regulatory genes can have substantial biological impact, especially in neurodegenerative contexts. Therefore, these thresholds were chosen to balance sensitivity and specificity, enabling the identification of subtle but consistent transcriptional changes associated with disease comorbidity, while minimizing false positives through FDR correction. Similar thresholds have been used in published transcriptomic studies of comorbid diseases (Qi et al., 2023). Differential expression results were visualized using volcano plots.

Weighted gene co-expression network analysis (WGCNA)

To identify gene modules with similar co-expression patterns and to investigate hub genes within the network and their association with disease status (T2DM or AD), we performed weighted gene co-expression network analysis (WGCNA) separately on each dataset (GSE109887 for AD and GSE104674 for T2DM) to account for tissue-specific expression profiles. The soft-thresholding power (β) was determined using the scale-free topology criterion (scale-free R² > 0.85) to ensure that the resulting adjacency matrix approximated a scale-free network. The adjacency matrix was constructed by raising the absolute value of the Pearson correlation coefficient between gene pairs to the power β (a_ij = |cor(i, j)|β). A power of β = 4 satisfied this criterion for both datasets. The topological overlap matrix (TOM) was then derived from this adjacency matrix to measure network interconnectedness. A hierarchical clustering dendrogram was generated based on TOM-based dissimilarity (1-TOM). Modules were identified using the dynamic tree cutting algorithm, and those with highly correlated module eigengenes (height < 0.25) were merged. The correlation between each module eigengene and disease status (AD or T2DM vs control) was calculated using Pearson’s correlation coefficient, and modules significantly associated (P < 0.05) with disease status were selected for further analysis.

Shared gene identification

DEARGs from each disease dataset were intersected with genes from WGCNA modules significantly correlated with the respective disease status, resulting in a set of shared ARGs potentially involved in the pathogenesis of both AD and T2DM.

Functional enrichment analysis of DEARGs

Functional enrichment analysis of the 33 shared genes for KEGG pathways and GO terms was performed using the R package “clusterProfiler”. GO enrichment analysis was performed across three categories: biological process (BP), cellular component (CC), and molecular function (MF). Pathways or terms with an FDR-adjusted corrected P value < 0.05 were considered statistically significant.

PPI network construction and differential expression analysis of hub DEARGs

A PPI network was constructed to identify key hub ARGs shared between AD and T2DM. The shared DEARGs were input into the STRING database to generate the PPI network. Hub genes within this network were identified based on nod degree, which represents the number of direct connections a node (protein) has within the network. A higher degree values indicate greater potential functional significance. Genes with a degree >3 were defined as candidate hub genes. Among them, the top 12 genes exhibiting the highest degree values were ranked and designated as key hub genes. To investigate differences in their expression, the expression levels of these key hub DEARGs were compared across normal, AD, and T2DM groups, and the results were visualized using boxplots.

External datasets validation and prediction accuracy assessment of hub genes

To validate the hub genes, we analyzed their expression in two independent datasets for AD and T2DM (GSE122063 and GSE64998) and the results were visualized using boxplots. The diagnostic performance of the hub genes was evaluated using ROC curve to differentiate AD/T2DM from normal controls and calculate the area under the curve (AUC) to evaluate their diagnostic value.

Validation of hub genes in animal models

To investigate the expression of hub genes in the brains of T2DM models, 6-week-old male db/db mice (n = 8) and age-matched male littermate db/m controls (n = 8) were purchased from the SiPeiFu (Beijing) Biotechnology Co., Ltd. Upon arrival, the mice underwent a 1-week acclimatization period followed by 5 weeks of experimental feeding, resulting in a total age of 11 weeks at the study endpoint. The mice were housed under specific-pathogen-free (SPF) conditions with a 12 h/12 h light/dark cycle, controlled temperature (20 ± 2 °C), and relative humidity (55%). Bedding was regularly replaced, and cages were routinely cleaned and disinfected. Four mice were housed in each cage to meet social needs, with environmental enrichment provided in the form of nesting materials and shelters. After 5 weeks of feeding, body weight and fasting blood glucose (FBG) levels were measured to confirm hyperglycemia in db/db mice compared to controls. FBG levels were determined through tail vein sampling. In accordance with institutional guidelines, no analgesia was administered due to the minimally invasive nature of the procedure. At the end of the study, mice were euthanized by intraperitoneal injection of entobarbital sodium (50 mg/kg body weight) and cervical dislocation. Brains were immediately harvested, dissected on an ice-cooled board, snap-frozen in liquid nitrogen, and stored at −80 °C for qRT-PCR analysis. All animal procedures were approved by the Animal Ethics Committee of Capital Medical University (AEEI-2023-296) and performed according to the relevant ethical standards.

Morris water maze test

The Morris water maze (MWM) test was conducted to evaluate spatial learning and memory abilities. Briefly, the MWM test was conducted in a circular black pool filled with water and a hidden platform. The pool was surrounded by black curtains to eliminate external visual cues and was divided into four quadrants: northeast, northwest, southeast, and southwest. During the acquisition phase, mice underwent four trials per day for 5 consecutive days. Each trial ended when the mouse located the hidden platform or after 60 s had elapsed. Escape latency (the time taken to locate the platform) and swim paths were recorded using an overhead camera. On the final day, a probe trial was conducted by removing the platform, and mice were allowed to swim for 60 s. The number of platform crossings and the time spent in the target quadrant were recorded using a video tracking system (JX business, Shanghai, China).

Quantitative qRT-PCR analysis

After the MWM test, all mice were sacrificed by cervical dislocation. The brains were immediately harvested, dissected on an ice-cooled board, snap-frozen in liquid nitrogen, and stored at −80 °C for quantitative qRT-PCR analysis. Total RNA was extracted using TRIzol reagent (Invitrogen, Waltham, MA, USA), and its purity and concentration were evaluated by measuring absorbance at 260 and 280 nm using an ultramicro UV-visible spectrophotometer. Complementary DNA (cDNA) was synthesized from the RNA using the Hifair® III 1st Strand cDNA Synthesis SuperMix (Yeason, Shanghai, China), according to the manufacturer’s instructions. The reverse transcription reaction was carried out at 55 °C for 15 min, followed by 85 °C for 5 s. Quantification of mRNA expression was performed through qRT-PCR using SYBR Green PCR kit (SAIPUBIO, Shenzhen, China) on a FQD-96A real-time PCR system (Bio-Rad, Hercules, CA, USA). All primer sequences were synthesized by Sangon Biotech (Shanghai, China). GAPDH was used as the internal reference gene for normalization. Relative gene expression levels were calculated using the 2^−ΔΔCt method. Primer sequences are shown in Table 1.

Statistical analysis

Statistical analyses were carried out using SPSS 29.0 (SPSS, Inc., Chicago, IL, USA) and GraphPad Prism 9.0 (GraphPad Software, San Diego, CA, USA). Data are presented as mean ± SD. Data were compared between two groups using the independent samples t-tests. Escape latency in the MWM test was analyzed using repeated-measures ANOVA. A two-tailed P < 0.05 was considered significant.

Identification of DEARGs

In this study, four GEO datasets (GSE109887, GSE104674, GSE122063, and GSE64998) were analyzed, and their information is summarized in Table 2. ARGs were compiled from multiple databases, including HADb, the Autophagy Database, HAMdb, MsigDB, and GeneCards, yielding a total of 2,181 ARGs. Differential expression analysis was performed using the thresholds of |log₂FC (fold change) | > 0.3 and FDR-adjusted P value < 0.05. Based on these criteria, 2,416 DEARGs were identified between brain tissue samples from AD patients and normal controls using the GSE109887 dataset (Fig. 1A). Similarly, 1,882 DEARGs were obtained in subcutaneous adipose tissue samples from T2DM patients and normal controls using the GSE104674 dataset (Fig. 1B).

Construction of co-expression modules by WGCNA in AD and T2DM

WGCNA was performed on AD samples from the GSE109887 dataset to identify gene co-expression modules associated with AD status (Fig. 2). The soft-thresholding power was determined based on the criterion of achieving a scale-free topology with a fit index (R²) > 0.85, resulting in an optimal power value of 4. With this threshold, a weighted co-expression network was constructed using a one-step approach (Fig. 2A). Most genes exhibited a low degree of connectivity, while a few genes displayed high connectivity, consistent with a scale-free network topology indicative of hub genes. The linear regression analysis yielded R² values of 0.86 and 0.84, with slopes of −1.08 and −1.52, respectively, further supporting the scale-free nature of the constructed network (Fig. 2B). Hierarchical clustering and module merging (height cut-off < 0.25) yielded five distinct gene modules (Fig. 2C). Subsequent correlation analysis between module eigengenes and AD clinical traits demonstrated that the blue, brown, and green modules were significantly associated with AD status (P < 0.05, Fig. 2D). Collectively, these three modules encompassed 800 genes, which were selected for downstream analysis as potential AD-related biomarkers. Figure 2 provides a comprehensive visualization of the WGCNA workflow and outcomes for AD.

In the T2DM analysis, WGCNA was applied to the GSE104674 dataset to identify co-expression networks associated with the pathogenesis of T2DM (Fig. 3). Consistent with the approach used for AD, the soft-thresholding power was determined based on the scale-free topology criterion (R² > 0.85), yielding an optimal power value of 4 (Figs. 3A and 3B). The network was constructed with this power, and hierarchical clustering with a height cut-off of <0.25 for module merging produced seven gene modules (Fig. 3C). Correlation analysis with T2DM-associated clinical traits revealed that five modules—blue, brown, yellow, green, and red—were significantly correlated with T2DM status (P < 0.05, Fig. 3D). These modules collectively contained 913 genes, which were considered potential T2DM-related markers. The detailed outcomes are depicted in Fig. 3, which captures the WGCNA process for T2DM.

Identification of shared genes

To investigate the autophagy-related co-pathogenesis of AD and T2DM, we conducted an intersection analysis between the previously mentioned DEARGs and the genes identified by WGCNA. As shown in Fig. 4, a total of 33 genes (ANXA5, BAG3, BNIP3, CCL2, CCND1, CDKN1A, COX5A, DNM3, FBXW7, GABARAPL1, GLS, GSTP1, HSPB1, HSPB8, KIF5B, MAP1B, MAP2K1, MET, NUAK1, PFKFB3, PHGDH, PLK2, PRKCH, RAB6B, RPN2, S100A4, TBC1D9, TNS3, TSPO, TUBB2A, TXNIP, UCHL1, YAP1) were identified. We hypothesize that these genes may play critical roles in the autophagy-related pathogenesis of both AD and T2DM and may represent a set of shared molecular targets (Fig. 4).

Function annotation of shared genes

To investigate the biological roles of the 33 shared genes, GO and KEGG enrichment analyses were performed. As shown in Figs. 5A–5C, GO enrichment analysis revealed that the shared genes were primarily enriched in autophagy-related BP, including “autophagy,” “process utilizing autophagic mechanism,” and “regulation of autophagy” (Fig. 5A). In terms of CC, the genes were mainly associated with the “microtubule” and “neuron projection cytoplasm” (Fig. 5B). For MF, enrichment was observed in “kinase regulator activity” and “ubiquitin protein ligase binding” (Fig. 5C). KEGG pathway analysis further indicated that these genes were predominantly enriched in the FoxO signaling pathway and pathways related to various types of cancer (Fig. 5D).

PPI network construction and identification of hub genes

The 33 shared genes between AD and T2DM were used to construct a PPI network via the STRING database. As indicated in Fig. 6A, the resulting network was visualized using Cytoscape. Based on node degree ranking, 12 genes (ANXA5, CCND1, MAP2K1, BAG3, BNIP3, HSPB1, CCL2, FBXW7, MET, YAP1, CDKN1A and PFKFB3) were identified as hub genes for further analysis (Fig. 6B). Moreover, the expression of these hub DEARGs were evaluated in the GSE109887 (AD) and GSE104674 (T2DM) datasets, comparing patient samples with controls (Figs. 7A and 7B). Seven hub genes, including ANXA5, BAG3, CCL2, CCND1, CDKN1A, HSPB1, and YAP1, were significantly upregulated in both the brain tissues of AD patients and the subcutaneous adipose tissues of T2DM patients. In contrast, BNIP3 and MET were significantly downregulated in both conditions, showing consistent expression patterns across AD and T2DM. However, three hub genes exhibited divergent expression trends between the two diseases. Specifically, PFKFB3 gene was significantly upregulated in the brain tissues of AD patients but markedly downregulated in the adipose tissues of T2DM patients. In contrast, FBXW7 and MAP2K1 were significantly downregulated in AD but upregulated in T2DM, demonstrating opposite expression patterns between the two conditions.

External datasets validation and diagnostic efficacy assessment of hub genes

To validate the identified hub genes, we analyzed two independent datasets: GSE122063 (AD) and GSE64998 (T2DM). To analyze differential expression patterns of hub genes (Figs. 8A and 8B). Differential expression patterns of the 12 previously identified hub genes were examined. Expression levels were compared between brain tissues of AD patients and controls, as well as between subcutaneous adipose tissues of T2DM patients and controls. The analysis revealed that ANXA5, BAG3, and CDKN1A remained significantly upregulated, while MET remained downregulated in both AD and T2DM patients. Conversely, PFKFB3 exhibited an opposite expression pattern between the two. However, BNIP3 and FBXW7 were significantly downregulated, and CCL2, HSPB1, and YAP1 were upregulated only in AD patients. In contrast, CCND1 was significantly upregulated only in T2DM patients.

The diagnostic performance of five hub ARGs (ANXA5, BAG3, CDKN1A, MET, and PFKFB3), which demonstrated significant differential expression in both the training and validation sets, was evaluated using ROC curve analysis. Among these genes, ANXA5, BAG3, and CDKN1A were significantly upregulated in both AD and T2DM patients, while MET was significantly downregulated in both groups, suggesting that the expression of these four genes could serve as shared diagnostic biomarkers. In contrast, PFKFB3 was significantly upregulated in AD patients but downregulated in T2DM patients, indicating that its expression may serve as a biomarker for diagnosing a single disease and may not be suitable for simultaneous diagnosis of both conditions. The area under the curve (AUC), 95% confidence intervals (CI), and P-values were calculated to quantify diagnostic accuracy. The results are visualized in Fig. 9 (training set) and Fig. 10 (validation set). In AD patients, the area under the curve (AUC) values for ANXA5, BAG3, CDKN1A, MET, and PFKFB3 were 0.681, 0.698, 0.736, 0.762, and 0.794 in the training set (Figs. 9A–9E, all P < 0.05) and 0.715, 0.820, 0.685, 0.884 and 0.845 in validation set (Figs. 10A–10E, all P < 0.05). For T2DM patients, the AUC values for ANXA5, BAG3, CDKN1A, MET, and PFKFB3 were 0.901, 0.913, 0.722, 0.87, and 0.866 in the training set (Figs. 9F–9J, all P < 0.05) and 0.857, 0.796, 0.918, 0.867, and 0.796 in validation set (Figs. 10F–10J, all P < 0.05). AUC values between 0.7 and 0.9 indicate good diagnostic accuracy, while values between 0.5 and 0.7 indicate moderate diagnostic accuracy. Based on this criterion, the five hub genes demonstrated excellent diagnostic performance for T2DM across both datasets (0.9 > AUC > 0.7). In AD, only MET and PFKFB3 consistently showed good diagnostic performance, whereas the remaining genes exhibited moderate accuracy.

Validation of candidate hub genes in animal models

Following the identification of 33 shared DEARGs and 12 hub genes, we applied a series of criteria to select the most robust candidate genes for further validation. The selection of five hub genes (ANXA5, BAG3, CDKN1A, MET, PFKFB3) was based on their initial hub status (degree > 3 in the PPI network), consistent dysregulation across both training and validation datasets (P < 0.05 for Wilcoxon test) and sustained diagnostic utility (AUC > 0.6 in ROC analysis) in both AD and T2DM contexts across all datasets. We next verified the expression of the five identified hub genes in brain tissue samples from T2DM animal models (db/db mice) and corresponding controls (db/m mice) (Figs. 11A and 11B). Behavioral assessments using MWM revealed that learning and memory abilities were significantly impaired in the T2DM model compared to controls, as evidenced by reduced time spent in the target quadrant and fewer platform crossings during both training and testing days (Figs. 12A–12E). To further examine gene expression, we conducted qRT-PCR on brain tissues from both groups (Fig. 13). The results indicated that ANXA5 and BAG3 were significantly upregulated in db/db mice, whereas MET was significantly downregulated, relative to db/m controls (Figs. 13A–13C). In contrast, CDKN1A and PFKFB3 expression did not differ significantly between the two groups (Figs. 13D and 13E). These findings partially diverge from bioinformatics results: CDKN1A was consistently upregulated in both AD and T2DM, while PFKFB3 exhibited disease-specific expression patterns—upregulated in AD but downregulated in T2DM. These discrepancies may be attributed to interspecies differences in autophagic regulation or the limitations of monodisease models that do not account for tissue-specific pathophysiological contexts.