Identification and validation of hub genes and potential drugs involved in osteoarthritis through bioinformatics analysis

Abstract

Purpose: Osteoarthritis (OA) is a common degenerative disease, which still lacks specific therapeutic drugs. Synovitis is one of the most important pathological process in OA. Therefore, we aim to identify and analyze the hub genes and their related networks of OA synovium with bioinformatics tools to provide theoretical basis for potential drugs.

Materials and methods: Two datasets were obtained from GEO. DEGs and hub genes of OA synovial tissue were screened through Gene Ontology (GO) annotation, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment as well as protein—protein interaction (PPI) network analysis. Subsequently, the correlation between expression of hub genes and ferroptosis or pyroptosis was analyzed. CeRNA regulatory network was constructed after predicting the upstream miRNAs and lncRNAs. The validation of hub genes was undertook through RT-qPCR and ELISA. Finally, potential drugs targeting pathways and hub genes were identified, followed by the validation of the effect of two potential drugs on OA.

Results: A total of 161 commom DEGs were obtained, of which 8 genes were finally identified as hub genes through GO and KEGG enrichment analysis as well as PPI network analysis. Eight genes related to ferroptosis and pyroptosis respectively were significantly correlated to the expression of hub genes. 24 miRNAs and 69 lncRNAs were identified to construct the ceRNA regulatory network. The validation of EGR1, JUN, MYC, FOSL1, and FOSL2 met the trend of bioinformatics analysis. Etanercept and Iguratimod reduced the secretion of MMP-13 and ADAMTS5 of fibroblast-like synoviocyte.

Conclusion: EGR1, JUN, MYC, FOSL1, and FOSL2 were identified as hub genes in the development of OA after series of bioinformatics analysis and validation. Etanercept and Iguratimod seemed to have opportunities to be novel drugs for OA.

1 Background

Osteoarthritis (OA) is a common degenerative disease in middle-aged and elderly people all over the world, of which the prevalence has gradually increased due to the aging population and the trend of the overweight (Loeser et al., 2016; Roemer et al., 2022). According to the data from World Health Organization (WHO) in 2019, about 250 million people worldwide suffered from OA (Hunter and Bierma-Zeinstra, 2019). Synovitis is one of the most important pathological manifestations in the occurrence and development of OA, running through the whole process (Sharma, 2021). Studies have shown that synovial lesions usually occur earlier than cartilage, which can be detected by MRI in several patients with small joint injury on X-ray. In addition, low-grade synovitis contributes to radiographic and pain progression (Sanchez-Lopez et al., 2022). The pathological changes of synovitis are complex and diverse, leading to the lack of specific treatment for OA. Therefore, it is necessary to explore the pathogenesis and diagnostic markers of OA from the perspective of synovium to find the therapeutic targets of OA, alleviate symptoms and promote the prognosis.

Ferroptosis is a programmed cell death mode caused by abnormal oxidation, which is regulated by glutathione peroxidase 4 (GPX4), with characteristics of iron-dependent accumulation of lipid peroxide. Kennish et al. (2014) showed that the iron level in serum of OA patients was positively correlated with the severity of OA, suggesting the existence of abnormal iron homeostasis, but its effect in synovium has not been reported. Pyroptosis is another mode of programmed cell death, which occurs when pattern recognition receptors (PRRS) induce the activation of cystine aspartic protein 1 (caspase-1) or caspase-11, which is characterized by the destruction of cell membrane and the release of cytokines. Pyroptosis has been shown to be involved in synovits, that is, Interleukin-1β (IL-1β) in the inflammatory environment derives from synovium rather than cartilage (Borgonio et al., 2014).

Bioinformatics tools have been widely used to process microarray data to determine differentially expressed genes (DEGs) and conduct various analysis. In this study, combining with bioinformatics analysis and verification, we explored and screened the hub genes of OA synovial tissue, and discussed the correlation between them and ferroptosis or pyroptosisy, then constructed the upstream ceRNA regulatory network to evaluate the reliability of these genes as the prediction or treatment targets of OA. Finally, we identified and validated potential drugs targeting pathways and hub genes through Comparative Toxicogenomics Database (CTD), Drugbank and Drugs-Genes interaction (DGI) database.

2 Materials and methods

2.1 Data collection

Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) was used as our data source (Barrett et al., 2013), where “Osteoarthritis” was entered as the keyword in the search box for detection. The selected data set included the gene expression array of osteoarthritis and normal synovial tissue of human samples. Finally, GSE55235 (GPL96 platform), GSE55457 (GPL96 platform) were determined as the data sets of this study.

2.2 Data preprocessing and identification of DEGs

The data of GSE55235, GSE55457 were downloaded with the format of MINiML. The mRNA expression data of OA and normal synovial tissue in each data set were analyzed by limma software package of R software (Yu et al., 2012). P-value was analyzed in GEO to correct false-positive results, which <0.05 and |log2 (fold change)| > 1 was defined as the threshold. DEGs were then obtained, and were visualized by volcano map and heat map. common DEGs were identified after the intersection of the two datasets.

2.3 Enrichment analysis of common DEGs and PPI network construction

Gene Ontology (GO) annotation of the common DEGs was performed by the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 2021) (Sherman et al., 2022), including biological process (BP), cellular component (CC), and molecular function (MF). p < 0.05 was determined as a significant margin for all analysis. Column chart was plotted by “http://www.bioinformatics.com.cn,” a free online platform for data analysis and visualization. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the common DEGs was carried out by R package “clusterprofiler” (Version 4). Hypergeometric test was performed to evaluate the significance of pathway–pathway association: phyper (k-1, M, N-M, n, lower. tail = F). Metascape (http://metascape.org) was used for protein - protein interaction (PPI) network analysis of common DEGs (Zhou et al., 2019). PPI analysis was performed using the following databases: STRING, BioGrid (Kondo et al., 2019), OmniPath and InWeb_IM (Li et al., 2017). Physical score >0.132 was selected to be the standard of gene screening. In addition, molecular complexity detection (MCODE) was used to identify densely connected network components and obtain the hub genes of OA (Bader and Hogue, 2003). The network was visualized using Cytoscape.

2.4 Correlation between OA hub genes and ferroptosis and pyroptosis related genes

19 genes related to ferroptosis and 21 genes related to pyroptosis contained in microarray data were selected for correlation analysis with the hub genes (Liu et al., 2020a; Wu et al., 2021). Due to the small sample size of the three groups of data related to ferroptosis, the samples were combined and normalized (Zhang et al., 2019). The data was standardized using “normalize.quantiles” function in the preprocessCore package of R, and was evaluated through the box diagram. The batch effect of data was evaluated by comparing the visual PCA diagram before and after removal. The polygenic correlation map was displayed by pheatmap package of R. Spearman’s correlation analysis was used to describe the correlation between quantitative variables without normal distribution. p < 0.05 was considered statistically significant.

2.5 Analysis of ceRNA regulatory network of hub genes

ENCORI (Li et al., 2014) (https://starbase.sysu.edu.cn) and TargetScan (McGeary et al., 2019) (https://www. targetscan. org/vert_80/) were used to predict miRNAs which regulated hub genes. The first three reliable miRNAs were selected after the intersection of the prediction. LncBase database v3.0 (Paraskevopoulou et al., 2016) (https://diana.e-ce.uth.gr/lncbasev3/interactions) was used to predict lncRNAs which regulated miRNAs above. The first three reliable lncRNAs were selected for each miRNA. Finally, the ceRNA regulatory network of hub genes was constructed.

2.6 Extraction of human fibroblast-like synoviocyte (FLS)

The study was designed according to the Declaration of Helsinki, and was approved by Ethic Committee of the Second Affiliated Hospital of Harbin Medical University (KY 2021-256). Informed consent was obtained from each donor. Synovium of 3 OA patients (age 54–70 years, Kellgren-Lawrence grade 4) that underwent total joint arthroplasty (TKA) and 3 patients (age 56–68 years) that underwent meniscectomy without OA were obtained at the time after surgery. All patients were confirmed without Rheumatoid Arthritis, acute trauma, tumor or infection of knee joint. Briefly, synovium was cut into pieces at the final size about 0.5 mm*0.5 mm, and put into 0.1% type I collagenase (Biosharp, China, BS163). α-MEM medium (Cytiva, United States, SH30265.01) was added with 10% foetal bovine serum (ExCell Bio, China, FSD500) after 2 h. Primary cells could be seen climbing out after about 3–5 days. FLS of passage 6-8 (P6-8) were used in this study. 10 ng/ml of IL-1β (PEPROTECH, United States, 200-01B) were used to stimulate FLS of OA groups for 48 h in order to imitate the environment of OA, while complete medium was added into the FLS of control group. 10 μg/ml of Etanercept and Iguratimod (MedChem Express, China, HY-108847, HY-17009) were added along with IL-1β to FLS for the following test.

2.7 Screening of potential drugs for OA

CTD database (version 16766M) was used to search drugs targeting the KEGG pathways above. Drugs targeting hub genes were searched and compared using Drugbank (version 5.1.8) database and DGI database (version 4.2.0).

2.8 RT q-PCR

Trizol (Beyotime, China, R0016) was used to extract total RNA from FLS. After determining the concentration, 2 µg of total RNA was used to synthesize cDNA through cDNA synthesis kit. SYBR Green (ES Science, China, QP002) was used for qRT-PCR according to the instructions. The primer sequence of 8 hub genes and GAPDH were listed in Table 1. The mRNA level of a specific gene was calculated as 2-^ΔΔ Ct and normalized to GAPDH.

2.9 Elisa analysis

The supernatant was collected and centrifuged at 1,000 × g at 4°C for 10 min, then was added to 96-well plates coverd by the antibody of each protein (Elabscience, United States). After incubated for 90 min at 37°C, Biotinylated Detection Ab was added to each well. After 60 min, the plates were washed for 3 times. HRP Conjugate was then added to each well and incubated for 30 min at 37°C in the dark. 90 μL of Substrate Reagent and 50 μL Stop Solution were then added. The optical density (OD value) of each well were determined at once by microplate reader at 450 nm.

2.10 Detection of intracellular ROS

FLSs of each group were seeded at a density about 2000 cells per well in 96-well plates. The medium was replaced by H2DCFDA (Biosharp, China, BL714A) working solution with the concentration of 10 μM. FLSs were then incubated at 37°C for 2 h in the dark and washed twice to fully remove the probes outside the cells. OD value of each well was determined at once by microplate reader at 525 nm.

2.11 Detection of lipid peroxidation

FLSs of each group were broken by ultrasonic cell crusher and centrifuged at 3,500 xg for 10 min at 4°C. MDA standard solution (Solarbio, China, BC0025) and the supernatant to be tested with chromogenic agent were boiled, centrifuged, and added to a 96-well plate respectively. The absorbance was measured at 450 nm, 520 nm, and 600 nm for detection by microplate reader. Then the concentration of MDA was calculated according to the instruction.

2.12 Cell viability assay

FLSs were seeded in 96-well plates at a density at about 2000 cells per well. After adding 10 ul of CCK-8 reagent (Beyotime, China, C0037) at 24 h or 48 h respectively, the plates were put into the incubator for 30 min away from light. The absorbance of each well was measured at 450 nm by microplate reader.

2.13 Statistical analysis

All the data of validation was presented as the means ± standard deviation (SD). Statistical analysis was performed using GraphPad Prism version 6.02. Differences in numerical data between two groups were determined by t-test, while four groups were determined by One-way ANOVA followed by a Bonferroni post-hoc test. p < 0.05 was defined as statistically significant.

3 Results

3.1 Data collection and identification of DEGs

The flow chart of the study design was showed in Figure 1. A total of 20 OA synovium samples and 20 normal synovium samples were obtained from two GEO database, including 10 OA samples and 10 normal samples from GSE55235 dataset and 10 OA samples and 10 normal samples from GSE55457 dataset. Relevant clinical data were not provided in the two datasets above. After screening the differential expression of mRNA between OA and normal synovium samples, 595 up-regulated DEGs and 346 down-regulated DEGs from GSE55235 as well as 175 up-regulated DEGs and 385 down-regulated DEGs from GSE55457 were obtained (Figures 2A, B). After the intersection, a total of 161 common DEGs were obtained (Figure 2C; Table 2).

FIGURE 1

FIGURE 2

3.2 Enrichment analysis and protein-protein interaction analysis of DEGs in OA

In order to further explore the function of DEGs in OA, GO annotation and KEGG enrichment analysis were performed subsequently. As the results, 156 GO terms with p < 0.05 were screened out (S1). Top 5 BP, 5 CC and 4 MF terms were showed in Figure 3A. It showed that the most significant BP of 161 commom DEGs was “positive regulation of transcription from RNA polymerase II promoter”, the most significant CC was “extracellular space”, and the most significant MF was “transcriptional activator activity, RNA polymerase II transcription regulatory region sequence-specific binding”. As the result of KEGG enrichment analysis, 38 pathways with p < 0.05 were found (S2), and top 15 pathways were showed in Figure 3B. It showed that these genes were involved in pathways including tumor necrosis factor (TNF) signaling pathway, osteoclast differentiation, IL-17 signaling pathway, nuclear factor kappa-B (NF-κB) signaling pathway and so on. In addition, as the result of PPI analysis combined with MCODE component algorithm, top 3 modules with highest MCODE score were showed in Figures 4A, B. Finally, 8 hub genes in the biggest module were identified as hub genes, which were ATF3, EGR1, FOSB, FOSL1, FOSL2, JUN, JUNB, and MYC (Table 3).

FIGURE 3

FIGURE 4

3.3 Correlation analysis between OA hub genes and ferroptosis related genes

We combined and normalized the data of 2 groups as described above, then evaluated the data standardization and batch effect through the box diagram (Figure 5A) and visual PCA diagram (Figure 5B), respectively. 19 ferroptosis related genes and 21 pyroptosis related genes were selected for correlation analysis and displayed by heat map. It showed that in the analysis of ferroptosis, ALOX15, CISD1, SAT1, and TFRC were positively correlated to hub genes, while ATP5MC3, GPX4, HSPB1, and MT1G were negatively correlated (Figure 5C). In the analysis of pyroptosis, Caspase-6, ELANE, GSDMB, IL-6 and NLRP1 were positively correlated to hub genes, while GPX4, NOD1 and PYCARD were negatively correlated (Figure 5D).

FIGURE 5

3.4 Analysis of ceRNA regulatory network of OA hub genes

A total of 24 reliable miRNAs that regulated mRNA of 8 hub genes were obtained, after the intersection of the prediction, using ENCORI and TargetScan database. Then, a total of 69 reliable lncRNAs were obtained after the analysis under LncBase v3.0 database and selection, which regulated miRNAs above. Finally, we constructed the ceRNA regulatory network of hub genes (Figure 6).

FIGURE 6

3.5 Validation of hub genes

qRT-PCR and ELISA were used to evaluate the mRNA and protein expression level of hub genes in FLS of control group and OA group. The results showed that the mRNA expression levels of JUN, MYC, EGR1, FOSL1, and FOSL2 of FLS in OA group was significantly lower than that in normal FLS (p < 0.001, p < 0.01, p < 0.01, p < 0.05, p < 0.05, respectively), meeting the trend of bioinformatics analysis (Figures 7A–E). The expression of JUNB and FOSB had no significant difference with control group (Figures 7F, G), while the expression of ATF3 was higher than that in normal FLS (p < 0.001, Figure 7H), which was opposite to bioinformatics analysis. The result of ELISA showed that the expression of EGR1, JUN, MYC were all significantly lower in the supernatant of FLS in OA group, same to the result of qRT-PCR (p < 0.0001, Figure 7I).

FIGURE 7

3.6 Identification of potential drugs for OA

Drugs targeting 8 hub genes and top 4 KEGG pathways were identified by CTD, DrugBank and DGI database. As the result, 85 drugs targeting 4 of hub genes were obtained, among which Quinapril targeted MYC and JUN at the same time (Figure 8A). 285 drugs were found targeting KEGG pathways, among which Elsubrutinib, Lenalidomide, Lenercept, Mifamurtide, Opinercept, Rebimastat, Tasonermin, Thalidomide and Tibulizumab targeted two pathways. Etanercept and Iguratimod targeted three pathways, which were selected for following validation (Table 4; Figure 8B).

FIGURE 8

3.7 Etanercept and Iguratimod protect FLS of OA from inflammation and cartilage degeneration

To verify the effect of Etanercept and Iguratimod on FLS, we co-treated FLS with IL-1β. By the test of cell viability, we found that Etanercept can significantly inhibit the growth rate of FLS compared with OA group (p < 0.01, Figure 9A). The results of intracellular ROS and MDA detection showed that the levels of these two species in Etanercept and Iguratimod groups were significantly lower than those of OA group (p < 0.01, p < 0.05, Figures 9B, C). At the same time, the mRNA expression of EGR1, JUN and MYC in Etanercept and Iguratimod group were significantly higher than those in OA group (p < 0.001, p < 0.01 for EGR1, p < 0.001, p < 0.05 for JUN and MYC, Figures 9D–F). Moreover, the level of MMP-13 and ADAMTS5 in the supernatant of FLS was significantly lower in Etanercept and Iguratimod group than those in OA group (p < 0.01 for MMP-13, p < 0.001, p < 0.01 for ADAMTS5, Figures 9G, H).

FIGURE 9

4 Discussion

As one of the major chronic diseases endangering the middle-aged and elderly, the affected population of OA is also showing a trend of younger age (Mahmoudian et al., 2021). The compliance of non-drug treatment such as kinesitherapy and physiotherapy is usually inexact. The operations for KOA mainly includes arthroscopic debridement, high tibial osteotomy (HTO), unicondylar knee arthroplasty (UKA) and TKA (Zhang et al., 2020), which still have the problems of high risk and cost. The gathered evidence suggests that mononuclear infiltration and over expression of inflammatory mediators in synovium are seen in early OA and predate radiographic damage in OA (Sokolove and Lepus, 2013). Synovium and synovial fluid are the main contributors to inflammation that can secrete key cytokines, most of which are the main regulators of matrix metalloproteinases (MMP), resulting in the loss of cartilage (Kulkarni et al., 2021). Therefore, in-depth study of the molecular mechanism of synovitis is of great significance to the prevention and treatment for OA.

In this study, we screened out 8 hub genes in synovim that may be involved in the progression of OA through a series of bioinformatics methods. EGR1 is potentially involved in postnatal bone biology and implicated in the regulation of osteoblastic cell growth and differentiation. Studies have shown that when EGR1 was inhibited, the bone mass of the limbs decreased in mice, with lower bone volume fraction and mineral density, same to our study (Reumann et al., 2011). However, over-expression of EGR1 regulated the expression of KLF5 and β-catenin signaling pathway, leading to the acceleration of cartilage hypertrophy and degeneration (Sun et al., 2019). MYC encodes transcription factors to regulate transcriptional activity and cell proliferation, growth, and had also been found to regulate OA process through multiple pathways. Wu et al. (O'Donnell et al., 2005) found that the effects of miR-24 on OA chondrocytes may be achieved by targeting MYC and further regulating the MAPK signaling pathway. The silence of c-MYC could promote proliferation of rat chondrocyte (Zou et al., 2018). ATF3 deficiency in chondrocytes had been reported to alleviates OA development (Iezaki et al., 2016). This effect may be achieved by ATF3 directly affecting the expression of MMP-13 thus reducing cartilage loss (Chan et al., 2017). However, There is little study on the mechanism of ATF3 in OA synoviocyte, except for the bioinformatics analysis.

The rest five hub genes, which were FOSB, FOSL1, FOSL2, JUN and JUNB, belonged to the AP-1 family. AP-1 was thought to be critically involved in the pathogenesis of arthritis due to the binding activity for its cognate recognition sites in the promoters of inflammatory cytokines and matrix-degrading enzymes (Huber et al., 2019). Study showed that mitochondrial dysfunction triggered a catabolic response in chondrocytes via activation of the JNK/AP1 pathway (Ansari et al., 2020). On the other hand, JUNB-FBXO21-ERK axis promoted cartilage degeneration in OA by inhibiting autophagy (Lin et al., 2021).

As independent forms of programmed death, ferroptosis and pyroptosis have been shown to play essential roles in the pathological processes of tumor, Alzheimer’s disease, cerebral hemorrhage, ischemia-reperfusion injury, OA and so on (Kenny et al., 2019; Liang et al., 2019; Liu et al., 2020b). In this study, we selected several ferroptosis and pyroptosis related genes. Although these genes were not differentially expressed in the two datasets, analyzing their relationship with hub genes will be still beneficial to explore the mechanism of FLS participating in OA. We found that FOSL1 and JUNB were associated with 3 of the 19 ferroptosis related genes. Overexpression of c-JUN inhibits ferroptosis induced by erastin in Schwann cells to promotes the rehabilitation of facial nerve function (Gao et al., 2022). Vertically, the strongest two genes related to hub genes were CISD1 and HSPB1. CISD1, an iron-containing outer mitochondrial membrane protein, inhibits ferroptosis by protection against mitochondrial lipid peroxidation (Yuan et al., 2016). The phosphorylation of HSPB1 at the Ser-15 site induced by Erastin is the key of the protective response to ferroptosis stress (Wang et al., 2022). In addition, over-expression of HSPB1 attenuates ferroptosis in rats through promoting G6PD expression (Dai and Hu, 2022). Therefore, we believe that CISD1 and HSPB1 are worth studying in the future to explore the mechanism of ferroptosis in OA synovium. In addition, we found that JUN was related to the most pyroptosis related genes. Study had showed that the Jun N-terminal kinases (JNK) pathway is largely upstream of the NLRP3 inflammasome, which exerts a crucial regulatory impact on microglia pyroptosis and inflammatory responses (He et al., 2021). However, the relationship between JUN and pyroptosis of FLS still needs further study. As a key regulator of innate immunity and inflammasome activation, caspase-6 promotes the activation of pyroptosis (Zheng et al., 2020).

Studies have found the significance of miRNAs and lncRNAs on OA. Xu et al. (2021) showed that SNHG7 ameliorated the development of OA by suppressing apoptosis through miR-214-5p-PPARGC1B-PPARγ axis. Wang et al. (2019) found that lncRNA FOXD2 Adjacent Opposite Strand RNA 1 (FOXD2-AS1) served as the protector for OA patients by inducing chondrocyte proliferation. The downregulation of lncRNA LOC101928134, which acts as a promoter of OA, can block the process of OA (Yang et al., 2019). Here, we built a ceRNA network after predicting the upstream miRNAs and lncRNAs, in order to provide more possible early biomarkers for diagnosis, or targeting drugs for OA patients.

Among the 11 drugs targeting two or more than two KEGG pathways, 10 acted on NF-κB signaling pathway, 9 acted on TNF signaling pathway, 3 acted on Osteoclast differentiation pathway, and 2 acted on IL-17 signaling pathway, suggesting that NF-κB was still the key point for OA treatment. Moreover, Etanercept and Iguratimod related to three pathways. Etanercept, a soluble fusion protein that binds TNF-α, has been proved to be an effective choice targeting several inflammatory diseases, especially rheumatoid arthritis (Graudal et al., 2015). Recent study showed that Etanercept lead to the repairment of cartilage with a special scaffold (Campos et al., 2022), and the reduction of MMP-3 (Kroon et al., 2020), suggesting the potential effect on OA. The effect of Iguratimod on OA still needs further study. Our study proved that both Etanercept and Iguratimod could reduce level of intracellular ROS and MDA of FLS as well as level of MMP-13 and ADAMTS5 in the supernatant. In addition, Etanercept could significantly alleviate the abnormal proliferation of FLS under OA environment. Drugs targeted to the hub genes can be divided into antineoplastic, antipsychotic, anti-inflammatory, anti-hyperglycemia, anti-hyperlipidemia, etc. Their roles in OA still needs to be explored.

Our study still had some limitations. First, we obtained gene arrays in synovial tissues of OA patients from GEO database. However, due to the lack of information of patients, it was hard to correlate the DEGs and hub genes obtained in this study with the severity or grade of OA. Second, the sample size of each group was not so large, which was needed for further research.

In summary, our study aimed to identify key genes involved in the pathophysiology of OA. 161 common DEGs and 8 hub genes were screened through GO analysis, KEGG pathway analysis, PPI network construction as well as MCODE, which may become potential targeting clinical diagnosis and treatment of OA. Furthermore, 5 of 8 genes met the similar expression trend with our result through the validation. Subsequently, we discussed the relationship between hub genes and key genes of ferroptosis and pyroptosis. Moreover, miRNAs and lncRNAs were identified to construct the ceRNA regulatory network of hub genes. Finally, we found that Etanercept and Iguratimod, as top two of potential drugs targeting KEGG pathways, had the protective effect on FLS in the OA environment.

5 Conclusion

EGR1, JUN, MYC, FOSL1, and FOSL2 were identified and validated as hub genes in the development of OA after series of bioinformatics analysis. They may have effect on OA development through different kinds of pathways, including the process of ferroptosis and pyroptosis. Etanercept and Iguratimod seemed to have stronger opportunities to be novel drugs for OA.

References

• 1

  AnsariM. Y.AhmadN.VoletiS.WaseS. J.NovakK.HaqqiT. M.et al (2020). Mitochondrial dysfunction triggers a catabolic response in chondrocytes via ROS-mediated activation of the JNK/AP1 pathway. J. Of Cell Sci.133 (22), jcs247353. 10.1242/jcs.247353

  • CrossRef
  • Google Scholar

• 2

  BaderG. D.HogueC. V. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinforma.4, 2. 10.1186/1471-2105-4-2

  • CrossRef
  • Google Scholar

• 3

  BarrettT.WilhiteS. E.LedouxP.EvangelistaC.KimI. F.TomashevskyM.et al (2013). NCBI GEO: Archive for functional genomics data sets-update. Nucleic Acids Res.41, D991–D995. 10.1093/nar/gks1193

  • CrossRef
  • Google Scholar

• 4

  BorgonioC.González-HuertaN. C.Romero-CórdobaS.Hidalgo-MirandaA.Miranda-DuarteA. (2014). Altered expression of circulating microRNA in plasma of patients with primary osteoarthritis and in silico analysis of their pathways. PloS one9 (6), e97690. 10.1371/journal.pone.0097690

  • CrossRef
  • Google Scholar

• 5

  CamposY.FuentesG.AlmirallA.QueI.SchomannT.ChungC. K.et al (2022). The incorporation of Etanercept into a porous tri-layer scaffold for restoring and repairing cartilage tissue. Pharmaceutics14 (2), 282. 10.3390/pharmaceutics14020282

  • CrossRef
  • Google Scholar

• 6

  ChanC.MacdonaldC.LitherlandG.WilkinsonD. J.SkeltonA.Europe-FinnerG. N.et al (2017). Cytokine-induced MMP13 expression in human chondrocytes is dependent on activating transcription factor 3 (ATF3) regulation. J. Of Biol. Chem.292 (5), 1625–1636. 10.1074/jbc.M116.756601

  • CrossRef
  • Google Scholar

• 7

  DaiY.HuL. (2022). HSPB1 overexpression improves hypoxic-ischemic brain damage by attenuating ferroptosis in rats through promoting G6PD expression. J. Of Neurophysiology128 (6), 1507–1517. 10.1152/jn.00306.2022

  • CrossRef
  • Google Scholar

• 8

  DhillonA. S.TulchinskyE. (2015). FRA-1 as a driver of tumour heterogeneity: A nexus between oncogenes and embryonic signalling pathways in cancer. Oncogene34, 4421–4428. 10.1038/onc.2014.374

  • CrossRef
  • Google Scholar

• 9

  GaoD.HuangY.SunX.YangJ.ChenJ.HeJ. (2022). Overexpression of c-Jun inhibits erastin-induced ferroptosis in Schwann cells and promotes repair of facial nerve function. J. Of Cell. And Mol. Med.26 (8), 2191–2204. 10.1111/jcmm.17241

  • CrossRef
  • Google Scholar

• 10

  GraudalN.JürgensG.CombeB. (2015). Etanercept tapering in rheumatoid arthritis. N. Engl. J. Of Med.372 (5), 489–490. 10.1056/NEJMc1414787

  • CrossRef
  • Google Scholar

• 11

  GuoJ.ShenS.LiuX.RuanX.ZhengJ.LiuY.et al (2019). Role of linc00174/miR-138-5p (miR-150-5p)/FOSL2 feedback loop on regulating the blood-tumor barrier permeability. Mol. Ther. Nucleic Acids.18, 1072–1090. 10.1016/j.omtn.2019.10.031

  • CrossRef
  • Google Scholar

• 12

  HeJ.LiuT.LiY.MiX.HanD.YangN.et al (2021). JNK inhibition alleviates delayed neurocognitive recovery after surgery by limiting microglia pyroptosis. Int. Immunopharmacol.99, 107962. 10.1016/j.intimp.2021.107962

  • CrossRef
  • Google Scholar

• 13

  HuberR.KirstenH.NäkkiA.PohlersD.ThudeH.EidnerT.et al (2019). Association of human FOS promoter variants with the occurrence of knee-osteoarthritis in a case control association study. Int. J. Mol. Sci.20 (6), 1382. 10.3390/ijms20061382

  • CrossRef
  • Google Scholar

• 14

  HunterD. J.Bierma-ZeinstraS. (2019). Osteoarthr. Lancet393 (10182), 1745–1759. 10.1016/S0140-6736(19)30417-9

  • CrossRef
  • Google Scholar

• 15

  IezakiT.OzakiK.FukasawaK.InoueM.KitajimaS.MunetaT.et al (2016). ATF3 deficiency in chondrocytes alleviates osteoarthritis development. J. Of Pathology239 (4), 426–437. 10.1002/path.4739

  • CrossRef
  • Google Scholar

• 16

  JamesC. G.WoodsA.UnderhillT. M.BeierF. (2006). The transcription factor ATF3 is upregulated during chondrocyte differentiation and represses cyclin D1 and A gene transcription. BMC Mol. Biol.7, 30. 10.1186/1471-2199-7-30

  • CrossRef
  • Google Scholar

• 17

  KennishL.AtturM.OhC.KrasnokutskyS.SamuelsJ.GreenbergJ. D.et al (2014). Age-dependent ferritin elevations and HFE C282Y mutation as risk factors for symptomatic knee osteoarthritis in males: A longitudinal cohort study. BMC Musculoskelet. Disord.15, 8. 10.1186/1471-2474-15-8

  • CrossRef
  • Google Scholar

• 18

  KennyE. M.FidanE.YangQ.AnthonymuthuT. S.NewL. A.MeyerE. A.et al (2019). Ferroptosis contributes to neuronal death and functional outcome after traumatic brain injury. Crit. Care Med.47 (3), 410–418. 10.1097/CCM.0000000000003555

  • CrossRef
  • Google Scholar

• 19

  KondoK.OhigashiI.TakahamaY. (2019). Thymus machinery for T-cell selection. Int. Immunol.31 (3), 119–125. 10.1093/intimm/dxy081

  • CrossRef
  • Google Scholar

• 20

  KroonF.Bay-JensenA. C.WittoekR.VerbruGGenG.SmolenJ. S.KloppenburgM.et al (2020). Etanercept therapy leads to reductions in matrix metalloproteinase-3 in patients with erosive hand osteoarthritis. Scand. J. Of Rheumatology49 (2), 167–168. 10.1080/03009742.2019.1657493

  • CrossRef
  • Google Scholar

• 21

  KulkarniP.MartsonA.VidyaR.ChitnavisS.HarsulkarA. (2021). Pathophysiological landscape of osteoarthritis. Adv. Clin. Chem.100, 37–90. 10.1016/bs.acc.2020.04.002

  • CrossRef
  • Google Scholar

• 22

  LiJ. H.LiuS.ZhouH.QuL. H.YangJ. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res.42, D92–D97. 10.1093/nar/gkt1248

  • CrossRef
  • Google Scholar

• 23

  LiT.WernerssonR.HansenR. B.HornH.MercerJ.SlodkowiczG.et al (2017). A scored human protein-protein interaction network to catalyze genomic interpretation. Nat. METHODS14 (1), 61–64. 10.1038/nmeth.4083

  • CrossRef
  • Google Scholar

• 24

  LiangC.ZhangX.YangM.DongX. (2019). Recent progress in ferroptosis inducers for cancer therapy. Adv. Mater.31 (51), e1904197. 10.1002/adma.201904197

  • CrossRef
  • Google Scholar

• 25

  LichtA. H.PeinO. T.FlorinL.HartensteinB.ReuterH.ArnoldB.et al (2006). JunB is required for endothelial cell morphogenesis by regulating core-binding factor beta. J. Of Cell Biol.175 (6), 981–991. 10.1083/jcb.200605149

  • CrossRef
  • Google Scholar

• 26

  LinZ.MiaoJ.ZhangT.HeM.WangZ.FengX.et al (2021). JUNB-FBXO21-ERK axis promotes cartilage degeneration in osteoarthritis by inhibiting autophagy. Aging Cell20 (2), e13306. 10.1111/acel.13306

  • CrossRef
  • Google Scholar

• 27

  LiuP.FengY.LiH.ChenX.WangG.XuS.et al (2020). Ferrostatin-1 alleviates lipopolysaccharide-induced acute lung injury via inhibiting ferroptosis. Cell Mol. Biol. Lett.25 (10), 10. 10.1186/s11658-020-00205-0

  • CrossRef
  • Google Scholar

• 28

  LiuZ.ZhaoQ.ZuoZ. X.YuanS. Q.YuK.ZhangQ.et al (2020). Systematic analysis of the aberrances and functional implications of ferroptosis in cancer. iScience23 (7), 101302. 10.1016/j.isci.2020.101302

  • CrossRef
  • Google Scholar

• 29

  LoeserR. F.CollinsJ. A.DiekmanB. O. (2016). Ageing and the pathogenesis of osteoarthritis. Nat. Rev. Rheumatol.12 (7), 412–420. 10.1038/nrrheum.2016.65

  • CrossRef
  • Google Scholar

• 30

  LuH.HouG.ZhangY.DaiY.ZhaoH. (2014). c-Jun transactivates Puma gene expression to promote osteoarthritis. Mol. Med. Rep.9 (5), 1606–1612. 10.3892/mmr.2014.1981

  • CrossRef
  • Google Scholar

• 31

  MahmoudianA.LohmanderL. S.MobasheriA.EnglundM.LuytenF. P. (2021). Early-stage symptomatic osteoarthritis of the knee - time for action. Nat. Rev. Rheumatol.17 (10), 621–632. 10.1038/s41584-021-00673-4

  • CrossRef
  • Google Scholar

• 32

  McGearyS. E.LinK. S.ShiC. Y.PhamT. M.BisariaN.KelleyG. M.et al (2019). The biochemical basis of microRNA targeting efficacy. Science366 (6472), eaav1741. 10.1126/science.aav1741

  • CrossRef
  • Google Scholar

• 33

  O'DonnellK. A.WentzelE. A.ZellerK. I.DangC. V.MendellJ. T. (2005). c-Myc-regulated microRNAs modulate E2F1 expression. Nature435 (7043), 839–843. 10.1038/nature03677

  • CrossRef
  • Google Scholar

• 34

  Papoudou-BaiA.HatzimichaelE.BarboutiA.KanavarosP. (2017). Expression patterns of the activator protein-1 (AP-1) family members in lymphoid neoplasms. Clin. And Exp. Med.17 (3), 291–304. 10.1007/s10238-016-0436-z

  • CrossRef
  • Google Scholar

• 35

  ParaskevopoulouM. D.VlachosI. S.KaragkouniD.GeorgakilasG.KanellosI.VergoulisT.et al (2016). DIANA-LncBase v2: Indexing microRNA targets on non-coding transcripts. Nucleic Acids Res.44 (D1), D231–D238. 10.1093/nar/gkv1270

  • CrossRef
  • Google Scholar

• 36

  ReumannM. K.StrachnaO.YagermanS.TorrecillaD.KimJ.DotyS. B.et al (2011). Loss of transcription factor early growth response gene 1 results in impaired endochondral bone repair. BONE49 (4), 743–752. 10.1016/j.bone.2011.06.023

  • CrossRef
  • Google Scholar

• 37

  RoemerF. W.GuermaziA.HannonM. J.FujiiT.OmoumiP.HunterD. J.et al (2022). Presence of magnetic resonance imaging- defined inflammation particularly in overweight and obese women increases risk of radiographic knee osteoarthritis: The POMA study. Arthritis Care & Res.74 (8), 1391–1398. 10.1002/acr.24568

  • CrossRef
  • Google Scholar

• 38

  Sanchez-LopezE.CorasR.TorresA.LaneN. E.GumaM. (2022). Synovial inflammation in osteoarthritis progression. Nat. Rev. Rheumatol.18 (5), 258–275. 10.1038/s41584-022-00749-9

  • CrossRef
  • Google Scholar

• 39

  SharmaL. (2021). Osteoarthritis of the knee. N. Engl. J. Of Med.384 (1), 51–59. 10.1056/NEJMcp1903768

  • CrossRef
  • Google Scholar

• 40

  ShermanB. T.HaoM.QiuJ.JiaoX.BaselerM. W.LaneH. C.et al (2022). David: A web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res.50, W216–W221. 10.1093/nar/gkac194

  • CrossRef
  • Google Scholar

• 41

  SokoloveJ.LepusC. M. (2013). Role of inflammation in the pathogenesis of osteoarthritis: Latest findings and interpretations. Ther. Adv. Musculoskelet. Dis.5 (2), 77–94. 10.1177/1759720X12467868

  • CrossRef
  • Google Scholar

• 42

  SunX.HuangH.PanX.LiS.XieZ.MaY.et al (2019). EGR1 promotes the cartilage degeneration and hypertrophy by activating the Krüppel-like factor 5 and β-catenin signaling. Biochimica Biophysica Acta-Molecular Basis Of Dis.1865 (9), 2490–2503. 10.1016/j.bbadis.2019.06.010

  • CrossRef
  • Google Scholar

• 43

  TanelA.FonsecaS. G.Yassine-DiabB.BordiR.ZeidanJ.ShiY.et al (2009). Cellular and molecular mechanisms of memory T-cell survival. Expert Rev. Vaccines8 (3), 299–312. 10.1586/14760584.8.3.299

  • CrossRef
  • Google Scholar

• 44

  WangL.WuS.HeH.AiK.XuR.ZhangL.et al (2022). CircRNA-ST6GALNAC6 increases the sensitivity of bladder cancer cells to erastin-induced ferroptosis by regulating the HSPB1/P38 axis. Lab. Investig.102 (12), 1323–1334. 10.1038/s41374-022-00826-3

  • CrossRef
  • Google Scholar

• 45

  WangY.CaoL.WangQ.HuangJ.XuS. (2019). LncRNA FOXD2-AS1 induces chondrocyte proliferation through sponging miR-27a-3p in osteoarthritis. Artif. Cells Nanomedicine Biotechnol.47 (1), 1241–1247. 10.1080/21691401.2019.1596940

  • CrossRef
  • Google Scholar

• 46

  WuJ. B.ZhuY.LuoM. M.LiL. (2021). Comprehensive analysis of pyroptosis-related genes and tumor microenvironment infiltration characterization in breast cancer. Front. Immunol.12, 748221. 10.3389/fimmu.2021.748221

  • CrossRef
  • Google Scholar

• 47

  XuJ.PeiY.LuJ.LiangX.LiY.WangJ.et al (2021). LncRNA SNHG7 alleviates IL-1β-induced osteoarthritis by inhibiting miR-214-5p-mediated PPARGC1B signaling pathways. Int. Immunopharmacol.90, 107150. 10.1016/j.intimp.2020.107150

  • CrossRef
  • Google Scholar

• 48

  YangD. W.ZhangX.QianG. B.JiangM. J.WangP.WangK. Z.et al (2019). Downregulation of long noncoding RNA LOC101928134 inhibits the synovial hyperplasia and cartilage destruction of osteoarthritis rats through the activation of the Janus kinase/signal transducers and activators of transcription signaling pathway by upregulating IFNA1. J. Of Cell. Physiology234 (7), 10523–10534. 10.1002/jcp.27730

  • CrossRef
  • Google Scholar

• 49

  YuG.WangL.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics-A J. Of Integr. Biol.16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 50

  YuanH.LiX.ZhangX.KangR.TangD. (2016). CISD1 inhibits ferroptosis by protection against mitochondrial lipid peroxidation. Biochem. And Biophysical Res. Commun.478 (2), 838–844. 10.1016/j.bbrc.2016.08.034

  • CrossRef
  • Google Scholar

• 51

  ZhangX.ZhangW.JiangY.LiuK.RanL.SongF. (2019). Identification of functional lncRNAs in gastric cancer by integrative analysis of GEO and TCGA data. J. Of Cell. Biochem.120 (10), 17898–17911. 10.1002/jcb.29058

  • CrossRef
  • Google Scholar

• 52

  ZhangZ.HuangC.JiangQ.ZhengY.LiuY.LiuS.et al (2020). Guidelines for the diagnosis and treatment of osteoarthritis in China (2019 edition). Ann. Transl. Med.8 (19), 1213. 10.21037/atm-20-4665

  • CrossRef
  • Google Scholar

• 53

  ZhengM.KarkiR.VogelP.KannegantiT. D. (2020). Caspase-6 is a key regulator of innate immunity, inflammasome activation, and host defense. Cell181 (3), 674–687. 10.1016/j.cell.2020.03.040

  • CrossRef
  • Google Scholar

• 54

  ZhouY.ZhouB.PacheL.ChangM.KhodabakhshiA. H.TanaseichukO.et al (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun.10 (1), 1523. 10.1038/s41467-019-09234-6

  • CrossRef
  • Google Scholar

• 55

  ZouJ.LiX.ShiZ.XueJ. F. (2018). Effects of C-myc gene silencing on interleukin-1β-induced rat chondrocyte cell proliferation, apoptosis and cytokine expression. J. Of Bone And Mineral Metabolism36 (3), 286–296. 10.1007/s00774-017-0845-4

  • CrossRef
  • Google Scholar