Introduction
Osteoarthritis (OA) is a degenerative joint disease with joint pain and dysfunction as the main symptoms, and is one of the most common causes of chronic pain and disability in the elderly ( Yuan et al., 2020 ; Molnar et al., 2021 ) which has attracted wide attention due to its high prevalence and high disability rate. According to statistics, about 300 million people around the world suffer from OA ( Abramoff and Caldera, 2020 ), which is particularly common in middle-aged and elderly people. The prevalence rate can reach 50% in people over 60 years old, and can reach 80% in people over 75 years old, and the disability rate can reach 53% ( Quicke et al., 2022 ), which not only impairs patients’ quality of life, but also brings heavy burden to families and society. It is understood that genetics, age, obesity, trauma, inflammation, and metabolic changes are all considered as risk factors for OA ( Chaganti and Lane, 2011 ). However, so far, the etiology and pathogenesis of OA have not been fully elucidated, and current treatment strategies, including oral non-steroidal anti-inflammatory drugs (NSAIDs), intra-articular injection, and joint arthroplasty, mainly focus on pain alleviation. However, pharmacological interventions are associated with gastrointestinal, cardiovascular, and other adverse reactions, whereas surgical treatment carries a risk of postoperative infection ( Bourget-Murray et al., 2023 ). Therefore, identifying diagnostic and therapeutic targets for OA continues to be a key focus of ongoing research. Our study aims to advance OA treatment by elucidating novel biomarkers and effective drug targets.
Amino acids (AAs) are small molecule metabolites, which not only participate in the construction of peptides and proteins, but also play an important role in immune function, cytokine secretion, and response. Some studies using ultra-high performance liquid chromatography-quadrupole time-of-flight mass spectrometry found that there were differences in AAs spectra between healthy controls and OA patients, and arginine and proline metabolism were the most significantly affected metabolic pathways ( Chen et al., 2018 ). A recent multi-omics analysis of synovial tissue and synovial fluid has demonstrated localized dysregulation of arginine metabolism in the osteoarthritic joint, as evidenced by significantly elevated arginine levels in the synovial membrane and enrichment of the arginine biosynthesis pathway in the synovial fluid, a pattern of local accumulation and enhanced biosynthetic activity that promotes cartilage degradation and synovial inflammation ( Ge et al., 2025 ; Wang et al., 2026 ). Proline, the other key amino acid in this pathway, is essential for collagen synthesis and the maintenance of extracellular matrix stability ( Karna et al., 2020 ). A recent study has demonstrated that the circular RNA circUbqln1 may impair proline metabolism in OA by modulating PRODH activity, thereby inhibiting collagen synthesis and compromising the structural integrity of articular cartilage ( Feng et al., 2024 ). These findings indicate that arginine accumulation and pathway enrichment in synovial tissues promote cartilage degradation and synovial inflammation, while proline dysregulation inhibits collagen synthesis and compromises cartilage structural integrity. It has been demonstrated that arginine and proline metabolism play an important role in the occurrence and development of OA. Therefore, the identification of key biomarkers and potential therapeutic targets related to arginine and proline metabolism can provide positive value for the diagnosis and treatment of OA.
Fortunately, in recent years, the development of transcriptome sequencing technology and the application of bioinformatics in the medical field have provided new approaches for studying the molecular mechanisms of various diseases. We collected synovial tissue data from healthy individuals and OA patients. Through the differential expression analysis of transcriptome sequencing and weighted gene co-expression network analysis (WGCNA), candidate genes associated with arginine and proline metabolism were identified and functionally characterized. Multiple bioinformatics methods were further employed to screen and pinpoint hub genes, enabling efficient determination of the functions and roles of differentially expressed genes. These findings were subsequently validated experimentally. Additionally, the pharmacological activities of two potential OA drug targets were verified via drug prediction, molecular docking, and molecular dynamics simulations, offering novel directions and potential therapeutic targets for the prevention and treatment of OA ( Figure 1 ).
Materials and methods
Synovial tissue sample collection and sequencing
Human synovial tissue samples were collected from OA patients and healthy controls from the Fuzhou Second General Hospital. The training set included 12 OA patient samples and 10 normal control samples, all from human donors. This research was approved by the Hospital Ethics Committee (Ethics Number: 2021014). All patients have signed the informed consent form. All operations were carried out strictly in accordance with the relevant guidelines and regulations of the hospital.
As per the manufacturer’s protocol, TRIzol reagent (Invitrogen, CA, USA) was utilized to separate and purify the total RNA from synovial tissue samples. RNA quantity and purity were evaluated with a NanoDrop ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, USA), and its integrity was tested via a Bioanalyzer 2,100 system (Agilent, CA, USA). Validation was further confirmed by agarose gel electrophoresis. For downstream analyses, samples were required to have a concentration over 50 ng/μL, RIN values above 7.0, OD260/280 ratios higher than 1.8, and total RNA amounts greater than 1 μg. After two rounds of purification with oligo (dT) magnetic beads (Dynabeads Oligo (dT), Thermo Fisher, USA), polyadenylated mRNA was isolated. The mRNA fragmentation was achieved at 94 °C for 5–7 min via the NEBNext® Magnesium RNA Fragmentation Module. The fragmented RNA was reverse transcribed into cDNA with the help of SuperScript™ II Reverse Transcriptase (Invitrogen). Double-stranded cDNA was synthesized via E. coli DNA polymerase I and RNase H (NEB), incorporating dUTP Solution (Thermo Fisher) to generate blunt-ended fragments. A single adenine base was added to the 3′ends to facilitate ligation with adapters containing a complementary thymine base. The library was purified and size-selected with the help of magnetic beads. Subsequently, double-stranded DNA was treated via UDG enzyme (NEB) to remove uracil residues. The DNA library underwent PCR amplification under the following conditions: 3-min denaturation at 95 °C, 8 cycles of 98 °C for 15 s, 60 °C for 15 s, and 72 °C for 30 s, finishing with a 5-min extension at 72 °C. A library with an average fragment size of 300 ± 50 bp was generated. Paired-end sequencing was subsequently performed on the Illumina Novaseq™ 6,000 platform (LC Bio Technology CO., Ltd., Hangzhou., China) with a 150 bp read length mode according to the manufacturer’s standard protocol. Raw sequencing reads were generated and archived in FASTQ format. Quality of the data was assessed with FastQC (v. 0.11.9) to remove low-quality reads ( Brown et al., 2017 ). Phred quality scores were used to evaluate base call accuracy, with a score ≥30 defined as acceptable, and nucleotide distribution was analyzed to detect any AT/CG bias. Clean reads were subsequently mapped to the human reference genome, which was downloaded via the Ensembl database ( https://asia.ensembl.org/info/about/species.html ), using Hierarchical Indexing for Spliced Alignment of Transcripts 2 (Hisat2). Transcripts assembly were carried out with Cufflinks (v. 2.2.1) to produce count data ( Trapnell et al., 2012 ). To assess sample clustering and distinguish between OA and control groups, principal component analysis (PCA) was executed with the help of “scatterplot3d” (v. 0.3–42, Li et al., 2022 ).
Gene expression dataset screening and processing
For the testing set, the GSE 1919 ( GPL91 ) dataset via the GEO database ( https://www.ncbi.nlm.nih.gov/geo/ ) was used, which encompassed 5 synovial tissue normal control samples and 5 OA patient samples. Additionally, 54 APRGs were extracted via the Molecular Signatures Database (MSigDB) ( https://www.gsea-msigdb.org/gsea/msigdb/ ), specifically from the KEGG arginine and proline metabolism gene set ( Wang et al., 2023 ). Gene expression data were subjected to Median Ratio Normalization using the DESeq2 (v. 3.50.3, Pang et al., 1989 ) package to ensure the comparability of gene expression levels among samples. Once the data were quality-controlled, raw counts were normalized for library size using the DESeq2 estimate size factors function (median-of-ratios method), and the normalized counts were used for downstream analysis.
Differential expression analysis
The DESeq2 package was employed to identify DEGs by comparing OA and control samples within the training set [|log2fold change (FC)| > 0.5 and p < 0.05]. Based on log2FC, most 10 upregulated and downregulated genes were established and ranked in descending order, subsequently highlighted in the volcano plot via the “ggplot2” R package (v. 3.3.6, Gustavsson et al., 2022 ), and a heatmap was built by the “ComplexHeatmap” R package (v. 2.14.0, Gu and Hübschmann, 2022 ).
Weighted gene Co-expression network analysis (WGCNA)
After identifying DEGs, the intersection of APRGs and DEGs was taken to pinpoint the differentially expressed APRGs in OA versus control groups. Using all samples within the training set, single-sample Gene Set Enrichment Analysis (ssGSEA) scores for the differentially expressed APRGs were processed via “GSVA” (v. 1.42.0, Hänzelmann et al., 2013 ). A Wilcoxon test was adopted to compare differences in APRG scores across the OA and control cohorts (p < 0.05). Subsequently, WGCNA was performed on the APRGs scores using the “WGCNA” R package (v. 1.71, Langfelder and Horvath, 2008 ) for building a co-expression network and detecting modules within the training set. First, hierarchical clustering was performed on the training set using the hclust function to identify potential deviant samples. Outliers were omitted to safeguard the accuracy of upcoming analyses. After outlier removal, the ssGSEA scores of APRGs for each sample were introduced, and hierarchical clustering was re-executed via the hclust function. To construct a co-expression network, the pickSoftThreshold function was harnessed to determine optimal soft-threshold power, ensuring that the network maximally approximated a scale-free topology. A soft-threshold parameter was selected from a range of 1 to 30, with the goal of identifying the appropriate threshold for gene-gene correlations. The R 2 threshold was set to 0.85, ensuring a scale-free network with lower mean connectivity. To better represent gene interaction strength and similarity, a soft-threshold power was chosen to maintain a scale-free topology, refining the accuracy of the adjacency matrix and topological overlap matrix (TOM). Based on the TOM, a gene dendrogram was generated using hierarchical clustering to group genes into modules. The dynamic tree-cutting algorithm identified modules, specifying a minimum size of 100 genes. Finally, modules were merged based on a height cutoff of 0.4, resulting in candidate modules for further analysis.
Identification of candidate genes
Further correlation analysis was performed between the ssGSEA scores of APRGs and the gene modules. Key modules with the strongest positive and negative correlations (|cor| > 0.3, p < 0.05) were opted for deeper exploration. Subsequently, module membership (MM) and gene significance (GS) were calculated for key module genes. MM represents the correlation between genes and their respective modules, while GS reflects the correlation between genes and the trait of interest. Scatter plots were rendered via “ggplot2”. Key module genes were those with |MM| above 0.9 and |GS| above 0.6. A Venn diagram was generated via “ggVenn” (v. 1.2.2, Zheng et al., 2022 ) to identify the intersection between DEGs and APRGs-related key module genes, overlapping genes were pinpointed as candidate genes.
Functional analysis of candidate genes
Then, candidate genes underwent gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis (p < 0.05) via the “clusterProfiler” R package (v. 4.7.1.003, Wu et al., 2021 ), with the human gene set from the “org.Hs.e.g.,.db” as the background. GO analysis was performed across 3 categories: biological process (BP), cellular component (CC), and molecular function (MF), and 5 most significant genes from each section enriched pathways were rendered via “GOplot” (v. 1.0.2, Walter et al., 2015 ). Similarly, the KEGG pathway analysis results were also visualized with the “GOplot” R package. To explore PPI, the candidate genes were uploaded to the search tool for the retrieval of interacting the genes (STRING) database ( https://cn.string-db.org/ ), the PPI network was established with a confidence score threshold of 0.40. Finally, PPI networks were visualized via Cytoscape software (v. 3.9.1, Shannon et al., 2003 ). The top 40 candidate genes were selected from five algorithms individually: Degree, Density of Maximum Neighborhood Component (DMNC), Closeness, Clustering Coefficient, as well as Edge Percolated Component (EPC). The intersecting genes from these selections were used for subsequent analysis.
Screening for key biomarkers
The intersecting genes were subjected to machine learning using random forest analysis by “randomForest” (v. 4.7–1.1, Alderden et al., 2018 ) to identify feature genes. The top 5 feature genes, ranked by their importance based on the mean decrease in the Gini index, were selected for further experiments. To assess their ability to differentiate OA samples from control samples, receiver operating characteristic (ROC) curves were implemented through “pROC” (v. 1.18.5, Robin et al., 2011 ), and area under the curve (AUC) was quantified on both databases. Feature genes with AUC values greater than 0.7 were evaluated for expression differences between OA and control in both databases using the Wilcoxon test (p < 0.05). Box plots were generated using “ggplot2” to visualize differential expression (p < 0.05). Genes with marked differential expression between groups and consistent trends across both datasets were identified as hub genes.
Creation and evaluation of nomograms based on hub genes
In the training set, a nomogram model was built from the expression levels of hub genes, with OA samples considered as the outcome event. The “rms” R package (v. 6.3–0, Sachs, 2017 ) was used to construct this nomogram. To evaluate the model’s predictive ability, calibration curves were generated via “rms”. A p-value >0.05 and a slope close to 1 indicated higher predictive accuracy. Additionally, ROC curves rendered through “pROC”, with an AUC >0.7 considered indicative of good model accuracy.
Gene set enrichment analysis (GSEA)
GSEA was executed with the help of gene sets from the MSigDB database, with the KEGG background gene set from “org.Hs.e.g.,.db”. For each of the hub genes in the training set, correlations were calculated across hub genes and all remaining genes to generate a prioritized catalog of associated genes. Then, GSEA was conducted based on this ranking, with the significance threshold set at p. adjust <0.05 as well as |Normalized Enrichment Score (NES)| > 1.
Immune infiltration assessment
In the training set, enrichment scores of 64 immune cell types were employed to highlight variation between OA and control samples via the xCell algorithm implemented through the “xCell” R package (v. 1.1.0, Aran et al., 2017 ). Next, the Wilcoxon test (p < 0.05) was used to identify significantly different immune cell types between the OA and control groups. To gain more insight into the interaction between hub genes and significantly different immune cells, Spearman correlation analysis was done on all samples in the training set via “psych” (v. 2.2.9, Robles-Jimenez et al., 2021 ). |cor| > 0.3 and p < 0.05 were regarded as evidence of a correlation, and a moderate correlation (|cor| > 0.38) was considered statistically significant.
Drug prediction, molecular docking, and molecular dynamics simulation
To evaluate the clinical potential of hub genes, the Drug Signatures Database (DSigDB) ( http://tanlab.ucdenver.edu/DSigDB , accessed on 10 November 2025) within EnrichR (v. 3.2, Sheng et al., 2025 ) was utilized in this study to predict candidate compounds interacting with hub genes. The compound with the highest combined score for each key gene was selected as the ligand for molecular docking analysis. The 3D structures of hub genes and their corresponding compounds were retrieved from the RCSB PDB database ( https://www.rcsb.org/ ) and PubChem database ( https://pubchem.ncbi.nlm.nih.gov/ ), respectively. Molecular docking was performed using the CB-Dock online server ( https://cadd.labshare.cn/cb-dock/php/index.php ), a tool implementing blind docking based on AutoDock Vina that can automatically identify binding pockets and calculate binding free energy. To verify the stability of the docked conformations, human-derived protein sequences were further obtained from NCBI, and 100 ns molecular dynamics (MD) simulations were conducted using GROMACS 2024.2 software (the simulated systems included MYOM2-MELPHALAN and TCAP-Testosterone Enanthate). In the MD simulations, the AMBER14SB force field and GAFF force field were employed to generate parameters and topology files for proteins and small-molecule ligands, respectively. The system was placed under periodic boundary conditions, and the simulation box was solvated with TIP3P water molecules. To maintain system electroneutrality, 0.15 mol/L Na + and Cl − ions were added. Energy minimization was accomplished using the steepest descent method. Subsequently, two-stage pre-equilibration was performed: first, 100 ps temperature equilibration at 300 K under the NVT ensemble, followed by 100 ps pressure equilibration at 1 bar under the NPT ensemble. MD simulations were performed to evaluate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), total energy, and number of hydrogen bonds of the complexes, respectively.
RT-qPCR analysis of synovial tissue in human OA
The Real-time quantitative polymerase chain reaction (RT-qPCR) experiment set included OA patient samples and normal control samples, all self-collected from human donors. Subsequently, the synovial tissues collected from clinical samples were minced with bone forceps and then pulverized into fine powder using a mortar pre-cooled with liquid nitrogen. Total RNA was extracted from each tissue sample following the protocol of the TRIzol kit, and the RNA concentration was quantified. Using the reverse transcriptase HiScript Q RT SuperMix for qPCR (+gDNA wiper), the reverse transcription system was prepared as per the manufacturer’s instructions, residual genomic DNA was removed, and complementary DNA (cDNA) was synthesized via reverse transcription. Subsequently, 2 mL cDNA was incorporated into the ChamQ SYBR qPCR Master Mix reaction system for RT-qPCR. The thermal cycling conditions comprised denaturation at 95 °C for 30 s, annealing at 95 °C for 10 s, and extension at 60 °C for 30 s for a total of 40 cycles. Following the acquisition of CT values for MYOM2 and TCAP genes in each group, the 2 −ΔΔCT method was used to determine the relative expression levels of each gene. Primer sequences are detailed in Table 1 .
TABLE 1
The gene primer sequences of MYOM2 and TCAP.
| Gene | Primer sequence (5'→3′) | Length (bp) | |
| MYOM2 | Forward primer Reverse primer | ATGGCGTACACACACTGGAG ACCGAACGGAATGGTTCCTC | 151 |
| TCAP | Forward primer Reverse primer | CATGAGACCTACCACCAGCA ATCTTGGCAGGGGTGAAGAT | 155 |
| GAPDH | Forward primer Reverse primer | GGTGTGAACCATGAGAAGTATGA GAGTCCTTCCACGATACCAAAG | 123 |
Results
APRG-associated key modules reveal dysregulated genes in OA
After quality control, no AT/CG separation was observed ( Supplementary Figure S1A ). The transcriptome data from self-collected synovial tissue samples showed a mapping rate above 90% for all samples, confirming high data quality ( Table 2 ). PCA revealed clear separation between OA and control samples, indicating distinct transcriptomic profiles in two groups ( Supplementary Figure S1B ). Differential expression analysis comparing OA and control samples revealed 2,227 DEGs. Specifically, 1,435 genes with increased expression and 792 with decreased expression in OA were identified ( Figures 2A,B ). The 54 APRGs were obtained from the MSigDB. By intersecting 54 APRGs with the 2,227 DEGs, 8 differentially expressed APRGs were identified comparing OA and control groups ( Figure 2C ). The ssGSEA scores of APRGs were markedly elevated in the OA group relative to the control group (p < 0.05), suggesting a strong correlation between APRG expression and OA ( Figure 2D ). Outlier samples in this study were excluded. Following sample exclusion, ssGSEA scores of APRG were recalculated for each sample, and clustering analysis was performed ( Supplementary Figure S1C,D ). A soft threshold of R2 = 0.85 was established, and a power value of 24 was chosen to achieve a scale-free network, with a low mean connectivity. Thus, the optimal soft threshold (β) was set to 24 ( Supplementary Figure S1E ). Correlation analysis identified 17 modules associated with the ssGSEA scores of APRGs ( Supplementary Figure S1F ). Among these, the MEblack and MEbrown modules were identified as key modules, with MEblack showing the top positive correlation (cor = 0.83, p < 0.001) and MEbrown showing the largest negative correlation (cor = −0.57, p < 0.05) with APRG scores ( Figure 2E ). A total of 120 key module genes for further analysis, with 105 black module and 15 brown module genes, with MEblack and MEbrown modules showing strong module membership–trait correlations (MEblack: cor = 0.87; MEbrown: cor = 0.37) ( Figure 2F ).
TABLE 2
Comparison rate of each sample.
| Sample | Total_reads | Mapped (%) |
| OA2 | 35,827,848 | 90.31% |
| OA3 | 36,794,650 | 93.55% |
| OA4 | 36,761,617 | 92.76% |
| OA5 | 33,206,496 | 94.10% |
| OA6 | 37,609,448 | 93.93% |
| OA7 | 37,404,510 | 92.16% |
| OA8 | 36,850,128 | 93.79% |
| OA9 | 38,197,749 | 94.08% |
| OA10 | 38,467,259 | 93.31% |
| OA11 | 35,937,128 | 90.11% |
| OA12 | 34,028,720 | 91.72% |
| OA13 | 40,332,314 | 93.48% |
| control1 | 36,567,931 | 93.80% |
| control2 | 32,318,946 | 93.52% |
| control3 | 37,786,353 | 94.26% |
| control4 | 36,040,191 | 94.58% |
| control5 | 31,607,790 | 93.77% |
| control6 | 31,087,738 | 93.23% |
| control7 | 38,281,588 | 94.84% |
| control8 | 38,171,818 | 94.14% |
| control10 | 34,778,731 | 93.58% |
| control11 | 48,291,671 | 93.17% |
Candidate genes converge on muscle-related pathways
By intersecting the 2,227 DEGs, 120 key module genes, 96 genes were identified and designated as candidate genes ( Figure 3A ). GO enrichment analysis identified 302 terms, including 209 BP, 41 CC, 52 MF. BP terms were mainly enriched in muscle-related processes (muscle organ development, muscle system process, and muscle tissue development); CC terms in contractile structures (contractile fiber, myofibril, and sarcomere); and MF terms in actin binding and structural constituents of muscle ( Figure 3B ). The KEGG pathway enrichment analysis identified 15 significantly enriched pathways, including cytoskeleton in muscle cells, motor proteins, cardiac muscle contraction, adrenergic signaling in cardiomyocytes, and hypertrophic cardiomyopathy ( Figure 3C ). The PPI network was composed of 77 nodes and 890 edges. Using overlap of top 40 genes from each of the five algorithms, 23 intersecting genes were identified ( Figure 3D i–ii).
MYOM2 and TCAP serve as hub diagnostic biomarkers in OA
A random forest analysis was performed on the 23 intersecting genes, and five feature genes (TNNC1, CASQ1, MYL3, MYOM2, and TCAP) were finally screened out ( Figure 4A ). The full variable importance plot is provided in Supplementary Figure S2 . In the training set, all five feature genes had an AUC greater than 0.7 ( Figure 4B ). In the testing set, CASQ1 was not identified, TNNC1 had AUC values less than 0.7, and three feature genes (MYL3, MYOM2, and TCAP) had AUC values greater than 0.7 ( Figure 4C ). Therefore, these three feature genes (MYL3, MYOM2, and TCAP) were selected for expression validation. Of these 3 feature genes, MYL3 did not show significant expression in the testing set. However, MYOM2 and TCAP were significantly differentially expressed in both sets (p < 0.05), with consistent expression trends ( Figure 4D ). Therefore, MYOM2 and TCAP were confirmed as OA hub genes.
MYOM2 and TCAP nomogram showed promising predictive performance for OA risk
Following the expression analysis of MYOM2 and TCAP, a nomogram was assembled to predict OA risk ( Figure 5A ). The calibration curve was used to further verify the accuracy of the model, with intergroup comparisons performed between uncorrected data and bias-corrected data. The results showed that p = 0.643 > 0.05, suggesting no statistically significant difference between the two groups. This finding confirms favorable result consistency before and after bias correction, and showed the satisfactory calibration capability and reliability of the nomogram ( Figure 5B ). Additionally, the ROC curve yielded an AUC of 0.875, showing a strong fit and high predictive performance for the nomogram ( Figure 5C ). These results suggested that MYOM2 and TCAP could be dependable hub genes for predicting OA risk.
MYOM2 and TCAP are enriched in neurodegenerative and cardiac pathways
GSEA revealed the functions of TCAP, showing its involvement in 5 signaling pathways ( Figure 6A ), while MYOM2 was enriched in 5 signaling pathways ( Figure 6B ). Both genes were commonly associated with pathways such as Alzheimer’s disease, Parkinson’s disease, cardiac muscle contraction, oxidative phosphorylation, and Huntington’s disease. The enrichment with MYOM2 and TCAP in these pathways suggested their potential role in OA, possibly by regulating cellular metabolism, oxidative stress response, and muscle function, thereby influencing the onset and progression of OA.
MYOM2 and TCAP correlate with skeletal muscle cells and endothelial cells in OA
A heatmap was generated to visualize the enrichment scores of 64 immune infiltrating cells in OA samples and normal controls ( Figure 7A ). Among these, 14 immune cells exhibited significant differences between the OA and control groups, and were defined as differential immune cells, such as common myeloid progenitors (CMP), endothelial cells, erythrocytes, fibroblasts, hematopoietic stem cells (HSC), immature dendritic cells (iDC), and lymphatic endothelial cells ( Figure 7B ). Notably, smooth muscle, skeletal muscle, pDC, memory B cells, and myocytes were more enriched in the OA group. Further correlation analysis revealed that TCAP had the strongest negative connection to endothelial cells (cor = −0.66, p < 0.001) and the strongest positive relationship with skeletal muscle cells (cor = 0.88, p < 0.001). Similarly, MYOM2 showed a substantial positive link to skeletal muscle cells (cor = 0.92, p < 0.001) and a notable negative connection with endothelial cells (cor = −0.54, p < 0.01) ( Figure 7C ; Table 3 ).
TABLE 3
Correlation results of immune cells.
| Gene | Cell | Correlation | P value |
| MYOM2 | CMP | −0.339,356,296 | 0.122327878 |
| TCAP | CMP | −0.178994918 | 0.425438244 |
| MYOM2 | Endothelial cells | −0.543349358 | 0.008963843 |
| TCAP | Endothelial cells | −0.664783986 | 0.000737658 |
| MYOM2 | Erythrocytes | 0.305978423 | 0.166091029 |
| TCAP | Erythrocytes | 0.269235672 | 0.225648839 |
| MYOM2 | Fibroblasts | −0.523016118 | 0.012503075 |
| TCAP | Fibroblasts | −0.37616494 | 0.084453917 |
| MYOM2 | HSC | −0.479616871 | 0.023899026 |
| TCAP | HSC | −0.500568258 | 0.017656713 |
| MYOM2 | iDC | −0.410502541 | 0.057741205 |
| TCAP | iDC | −0.304347826 | 0.16847502 |
| MYOM2 | Endothelial cells | −0.424738791 | 0.048797687 |
| TCAP | Endothelial cells | −0.551821541 | 0.007756195 |
| MYOM2 | Memory B-cells | 0.19566495 | 0.382849291 |
| TCAP | Memory B-cells | 0.114185704 | 0.612873626 |
| MYOM2 | MEP | 0.026546174 | 0.90665028 |
| TCAP | MEP | 0.220276767 | 0.324599577 |
| MYOM2 | Mv endothelial cells | −0.477131564 | 0.024743739 |
| TCAP | Mv endothelial cells | −0.575381141 | 0.005083271 |
| MYOM2 | Myocytes | 0.716310912 | 0.000177055 |
| TCAP | Myocytes | 0.784827608 | 1.53E-05 |
| MYOM2 | pDC | 0.490472102 | 0.020475709 |
| TCAP | pDC | 0.499026848 | 0.018065556 |
| MYOM2 | Skeletal muscle | 0.923240189 | 9.31E-10 |
| TCAP | Skeletal muscle | 0.879631579 | 6.95E-08 |
| MYOM2 | Smooth muscle | −0.105851193 | 0.639193945 |
| TCAP | Smooth muscle | −0.183046538 | 0.414859631 |
MELPHALAN and Testosterone Enanthate show stable binding affinity to MYOM2 and TCAP
DSigDB is a specialized database dedicated to associating drugs/compounds with target genes and can be applied to GSEA. It contains a total of 22,527 gene sets, covering 17,389 unique compounds and 19,531 genes. Constructed based on data of quantitative drug-induced inhibition and gene expression alterations, the database is classified into four datasets. In the present study, this database was utilized to carry out relevant analyses including drug-gene association analysis and drug repositioning prediction. Analysis via the DSigDB (v. 1.0) predicted 10 and 3 compounds that interact with MYOM2 and TCAP, respectively, with no shared compounds identified between the two. Among these, MELPHALAN (for MYOM2) and Testosterone Enanthate (for TCAP) exhibited the highest combined scores and were thus selected as ligands for subsequent analyses ( Figure 8A ). Molecular docking results showed favorable binding capacities of these ligands to their respective targets, with binding free energies of −6.0 kcal/mol and −5.9 kcal/mol, respectively ( Figure 8B ). To validate the stability of the docked conformations, 100 ns molecular dynamics simulations were performed. RMSD analysis revealed that the TCAP-Testosterone Enanthate complex displayed minimal structural fluctuations throughout the simulation (RMSD maintained within 0–2.5 nm), while the MYOM2-MELPHALAN complex also exhibited similar stability during the 0–80 ns period ( Figure 8C ). RMSF analysis indicated that the MYOM2-MELPHALAN complex possessed higher flexibility in the 100–125 and 160–180 residue segments. In contrast, the TCAP-Testosterone Enanthate complex showed significant fluctuations in the 0–60 residue region, suggesting this segment may serve as the ligand-binding interface ( Figure 8D ). Both complexes exhibited low total energy with gentle fluctuations, indicating stable energy states of the systems ( Figure 8E ). The number of hydrogen bonds showed periodic changes during the simulation, with the TCAP-Testosterone Enanthate and MYOM2-MELPHALAN complexes maintaining 0–2 and 0-4 hydrogen bonds, respectively. This reflects a dynamic equilibrium at the binding interface and good elasticity of the hydrogen bond network ( Figure 8F ). Collectively, the results of molecular docking and dynamics simulations support that both the MYOM2-MELPHALAN and TCAP-Testosterone Enanthate complexes possess high structural stability.
RT-qPCR shows upregulation of MYOM2 and TCAP in OA
RT-qPCR was performed to compare the mRNA expression levels of MYOM2 and TCAP between the control group and osteoarthritis model group. The results showed that there were significant differences in the expression of these two genes (p < 0.05). Specifically, the expression level of MYOM2 was significantly increased in the osteoarthritis model group compared with the control group (p < 0.05) ( Figure 9A ). Similarly, the expression of TCAP was also remarkably upregulated in the model group relative to the control group (p < 0.05) ( Figure 9B ). Collectively, these findings further verified the authenticity and reliability of the bioinformatics analysis results. In addition, it also suggests that MYOM2 and TCAP, the key target genes of arginine and proline metabolism, are closely associated with the pathogenesis of OA.
Discussion
OA is a chronic degenerative joint disease with joint pain, swelling, and even joint mobility disorders. Its course is long, the disease is prolonged and difficult to cure, and the pathogenesis is not clear. Recent studies have shown that dysregulated arginine and proline metabolism within the joint is closely associated with OA pathogenesis. Metabolomic profiling of human OA synovial fluid has revealed significantly upregulated arginine and proline metabolism compared with healthy controls ( Carlson et al., 2018 ). In OA cartilage from both human patients and mouse models, the arginine-metabolizing enzyme arginase II is specifically upregulated in chondrocytes; overexpression of arginase II in mouse joint tissues causes OA pathogenesis, whereas genetic ablation of arginase II abolishes surgically induced OA ( Choi et al., 2019 ). Concurrently, multi-omics analysis of paired OA synovial tissue and fluid has detected enrichment of proline-related metabolic pathways and elevated hydroxyproline levels, reflecting active collagen degradation ( Ge et al., 2025 ). These findings indicate that aberrant arginine and proline metabolism promotes synovial inflammation and cartilage destruction in OA. However, the specific genes through which this metabolic dysregulation drives OA pathogenesis have not been fully characterized. By using transcriptomics sequencing data and integrating with comprehensive gene expression databases for bioinformatics analysis, it has been widely applied to discover new therapeutic targets for diseases. Therefore, the focus of our research was to determine the diagnostic markers and potential therapeutic targets related to arginine and proline metabolism in OA. Through a systematic analysis, we unexpectedly discovered that MYOM2 and TCAP are key biomarkers and therapeutic targets related to arginine and proline metabolism in OA, ultimately influencing the progression of OA.
In this study, clinical synovial tissues from OA patients and normal control groups were collected for sequencing, and the differential analysis of sequencing data training sets of 12 OA patient samples and 10 normal control samples was performed. A total of 2,227 DEGs were screened out, including 1,435 upregulated genes and 792 downregulated genes, which indicated that OA was involved in multiple complex biological processes ( Chen and Zhang, 2024 ). In the GO and KEGG enrichment analysis, 96 candidate genes of arginine and proline metabolism were found to be related to a series of biological processes and signaling pathways of muscle tissues and functions, especially key functions such as muscle organ development, muscle system process, and muscle tissue development. This indicates that muscle-specific kinase (MUSK)-related genes and their functions may play an important role in maintaining the motor function and structural stability of joints during OA progression, and also suggests that abnormal muscle function may be an important factor affecting the function and structural stability of joints during OA progression ( Shorter et al., 2019 ; Krishnasamy et al., 2018 ). These results not only deepen our understanding of OA pathogenesis, but also suggest that new therapeutic strategies can be developed in the future by modulating muscle and related signaling pathways.
This study found that MYOM2 and TCAP are the key biomarkers of arginine and proline metabolism in OA, which play an important role in the occurrence and development of OA. The TCAP gene plays a key role in maintaining muscle structure and function. The TCAP protein encoded by the gene interacts with the giant Titin to regulate the elasticity and stability of muscle fibers ( Zhao et al., 2024 ). TCAP not only participates in the regulation of muscle contraction, but also plays an important role in cytoskeleton signaling and mechanical force sensing ( Lv et al., 2023 ; Knöll et al., 2011 ). For patients with OA, the injury and degeneration of muscle tissue around the joint is a common phenomenon, and the loss of muscle function may further aggravate the wear and inflammation of the joint ( Yang et al., 2023 ). Abnormal expression or dysfunction of TCAP gene may accelerate the development of OA by affecting the structural integrity and ability of muscle contraction, thereby leading to decreased joint stability ( Hodge et al., 2019 ). Therefore, TCAP gene plays an important role in the maintenance of joint and muscle function, and its regulatory role in OA deserves further study in order to provide new molecular targets for the treatment of OA. MYOM2 is essential for maintaining the structural integrity of muscle cells and responding to mechanical stress ( Woods et al., 2008 ). The MYOM2 encoded by this gene is a muscle-specific structural protein located at the M-line of muscle fibers and mainly involved in regulating the mechanical stability and elasticity of muscles ( Wang et al., 2024 ). MYOM2 is critical in maintaining the structural integrity of muscle cells and responding to mechanical stress. In OA, the muscle tissue around the joint is often damaged or degraded, and the weakening of muscle function may further aggravate the degenerative changes and inflammatory response of the joint. The abnormal expression or dysfunction of MYOM2 gene may affect the structural stability of muscle fibers, leading to the dysfunction of muscles around the joint. Thereby it reduces the ability of the joint to support and accelerates cartilage damage and joint degradation ( Liu et al., 2019 ). The MYOM2 gene plays an important role in regulating muscle and joint stability, and its potential regulatory mechanisms may provide new molecular targets for OA intervention.
The shared enrichment of TCAP and MYOM2 in neurodegenerative disease pathways suggests that these hub genes may participate in neurobiological processes extending beyond the joint. Emerging evidence indicates that OA involves not only local joint degeneration but also peripheral and central nervous system sensitization. Within OA joints, nociceptive afferent neurons are sensitized by neuropeptides, neurotrophins, and pro-inflammatory cytokines released from inflamed synovium and degrading cartilage, leading to chronic pain signaling ( Batchelor et al., 2025 ). Beyond peripheral sensitization, a growing body of literature has established epidemiological and mechanistic links between OA and central neuroinflammation. A recent review demonstrated that OA-associated systemic inflammation can compromise blood-brain barrier integrity and exacerbate neuroinflammatory processes, while shared mechanisms such as oxidative stress and mitochondrial dysfunction connect OA with neurodegenerative conditions including Alzheimer’s and Parkinson’s diseases ( Rabie et al., 2025 ). The co-enrichment of MYOM2 and TCAP in Alzheimer’s, Parkinson’s, and Huntington’s disease pathways identified in our GSEA analysis may therefore reflect these shared cellular stress and metabolic dysregulation mechanisms at the systems level. Importantly, central sensitization and neuropathic-like pain have been recognized as key drivers of chronic OA pain, with studies showing that central sensitization-related signs are present in a substantial proportion of OA patients and are independent predictors of poor functional outcomes ( Salaffi et al., 2025 ). This has therapeutic implications, as interventions targeting central pain mechanisms rather than peripheral nociception alone are increasingly recognized as necessary components of comprehensive OA management. Furthermore, the enrichment of TCAP and MYOM2 in cardiac muscle contraction pathways, together with their predominant expression in skeletal muscle and heart ( GTEx Consortium., 2020 ), raises the possibility that these genes participate in a broader musculoskeletal-cardiovascular regulatory network that may be coordinately dysregulated in OA, a disease known to share risk factors and comorbidities with cardiovascular disorders. Further investigation into whether MYOM2 and TCAP modulation affects neural sensitization pathways or muscle-joint-neural crosstalk may reveal novel therapeutic strategies targeting both structural joint pathology and chronic pain in OA. Analysis of immune cell infiltration also revealed that TCAP and MYOM2 showed significant correlations with endothelial cells and skeletal muscle cells, which was consistent with the GSEA enrichment results. These findings further support the potential involvement of these two genes in skeletal muscle and cardiovascular functions during OA pathogenesis, thereby providing new directions for mechanistic investigations of this disease.
This study collected clinical samples from OA patients and normal controls, and performed quantitative RT-qPCR analysis on key targets. The results showed that MYOM2 and TCAP were significantly upregulated in OA synovial tissues, further suggesting the important roles of the screened MYOM2 and TCAP in OA progression. While these findings reinforce the potential of MYOM2 and TCAP as biomarkers for OA, their specificity requires further validation. Although our GSEA analysis revealed that these genes are co-enriched with Alzheimer’s and Huntington’s disease pathways, GTEx data demonstrate that both MYOM2 and TCAP are predominantly expressed in skeletal muscle and heart, with minimal expression detected in the central nervous system ( GTEx Consortium., 2020 ). This pathway-level co-enrichment more likely reflects shared cellular stress and metabolic dysfunction mechanisms rather than genuine expression of these genes in neural tissues. Another consideration is whether MYOM2 and TCAP are similarly dysregulated in other inflammatory joint diseases, as this would directly affect their diagnostic specificity. Currently, publicly available datasets do not provide direct comparisons of MYOM2 or TCAP expression in synovial tissues from patients with other acute inflammatory arthropathies. However, expression data from non-synovial sources offer informative contrasts. RNA sequencing of peripheral blood mononuclear cells from patients with ankylosing spondylitis revealed that MYOM2 is significantly downregulated compared with healthy controls ( Huang et al., 2020 ), a direction of change opposite to the upregulation observed in OA synovial tissue in the present study. This divergent expression pattern suggests that MYOM2 dysregulation is disease-context-dependent rather than a universal feature of inflammatory arthritis. Furthermore, a single nucleotide polymorphism of MYOM2 has been reported to be significantly associated with responsiveness to etanercept treatment in ankylosing spondylitis patients ( Liu et al., 2019 ). Collectively, the existing evidence supports the potential of MYOM2 and TCAP as relatively specific indicators of OA-related joint pathology. Systematic profiling of these genes across synovial tissues from multiple inflammatory joint diseases remains an important direction for future validation.
The study also predicted potential therapeutic agents targeting two hub genes. The analysis revealed that MYOM2 may potentially interact with 10 small-molecule drugs, among which MELPHALAN achieved the highest comprehensive score for MYOM2. Molecular docking and molecular dynamics simulations showed its strong target binding affinity and structural stability. Clinically, MELPHALAN is a nitrogen mustard alkylating agent used in the treatment of multiple myeloma and ovarian cancer. Its antitumor activity is mediated through DNA interstrand cross-linking at the N7 position of guanine, leading to inhibition of DNA replication and induction of apoptosis in rapidly dividing cells ( Esma et al., 2017 ). Notably, the biological effects of MELPHALAN are dose-dependent: at concentrations approximately 100-fold lower than cytostatic doses, MELPHALAN loses its DNA-damaging activity and instead disrupts signal transduction through TNF-α and IL-2β receptors, thereby exerting anti-inflammatory effects independent of its alkylating properties, as evidenced by reduced severity of experimental colitis in mice following non-cytotoxic dosing ( Pukhalsky et al., 2006 ); conversely, at cytotoxic doses MELPHALAN induces oxidative stress, inflammation, and apoptosis through modulation of pro-inflammatory cytokines and apoptotic markers ( Panghal et al., 2023 ). These findings collectively establish that MELPHALAN can exert dose-dependent regulation of cell proliferation, apoptosis, and inflammatory signaling, suggesting that MYOM2 in non-oncological conditions such as OA may be involved in analogous regulatory pathways. The analysis also indicated that TCAP is associated with three drugs, primarily involved in hormonal regulation, anti-inflammatory effects, and metabolic pathways. Among them, Testosterone Enanthate scored the highest for TCAP. As a hormonal regulator, it may potentially modulate muscle function and promote tissue repair, which holds significant implications for the clinical treatment of OA ( Howard et al., 2022 ). The construction of a drug-gene interaction network revealed that regulating pathways associated with MYOM2 and TCAP could lead to the development of novel therapeutic strategies for OA. The association between MYOM2 and chemotherapeutic agents suggests that inhibiting specific cellular metabolic and proliferative pathways may influence the pathological progression of OA. Meanwhile, the link between TCAP and anti-inflammatory or hormone-regulating drugs implies that modulating joint inflammation and tissue repair through these pathways could offer more targeted treatment options for patients. We employed computational drug prediction and molecular binding validation to assess the pharmacological relevance of these targets. These findings provide compelling new leads for the development of more effective OA therapies, with the potential to reduce drug development costs and advance personalized medicine. This study contributes to the field by highlighting the significance of these promising therapeutic targets in OA treatment.
Although our study provided valuable insights, several limitations should be acknowledged. Our basic experiments focused on validation at the gene expression level and did not include functional assays. Furthermore, the molecular docking and dynamics simulations conducted in this study rely heavily on the quality of the three-dimensional structures of the proteins and ligands. Although such computational analyses can offer important clues for drug target identification, they are insufficient to fully predict the actual effects of these targets in complex physiological environments. Subsequent validation—from in vitro and in vivo studies to eventual clinical trials—is essential to confirm the therapeutic potential of these targets. Such efforts will help to better understand the roles of arginine and proline metabolism in OA and facilitate the clinical translation of related research findings.
In conclusion, this study collected clinical samples for transcriptome sequencing and combined it with multiple bioinformatics techniques to identify TCAP and MYOM2 as hub genes involved in arginine and proline metabolism in OA. These genes are key biomarkers and potential therapeutic targets, and these targets were further validated through in vitro experiments. Furthermore, pharmacological prediction, molecular docking, and molecular dynamics simulations not only verified the feasibility of these key targets at the pharmacological level but also screened promising lead compounds. These findings provide candidate drugs and a theoretical framework for interventional therapy of OA, holding significant potential for accelerating the development of cost-effective and highly efficient personalized treatment strategies. This study may offer new directions for precision medicine in OA.
Acknowledgements
We sincerely express our gratitude to the team led by Shengjian Weng from the Orthopedics Department of the Fuzhou Second General Hospital for generously providing us with human synovial tissue for this research.
Edited by: Felipe Eltit Guersetti , Vancouver Prostate Centre, Canada
Reviewed by: Robyn Birch , University of British Columbia, Canada
Mauricio Quezada , Finis Terrae University, Chile
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material .
Ethics statement
The studies involving humans were approved by Fuzhou Second General Hospital Medical Ethics Committee (Ethics Number:2021014). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
X-HC: Writing – review and editing, Writing – original draft, Conceptualization, Methodology. JL: Writing – original draft, Data curation, Conceptualization, Writing – review and editing. LQ: Investigation, Writing – original draft. YZ: Validation, Writing – original draft. Z-MZ: Writing – review and editing, Software. PC: Writing – review and editing, Methodology. Y-XS: Writing – original draft. J-MG: Project administration, Writing – review and editing, Supervision, Writing – original draft. S-JW: Writing – review and editing, Project administration, Supervision, Writing – original draft.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2026.1795267/full#supplementary-material
References
Abramoff B. Caldera F. E. ( 2020 ). Osteoarthritis: pathology, diagnosis, and treatment options . Med. Clin. North Am. 104 ( 2 ), 293 – 311 . 10.1016/j.mcna.2019.10.007 32035570
Alderden J. Pepper G. A. Wilson A. Whitney J. D. Richardson S. Butcher R. ( 2018 ). Predicting pressure injury in critical care patients: a machine-learning model . Am. J. Crit. Care 27 ( 6 ), 461 – 468 . 10.4037/ajcc2018525 30385537 PMC6247790
Aran D. Hu Z. Butte A. J. ( 2017 ). xCell: digitally portraying the tissue cellular heterogeneity landscape . Genome Biol. 18 ( 1 ), 220 . 10.1186/s13059-017-1349-1 29141660 PMC5688663
Batchelor V. Perry T. A. Cader M. Z. Vincent T. L. ( 2025 ). Peripheral neuronal sensitization and neurovascular remodelling in osteoarthritis pain . Nat. Rev. Rheumatol. 21 ( 9 ), 526 – 545 . 10.1038/s41584-025-01280-3 40796685
Bourget-Murray J. Piroozfar S. Smith C. Ellison J. Bansal R. Sharma R. ( 2023 ). Annual incidence and assessment of risk factors for early-onset deep surgical site infection following primary total knee arthroplasty in osteoarthritis . Bone Jt. J. 105-B ( 9 ), 971 – 976 . 10.1302/0301-620X.105B9 37654121
Brown J. Pirrung M. McCue L. A. ( 2017 ). FQC dashboard: integrates FastQC results into a web-based, interactive, and extensible FASTQ quality control tool . Bioinformatics 33 ( 19 ), 3137 – 3139 . 10.1093/bioinformatics/btx373 28605449 PMC5870778
Carlson A. K. Rawle R. A. Adams E. Greenwood M. C. Bothner B. June R. K. ( 2018 ). Application of global metabolomic profiling of synovial fluid for osteoarthritis biomarkers . Biochem. Biophys. Res. Commun. 499 ( 2 ), 182 – 188 . 10.1016/j.bbrc.2018.03.117 29551687 PMC5890938
Chaganti R. K. Lane N. E. ( 2011 ). Risk factors for incident osteoarthritis of the hip and knee . Curr. Rev. Musculoskelet. Med. 4 ( 3 ), 99 – 104 . 10.1007/s12178-011-9088-5 21808997 PMC3261259
Chen S. Zhang C. ( 2024 ). Effects of ischaemic post-conditioning on eccentric exercise-induced muscle damage . Biol. Sport 41 ( 2 ), 27 – 35 . 10.5114/biolsport.2024.129483 38524812 PMC10955728
Chen R. Han S. Liu X. Wang K. Zhou Y. Yang C. ( 2018 ). Perturbations in amino acids and metabolic pathways in osteoarthritis patients determined by targeted metabolomics analysis . J. Chromatogr. B Anal. Technol. Biomed. Life Sci. 1085 , 54 – 62 . 10.1016/j.jchromb.2018.03.047 29631251
Choi W. S. Yang J. I. Kim W. Kim H. E. Kim S. K. Won Y. ( 2019 ). Critical role for arginase II in osteoarthritis pathogenesis . Ann. Rheum. Dis. 78 ( 3 ), 421 – 428 . 10.1136/annrheumdis-2018-214282 30610061 PMC6390026
Esma F. Salvini M. Troia R. Boccadoro M. Larocca A. Pautasso C. ( 2017 ). Melphalan hydrochloride for the treatment of multiple myeloma . Expert Opin. Pharmacother. 18 ( 11 ), 1127 – 1136 . 10.1080/14656566.2017.1349102 28658983
Feng N. Ye Y. Pan Y. Kuang B. Du Y. Geng N. ( 2024 ). The circUbqln1, regulated by XBP1s, interplays with 14-3-3ζ to inhibit collagen synthesis and promote osteoarthritis by controlling PRODH activity and proline metabolism . J. Adv. Res. 66 , 267 – 284 . 10.1016/j.jare.2024.01.007 38219870 PMC11674786
Ge M. Sun W. Xu T. Yang R. Zhang K. Li J. ( 2025 ). Multi-omics analysis of synovial tissue and fluid reveals differentially expressed proteins and metabolites in osteoarthritis . J. Transl. Med. 23 ( 1 ), 285 . 10.1186/s12967-025-06310-y 40050855 PMC11887128
GTEx Consortium ( 2020 ). The GTEx consortium atlas of genetic regulatory effects across human tissues . Science 369 ( 6509 ), 1318 – 1330 . 10.1126/science.aaz1776 32913098 PMC7737656
Gu Z. Hübschmann D. ( 2022 ). Make interactive complex heatmaps in R . Bioinformatics 38 ( 5 ), 1460 – 1462 . 10.1093/bioinformatics/btab806 34864868 PMC8826183
Gustavsson E. K. Zhang D. Reynolds R. H. Garcia-Ruiz S. Ryten M. ( 2022 ). Ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2 . Bioinformatics 38 ( 15 ), 3844 – 3846 . 10.1093/bioinformatics/btac409 35751589 PMC9344834
Hänzelmann S. Castelo R. Guinney J. ( 2013 ). GSVA: gene set variation analysis for microarray and RNA-seq data . BMC Bioinforma. 14 , 7 . 10.1186/1471-2105-14-7 23323831 PMC3618321
Hodge B. A. Zhang X. Gutierrez-Monreal M. A. Cao Y. Hammers D. W. Yao Z. ( 2019 ). MYOD1 functions as a clock amplifier as well as a critical co-factor for downstream circadian gene expression in muscle . Elife 8 , e43017 . 10.7554/eLife.43017 30789342 PMC6398978
Howard E. E. Shankaran M. Evans W. J. Berryman C. E. Margolis L. M. Lieberman H. R. ( 2022 ). Effects of testosterone on mixed-muscle protein synthesis and proteome dynamics during energy deficit . J. Clin. Endocrinol. Metab. 107 ( 8 ), e3254 – e3263 . 10.1210/clinem/dgac295 35532889
Huang D. Liu J. Cao Y. Wan L. Jiang H. Sun Y. ( 2020 ). RNA sequencing for gene expression profiles in peripheral blood mononuclear cells with ankylosing spondylitis RNA . Biomed. Res. Int. 2020 , 5304578 . 10.1155/2020/5304578 32596323 PMC7298317
Karna E. Szoka L. Huynh T. Y. L. Palka J. A. ( 2020 ). Proline-dependent regulation of collagen metabolism . Cell Mol. Life Sci. 77 ( 10 ), 1911 – 1918 . 10.1007/s00018-019-03363-3 31740988 PMC7228914
Knöll R. Linke W. A. Zou P. Miocic S. Kostin S. Buyandelger B. ( 2011 ). Telethonin deficiency is associated with maladaptation to biomechanical stress in the mammalian heart . Circ. Res. 109 ( 7 ), 758 – 769 . 10.1161/CIRCRESAHA.111.245787 21799151 PMC3664427
Krishnasamy P. Hall M. Robbins S. R. ( 2018 ). The role of skeletal muscle in the pathophysiology and management of knee osteoarthritis . Rheumatol. Oxf. 57 ( Suppl. l_4 ), iv124 . 10.1093/rheumatology/key039 29444285
Langfelder P. Horvath S. ( 2008 ). WGCNA: an R package for weighted correlation network analysis . BMC Bioinforma. 9 , 559 . 10.1186/1471-2105-9-559 19114008 PMC2631488
Li L. Shao C. Liu Z. Wu X. Yang J. Wan H. ( 2022 ). Comparative efficacy of Honghua class injections for treating acute ischemic stroke: a Bayesian network meta-analysis of randomized controlled trials . Front. Pharmacol. 13 , 1010533 . 10.3389/fphar.2022.1010533 36249799 PMC9554475
Liu J. Zhu Q. Han J. Zhang H. Li Y. Ma Y. ( 2019 ). IgG Galactosylation status combined with MYOM2-rs2294066 precisely predicts anti-TNF response in ankylosing spondylitis . Mol. Med. 25 ( 1 ), 25 . 10.1186/s10020-019-0093-2 31195969 PMC6567531
Lv X. Lin F. Wu W. Wang H. Luo Y. Wang Z. ( 2023 ). The clinical features and TCAP mutation spectrum in a Chinese cohort of patients with limb-girdle muscular dystrophy R7 . Hum. Mol. Genet. 32 ( 15 ), 2502 – 2510 . 10.1093/hmg/ddad090 37216648
Molnar V. Matišić V. Kodvanj I. Bjelica R. Jeleč Ž. Hudetz D. ( 2021 ). Cytokines and chemokines involved in osteoarthritis pathogenesis . Int. J. Mol. Sci. 22 ( 17 ), 9208 . 10.3390/ijms22179208 34502117 PMC8431625
Pang X. P. Hershman J. M. Chung M. Pekary A. E. ( 1989 ). Characterization of tumor necrosis factor-alpha receptors in human and rat thyroid cells and regulation of the receptors by thyrotropin . Endocrinology 125 ( 4 ), 1783 – 1788 . 10.1210/endo-125-4-1783 2551628
Panghal A. Kumar V. Jena G. ( 2023 ). Melphalan induced germ cell toxicity and dose-dependent effects of β-aminoisobutyric acid in experimental rat model: role of oxidative stress, inflammation and apoptosis . J. Biochem. Mol. Toxicol. 37 ( 8 ), e23374 . 10.1002/jbt.23374 37086025
Pukhalsky A. Shmarina G. Alioshkin V. Sabelnikov A. ( 2006 ). Alkylating drugs applied in non-cytotoxic doses as a novel compounds targeting inflammatory signal pathway . Biochem. Pharmacol. 72 ( 11 ), 1432 – 1438 . 10.1016/j.bcp.2006.03.008 16620792
Quicke J. G. Conaghan P. G. Corp N. Peat G. ( 2022 ). Osteoarthritis year in review 2021: epidemiology and therapy . Osteoarthr. Cartil. 30 ( 2 ), 196 – 206 . 10.1016/j.joca.2021.10.003 34695571
Rabie M. A. Madry H. Cucchiarini M. El-Sayed N. S. ( 2025 ). The brain-joint axis: links between osteoarthritis and neurodegenerative disorders in aging . J. Adv. Res. S2090-1232 ( 25 ). 10.1016/j.jare.2025.10.023 41130360
Robin X. Turck N. Hainard A. Tiberti N. Lisacek F. Sanchez J. C. ( 2011 ). pROC: an open-source package for R and S+ to analyze and compare ROC curves . BMC Bioinforma. 12 , 77 . 10.1186/1471-2105-12-77 21414208 PMC3068975
Robles-Jimenez L. E. Aranda-Aguirre E. Castelan-Ortega O. A. Shettino-Bermudez B. S. Ortiz-Salinas R. Miranda M. ( 2021 ). Worldwide traceability of antibiotic residues from livestock in wastewater and soil: a systematic review . Animals (Basel) 12 ( 1 ), 60 . 10.3390/ani12010060 35011166 PMC8749557
Sachs M. C. ( 2017 ). plotROC: a tool for plotting ROC curves . J. Stat. Softw. 79 , 2 . 10.18637/jss.v079.c02 30686944 PMC6347406
Salaffi F. Carotti M. Farah S. Ciccullo C. Gigante A. P. Bandinelli F. ( 2025 ). A mediation appraisal of neuropathic-like symptoms, pain catastrophizing, and central sensitization-related signs in adults with knee Osteoarthritis-A cross-sectional study . J. Pers. Med. 15 ( 1 ), 22 . 10.3390/jpm15010022 39852214 PMC11767182
Shannon P. Markiel A. Ozier O. Baliga N. S. Wang J. T. Ramage D. ( 2003 ). Cytoscape: a software environment for integrated models of biomolecular interaction networks . Genome Res. 13 ( 11 ), 2498 – 2504 . 10.1101/gr.1239303 14597658 PMC403769
Sheng L. Liu X. Zhang H. Qian X. Gu Y. Zhu W. ( 2025 ). Identification of biomarkers related to neutrophil extracellular traps and potential therapeutic drugs for rheumatoid arthritis using computational analysis . Eur. J. Med. Res. 30 ( 1 ), 965 . 10.1186/s40001-025-03213-y 41088461 PMC12522263
Shorter E. Sannicandro A. J. Poulet B. Goljanek-Whysall K. ( 2019 ). Skeletal muscle wasting and its relationship with osteoarthritis: a mini-review of mechanisms and current interventions . Curr. Rheumatol. Rep. 21 ( 8 ), 40 . 10.1007/s11926-019-0839-4 31203463 PMC6571089
Trapnell C. Roberts A. Goff L. Pertea G. Kim D. Kelley D. R. ( 2012 ). Differential gene and transcript expression analysis of RNA-Seq experiments with TopHat and Cufflinks . Nat. Protoc. 7 ( 3 ), 562 – 578 . 10.1038/nprot.2012.016 22383036 PMC3334321
Walter W. Sánchez-Cabo F. Ricote M. ( 2015 ). GOplot: an R package for visually combining expression data with functional analysis . Bioinformatics 31 ( 17 ), 2912 – 2914 . 10.1093/bioinformatics/btv300 25964631
Wang Z. Zhang J. Shi S. Ma H. Wang D. Zuo C. ( 2023 ). Predicting lung adenocarcinoma prognosis, immune escape, and pharmacomic profile from arginine and proline-related genes . Sci. Rep. 13 ( 1 ), 15198 . 10.1038/s41598-023-42541-z 37709932 PMC10502151
Wang X. Zhang J. Su J. Huang T. Lian L. Nie Q. ( 2024 ). Genome-wide mapping of the binding sites of myocyte enhancer factor 2A in chicken primary myoblasts . Poult. Sci. 103 ( 10 ), 104097 . 10.1016/j.psj.2024.104097 39094502 PMC11345569
Wang S. Ren Q. Huang Y. Ma Y. Zhang X. Liu Y. ( 2026 ). Emerging roles of arginine metabolism in skeletal health and disease . Metabolism 175 , 156451 . 10.1016/j.metabol.2025.156451 41270967
Woods S. Farrall A. Procko C. Whitelaw M. L. ( 2008 ). The bHLH/Per-Arnt-Sim transcription factor SIM2 regulates muscle transcript myomesin2 via a novel, non-canonical E-box sequence . Nucleic Acids Res. 36 ( 11 ), 3716 – 3727 . 10.1093/nar/gkn247 18480125 PMC2441813
Wu T. Hu E. Xu S. Chen M. Guo P. Dai Z. ( 2021 ). clusterProfiler 4.0: a universal enrichment tool for interpreting omics data . Innovation (Camb) 2 ( 3 ), 100141 . 10.1016/j.xinn.2021.100141 34557778 PMC8454663
Yang J. Liu P. Wang S. Jiang T. Zhang Y. Liu W. ( 2023 ). Causal relationship between sarcopenia and osteoarthritis: a bi-directional two-sample mendelian randomized study . Eur. J. Med. Res. 28 ( 1 ), 327 . 10.1186/s40001-023-01322-0 37689698 PMC10492359
Yuan C. Pan Z. Zhao K. Li J. Sheng Z. Yao X. ( 2020 ). Classification of four distinct osteoarthritis subtypes with a knee joint tissue transcriptome atlas . Bone Res. 8 ( 1 ), 38 . 10.1038/s41413-020-00109-x 33298863 PMC7658991
Zhao Y. Liang J. Liu X. Li H. Chang C. Gao P. ( 2024 ). Tcap deficiency impedes striated muscle function and heart regeneration with elevated ROS and autophagy . Biochim. Biophys. Acta Mol. Basis Dis. 1870 ( 8 ), 167485 . 10.1016/j.bbadis.2024.167485 39226992
Zheng Y. Gao W. Zhang Q. Cheng X. Liu Y. Qi Z. ( 2022 ). Ferroptosis and autophagy-related genes in the pathogenesis of ischemic cardiomyopathy . Front. Cardiovasc. Med. 9 , 906753 . 10.3389/fcvm.2022.906753 35845045 PMC9279674