Identification of osteoarthritis-characteristic genes and immunological micro-environment features through bioinformatics and machine learning-based approaches
BMC Medical Genomics volume 16, Article number: 236 (2023)
Osteoarthritis (OA) is a multifaceted chronic joint disease characterized by complex mechanisms. It has a detrimental impact on the quality of life for individuals in the middle-aged and elderly population while also imposing a significant socioeconomic burden. At present, there remains a lack of comprehensive understanding regarding the pathophysiology of OA. The objective of this study was to examine the genes, functional pathways, and immune infiltration characteristics associated with the development and advancement of OA.
The Gene Expression Omnibus (GEO) database was utilized to acquire gene expression profiles. The R software was employed to conduct the screening of differentially expressed genes (DEGs) and perform enrichment analysis on these genes. The OA-characteristic genes were identified using the Weighted Gene Co-expression Network Analysis (WGCNA) and the Lasso algorithm. In addition, the infiltration levels of immune cells in cartilage were assessed using single-sample gene set enrichment analysis (ssGSEA). Subsequently, a correlation analysis was conducted to examine the relationship between immune cells and the OA-characteristic genes.
A total of 80 DEGs were identified. As determined by functional enrichment, these DEGs were associated with chondrocyte metabolism, apoptosis, and inflammation. Three OA-characteristic genes were identified using WGCNA and the lasso algorithm, and their expression levels were then validated using the verification set. Finally, the analysis of immune cell infiltration revealed that T cells and B cells were primarily associated with OA. In addition, Tspan2, HtrA1 demonstrated a correlation with some of the infiltrating immune cells.
The findings of an extensive bioinformatics analysis revealed that OA is correlated with a variety of distinct genes, functional pathways, and processes involving immune cell infiltration. The present study has successfully identified characteristic genes and functional pathways that hold potential as biomarkers for guiding drug treatment and facilitating molecular-level research on OA.
Osteoarthritis (OA) is a highly prevalent musculoskeletal condition on a global scale, and it stands as the primary contributor to disability . OA predominantly impacts the knee joint, and its occurrence is associated with various factors, such as joint injury, joint dysplasia, advancing age, obesity, and genetic predisposition . The quality of life of middle-aged and older individuals can be adversely affected by chronic pain and functional impairment, which can also impose a substantial economic burden on society . OA impacts various anatomical structures within the periarticular region. The process of irreversible articular cartilage destruction may initiate prior to the manifestation of clinical symptoms and radiographic abnormalities . Despite numerous investigations, the pathophysiology and therapeutic options for OA remain elusive. Therefore, the development of innovative biomarkers capable of detecting the degradation of articular cartilage is imperative for the diagnosis and management of OA.
Utilizing high-throughput gene expression profiling, the Weighted Gene Co-expression Network Analysis (WGCNA) method can be employed to discern gene co-expression networks within diseases . This approach facilitates the identification of gene modules that exhibit significant connectivity as well as, the identification of core genes within these modules .
The flowchart for this study is depicted in Fig. 1. First, gene expression profiles of OA cartilage tissues were acquired from the Gene Expression Omnibus (GEO) database. Subsequently, the GSE114007 and GSE57218 datasets were selected and retrieved for further analysis. To identify the differentially expressed genes (DEGs) in healthy and OA-affected cartilage, the combined and corrected datasets were subjected to differential analysis using the limma package in the R software. Subsequently, functional enrichment analyses and the construction of protein-protein interaction networks were performed for the DEGs. The combined dataset was then subjected to WGCNA analysis in order to identify the gene module that is most pertinent to the disease, as well as the core genes within the key module. The core genes derived from WGCNA were compared with DEGs to identify genes that were present in both datasets. Subsequently, the lasso algorithm was employed to acquire the OA-characteristic genes, which were subsequently validated using the GSE169077 dataset. The validation dataset GSE169077 was obtained from the GEO database as well. The present study also encompassed an examination of immune cell infiltration within the samples in order to ascertain the differential infiltration of immune cells in healthy and OA-affected cartilage. Furthermore, the relationship between the OA-characteristic genes and the immune cells was determined via correlation analysis.
This study employed bioinformatics and machine learning techniques to identify the signaling pathways and OA-characteristic genes. In addition, we investigated the relationship between these characteristic genes and immune cells. Furthermore, the identified OA-characteristic genes were validated using additional datasets. Thus, the current study establishes a foundation for understanding the pathogenesis and progression of OA, offering valuable insights for the development of gene-targeted therapies.
Materials & methods
Downloading and merging of data sets
The GEO database https://www.ncbi.nlm.nih.gov/GEO/), is a globally recognized public repository that offers free access to high-throughput datasets. The OA cartilage-related gene expression datasets GSE114007, GSE57218, and GSE169077 (Table 1) were obtained through a comprehensive search and subsequent download. It is important to note that GSE114007 and GSE57218 were utilized as training group datasets, while GSE169077 served as the validation group dataset for the present study. The expression matrix of the training group datasets was consolidated, and batch effects were mitigated using the SVA package in R version 4.21 software. As the study was conducted utilizing a publicly available database and did not involve the use of animal or human subjects, it was not necessary to obtain approval from an institutional review board.
Identification of DEGs
We analyzed gene expression between normal and OA articular cartilage using the limma package in the R programming language based on Bayesian calculations of t-value, f-value, and logarithmic dominance. Correspondingly, we identified the DEGs that between the two groups that met the specified conditions (|logFC| ≥1 and adj. P. Val < 0.05). The ggplot2 and pheatmap packages are then utilized to generate a volcanic map and heat map, respectively, in order to visually represent and elucidate the obtained outcomes.
The functionally enrichment and PPI network construction of DEGs
The Clusterprofiler package in the R programming language was utilized to perform Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) [7,8,9], and Gene Set Enrichment Analysis (GSEA) enrichment analyses. In addition, the PPI network was constructed using the String database (https://cn.string-db.org), with a minimum connectivity score threshold of 0.4.
Weighted Gene Co-expression Network Analysis (WGCNA)
The WGCNA is a bioinformatics approach utilized for the identification of gene modules associated with diseases. This approach involves the construction of gene co-expression networks and the subsequent identification of significant pathogenic pathways or potential therapeutic targets . The genes that exhibited significant correlations were organized into modules through the application of clustering and dynamic cropping techniques. In order to ascertain gene modules that are linked to OA, we applied a multiplication factor of 60 to the quantity of RNAs within the module set and established a threshold setting of 0.25 as the cut height. The significance of each module was visually represented. The key module, which is most closely linked to OA, was selected for analysis. Core genes within this module were identified based on criteria including a gene significance value > 0.5 and a gene module correlation > 0.80.
Screening of OA characteristic genes
The identification of overlapping genes was performed by utilizing the Venn package of R version 4.21 software, which involved the integration of DEGs and module-key genes. The LASSO regression model can effectively compress the variable coefficients in a regression model by employing a penalty function. The process of identifying the minimum classification error value (λ) has the potential to lead to the reduction of variable dimensionality. The identification of OA-characteristic genes was performed by screening the acquired overlapping genes, followed by the construction of a lasso regression model using the R programming language.
Validation of OA characteristic genes in training and testing groups
A comparative analysis was performed to assess the disparity in gene expression in healthy and OA-affected cartilage. The statistical analysis was performed utilizing the limma package in the R programming language. A significance level of 0.05 was adopted for determining statistical significance. The “proc” package of R version 4.21 software was employed to generate receiver operating characteristic (ROC) curves and calculate the area under the curve (AUC) values for healthy and OA-affected samples to investigate the classification effect of the disease-characteristic genes.
Immune cell infiltration analysis
Single sample gene set enrichment analysis (ssGSEA) was used to evaluate the scores of 28 immune cells to obtain a score file representing the immune cell composition in each individual sample. The immune cell scores of individuals in OA and healthy groups were subjected to statistical analysis using t-tests. The resulting data was then visually represented through the use of violin plots. The association between the disease-characteristic genes and immune cells was examined using the correlation test function in the R software. The results were then visualized using the ggplot package. Subsequently, a correlation analysis was conducted to examine the relationship between each dysregulated immune cell, and the ensuing outcomes were graphically represented.
Data preprocessing and differential expression gene screening
Following the integration of the datasets GSE57218 and GSE114007, a comprehensive analysis of gene expression profiling was conducted. Based on the principles of differential analysis, a total of 80 DEGs were identified, with 65 genes exhibiting upregulation and 15 genes displaying downregulation (Fig. 2A). Subsequently, a heat map representing the DEGs was generated (Fig. 2B).
Enrichment analysis of differentially expressed genes
The DEGs were subjected to analysis using the clusterprofiler package in the R programming language in order to obtain enrichment results for the GO and KEGG pathways. A comprehensive analysis revealed the identification of 93 biological processes (BP), 14 cellular components (CC), 25 molecular functions (MF), and 10 signaling pathways, all of which were found to be statistically significant at a significance level of P < 0.05. The outcomes of BP are primarily associated with the organization of the extracellular matrix (ECM), the response to nutrients, the response to nutrient availability, and the organization of extracellular structures, etc. The top 10 enrichment outcomes for MF and CC are shown in Fig. 3A. According to the KEGG results, DEGs are predominantly involved in the ECM-receptor interaction, PI3K-Akt, protein digestion and absorption, focal adhesion, apelin, relaxin, p53, and other signaling pathways (Fig. 3B).
GSEA analysis and ppi network construction
The results of the GSEA analysis are presented in Table 2. The gene set of the OA group exhibits significant enrichment in datasets associated with T cells and B cells (Fig. 4A, B). In addition, as shown in Fig. 4C, the PPI network reveals the connections between DEGs and pathways associated with OA.
Weighted Gene Co-expression Network Analysis (WGCNA)
The WGCNA was conducted on the merged data set, and a value of β = 7 was selected in order to achieve a suitable fit of the gene distribution to a scale-free network. Sample trees and soft thresholds were constructed based on connectivity degrees, as depicted in Fig. 5A and B, respectively. The vertical axes of the graph depict the scale-free topology fit index R^2 and mean connectivity, while the relationship between these variables and the soft threshold is illustrated. The observed negative correlation between the variables K and P (k) (correlation coefficient = 0.89) suggests that selecting β = 7 satisfies the requirements for constructing a gene scale-free network (Fig. 5 C, D).
Nine modules were identified in the gene expression tree and module eigen correlation heatmap following the conversion of the gene expression matrix into an adjacency and topological matrix (Fig. 6A, B). A score was assigned to each module to determine the importance of gene signatures (Fig. 6C). Moreover, the green module, which received the highest score, was identified as the key module. The three central genes of the principal module were examined (Fig. 6D) using the criteria of gene importance > 0.5 and gene module correlation > 0.8.
Screening of OA characteristic genes
Three overlapping genes were identified by intersecting the module core genes with DEGs using the R software (Fig. 7A). Following the application of the lasso regression algorithm, three distinct genes exhibiting substantial discriminatory capabilities between healthy and OA-affected samples were successfully identified (Fig. 7B, C).
Validation of OA characteristic genes in training and testing groups
Based on the findings obtained by our research team, it was observed that the OA samples exhibited elevated expression levels of HtrA1 and Tspan2 compared to the healthy samples within the training group (Fig. 8A, B). Conversely, a notable down-regulation of Herc5 was observed in the OA-affected samples (Fig. 8C). In order to enhance our comprehension of the diagnostic efficacy of HtrA1, Tspan2, and Herc5, we conducted ROC analyses on these three biomarkers. The AUC values for the three OA characteristics examined in this study are > 0.85, indicating a robust ability to differentiate OA samples from normal samples (Fig. 8D, E, and F).
The GSE169077 dataset was selected for the purpose of validating the expression of specific genes associated with OA, viz. Herc5, Tspan2, and HtrA1 in both healthy and OA-affected cartilage. This validation was conducted through the analysis of differential gene expression (Fig. 9A, B, and C) and through ROC analysis (Fig. 9D, E, and F). The research team observed a notable increase in the expression levels of HtrA1 and Tspan2, while the expression levels of Herc5 were observed to be lower in cartilage associated with OA. The findings exhibited congruity with the outcomes observed in the training cohort. Based on our research findings, it can be concluded that HtrA1, Tspan2, and Herc5 possess significant potential as biomarkers for the diagnosis of OA.
Infiltrated immune cells were correlated with Tspan2, HtrA1
The characteristics of immune cells were also investigated by employing the cryptosort method. The study yielded data on immune cell infiltration, specifically the quantification of 28 distinct immune cell types within cartilage samples. Figure 10 (A) displays the outcomes of the comparative analysis conducted between the OA and control groups. Dysregulated levels were observed in various cell types within the OA samples, including activated dendritic cells, immature dendritic cells, Th1 cells, Th2 cells, Th17 cells, T follicular helper cells, Gamma delta T cells, regulatory T cells, central memory CD8 T cells, activated B cells, memory B cells, eosinophils, myeloid-derived suppressor cells (MDSC), and CD56dim natural killer cells. Figure 10 (B, C) illustrate the correlation between immune cells and OA-characteristic genes (Herc5, Tspan2, and HtrA1). Additionally, it depicts the interrelationship among these dysregulated immune cells. The relationship between Tspan2 and dysregulated immune cells involves the activation of B cells, memory B cells, γδT cells, and Th2 cells. Similarly, dysregulated immune cells associated with HtrA1 involve activated B cells, memory B cells, eosinophils, Th2 cells, and Th17 cells. The findings of our study, therefore, indicate that Tspan2 and HtrA1 exhibit interactions with various immune cells, thereby contributing to the progression of OA.
OA is a chronic degenerative joint disorder that manifests as joint pain, inflammation, rigidity, and progressive loss of mobility, severely affecting the quality of life for the middle-aged and elderly populace. In addition to synovitis, the primary characteristics of this condition encompass the degeneration and destruction of articular chondrocytes [1,2,3,4]. One of the primary roles of cartilage is to effectively absorb mechanical forces by means of its ECM . The ECM is a highly dynamic, complex molecular network primarily consisting of hyaluronic acid, collagen fibers, and proteoglycans [11, 12]. These constituents play a crucial role in preserving the normal structure and function of cartilage [11, 12]. The joint undergoes a series of molecular pathway dysfunctions during the progression of OA, leading to an imbalance in proteolysis. This imbalance ultimately leads to the degradation of cartilage’s structural integrity and biomechanical properties . Furthermore, these alterations take place prior to the degradation of cartilage .
Herein, enrichment analyses using GO and KEGG were performed on a set of 80 DEGs derived from samples of individuals with OA and healthy controls. Based on the results of the GO enrichment analysis, it is evident that DEGs exhibit a prominent presence in various biological processes, including extracellular matrix organization, extracellular structure organization, response to transforming growth factor beta, external encapsulating structure organization, response to nutrient levels, and other related processes. Furthermore, the findings of this study indicate that the DEGs were predominantly enriched in two cellular compartments: the extracellular matrix, which contains collagen, and the lumen of the endoplasmic reticulum. The MF indicated that the DEGs were implicated in various biological processes, including extracellular matrix structural constituents, growth factor binding, receptor-ligand activity, signaling receptor activator activity, and protease binding. In addition, the findings from the GO enrichment analysis indicated that the OA cartilage exhibited more pronounced alterations in the extracellular matrix and a greater degree of involvement in diverse biological processes compared to healthy cartilage. Moreover, the KEGG enrichment analysis indicates that the focal adhesion and PI3K-Akt signaling pathways exhibit the highest count of enriched genes. Adhesion is one of the major ECM receptor-interacting signaling pathways in which integrins play a crucial role. When ECM binds to integrins, focal adhesion kinase (FAK) in cells is phosphorylated. Following its activation, FAK participates in a diverse range of cellular processes, encompassing cell migration, cell differentiation, matrix remodeling, growth factor signaling, and cell survival [15, 16]. Thus, the integrin-FAK interaction was confirmed to be essential in preserving chondrocyte homeostasis . The PI3K-Akt signaling pathway is a very complex pathway comprising over 150 proteins, which plays a crucial role in various cellular processes such as inflammation, metabolism, cell cycle regulation, cell survival, and programmed cell death . These processes are essential for cellular homeostasis [19,20,21]. The PI3K-Akt signaling pathway exerts a suppressive influence on ECM catabolism. The activation of PI3K-Akt by TGF-b leads to a decrease in the expression of MMP13 . A separate study also demonstrated that IGF-1 activates PI3K and extracellular signal-regulated kinases (ERK), resulting in an increased expression of COL2A1 and the inhibition of MMP13 . Moreover, previous studies have confirmed the significance of the PI3K-Akt signaling pathway in the regulation of chondrocyte survival and apoptosis. Under various pathological conditions, the activation of signals has been observed to have the ability to hinder chondrocyte apoptosis, thereby restricting the advancement of OA . Additional pathways implicated in the pathogenesis of OA have been elucidated through pertinent investigations. These pathways encompass protein digestion and absorption, ECM-receptor interaction, Apelin, and the p53 signaling pathway. The aforementioned pathways indicate a strong correlation between the metabolism, inflammatory response, and apoptosis process of chondrocytes and OA [25,26,27,28]. It is imperative to acknowledge that certain pathways, such as the relaxin pathway, hold significant potential for research. Several studies have indicated that relaxin receptors are present in the ligament, cartilage, and synovium of the trapeziometacarpal joint (TM), thereby exerting an influence on the stability of the joint. The potential reason for this phenomenon may be attributed to its involvement in the regulation of collagen, as suggested by previous research . Furthermore, the GSEA revealed a notable enrichment of OA cells in datasets associated with T-cells and B-cells, surpassing the enrichment observed in normal chondrocytes. This suggests a significant impact of inflammation on the progression of OA. Moreover, the findings of extensive enrichment analyses suggest that various pathways, including inflammation, metabolism, apoptosis, and others, may play a role in the processes associated with OA. Furthermore, molecules that have the ability to participate in multiple pathways concurrently may hold promise as potential targets for therapeutic interventions.
Three genes, viz. HtrA1, Tspan2, and Herc5, were identified as having the strongest association with OA through the utilization of two machine learning algorithms. Herc5, a protein with multiple domains and a molecular weight of 114 kDa, belongs to the HEC E3 ubiquitin ligase and RCC1 superfamily. It serves as the primary E3 ligase responsible for the modification of ISG15 [30, 31]. According to reports, there is evidence suggesting that ISG15 has the ability to interact with p53. Additionally, the down-regulation of Herc5 expression has been observed to indirectly facilitate the degradation of p53, resulting in an increase in apoptosis. This process has been well-studied in cancer . Recent research has suggested that Herc5 may also be influenced by IL-1β and TNF-α through distinct signal transduction pathways . Furthermore, the downregulation of Herc5 triggers the activation of the IL-17 A pathway, a crucial component in the inflammatory response, leading to an upregulation of CXCL13, CXCL15, and CXCL16 expression . According to our investigation, it was observed that the expression of Herc5 was reduced in chondrocytes affected by OA. This particular modification could potentially be linked to the occurrence of apoptosis and the inflammatory reaction in individuals with OA. While the immunological importance of Herc5 has been extensively validated, its precise involvement in the pathogenesis of OA remains unclear and necessitates additional research. Tetraspanin 2 (Tspan2), which belongs to the tetraspanin superfamily , has been identified on the plasma membrane as well as intracellular organelle membranes in a wide range of cell and tissue types. This regulatory mechanism plays a crucial role in the regulation of cancer, the immune system, and infectious diseases [36, 37]. Tspan2 is responsible for the regulation of protein transportation, specifically its partner proteins such as CXCL-12 and CXCR-4. Additionally, Tspan2 is linked to the expression of IL-13 and IL-10 . Tspan2 has the ability to interact with PIK3-R3 and contribute to the PIK3 signaling pathway . Furthermore, Tspan2, being a cytokine with proapoptotic properties, has the ability to induce cell apoptosis through the activation of the JNK, Wnt, and Akt signaling pathways . The results of our analysis indicate that there was an increase in tspan2 expression in OA cartilage when compared to samples from healthy individuals. Therefore, it is postulated that Tspan2 potentially participates in various chemokine, metabolic, and apoptotic pathways in the progression of OA; however, additional research is required. HtrA1 serves as the primary protease in human OA chondrocytes, effectively initiating catabolic pathways that lead to the deterioration of cartilage integrity [39, 40]. This degradation is facilitated through the cleavage of the fibronectin fragment within the ECM . Furthermore, aside from its intrinsic metabolic function, the indirect activation of HtrA1 and subsequent protein cleavage can be facilitated by the activation of other proteases or the inhibition of HtrA1 protease inhibitors . Moreover, previous studies have demonstrated that the expression of HtrA1 is increased during the initial stages of cartilage degeneration or in cases of minor damage. This finding implies that the expression of HtrA1 takes place during the initial phases of OA . The findings of our study, thus, provide evidence that the upregulation of HtrA1 in chondrocytes affected by OA substantiates its involvement in the development of the disease.
Multiple studies have provided evidence to support the notion that the infiltration of immune cells plays a substantial role in both the occurrence and advancement of OA . The assessment and characterization of the variety of immune cells that infiltrate the affected area are of utmost importance in understanding the underlying molecular mechanisms of OA and developing novel targets for immunotherapy. Nevertheless, the precise functions and contributions of different immune cells within the microenvironment linked to OA remain unclear. As a result, we conducted an immune cell infiltration analysis on the samples. The findings indicate that there is dysregulation in the infiltration levels of various immune cell types, including activated dendritic cells, immature dendritic cells, Th1 cells, Th2 cells, Th17 cells, Tfh cells, γδT cells, Treg cells, central memory CD8 T cells, activated B cells, memory B cells, eosinophils, CD56dim natural killer cells, and MDSCs, within OA cartilage. Dendritic cells (DCs) consist predominantly of both immature and mature cells, functioning as antigen-presenting cells [43, 44]. According to the literature, it has been observed that immature dendritic cells have the ability to secrete both inflammatory cytokines and immunosuppressive cytokines . Conversely, mature dendritic cells have been found to secrete inflammatory cytokines as well as immunostimulatory cytokines . Cytokines, including TNF-α, IL-1β, IFN, and IL-6, which are secreted by DCs, play a substantial role in the development of OA . Furthermore, several studies have provided evidence that DCs can enhance the levels of MMP-1 in synovial fluid . B lymphocytes originate from pluripotent stem cells found in the bone marrow and are closely associated with OA. Patients diagnosed with OA commonly exhibit a notable presence of B lymphocytes within the synovium, and there is often a positive association between the extent of this infiltration and the level of synovial inflammation . However, scholarly investigations have revealed that modifications occurring in the matrix at the chondro-osseous junction lead to a reduction in the generation of B lymphocytes and the manifestation of B lymphocyte cytokines . Furthermore, the regulation of B lymphocyte levels is influenced by osteoclasts, and a reduction in B lymphocytes is observed as a symptom of osteosclerosis . Herein, we observed that there was a decrease in the quantity of B lymphocytes within the cartilage affected by OA, potentially attributable to changes in the composition of the OA cartilage matrix. Given the intricate and multifaceted nature of OA, it is imperative to conduct additional research to explore the involvement of B lymphocytes in OA tissues. T cells, which constitute the predominant population of synovial infiltration in patients with OA, also exert a pivotal influence on the pathogenesis of this condition. The observation of CD4 + and CD8 + T cells within synovial aggregates is frequently reported . In this study, the predominant immune cell populations found within the infiltrated cartilage of OA are regulatory T cells (Treg) and CD4 + T cell-derived subsets, including Th1, Th2, Th17, and helper follicular cells. T regulatory (Treg) cells are derived from undifferentiated T cells and possess the ability to modulate the release of anti-inflammatory cytokines, such as interleukin-10 (IL-10), as well as cytokine receptors . The occurrence of inflammation is frequently accompanied by an elevation in the Th17/Treg ratio . Nevertheless, several studies have indicated that although the secretion of IL-10 by Tregs in the peripheral blood of patients with OA is low, there is an increase in the quantity of Treg cells, potentially linked to a diminished Treg cell reaction . The present study revealed an observed augmentation in the quantity of Treg within OA cartilage. However, additional investigation is required to elucidate the precise contribution of Treg in the pathogenesis of OA. Th cells and Tfh cells originate from CD4 + T cells, and various Th cell phenotypes are associated with different cytokines. The Th1 phenotype is associated with the synthesis of proinflammatory cytokines, which play a crucial role in the immune response against intracellular pathogens (as IFN-γ and IL-1β). The Th2 phenotype plays a crucial role in facilitating humoral immunity, regulating inflammatory processes, and contributing to the wound healing response. The primary cytokines associated with the Th2 subset are interleukin-4 (IL-4) and IL-13. The Th17 phenotype has the capability to generate IL-17, which in turn triggers an autoimmune response and stimulates the activation of neutrophils . Follicular helper T cells have the capability to express a range of genes associated with inflammation and can additionally stimulate B cells to generate immunoglobulin . An increasing body of scholarly literature indicates that TFH cells could potentially influence the severity of autoimmune diseases [including rheumatoid arthritis (RA]. However, the specific mechanism by which TFH cells operate in the context of OA remains uncertain. In addition to the aforementioned infiltration, there exists a presence of atypical T cells, specifically γδT cells. Several studies have provided evidence suggesting that γδT cells play a substantial role in the development and progression of RA . While the precise mechanism by which non-traditional T cells contribute to OA remains unclear, it would be advantageous to explore the importance of these cells in the context of OA. Furthermore, it has been established that the dysregulation of CD56 natural killer cells, MDSCs, and eosinophils contributes to the onset and advancement of OA [56,57,58].
In addition, we investigated the potential association between genes characteristic of OA and dysregulated immune cells within the chondrocytes of patients diagnosed with OA. Tspan2 exhibits a correlation with γδT cells, Th2 cells, activated B cells, and memory B cells, whereas HtrA1 demonstrates an association with Th2 cells, Th17 cells, activated B cells, and memory B cells. This serves as the foundation for the involvement of OA-characteristic genes in the process of immune cell infiltration.
Nonetheless, our study does exhibit certain limitations. Firstly, given that OA exerts its influence across a spectrum of tissues, confining our assessment of immune cell infiltration solely to cartilage tissue possesses inherent constraints in portraying the complete engagement of immune cells in OA. Secondly, the correlation analysis of OA characteristic genes and infiltrating immune cells can only suggest a correlation, not reveal their causal relationship.
Our investigation has revealed that osteoarthritis (OA) entails a multifaceted process influenced by numerous signaling pathways. The characteristic OA-related genes we identified, namely Herc5, Tspan2, and HtrA1, exhibit involvement in a diverse array of biological processes encompassing inflammatory responses, metabolism, and apoptosis. Furthermore, these genes are intricately linked to immune cell functions. This study not only serves as a wellspring of insight into unraveling the molecular underpinnings accountable for cartilage damage in OA but also furnishes a conceptual and practical foundation for the early identification and targeted therapeutic interventions of osteoarthritis.
Availability of data and materials
All data and materials for this article are derived from the geo database and are genuinely available.The detailed website is as follows:
Neogi T. The epidemiology and impact of pain in osteoarthritis. Osteoarthritis Cartilage. 2013;21(9):1145–53. https://doi.org/10.1016/j.joca.2013.03.018. PMID: 23973124; PMCID: PMC3753584.
Yelin E, Weinstein S, King T. The burden of musculoskeletal diseases in the United States. Semin Arthritis Rheum. 2016;46(3):259–60. https://doi.org/10.1016/j.semarthrit.2016.07.013. Epub 2016 Jul 26. PMID: 27519477.
Driban JB, Harkey MS, Barbe MF, Ward RJ, MacKay JW, Davis JE, Lu B, Price LL, Eaton CB, Lo GH, McAlindon TE. Risk factors and the natural history of accelerated knee osteoarthritis: a narrative review. BMC Musculoskelet Disord. 2020;21(1):332. https://doi.org/10.1186/s12891-020-03367-2. PMID: 32471412; PMCID: PMC7260785.
Hunter DJ, Bierma-Zeinstra S. Osteoarthritis. Lancet. 2019;393(10182):1745–59. https://doi.org/10.1016/S0140-6736(19)30417-9. PMID: 31034380.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. https://doi.org/10.1186/1471-2105-9-559. PMID: 19114008; PMCID: PMC2631488.
Lu X, Fan Y, Li M, Chang X, Qian J. HTR2B and SLC5A3 are specific markers in age-related osteoarthritis and involved in apoptosis and inflammation of Osteoarthritis Synovial cells. Front Mol Biosci. 2021;8:691602. https://doi.org/10.3389/fmolb.2021.691602. PMID: 34222340; PMCID: PMC8241908.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. https://doi.org/10.1093/nar/28.1.27. PMID: 10592173; PMCID: PMC102409.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51. https://doi.org/10.1002/pro.3715. Epub 2019 Sep 9. PMID: 31441146; PMCID: PMC6798127.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92. https://doi.org/10.1093/nar/gkac963. PMID: 36300620; PMCID: PMC9825424.
Heinegård D, Saxne T. The role of the cartilage matrix in osteoarthritis. Nat Rev Rheumatol. 2011;7(1):50–6. https://doi.org/10.1038/nrrheum.2010.198. Epub 2010 Nov 30. PMID: 21119607.
Carballo CB, Nakagawa Y, Sekiya I, Rodeo SA. Basic Science of Articular Cartilage. Clin Sports Med. 2017;36(3):413–425. https://doi.org/10.1016/j.csm.2017.02.001. Epub 2017 Apr 26. PMID: 28577703.
Guilak F, Nims RJ, Dicks A, Wu CL, Meulenbelt I. Osteoarthritis as a disease of the cartilage pericellular matrix. Matrix Biol. 2018;71–72:40–50. https://doi.org/10.1016/j.matbio.2018.05.008. Epub 2018 May 22. PMID: 29800616; PMCID: PMC6146061.
Lotz M, Loeser RF. Effects of aging on articular cartilage homeostasis. Bone. 2012;51(2):241–8. https://doi.org/10.1016/j.bone.2012.03.023. Epub 2012 Mar 28. PMID: 22487298; PMCID: PMC3372644.
Onuora S, Osteoarthritis. Cartilage matrix stiffness regulates chondrocyte metabolism and OA pathogenesis. Nat Rev Rheumatol. 2015;11(9):504. https://doi.org/10.1038/nrrheum.2015.107. Epub 2015 Aug 4. PMID: 26241187.
Gao Y, Liu S, Huang J, Guo W, Chen J, Zhang L, Zhao B, Peng J, Wang A, Wang Y, Xu W, Lu S, Yuan M, Guo Q. The ECM-cell interaction of cartilage extracellular matrix on chondrocytes. Biomed Res Int. 2014;2014:648459. https://doi.org/10.1155/2014/648459. Epub 2014 May 18. PMID: 24959581; PMCID: PMC4052144.
Loeser RF. Integrins and chondrocyte-matrix interactions in articular cartilage. Matrix Biol. 2014;39:11–6. https://doi.org/10.1016/j.matbio.2014.08.007. Epub 2014 Aug 25. PMID: 25169886; PMCID: PMC4699681.
Cheng K, Xia P, Lin Q, Shen S, Gao M, Ren S, Li X. Effects of low-intensity pulsed ultrasound on integrin-FAK-PI3K/Akt mechanochemical transduction in rabbit osteoarthritis chondrocytes. Ultrasound Med Biol. 2014;40(7):1609–18. Epub 2014 Apr 16. PMID: 24742749.
Jafari M, Ghadami E, Dadkhah T, Akhavan-Niaki H. PI3k/AKT signaling pathway: erythropoiesis and beyond. J Cell Physiol. 2019;234(3):2373–85. https://doi.org/10.1002/jcp.27262. Epub 2018 Sep 7. PMID: 30192008.
Tang F, Wang Y, Hemmings BA, Rüegg C, Xue G. PKB/Akt-dependent regulation of inflammation in cancer. Semin Cancer Biol. 2018;48:62–9. https://doi.org/10.1016/j.semcancer.2017.04.018. Epub 2017 May 2. PMID: 28476657.
Sun K, Luo J, Guo J, Yao X, Jing X, Guo F. The PI3K/AKT/mTOR signaling pathway in osteoarthritis: a narrative review. Osteoarthritis Cartilage. 2020;28(4):400–9. https://doi.org/10.1016/j.joca.2020.02.027. Epub 2020 Feb 18. PMID: 32081707.
Yao X, Zhang J, Jing X, Ye Y, Guo J, Sun K, Guo F. Fibroblast growth factor 18 exerts anti-osteoarthritic effects through PI3K-AKT signaling and mitochondrial fusion and fission. Pharmacol Res. 2019;139:314–24. Epub 2018 Sep 28. PMID: 30273654.
Zhang Q, Lai S, Hou X, Cao W, Zhang Y, Zhang Z. Protective effects of PI3K/Akt signal pathway induced cell autophagy in rat knee joint cartilage injury. Am J Transl Res. 2018;10(3):762–70. PMID: 29636866; PMCID: PMC5883117.
Zhang M, Zhou Q, Liang QQ, Li CG, Holz JD, Tang D, Sheu TJ, Li TF, Shi Q, Wang YJ. IGF-1 regulation of type II collagen and MMP-13 expression in rat endplate chondrocytes via distinct signaling pathways. Osteoarthritis Cartilage. 2009;17(1):100–6. https://doi.org/10.1016/j.joca.2008.05.007. Epub 2008 Jul 1. PMID: 18595745.
Cravero JD, Carlson CS, Im HJ, Yammani RR, Long D, Loeser RF. Increased expression of the Akt/PKB inhibitor TRB3 in osteoarthritic chondrocytes inhibits insulin-like growth factor 1-mediated cell survival and proteoglycan synthesis. Arthritis Rheum. 2009;60(2):492–500. https://doi.org/10.1002/art.24225. PMID: 19180501; PMCID: PMC2637941.
Chang TK, Wang YH, Kuo SJ, Wang SW, Tsai CH, Fong YC, Wu NL, Liu SC, Tang CH. Apelin enhances IL-1β expression in human synovial fibroblasts by inhibiting mir-144-3p through the PI3K and ERK pathways. Aging. 2020;12(10):9224–39. Epub 2020 May 18. PMID: 32420902; PMCID: PMC7288923.
Wang C, Li N, Liu Q, Su L, Wang S, Chen Y, Liu M, Lin H. The role of circRNA derived from RUNX2 in the serum of osteoarthritis and its clinical value. J Clin Lab Anal. 2021;35(7):e23858. https://doi.org/10.1002/jcla.23858. Epub 2021 Jun 24. PMID: 34165827; PMCID: PMC8274987.
Xu M, Feng M, Peng H, Qian Z, Zhao L, Wu S. Epigenetic regulation of chondrocyte hypertrophy and apoptosis through Sirt1/P53/P21 pathway in surgery-induced osteoarthritis. Biochem Biophys Res Commun. 2020;528(1):179–85. Epub 2020 Jun 1. PMID: 32499111.
Gu HY, Yang M, Guo J, Zhang C, Lin LL, Liu Y, Wei RX. Identification of the biomarkers and pathological process of Osteoarthritis: weighted gene co-expression network analysis. Front Physiol. 2019;10:275. https://doi.org/10.3389/fphys.2019.00275. PMID: 30941059; PMCID: PMC6433881.
Clifton KB, Rodner C, Wolf JM. Detection of relaxin receptor in the dorsoradial ligament, synovium, and articular cartilage of the trapeziometacarpal joint. J Orthop Res. 2014;32(8):1061-7. https://doi.org/10.1002/jor.22640. Epub 2014 May 3. PMID: 24797570.
Hadjebi O, Casas-Terradellas E, Garcia-Gonzalo FR, Rosa JL. The RCC1 superfamily: from genes, to function, to disease. Biochim Biophys Acta. 2008;1783(8):1467–79. https://doi.org/10.1016/j.bbamcr.2008.03.015. Epub 2008 Apr 10. PMID: 18442486.
Wong JJ, Pung YF, Sze NS, Chin KC. HERC5 is an IFN-induced HECT-type E3 protein ligase that mediates type I IFN-induced ISGylation of protein targets. Proc Natl Acad Sci U S A. 2006;103(28):10735–40. https://doi.org/10.1073/pnas.0600397103. Epub 2006 Jun 30. PMID: 16815975; PMCID: PMC1484417.
Wang Y, Ding Q, Xu T, Li CY, Zhou DD, Zhang L. HZ-6d targeted HERC5 to regulate p53 ISGylation in human hepatocellular carcinoma. Toxicol Appl Pharmacol. 2017;334:180–91. Epub 2017 Sep 15. PMID: 28919514.
Mathieu NA, Paparisto E, Barr SD, Spratt DE. HERC5 and the ISGylation Pathway: critical modulators of the antiviral Immune response. Viruses. 2021;13(6):1102. https://doi.org/10.3390/v13061102. PMID: 34207696; PMCID: PMC8228270.
Xue F, Higgs BW, Huang J, Morehouse C, Zhu W, Yao X, Brohawn P, Xiao Z, Sebastian Y, Liu Z, Xia Y, Shen D, Kuziora M, Dong Z, Han H, Gu Y, Gu J, Xia Q, Yao Y. HERC5 is a prognostic biomarker for post-liver transplant recurrent human hepatocellular carcinoma. J Transl Med. 2015;13:379. https://doi.org/10.1186/s12967-015-0743-2. PMID: 26653219; PMCID: PMC4676172.
Reynolds JL, Mahajan SD. Transmigration of tetraspanin 2 (Tspan2) siRNA Via Microglia Derived Exosomes across the blood brain barrier modifies the production of Immune Mediators by Microglia cells. J Neuroimmune Pharmacol. 2020;15(3):554–63. https://doi.org/10.1007/s11481-019-09895-6. Epub 2019 Dec 10. PMID: 31823250; PMCID: PMC7282939.
Grossmann A, Benlasfer N, Birth P, Hegele A, Wachsmuth F, Apelt L, Stelzl U. Phospho-tyrosine dependent protein-protein interaction network. Mol Syst Biol. 2015;11(3):794. https://doi.org/10.15252/msb.20145968. PMID: 25814554; PMCID: PMC4380928.
Zhao J, Wu W, Zhang W, Lu YW, Tou E, Ye J, Gao P, Jourd’heuil D, Singer HA, Wu M, Long X. Selective expression of TSPAN2 in vascular smooth muscle is independently regulated by TGF-β1/SMAD and myocardin/serum response factor. FASEB J. 2017;31(6):2576–91. Epub 2017 Mar 3. PMID: 28258189; PMCID: PMC5434656.
Hwang IH, Park J, Kim JM, Kim SI, Choi JS, Lee KB, Yun SH, Lee MG, Park SJ, Jang IS. Tetraspanin-2 promotes glucotoxic apoptosis by regulating the JNK/β-catenin signaling pathway in human pancreatic β cells. FASEB J. 2016;30(9):3107–16. https://doi.org/10.1096/fj.201600240RR. Epub 2016 May 31. PMID: 27247127; PMCID: PMC5001516.
Bhutada S, Li L, Willard B, Muschler G, Piuzzi N, Apte SS. Forward and reverse degradomics defines the proteolytic landscape of human knee osteoarthritic cartilage and the role of the serine protease HtrA1. Osteoarthritis Cartilage. 2022;30(8):1091–102. Epub 2022 Mar 24. PMID: 35339693.
Ding L, Guo D, Homandberg GA. The cartilage chondrolytic mechanism of fibronectin fragments involves MAP kinases: comparison of three fragments and native fibronectin. Osteoarthritis Cartilage. 2008;16(10):1253–62. https://doi.org/10.1016/j.joca.2008.02.015. Epub 2008 Apr 18. PMID: 18396067.
Hernandez PA, Wells J, Usheva E, Nakonezny PA, Barati Z, Gonzalez R, Kassem L, Henson FMD. Early-Onset Osteoarthritis originates at the chondrocyte level in hip dysplasia. Sci Rep. 2020;10(1):627. https://doi.org/10.1038/s41598-020-57431-x. PMID: 31953438; PMCID: PMC6969105.
Zhao C. Identifying the hub gene and immune infiltration of osteoarthritis by bioinformatical methods. Clin Rheumatol. 2021;40(3):1027–37. https://doi.org/10.1007/s10067-020-05311-0. Epub 2020 Aug 12. PMID: 32785809.
Reis e Sousa C. Dendritic cells in a mature age. Nat Rev Immunol. 2006;6(6):476–83. https://doi.org/10.1038/nri1845. PMID: 16691244.
Audiger C, Rahman MJ, Yun TJ, Tarbell KV, Lesage S. The importance of dendritic cells in maintaining Immune Tolerance. J Immunol. 2017;198(6):2223–31. https://doi.org/10.4049/jimmunol.1601629. PMID: 28264998; PMCID: PMC5343761.
Garg AD, Kaczmarek A, Krysko O, Vandenabeele P, Krysko DV, Agostinis P. ER stress-induced inflammation: does it aid or impede disease progression? Trends Mol Med. 2012;18(10):589–98. https://doi.org/10.1016/j.molmed.2012.06.010. Epub 2012 Aug 8. PMID: 22883813.
Kriegova E, Manukyan G, Mikulkova Z, Gabcova G, Kudelka M, Gajdos P, Gallo J. Gender-related differences observed among immune cells in synovial fluid in knee osteoarthritis. Osteoarthritis Cartilage. 2018;26(9):1247–56. https://doi.org/10.1016/j.joca.2018.04.016. Epub 2018 May 19. PMID: 29753948.
Hirohata S, Nagai T, Asako K, Tomita T, Yoshikawa H. Induction of type B synoviocyte-like cells from plasmacytoid dendritic cells of the bone marrow in rheumatoid arthritis and osteoarthritis. Clin Immunol. 2011;140(3):276–83. Epub 2011 Apr 20. PMID: 21550856.
Da RR, Qin Y, Baeten D, Zhang Y. B cell clonal expansion and somatic hypermutation of Ig variable heavy chain genes in the synovial membrane of patients with osteoarthritis. J Immunol. 2007;178(1):557 – 65. https://doi.org/10.4049/jimmunol.178.1.557. PMID: 17182596.
Sweeney E, Roberts D, Jacenko O. Altered matrix at the chondro-osseous junction leads to defects in lymphopoiesis. Ann N Y Acad Sci. 2011;1237:79–87. https://doi.org/10.1111/j.1749-6632.2011.06227.x. PMID: 22082369.
Mansour A, Anginot A, Mancini SJ, Schiff C, Carle GF, Wakkach A, Blin-Wakkach C. Osteoclast activity modulates B-cell development in the bone marrow. Cell Res. 2011;21(7):1102–15. https://doi.org/10.1038/cr.2011.21. Epub 2011 Feb 15. PMID: 21321604; PMCID: PMC3193501.
Li YS, Luo W, Zhu SA, Lei GH. T cells in Osteoarthritis: alterations and Beyond. Front Immunol. 2017;8:356. https://doi.org/10.3389/fimmu.2017.00356. PMID: 28424692; PMCID: PMC5371609.
Miyara M, Sakaguchi S. Natural regulatory T cells: mechanisms of suppression. Trends Mol Med. 2007;13(3):108–16. https://doi.org/10.1016/j.molmed.2007.01.003. Epub 2007 Jan 24. PMID: 17257897.
Paradowska-Gorycka A, Wajda A, Romanowska-Próchnicka K, Walczuk E, Kuca-Warnawin E, Kmiolek T, Stypinska B, Rzeszotarska E, Majewski D, Jagodzinski PP, Pawlik A. Th17/Treg-Related transcriptional factor expression and Cytokine Profile in patients with rheumatoid arthritis. Front Immunol. 2020;11:572858. https://doi.org/10.3389/fimmu.2020.572858. PMID: 33362761; PMCID: PMC7759671.
Li S, Wan J, Anderson W, Sun H, Zhang H, Peng X, Yu Z, Wang T, Yan X, Smith W. Downregulation of IL-10 secretion by Treg cells in osteoarthritis is associated with a reduction in Tim-3 expression. Biomed Pharmacother. 2016;79:159–65. Epub 2016 Feb 23. PMID: 27044824.
Pemmari A, Leppänen T, Hämäläinen M, Moilanen T, Moilanen E. Chondrocytes from Osteoarthritis Patients adopt distinct phenotypes in response to Central TH1/TH2/TH17 cytokines. Int J Mol Sci. 2021;22(17):9463. https://doi.org/10.3390/ijms22179463. PMID: 34502384; PMCID: PMC8431052.
Zhang L, Kirkwood CL, Sohn J, Lau A, Bayers-Thering M, Bali SK, Rachala S, Marzo JM, Anders MJ, Beier F, Kirkwood KL. Expansion of myeloid-derived suppressor cells contributes to metabolic osteoarthritis through subchondral bone remodeling. Arthritis Res Ther. 2021;23(1):287. https://doi.org/10.1186/s13075-021-02663-z. PMID: 34784965; PMCID: PMC8594239.
Qin Y, Li J, Zhou Y, Yin C, Li Y, Chen M, Du Y, Li T, Yan J. Apolipoprotein D as a potential biomarker and construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by weighted gene coexpression network analysis. Cartilage. 2021;13(1suppl):1702S–17. https://doi.org/10.1177/19476035211053824. Epub 2021 Oct 31. PMID: 34719950; PMCID: PMC8808834.
Jaime P, García-Guerrero N, Estella R, Pardo J, García-Álvarez F, Martinez-Lostao L. CD56+/CD16- natural killer cells expressing the inflammatory protease granzyme A are enriched in synovial fluid from patients with osteoarthritis. Osteoarthritis Cartilage. 2017;25(10):1708–18. Epub 2017 Jun 29. PMID: 28668542.
We thank each author for his contributions.
There is no financial support for this article.
Ethics approval and consent to participate
Since the data in this study are all derived from public databases and do not involve human or animal experiments, no ethical approval and consent to participate were required.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Da, Z., Guo, R., Sun, J. et al. Identification of osteoarthritis-characteristic genes and immunological micro-environment features through bioinformatics and machine learning-based approaches. BMC Med Genomics 16, 236 (2023). https://doi.org/10.1186/s12920-023-01672-y