Exploration of effective biomarkers and infiltrating Immune cells in Osteoarthritis based on bioinformatics analysis

Abstract

Osteoarthritis (OA) is a multi-factorial chronic joint disease mainly identified by synovial inflammation, cartilage damage, and degeneration. Our study applied bioinformatics analysis to uncover the immunity in OA and tried to explore the underlying immune-related molecular mechanism. First, OA-related gene-expression profiling data were retrieved from GEO database. Then, we analysed a series of datadata with using the xCell algorithm, GEO2R, enrichment analysis of SangerBox website, CytoHubba, ROC logistic regression and correlation analysis. Finally, Nine infiltrating immune cells with differential abundance between OA and normal samples were obtained. There were 42 IODEGs in OA, and their functions were associated with immune cells and corresponding biological processes. Moreover, 5 hub genes, including GREM1, NRP1, VEGFA, FYN and IL6R, were identified. Correlation analysis demonstrated that NRP1 was negatively associated with NKT cells, NRP1 and GREM1 were positively associated with aDC, VEGFA was positively associated with CD8+ naïve T cells, while VEGFA, FYN and IL6R were negatively associated with Macrophages M1. The 5 hub genes could be employed as effective diagnostic biomarkers for OA. In addition, they may participate in OA pathogenesis via interactions with infiltrating immune cells.

Keywords:

• Osteoarthritis
• immune infiltration
• differentially expressed genes
• bioinformatics analysis
• immune cells

Introduction

Osteoarthritis (OA) is a chronic degenerative joint disease of unknown aetiology. It is identified by articular cartilage degeneration and secondary osteoproliferation, and it clinically manifests progressive joint pain that leads to joint stiffness, limitation and even loss of function [1]. On a global scale, OA is the most common amongst chronic joint diseases with the morbidity gradually increasing with age. As of 2017, OA affected 0.303 billion people, posing a major threat for human health and creating an enormous economic burden [2]. Nonetheless, the pathogenesis of OA remains an open issue. To develop better preventive and treatment strategies for OA, it is necessary to clarify the potential mechanisms involved in occurrence and development of this medical condition.

OA is a highly complicated disease, and it is caused by an interplay between multiple factors, such as age, occupation, obesity, metabolic factors, and immune factors. Notably, immune factors have important implications throughout the course of OA [3]. Research revealed that the infiltrating macrophages, neutrophils and lymphocytes within the synovial tissue are important in regulation of the chronic inflammatory process in OA [4]. via producing a large amount of cytokine and chemokine, these infiltrating cells stimulate chondrocytes to produce matrix-degrading enzymes, accelerating cartilage breakdown and promoting the development of OA [5]. In recent years, macrophage polarization and related cytokines have been extensively investigated in OA and reported to have a close relationship with the synovial and cartilage inflammation and responses to injury [6]. In addition, the persistent articular cartilage damage and repair that are induced by chronic inflammation are correlated to multiple immune cells, including T lymphocytes, B lymphocytes, NK cells, and dendritic cells, etc. [7–9]. In this context, there is an urgent need of systemic methods that can estimate the role of immune cells in OA and explore immune-related diagnostic biomarkers to help uncover the immune mechanisms in OA.

With the rapid development and extensive application of microarray sequencing, bioinformatics has been applied in a growing number of studies for identifying novel biomarkers for various diseases, exploring the involved epigenetic alterations, and unravelling the underlying molecular mechanisms. Existing bioinformatics research demonstrated that ZNF160 could be employed as a biomarker for diagnosis of OA [10]. Infiltrating immune cells are critical for OA, but only a minority of bioinformatics studies explored their relationship. xCell algorithm is a gene signature-based method that performs cell type enrichment analysis from gene expression data for 64 immune and stromal cell types, followed by validation with extensive computer simulations and immunophenotyping [11]. Therefore, it may allow for a systemic analysis of the immune cell infiltration in OA tissues.

In the present study, we firstly applied the xCell algorithm to explore the difference in immune infiltration between OA and healthy tissues. Subsequently, a series of bioinformatics analyses were used to screen out OA-related differentially expressed genes (DEGs) and investigate potential pathological mechanism of OA, including differential expression analysis, GO and KEGG enrichment analysis, protein-protein interaction (PPI) network analysis, Hub gene screening, ROC logistic regression. Finally, Spearman method was employed to analyze the relationship between the candidate biomarkers and significant infiltrating immune cells. In all, this study may help identify new biomarkers for diagnosis of OA and provide new insights into the diagnosis and immunotherapy of OA.

Methods and materials

Geo data acquisition and processing

OA-related gene expression profiling data were retrieved from GEO database (https://www.ncbi.nlm.nih.gov/geo/) based on the following criteria: (1) Homo sapiens; (2) RNA expression profiling by microarray; (3) OA and normal tissue samples. As a result, 3 eligible datasets were retrieved, including the GSE169077 [12] used as the training set and the GSE55235 [13] and GSE55457 [13] used as the test sets. GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r) was applied for data normalization and follow-up analysis. Basic information of selected datasets acquired from GEO was shown in Table 1. Immune-related genes were retrieved from ImmPort [14] database (http://www.immport.org).

Table 1. Basic Information of selected datasets acquired from GEO.

Datasets	Platform	Species	Total	OA sample	Normal sample	Type
GSE169077	GPL96	Homo sapiens	11	6	5	Array
GSE55235	GPL96	Homo sapiens	20	10	10	Array
GSE55457	GPL96	Homo sapiens	20	10	10	Array

Estimation of immune cell infiltration

xCell [11] algorithm (https://xcell.ucsf.edu/) is a gene signature-based method that confers a relative infiltration score to 64 immune and stromal cell types according to the gene expression data measured by microarray sequencing. Here, the xCell algorithm was applied to estimate the specific infiltration of immune cells in OA and normal tissues. Principal component analysis (PCA) was performed to confirm whether the immune infiltrates between the groups were statistically significant. The infiltrating immune cells with p < .05 were selected and processed for correlation analysis using Spearman method.

Screening of immune-related DEGs

R packages “GEOquery” [15] and “limma” [16] of the interactive online tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r) were applied for data normalization and screening of DEGs (|fold change| ≥1 and adj. p < .05). The immune-related and osteoarthritis-related differentially expressed genes (IODEGs) were obtained by intersecting the DEGs with the immune-related genes using the Draw Venn Diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). The resulting DEGs and IODEGs were visualized as Heatmaps.

Functional enrichment analysis

Gene Ontology (GO) was employed to analyze the IODEGs from three aspects: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) [17]. Kyoto Encyclopaedia of Genes and Genomes (KEGG) [18] was applied to explore the IODEG-related pathways and high-level biological functions. The background gene sets selected the GO annotations from SangerBox platform and the latest KEGG pathway gene annotations that were obtained from KEGG rest API (https://www.kegg.jp/kegg/rest/keggapi.html). Enrichment analysis was fulfilled with SangerBox platform [19]. The maximum and minimum gene set sizes were set to 5,000 and 5, respectively, and the threshold of statistical significance was set to p-Value < .05 and FDR < 0.25.

Identification of Hub genes

PPI of IODEGs were analyzed with STRING (https://cn.string-db.org), and the PPIs with a score >0.4 were visualized as a network using Cytoscape (v3.9.1) [20]. The top 15 nodes (genes) were scored with 8 algorithms as implemented by plug-in CytoHubba of Cytoscape, including Maximal Clique Centrality (MCC), Degree, Edge Percolated Component (EPC), BottleNeck, EcCentricity, Closeness, Radiality, and Betweenness. Hub genes were identified by R package “UpSet” [21] of SangerBox platform and validated using the GSE55235 and GSE55457 gene sets.

Assessment of the diagnostic effectiveness of candidate biomarkers

ROC logistic regression as implemented by SPSS 20.0 was used to validate the effectiveness of the candidate biomarkers for diagnosis of OA in the GSE55235 and GSE55457 gene sets, with the area under the ROC curve (ROC-AUC) used as the criterion.

Relationship between diagnostic biomarkers and infiltrating immune cells

To gain more insight into the potential mechanism of action of infiltrating immune cells in OA, the link between the candidate diagnostic biomarkers and the infiltrating immune cells in OA was analyzed using the Spearman method (p < .05). The results were visualized with R package “ggplot2” [22] (https://sourceforge.net/projects/ggplot2.mirror/) of SangerBox platform.

Data analysis and visualization

Statistical analysis and data visualization were conducted with RStudio 3.6.3, SPSS 20.0, and Sangerbox3.0 (http://vip.sangerbox.com/home.html). Spearman method was applied for correlation analysis. All data were shown as mean ± standard deviation. A p < .05 was interpreted as statistically significant.

Results

Immune cell infiltration between OA and normal samples

Immune cell infiltration in OA and normal tissues was firstly analysed based on GEO data. PCA analysis demonstrated that the difference in immune cell infiltration between the two sample types was statistically significant (Figure 1A). The specific infiltration of 64 immune cell populations was analysed by the xCell algorithm and visualized as Heatmap (Figure 1B) and Box plot (Figure 1C). As analyzed, the infiltrating immune cells with the most significant differences between OA and normal tissues were dominated by aDC, Astrocytes, CD8+ naïve T cells, CMP, Macrophages, Macrophages M1, Neutrophils, NKT, and Th1 cells. The corresponding proportions of the nine immune cell populations were calculated (Figure 1D). Correlation analysis demonstrated that Astrocytes had the strongest positive relationship with CMP (r = 0.87), while NKT exhibited the strongest negative relationship with Neutrophils (r = −0.77) (Figure 1E). In comparison to the normal tissues, the OA tissues had significantly abundant aDC, Astrocytes, CMP, Macrophages, Macrophages M1, and Neutrophils (Figure 2A–F), but less abundant CD8+ naïve T cells, NKT, and Th1 cells (Figure 2G–I).

Figure 1. Immune cell infiltration in OA and normal tissues. (A) PCA of 64 immune cell types, and corresponding abundance calculated by xCell algorithm; (B) Heatmap of the 64 immune cell types; (C) Box plot of the 64 immune cell types; (D) Relative percentage of the 9 differential infiltrating immune cell types (The X-axis represents 5× normal samples and 6× OA samples, and the Y-axis represents the relative percentage of the 9 immune cell types in each sample; The coloured bars on the right represent the type of immune cells); (E) Heatmap of the 9 differential infiltrating immune cell types.

Figure 2. Up-regulated and down-regulated infiltrating immune cells in OA. Red represents the expression in normal tissues, and blue represents the expression in OA tissues. Up-regulated infiltrating immune cells in OA included (A) aDC; (B) CMP; (C) Macrophages M1; (D) Astrocytes; (E) Macrophages; (F) Neutrophils; while down-regulated infiltrating immune cells were (G) CD8+ naïve T cells;(H) NKT; (I) Th1 cells.

Iodegs between OA and normal tissues

In the GSE169077 dataset, a total of 382 DEGs were screened out between OA and normal tissues, including 166 up-regulated and 216 down-regulated genes in OA samples (Figure 3A, B). In the meantime, 1,793 immune-related genes were retrieved from the ImmPort database. IODEGs were obtained by intersecting the 382 DEGs with the 1,793 immune-related genes (Figure 3C–D). As a result, 42 IODEGs were acquired, including 23 up-regulated and 19 down-regulated genes (Table 2).

Figure 3. OA-related DEGs and IODEGs. (A) Volcano Plot of OA-related DEGs (VEGFA, FYN, IL6R, NRP1, and GREM1 were diagnostic biomarkers for OA); (B) Heatmap of OA-related DEGs; (C) Venn diagram of OA-related DEGs and immune-related genes retrieved from the ImmPort database; (D) Heatmap of 42 OA-related IODEGs.

Table 2. IODEGs between OA and normal groups.

IODEGs	Name	Number
Up-regulated	HLA-DRA, BMP1, TNC, GREM1, NRP1, OGN C3AR1, PDGFC, ANGPTL2, ANOS1, INHBA PENK, GDF5, SEMA3C, PTGER3, FGF1 HLA-DRB1, NGF, IL21R, PDGFRL, LTBP1, S100A10, PTN	23
Down-regulated	ADM, STC2, ANGPTL4, RARA, VEGFA FYN, APOD, CBLB, NOD1, NAMPT GBP2, IL6R, RBP4, LEP, ACKR3, NCK1 LMBR1L, HLA-E, MICA	19

IODEGs were set as |fold change|≥1 and adj. p-value <.05.

IODEGs: immune-related and OA-related differentially expressed genes.

Go and KEGG pathway analysis for IODEGs

GO terms of the IODEGs were represented by BP, CC, and MF. With respect to BP, the IODEGs were mainly enriched in immune cell-related pathways, including cell surface receptor signalling pathway, regulation of immune system process, regulation of cell adhesion, and positive regulation of immune system process (Figure 4A and B); As for CC, the IODEGs were enriched in MHC protein complex, extracellular region, and extracellular matrix (Figure 4C and D); while for MF, the IODEGs were mostly enriched in signalling receptor binding, receptor ligand activity, receptor regulator activity, T cell receptor binding, and G protein-coupled glutamate receptor binding (Figure 5A and B). The most enriched KEGG pathways of the IODEGs were involved in Th17 cell differentiation, Cytokine-cytokine receptor interaction, TGF-beta signalling pathway, PI3K-Akt signalling pathway, Antigen processing and presentation, Rheumatoid arthritis, T cell receptor signalling pathway, AND Jak-STAT signalling pathway (Figure 5C and D). Collectively, the function of IODEGs was closely linked to immune cells.

Figure 4. GO terms of the IODEGs from aspects of BP and CC. (A,B) Bubble Chart and Circle Graph of GO terms of the IODEGs from aspect of BP; (C,D) Bubble Chart and Circle Graph of GO terms of the IODEGs from aspect of CC.

Figure 5. GO terms of the IODEGs from aspect of MF and KEGG pathways. (A,B) Bubble Chart and Circle Graph of GO terms of the IODEGs from aspect of MF; (C,D) Bubble Chart and Circle Graph of KEGG pathways of the IODEGs.

Identification and visualization of Hub genes

A PPI network of the IODEGs was generated with the STRING database and visualized by the Cytoscape (Figure 6A,B). With the 8 algorithms as implemented by the plug-in CytoHubba of the Cytoscape, the top 15 nodes (genes) were scored and visualized (Figure 7A–H). Further screening with R package “UpSet” resulted in 8 Hub genes, including VEGFA, FYN, and IL6R which showed down-regulated expression in OA, and NGF, PTN, FGF1, NRP1, and GREM1 which showed up-regulated expression. Their expression in the GSE169077 microarray dataset was visualized as Heatmap (Figure 8A,B). Expression of the 8 Hub genes was validated in the GSE55235 and GSE55457 datasets and visualized as Boxplot and Heatmap (Figure 8C,D). Consistently, VEGFA, FYN, and IL6R had significantly lower expression in OA tissues than normal tissues, while NRP1 and GREM1 reversely had remarkably higher expression in OA tissues. It could be inferred that VEGFA, FYN, IL6R, NRP1 and GREM1 could be employed as candidate biomarkers for diagnosis of OA.

Figure 6. PPI network establishment. (A) PPI of the 42 IODEGs analyzed with the STRING database; (B) PPI network visualized by the Cytoscape, where red nodes represent up-regulated genes and green nodes represent down-regulated genes.

Figure 7. Screening of Hub genes based on the 8 algorithms as implemented by the plug-in CytoHubba. (A) MCC; (B) Degree; (C) EPC; (D) EcCentricity; (E) BottleNeck; (F) Closeness; (G) Radiality; (H) Betweenness. The rank of protein is indicated by colour intensity (from yellow to red), i.e. a higher rank of protein is indicated by a darker red colour.

Figure 8. Identification and validation of Hub genes. (A) Identification of Hub genes with R package “UpSet”; (B) Heatmap of the 8 Hub genes in the GSE169077 datasets; (C,D) Boxplot and Heatmap of the 8 Hub genes in the GSE55235 and GSE55457 datasets.

Assessment of the diagnostic effectiveness of candidate biomarkers

To validate the effectiveness of the 5 Hub genes in diagnosis of OA, ROC logistic regression analysis was performed based on the GSE55235 and GSE55457 datasets. AUC value >0.700 is considered as effective and highly reliable. It was found that the AUC values of VEGFA, FYN, GREM1, NRP1, and IL6R were 0.947, 0.795, 0.768, 0.775, and 0.776, respectively (Figure 9A–E), indicating good sensitivity and specificity of the 5 Hub genes in diagnosis of OA.

Figure 9. Validation of the 5 Hub genes in OA diagnosis in the GSE55235 and GSE55457 datasets with ROC curve analysis. (A) VEGFA; (B) FYN; (C) GREM1; (D) NRP1; (E) IL6R.

Relationship between diagnostic biomarkers and infiltrating immune cells in OA

Relationship between the 5 candidate diagnostic biomarkers (VEGFA, FYN, GREM1, NRP1, and IL6R) and 9 differential infiltrating immune cell types (aDC, Astrocytes, CD8+ naïve T cells, CMP, Macrophages, Macrophages M1, Neutrophils, NKT, and Th1 cells) was analyzed, and the significant relationship pairs were identified by r > 0.70 and p < .05 (Figure 10A). VEGFA was significantly positively associated with CD8+ naïve T cells (r = 0.852, p < .001) (Figure 10B), while VEGFA (r = −0.840, p = .001), FYN (r = −0.765, p = .006), and IL6R (r = −0.732, p = .010) were negatively associated with Macrophages M1 (Figure 10C–E). In addition, NRP1 was negatively associated with NKT (r = −0.764, p = .009) (Figure 10F), while NRP1 (r = 0.700, p = .021) and GREM1 (r = 0.764, p = .009) were positively associated with aDC (Figure 10G,H).

Figure 10. Relationship between candidate diagnostic biomarkers and differential infiltrating immune cell types. (A) Matrixplot of correlation coefficients; (B) VEGFA-CD8+ naïve T cells; (C) VEGFA-Macrophages M1; (D) FYN-Macrophages M1; (E) IL6R-Macrophages M1; (F) NRP1-NKT; (G) NRP1-aDC; (H) GREM1-aDC.

Discussion

OA is a chronic irreversible degenerative disease that might be driven by persistent inflammation, which is resulted from congenital and adaptive immune responses and can cause irreversible synovial damage and cartilage breakdown leading to joint dysfunction [23,24]. Multiple studies have noted that activation and interaction of infiltrating immune cells are critical for regulation of synovitis formation and persistent cartilage damage in OA. Therefore, regulation of immune response is necessary for OA treatment [25–27]. Nonetheless, no effective biomarkers for OA diagnosis are available at present due to the complicated pathogenesis of OA. High-throughput sequencing technique allows for better exploring the mechanisms underlying the occurrence and development of OA as well as looking for new effective diagnostic biomarkers. In this study, we combined the high-throughput sequencing with bioinformatics to analyse the role of infiltrating immune cells in OA and identify new effective diagnostic biomarkers.

At first, the GSE169077 dataset was obtained for analysis of infiltrating immune cells. We found that Astrocytes, M1 Macrophages, Macrophages, Common Myeloid Progenitors (CMP), Neutrophils, and aDC were more abundant in OA samples than normal samples, while CD8+ naïve T cells, Th1 cells, and NKT were less abundant in OA samples. It has been proven that macrophages and the various factors they secrete are important in regulation of joint synovial homeostasis and cartilage functioning [28]. Upon dysregulation of macrophage function in synovial homeostasis, more macrophages will be recruited to the joint synovium resulting in exacerbated synovial inflammation. In the meantime, there is a trend of macrophage polarization towards the M1 phenotype in OA, which also plays a pro-inflammatory role [29,30]. Some studies found that suppression of macrophage M1 polarization could slow the progression of OA, and that promotion of polarization towards the M2 phenotype via the NF-κB pathway was conducive to recovery from OA [31–33]. M1 macrophages are pro-inflammatory that can secrete cytokines to promote the recruitment of neutrophils and the proliferation and activation of synovial fibroblasts. The neutrophils and fibroblasts exert synergistic effects to potentiate persistent development of inflammation [34]. T-cell polarization, especially CD4+ T-cell differentiation to Th1 cell, and the abnormal immune response it mediates are also involved in the occurrence and development of OA [35]. It was report that CD8+ naïve T cells and NKT cells significantly changed in the peripheral blood, synovial fluid and synovium in OA patients and played a role in OA onset [36,37]. Dendritic cells (DC) are heterogeneous, specialized antigen presenting cells that are vital in the congenital and adaptive immunity. Mature DCs can induce the production of a large amount of inflammatory factors through increasing the expression of Toll-like receptors, promoting T-lymphocyte proliferation and subsequently leading to synovitis and cartilage degeneration [38]. Our findings are basically consistent with the studies mentioned above. However, role of these immune cells in OA was only highlighted by bioinformatics, requiring further experimental validation.

To obtain more effective diagnostic biomarkers, 382 DEGs between OA and normal samples were obtained, including 42 IODEGs. GO and KEGG pathway analyses demonstrated that the IODEGs were mainly involved in inflammatory and immune responses, suggesting that the function of these IODEGs might be associated with immune cells. Subsequently, 5 Hub genes, including GREM1, NRP1, VEGFA, FYN and IL6R, were identified, and as validated by the test sets and ROC logistic regression analysis, they had good diagnostic performance. This finding made our study results more reliable. We speculated that the 5 Hub genes may play an important regulatory role in pathogenesis of OA.

GREM1 (Gremlin 1, DAN Family BMP Antagonist) can modulate the bone morphogenetic protein (BMP) signalling pathway by serving as a BMP antagonist. Previous research showed significant up-regulation of GREM1 in chondrocytes in human OA, which affected tissue remodelling by regulating the BMP activity. Additionally, it was also validated in animal models that GREM1 was positively associated with disease severity [39,40]. Another study revealed that the persistent inflammation and immune microenvironment in OA could also lead to GREM1 elevation, resulting in synovial hyperplasia and long-term inflammation [41]. Moreover, the latest mechanistic study reported that the dysregulation of circ_SPG11 in OA could promote the IL-1β-induced chondrocyte apoptosis and extracellular matrix (ECM) degradation via increasing GREM1 expression [42]. NRP1 (Neuropilin 1) is the binding receptor of Semaphorin 3 A (SEMA3A), while their binding is critical to the immunity and bone metabolism [43]. Early research showed that both NRP1 and SEMA3A were highly expressed in OA cartilage, and their synergistic effects were suppressive for the effect of VEGF165 in promoting migration of chondrocytes in early OA [44]. In addition, high expression levels of NRP1 and SEMA3A could also inhibit the early-stage inflammation in OA [45]. Other scholars found that persistent low-grade inflammation in OA could lead to persistent high expression of NRP1 and its ligand protein upon IL-1β and TNF-α stimulation, subsequently promoting apoptosis in chondrocytes through the PI3K/Akt signalling pathway [46]. Combining the studies, more experiments are warranted to clarify the role of NRP1 in OA. VEGFA (Vascular Endothelial Growth Factor A) is a homodimeric vasoactive glycoprotein critical to angiogenesis and important for inflammatory response [47]. It was reported that VEGFA was an important regulator of angiogenesis during the process of endochondral bone formation. It could promote the formation and maturation of blood vessels with the cartilage, regulate the activity and maintain the survival of osteoblasts, thereby potentiating the formation of the new bone [48–50]. In addition, under the regulation of upstream miR-214-3p and miR-140-5p, VEGFA could also participate in the inflammatory response in OA through reducing the cartilage and synovial damages in OA [51,52]. FYN (FYN Proto-Oncogene, Src Family Tyrosine Kinase) is one of the tyrosine-specific phosphotransferases. It undergoes autophosphorylation upon stimulation by multiple transmembrane receptors and plays a role in biological processes, such as inflammation and cell survival, by mediating various signalling pathways [53]. A recent study suggested that effective reduction in FYN and β-catenin protein expression could decease the level of enzymes involved in ECM catabolism within the bone, thereby delaying the exacerbation of OA [54]. IL-6R is the receptor protein of pro-inflammatory cytokine IL-6, and it plays an important regulatory role in inflammatory and immune responses. Research revealed that application of Let-7a inhibitor in lipopolysaccharide-induced chondrocytes could significantly reduce the expression of inflammation-related factors and increase the expression of IL-6R to promote chondrocyte viability and then decrease apoptosis [55]. In addition, studies have also confirmed the related pathways of 5 Hub genes and their related effects. It was found that exosomal miR-150 derived from BMSCs inhibited TNF-α-mediated osteoblast apoptosis by regulating GREM1/NF-κB signalling pathway, thereby improving the symptoms of osteonecrosis of the femoral head [56]. In breast cancer metastasis experiments, FOXO3a-driven ectopic expression of miR-29b-2/miR-338 could inhibit EMT, migration, invasion and in vivo metastasis of breast cancer by regulating VEGFA/NRP1 signalling [57]. Yue et al. found that resveratrol coould also mediate the activation of NRF2 through the AKT/GSK3β/FYN signalling axis, significantly improving the osteoblast dysfunction induced by high glucose [58]. Further studies had shown that circ_0016347 could regulate the function and action of osteosarcoma cells by inactivating the JAK2/STAT3 signalling pathway via the miR-661/IL6R axis [59]. However, its concrete effects and correlation pathway mechanisms in the initiation and progression of OA need to be confirmed by rigorous experimental research. The current study identified that the 5 Hub genes had a potential role as a predictor of OA. Therefore, the 5 Hub genes could be considered as potential effective diagnostic biomarkers for OA.

Since both the 5 Hub genes and infiltrating immune cells identified before have a significant regulatory role in OA, we further analysed their relationship. We noted that VEGFA (r = −0.840, p = .001), FYN (r = −0.765, p = .006) and IL6R (r = −0.732, p = .010) were negatively associated with M1 macrophage. Existing studies proved that the VEGF165b secreted by macrophage would competitively inhibit VEGFA expression and downstream pathway activation, inducing more M1 macrophages and further suppressing angiogenesis and perfusion recovery of ischaemic muscle [60]. Likewise, activation of FYN-related pathways could effectively improve the Nrf2-mediated redox homeostasis, thereby inhibiting macrophage infiltration [61]. A study on breast cancer showed that suppression of pro-oncogenic gene MCT-1 up-regulated miR-34a expression, which inhibited IL-6R expression and promoted M1 polarization [62]. This study also noted that NRP1 was positively associated with aDC (r = 0.700, p = .021) but negatively associated with NKT (r = −0.764, p = .009). In an environment rich in vascular endothelial factors, high NRP1 expression activates downstream related pathways, promoting the maturation and activation of DCs [63]. In addition, NRP1 as the surface marker of Tregs is conducive to long-term interactions between Tregs and immature DCs [64]. Stradner MH et al. [65]found that Id2 deficiency could induce a hyporesponsive state of iNKT cells and then maintain NRP1 at a high level, thereby playing an anti-inflammatory role. Other relationship pairs we found here were first reported as potentially linked to the occurrence and development of OA, which provides new insights into relevant research fields. Collectively, we reasoned that the 5 diagnostic biomarkers may be a part of OA pathogenesis via interacting with related immune cells. Nevertheless, further experiments are in demand for validation.

In spite of the findings, this study still has some limitations. First, this is a retrospective study, and prospective studies with more samples and methods are required to validate the findings of the current study. Second, clinical information is not covered in this study. Besides, the upstream and downstream regulatory pathways of the identified biomarkers need to be further investigated. Finally, in-vivo and ex-invivo experiments, including experiments in human samples and animals, should be performed before clinical translation of the findings of the study.

To sum up, this study made a detailed analysis on OA-related immune cell infiltration and identified 5 potentially effective diagnostic biomarkers for OA (GREM1, NRP1, VEGFA, FYN and IL6R). Through a series of bioinformatics analyses, we established multiple biomarker-infiltrating immune cell relationship pairs, including VEGFA-CD8+ naïve T cells, VEGFA-Macrophages M1, FYN-Macrophages M1, IL6R-Macrophages M1, NRP1-NKT, NRP1-aDC, and GREM1-aDC. We believed that these relationship pairs played their role during the occurrence and development of OA via participating in immune responses. In all, the current study provides new insights into the immunomodulatory mechanism of OA and also a new direction for research into the diagnosis and treatment of OA.

Author contributions

Jiachen Peng, Tao Zhang and Piaotao Cheng had substantial contributions to the conception or design of the experiments and supervised the study. Piaotao Cheng, Shouhang Gong, Caopei Guo, Ping Kong, Chencheng Li and Chengbing Yang were involved in the acquisition, analysis or interpretation of data for the work. Jiachen Peng, Tao Zhang, Piaotao Cheng, Shouhang Gong, Caopei Guo, Ping Kong, Chencheng Li and Chengbing Yang drafted the work and revised it critically for important intellectual content and finally approved the version to be published. All authors agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Acknowledgments

We acknowledge the GEO database and SangerBox website for providing their platforms as well as their contributors for uploading meaningful datasets.

Disclosure statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Data availability statement

The GSE169077, GSE55235 and GSE55457 datasets were publicly available and obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The names of the repository and accession number(s) can be found in the online repositories.

Additional information

Funding

This work was supported by grants from the National Natural Science Foundation of China [No.81760400].