Data were obtained from the Gene Expression Omnibus (GEO) database. ¹ The GSE26168 dataset (platform GPL6883) served as the training set, comprising 9 T2DM whole blood samples and 8 healthy controls (7 impaired fasting glucose samples were excluded) [17]. The GSE95849 dataset (platform GPL22448) was used as the validation set, including 6 T2DM peripheral blood samples and 6 healthy controls (6 diabetic peripheral neuropathy samples were excluded) [18]. The GSE18309 dataset (platform GPL570) served as the MCI training set, containing 3 MCI whole blood samples and 3 healthy controls (3 AD samples were excluded) [19]. All datasets were downloaded from GEO as preprocessed Series Matrix files (background-corrected and normalized). Data were checked for distribution and log2-transformed if needed; probes with missing values were removed, and the probe with the highest mean expression was retained per gene. The analyses were performed separately for each dataset without cross-platform merging, and integration was conducted at the gene-symbol level by intersection, avoiding batch effects.

Additionally, 3,632 uric acid metabolism-related genes were retrieved from the Genecards database ² using the keyword “Uric Acid” [20] (Supplementary Table 1). To assess specificity, a biologically constrained sensitivity analysis was performed by intersecting the GeneCards gene set with KEGG purine metabolism genes and GO terms related to uric acid metabolism (e.g., inflammatory response, platelet activation, and cell adhesion), yielding 585 refined genes (Supplementary Table 1). The broad GeneCards-derived gene set was used for primary analysis to maximize sensitivity, while the constrained gene set served as a confirmatory sensitivity analysis.

Differential expression analysis was performed on the T2DM training set (GSE26168) using the limma package (v3.54.0) [21]. T2DM-DEGs were defined with thresholds of |log₂Fold Change (FC)| > 1 and p < 0.05. Volcano plots were generated with ggplot2 (v3.4.1) [22], annotating the top 10 upregulated and downregulated genes. Heatmaps were constructed using ComplexHeatmap (v2.14.0) [23] and circlize package (v0.4.15) [24] to visualize expression patterns.

Simultaneously, to obtain the DEGs between MCI and control groups in the training set GSE18309, this study performed differential analysis using the limma package (v3.54.0) on MCI and control samples. MCI-DEGs were defined with thresholds of |log₂Fold Change (FC)| > 0.5 and p < 0.05 to ensure statistical and biological significance. The ggplot2 package (v3.4.1) was used to generate a volcano plot, in which the top 10 most significantly up-regulated and top 10 most significantly down-regulated genes (sorted by p value from low to high) were annotated. Additionally, a heatmap was constructed with the ComplexHeatmap package (v2.14.0) to visualize the expression patterns of the top 10 up-regulated and top 10 down-regulated genes.

In order to obtain candidate genes associated with T2DM, MCI, and uric acid metabolism, the VennDiagram package (v1.7.1) [25] was used to intersect the sets of T2DM-DEGs, MCI-DEGs, and uric acid metabolism-related genes, and this study defined the intersection genes as candidate genes. As a sensitivity analysis, the same intersection procedure was repeated using the 585 biologically constrained uric acid-related genes to evaluate the robustness of candidate gene identification against gene set definition breadth.

Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to determine the potential biological functions and mechanisms associated with the candidate genes, employing the clusterProfiler package (v4.2.2) [26] (p < 0.05). The ggplot2 package (v3.4.1) was used to visualize the top 10 GO enrichment terms sorted by p value from low to high. Among them, GO comprised 3 functional dimensions, namely, molecular function (MF), cellular component [27], and biological process (BP). KEGG pathway enrichment revealed the classic immune pathways and cellular signaling mechanisms in which these candidate genes were involved.

Given the small candidate gene set (n = 11), the protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database (; ³ interaction score >0.15), not for hub screening, but to provide a protein-level interaction context supporting the functional coherence of the genes prior to machine learning selection. The Cytoscape package (v3.8.2) [28] was used to visualize the PPI network results. Robustness was assessed via sensitivity analysis at a medium confidence threshold (score >0.4), followed by GeneMANIA ⁴ co-expression networking to validate the functional hub status of the candidate genes.

In order to further identify candidate biomarkers associated with T2DM and MCI, this study performed the least absolute shrinkage and selection operator (LASSO) and extreme gradient boosting (XGBoost) algorithm on the GSE26168. Specifically, the glmnet package (v4.1-8) [29] was used to perform LASSO regression analysis. The variables were determined at minimum lambda (10-fold cross-validation), and genes with non-zero regression coefficients were selected as LASSO candidate characteristic genes. To maximize data utility and test feature selection robustness for our small sample (n = 17), we also used leave-one-out cross-validation (LOOCV) to assess the robustness of feature selection against different data partitioning strategies. Given the small sample size (n = 17, p = 11), standard logistic regression was not applied due to complete separation and unstable coefficient estimates observed in preliminary analyses. Instead, LASSO was preferred as it effectively reduces overfitting through L1 regularization.

Meanwhile, the xgboost package (v2.0.3.1) [30] was employed to extract the core metric Gain, which measured feature importance. And the top-ranked genes based on Gain were selected as the XGBoost candidate characteristic genes. In addition, the ggvenn package (v0.1.9) [31] was used to identify the intersection of the candidate characteristic genes from 2 algorithms and defined these overlapping genes as the characteristic genes. For methodological robustness, univariate logistic regression was conducted as a sensitivity analysis.

To evaluate the ability of characteristic genes to distinguish T2DM and control samples and avoid circularity, area under curve (AUC)-based validation of the identified characteristic genes was performed with pROC package (v1.18.0) [32] in the independent external dataset GSE95849 (n = 12), which was not used in the feature selection step. The training set GSE26168 AUC values are also reported for reference only. Subsequently, characteristic genes that exhibited AUC >0.7 in both the training and validation sets were selected as candidate biomarkers. Gene expression analysis was subsequently conducted on these candidate biomarkers by comparing T2DM and control samples from the GSE26168 and GSE95849 datasets. Based on the candidate biomarkers, the Wilcoxon test was used to compare the expression differences of the candidate biomarkers between the 2 groups (p < 0.05), and box plots were drawn using the ggplot2 package (v3.4.1). Candidate biomarkers that showed significant differences between the 2 groups and consistent expression trends were selected as biomarkers.

To further investigate the relevant signaling pathways and biological information associated with biomarkers in T2DM, this study downloaded the reference gene set (c2. cp.kegg.v7.4. symbols.gmt) from the molecular signatures database (MsigDB, ⁵ ). Based on biomarkers and all other genes in the dataset, the Spearman correlation coefficient was calculated using the psych package (v2.1.6) [33]. The genes were then sorted from largest to smallest using the correlation coefficient as the sorting criterion. Subsequently, GSEA was performed employing the clusterProfiler package (v4.2.2), with 1,000 permutations and thresholds of |Normalised enrichment score (NES)| > 1, p < 0.05. The same GSEA procedure was additionally applied to the MCI dataset (GSE18309) using HP and ITGB3 as core genes, to explore pathway enrichment patterns in the MCI context and compare them with T2DM findings.

To assess differences in immune cell infiltration between the T2DM and control samples, the GSVA package (v1.50.0) [34] was used to quantify the relative abundance of 28 immune cell types across all samples in the training set GSE26168. The infiltration abundance of immune cells was then visualized using the pheatmap package (v1.0.12) [35]. Moreover, the Wilcoxon test was conducted to identify immune cells exhibiting notable disparities between T2DM and control groups, with p-values adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR); adjusted p < 0.05 was considered significant. The ggplot2 package (v3.4.1) was used for box plot visualization. To analyze the correlations between differential immune cells, Spearman correlation analysis was performed using the psych package (v2.1.6). Furthermore, to further investigate the relationship between biomarker expression and immune cell infiltration, the psych package (v2.1.6) was employed to calculate Spearman correlation coefficient between biomarkers and the immune cell enrichment scores. Analyses were conducted with the thresholds of |correlation coefficient(cor)| > 0.3 and BH-adjusted p < 0.05.

To gain deeper insights into the correlations among biomarkers, Spearman correlation analysis was performed using the psych package (v2.1.6) based on the T2DM and control samples from the training set GSE26168. Correlations were identified with p < 0.05 and |cor| > 0.3. To assess immune microenvironment changes in MCI, the same immune infiltration analysis pipeline was independently applied to the GSE18309 dataset, with BH-FDR correction applied to all comparisons.

To investigate the distribution of biomarkers across tissues, the Human Proteome Atlas (HPA, ⁶ ) database was used to explore the expression levels of biomarkers in different tissues. Visualization of biomarker expression levels across tissues was achieved using histograms. The top 5 tissue distribution data for each biomarker were imported into Cytoscape software (v3.8.2) to generate tissue-biomarker network diagrams.

To further identify the specific subcellular localization of the biomarkers, this study used the Genecards ² database to obtain the subcellular distribution information of the biomarkers. The list of subcellular localization data with confidence = 5 was selected.

First, to explore potential therapeutic drugs, this study conducted Connectivity Map (CMap) analysis to predict small-molecule compounds for T2DM. Specifically, drug features were extracted from the Connectivity Map database. ⁷ Biomarkers served as input data. The final analysis results were assigned a score ranging from −100 to 100. Negative scores indicated the ability to reverse gene expression, signifying potential therapeutic significance. CMap profiles are mainly derived from cancer cell lines; thus, drug predictions are purely computational and hypothesis-generating, requiring validation in disease-relevant models before therapeutic conclusions can be made. Simultaneously, the biomarkers were imported into the Coremine Medical database, ⁸ where relevant Traditional Chinese Medicine (TCM) information pertaining to the biomarkers was downloaded. Potential therapeutic traditional Chinese medicines were identified at a significance level of p < 0.05.

The expression levels of these biomarkers were verified through the application of RT-qPCR. Samples of blood were collected from patients diagnosed with T2DM and healthy control subjects at Anhui Public Health Clinical Center,The First Affiliated Hospital of Anhui Medical University, North District, with each sample representing an independent biological replicate (n = 5 per group). Prior to sample collection, informed consent was duly obtained from each participant, and the study protocol had been approved by the Ethics Committee of Anhui Public Health Clinical Center(Approval No. PJ-YX2025-044). Detailed clinical characteristics of all participants, including age, sex, BMI, HbA1c, and current medications, are presented in Supplementary Table 12.

For the extraction of total RNA, the TRIzol (R401-01) reagent (produced by Novozyme Biotech Co., Ltd., located in Nanjing, China) was employed. The extraction process was meticulously conducted in strict accordance with the manufacturer - provided protocol. Subsequently, the reverse transcription of the extracted total RNA into cDNA was executed using the HP All-in-one qRT Master Mix II RT203-Ver.1 for qPCR, which was manufactured by Kunming Yungeng Biotechnology Co., Ltd. The reverse transcription procedure was carried out precisely as per the instructions provided by the supplier. The RT-qPCR reactions were performed with the utilization of the ® 2×Universal Blue SYBR Green qPCR Master Mix (supplied by Servicebio, based in Wuhan, China). The primer sequences utilized for the polymerase chain reaction amplification were comprehensively detailed in Supplementary Table 2. To exclude gDNA contamination, cDNA synthesis was performed using a kit with integrated gDNA digestion (Yeasen, Cat. No. 11141ES60). Amplification curves displayed typical sigmoidal kinetics; despite the absence of formal efficiency calculation, specificity was supported by sigmoidal profiles and single-peak melt curves, indicating no non-specific amplification or primer-dimer formation. Although SELP primers were initially designed against a predicted transcript (XM accession), NCBI Primer-BLAST confirmed the amplicon resides within a conserved coding region identical to the curated reference (NM_003005), ensuring valid mRNA-level detection. Among various candidate genes, GAPDH was specifically selected as the internal reference gene for normalization purposes. Each sample was assayed in technical triplicate, and the mean Ct value was used for subsequent calculations. The relative expression levels of the biomarkers were calculated by implementing the 2^−ΔΔCT method, as previously described and validated in a published study [36]. After the completion of data acquisition from the RT-qPCR experiments, the resultant data were initially exported to Microsoft Excel for preliminary organization. Subsequently, the data underwent in - depth analysis and visualization using GraphPad Prism 10.1.2 software ⁹ (p < 0.05), which is a widely recognized and highly utilized tool in the field of scientific data analysis and graphical representation.

Serum concentrations of HP, ITGβ1, and SELP in T2DM patients (n = 11) and healthy controls (n = 10) were measured using commercially available enzyme-linked immunosorbent assay (ELISA) kits (Fenxibio, Beijing, China): HP (Cat# FXs202379), ITGβ1 (Cat# FXs207165), and SELP (Cat# FXs205897). All assays were performed in accordance with the manufacturer’s instructions.

All bioinformatic analyses were performed in the R language (v4.2.2). Data from different cohorts were compared using the Wilcoxon test. The Spearman method was utilized for correlation analysis. Differences were regarded as significant at p < 0.05. Post hoc power analysis was performed based on qPCR expression data using a two-sample t-test framework (two-sided, α = 0.05). Effect sizes (Cohen’s d) were calculated from the mean and standard deviation of each group.

Principal component analysis (PCA) confirmed acceptable within-dataset sample distributions across all three datasets, with no outlier samples detected (Supplementary Figure 1). Differential expression analysis identified 2,434 T2DM-DEGs, with 1,220 upregulated and 1,214 downregulated genes in T2DM versus controls. The top differentially expressed genes included APOL4, GRIN2B, and TRPC5 (downregulated), as well as CES4 and IL6R (upregulated) (Figure 1A). The heatmap similarly demonstrated expression patterns of these genes in T2DM and control groups (Figure 1B), indicating differential activation or suppression in T2DM. Similarly, a total of 655 MCI-DEGs was identified between MCI and control groups. Compared to the control samples, there were 350 up-regulated genes (such as FOLR3, OLFML1, and RPS11P1) and 305 down-regulated genes (such as LINC00339 and LOC158434) in the MCI samples (Figures 1C,D). Finally, a total of 11 intersection genes between the 2,434 T2DM-DEGs, 655 MCI-DEGs, and 3,632 uric acid metabolism-related genes were identified as candidate genes (Figure 1E). As a sensitivity analysis using the biologically constrained 585-gene set, intersection with T2DM-DEGs and MCI-DEGs yielded five genes (FSTL1, HP, ITGB3, SELP, and THBD), all overlapping with the original 11 candidates (Supplementary Figure 2A).

The functional enrichment analysis of these candidate genes revealed associations with 177 GO terms (p < 0.05). Among these, the BPs mainly involved the processes of platelet activation and blood coagulation. The CCs primarily included the components of platelet alpha granule membrane, external side of plasma membrane, and haptoglobin-hemoglobin complex. In terms of MFs, the candidate genes were mainly enriched in the functions of heparin binding and organic acid binding (Figure 1F; Supplementary Table 3). Additionally, 18 KEGG pathways were identified, including Fluid shear stress and atherosclerosis, Neutrophil extracellular trap formation, and Rap1 signaling pathway (p < 0.05) (Figure 1G; Supplementary Table 4). These results suggested that the candidate genes might be involved in immune system, signal transduction and metabolism of cofactors and vitamins. The PPI network (11 nodes, 14 edges) provided protein-level interaction context (Figure 1H); sensitivity analysis at a higher threshold (>0.4) and GeneMANIA co-expression analysis consistently identified HP, SELP, ITGB3, THBD, and FSTL1 as functionally coherent hubs (Supplementary Figures 2B,2C), supporting the biological relevance of the candidate gene set. The results showed these genes might influence each other directly or indirectly, and might participate in the biological processes linking the uric acid metabolism to T2DM and MCI.

In the LASSO regression analysis, the optimal value of lambda. min = 0.03627 and 5 LASSO candidate characteristic genes including HP, ITGB3, TAT, FSTL1, and SELP were ultimately identified (Figure 2A). Notably, the LOOCV feature selection results were identical to those obtained from the original 10-fold cross-validation (Supplementary Table 5). Through the XGBoost model and Gain-based ranking, 3 XGBoost candidate characteristic genes (HP, ITGB3, and SELP) were selected from the candidate genes (Figure 2B). Subsequently, this study intersected the genes obtained from the mentioned above machine learning algorithms and finally found that 3 genes (HP, ITGB3, and SELP) could be indicated by all the algorithms, meaning that these genes could be used as characteristic genes (Figure 2C). As a sensitivity analysis, univariate logistic regression across 11 candidate genes (Supplementary Table 6) indicated that HP, ITGB3, and SELP showed the largest effect sizes. Stepwise logistic regression selected HP and TAT, but complete separation resulted in unstable estimates. These results provide supportive evidence for the robustness of the LASSO/XGBoost-based feature selection.

To avoid circularity, biomarker performance was primarily evaluated in the independent GSE95849 dataset. HP (AUC = 1.000, 95% CI: 1.000–1.000), ITGB3 (AUC = 0.917, 95% CI: 0.759–1.000), and SELP (AUC = 0.917, 95% CI: 0.738–1.000) all demonstrated strong discriminatory ability (AUC >0.7). However, given the small sample size, the result AUC = 1.0 should be interpreted with caution as it may reflect potential overestimation of performance. The results observed in the training dataset (AUC: 0.806–0.986) are provided for reference only and showed a consistent trend (Figure 2D). Accordingly, these 3 characteristic genes were identified as candidate biomarkers. Next, the expression of HP, ITGB3, and SELP showed significant differences in both 2 datasets, and all of them showed a significant upward trend in the T2DM group (Figure 2E). Finally, these 3 candidate biomarkers were identified as biomarkers. The results reinforced the correlation between biomarkers and disease states, as well as their potential biological significance.

Given the limited sample size of the training dataset (n = 17), the following GSEA results should be interpreted as exploratory. The GSEA observed that pathways associated with the expression of HP and ITGB3 were 7 and 4, respectively (Supplementary Tables 7, 8). Interestingly, SELP was not significantly enriched in pathways under the thresholds of |NES| > 1 and p < 0.05. Specifically, HP was remarkably enriched in the pathways of ribosome, primary immunodeficiency, and type I diabetes mellitus, while ITGB3 was enriched in glycosaminoglycan biosynthesis heparan sulfate, adherens junction, and ribosome (Figures 3A,B). GSEA in the GSE18309 dataset showed that HP was enriched in proteasome, ubiquitin-mediated proteolysis, neuroactive ligand–receptor interaction, and ECM-receptor interaction pathways, while ITGB3 was enriched in ECM-receptor interaction and neuroactive ligand-receptor interaction. Both T2DM and MCI shared enrichment in vascular and neural signaling pathways, but enrichment was weaker and less extensive in MCI, consistent with its earlier disease stage (Supplementary Figures 3A,B).

Based on the samples in GSE26168, the relative composition of all cell types in each sample was shown in Figure 4A. The results showed that compared with the control group, Neutrophils constituted the highest proportion of cells in T2DM group, suggesting that the composition of immune cells in the immune microenvironment changed under the T2DM condition (p < 0.05). Following Benjamini-Hochberg FDR correction, 2 immune cell types exhibited significant differences between T2DM and control samples: immature B cells (d = 0.509, 95% CI: 0.263–0.738, FDR = 0.014) and plasmacytoid dendritic cells (d = 0.376, 95% CI: 0.186–0.457, FDR = 0.009) (Figure 4B). Correlation analysis revealed a significant positive correlation between immature B cells and plasmacytoid dendritic cells (cor = 0.66, FDR-adjusted p < 0.05) (Figure 4C). Additionally, Biomarker-immune cell correlation analysis demonstrated that HP showed significant positive correlations with plasmacytoid dendritic cells (cor = 0.67, 95%CI: 0.3291–0.8713, FDR = 0.0073) and immature B cells (cor = 0.62, 95%CI: 0.255–0.816, FDR = 0.0114); ITGB3 exhibited the strongest positive correlations with plasmacytoid dendritic cells (cor = 0.67, 95%CI: 0.353–0.819, FDR = 0.0073) and immature B cells (cor = 0.76, 95%CI: 0.432–0.902, FDR = 0.0021); SELP showed no significant correlation with either differential immune cell type (all FDR-adjusted p > 0.05) (Figure 4D).

Notably, biomarkers HP and ITGB3 showed a significant positive correlation (cor = 0.73, p < 0.05), while ITGB3 and SELP demonstrated a significant correlation (cor = 0.56, p < 0.05). However, no significant correlation was observed between HP and SELP (p > 0.05) (Figure 4E). The significant correlation between the biomarkers and the abundance of immune cells implied that HP, ITGB3, and SELP might have potential roles in regulating the immune microenvironment of T2DM.

The same immune infiltration analysis in the MCI dataset (GSE18309) showed a trend toward increased monocytes, macrophages, and neutrophils; however, no cell type remained significant after BH-FDR correction. Compared to T2DM, the weaker immune signal in MCI may reflect a less advanced stage of immune activation (Supplementary Figure 3C).

In the tissue-biomarker network, neither SELP nor ITGB3 exhibited tissue specificity, while HP demonstrated extremely high tissue specificity, being expressed exclusively in the liver (Figures 5A–C). Specifically, SELP was primarily distributed in the urinary bladder, lung, and appendix tonsil, whereas ITGB3 was expressed in tissues such as the thyroid gland, seminal vesicle, and kidney. Both SELP and ITGB3 were co-localized in smooth muscle tissue (Figure 5D).

Furthermore, subcellular localization results indicated that all biomarkers, HP, ITGB3, and SELP, were primarily localized to the extracellular (Confidence = 5). Besides, both SELP and ITGB3 are predominantly localized to the plasma membrane, consistent with their well-characterized roles as membrane-associated proteins. (Supplementary Table 9). Given the lack of UniProt support for GeneCards’ nuclear annotation of ITGB3, only its classical plasma membrane localization was considered to maintain functional and biological accuracy.

First, a total of 2,429 drug-biomarker interactions were predicted in the CMap database. Table 1 presents the top 10 and bottom 10 candidate compounds (by scores). For example, drugs such as calcifediol (vitamin D receptor agonist) and meclofenamic acid (cyclooxygenase inhibitor) exhibited high negative scores, indicating suggesting their gene expression signatures may computationally oppose the disease-associated transcriptomic pattern. TCM associations predicted by the Coremine Medical database are summarized in Supplementary Table 10 as exploratory text-mining outputs.

RT-qPCR analysis indicated that, when juxtaposed with the control group, the expression levels of HP (P = 0.0050), ITGB3 (P = 0.0038), and SELP (P = 0.0358) were markedly increased in individuals diagnosed with T2DM (Figures 6A–C). These findings were corroborated at the protein level: serum concentrations of HP, ITGβ1, and SELP were similarly elevated in T2DM patients (n = 11) compared with healthy controls (n = 10), as determined by ELISA (all P < 0.05; Figures 6D–F). These observations aligned closely with the trends unearthed by our prior bioinformatics analyses, thereby providing additional validation for the relevance and reliability of these biomarkers in the pathophysiological context of T2DM. Effect sizes (Cohen’s d) for HP (2.46), ITGB3 (2.52), and SELP (1.56) indicated substantial biological impacts (Supplementary Table 11). Based on the smallest observed effect size (d ≈ 1.6), a requirement of n = 6 per group (80% power, α = 0.05). Our study utilized n = 5 per group, aligning closely with this calculated threshold for robust statistical inference.