Comprehensive analysis of pivotal biomarkers, immune cell infiltration and therapeutic drugs for steroid-induced osteonecrosis of the femoral head

ABSTRACT

Steroid-induced osteonecrosis of the femoral head (SONFH) is a progressive disease that leads to an increased disability rate. This study aimed to ascertain biomarkers, infiltrating immune cells, and therapeutic drugs for SONFH. The gene expression profile of the GSE123568 dataset was downloaded from the Gene Expression Omnibus (GEO) database. The differentially expressed genes (DEGs) were identified using the NetworkAnalyst platform. Functional enrichment, protein-protein interaction network (PPI), and module analyses were performed using Metascape tools. An immune cell abundance identifier was used to explore immune cell infiltration. Furthermore, hub genes were identified based on maximal clique centrality (MCC) evaluation using cytoHubba application and confirmed by qRT-PCR using clinical samples. Finally, the L1000 platform was used to determine potential drugs for SONFH treatment. The SONFH mouse model was used to determine the therapeutic effects of aspirin. In total, 429 DEGs were identified in SONFH samples. Functional enrichment analysis showed that these DEGs were enriched in myeloid leukocyte activation and osteoclast differentiation processes. A set of nine immune cell types was confirmed to be markedly different between the SONFH and control samples. All 10 hub genes were significantly highly expressed in the serum of SONFH patients, as shown by qRT-PCR. Finally, the therapeutic effect of aspirin on SONFH was examined in animal experiments. Taken together, our data revealed the hub genes and infiltrating immune cells in SONFH, and we also screened potential drugs for use in SONFH treatment.

KEYWORDS:

• SONFH
• immune cell infiltration
• biomarker
• therapeutics

Introduction

Glucocorticoids are widely used to treat connective tissue and inflammatory disorders, such as systemic lupus erythematosus, Sjogren’s syndrome, rheumatoid arthritis, and severe acute respiratory syndrome [1,2]. However, one of the most common complications following administration of steroids is steroid-induced osteonecrosis of the femoral head (SONFH) [3]. Recent studies have shown that glucocorticoids are the primary risk factor for nontraumatic ONFH in both China and Japan [4,5]. SONFH is a progressive disease, and when untreated, it can spark the devastation and dysfunction of hip joints and eventually impair quality of life. As a consequence of this disease, more than half of the patients require artificial joint replacement. Regretfully, the pathogenesis of this disease is unclear. Therefore, exploring the molecular etiological factors and discovering novel biomarkers for early diagnosis and individualized therapy are urgently needed.

Apoptosis, necrosis, pyrolysis of osteoblasts and osteocytes, intravascular fat embolism, elevated intraosseous pressure, circulatory impairment, coagulation disorders, and genetic polymorphisms have been confirmed to play crucial roles in the pathogenesis of SONFH [6–10]. Several studies have focused on the role of inflammatory reaction in the development of ONFH [11,12]. It has been found that abnormal immune cell infiltration, inflammatory pathway activation, and cytokine release might be associated with the pathophysiology of ONFH [13,14]. However, few studies have investigated immune cell infiltration of SONFH. Thus, exploring the distribution of immune cells infiltrating the femoral head in SONFH may contribute to the clarification of the pathogenesis of SONFH.

Bioinformatics is an emerging interdisciplinary field involving molecular biology and information science [15]. A large number of researchers have performed chip sequencing or high-throughput sequencing experiments and uploaded the results to openly accessed repositories, such as the ArrayExpress and Gene Expression Omnibus (GEO) database [16,17]. In our study, we downloaded the primitive microarray dataset GSE123568 from the GEO database and analyzed it using bioinformatics tools. Moreover, immune cell infiltration in SONFH samples was evaluated using the ImmuCellAI method, which is applied to evaluate the abundance of 24 immune cell types. Moreover, the expression levels of key genes were examined using qRT-PCR. Additionally, drugs with potential therapeutic effects on SONFH were predicted based on the differentially expressed genes (DEGs) identified. Ultimately, the effects of aspirin were validated in SONFH treatment using a glucocorticoid-mediated mouse model of ONFH.

Thus, the aim of our research was to identify crucial biomarkers of DEGs and immune infiltration in SONFH, and explore drugs with potential therapeutic effects on SONFH. Eventually, novel diagnostic and therapeutic measures should be established for SONFH.

Materials and methods

Microarray data and identification of DEGs

The microarray data of the GSE123568 dataset were downloaded from the GEO database (GPL15207 platform, Affymetrix Human Gene Expression Array). GES123568 comprises 40 samples of 30 peripheral serum samples obtained from SONFH patients and 10 peripheral serum samples from healthy patients.

NetworkAnalyst (https://www.networkanalyst.ca/), an online tool [18], was used to explore the DEGs by comparing the expression values between SONFH patient samples and normal control samples. In the present study, an adjusted p value less than 0.01 and log FC value greater than 1.0, were selected as the cutoff criteria to identify DEGs. Additionally, we generated a volcano map based on adjusted p and log FC values. Moreover, we applied a hierarchical clustering analysis method to evaluate the DEGs (50 genes are shown) between SONFH samples and control samples, and the results were visualized using NetworkAnalyst.

Function enrichment analysis and PPI network construction

Metascape (http://metascape.org) is a free-of-charge, data-rich, and well-maintained online tool for gene annotation and analysis [19]. It can provide gene enrichment information according to biological categories, including biological process (BP), cellular component (CC), molecular function (MF), and pathway enrichment annotation to deliver exhaustive and elaborate information on each gene. Therefore, in this study, Gene Ontology (GO) terms, consisting of BP, CC, MF, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, were used to evaluate enriched genes in Metascape. Only terms with p < 0.01, a minimum count of 3, and an enrichment factor > 1.5 were considered significant. Protein – protein interaction (PPI) enrichment analysis was performed using Metascape. Furthermore, the molecular complex detection (MCODE) algorithm was used to investigate highly related network units.

Immune cell infiltration as determined via ImmuCellAI analysis

The ImmuCellAI (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) algorithm [20] was used to assess the abundance of 24 immune cell types in SONFH and normal samples. These immune cells include 18 T cell subsets and six other critical immune cells, including B cells, macrophages, monocytes, neutrophils, dendritic cells (DCs), and natural killer (NK) cells [17]. We first uploaded the gene expression files (GSE123568) obtained previously. Then, immune cell abundance in the groups was set up for analysis selection. The abundance of each immune cell type in the two groups was calculated using this method.

Hub genes identification and validation in vitro

Hub genes (10 most expressed genes) were obtained based on maximal clique centrality (MCC) evaluation by cytoHubba, a plugin in Cytoscape software [21]. Peripheral serum samples from six SONFH patients and six healthy control adults were harvested for qRT-PCR validation. The study protocol was approved by the Research Ethics Commission of Wuhan Union Hospital, and informed consent was obtained from all participants.

Quantitative real-time-PCR (qRT-PCR)

The mRNA transcripts were quantified by qRT-PCR using AceQ® Universal SYBR qPCR Master Mix (Vazyme Biotech, Q511-02) and a quantitative real-time PCR system (StepOnePlus, ABI, USA). The amplification conditions were set as follows: 94°C for 4 min, followed by 35 cycles at 94°C for 20 s, 60°C for 30 s, and 72°C for 30 s. Then, a melting curve was established to obtain the experimental data. GAPDH was used as the reference gene, and all experiments were conducted in triplicates. Relative target gene expression levels were calculated using the 2^−ΔΔCt method. The primers used in this study are listed in Supplementary Table S1.

Prediction of potential therapeutic drugs

A connectivity map (https://clue.io/query) was used to identify potential therapeutic drugs for SONFH based on uploaded DEG data. This connectivity map, an online tool, is used to compare queried signatures with a gene expression profile database of several cell lines treated with more than 2000 compounds, most of which have been approved by the U.S. Food and Drug Administration [22]. Compounds with therapeutic potential usually affect biological processes opposite of the disease. Thus, these compounds may have certain effects.

Animals and experimental procedures

The research protocol was approved by the Committee on Ethics of Animal Experiments of the Wuhan Union Hospital. Male BALB/c mice weighing 18.2–26.1 g (8 weeks) were purchased from the Laboratory Animal Center of Huazhong University of Science and Technology (Wuhan, China). They were then randomly divided into three groups: the control group, glucocorticoid (GC) group, and GC+ aspirin group, with 5 mice in each group. The SONFH model was induced by treating the mice with 2 mg/ml dexamethasone (Dex, National Medicine Standard H41020055, Zhengzhou Cheuk-Fung Pharmaceutical Company, China) administered in drinking water for 8 weeks. For the GC+ aspirin treatment group, mice were treated with aspirin (10 mg/kg) daily via gavage. The diagnosis of osteonecrosis was established based on the presence of empty lacunae or osteocytes with pyknotic nuclei in the bone trabeculae, accompanied by surrounding bone marrow necrosis [23]. At the end of the treatment, mice were sacrificed and samples of the femoral head were obtained for RNA extraction and histological examination. For histological examination, femoral head samples were decalcificated, paraffin embedded, sectioned and stained with hematoxylin and eosin or TRAP staining kit using standard procedures [24,25]. For RNA extraction, femoral head samples were collected, immediately snaped frozen in liquid nitrogen, and RNA was extracted as described above.

Statistical analysis

All results are presented as the mean ± standard deviation (SD) from at least three independent experiments. Statistical analysis was carried out using SPSS software (version 23.0). Statistical significance was determined using a two-tailed Student’s t-test when two groups were compared, and one-way ANOVA was followed by Bonferroni’s post hoc test when comparing more than two groups. GraphPad Prism 7 software was used for drawing the diagrams. A p value lower than 0.05, was considered to represent a statistically significant difference.

Results

Identification of DEGs

Following the search keywords ‘osteonecrosis of the femoral head,’ ‘glucocorticoid’ and ‘expression profiles by array’ we obtained the GSE123568 dataset from the GEO database. A total of 429 DEGs (Supplementary Table S2), including 214 upregulated genes and 215 downregulated genes, were identified in the SONFH samples compared with the normal controls using NetworkAnalyst. A volcano diagram and heat map of the DEGs are shown in Figure 1.

Figure 1. Identification of DEGs in SONFH. (a) Hierarchical clustering heat map of SONFH samples and control samples based on identified DEGs (50 genes shown). (b) Volcano plot of the DEGs between SONFH samples and control samples

Functional enrichment analysis and PPI network construction

The functions of the DEGs were assessed by conducting GO and KEGG analyses in Metascape. GO functional enrichment analysis in the BP category revealed that the DEGs were mainly involved in myeloid leukocyte activation, defense response to other organisms, positive regulation of cytokine production, and myeloid cell differentiation (Figure 2a). In the CC category, the DEGs were mainly enriched in the secretory granule membrane, side of membrane, lytic vacuole, and secretory granule lumen (Figure 2b). In addition, in the MF category, the DEGs were predominantly enriched in pattern recognition receptor activity, lipid binding, oxidoreductase activity, and cytokine receptor activity (Figure 2c). KEGG pathway analysis also demonstrated that the DEGs were notably enriched in osteoclast differentiation, leishmaniasis, chemokine signaling pathway, and cytokine-cytokine receptor interaction (Figure 2d). In addition, to better explain the correlation among each DEG, a Metascape protein–protein interaction (PPI) enrichment analysis was performed. The PPI network and MCODE components in the gene list are illustrated in Figure 2e, f. The six most significant MCODE components were extracted from PPIs. After pathway and process enrichment analysis, which was independently performed for each MCODE component, the results revealed that the biological functions were chiefly associated with protein refolding, GPCR ligand binding, MyD88-dependent Toll-like receptor signaling pathway, G alpha (q) signaling events, negative regulation of cell proliferation, and regulation of osteoblast differentiation (Table 1).

Table 1. The six most significant MCODE components were extracted from the PPI

MCODE_ALL	R-HSA-373076 Class A/1 (Rhodopsin-like receptors); R-HSA-500792 GPCR ligand binding; GO:0002274 myeloid leukocyte activation
MCODE_1	GO:0042026 protein refolding; GO:0019320 hexose catabolic process; hsa00010 Glycolysis/Gluconeogenesis
MCODE_2	R-HSA-373076 Class A/1 (Rhodopsin-like receptors); R-HSA-418594 G alpha (i) signaling events; R-HSA-500792 GPCR ligand binding
MCODE_3	GO:0002755 MyD88-dependent toll-like receptor signaling pathway; GO:0032757 positive regulation of interleukin-8 production; GO:0032677 regulation of interleukin-8 production
MCODE_4	R-HSA-416476 G alpha (q) signaling events; R-HSA-373076 Class A/1 (Rhodopsin-like receptors); R-HSA-500792 GPCR ligand binding
MCODE_5	GO:0008285 negative regulation of cell proliferation
MCODE_6	GO:0045669 positive regulation of osteoblast differentiation; GO:0045778 positive regulation of ossification; GO:0045667 regulation of osteoblast differentiation

Figure 2. Function enrichment analysis and PPI network construction. (a) GO functional enrichment analysis of BP. (b) GO functional enrichment analysis of CC. (c) GO functional enrichment analysis of MF. (d) KEGG pathway analyses results of the DEGs. (e) The PPI network of the DEGs. (f) The six most significant modules were obtained from the PPI network

Immune cell infiltration analysis

Due to scientific and technological limitations, the pattern of immune cell infiltration in SONFH has not been completely clarified, particularly in subgroups consisting of few cells. By utilizing the ImmuCellAI algorithm, we first investigated the differences in immune cell infiltration among 24 subpopulations of immune cells in these two groups (Figure 3a, b). Compared with the control group, the SONFH group displayed an elevated percentage of naive CD8 + T cells and macrophages, whereas the proportions of exhausted T cells, induced regulatory T cells (iTregs), T helper 1 cells (Th1 cells), effector memory T cells, mucosal-associated invariant T cells (MAIT cells), DCs, and gamma delta T cells were lower (Figure 3c, p < 0.05). Furthermore, we investigated the correlation between these types of immune cells. Macrophages were highly correlated with monocytes (Pearson correlation = 0.86), indicating that the functions of macrophages and monocytes may be synergistic (Figure 3d).

Figure 3. The differences in immune cell infiltration between SONFH and normal controls. (a) The relative percentage of 24 subpopulations of immune cells in 40 samples from the GSE123568 dataset. (b) Correlation analyses of the 24 subpopulations of immune cells estimated in 40 samples from GSE123568. (c) Heat map of the 24 subpopulations of immune cells estimated in 40 samples from GSE123568. (d) The differences in immune cell infiltration between SONFH patients and normal controls (*p values < 0.05 were considered statistically significant)

Hub gene identification and validation in vitro

Based on the MCC centrality evaluation by cytoHubba, the 10 most highly expressed hub genes CCR1, CCR2, CCR3, FPR2, CXCR2, CXCR1, CXCL5, PF4, HCAR2, and P2RY13, were identified (Figure 4a). We then detected the expression of these 10 hub genes in six SONFH and six control samples via qRT-PCR. The results illustrated that the comparative expression levels of all 10 hub genes were in line with the microarray hybridization results (Figure 4b).

Figure 4. The 10 most highly expressed hub genes in the PPI network as confirmed with RT-PCR. (a) The 10 most highly expressed hub genes in the PPI network as determined by cytoHubba. (b) RT-PCR validation of the hub genes in SONFH and normal controls. All experiments were performed in triplicate, and the results are presented as the means ± SD. (*p < 0.05)

Prediction of potential therapeutic drugs for SONFH

To identify potentially therapeutically valuable compounds that may be the most suitable for targeting DEGs, we uploaded DEGs to the connectivity map L1000 platform, a tool used to compile gene expression profiles associated with a variety of therapeutic compounds. The 10 compounds most closely associated with DEGs are listed in Table 2 with their gene targets, including phylloquinone, cholic acid, MRS-1220, bucladesine, isotretinoin, doxercalciferol, betahistine, aspirin, maraviroc, and SR-27897. Among the target genes, we found that PTGS2, a gene target of aspirin, was significantly upregulated in the GSE123568 dataset, indicating a valuable role for aspirin in SONFH treatment.

Table 2. The top 10 compounds with activity against SONFH as predicted via connectivity map

Score	Name	Description	Target
−99.68	Phylloquinone	Vitamin K	BGLAP, GGCX
−99.58	cholic-acid	Bile acid	ADH1C, CES1, COX4I1, COX5A, COX5B, COX6A2, COX6B1, COX6C, COX7A1, COX7B, COX7C, COX8A, ESRRG, FABP6, FECH, GPBAR1, MT-CO1, MT-CO2, MT-CO3, PLA2G1B
−99.58	MRS-1220	Adenosine receptor antagonist	ADORA2B, ADORA3
−99.54	Bucladesine	Adenosine receptor agonist	PRKACA
−99.51	Isotretinoin	Retinoid receptor agonist	CYP2B6, CYP2C19, CYP2C8, CYP3A5, CYP3A7, NR2C2, PPARD, RARA, RARB, RARG, RORB
−99.22	Doxercalciferol	Vitamin D receptor agonist	VDR
−98.8	Betahistine	Histamine receptor agonist	HRH1, HRH3
−98.8	Aspirin	Cyclooxygenase inhibitor	AKR1C1, ASIC3, EDNRA, HSPA5, IKBKB, NFKB1, NFKB2, NFKBIA, PRKAA1, PRKAA2, PRKAB1, PRKAB2, PRKAG1, PRKAG2, PRKAG3, PTGS1, PTGS2, RPS6KA3, TP53
−98.73	Maraviroc	CC chemokine receptor antagonist	CCR5, CYP3A5
−98.38	SR-27897	CCK receptor antagonist	CCKAR

A mouse model of glucocorticoid (GC)-induced ONFH was generated to further explore the role of aspirin in the treatment of SONFH (Figure 5a). Compared with mice in the GC treatment group, GC+ aspirin treatment ameliorated GC-induced ONFH in mice by increasing bone volume/total volume (BV/TV) and reducing the empty lacunae rate (Figure 5b-f). However, aspirin treatment did not inhibit osteoclastogenesis in GC-induced ONFH in mice. As shown in Figure 5g-I, compared with that in the GC treatment group, the number of osteoclasts in the femoral head tissue in the GC+ aspirin treatment group was not significantly reduced by TRAP staining. Furthermore, relative Runx2 mRNA expression in the femoral head tissue was significantly increased by aspirin treatment (Figure 5j). Relative RANKL mRNA expression in the femoral head tissue was significantly decreased by aspirin treatment, compared with the GC treatment group (Figure 5k). Since PTGS2 was predicted to be the target gene of aspirin, PTGS2 mRNA expression was detected by qRT-PCR in different groups of mice. The results revealed that aspirin treatment reversed GC treatment-induced increase in PTGS2 expression (Figure 5l). Taken together, these results indicate that aspirin may be a drug for the treatment of SONFH.

Figure 5. Identification of the therapeutic effect of aspirin in GC-induced ONFH. (a) Schematic diagram of the treatment of different groups of mice. H&E staining of representative femoral heads in each group. (b) Bone volume/total volume ration (BV/TV) in each group. (c) Empty lacunae rate in each group. (d-f) H&E staining of representative femoral heads in each group (scale bar = 100 μm). (g-i) TRAP staining of representative femoral heads in each group (scale bar = 100 μm). (j-l) Relative Runx2, RANKL and PTGS2 mRNA expression in femoral heads in each group. All experiments were performed in triplicate, and the results are presented as the means ± SD. (*p < 0.05 versus the control; #p < 0.05 versus the GC)

Discussion

The medical history of glucocorticoid use and related clinical manifestations, such as pain in the hip joint, and imaging examinations, including those on plain film, bone scintillation, and magnetic resonance imaging, are often used for the diagnosis of SONFH [26]. However, symptoms may not appear in the early stages of the disease, but rather in the middle and late stages. Because the disease cannot be diagnosed early, the hip preservation treatment is delayed. Although magnetic resonance imaging is often used as a method for the early diagnosis of this disease, it is not suitable for screening all patients receiving GC therapy because of the high cost and time commitment required for obtaining scans. Thus, identifying key circulating markers is essential for the early diagnosis of SONFH. In recent years, some pivotal biomarkers for SONFH have been explored, such as type I collagen cross-linked C-telopeptide and amino terminal telopeptide of procollagen type I [27]. However, both two genes were found to be associated with the development of osteoporosis, which limited their diagnostic value [28]. Therefore, it is necessary to identify promising biomarkers for SONFH to improve its diagnosis and treatment. In our study, we attempted to identify potentially crucial genes related to SONFH. A bioinformatics assessment of GO and KEGG enrichment analyses was performed, and pivotal hub genes were identified in our study. Moreover, we adopted the ImmuCellAI tool to identify immune cell infiltration in SONFH. The results showed a remarkable difference in immune cell infiltration between the SONFH and control samples.

In the present study, we discovered that the DEGs were mostly enriched in myeloid leukocyte activation, positive regulation of cytokine production, myeloid cell differentiation, and cytokine receptor activity. Previous studies have verified that the cytokines IL-10, IL-12, and TNF-α are related to the development of ONFH [29]. In addition, osteoprotegerin, receptor activator of nuclear factor-κB and receptor activator of nuclear factor-κB ligand (RANKL), a cytokine/cytokine receptor signal axis, have been confirmed to regulate the balance between osteogenesis and osteoclastogenesis, which has also been found to be associated with the pathophysiology of ONFH [30]. Additionally, the outcomes of KEGG pathway analysis revealed that DEGs were notably enriched in osteoclast differentiation, chemokine signaling pathways, and cytokine-cytokine receptor interactions. Moreover, six of the most significant MCODE components extracted from the PPI network were also enriched in the MyD88-dependent Toll-like receptor signaling pathway and in cell proliferation and osteoblast differentiation regulatory pathways. Two studies demonstrated that the immune response associated with the Toll-like receptor 4 signaling pathway can lead to SONFH [31,32]. In addition, abnormal osteoclast differentiation of marrow mesenchymal stem cells (MSCs) has been confirmed to be significantly associated with the development of SONFH [33–35]. Thus, we hypothesize that the immune response associated with cytokine and chemokine signaling pathways and osteoclast and osteoblast differentiation may underlie the potential pathophysiology of SONFH.

Ten crucial hub genes related to SONFH were identified from the PPI network, namely, CCR1, CCR2, CCR3, CXCR1, FPR2, HCAR2, PF4, CXCR2, CXCL5, and P2RY13, using the Cytoscape plugin named CytoHubba. In addition, we collected six serum samples from patients with SONFH and six control serum samples to explore the expression of 10 key pivotal genes and found similar trends with their expression levels in the GEO datasets, which confirmed the reliability of our findings. Previous studies have indicated the role of some hub genes in bone and joint diseases. For instance, CCR2 participates in chondrocyte degradation and leads to the progression of osteoarthritis [36,37]. In addition, it has been demonstrated that overexpression of CCR3 increases the recruitment of circulating monocytes in bone, enhances monocyte differentiation into osteoclasts, and eventually promotes the development of osteoporosis [38]. Furthermore, although CXCR1 and CXCR2 are members of a family of chemokine receptors, they seem to play opposite roles in osteoarthritis. Some researchers have pointed out that CXCR1/2 signaling can maintain cartilage homeostasis by increasing extracellular matrix production and inhibiting chondrocyte apoptosis [39,40]. As for PF4 and CXCL5, PF4 is also called CXCL4, yet their role has not been reported. CXCL5 inhibits osteoclastogenesis of MSCs [41]. However, the critical role of several key genes screened in SONFH remains unclear, such as CCR1, FPR2, HCAR2, and P2RY13. Further research is required to determine the exact roles of these hub genes in SONFH. We believe that these results can provide useful knowledge for improving the accuracy of diagnosis and treatment in patients with SONFH.

From the analysis of immune cell infiltration, we also found significant differences in the abundance of several immune cells, such as naive CD8 + T cells and macrophages. The effects of several types of immune cells in SONFH have also been explored. For example, previous studies have shown a high correlation between increased numbers of macrophages and glucocorticoid administration in patients with SONFH, which demonstrated that macrophage infiltration might be one of the causes of SONFH [42]. In addition, inflammatory infiltration of macrophages and CD4 + T cells was also observed in patients with SONFH [14]. However, the role of naive CD8 + T cell infiltration in SONFH requires further study. We also investigated the correlation between these immune cells. Macrophages were highly correlated with monocytes. Immune cell infiltration analysis showed that monocyte levels were higher in the SONFH group than in the control group, although the p value was slightly greater than 0.05 (p = 0.053). These results reveal that macrophages and monocytes have a synergistic effect in the development of SONFH. Many researchers have reported that macrophages and monocytes are precursors of osteoclasts [43,44]. Osteoclast differentiation-mediated bone resorption is one of the causes of femoral head collapse [45]. Hence, we speculate that high-dose glucocorticoids might activate macrophage and monocyte systems in the blood, which will infiltrate the femoral head tissues, increasing local osteoclast differentiation and leading to the destruction of the balance between bone formation and bone resorption.

In the past few decades, conservative treatment of SONFH has mainly focused on preventing collapse of the femoral head. Administration of bisphosphonate was reported to prevent collapse in early stage ONFH in a small group of patients [46]. However, the prolonged anti-collapse effects of bisphosphonates have not been elucidated. Once the femoral head progresses and collapses, surgery is needed, including hip-preserving and hip replacement surgeries. However, surgical treatment often results in greater trauma and higher costs. Therefore, it is necessary to identify novel agents and targets for the early treatment of SONFH. In this study, we screened 10 potential drugs for the treatment of SONFH. Most of them have been proven to be associated with the maintenance of normal bone metabolism and have potential therapeutic effects in SONFH. Considering the target genes of these compounds, aspirin, which targets PTGS2, attracted our attention. Because of its safety, aspirin is often used clinically as an antipyretic, analgesic, and anti-inflammatory drug. Previous studies have demonstrated its potential effects on postmenopausal osteoporosis [47]. Aspirin can inhibit osteoclast formation and improve osteogenesis by affecting a variety of biological pathways, such as inhibiting the activities of COX2, COX1, and PGE2 [48]. Furthermore, some studies have shown that aspirin may be a simple and effective treatment option to delay the progression of early ONFH in patients [49]. In this study, a mouse model of Dex-induced ONFH was developed to further explore the role of aspirin in SONFH treatment, and the results revealed that aspirin could protect against GC-induced SONFH, indicating the therapeutic effect of aspirin in SONFH. However, our results demonstrated that aspirin did not ameliorate SONFH by inhibiting osteoclastogenesis.

Although several critical hub genes, types of infiltrating immune cells, and potential therapeutic agents for SONFH were identified in our study, there are some limitations to our study. First, the sample size was relatively small, which may have led to low statistical power for examining gene expression differences between the two groups. With the establishment of the steroid-induced femoral head necrosis specimen library, we will collect more specimens for testing related indicators in the future to reduce the error caused by the small sample size. Second, although we confirmed in a mouse model that aspirin could inhibit the occurrence of SONFH, follow-up clinical trials are still needed to confirm its clinical efficacy. Finally, further studies on the molecular mechanism by which aspirin inhibits SONFH are needed.

Conclusion

In summary, 429 DEGs were identified between the SONFH and the control groups. Analysis of the biological functions and pathways of these DEGs contributed to a better understanding of the molecular mechanism of SONFH. Ten pivotal hub genes in SONFH were identified, and the difference in infiltrating immune cells between SONFH and normal control samples was explored. In addition, considering the DEGs identified, we successfully predicted that 10 potential drugs might have therapeutic effects on SONFH. Finally, the role of aspirin in SONFH treatment was verified using a mouse model of DEX-induced ONFH. Overall, our findings might provide novel clues for investigating the pathogenesis and treatment of SONFH.

Highlights

1. Ten hub genes and infiltrating immune cells were identified in SONFH.

2. Myeloid leukocyte activation may be the causes of SONFH.

3. Aspirin may be a novel treatment for SONFH.

Authors’ contributions

YF and WX conceived and designed the study. BW and SG performed the experiments. LH, WS, and ZL performed the data analysis. BW wrote the manuscript. ZZ and YZ helped in designing the figures. XW and YM contributed to manuscript editing. YF supervised the study. All authors and participants reviewed and approved the final manuscript.

Abbreviations

SONFH, steroid-induced osteonecrosis of the femoral head; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes; PPI, protein-protein interaction network; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCODE, molecular complex detection; MCC, maximal clique centrality; GC, glucocorticoid; DEX, dexamethasone; MSCs, marrow mesenchymal stem cells; RANKL, receptor activator of nuclear factor-κB ligand.

Ethical statement

The present study was reviewed and approved by the Research Ethics Commission of Wuhan Union Hospital.

Data availability statement

All the data in our study will be made available by the authors upon reasonable request. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123568

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (No. 81974337, No. 81672155) and the Project of Scientific Research Plan of Wuhan Municipal Health Commission (WX21Q19).