PDF
Citation
Print
INTRODUCTION
Osteoarthritis (OA) is the most prevalent joint disease, causing joint pain, swelling, and permanent damage [1]. It is characterized by pathological changes in cartilage, bone, synovium, ligaments, muscles, and periarticular fat, resulting in joint dysfunction, pain, stiffness, limited function, and loss of function. Risk factors include age (33% of individuals over 75 years old have symptoms and radiographic knee OA), gender, obesity, genetics, and severe joint injuries [1]. The worldwide prevalence of OA increased by 113.25%, from 247.51 million cases in 1990 to 527.81 million cases in 2019 [2]. In recent decades, OA has become a significant public health concern globally, highlighting the importance of prevention and early treatment in alleviating its growing burden. The natural aging process is a major contributing factor to the inability of joints to repair themselves and maintain their optimal health. A study has shown a correlation between chondrocytes aging and an increase in the number of aging chondrocytes within OA cartilage [3]. These aging cells release mediators that damage joint tissue [4].
Telomere shortening and damage are well-established causes of cellular aging and senescence [5]. Several human diseases associated with normal aging are linked to accelerated telomere dysfunction, making this a crucial area of study in human pathology [6]. Previous research has indicated a potential connection between OA and telomeres, suggesting that telomere-induced cellular aging may play a significant role in the OA process. Chen et al. [7] discovered that metal exposure can accelerate cellular aging by impacting telomere length (TL), contributing to systemic aging and ultimately exacerbating the OA process. Guillén et al. [8] found that in patients with knee OA, telomeres in peripheral blood leukocytes decay more rapidly, which may reflect systemic accelerated aging and worsen cartilage degeneration. Additionally, the severity of pain in OA patients has been linked to shorter TL in white blood cells [9]. Lou et al. [10] demonstrated that paclitaxel can inhibit TNF-α-induced increase in telomerase activity and expression of p16 and p53, thereby potentially serving as a safe therapeutic option for preventing OA-related chondrocyte aging and oxidative stress. Telomeric silencing 1-like (DOT1L) histone methyltransferase disruptor acts as the primary protector of cartilage health by limiting excessive activation of the Wnt pathway. Cartilage-specific homozygous DOT1L knockout mice display severe growth phenotypes and perinatal death, indicating the importance of maintaining or enhancing DOT1L activity during aging or after injury to prevent the onset and progression of OA [11]. While it is evident that telomeres play a crucial role in the onset and progression of OA, the underlying mechanisms are not yet fully understood.
This study used bioinformatics to screen for differentially expressed telomere related genes (DETMRGs) in OA employing least absolute shrinkage and selection operator (LASSO), support vector machine-recursive feature elimination (SVM-RFE), and RandomForest algorithms. The diagnostic value was assessed using receiver operating characteristic (ROC) curves, gene expression, and nomogram. Additionally, CIBERSORT and xcell algorithms were applied to investigate the correlation between feature genes and immune infiltration, while the impact of feature genes on phenotype was evaluated through gene set enrichment analysis (GSEA). Finally, potential drugs targeting feature genes were predicted, offering a new avenue for the early diagnosis and treatment of OA.
METHODS
Datasets
GSE51588, GSE12021 and GSE55457 microarray datasets of OA, GSE77298 and GSE1919 microarray datasets of rheumatoid arthritis (RA) were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) data repository. GSE51588 dataset contains 10 normal samples and 40 OA samples; GSE12021 dataset contains 9 normal samples and 10 OA samples; GSE55457 dataset contains 10 normal samples and 10 OA samples. GSE12021 and GSE55457 datasets were merged as the validation set. GSE77298 dataset contains 7 normal samples and 14 RA samples; GSE1919 dataset contains 5 normal samples and 5 RA samples. Feature genes related to telomeres were acquired from the data repository (http://www.cancertelsys.org/telnet/) [12].
Analysis of DETMRGs
Gene differential analysis of OA between normal and OA samples was conducted using limma package (|log fold change (FC)|> 1, adjust p-value < 0.05), then we drew volcano and heatmaps of differentially expressed genes (DEGs), with the heat map only displaying the top 20 DEGs in |log FC|. DETMRGs were obtained by intersection of DEGs with DETMRGs. We then calculated the correlation between the DETMRGs and drew a heatmap. Gene ontology (GO) [13] and Kyoto encyclopedia of genes and genomes (KEGG) [14] pathway enrichment analysis of the DETMRGs were acquired by using the R software clusterProfile package (p < 0.05). A protein and protein interaction (PPI) network of DETMRGs (confidence > 0.4) was constructed by using the STRING database, then visualized it using Cytoscape v3.10.2 [15,16].
Screening feature genes based on machine learning
To identify feature genes pertinent to our research objectives and to streamline model complexity, we utilized three distinct feature selection techniques: LASSO regression, SVM-RFE analysis, and Random Forest analysis. Each method underwent evaluation through 10-fold cross-validation to ascertain the model’s stability and generalizability.
For LASSO regression analysis, we employed the R package glmnet, using candidate gene expression data as input variables, with the target variable indicating whether the sample is an OA sample. We implemented 10-fold cross-validation to determine the optimal penalty parameter λ, ensuring a balanced model that considers the trade-off between error and feature count. This process allowed us to filter out genes most significantly associated with the target classification while eliminating genes that exhibited strong correlation or minimal contribution to the model, thereby reducing complexity.
For SVM-RFE regression analysis, we utilized the R package caret. SVM-RFE is a feature selection approach that iteratively removes features that contribute less to classification based on the weight information derived from SVM. We trained candidate genes using the SVM model and computed feature weights for each gene. In each iteration, the gene with the least contribution was discarded, and the training and evaluation process was repeated until the feature set stabilized. We assessed feature stability and model performance through 10-fold cross-validation.
To construct a Random Forest model, we employed the R package randomForest, training the model with candidate gene expression data and calculating the Mean Decrease Gini value for each gene. This metric indicates the significance of genes in classification tasks, with higher values denoting greater contributions to the model. We selected the top 10 genes based on the Mean Decrease Gini ranking as critical features. Finally, we utilized 10-fold cross-validation to evaluate the model's universality in screening gene sets, ensuring the robustness and accuracy of the classification model.
GSEA
Pathway enrichment analysis on feature genes was performed by R package clusterProfiler and the data repository (c2.cp.v2023.2.Hs.symbols.gmt) was used as the reference dataset. To minimize false positive outcomes, we established a screening threshold: q.value < 0.25 and |NES| > 1. The top 5 pathways with the most significant p-value were displayed.
ROC curves and expression level validation
The ROC curve serves as a tool to assess the predictive efficacy of feature gene expression levels in differentiating outcome categories between training and validation datasets. This analysis was conducted using the R package pROC. The expression levels of feature genes, selected through machine learning techniques, were utilized as independent variables within the logistic regression model. The model was trained on the training dataset, with OA samples designated as the outcome variable. For each feature gene, the predicted probabilities were plotted against the actual binary outcomes to generate an ROC curve. The area under the curve (AUC) was computed to gauge the effectiveness of the expression levels in distinguishing between the two categories; a higher AUC indicates a greater discriminative capacity. ROC curves were constructed for both the training set and the validation set to verify that the selected feature genes maintain consistent predictive power across different datasets. To assess the significance of the differences in feature gene expression levels between the two outcome groups, we employed the Wilcoxon test. This non-parametric test was chosen due to its lack of assumptions regarding the normality of data, making it appropriate for comparing the medians of two independent groups. Wilcoxon tests were conducted for each feature gene in both the training and validation sets, with p-values calculated to determine the statistical significance of the differences in expression levels.
Immune cell infiltration analysis
To assess immune cell infiltration and accumulation in OA and normal tissue samples, we employed two well-established algorithms: CIBERSORT and xCell. These algorithms estimate the relative abundance of different immune cell types in tissue samples using gene expression data. Both methodologies were implemented via the IOBR package in R, offering an efficient interface for these analyses.
The CIBERSORT algorithm is a deconvolution technique that utilizes gene expression data to estimate the proportions of different immune cell types within a mixed sample. It operates by comparing the expression profile of the sample against a reference gene expression matrix of known immune cell types. For each sample analyzed, the output includes estimated percentages of various immune cell populations. We utilized the default LM22 feature matrix, specifically tailored for human immune cell populations, for our analysis. The xCell algorithm represents another gene set-based quantification method for immune cell types, leveraging gene expression data to estimate the relative abundance of 64 immune and stromal cell types. In contrast to CIBERSORT, which relies on linear regression, xCell applies a gene set enrichment approach to evaluate immune cell infiltration. We applied the xCell method to further corroborate the findings from CIBERSORT and to ensure the robustness of our immune infiltration analysis. For both algorithms, immune cell infiltration scores were generated for each sample, reflecting the relative abundance of each immune cell type present. The IOBR package was utilized to automate the processing and generation of these scores. Subsequently, we performed correlation analysis to investigate the relationship between differential immune cell populations and the characteristic genes of interest. This step aimed to elucidate the connections between specific immune cell types and the expression of OA-related feature genes. We calculated the correlations between differential immune cell types and the associations between immune cell types and characteristic genes, employing the ggcorrplot package in R to visualize the correlation results.
Construction of diagnostic model
We utilized feature genes identified through machine learning techniques as inputs to develop a diagnostic model. The purpose of constructing these diagnostic models was to assess the predictive capability of the selected feature genes in differentiating OA from normal samples. This process includes the creation of column charts, ROC curve analysis, and decision curve analysis (DCA). The rms package in R was employed to create a nomogram, which estimated the probability of specific outcomes, such as the likelihood of OA, based on the values of the selected feature genes. The rms package was also used to fit logistic regression models with the characteristic genes and to generate nomographs using the model’s coefficients. The nomograph was developed by assigning scores to each feature gene according to its regression coefficient, and the total score was then utilized to predict outcomes. To assess the diagnostic model’s performance, ROC curves were generated using the pROC package. The ROC curve was essential for evaluating the model’s sensitivity and specificity across various thresholds. Additionally, to further analyze the clinical utility of the diagnostic model, DCA was conducted using the RMDA package. DCA evaluated the balance between false positives and false negatives at different decision thresholds, providing a comprehensive understanding of the net benefits of the diagnostic models.
Drug prediction
To forecast potential pharmaceuticals linked to specific genes, we utilized the Drug Gene Interaction Database (DGIdb, https://www.dgidb.org/), a comprehensive platform that consolidates drug-gene interaction information from diverse sources. The DGIdb interface offers tools for pinpointing existing medications or compounds that engage with target genes, which may act as prospective therapeutic agents for conditions like OA. The feature genes selected from the diagnostic model were entered into the DGIdb database. We extracted data to compile a list of drugs anticipated to interact with each feature gene. To enhance our understanding of the connections between characteristic genes and potential pharmaceuticals, we developed interaction networks for each gene of interest. These networks illustrate the direct or indirect interactions between each drug and the genes. By employing Cytoscape (v3.10.2) software, we visualized the relationships between characteristic genes and candidate drugs, emphasizing potential treatment avenues. These networks allow us to determine whether certain drugs engage with multiple characteristic genes or if specific genes are targeted by various drugs.
Construction of competing endogenous RNA (ceRNA) network
The development of a competitive ceRNA network is focused on investigating the regulatory relationships among lncRNA, miRNA, and differentially expressed target mRNA. According to the ceRNA hypothesis, lncRNA can modulate mRNA expression by sequestering miRNA, thus inhibiting miRNA's interaction with its target mRNA. Initially, the multiMiR package in R was utilized to identify miRNAs that may engage with mRNA. This package consolidates miRNA target interaction data from various databases, yielding a comprehensive list of predicted miRNAs capable of regulating mRNA. The interaction data between lncRNA and miRNA is sourced from the Encyclopedia of RNA Interactions (ENCORI) database (https://rnasysu.com/encori/), which offers experimentally validated interactions between lncRNAs and miRNAs. To ensure the reliability of lncRNA-miRNA interactions, we filtered lncRNAs based on the criterion of clipExpNum > 10. After identifying the miRNAs that interact with mRNA and lncRNA, we constructed a ceRNA network. In this network, mRNA, lncRNA, and miRNA are depicted as nodes, while the edges illustrate their interactions. The ceRNA network was visualized using Cytoscape v3.10.2.
Quantitative real time polymerase chain reaction (qRT-PCR) validation of the feature genes
The human chondrocyte cells (CP-H107, Procell) and human chondrocytes cells from patients with OA (402OAK-05a, BioVector NTCC Inc.) were cultured in DMEM/F12 medium (Corning). The medium was supplemented with 10% (v/v) fetal bovine serum (Sigma-Aldrich) and 1% (v/v) penicillin/streptomycin (Procell). The cells were then incubated at 37°C in an environment enriched with 5% CO2.
To detect the expression levels of feature genes predicted by bioinformatics in both normal human chondrocyte cells and OA human chondrocyte cells. Intracellular RNA was isolated from these cells for analysis. Gene amplification was performed according to the manufacturer’s instructions (GeniuScript III Select RT Kit for qPCR, uGreener Flex qPCR 2X Mix, U&G Bio). Primer sequences used for qRT-PCR were as follows: PGD forward 5’- CCTATGAACTCTTGGCCAAAC-3’, reverse 5’-AAACTGTGTGTCTTCCTCAGG-3’; SLC7A5 forward 5’-GGGAAGGGTGATGTGTCC AAT-3’, reverse 5’-AGGCCGCTGTATAATGCCAG-3’; TKT forward 5’-CCCGAAACAAGC TTTCACCG-3’, reverse 5’-TAGACTCGGTAGCTGGCCTT-3’; GAPDH forward 5’-ACATCGCTCAGACACCATG-3’, reverse 5’-TGTAGTTGA GGTCAATGAAGGG-3’. Gene expression levels were quantified as relative fold changes, normalized against GAPDH expression.
RESULTS
Identification of DETMRGs and their functional enrichment
We used GSE51588 microarray datasets, which contain totally 10 normal samples and 40 OA samples, to identify DEGs (Fig. 1A, B). DETMRGs were identified as the overlapping genes between DEGs and telomere related genes (TMRGs), resulting in a total of 34 DETMRGs (Fig. 1C). The PPI network, comprising 25 nodes and 41 edges (Fig. 1D), demonstrated that DETMRGs interact closely at the protein level.
The correlations between 34 DETMRGs are shown in the correlation matrix (Fig. 2A), glycine decarboxylase (GLDC) is highly positively correlated with RHBDL2, DDX21 with LDHA, HIST1H2AL with HIST1H3H, GATA1 with KLF1 (r > 0.85, p < 0.001); HOXA3 is significantly negatively correlated with SLC7A5 (r < –0.7, p < 0.001). Furthermore, ceRNA network was constructed with 612 interactions, including 186 miRNAs, 37 lncRNAs and 16 mRNAs (Fig. 2B). To investigate the biological functions and pathways associated with DETMRGs, we performed GO and KEGG enrichment analysis (Fig. 2C, D). Central carbon metabolism in cancer, carbon metabolism and Epstein-Barr Virus Infection etc. frequently presented in the KEGG enrichment results (Fig. 2C). On the other side, GO enrichment collected 53 results, including 45 biological processes, 2 cellular components, and 6 molecular functions. Notably, biological functions such as erythrocyte differentiation, erythrocyte homeostasis, myeloid cell homeostasis, myeloid cell differentiation, homeostasis of number of cells and response to molecule of bacterial etc. were mainly enriched in the GO enrichment results (Fig. 2D).
Selection of diagnostic biomarkers from DETMRGs for OA
We used three machine learning algorithms to improve the accuracy of OA diagnostic biomarkers from DETMRGs. First, we identified 12 predictive genes using LASSO regression analysis (Fig. 3A, B). Secondly, SVM-RFE selected 3 features (Fig. 3C). Next, we applied RandomForest algorithm to calculate the relationship between the error rate and the number of trees (Fig. 3D), and ranked 30 genes with a relative importance score > 0.2 based on MeanDecreaseGini (Fig. 3E). Finally, we used the top ten feature genes selected by RandomForest screening to intersect with the feature genes of LASSO and SVM-RFE, resulting in three feature genes: PGD, SLC7A5, and TKT (Fig. 3F).
Diagnostic value and expression pattern of the feature genes for OA
We next sought to verify the diagnostic value of the feature genes PGD, SLC7A5, and TKT for OA. The AUC values of PGD, SLC7A5, and TKT were 0.973, 0.985, and 0.995 (Fig. 4A) in training set GSE51588. Meanwhile, the expression of PGD, SLC7A5, and TKT was significantly lower compared to the normal control group (Fig. 4B). The AUC values for the three feature genes exceed 0.8, suggesting that these genes demonstrate a high level of accuracy in the diagnosis of OA. On the other hand, the AUC values of PGD, SLC7A5, and TKT were 0.837, 0.805 and 0.913 in the validation set GSE51588 (Fig. 4C), and their expression levels were also significantly lower compared to the normal control group (Fig. 4D). We examined the gene expression profiles of PGD, SLC7A5, and TKT within the RA datasets. In the GSE77298 dataset, levels of PGD, SLC7A5, and TKT in RA cases exceeded those observed in the normal control group. Conversely, GSE1919 showed no significant differences in PGD and SLC7A5 levels between RA and control subjects, although TKT levels were notably lower than those in the normal controls (Supplementary Fig. 1). The inconsistency across these datasets indicates that PGD, SLC7A5, and TKT do not serve as reliable characteristic genes for RA. Nonetheless, the role of PGD, SLC7A5, and TKT in identifying OA remained specific. To predict the risk of OA in patients, we constructed a diagnostic nomogram based on the three feature genes PGD, SLC7A5, and TKT illustrating the expression scores of each gene (Fig. 4E). The AUC value of nomogram was 1, which was better than PGD, SLC7A5, and TKT (Fig. 4F). To obtain a more applicable diagnostic model for OA, we constructed a DCA (Fig. 4G) and clinical calibration curves (Fig. 4H). Overall, the results indicated the diagnostic model for OA constructed from the three feature genes was highly predictive.
Immune cells infiltration and its association with feature genes in OA
The infiltration of immune cells in the immune microenvironment is crucial for the diagnosis and treatment of OA. We used CIBERSORT algorithm to evaluate the infiltration of 22 immune cell types in OA and normal samples (Fig. 5A). Then, we analyzed the differences in the 22 immune cells types between OA and control samples. The results revealed significant differences in 15 immune cells types, including B cell naïve, plasma cells, T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells regulatory (Tregs), NK cells activated, monocytes, macrophage M1, macrophage M2, dendritic cells (DC) resting, DC activated, mast cell resting, mast cell activated, and neutrophils (Fig. 5B). In order to determine the correlation between immune cells, Spearman was used to analyze the correlation of immune cells with differences between OA and normal samples (Fig. 5C). To explore the correlation between the feature genes PGD, SLC7A5, and TKT and immune cells, we used Spearman to analyze the correlation between PGD, SLC7A5, and TKT and differentially expressed immune cells. The results indicated that PGD, SLC7A5, and TKT were significantly positively correlated with mast cell resting and neutrophils (Fig. 5D). Results with |correlation| > 0.3 and p < 0.05 were considered significant difference.
Afterwards, we utilized the xCell algorithm to assess the infiltration of 22 immune cells in OA and normal samples (Supplementary Fig. 2). Subsequently, we analyzed the differences in the infiltration of 22 immune cells between OA and control samples. The findings revealed significant differences in 8 immune cells, namely astrocytes, CD4+ memory T-cells, CD8+ Tem, Class-switched memory B-cells, common lymphoid precursor (CLP), common myeloid progenitor, DC, and eosinophils (Supplementary Fig. 3). We conducted a Spearman analysis to examine the relationship between PGD, SLC7A5, and TKT, and differentially expressed immune cells. The results indicated that PGD, SLC7A5, and TKT were significantly positively correlated with CLP and eosinophils (Supplementary Fig. 4). Results with |correlation| > 0.3 and p < 0.05 were considered significant.
GSEA of feature genes and prediction of potential drug target
We utilized GSEA to identify the significant pathways for each gene (Fig. 6A-C). TKT was found to be positively enriched in five pathways: neutrophil degranulation, retinoblastoma gene in cancer, antimicrobial peptides, translation, and cell cycle checkpoints (Fig. 6A). SLC7A5 also showed positive enrichment in neutrophil degranulation, cell cycle checkpoints, antimicrobial peptides, cell cycle, and retinoblastoma gene in cancer (Fig. 6B). Conversely, PGD showed positive enrichment in two pathways (cellular responses to stimuli and neutrophil degranulation) and negative enrichment in three pathways (olfactory signaling pathway, olfactory transduction, and sensory perception) (Fig. 6C). The results demonstrated a positive correlation of PGD, SLC7A5, and TKT with neutrophil degranulation, aligning with the CIBERSORT algorithm findings. This suggested that downregulating PGD, SLC7A5, and TKT may reduce neutrophil infiltration and impair neutrophil functionality. Finally, our predictions identified phenobarbital, pralmorelin, and penicillamine as potential drugs targeting PGD, while gabapentin and melphalan were identified as candidates for targeting SLC7A5 (Fig. 6D).
Validation of feature genes by RT-qPCR
We performed RT-qPCR assays to assess the expression levels of PGD, SLC7A5, and TKT in both normal human chondrocytes (CP-H107) and those derived from patients with OA (402OAK-05a). The findings indicated that, compared to CP-H107, the expression levels of PGD, SLC7A5, and TKT in 402OAK-05a were significantly downregulated (Fig. 7). This observation aligns with the results obtained from our bioinformatics analysis. PGD, SLC7A5, and TKT serve as reliable telomere-related diagnostic biomarkers for OA and possess significant research potential.
DISCUSSION
OA is a prevalent chronic joint disease that poses a significant threat to human health, particularly among the elderly population. The heterogeneous nature of OA leads to diverse manifestations among patients, tissues, and cells. Consequently, various approaches have been proposed for diagnosing and treating OA, including targeting cell pyroptosis [17], ferroptosis [18,19], and autophagy [20-23], as well as investigating their individual roles in OA and potential as molecular targets for treatment [1]. Constant research into new approaches and methods is necessary to gain a deeper understanding of the pathological mechanisms of OA, thereby facilitating more accurate diagnosis and precise treatment. Recent years have witnessed significant advancements in the screening and identification of diagnostic biomarkers for OA through bioinformatics, particularly focusing on telomere-related genes as diagnostic biomarkers, highlighting their critical importance.
The GO analysis results suggest that DETMRGs are intricately linked to the differentiation and homeostasis of red blood cells, bone marrow cells, and cell quantity. Demirci Yildirim and Sari [24] reported that during OA, there may be an increase in positive acute phase markers such as erythrocyte sedimentation rate, C-reactive protein, complement levels, mild leukocytosis, and thrombocytosis. Mariano et al. [25] observed a reduction in the expression of red blood cell antioxidant enzymes in the synovial fluid and blood of OA patients. OA is characterized by low-grade chronic inflammation, and abnormal red blood cell function can lead to the production of reactive oxygen species (ROS). Oxidative stress and inflammation are closely linked, and both contribute to joint disease. The KEGG analysis results show that carbon metabolism is closely associated with OA, and previous studies have demonstrated that several pathways of central carbon metabolism are critical for maintaining cartilage homeostasis [26]. Dysregulation energy metabolism can impair chondrocyte function in OA. OA chondrocytes have been found to exhibit increased basal extracellular acidification rate and lactate production compared to non-OA chondrocytes [27]. The previous research findings are consistent with our analysis, indicating that the occurrence of OA is closely linked to the differentiation and homeostasis of red blood cells and bone marrow cells, as well as carbon metabolism.
PPI analysis and functional similarity results have revealed that certain genes, such as GLDC, HOXA3, LDHA, and TKT, play prominent roles in the network. Li et al. [28] discovered that homeobox A3 (HOXA3) is a target gene for miR-10a-5p, and that miR-10a-5p promotes the progression of OA by inhibiting the expression of HOXA3. These research findings unveiled the regulatory mechanism of miR-10a-5p, offering a potential new therapeutic target for OA. The analysis results of the ecRNA network have shown a close relationship between HOXA3 and hsa-miR-4295, hsa-miR-301a-3p, hsa-miR-3666, hsa-miR-301b-3p, hsa-miR-454-3p, hsa-miR-130b-3p, hsa-miR-130a-3p, hsa-miR-10b-5p, as well as hsa-miR-10a-5p, which align with previous results. This suggests that other miRNAs related to HOXA3 may also hold potential for OA treatment, warranting further investigation in the future. Arra et al. [29] found that under inflammatory conditions, chondrocytes undergo metabolic changes regulated by NF-κB activation, leading to reprogramming of cellular metabolism towards glycolysis and lactate dehydrogenase A (LDHA). They identified LDHA as a potential therapeutic target for OA treatment. Gao et al.’s study [30] demonstrated that Hsa-piR-019914 reduces LDHA-dependent ROS production, decreases inflammation-mediated chondrocyte apoptosis, maintains cell proliferation and colony formation by targeting LDHA expression, and has a protective effect on chondrocytes in vitro. Furthermore, the analysis results of the ceRNA network indicate that LDHA is a target gene of hsa-miR-33a-5p, hsa-miR-33b-5p, hsa-miR-34a-5p, hsa-miR-449a, hsa-miR-449b-5p and hsa-miR-34c-5p suggesting that regulating these miRNAs could offer a novel therapeutic approach for targeting LDHA in OA treatment. Finally, current literature lacks reports of GLDC and thiamine-dependent enzyme (TKT) as targets for OA treatment, highlighting the need for further research in these areas.
To ensure the reliability of our results, we integrated the outputs of three machine learning algorithms (SVM-RFE, LASSO, RandomForest) and ultimately identified three feature genes: PGD, SLC7A5, and TKT. PGD, a phosphogluconate dehydrogenase, is the second dehydrogenase in the pentose phosphate pathway [31]. While Wang et al. [32] suggested PGD as a potential biomarker for OA diagnosis, research on its therapeutic applications in OA remains limited. However, the study conducted by Batsios et al. [33] indicates that the expression of telomerase reverse transcriptase (TERT) in oligodendroglioma correlates with the upregulation of PGD. TERT plays a vital role in preserving telomeres, which is essential for the immortality of various cancers. Consequently, PGD downregulation may affect chondrocyte aging and apoptosis in OA patients through modulation of TERT expression. SLC7A5 is a cytoplasmic membrane component and an amino acid transporter protein involved in the decarboxylation reactions of various amino acids and their transportation into the cytoplasm [34,35]. The research conducted by He et al. [36] demonstrated that the downregulation of the amino acid transporter SLC7A5 can result in compromised cellular amino acid absorption, which subsequently inhibits mTOR signaling and hampers the proliferation of BeWo cells, while also promoting autophagy and apoptosis. Consequently, the downregulation of SLC7A5 may similarly facilitate autophagy and apoptosis in chondrocytes of OA patients by influencing the mTOR signaling pathway. Moreover, solute carrier transporters play a crucial role in maintaining joint tissue homeostasis and regulating inflammatory processes in both RA and OA [37]. These findings support our prediction that SLC7A5 could be a promising target for OA diagnosis and treatment. TKT, a thiamine-dependent enzyme, directs excess sugar phosphates into glycolysis through the pentose phosphate pathway [38]. Currently, there are no research reports connecting TKT with OA. In conclusion, the research on PGD, SLC7A5, and TKT in the context of OA is limited, indicating a significant gap that warrants further investigation in the future.
The results of immune cell infiltration related to feature genes indicate a reduction in mast cells, neutrophils, CLP, and eosinophils infiltration in OA samples. Zhao et al.’s [39] study revealed a correlation between mast cells presence in the synovium of OA and disease severity. Additionally, they found that pharmacological blockade of histamine, a key mast cell mediator, has the potential to improve OA disease outcomes. Similarly, studies have shown that mast cells can exhibit anti-inflammatory effects by secreting mediators that inactivate pro-inflammatory cytokines, including interleukin-6 [40]. Furthermore, research has shown that neutrophils play a role in OA through multiple mechanisms, including tissue degeneration, osteophyte development, and the release of inflammatory cytokines and chemokines via neutrophil elastase [41]. This aligns with the findings from our CIBERSORT algorithm analysis. Neutrophils significantly contribute to early host defense mechanisms due to their ability to migrate toward inflamed tissues [42]. The cytoskeleton is essential for neutrophils function; thus, SLC7A5, as a cytoskeletal component, is closely linked to these cells. A decrease in neutrophils could explain the downregulation of SLC7A5 observed in OA. Neutrophil energy primarily relies on glycolysis [43], with both PGD and TKT serving as key enzymes in this metabolic pathway, also participating in the pentose phosphate pathway. Consequently, PGD and TKT might influence neutrophil proliferation by modulating their glycolytic process. Discovery the role of neutrophils in OA has revealed potential new therapeutic targets and diagnostic methods. However, no reports currently exist on the study of CLP and eosinophils in OA, highlighting the need for further research to establish their relationship with the condition. These results suggest that PGD, SLC7A5, and TKT may promote OA pathogenesis by regulating the activity of these cells, a promising area for future investigation. Mast cells, neutrophils, CLP, and eosinophils may indeed become key cells for future OA treatment targets.
The GSEA results indicated that PGD, SLC7A5, and TKT were positively correlated with neutrophil degranulation. These results align with those from our CIBERSORT analysis. Furthermore, immune cell infiltration results revealed a positive correlation between the levels of PGD, SLC7A5, and TKT and neutrophil infiltration levels. The findings indicated that OA patients, downregulation of PGD, SLC7A5, and TKT may correlate with reduced neutrophil levels. Given that PGD, SLC7A5, and TKT are associated with telomere function, their decreased expression reflects neutrophil senescence and apoptosis. This insight suggests that enhancing the expression of PGD, SLC7A5, and TKT in neutrophils could promote their proliferation and maintenance, potentially offering novel therapeutic strategies for alleviating and managing OA.
Potential drugs for PGD, SLC7A5, and TKT were predicted, with no potential drugs identified for TKT, and phenobarbital, pralmorelin, and penicillamine identified for treating OA in PGD. However, phenobarbital and pralmorelin have been rarely reported for OA treatment. On the other hand, penicillamine, as a slow-release drug, has shown promise in reducing imaging progression and improving clinical and biochemical indicators in RA [44]. It also lowers urinary hydroxyl pyridine cross-linking levels and reduces the index of bone involvement in arthritis diseases [45]. Additionally, gabapentin and melphalan have been predicted as potential drugs for treating OA associated with SLC7A5, with gabapentin being effective in alleviating pain in knee OA [46]. Melphalan is an effective DNA alkylator with anti-tumor activity, although no reports currently exist regarding its use in OA treatment. Overall, research on drugs related to PGD, SLC7A5, and TKT in OA is limited, and further studies and evaluation of their effects could provide new directions for the clinical management of OA in the future.
In conclusion, although the integrated bioinformatics analysis has allowed us to pinpoint telomere-related genes PGD, SLC7A5, and TKT as potential diagnostic biomarkers for OA, revealed their connection to immune cell infiltration in OA, and predicted potential therapeutic drugs, further in vivo and in vitro experiments are necessary to validate our findings. In summary, these findings may contribute to our understanding of the role of telomeres in the onset of OA and provide valuable insights for its diagnosis and treatment.
SUPPLEMENTARY MATERIALS
Supplementary data including four figures can be found with this article online at https://doi.org/10.4196/kjpp.24.322
XML Download