This study was approved by the ethics committee of the First Affiliated Hospital of Guizhou University of Traditional Chinese Medicine (Number: KL2024-017). All participants provided written informed consent prior to the study. A total of 10 newly diagnosed patients with T2DM, aged 40–65 years (mean age: 47.23 ± 7.41 years), and 10 healthy individuals, aged 40–65 years (mean age: 47.23 ± 7.41 years), were included in this study. The diagnosis of T2DM was established based on the criteria set by the American Diabetes Association (ElSayed et al., 2023). Participants were recruited from individuals attending our hospital. Participants with acute or chronic inflammatory diseases, cardiovascular disease, uncontrolled hypertension, type 1 diabetes mellitus (T1DM), gestational diabetes, smoking habits, or alcohol consumption were excluded.

Peripheral blood mononuclear cells (PBMCs) were isolated from heparinized whole blood samples immediately using Ficoll-Hypaque density-gradient centrifugation under sterile conditions as previously conducted (Arab Sadeghabadi et al., 2022). After washing the cells twice with phosphate-buffered saline, isolated PBMCs were suspended in RPMI 1640 medium. Then PBMCs were counted, 5 × 10⁶ cells/well plated in a 6-well plate containing RPMI 1640 supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin for 24 h. Total RNA (500 ng) was extracted from PBMCs with TRIzol reagent (Nordic Bioscience, Beijing, China) and reverse into cDNA (10 uL) by reverse transcription kit (TakaRa, PrimeScript™ RT Master Mix, Cat No: RR036Q). The primer sequences of the specific genes (PPP1CA and CTSD) are listed in Table 1. Real-time quantitative PCR (RT-qPCR) for these mRNAs was performed and analyzed using cDNA and SYBR Green PCR Master Mix (Nordic Bioscience, Beijing, China, dilution: 5X). RT-qPCR was performed on a 7,500 fast real-time PCR system (Applied Biosystems, Foster City, CA, United States). The relative amounts of mRNA were determined based on 2^−ΔΔCT calculations with GAPDH as a control reference.

Gene expression profiles for training cohort (GSE15932) were extracted from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (GPL570 platform), which comprised peripheral blood samples from 8 T2DM patients and 8 controls. The GSE26168 dataset (GPL6883 platform), containing transcriptomic data from 9 T2DM patients and 8 healthy controls, was downloaded from the same database for external validation. A total of 1,519 asparagine-related genes (ARGs) were identified by merging two independent sources: 326 genes from Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) and 1,322 genes from GeneCards (https://www.genecards.org/) with relevance scores >6 (Li et al., 2024) (Supplementary Table S1).

Differentially expressed genes (DEGs) between T2DM samples and controls were identified from the training cohort (GSE15932) using the R package “limma” (v 3.54.0) (Ritchie et al., 2015) (|log₂ fold change (FC)| >0.5, P < 0.05). To visualize the distribution of DEGs, a volcano plot was created with “ggplot2” package (v 3.4.4) (Gustavsson et al., 2022). The top 10 upregulated/downregulated DEGs (ranked by |log₂FC|) were labeled on the volcano plot. Additionally, a heatmap was constructed using the R package “Complex Heatmap” (v 2.14.0) (Gu et al., 2016) to display the expression patterns of these top 20 DEGs across all samples. To obtain candidate genes associated with both asparagine and T2DM pathogenesis, the intersection of DEGs and ARGs was implemented with “Venn Diagram” package (v 1.7.3).

To shed light on the biological functions and signaling pathways associated with the candidate genes, GO and KEGG pathway enrichment analyses were executed with “cluster Profiler” package (v 4.7.1.003) (Wu et al., 2021). The GO analysis categorized gene functions into biological process (BP), cellular component (CC), and molecular function (MF). KEGG pathway analysis was conducted to detect metabolic and signal transduction pathways linked to the candidate genes. For both analyses, a significance threshold of adjusted P < 0.05 was implemented. Enrichment results were ranked by adjusted p-values in ascending order, and the top 5 most significant terms from each category were obtained for visualization.

PPI network for the candidate genes was analyzed based on STRING database (https://string-db.org/) (confidence score >0.4). The resulting interaction data were visualized with Cytoscape (v 3.9.1). To identify candidate feature genes within the PPI network, four topological algorithms from the Cytoscape plugin “Cytohubba”—Maximum Clique Centrality (MCC), Edge Percolated Component (EPC), Closeness, and Betweenness—were applied. Genes with weights above the median threshold (i.e., exceeding the median value for MCC, EPC, Closeness, and Betweenness scores) were retained. The shared genes identified by all four algorithms were determined as candidate feature genes using the R package “ggvenn” (v 0.1.9) (Gao et al., 2021).

To detect feature genes associated with asparagine in T2DM, the Boruta algorithm was implemented using the R package “Boruta” (v 8.0.0) (Zhou et al., 2023). The training cohort (GSE15932) was subjected to iterative feature importance analysis with a significance threshold of p = 0.001 and a maximum iteration parameter of maxRuns = 100. Genes classified as “Confirmed” by the Boruta algorithm—indicating statistically significant importance in distinguishing T2DM from controls—were retained for downstream analysis. A random forest classifier was constructed using “random Forest” package (v 4.7.1.1) (Alderden et al., 2018) to further refine feature selection. The model utilized 5-fold cross-validation, with the out-of-bag (OOB) error rate serving as the performance metric to determine the optimal number of decision trees (ntree). The ntree value corresponding to the lowest OOB error rate was selected to finalize the model. Genes ranked in the top 50% of mean decrease Gini importance scores were defined as feature genes. Support Vector Machine Recursive Feature Elimination (SVM-RFE) was executed employing “caret” package (v 6.0-93). The algorithm iteratively removed the least important features from the candidate gene set based on linear kernel SVM weights. A 5-fold cross-validation strategy was applied to evaluate classification accuracy, and the gene subset achieving the lowest error rate was identified. The overlapping genes across all three machine learning methods were defined as feature genes.

To validate the expression consistency of feature genes across cohorts, Wilcoxon rank-sum tests were performed separately in the training cohort (GSE15932) and validation cohort (GSE26168). Genes exhibiting statistically significant differential expression in both cohorts (P < 0.05), with consistent trends between T2DM and controls, were defined as biomarkers. To assess the diagnostic potential of these biomarkers, receiver operating characteristic (ROC) curves were generated. The area under the curve (AUC) was calculated to evaluate discriminatory ability between T2DM samples and controls, with AUC values interpreted as follows: AUC >0.7 indicated strong discriminatory potential. The ROC analysis was performed using the R package “pROC” (v 1.18.0).

To delve into biological pathways associated with the biomarker in T2DM, GSEA was implemented with R “clusterProfiler” package (v 4.7.1.003) (|normalized enrichment score (NES)| >1, adjusted P < 0.05). Spearman correlation coefficients (cor) between the identified biomarkers and all other genes in the training cohort (GSE15932) were calculated with “psych” package (v 2.1.6). Genes were ordered by the cor values. The reference gene set (h.all.v2024.1.Hs.symbols.gmt) was acquired from the MSigDB.

To quantify the infiltration levels of 28 immune cell types in peripheral blood samples from training cohort (GSE15932) (Zhang Y. et al., 2023), ssGSEA was executed employing “GSVA” package (v 1.46.0) (Hänzelmann et al., 2013). A heatmap depicting immune cell infiltration scores across T2DM samples and controls was generated with “pheatmap” package (v 1.0.12). Wilcoxon tests were used to contrast immune cell enrichment scores between T2DM samples and controls. Immune cell types with adjusted P < 0.05 were defined as differentially infiltrated immune cells. Boxplots visualizing these differences were constructed with “ggplot2” package (v 3.4.4). Cor values between differentially infiltrated immune cells were calculated (|cor| >0.3, P < 0.05). A correlation heatmap was created with “ggcorrplot” package (v 0.1.4). Spearman correlation analysis was executed using “psych” package (v 2.1.6) to assess relationships between differentially infiltrated immune cells and biomarkers. Cor values and p-values were calculated (|cor| >0.3, P < 0.05).

Protein sequences corresponding to the identified biomarkers were acquired from the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/). The subcellular localization of biomarker-associated mRNAs was predicted based on mRNALocater (http://bio-bigdata.cn/mRNALocater/), a machine learning-based database specialized in eukaryotic mRNA localization. The distribution of predicted subcellular localizations was visualized with “ggplot2” package (v 3.4.4).

Functional associations between the identified biomarkers and other genes were analyzed based on the GeneMANIA database (https://genemania.org/), which integrates heterogeneous genomic and proteomic data, including protein-protein interactions, co-expression networks, genetic interactions, and pathway associations. Biomarkers were uploaded to GeneMANIA, and the 20 genes with the highest functional similarity scores (based on weighted interaction evidence) were selected. Enriched pathways and biological functions associated with the biomarker-top 20 gene network were identified (p < 0.05, false discovery rate (FDR) <0.05). The top 7 most significant pathways and functions were retained for visualization.

Transcription factors (TFs) regulating the identified biomarkers were retrieved from ChIPBase database (http://rna.sysu.edu.cn/chipbase/). In addition, microRNA (miRNA) targeting the biomarkers was predicted using MicroCosm database (https://tools4mirs.org/software/mirna_databases/microcosm-targets/). Subsequently, long non-coding RNAs (lncRNAs) regulating identified miRNAs were predicted employing ENCORI database (https://rnasysu.com/encori/). Interactions with a number of supporting CLIP-seq experiments (clipExpNum) >20 were retained to ensure high-confidence lncRNA-miRNA pairs. A lncRNA-miRNA-mRNA regulatory network was created by integrating miRNA-mRNA and lncRNA-miRNA interaction pairs. The networks were visualized with Cytoscape software (v 3.9.1).

Potential therapeutic compounds targeting the biomarkers were detected using Comparative Toxicogenomics Database (CTD, https://ctdbase.org/). Direct and indirect drug-target interactions were retrieved (combined score >50). A network was constructed for visualization employing Cytoscape software (v 3.9.1). 3D structures of the top-ranked drugs (ranked by combined score) were downloaded in SDF format from NCBI PubChem Compound database (https://pubchem.ncbi.nlm.nih.gov/). Crystal structures of biomarker encoding proteins were obtained from UniProt database (https://www.uniprot.org/) in PDB format. Molecular docking was performed using CB-Dock2 (https://cadd.labshare.cn/cb-dock2/). The optimal binding conformation of the protein and drug was visualized using PyMOL software (v 3.1) (Seeliger and de Groot, 2010).

Bioinformatics analyses were executed utilizing the R software (v 4.2.2). Comparative analysis of data from different groups was carried out using the Wilcoxon test (P < 0.05).

Transcriptome profiling revealed 1,093 DEGs between T2DM samples and controls (|log₂FC| >0.5, P < 0.05), with 692 genes upregulated and 401 downregulated in T2DM samples (Figures 1A,B). Intersection analysis of these DEGs with 1,519 ARGs identified 90 candidate genes potentially involved in T2DM pathogenesis (Figure 1C). Functional annotation through GO enrichment analysis demonstrated significant associations across three categories. Candidate genes were predominantly enriched in some BP, such as glucose metabolic process, hexose metabolic process, and monosaccharide metabolic process (Figure 1D). CC involved vesicle-related structures such as cytoplasmic vesicle lumen, endocytic vesicle lumen, and secretory granule lumen (Figure 1E). For MF, candidate genes were significantly enriched in alpha-mannosidase activity, insulin receptor binding, and mannosidase activity (Figure 1F; Supplementary Table S2). KEGG pathway analysis identified 7 significantly enriched metabolic and disease pathways. Notably, the candidate genes were associated with Type II diabetes mellitus, diabetic cardiomyopathy, and Non−alcoholic fatty liver disease (Figure 1G; Supplementary Table S3). These findings strongly indicate that the regulation of asparagine could play a pivotal role in orchestrating the progression of T2DM through specific candidate genes.

PPI network analysis of 80 candidate genes (excluding 10 outlier proteins) revealed 813 interaction pairs, with SOD2, TGFB1, and MMP9 exhibiting the highest connectivity (Figure 2A). Subsequent integration of four algorithmic approaches identified 25 candidate feature genes (Figure 2B). Boruta algorithm further refined this set to 16 robust features, including APOA1, BRCA1, CTSD, FECH, GFM1, HBB, IGFBP3, MAPK3, PKM, PPP1CA, SOD2, TGFB1, TIMP1, and TTR (Figure 2C). Permutation feature importance with optimal ntree = 6 (minimum error rate) prioritized 5 feature genes, including CTSD, G6PD, PKM, PPP1CA, and TIMP1 (Figure 2D). The intersection of these three algorithmic outputs yielded 4 shared feature genes (PPP1CA, CTSD, PKM, TIMP1) (Figure 2E). In the training cohort (GSE15932), all four feature genes exhibited significant upregulation in T2DM samples compared to controls (P < 0.01). Subsequent validation in the independent cohort (GSE26168) revealed that PPP1CA and CTSD maintained consistent overexpression patterns in T2DM samples (P < 0.05). In contrast, TIMP1 showed nonsignificant expression trends, while PKM was excluded from analysis due to its absence in the validation cohort (Figure 2F). To further assess the discriminatory potential of PPP1CA and CTSD, ROC analysis was performed. In the training cohort (GSE15932), PPP1CA achieved an AUC of 0.969 and CTSD an AUC of 0.984, both exceeding the threshold of 0.9 indicative of strong diagnostic utility. In the validation cohort (GSE26168), PPP1CA achieved an AUC of 0.806 and CTSD an AUC of 0.875, where both genes again demonstrated AUC values >0.8 (Figure 2G). Then, mRNA expression of PPP1CA and CTSD was validated by RT-qPCR. Results indicated PPP1CA and CTSD were significantly upregulated in T2DM samples (P < 0.001, Figure 2H). These results establish PPP1CA and CTSD as robust biomarkers for T2DM.

To elucidate the biological pathways associated with PPP1CA and CTSD, GSEA was systematically performed. PPP1CA demonstrated significant enrichment in 8 hallmark pathways (|NES| >1, adjusted P < 0.05), with prominent involvement in MYC targets V1, E2F targets, G2M checkpoint, and myogenesis (Figure 3A; Supplementary Table S4). For CTSD, significant pathway enrichment was observed in myogenesis, G2M checkpoint, E2F targets, MYC targets V1 (Figure 3B; Supplementary Table S5). These results collectively highlight distinct yet complementary functional roles of PPP1CA and CTSD in T2DM progression.

Immune infiltration analysis unveiled significant alterations in immune cell composition between T2DM samples and controls. Activated dendritic cells (DCs, P < 0.01), mast cells (P < 0.05), and myeloid-derived suppressor cells (MDSCs) (P < 0.05) exhibited elevated infiltration levels in T2DM samples compared to controls, while immature DCs showed marked depletion (P < 0.001, Figures 4A,B). Correlation analysis demonstrated a strong positive association between mast cells and activated DCs (cor = 0.49, P < 0.05), whereas MDSCs inversely correlated with immature DCs (cor = −0.54, P < 0.05) (Figure 4C). Notably, PPP1CA and CTSD displayed distinct immunomodulatory associations. PPP1CA was negatively correlated with immature DCs infiltration (cor = −0.61, P < 0.05) but positively correlated with activated DCs (cor = 0.65, P < 0.01) and MDSCs (cor = 0.62, P < 0.01). CTSD showed positive correlations with activated DCs (cor = 0.58, P < 0.05) and mast cells (cor = 0.53, P < 0.05), alongside a negative correlation with immature DCs (cor = −0.54, P < 0.05) (Figure 4D; Supplementary Table S6). These findings collectively indicated systemic remodeling of the immune microenvironment in T2DM. The biomarker-immunocyte correlations suggested that PPP1CA and CTSD might orchestrate immune-metabolic crosstalk, potentially linking asparagine metabolism to T2DM progression.

Subcellular localization analysis revealed distinct compartmentalization patterns for the biomarkers. CTSD exhibited balanced distribution across cytoplasmic matrix, extracellular space, nucleus, and plasma membrane, with moderate presence in lysosomes, endoplasmic reticulum, Golgi apparatus, endosomes, cytoskeleton, and mitochondria, and minimal peroxisomal localization. In contrast, PPP1CA demonstrated predominant localization in cytoplasmic matrix, nucleus, and plasma membrane, followed by extracellular space and endoplasmic reticulum (Figure 5A). GeneMANIA network analysis revealed PPP1CA and CTSD primarily interacted with co-expressed genes through physical interactions, co-expression, predicted, and co-localization. Functional enrichment of these interaction partners demonstrated significant involvement in phosphatase-related processes, such as regulation of dephosphorylation, regulation of protein dephosphorylation, phosphatase regulator activity, protein dephosphorylation, protein serine/threonine phosphatase complex formation, protein phosphatase regulator activity, and phosphatase complex assembly (Figure 5B). PPP1CA demonstrated binding capacity with 69 TFs, while CTSD interacted with 38 TFs. Notably, RUNX1T1, MYC, and IRF1 exhibited dual regulatory control over both PPP1CA and CTSD, suggesting potential shared transcriptional mechanisms in T2DM pathophysiology (Figure 5C). MicroRNA interaction profiling identified 75 miRNA-mRNA regulatory pairs, with PPP1CA targeted by 39 miRNAs and CTSD by 36 miRNAs. The conserved miRNAs hsa-miR-885-3p, hsa-miR-675-5p, and hsa-miR-940 emerged as co-regulators of both biomarkers, indicating potential post-transcriptional coordination (Figure 5D). Exploration of competing endogenous RNA networks revealed lncRNA-mediated regulatory axes. The lncRNA AL355075.4 was predicted to simultaneously regulate both PPP1CA and CTSD through competitive binding to hsa-miR-330-5p, forming the AL355075.4-hsa-miR-330-5p-PPP1CA and AL355075.4-hsa-miR-330-5p-CTSD regulatory axes (Figure 5E). Collectively, these results delineated multi-layered regulatory networks governing PPP1CA and CTSD expression, encompassing transcriptional, post-transcriptional, and competing endogenous RNA-mediated mechanisms, which might critically influence T2DM pathogenesis.

To explore therapeutic candidates targeting the identified biomarkers, pharmacological profiling was performed. A total of 103 compounds were predicted to interact with PPP1CA and CTSD, with 72 compounds specifically targeting PPP1CA and 31 targeting CTSD. Among these, norcantharidin was identified as a specific ligand for PPP1CA, while naringenin exhibited exclusive targeting toward CTSD. Notably, sarin, benzo[a]pyrene, and hydrogen peroxide demonstrated dual targeting capabilities for both biomarkers (Figure 6A). Molecular docking simulations were conducted to validate binding interactions between the biomarkers and their top-scoring ligands. PPP1CA showed the strongest affinity for norcantharidin, with a binding energy of −6.2 kcal/mol (Table 2; Figure 6B), with key interactions with pocket 2 (primarily involves residues on protein chain B, such as PRO24, TYR69, and ASP71). Similarly, CTSD exhibited optimal structural compatibility with naringenin, achieving a binding energy of −7.4 kcal/mol (Table 2; Figure 6C), with key interactions with pocket 1 (primarily involves residues on protein chains C and D, including TYR78, HIS77, ASP33, ILE124, PHE126). The favorable binding energies observed in molecular docking supported the therapeutic potential of these compounds as modulators of PPP1CA and CTSD activity in T2DM pathogenesis.