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

Datasets relevant to AD and T2DM were systematically curated from the Gene Expression Omnibus (GEO) database () 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)for AD (training set), consisting of mRNA expression profiles from brain tissue of 46 AD patients and 32 healthy controls; (2)for T2DM (training set): consisting of mRNA expression profiles of subcutaneous adipose tissue from 24 T2DM patients and 24 healthy controls; (3)for AD (validation set), consisting of mRNA expression profiles from brain tissue of 44 AD patients and 56 healthy controls; (4)for T2DM (validation set): mRNA expression profiles from liver tissue of 14 T2DM patients and seven healthy controls. https://www.ncbi.nlm.nih.gov/geo/↗ GSE109887 GSE104674 GSE122063 GSE64998

ARGs were obtained on October 28, 2021, from multiple databases, including HADb (), AUTOPHAGY DATABASE (), HAMdb (), and MsigDB (). For MsigDB, 19 gene sets closely associated with autophagy were selected using "autophagy" as the keyword, encompassing a total of 557 genes. Additionally, GeneCards () 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 (). http://www.autophagy.lu/↗ http://autophagy.info↗ http://hamdb.scbdd.com/↗ https://www.gsea-msigdb.org/gsea/msigdb/index.jsp↗ http://www.genecards.org/↗ Table S1

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 thedataset, and the "edgeR" package for the T2DM and control groups in thedataset. 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-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 |logFC| > 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 (., 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 (). Differential expression results were visualized using volcano plots. GSE109887 GSE104674 2 2 P e.g [Qi et al., 2023]

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 (for AD andfor 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= |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 T2DMcontrol) was calculated using Pearson's correlation coefficient, and modules significantly associated (< 0.05) with disease status were selected for further analysis. GSE109887 GSE104674 2 ij vs P

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 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 correctedvalue < 0.05 were considered statistically significant. P

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.

To validate the hub genes, we analyzed their expression in two independent datasets for AD and T2DM (and) 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. GSE122063 GSE64998

To investigate the expression of hub genes in the brains of T2DM models, 6-week-old male db/db mice (= 8) and age-matched male littermate db/m controls (= 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. n n

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).

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).was used as the internal reference gene for normalization. Relative gene expression levels were calculated using the 2method. Primer sequences are shown in. GAPDH −ΔΔCt Table 1

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 < 0.05 was considered significant. P

In this study, four GEO datasets (,,, and) were analyzed, and their information is summarized in. 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 |logFC (fold change) | > 0.3 and FDR-adjustedvalue < 0.05. Based on these criteria, 2,416 DEARGs were identified between brain tissue samples from AD patients and normal controls using thedataset (). Similarly, 1,882 DEARGs were obtained in subcutaneous adipose tissue samples from T2DM patients and normal controls using thedataset (). GSE109887 GSE104674 GSE122063 GSE64998 GSE109887 GSE104674 Table 2 Fig. 1A Fig. 1B 2 P

WGCNA was performed on AD samples from thedataset to identify gene co-expression modules associated with AD status (). 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 (). 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 Rvalues 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 (). Hierarchical clustering and module merging (height cut-off < 0.25) yielded five distinct gene modules (). Subsequent correlation analysis between module eigengenes and AD clinical traits demonstrated that the blue, brown, and green modules were significantly associated with AD status (< 0.05,). Collectively, these three modules encompassed 800 genes, which were selected for downstream analysis as potential AD-related biomarkers.provides a comprehensive visualization of the WGCNA workflow and outcomes for AD. GSE109887 Fig. 2 Fig. 2A Fig. 2B Fig. 2C Fig. 2D Figure 2 2 2 P

In the T2DM analysis, WGCNA was applied to thedataset to identify co-expression networks associated with the pathogenesis of T2DM (). 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 (and). 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 (). Correlation analysis with T2DM-associated clinical traits revealed that five modules—blue, brown, yellow, green, and red—were significantly correlated with T2DM status (< 0.05,). These modules collectively contained 913 genes, which were considered potential T2DM-related markers. The detailed outcomes are depicted in, which captures the WGCNA process for T2DM. GSE104674 Fig. 3 Figs. 3A 3B Fig. 3C Fig. 3D Fig. 3 2 P

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, a total of 33 genes () 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 Fig. 4 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

To investigate the biological roles of the 33 shared genes, GO and KEGG enrichment analyses were performed. As shown in–, 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" (). In terms of CC, the genes were mainly associated with the "microtubule" and "neuron projection cytoplasm" (). For MF, enrichment was observed in "kinase regulator activity" and "ubiquitin protein ligase binding" (). KEGG pathway analysis further indicated that these genes were predominantly enriched in the FoxO signaling pathway and pathways related to various types of cancer (). Figs. 5A 5C Fig. 5A Fig. 5B Fig. 5C Fig. 5D

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

To validate the identified hub genes, we analyzed two independent datasets:(AD) and(T2DM). To analyze differential expression patterns of hub genes (and). 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, andremained significantly upregulated, whileremained downregulated in both AD and T2DM patients. Conversely,exhibited an opposite expression pattern between the two. However,andwere significantly downregulated, and, andwere upregulated only in AD patients. In contrast,was significantly upregulated only in T2DM patients. GSE122063 GSE64998 Figs. 8A 8B ANXA5, BAG3 CDKN1A MET PFKFB3 BNIP3 FBXW7 CCL2, HSPB1 YAP1 CCND1

The diagnostic performance of five hub ARGs (,,,, and), which demonstrated significant differential expression in both the training and validation sets, was evaluated using ROC curve analysis. Among these genes,, andwere significantly upregulated in both AD and T2DM patients, whilewas significantly downregulated in both groups, suggesting that the expression of these four genes could serve as shared diagnostic biomarkers. In contrast,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-values were calculated to quantify diagnostic accuracy. The results are visualized in(training set) and(validation set). In AD patients, the area under the curve (AUC) values for,, andwere 0.681, 0.698, 0.736, 0.762, and 0.794 in the training set (–, all< 0.05) and 0.715, 0.820, 0.685, 0.884 and 0.845 in validation set (–, all< 0.05). For T2DM patients, the AUC values for,, andwere 0.901, 0.913, 0.722, 0.87, and 0.866 in the training set (–, all< 0.05) and 0.857, 0.796, 0.918, 0.867, and 0.796 in validation set (–, all< 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, onlyandconsistently showed good diagnostic performance, whereas the remaining genes exhibited moderate accuracy. ANXA5 BAG3 CDKN1A MET PFKFB3 ANXA5, BAG3 CDKN1A MET PFKFB3 P ANXA5, BAG3, CDKN1A MET PFKFB3 P P ANXA5, BAG3, CDKN1A MET PFKFB3 P P MET PFKFB3 Fig. 9 Fig. 10 Figs. 9A 9E Figs. 10A 10E Figs. 9F 9J Figs. 10F 10J

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 (,,,,) was based on their initial hub status (degree > 3 in the PPI network), consistent dysregulation across both training and validation datasets (< 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) (and). 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 (–). To further examine gene expression, we conducted qRT-PCR on brain tissues from both groups (). The results indicated thatandwere significantly upregulated in db/db mice, whereaswas significantly downregulated, relative to db/m controls (–). In contrast,andexpression did not differ significantly between the two groups (and). These findings partially diverge from bioinformatics results:was consistently upregulated in both AD and T2DM, whileexhibited 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. ANXA5 BAG3 CDKN1A MET PFKFB3 P ANXA5 BAG3 MET CDKN1A PFKFB3 CDKN1A PFKFB3 Figs. 11A 11B Figs. 12A 12E Fig. 13 Figs. 13A 13C Figs. 13D 13E