Machine learning-based endoplasmic reticulum-related diagnostic biomarker and immune microenvironment landscape for osteoarthritis

Tingting Liu; Xiaomao Li; Mu Pang; Lifen Wang; Ye Li; Xizhe Sun

doi:10.18632/aging.205611

Research Paper Volume 16, Issue 5 pp 4563—4578

Machine learning-based endoplasmic reticulum-related diagnostic biomarker and immune microenvironment landscape for osteoarthritis

Tingting Liu^1, , Xiaomao Li^2, , Mu Pang^3, , Lifen Wang^1, , Ye Li^4, , Xizhe Sun^1, ,

¹ Research Center for Drug Safety Evaluation of Hainan, Hainan Medical University, Haikou, Hainan 571199, China
² Jiangsu Food and Pharmaceutical Science College, Huaian, Jiangsu 223023, China
³ The Fourth Clinical Medical College of Guangzhou University of Chinese Medicine (Shenzhen Traditional Chinese Medicine Hospital), Shenzhen, Guangdong 518000, China
⁴ Chongqing Three Gorges Medical College, Chongqing 404120, China

Received: October 13, 2023 Accepted: January 23, 2024 Published: February 28, 2024

https://doi.org/10.18632/aging.205611
How to Cite

Copyright: © 2024 Liu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Osteoarthritis (OA) is the most common degenerative joint disease worldwide. Further improving the current limited understanding of osteoarthritis has positive clinical value.

Methods: OA samples were collected from GEO database and endoplasmic reticulum related genes (ERRGs) were identified. The WGCNA network was further built to identify the crucial gene module. Based on the expression profiles of characteristic ERRGs, LASSO algorithm was used to select key factors according to the minimum λ value. Random forest (RF) algorithm was used to calculate the importance of ERRGs. Subsequently, overlapping genes based on LASSO and RF algorithms were identified as ERRGs-related diagnostic biomarkers. In addition, OA specimens were also collected and performed qRT-PCR quantitative analysis of selected ERRGs.

Results: We identified four ERRGs associated with OA risk assessment through machine learning methods, and verified the abnormal expressions of these screened markers in OA patients through in vitro experiments. The influence of selected markers on OA immune infiltration was also evaluated.

Conclusions: Our results provide new evidence for the role of ER stress in the OA progression, as well as new markers and potential intervention targets for OA.

Introduction

Osteoarthritis (OA) is the most common degenerative joint disease worldwide, characterized by symptoms such as cartilage degradation, pain and limited movement, accompanied by a significant decline in the patient’s quality of life [1]. There are currently no pharmacological methods to inhibit disease progression or reduce cartilage damage. OA treatment strategies are limited to adjunctive therapy such as joint replacement surgery and physical therapy to improve function [2]. Further improving the current limited understanding of OA has positive clinical value.

Due to the heterogeneity of OA, individual etiology, clinical manifestations and response to treatment were inconsistent [3, 4]. These present great challenges for the prevention, diagnosis and prognosis prediction of OA. As objective, quantifiable features, biomarkers contribute to disease risk assessment by analyzing measurement reliability and summarizing biological, physiological or pathological pathways [5]. Currently, the only clinically used OA biomarkers are imaging markers. However, imaging markers cannot detect molecular alterations prior to the appearance of structural changes or provide potential therapeutic targets [6]. The predictive reliability of currently discovered OA biomarkers related to cartilage/bone structure, inflammation, and metabolism also requires further careful consideration [7–9]. At the same time, the reliability of multiple marker combinations in predicting OA severity and progression was significantly higher than that of a single marker [10, 11]. Therefore, there is an urgent need to screen for new meaningful diagnostic and prognostic markers.

Endoplasmic reticulum (ER) stress is caused by impaired protein folding capacity in ER, leading to accumulation of incorrectly folded proteins in ER, which adversely affects cell physiological function [12]. There is a lot of evidence showing the role of ER stress in the development of OA. Chondrodysplasia is usually caused by mutations in genes that code for cartilage components. This genetic mutation causes impaired synthesis or secretion of extracellular matrix (ECM) components that make up the main part of cartilage and aggregate in ER, ultimately leading to loss of ECM and disruption of chondrocyte homeostasis [13–15]. As two major risk factors for OA, aging and obesity both impair the function of key ER molecular chaperones, leading to improper protein folding and ER stress [16]. Thus, ER stress-induced chondrocyte death has been identified as a contributing factor to OA and a potentially reliable treatment strategy [17]. It is necessary to further explore the mechanism and clinical value of ER stress in OA.

The correlation between endoplasmic reticulum stress and immune infiltration has also been gradually clarified. Further activation of the unfolded protein response (UPR) during endoplasmic reticulum stress is characteristic of many autoimmune diseases [18]. Microenvironmental stress in the immune-infiltrating microenvironment, including hypoxia, reactive oxygen species, and pro-inflammatory cytokines, may increase the level of endoplasmic reticulum stress in dendritic cells (DC) and fibroblast-like synovial cells (FLS) in joints [19]. There is evidence that endoplasmic reticulum stress plays a role in the development process from B cells to plasma cells and the secretion of immunoglobulin [20, 21]. In addition, accumulation of misfolded proteins during ER stress can lead to increased MHC presentation on the cell surface, thereby increasing the chance of auto-reactive T cell activation [18]. Several pro-inflammatory cytokines have also been reported to act through endoplasmic reticulum stress processes [22, 23]. Therefore, the influence of endoplasmic reticulum stress on the immune infiltration of adaptive immune cells and related genes are worthy of further exploration.

In this study, we identified four endoplasmic reticulum related genes (ERRGs) associated with OA risk assessment through machine learning methods, and verified the abnormal expression of these screened markers in OA patients through in vitro experiments. At the same time, we evaluated the effect of markers on OA immune infiltration. Our results provide new evidence for the role of ER stress in the progression of OA, as well as new markers and potential intervention targets for OA.

Materials and Methods

Data collection and pre-processing

Three open access datasets of OA and HC samples were collected from the GEO database, including GSE51588, GSE98918 and GSE117999. The “limma” script was used to preprocess the raw data of each GEO dataset and “SVA” script was utilized to normalize the raw data and eliminate the batch effect of the three GEO datasets in the R language environment [24]. In this study, the endoplasmic reticulum related genes (ERRGs) were identified based on the GeneCards database [25] (Supplementary Table 1).

Differential expression analysis of transcriptome data and molecular pathways enrichment prediction

After the normalization of the transcriptome data of HC and OA samples, we conducted a differential expression analysis via “limma” script in the R environment. Under the selection criteria of p.adjust < 0.05, the DEGs between the OA and HC groups was identified. Based on the SRTING database, the potential interaction of ERRGs was revealed using the “string” script. Moreover, we used the “clusterProfiler” package to predict the potential molecular function and KEGG pathways of the ERRGs.

Construction of WGCNA network to identify the crucial gene module for OA

A WGCNA network model was established to identify the crucial gene module for OA. In first, according to the transcriptome data of all samples (HC samples and OA samples), we developed a clustering tree to exclude the abnormal samples. Next, based on the optimal soft threshold (β), we established a WGCNA network. After the exclusion filter of expression level was set at 0.5, the genes were divided into the different gene modules and the dynamic tree cut method to merge the similar gene modules. Pearson correlation was used to evaluate the relationship between each gene module and estimate the potential correlation between clinical trait and each gene module. Finally, the most relevant gene module was chosen for the subsequent analysis.

Generation of ERRG related diagnostic biomarker based on the machine learning algorithm

Two unique machine learning algorithms were utilized to identify the ERRGs related diagnostic biomarkers. Based on the expression profiler of the feature ERRGs related genes, a LASSO algorithm was performed to select the pivotal ERRG related genes according to the minimum λ value. Random forest (RF) algorithm was conducted to calculate the importance of each feature ERRGs related genes. Subsequently, the overlapping genes based on LASSO and RF algorithms were identified as the ERRGs related diagnostic biomarker.

Model diagnostic effectiveness evaluation of ERRGs related biomarker and nomogram

R package “rms” was used to establish the nomogram model based on the expression profiler of ERRGs related biomarkers. The formula for nomogram was: nomogram score = HSPA5 × −5.4 + UBL4A × 5.2 + ATF4 × −4.9 + PPP1R15A × −4.4. Utilizing the “pROC” script, the ROC curve was carried out to evaluate the model diagnostic effectiveness of the ERRGs related biomarkers for OA.

Characterization of immune microenvironment landscape

According to the 22 immune cell gene markers signature, we used the CIBERSORT estimation algorithm to evaluate the immune cell relative percent of OA and HC samples. In addition, Pearson correlation method was performed to estimate the potential correlation of ERRGs related biomarkers and 22 immune cells.

qRT-PCR quantitative analysis of ERRGs related biomarkers

Real-time fluorescence quantitative qPCR was performed to detect the expression levels of HSPA5, UBL4A, ATF4, and PPP1R15A. Frozen OA tissues were taken out at −80°C, and total RNA was extracted, and all the steps were strictly in accordance with the instructions of the RNA Extraction Kit (Qiagen, Germany). The RNA concentration of each sample was determined by Nanodrop 2000, and RNA was reverse transcribed into cDNA according to the Reverse Transcription Kit (TaKaRa, Japan). cDNA was then used as a template for amplification of HSPA5, UBL4A, ATF4 and PPP1R15A on a Bio-Rad CFX90 Real-time PCR instrument. qRT-PCR reaction conditions were as follows: 95°C pre-denaturation for 30 s; 95°C denaturation for 5 s, RT-PCR reaction conditions were as follows: 95°C pre-denaturation for 30 s; 95°C denaturation for 5 s, and 95°C denaturation for 5 s. The reaction was performed at the same time as the reaction. The reaction conditions were: pre-denaturation at 95°C for 30 s; denaturation at 95°C for 5 s; annealing at 56°C for 5 s; and extension at 65°C for 5 s. A total of 39 cycles were performed. The relative expression of target genes was analyzed by using 2^−ΔΔCt (ΔCt = Ct value of target gene - Ct value of internal reference), and GAPDH was used as internal reference (Supplementary Table 2).

Statistical analysis

All statistical analysis were performed under the R software 4.1.0 (https://cran.r-project.org/) and Perl language environment. Statistical comparison of data between the two groups was calculated using “limma” R package (Wilcox rank-sum test). p < 0.05 was considered to be statistically significant.

Data availability statement

The bioinformatic datasets presented in this study can be found in online repositories. The names of the online repositories and accession numbers can be found in the article.

Results

Data processing and differential expression analysis of OA and HC groups

We extracted a total of 25 HC samples and 25 OA samples from the GEO database to investigate the potential function of ERAGs in OA. Utilizing the “sva” software, the batch effect of each sample was removed and normalized (Figure 1A, 1B). According to the “limma” package, the differential expression analysis was performed with the analysis standards set at p.adjust < 0.05 and |FC|>1 (Figure 1C). The heatmap reveals the most significantly DEGs between the OA and HC groups (Figure 1D).

The workflows of data processing and differential expression analysis between HC and OA groups. (A, B) Data pre-processing of HC and OA samples in GEO database. (C) Identification of the DEGs between the HC and OA samples, the red dot indicates the up-DEGs and the blue dots indicate the down-DEGs in OA. The standard for selecting DEGs is set at p.adjust D) The expression analysis of DEGs between HC and OA groups.

Figure 1. The workflows of data processing and differential expression analysis between HC and OA groups. (A, B) Data pre-processing of HC and OA samples in GEO database. (C) Identification of the DEGs between the HC and OA samples, the red dot indicates the up-DEGs and the blue dots indicate the down-DEGs in OA. The standard for selecting DEGs is set at p.adjust < 0.05. (D) The expression analysis of DEGs between HC and OA groups.

Identification of the pivotal gene module associated with OA via WGCNA

A total of 25 HC and 25 OA samples were enrolled to identify the pivotal gene module for OA based on the WGCNA. Firstly, the samples were clustered to exclude the abnormal samples. Based on the filter condition of scale-free topology (R²) set at >0.85, we constructed a WGCNA network with the soft threshold (power) selected as β = 6 (Figure 2A). With the height of gene modules set at 0.25, the feature genes were divided into the 25 inimitable gene modules (Figure 2B). Correlation analysis suggested the stable independence between the 25 inimitable gene modules (Figure 2C). Thereafter, we further evaluated the association between the clinical trait and gene modules and the result of module-trait relationships illustrated that the brown gene module was positively associated with HC but negatively associated with OA (Figure 2D). The scatter plot of brown gene module illustrated a strong correlation between the module membership and gene significance (r = 0.92, p < 1e-200), and was selected for the subsequent analysis (Figure 2E).

Establishment of WGCNA model for selecting the pivotal gene module in OA. (A) Scale independence and mean connectivity. (B) The height of different gene modules and dynamic tree cut. (C) Potential association of 25 unique gene modules. (D) Correlation analysis of 25 unique gene modules and clinical features. (E) The relationship of module membership and gene significance in brown module.

Figure 2. Establishment of WGCNA model for selecting the pivotal gene module in OA. (A) Scale independence and mean connectivity. (B) The height of different gene modules and dynamic tree cut. (C) Potential association of 25 unique gene modules. (D) Correlation analysis of 25 unique gene modules and clinical features. (E) The relationship of module membership and gene significance in brown module.

Analysis of the pivotal DE-ERRGs and key molecular pathways

Utilizing the differential expression analysis and WGCNA network (brown module), we observed the 10 overlapping DE-ERRGs were considered as the pivotal DE-ERRGs for OA (Figure 3A). The PPI network revealed a strong relationship of 10 pivotal DE-ERRGs (Figure 3B). To further investigate the potential molecular function of pivotal DE-ERRGs, we performed the GO ang KEGG enrichment analysis. GO enrichment results suggested that the 10 pivotal DE-ERRGs was enriched in response to unfolded protein, response to topologically incorrect protein, smooth endoplasmic reticulum, mitochondrial outer membrane and chaperone binding, while the KEGG enrichment analysis revealed that protein processing in endoplasmic reticulum, amyotrophic lateral sclerosis and pathways of neurodegeneration- multiple diseases may medicated the role of DE-ERRGs in OA (Figure 3C, 3D).

Identification of pivotal DE-ERRGs and molecular pathway enrichment analysis. (A) Selection of pivotal DE-ERRGs via differential expression analysis and WGCNA. (B) PPI network analysis of 10 pivotal DE-ERRGs. (C, D) Molecular function analysis of 10 pivotal DE-ERRGs.

Figure 3. Identification of pivotal DE-ERRGs and molecular pathway enrichment analysis. (A) Selection of pivotal DE-ERRGs via differential expression analysis and WGCNA. (B) PPI network analysis of 10 pivotal DE-ERRGs. (C, D) Molecular function analysis of 10 pivotal DE-ERRGs.

Characteristic DE-ERRGs biomarker identification for OA

LASSO and RF machine learning algorithms were performed to identify the characteristic ERRGs-related biomarkers for OA. Based on the expression of 10 pivotal DE-ERRGs, the LASSO algorithm revealed the coefficients of each variable and 8 feature variables was selected with the minimum log lambda (λ = 8) (Figure 4A). According to the RF algorithm, the 10 pivotal DE-ERRGs were ranked according the variable importance and 5 important variables were identified (Figure 4B). By Venn diagram, 4 overlapping DE-ERRGs were considered as the feature biomarkers for OA, involving HSPA5, UBL4A, ATF4 and PPP1R15A (Figure 4C). The Pearson correlation analysis of 4 ERGGs feature biomarkers illustrated that HSPA5 was negatively associated with UBL4A, while was positively associated with ATF4 and PPP1R15A; UBL4A was negatively correlated with ATF4 and PPP1R15A; while ATF4 was positively linked with PPP1R15A (Figure 4D).

Feature ERRGs biomarkers identification using LASSO and RF machine learning algorithm. (A) LASSO algorithm for selecting feature ERRGs related biomarkers. (B) The importance ranking of 10 pivotal DE-ERRGs via RF algorithm. (C) Identification of feature ERRGs biomarkers via RF and LASSO machine learning algorithms. (D) Pearson correlation analysis of HSPA5, UBL4A, ATF4 and PPP1R15A.

Figure 4. Feature ERRGs biomarkers identification using LASSO and RF machine learning algorithm. (A) LASSO algorithm for selecting feature ERRGs related biomarkers. (B) The importance ranking of 10 pivotal DE-ERRGs via RF algorithm. (C) Identification of feature ERRGs biomarkers via RF and LASSO machine learning algorithms. (D) Pearson correlation analysis of HSPA5, UBL4A, ATF4 and PPP1R15A.

Model effectiveness assessment and nomogram construction of feature ERRGs biomarkers

We further evaluated the diagnostic effectiveness of the ERRGs related biomarkers in OA and the expression profile results suggested that the expression of HSPA5, ATF4 and PPP1R15A in HC group was greatly overexpressed, while the expression of UBL4A was significantly overexpressed in OA group (Figure 5A–5D). According to the expression profile of the four ERRGs related biomarkers, we established a newly nomogram model to evaluate the diagnostic effectiveness for OA (Figure 5E). The ROC analysis result indicated that the AUC of HSPA5, UBL4A, ATF4 and PPP1R15A was 0.805, 0.776, 0.882 and 0.955, respectively. Notably, the AUC of nomogram model was 0.987, which was higher than HSPA5, UBL4A, ATF4 and PPP1R15A, illustrating a satisfactory diagnostic ability of nomogram model for OA (Figure 5F).

Figure 5. Diagnostic effectiveness evaluation and nomogram construction based on the ERRGs related biomarkers. (A–D) The expression profile analysis of HSPA5, UBL4A, ATF4 and PPP1R15A in HC and OA groups. (E) Nomogram construction based on the four ERRGs related biomarkers. (F) Diagnostic effectiveness evaluation of HSPA5, UBL4A, ATF4, PPP1R15A and nomogram score.

Immune microenvironment characteristic and GSEA analysis

Utilizing the signature of 22 immune cell subtypes, we evaluated the proportion of 22 immune cell of each HC and OA sample by CIBERSORT algorithm (Figure 6A, 6B). The quantitative results of 22 immune cells indicated that the relative percent of plasma cells, NK cells activated, macrophages M1 and mast cells resting in OA group was significantly higher, whereas the relative percent of T cells CD4 memory resting, Dendritic cells activated and mast cells activated in HC group was remarkable up-regulated than the OA group. In HC group, the GSEA result suggested that adipocytokine signaling pathway, MAPK signaling pathway and NOD like receptor signaling pathway was up-regulated; however, the lysosome, allograft rejection and cell adhesion molecules cams was up-regulated in OA group (Figure 6C, 6D).

Estimation of immune microenvironment characteristic and KEGG related GSEA analysis. (A) The evaluation of 22 immune cell subtypes of HC and OA groups based on the CIBERSORT estimation algorithm. (B) The quantitative analysis of relative percent of 22 immune cell subtypes in HC and OA groups. (C, D) KEGG related pathway analysis in HC and OA groups based on the GSEA analysis.

Figure 6. Estimation of immune microenvironment characteristic and KEGG related GSEA analysis. (A) The evaluation of 22 immune cell subtypes of HC and OA groups based on the CIBERSORT estimation algorithm. (B) The quantitative analysis of relative percent of 22 immune cell subtypes in HC and OA groups. (C, D) KEGG related pathway analysis in HC and OA groups based on the GSEA analysis.

We further predicted the potential association of ERRGs related biomarkers and immune microenvironment characteristic using the Pearson correlation algorithm. As illustrated in Figure 7A–7D, we observed that ATF4 was positively associated with mast cells activated and negatively associated with mast cells resting, NK cells activated, plasma cells and B cells naïve; HSPA5 was positively correlated with dendritic cells activated and mast cells activated but negatively correlated with mast cells resting, plasma cells, NK cells activated and dendritic cells resting; PPP1R15A was positively correlated with T cells CD4 memory resting, dendritic cells activated and mast cells activated, while was negatively associated with mast cells resting, plasma cells and NK cells activated; UBL4A was positively correlated with dendritic cells resting and mast cells resting but negatively correlated with mast cells activated and dendritic cells activated.

Potential association analysis of ERRGs related biomarkers and immune microenvironment characteristic. The lollipop plot shows the association of immune microenvironment and (A) ATF4, (B) HSPA5, (C) PPP1R15A and (D) UBL4A.

Figure 7. Potential association analysis of ERRGs related biomarkers and immune microenvironment characteristic. The lollipop plot shows the association of immune microenvironment and (A) ATF4, (B) HSPA5, (C) PPP1R15A and (D) UBL4A.

In vitro validation of ERRGs related biomarkers

We further validated the expression profiler of the four ERRGs related biomarkers in the HC and OA samples. The expression results revealed that the HC group had higher expression of ATF4, HSPA5 and PPP1R15A, whereas the expression of UBL4A in OA group was significantly higher than HC group (Figure 8A–8D).

qRT-PCR analysis of the four ERRGs related biomarkers in OA and HC groups. The expression profiler of (A) ATF4, (B) HSPA5, (C) PPP1R15A and (D) UBL4A in HC and OA groups. *p **p ***p ****p

Figure 8. qRT-PCR analysis of the four ERRGs related biomarkers in OA and HC groups. The expression profiler of (A) ATF4, (B) HSPA5, (C) PPP1R15A and (D) UBL4A in HC and OA groups. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001.

Discussion

In this study, we screened four ERRGs with OA severity and progression, evaluated their effects on immune infiltration, and verified their abnormal expressions in OA patients through in vitro experiments.

Among the four ERRGs we screened, PPP1R15A had the most significant effect on OA. PPP1R15A is a key factor in the integrative stress response of mammals [26]. PPP1R15A mediates dephosphorylation of eIF2α [27, 28]. When eIF2α is phosphorylated, global protein synthesis is reduced, which is beneficial for cell survival and recovery [29]. In contrast, dephosphorylation of eIF2α allowed the cells to resume normal protein synthesis processes [30, 31]. Our results showed reduced PPP1R15A levels in patients with OA. PPP1R15A can direct the breakdown of the unfolded protein response (UPR) and lead to the restoration of normal ribosome activity. UPR is a direct result of endoplasmic reticulum (ER) stress [32]. Mice lacking PPP1R15A activity were healthier than wild-type controls, had improved insulin sensitivity and were more resistant to ER stress [33] and inhibited PPP1R15A to protect cells from ER stress-induced apoptosis [34]. At the same time, PPP1R15A’s promoting effect on Cd8+ T cells has also been reported [35], which is consistent with our observed immune infiltration results. Our data demonstrate the role of PPP1R15A in OA development from a new perspective.

UBL4A is essential for the mitochondrial fusion process under nutrient deprivation stress [36]. Meanwhile, as a ubiquitin ligase-associated protein, UBL4A is involved in proteasome degradation and mediates DNA damage signaling and cell death [37, 38]. UBL4A plays an important role in the development of immune dysfunction and subsequent abnormal bone metabolism. UBL4A contributes to the development of inflammatory diseases by regulating NF-CUMB signaling in macrophages and dendritic cells [39]. UBL4A knockout mice showed mild kyphosis and scoliosis with dysregulation of osteoblastogenesis and chondrogenesis [40]. UBL4A knockout mice resist collagen-induced arthritis by regulating the balance of Th1, Th17, and regulatory T cells in the T cell subpopulation [41]. We observed that the higher expression of UBL4A in OA patients may be related to the role of UBL4A in promoting inflammation and abnormal bone metabolism reported in the literature.

ATF4 is a member of the activating transcription factor (ATF)/cyclic adenosine phosphate Responsive element binding (CREB) family and plays a key role in the regulation of osteoblast function [42]. ATF4 accumulates in osteoblasts and is a specific activator of osteocalcin specific element 1 (OSE1), which is particularly active in osteoblasts [43, 44]. In addition to participating in amino acid metabolism, ATF4 plays an important role in type I collagen synthesis and transcriptional control of several major osteoblast genes [45, 46]. ATF4 is involved in various diseases related to bone metabolism, including Coffin-Lowry Syndrome [44]. In OA, AFT4 is associated with RUNt-associated transcription factor 2 (Runx2), which is essential for chondrocyte hypertrophy, and regulates osteoblast differentiation and chondro-development [47–49].

HSPA5, member A of the heat shock protein family (HSPA5), is a chaperone member mainly expressed in ER [50]. HSPA5 is involved in UPR process and promotes cell survival under ER stress [51]. Recently, HSPA5 was found to be a suppressor of the ferroptosis process [52, 53]. In OA, inhibition of HSPA5 expression can degrade GPX4, thus promoting ferroptosis in chondrocytes. However, upregulated HSPA5 expression could inhibit inflammatory damage and ferroptosis, thus alleviating OA progression [54]. Similarly, we observed that compared with the control group, the expression level of HSPA5 in OA patients was significantly reduced, which also suggested the role of HSPA5 and the ferroptosis involved in OA disease progression.

In the process of OA development, low-grade inflammatory processes formed by immune infiltration are involved. Our immunoinfiltration results showed that OA patients had significantly increased M1 macrophage infiltration levels compared with healthy controls. The number of macrophages in the synovial membrane of OA increased and was associated with OA disease progression and pain [55]. In addition, the disturbance of macrophage polarization in synovial tissue may contribute to the occurrence and progression of OA [56]. The degree of imbalance of the M1/M2 type of macrophage polarization is helpful in evaluating the severity of knee OA [57]. Further studies showed that the activated synovial macrophages in OA patients were mainly M1 macrophages, and the polarization of macrophages to M2 would weaken the development of OA, suggesting that the increase of M1 synovial macrophages may be the key reason for the deterioration of OA, which is consistent with the phenomenon we observed [58]. M1 macrophages can be induced by IFN-γ and TNF-α in vitro [59]. Because of its ability to produce pro-inflammatory cytokines including TNF-α and IL-1, M1 macrophages are known as pro-inflammatory macrophages [60]. The possible mechanism by which M1 macrophages affect OA is that M1 macrophages with increased content in OA synovium can secrete more cytokines. Subsequent chronic inflammation leads to cartilage degradation and osteophyte formation [61].

The results of DE-ERRGs pathway enrichment analysis showed that UPR was the most relevant pathway in biological process, suggesting the possible role of UPR in OA. UPR is induced by loading of unfolded or misfolded proteins that accumulate in ER and is intended to restore ER homeostasis by initiating apoptosis when ER stress persists [51]. UPR plays a crucial role in cartilage formation [62]. In chondrocytes and osteoblasts, bone morphogenetic protein 2 (BMP2) is an activator of UPR signaling [63]. Severe chondrodysplasia was observed in BMP2-knockout mice, accompanied by disturbances of growth plate chondrocytes [64]. In addition, the aforementioned ATF4 is also regulated by BMP2 and is involved in the expression of UPR transcription factors [63]. Targeting UPR has potential value in improving symptoms in cartilage-related diseases [65].

From targeted therapy to molecular classification and patient stratification to prognostic prediction, microarray has been shown to have important clinical applications in a variety of diseases including OA [66, 67]. In this study, we screened four ERRGs associated with OA progression through microarray analysis of a public database and verified their abnormal expression in OA patients. By analyzing the differentially expressed genes from different sources and involved pathways of OA, we can not only clarify the potential genetic mechanism of OA pathogenesis, but also explore potential OA targeted drug therapeutic targets [66, 68, 69]. Limited to the condition of public database, the selection bias of race and region existed in this study. In the future, further large-scale, multi-center analyses will better analyze the reliability and potential clinical application value of the selected targets. In addition, this study did not involve mechanism studies, and further molecular mechanism studies in the future can better elucidate the influence of different expressed genes on OA.

Supplementary Materials

Supplementary Table 1

Supplementary Table 2

Author Contributions

XS, YL and TL designed the study. XL, MP and LW collected and analyzed the data. XS, TL and XL verified experiment based on analytical results. TL wrote and edited the manuscript. All authors contributed to the article and approved the submitted version.

Conflicts of Interest

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.

Ethical Statement

No ethical approval was required for this study. Experiments used data available from open access datasets. This work was completed without further animal model or human participants.

Funding

The research was sponsored by Hainan Medical University and supported by the Education Department of Hainan Province, with the project number: Hnky2023-33. Furthermore, it received support from the Huai’an City Science and Technology Project (HABZ202226) and Chongqing Science and Technology Youth Project (KJQN202302721).

References

1. Lee SY, Wong PF, Jamal J, Roebuck MM. Naturally-derived endoplasmic reticulum stress inhibitors for osteoarthritis? Eur J Pharmacol. 2022; 922:174903. https://doi.org/10.1016/j.ejphar.2022.174903 [PubMed]
2. Buttgereit F, Burmester GR, Bijlsma JW. Non-surgical management of knee osteoarthritis: where are we now and where do we need to go? RMD Open. 2015; 1:e000027. https://doi.org/10.1136/rmdopen-2014-000027 [PubMed]
3. Costello CA, Rockel JS, Liu M, Gandhi R, Perruccio AV, Rampersaud YR, Mahomed NN, Rahman P, Randell EW, Furey A, Kapoor M, Zhai G. Individual participant data meta-analysis of metabolomics on sustained knee pain in primary osteoarthritis patients. Rheumatology (Oxford). 2023; 62:1964–71. https://doi.org/10.1093/rheumatology/keac545 [PubMed]
4. Javaid MK, Kiran A, Guermazi A, Kwoh CK, Zaim S, Carbone L, Harris T, McCulloch CE, Arden NK, Lane NE, Felson D, Nevitt M, and Health ABC Study. Individual magnetic resonance imaging and radiographic features of knee osteoarthritis in subjects with unilateral knee pain: the health, aging, and body composition study. Arthritis Rheum. 2012; 64:3246–55. https://doi.org/10.1002/art.34594 [PubMed]
5. Califf RM. Biomarker definitions and their applications. Exp Biol Med (Maywood). 2018; 243:213–21. https://doi.org/10.1177/1535370217750088 [PubMed]
6. Liu M, Haque N, Huang J, Zhai G. Osteoarthritis year in review 2023: metabolite and protein biomarkers. Osteoarthritis Cartilage. 2023; 31:1437–53. https://doi.org/10.1016/j.joca.2023.08.005 [PubMed]
7. Angelini F, Widera P, Mobasheri A, Blair J, Struglics A, Uebelhoer M, Henrotin Y, Marijnissen AC, Kloppenburg M, Blanco FJ, Haugen IK, Berenbaum F, Ladel C, et al. Osteoarthritis endotype discovery via clustering of biochemical marker data. Ann Rheum Dis. 2022; 81:666–75. https://doi.org/10.1136/annrheumdis-2021-221763 [PubMed]
8. Rockel JS, Layeghifard M, Rampersaud YR, Perruccio AV, Mahomed NN, Davey JR, Syed K, Gandhi R, Kapoor M. Identification of a differential metabolite-based signature in patients with late-stage knee osteoarthritis. Osteoarthr Cartil Open. 2022; 4:100258. https://doi.org/10.1016/j.ocarto.2022.100258 [PubMed]
9. Xie Z, Aitken D, Liu M, Lei G, Jones G, Cicuttini F, Zhai G. Serum Metabolomic Signatures for Knee Cartilage Volume Loss over 10 Years in Community-Dwelling Older Adults. Life (Basel). 2022; 12:869. https://doi.org/10.3390/life12060869 [PubMed]
10. Loef M, van de Stadt L, Böhringer S, Bay-Jensen AC, Mobasheri A, Larkin J, Lafeber FPJ, Blanco FJ, Haugen IK, Berenbaum F, Giera M, Ioan-Facsinay A, Kloppenburg M. The association of the lipid profile with knee and hand osteoarthritis severity: the IMI-APPROACH cohort. Osteoarthritis Cartilage. 2022; 30:1062–9. https://doi.org/10.1016/j.joca.2022.05.008 [PubMed]
11. Liem Y, Judge A, Li Y, Sharif M. Biochemical, clinical, demographic and imaging biomarkers for disease progression in knee osteoarthritis. Biomark Med. 2022; 16:633–45. https://doi.org/10.2217/bmm-2021-0579 [PubMed]
12. Tabas I, Ron D. Integrating the mechanisms of apoptosis induced by endoplasmic reticulum stress. Nat Cell Biol. 2011; 13:184–90. https://doi.org/10.1038/ncb0311-184 [PubMed]
13. Kung LH, Rajpar MH, Briggs MD, Boot-Handford RP. Hypertrophic chondrocytes have a limited capacity to cope with increases in endoplasmic reticulum stress without triggering the unfolded protein response. J Histochem Cytochem. 2012; 60:734–48. https://doi.org/10.1369/0022155412458436 [PubMed]
14. Hartley CL, Edwards S, Mullan L, Bell PA, Fresquet M, Boot-Handford RP, Briggs MD. Armet/Manf and Creld2 are components of a specialized ER stress response provoked by inappropriate formation of disulphide bonds: implications for genetic skeletal diseases. Hum Mol Genet. 2013; 22:5262–75. https://doi.org/10.1093/hmg/ddt383 [PubMed]
15. Gualeni B, Rajpar MH, Kellogg A, Bell PA, Arvan P, Boot-Handford RP, Briggs MD. A novel transgenic mouse model of growth plate dysplasia reveals that decreased chondrocyte proliferation due to chronic ER stress is a key factor in reduced bone growth. Dis Model Mech. 2013; 6:1414–25. https://doi.org/10.1242/dmm.013342 [PubMed]
16. Mazor M, Best TM, Cesaro A, Lespessailles E, Toumi H. Osteoarthritis biomarker responses and cartilage adaptation to exercise: A review of animal and human models. Scand J Med Sci Sports. 2019; 29:1072–82. https://doi.org/10.1111/sms.13435 [PubMed]
17. Rellmann Y, Eidhof E, Dreier R. Review: ER stress-induced cell death in osteoarthritic cartilage. Cell Signal. 2021; 78:109880. https://doi.org/10.1016/j.cellsig.2020.109880 [PubMed]
18. Hasnain SZ, Lourie R, Das I, Chen AC, McGuckin MA. The interplay between endoplasmic reticulum stress and inflammation. Immunol Cell Biol. 2012; 90:260–70. https://doi.org/10.1038/icb.2011.112 [PubMed]
19. Park YJ, Yoo SA, Kim WU. Role of endoplasmic reticulum stress in rheumatoid arthritis pathogenesis. J Korean Med Sci. 2014; 29:2–11. https://doi.org/10.3346/jkms.2014.29.1.2 [PubMed]
20. Gass JN, Gunn KE, Sriburi R, Brewer JW. Stressed-out B cells? Plasma-cell differentiation and the unfolded protein response. Trends Immunol. 2004; 25:17–24. https://doi.org/10.1016/j.it.2003.11.004 [PubMed]
21. Brewer JW, Hendershot LM. Building an antibody factory: a job for the unfolded protein response. Nat Immunol. 2005; 6:23–9. https://doi.org/10.1038/ni1149 [PubMed]
22. Xue X, Piao JH, Nakajima A, Sakon-Komazawa S, Kojima Y, Mori K, Yagita H, Okumura K, Harding H, Nakano H. Tumor necrosis factor alpha (TNFalpha) induces the unfolded protein response (UPR) in a reactive oxygen species (ROS)-dependent fashion, and the UPR counteracts ROS accumulation by TNFalpha. J Biol Chem. 2005; 280:33917–25. https://doi.org/10.1074/jbc.M505818200 [PubMed]
23. Nakajima S, Hiramatsu N, Hayakawa K, Saito Y, Kato H, Huang T, Yao J, Paton AW, Paton JC, Kitamura M. Selective abrogation of BiP/GRP78 blunts activation of NF-κB through the ATF6 branch of the UPR: involvement of C/EBPβ and mTOR-dependent dephosphorylation of Akt. Mol Cell Biol. 2011; 31:1710–8. https://doi.org/10.1128/MCB.00939-10 [PubMed]
24. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012; 28:882–3. https://doi.org/10.1093/bioinformatics/bts034 [PubMed]
25. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y, Kaplan S, Dahary D, Warshawsky D, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics. 2016; 54:1.30.1–33. https://doi.org/10.1002/cpbi.5 [PubMed]
26. Hollander MC, Zhan Q, Bae I, Fornace AJ Jr. Mammalian GADD34, an apoptosis- and DNA damage-inducible gene. J Biol Chem. 1997; 272:13731–7. https://doi.org/10.1074/jbc.272.21.13731 [PubMed]
27. Kidwell A, Yadav SPS, Maier B, Zollman A, Ni K, Halim A, Janosevic D, Myslinski J, Syed F, Zeng L, Waffo AB, Banno K, Xuei X, et al. Translation Rescue by Targeting Ppp1r15a through Its Upstream Open Reading Frame in Sepsis-Induced Acute Kidney Injury in a Murine Model. J Am Soc Nephrol. 2023; 34:220–40. https://doi.org/10.1681/ASN.2022060644 [PubMed]
28. Claes Z, Jonkhout M, Crespillo-Casado A, Bollen M. The antibiotic robenidine exhibits guanabenz-like cytoprotective properties by a mechanism independent of protein phosphatase PP1:PPP1R15A. J Biol Chem. 2019; 294:13478–86. https://doi.org/10.1074/jbc.RA119.008857 [PubMed]
29. Lu PD, Harding HP, Ron D. Translation reinitiation at alternative open reading frames regulates gene expression in an integrated stress response. J Cell Biol. 2004; 167:27–33. https://doi.org/10.1083/jcb.200408003 [PubMed]
30. Novoa I, Zeng H, Harding HP, Ron D. Feedback inhibition of the unfolded protein response by GADD34-mediated dephosphorylation of eIF2alpha. J Cell Biol. 2001; 153:1011–22. https://doi.org/10.1083/jcb.153.5.1011 [PubMed]
31. Novoa I, Zhang Y, Zeng H, Jungreis R, Harding HP, Ron D. Stress-induced gene expression requires programmed recovery from translational repression. EMBO J. 2003; 22:1180–7. https://doi.org/10.1093/emboj/cdg112 [PubMed]
32. Monkley S, Overed-Sayer C, Parfrey H, Rassl D, Crowther D, Escudero-Ibarz L, Davis N, Carruthers A, Berks R, Coetzee M, Kolosionek E, Karlsson M, Griffin LR, et al. Sensitization of the UPR by loss of PPP1R15A promotes fibrosis and senescence in IPF. Sci Rep. 2021; 11:21584. https://doi.org/10.1038/s41598-021-00769-7 [PubMed]
33. Patel V, Bidault G, Chambers JE, Carobbio S, Everden AJT, Garcés C, Dalton LE, Gribble FM, Vidal-Puig A, Marciniak SJ. Inactivation of Ppp1r15a minimises weight gain and insulin resistance during caloric excess in female mice. Sci Rep. 2019; 9:2903. https://doi.org/10.1038/s41598-019-39562-y [PubMed]
34. Lee JE, Morrison W, Hollien J. Hairy and enhancer of split 1 (HES1) protects cells from endoplasmic reticulum stress-induced apoptosis through repression of GADD34. J Biol Chem. 2018; 293:5947–55. https://doi.org/10.1074/jbc.RA118.002124 [PubMed]
35. Wang R, Zhang Y, Guo S, Pei S, Guo W, Wu Z, Wang H, Wang M, Li Y, Zhu Y, Meng LH, Lang J, Jin G, et al. Single-cell RNA sequencing reveals the suppressive effect of PPP1R15A inhibitor Sephin1 in antitumor immunity. iScience. 2023; 26:105954. https://doi.org/10.1016/j.isci.2023.105954 [PubMed]
36. Zhang H, Zhao Y, Yao Q, Ye Z, Mañas A, Xiang J. Ubl4A is critical for mitochondrial fusion process under nutrient deprivation stress. PLoS One. 2020; 15:e0242700. https://doi.org/10.1371/journal.pone.0242700 [PubMed]
37. Wang Q, Liu Y, Soetandyo N, Baek K, Hegde R, Ye Y. A ubiquitin ligase-associated chaperone holdase maintains polypeptides in soluble states for proteasome degradation. Mol Cell. 2011; 42:758–70. https://doi.org/10.1016/j.molcel.2011.05.010 [PubMed]
38. Krenciute G, Liu S, Yucer N, Shi Y, Ortiz P, Liu Q, Kim BJ, Odejimi AO, Leng M, Qin J, Wang Y. Nuclear BAG6-UBL4A-GET4 complex mediates DNA damage signaling and cell death. J Biol Chem. 2013; 288:20547–57. https://doi.org/10.1074/jbc.M112.443416 [PubMed]
39. Liu C, Zhou Y, Li M, Wang Y, Yang L, Yang S, Feng Y, Wang Y, Wang Y, Ren F, Li J, Dong Z, Chin YE, et al. Absence of GdX/UBL4A Protects against Inflammatory Diseases by Regulating NF-κB Signaling in Macrophages and Dendritic Cells. Theranostics. 2019; 9:1369–84. https://doi.org/10.7150/thno.32451 [PubMed]
40. Liang J, Li J, Fu Y, Ren F, Xu J, Zhou M, Li P, Feng H, Wang Y. GdX/UBL4A null mice exhibit mild kyphosis and scoliosis accompanied by dysregulation of osteoblastogenesis and chondrogenesis. Cell Biochem Funct. 2018; 36:129–36. https://doi.org/10.1002/cbf.3324 [PubMed]
41. Fu Y, Liu S, Wang Y, Ren F, Fan X, Liang J, Liu C, Li J, Ju Y, Chang Z. GdX/UBL4A-knockout mice resist collagen-induced arthritis by balancing the population of T_h1/T_h17 and regulatory T cells. FASEB J. 2019; 33:8375–85. https://doi.org/10.1096/fj.201802217RR [PubMed]
42. St-Arnaud R, Hekmatnejad B. Combinatorial control of ATF4-dependent gene transcription in osteoblasts. Ann N Y Acad Sci. 2011; 1237:11–8. https://doi.org/10.1111/j.1749-6632.2011.06197.x [PubMed]
43. Yang X, Karsenty G. ATF4, the osteoblast accumulation of which is determined post-translationally, can induce osteoblast-specific gene expression in non-osteoblastic cells. J Biol Chem. 2004; 279:47109–14. https://doi.org/10.1074/jbc.M410010200 [PubMed]
44. Yang X, Matsuda K, Bialek P, Jacquot S, Masuoka HC, Schinke T, Li L, Brancorsini S, Sassone-Corsi P, Townes TM, Hanauer A, Karsenty G. ATF4 is a substrate of RSK2 and an essential regulator of osteoblast biology; implication for Coffin-Lowry Syndrome. Cell. 2004; 117:387–98. https://doi.org/10.1016/s0092-8674(04)00344-7 [PubMed]
45. Elefteriou F, Ahn JD, Takeda S, Starbuck M, Yang X, Liu X, Kondo H, Richards WG, Bannon TW, Noda M, Clement K, Vaisse C, Karsenty G. Leptin regulation of bone resorption by the sympathetic nervous system and CART. Nature. 2005; 434:514–20. https://doi.org/10.1038/nature03398 [PubMed]
46. Yoshizawa T, Hinoi E, Jung DY, Kajimura D, Ferron M, Seo J, Graff JM, Kim JK, Karsenty G. The transcription factor ATF4 regulates glucose metabolism in mice through its expression in osteoblasts. J Clin Invest. 2009; 119:2807–17. https://doi.org/10.1172/JCI39366 [PubMed]
47. Hata K, Nishimura R, Ueda M, Ikeda F, Matsubara T, Ichida F, Hisada K, Nokubi T, Yamaguchi A, Yoneda T. A CCAAT/enhancer binding protein beta isoform, liver-enriched inhibitory protein, regulates commitment of osteoblasts and adipocytes. Mol Cell Biol. 2005; 25:1971–9. https://doi.org/10.1128/MCB.25.5.1971-1979.2005 [PubMed]
48. Tominaga H, Maeda S, Hayashi M, Takeda S, Akira S, Komiya S, Nakamura T, Akiyama H, Imamura T. CCAAT/enhancer-binding protein beta promotes osteoblast differentiation by enhancing Runx2 activity with ATF4. Mol Biol Cell. 2008; 19:5373–86. https://doi.org/10.1091/mbc.e08-03-0329 [PubMed]
49. Hirata M, Kugimiya F, Fukai A, Saito T, Yano F, Ikeda T, Mabuchi A, Sapkota BR, Akune T, Nishida N, Yoshimura N, Nakagawa T, Tokunaga K, et al. C/EBPβ and RUNX2 cooperate to degrade cartilage with MMP-13 as the target and HIF-2α as the inducer in chondrocytes. Hum Mol Genet. 2012; 21:1111–23. https://doi.org/10.1093/hmg/ddr540 [PubMed]
50. Lee AS. Glucose-regulated proteins in cancer: molecular mechanisms and therapeutic potential. Nat Rev Cancer. 2014; 14:263–76. https://doi.org/10.1038/nrc3701 [PubMed]
51. Hetz C. The unfolded protein response: controlling cell fate decisions under ER stress and beyond. Nat Rev Mol Cell Biol. 2012; 13:89–102. https://doi.org/10.1038/nrm3270 [PubMed]
52. Chen Y, Mi Y, Zhang X, Ma Q, Song Y, Zhang L, Wang D, Xing J, Hou B, Li H, Jin H, Du W, Zou Z. Dihydroartemisinin-induced unfolded protein response feedback attenuates ferroptosis via PERK/ATF4/HSPA5 pathway in glioma cells. J Exp Clin Cancer Res. 2019; 38:402. https://doi.org/10.1186/s13046-019-1413-7 [PubMed]
53. Zhu S, Zhang Q, Sun X, Zeh HJ 3rd, Lotze MT, Kang R, Tang D. HSPA5 Regulates Ferroptotic Cell Death in Cancer Cells. Cancer Res. 2017; 77:2064–77. https://doi.org/10.1158/0008-5472.CAN-16-1979 [PubMed]
54. Lv M, Cai Y, Hou W, Peng K, Xu K, Lu C, Yu W, Zhang W, Liu L. The RNA-binding protein SND1 promotes the degradation of GPX4 by destabilizing the HSPA5 mRNA and suppressing HSPA5 expression, promoting ferroptosis in osteoarthritis chondrocytes. Inflamm Res. 2022; 71:461–72. https://doi.org/10.1007/s00011-022-01547-5 [PubMed]
55. Benito MJ, Veale DJ, FitzGerald O, van den Berg WB, Bresnihan B. Synovial tissue inflammation in early and late osteoarthritis. Ann Rheum Dis. 2005; 64:1263–7. https://doi.org/10.1136/ard.2004.025270 [PubMed]
56. Yarnall BW, Chamberlain CS, Hao Z, Muir P. Proinflammatory polarization of stifle synovial macrophages in dogs with cruciate ligament rupture. Vet Surg. 2019; 48:1005–12. https://doi.org/10.1111/vsu.13261 [PubMed]
57. Liu B, Zhang M, Zhao J, Zheng M, Yang H. Imbalance of M1/M2 macrophages is linked to severity level of knee osteoarthritis. Exp Ther Med. 2018; 16:5009–14. https://doi.org/10.3892/etm.2018.6852 [PubMed]
58. Zhang H, Lin C, Zeng C, Wang Z, Wang H, Lu J, Liu X, Shao Y, Zhao C, Pan J, Xu S, Zhang Y, Xie D, et al. Synovial macrophage M1 polarisation exacerbates experimental osteoarthritis partially through R-spondin-2. Ann Rheum Dis. 2018; 77:1524–34. https://doi.org/10.1136/annrheumdis-2018-213450 [PubMed]
59. Weber A, Chan PMB, Wen C. Do immune cells lead the way in subchondral bone disturbance in osteoarthritis? Prog Biophys Mol Biol. 2019; 148:21–31. https://doi.org/10.1016/j.pbiomolbio.2017.12.004 [PubMed]
60. Zhao K, Ruan J, Nie L, Ye X, Li J. Effects of synovial macrophages in osteoarthritis. Front Immunol. 2023; 14:1164137. https://doi.org/10.3389/fimmu.2023.1164137 [PubMed]
61. Zhang H, Cai D, Bai X. Macrophages regulate the progression of osteoarthritis. Osteoarthritis Cartilage. 2020; 28:555–61. https://doi.org/10.1016/j.joca.2020.01.007 [PubMed]
62. Zhou Z, Lu J, Yang M, Cai J, Fu Q, Ma J, Zhu L. The mitochondrial unfolded protein response (UPR^mt) protects against osteoarthritis. Exp Mol Med. 2022; 54:1979–90. https://doi.org/10.1038/s12276-022-00885-y [PubMed]
63. Jang WG, Kim EJ, Kim DK, Ryoo HM, Lee KB, Kim SH, Choi HS, Koh JT. BMP2 protein regulates osteocalcin expression via Runx2-mediated Atf6 gene transcription. J Biol Chem. 2012; 287:905–15. https://doi.org/10.1074/jbc.M111.253187 [PubMed]
64. Shu B, Zhang M, Xie R, Wang M, Jin H, Hou W, Tang D, Harris SE, Mishina Y, O'Keefe RJ, Hilton MJ, Wang Y, Chen D. BMP2, but not BMP4, is crucial for chondrocyte proliferation and maturation during endochondral bone development. J Cell Sci. 2011; 124:3428–40. https://doi.org/10.1242/jcs.083659 [PubMed]
65. Hughes A, Oxford AE, Tawara K, Jorcyk CL, Oxford JT. Endoplasmic Reticulum Stress and Unfolded Protein Response in Cartilage Pathophysiology; Contributing Factors to Apoptosis and Osteoarthritis. Int J Mol Sci. 2017; 18:665. https://doi.org/10.3390/ijms18030665 [PubMed]
66. Liu W, Jiao Y, Tian C, Hasty K, Song L, Kelly DM, Li J, Chen H, Gu W, Liu S. Gene Expression Profiling Studies Using Microarray in Osteoarthritis: Genes in Common and Different Conditions. Arch Immunol Ther Exp (Warsz). 2020; 68:28. https://doi.org/10.1007/s00005-020-00592-4 [PubMed]
67. Zhang R, Yang X, Wang J, Han L, Yang A, Zhang J, Zhang D, Li B, Li Z, Xiong Y. Identification of potential biomarkers for differential diagnosis between rheumatoid arthritis and osteoarthritis via integrative genome-wide gene expression profiling analysis. Mol Med Rep. 2019; 19:30–40. https://doi.org/10.3892/mmr.2018.9677 [PubMed]
68. He P, Zhang Z, Liao W, Xu D, Fu M, Kang Y. Screening of gene signatures for rheumatoid arthritis and osteoarthritis based on bioinformatics analysis. Mol Med Rep. 2016; 14:1587–93. https://doi.org/10.3892/mmr.2016.5423 [PubMed]
69. Lin J, Wu G, Zhao Z, Huang Y, Chen J, Fu C, Ye J, Liu X. Bioinformatics analysis to identify key genes and pathways influencing synovial inflammation in osteoarthritis. Mol Med Rep. 2018; 18:5594–602. https://doi.org/10.3892/mmr.2018.9575 [PubMed]

Research Paper Volume 16, Issue 5 pp 4563—4578