We obtained our data from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). Using “type 2 diabetes mellitus” and “osteoarthritis” as search keywords, we selected the GSE26168 and GSE114007 datasets for our screening analysis. The selected dataset which contains more samples and comprehensive clinical information has good quality and meets the requirements of research analysis. The GSE26168 dataset consists of gene expression data from the blood tissue of nine patients with type 2 diabetes mellitus (T2D) and eight healthy controls. The GSE114007 dataset includes gene expression information from the cartilage tissue of 20 patients with osteoarthritis (OA) and 18 normal controls. To validate our findings, we also retrieved the GSE20966 (T2D) and GSE51588 (OA) datasets from the GEO database.

We utilized RStudio software (version 4.2.1) and the “Limma” R package to process the datasets and identify DEGs between the disease and control groups (18, 19). To ensure accuracy, we excluded probe sets without corresponding gene symbols and calculated averages for genes represented by multiple probe sets. Our significance threshold was set at p<0.05 and |logFC|>1. Additionally, we used the Venn diagram online analysis tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) to identify common DEGs between T2D and OA datasets.

To analyze the functional characteristics of genes shared by T2D and OA, we utilized RStudio software for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, using the “cluster Profiler” R package. GO is an internationally recognized gene functional classification system that provides a comprehensive and regularly updated vocabulary, as well as well-defined concepts for identifying genes and their products (20). It consists of three ontologies: molecular function, cellular component, and biological process. GO enrichment analysis allows for the identification of significantly enriched GO terms compared to the genome background, and associates DEGs with specific biological functions. KEGG, a principal public database for pathway-related information, enables the identification of pathways that are significantly enriched in DEGs compared to the genome-wide background (21).

Using the STRING online tool (http://string-db.org), we explored the interactions among common DEGs and constructed a Protein-Protein Interaction (PPI) network with extensive regulatory links (22). The STRING database facilitates the search for relationships between proteins, encompassing direct binding relationships and interconnected upstream and downstream regulatory pathways. Interactions with a confidence value above 0.4 were deemed statistically significant, with all other parameters set to default. We employed Cytoscape software (version 3.8.2) for PPI network analysis and visualization and utilized the Cytoscape plug-in, Molecular Complex Detection (MCODE), for identifying primary functional modules (setting k-core=2, degree cutoff=2, max depth=100, node score cutoff=0.2) (23, 24). Key genes were identified using Cytoscape’s CytoHubba plugin (25), and an upset diagram was used to filter and determine shared key genes. Additionally, we established a co-expression network for these key genes using GeneMANIA (http://www.genemania.org/), a tool renowned for uncovering the internal relationships within gene sets (26).

We utilized the GSE20966 (T2D) and GSE51588 (OA) datasets to validate the mRNA expression of key genes. This process aimed to definitively identify the key genes for T2D and OA. The GSE20966 dataset includes samples from 10 T2D patients and ten controls, while the GSE51588 dataset comprises 40 OA patient samples and ten control samples.

Quantitative real-time polymerase chain reaction (qRT-PCR) was employed to assess the mRNA expression of key genes in T2D (n=8, peripheral blood), OA (n=7, articular cartilage), T2D combined with OA (T2D+OA) samples (n=4,articular cartilage), and normal controls (n=9,peripheral blood). The inclusion and exclusion criteria for the clinical study were listed in Supplementary Table 1 . Clinical tissue samples were acquired from Panzhihua Central Hospital. For the study, we used TRIpure Total RNA Extraction Reagent (EP013, ELK Biotechnology) to extract RNA from the samples. The isolated total RNA was reverse-transcribed into complementary DNA (cDNA) using the EntiLink™ 1st Strand cDNA Synthesis Super Mix (EQ031, ELK Biotechnology). qRT-PCR was performed with the QuantStudio 6 Flex System PCR instrument (Life Technologies) and the EnTurbo™ SYBR Green PCR SuperMix kit. Each sample was prepared in triplicate wells. β-actin served as the internal reference gene, and the 2-ΔΔCt method was used for calculating mRNA fold changes and normalizing the data across all samples. The primer sequences for PCR are listed in Supplementary Table 2 .

For protein extraction, we initially rinsed 3 samples each of T2D peripheral blood, OA articular cartilage, T2D+OA articular cartilage, and peripheral blood from healthy controls with pre-chilled PBS buffer (ASPEN, AS1025). This process was repeated 2-3 times. Subsequently, RIPA total protein lysate (ASPEN, AS1004) was added in volume 10-20 times that of the tissue samples. The mixture was then centrifuged at 12,000 rpm at 4°C for 5 minutes, and the supernatant was collected. Protein concentration in the samples was determined using the BCA Protein Concentration Assay Kit (ASPEN, AS1086). For western blot (WB), SDS-PAGE electrophoresis was performed. An appropriate amount of 5×protein loading buffer (ASPEN, AS1011) was added to the samples, which were then incubated in a boiling water bath at 95-100°C for 5 minutes. Electrophoresis was conducted at 80V for the stacking gel and 120V for the separation gel. A PVDF membrane (Millipore, IPVH00010) was activated with methanol for 3 minutes before use. The transfer of proteins onto the membrane was carried out at a constant flow of 300mA.The membranes were immunoblotted with the diluted primary antibodies (ASPEN, AS1061) for overnight at 4°C. The blots were washed thrice with TBST. Afterward, the diluted secondary antibody was incubated at room temperature for 30 minutes. Excess antibody was washed off with TBST for 4 times. Immunoreactivity was detected using ECL western blot reagent. Finally, the signal bands were quantified by densitometry analysis using the AlphaEaseFC software processing system after scanning the blotted membrane.

To ascertain the diagnostic value of pivotal genes in OA, we utilized external datasets GSE20966 (T2D) and GSE51588 (OA) to validate clinical diagnosis prediction models. Employing RStudio software (version 4.2.1), we used the pROC package to generate Receiver Operating Characteristic (ROC) curves. Additionally, logistic regression analysis was conducted to evaluate the role of key genes in distinguishing between T2D, OA, and healthy individuals. This analysis facilitated the establishment of a clinical diagnosis prediction model based on the identified key genes.

Transcription factors (TFs) are proteins that bind to specific DNA sequences, forming complex regulatory systems to control gene expression. We utilized Network Analyst 3.0 (https://www.networkanalyst.ca/) to examine the interactions between key genes and transcription factors in T2D and OA (27). This analysis aimed to evaluate the impact of TFs on the expression and functional pathways of key genes. The transcription factor and gene target data were sourced from the ENCODE ChIP-seq database (28). We applied a peak intensity signal threshold of <500 and a regulatory potential score of <1 point, as predicted using the BETA negative algorithm. Furthermore, we employed Cytoscape to visually represent the TFs-mRNA regulatory network.

MicroRNAs (miRNAs) are endogenous, short, non-coding RNAs that play a crucial role in inhibiting or degradating target mRNAs. In order to better understand how gene expression is affected in different physiological and disease contexts, we analyzed the regulatory network between miRNAs and mRNAs (29). To do this, we utilized miRTarBase (version 8.0), a database specifically designed for predicting miRNA binding sites, to identify potential interactions between key genes and miRNAs. This involved conducting a cross-analysis of mRNA-miRNA binding and constructing a competing endogenous RNA (ceRNA) network. The results of this analysis were then visualized using Cytoscape software (version 3.8.0).

The Drug-Gene Interaction Database (DGIdb) (http://www.dgidb.org) is a comprehensive resource that provides information on interactions between drugs and genes (30). In this study, we utilized DGIdb to identify potential drugs that target pivotal genes, supplemented by candidate compounds sourced from existing literature for treating T2D and OA. We downloaded protein and ligand files for MMP9 and ANGPTL4 from the RCSB Protein Data Bank (https://www.rcsb.org/) in “PDB”. The PDB format focuses on storing 3D structural data of biological macromolecules such as proteins and nucleic acids, and contains detailed information of protein-ligand binding sites. It is more suitable for the screening of predictive gene drugs. The SDF format is mainly used to store and exchange detailed molecular structure information, including 2D or 3D coordinates of molecules, atomic types, bond information. It is suitable for occasions that require batch processing of large amounts of molecular data, such as virtual screening and compound library management (31). Therefore, after establishing the corresponding coordinates between the proteins and ligands, we converted the protein files to the “PDBQT” format. The “SDF” format ligands of candidate compounds were also converted into “PDBQT” format using the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (32). Molecular docking was performed using AutoDock (Vina 1.2.2), which calculated the binding affinity between receptors and ligands. A binding energy lower than −5.0 kcal/mol typically indicates strong binding activity to the target protein. We analyzed the docking results with the Protein-Ligand Interaction Profiler (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index) and visualized them using PyMOL (version 2.0). This allows us to demonstrate the ligand-receptor binding process through hydrogen bonds and amino acid residues.

All data were processed and analyzed using the GraphPad Prism software 8.0 (GraphPad Software, La Jolla, CA, USA), and results were expressed as mean ± standard deviation (SD). Significant differences between the groups were evaluated by unpaired student’s t-test. Unidirectional Analysis of Variance (ANOVA) with Tukey’s multiple comparison tests was utilized where appropriate in cases involving multiple groups. A p-value of less than 0.05 was regarded as significant.

We analyzed tissue samples from GEO, performing data on the datasets for T2D and OA. In the T2D dataset (GSE26168), we identified 1,746 DEGs. Similarly, in the OA dataset (GSE114007), we identified 2,910 DEGs. These DEGs for both T2D and OA were visually represented through volcano plots and heat maps ( Figures 2A–D ). A Venn diagram was utilized to highlight the overlap between the two datasets, revealing 209 common DEGs ( Figures 2E, F ). The results underscore the presence of many common genes shared between T2D and OA.

GO and KEGG enrichment analyses were conducted on the 209 common DEGs. The GO analysis ( Figure 2G ) revealed significant gene enrichment in various categories. In Biological Processes (BP), the most enriched pathways included ossification, ameboidal-type cell migration, wound healing, and regulation of epithelial cell proliferation. Within the Cellular Components (CC) category, collagen-containing extracellular matrix, presynapse, and neuron projection terminus were most prominent. Regarding Molecular Functions (MF), the top enriched functions were extracellular matrix structural constituent (ECM), glycosaminoglycan binding, and sulfur compound binding. KEGG analysis ( Figure 2H ) indicated substantial enrichment in pathways such as ECM-receptor interaction, PI3K-Akt signaling pathway, Wnt signaling pathway, intestinal immune network for IgA production, and protein digestion and absorption. These results strongly suggest that the common DEGs between T2D and OA play crucial roles in regulating chemokines, cytokines, and inflammatory responses. These genes are implicated in cellular processes like production, metabolism, and apoptosis and are significantly associated with the onset and progression of T2D and OA.

We conducted a comprehensive analysis of the DEGs shared by T2D and OA using the STRING online network tool and Cytoscape software to elucidate PPI. A universal PPI network was constructed with a minimum interaction score of 0.4, consisting of 207 nodes and 226 edges ( Figure 3A ). This complex PPI network showed that 96.86% of the genes had co-expression and 2.05% shared protein domains, and 1.08% had interaction responses. In this network, nodes represent proteins, and edges represent their interactions. The size and color intensity of a node indicates its degree value, while the thickness of an edge reflects the strength of the relationship between the proteins. Using the Cytoscape plug-in MCODE, we identified the most significant gene module ( Figure 3B ). Furthermore, we used the Cytoscape’s cytoHubba plug-in to determine the top 14 hub genes, which included MMP9, NGFR, CXCL12, CACNA1A, SERPINE1, CCR9, CALCA, NRCAM, MME, ANGPTL4, PTPRU, DLL4, SELP, and FAIM2 ( Figures 3C–E ). To further explore the functions of these hub genes, we analyzed their co-expression network and related functions using the GeneMANIA database ( Figure 3F ).

To validate the mRNA expression levels of the 14 identified hub genes, we utilized an additional T2D dataset and an OA dataset. The analysis of these datasets revealed significant expression differences in the hub genes. In the T2D external dataset (GSE20966), three hub genes-ANGPTL4, MMP9, and NRCAM-showed substantial differential expression compared to the normal control ( Figure 4A ). Similarly, in the OA external dataset (GSE51588), three hub gene-ANGPTL4, MMP9, and CACNA1A-exhibited significant expression differences ( Figure 4B ). By intersecting the results from both datasets, ANGPTL4 and MMP9 were conclusively identified as the pivotal genes shared between T2D and OA.

Firstly, we statistically analyzed the clinical data of the tested samples and established a baseline table of clinical information ( Table 1 ). In the limited clinical cases, we can find that patients in the T2D + OA group had the highest blood glucose and Glycated hemoglobin A1c (GHbA1c), and 5 patients who had poor long-term blood glucose control in T2D group were rated grade 1 or 2 for Kellgren-Lawrence (KL). Four patients in OA group were in the abnormal range of blood glucose. There were no significant differences in age, sex and bone mineral density among the groups. It suggested that hyperglycemia may promote the development of osteoarthritis. Secondly, As illustrated in the figure, the test results indicated significant differences in the mRNA expression levels of these genes among the different groups. Specifically, MMP9 expression was significantly higher in T2D patients (p=0.046), OA patients (P=0.009), and T2D + OA patients (p=0.048) compared to healthy individuals ( Figure 5A ). Similarly, ANGPTL4 showed significantly elevated mRNA expression levels in T2D patients (P=0.001), OA patients (p=0.013), and T2D + OA patients (p=0.045) relative to healthy controls ( Figure 5B ).

Combined with the clinical data above, we can find that blood glucose was also higher in T2D and T2D + OA patients with higher MMP9 and ANGPTL4 expression. In T2D patients who had not been diagnosed with OA, higher expressions of MMP9 and ANGPTL4 were associated with higher KL grades. It suggests that the expression of MMP9 and ANGPTL4 may be positively correlated with blood glucose and the severity of osteoarthritis.The details on clinical data above can be found in the Supplementary Table 3 .

We utilized western blot analysis to evaluate the protein levels of ANGPTL4 and MMP9 in tissue samples from T2D patients, OA patients, T2D+OA patients, and healthy controls. The results showed a significant increase in MMP9 protein in T2D patients (5.57 times higher, p < 0.001), OA patients (9.71 times higher, p < 0.001), and T2D+OA patients (3times higher, p = 0.046), when compared to healthy controls ( Figures 5C, E ). Similarly, ANGPTL4 protein level was elevated in in T2D patients (6 times higher, p=0.001), OA patients (3.6 times higher, p=0.036), and T2D+OA patients (10 times higher, p < 0.001) compared to the healthy group ( Figures 5D, E ). The original western blot gel was found in the Supplementary Figure 1 .

The method focused on assessing the area under curve (AUC) values to determine the sensitivity and specificity of these genes in diagnosing T2D and OA. In the T2D dataset, both ANGPTL4 and MMP9 had AUC values greater than 0.8 ( Figures 6A, B ), indicating a high level of diagnostic accuracy. Similarly, in the OA dataset, both genes had AUC values exceeding 0.75 ( Figures 6D, E ), confirming their strong diagnostic value for OA. Furthermore, we integrated ANGPTL4 and MMP9 with multiple markers to enhance prognosis prediction and developed a multi-marker diagnostic model using logistic regression analysis. The ROC results demonstrated that this multi-marker model effectively predicted the diagnosis of T2D (AUC = 0.980) and OA (AUC = 0.855) ( Figures 6C, F ), highlighting its potential clinical utility.

The results in the construction of an interaction network between key genes and TFs ( Figure 6G ). The network consists of 50 nodes and 57 edges. Notably, ANGPTL4 is influenced by 35 TF genes, while MMP is regulated by 22 genes. Our analysis revealed that different key genes are modulated by distinct TFs, with 7 TFs (ZBTB26, MXD3, DRAP1, GATAD1, SSRP1, TFAP4, MXD4) simultaneously regulating both key genes. Additionally, we used miRTarBase to predicted miRNA binding sites and establish an interaction regulatory network between key genes and miRNA ( Figure 6H ). This network consists of 24 nodes and 25 edges. In this configuration, ANGPTL4 is regulated by 4 TF genes and MMP9 by 21 genes. Our results indicate that different miRNAs regulate distinct key genes, with hsa-mir-29b-3p concurrently regulating both key genes.

Our research has identified 10 candidate compounds that interact with ANGPTL4 and MMP9, sourced from DGIdb ( Supplementary Table 4 ). Additionally, we have identified 10 potential therapeutic drugs for T2D combined with OA from relevant literature references (33–37). Using AutoDock, we have calculated the binding energies between the key target proteins MMP9 (PDB ID: 6ESM) and ANGPTL4 (PDB ID: 6U1U) and the candidate compounds. Upon comparing the binding energies to those of the original ligands, we have observed that in the MMP9 molecular docking group, two compounds exhibited binding energies equal to or greater than the original ligand of the protein. The drug with the highest binding affinity in this group was identified as Raloxifene ( Table 2 ). In the ANGPTL4 molecular docking group, 17 compounds had binding energies equal to or greater than that of the original protein-ligand. The candidate drug showing the greatest binding activity in this group was Ezetimibe ( Table 3 ). Notably, Raloxifene and S-3304 exhibited binding energies surpassing the original protein ligands of ANGPTL4 and MMP9. After analyzing the docking results with the best binding properties for MMP9 and ANGPTL4 through the Protein-Ligand Interaction Profiler website, we used PyMOL to visualize the molecular docking interactions ( Figures 6I–L ).