Bioinformatics analysis and identification of genes and molecular pathways in steroid-induced osteonecrosis of the femoral head

Background Steroid-induced osteonecrosis of the femoral head (ONFH) is a common hip joint disease and is difficult to be diagnosed early. At present, the pathogenesis of steroid-induced ONFH remains unclear, and recognized and effective diagnostic biomarkers are deficient. The present study aimed to identify potentially important genes and signaling pathways involved in steroid-induced ONFH and investigate their molecular mechanisms. Methods Microarray data sets GSE123568 (peripheral blood) and GSE74089 (cartilage) were obtained from the Gene Expression Omnibus database, including 34 ONFH samples and 14 control samples. Morpheus software and Venn diagram were used to identify DEGs and co-expressed DEGs, respectively. Besides, we conducted Kyoto Encyclopedia of Genome (KEGG) and gene ontology (GO) pathway enrichment analysis. We construct a protein-protein interaction (PPI) network through GEO2R and used cytoHubba to divide the PPI network into multiple sub-networks. Additionally, quantitative real-time polymerase chain reaction (qRT-PCR) was performed to verify the bioinformatics analysis results. Results A total of 118 intersecting DEGs were obtained between the peripheral blood and cartilage samples, including 40 upregulated genes and 78 downregulated genes. Then, GO and KEGG pathway enrichment analysis revealed that upregulated DEGs focused on the signaling pathways related to staphylococcus aureus infection, leishmaniasis, antigen processing, and presentation, as well as asthma and graft-versus-host disease. Downregulated genes were concentrated in the FoxO signaling pathway, AMPK signaling pathway, signaling pathway regulating stem cell pluripotency, and mTOR signaling pathway. Some hub genes with high interactions such as CXCR1, FPR1, MAPK1, FOXO3, FPR2, CXCR2, and TYROBP were identified in the PPI network. The results of qRT-PCR demonstrated that CXCR1, FPR1, and TYROBP were upregulated while MAPK1 was downregulated in peripheral blood of steroid-induced ONFH patients. This was consistent with the bioinformatics analysis. Conclusions The present study would provide novel insight into the genes and associated pathways involved in steroid-induced ONFH. CXCR1, FPR1, TYROBP, and MAPK1 may be used as potential drug targets and biomarkers for the diagnosis and prognosis of steroid-induced ONFH.


Introduction
Osteonecrosis of the femoral head (ONFH) is a common disease of the hip with a higher incidence in elderly patients, and the most common clinical symptom is severe pain [1]. There are about 150,000-200,000 new cases in China each year [2], and about 20,000 new cases in the USA each year [3]. The pathological features of steroidinduced ONFH include decreased blood supply to the hip joint, the collapse of the femoral head, and accumulation of microfractures without continuous remodeling [4,5]. Besides, ONFH can cause rapid destruction and dysfunction of the hip joint, and approximately 65-70% of patients with advanced ONFH require total hip replacement [6]. The pathogenesis and molecular mechanism of steroid-induced ONFH remain unclear, and there are few effective prevention and early treatment methods.
Studies have indicated that successful results can be obtained when joint-preserving surgery is used in the pre-collapse stage of steroid-induced ONFH [7]. After the femoral head necrosis collapses, cartilage wrinkles and the hip joint space gradually narrows. Meanwhile, the failure rate of joint-preserving treatment is high. Therefore, searching for new biomarkers in the blood is of great clinical significance for the diagnosis of steroidinduced ONFH, as well as the further early diagnosis and treatment of ONFH. Some genetic factors are considered key factors in the progress of steroid-induced ONFH, such as peroxisome proliferator-activated receptor γ, matrix metallopeptidase 9, and SMAD family member 3 and gremlin 1 [8][9][10]. Currently, microarray analysis using a high-throughput platform is promising and efficient for exploring the molecular mechanisms of diseases and identifying useful biomarkers for diagnosis and prognosis of diseases [11,12]. However, there is little comprehensive research on steroid-induced ONFH with high-throughput platforms. During the development of steroid-induced ONFH, the cartilage on the surface of the femoral head and the labrum of the acetabulum are significantly damaged. Degeneration of cartilage on the surface of the femoral head and tearing of the labrum enhance the instability of the hip joint and accelerate the development of steroid-induced ONFH [13]. Prevention and early treatment of hip cartilage damage may slow the development of steroid-induced ONFH and relieve hip dysfunction. However, there are few studies on the molecular mechanism of steroidinduced ONFH hip cartilage damage. Therefore, it is hypothesized in this paper that differentially expressed genes co-exist in the peripheral serum and cartilage of ONFH patients; this can be used for early diagnosis and observation of disease progression.
In this study, hub genes are screened out from the differences expressed genes (degrees) using bioinformatics methods based on the National Center for Biotechnology Information (NCBI http://www.ncbi.nlm.nih.gov/ geography/) from the Gene Expression Dataset GSE123568 (peripheral serum) and GSE74089 (hip articular cartilage) Comprehensive (GEO) database. Quantitative real-time PCR (qRT-PCR) was conducted to validate the gene expression profiling results based on an independent sera sample of eight patients with ONFH and twelve healthy controls. The present study aimed to identify potentially important genes and signaling pathways involved in ONFH and investigate their molecular mechanisms.

Microarray data source
The gene chip data related to femoral head necrosis (GSE123568, GSE74089) is downloaded in the GEO (gene expression omnibus) database (https://www.ncbi. nlm.nih.gov/geo/). GSE123568 chip data contains 10 normal patient serum samples and 30 steroid-induced femoral head necrosis serum samples. Besides, the GSE74089 chip data is composed of hip cartilage specimens from 4 NFH patients and 4 healthy controls. All the above data sets are based on GPL15207 and GPL13497 platform GEO database downloaded, and the information of the above two data sets are detailed in Table 1.
Differentially expressed genes analysis GEO 2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) is an online analysis program based on the R language provided by the GEO database (R 3.2.3 version) [14]. The Bioconductor R package in GEO2R is used to perform background correction, normalization, and expression value calculation on the GSE123568 and GSE74089 data; the limma 3.26.8 package is adopted to calculate the differential genes between the two groups and derive the GEO2R processed data. Set P< 0.01, adjust P value < 0.01, and the criteria for screening differential genes is that the expression change range is greater than or equal to twice (|log2 FC|≥1.0). Besides, log2FC≥1.0 indicates upregulation of gene expression; log2FC ≤− 1.0 suggests downregulation of gene expression. Finally, the differentially expressed genes of the femoral head necrosis group and the healthy control group, namely the differentially expressed genes (DEGs) of femoral head necrosis, were obtained. Using the online analysis website ClustVis (https://biit.cs.ut.ee/clustvis/) [15], the heat map and cluster analysis of the selected differential genes are completed, and the -log10 conversion is performed for the adjust P value in the data processed by GEO2R. Then, according to log2 FC, -log10 (adjust P value) was divided into groups of upregulated genome, downregulated genome, and non-statistically different genome. Next, the processed data is imported into GraphPad Prism 7 to draw a heat map. Finally, the Venn diagrams of the common upregulated and downregulated DEGs in the two data sets is illustrated using the online Venn diagram production website bioinfogp (https://bioinfogp. cnb.csic.es/tools/venny/index.html) [16].

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis
The online database DAVID (DAVID Bioinformatics Resources 6.8) was used to analyze the function and pathway enrichment of the upregulated and downregulated common differentially expressed genes obtained from the intersection of the two chip data sets to further explore the potential pathogenesis of the disease. The upregulated and downregulated differentially expressed genes are imported according to the set differential gene screening criteria from the data after GEO2R analysis and processing into the online enrichment analysis website DAVID (https://david.ncifcrf.gov/tools.jsp) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) signal pathway analysis [17,18]. The GO analysis mainly includes the cellular components (CC), molecular functions (MF), and biological processes (BP) of differential genes.

Protein-protein interaction (PPI) network analysis
The online analysis website STRING was used to construct a protein-protein interaction network (proteinprotein interaction, PPI). The screened upregulated and downregulated common differentially expressed genes were uploaded to the STRING online analysis website (https://string-db.org/); "highest confidence > 0.700" was set in the lowest interaction score, and the proteinprotein interaction result data were exported. Furthermore, the protein-protein interaction network (PPI network) of differentially expressed genes was obtained using Cytoscape 3.7.2 software [19], and the cytoHubba plug-in was adopted to screen Hub genes [20]. Meanwhile, three different algorithms (MNC, Degree, and MCC) were used to calculate the top 10 upregulated differentially expressed genes and the top 10 downregulated differentially expressed genes in their respective scores. Besides, the key genes are obtained by the intersection of the gene sets obtained by three different algorithms with the online analysis platform Draw Venn Diagram (http://bioinformatics. psb.ugent.be/webtools/ Venn/).

Quantitative real-time PCR analysis
Total RNA from peripheral serum samples was isolated with Trizol LS. The protocol for cDNA synthesis was reverse transcription at 37°C for 15 min and at 85°C for 5 s. The qPCR protocol was 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 30 s. The relative expression level was determined as targeting genes divided by β-Actin (ACTB). Relative miRNA expression was generated with the 2 −△△Cq method. Primers used in this study are presented in Table 2. All experiments were performed independently in triplicate.

Statistical analysis
GraphPad 8.0 was used to perform statistical analysis.
The results are presented as mean ± standard deviation (mean ± SD). Data conforming to normal distribution were compared using Student t test, while those with non-normally distributed were tested using Mann-Whitney U test. P < 0.05 indicates statistical significance. All experiments were performed in triplicate.

Identification of DEGs
This study screened out two chip data sets (GSE123568 and GSE74089). Among them, 30 serum samples of patients with steroid hormone-induced osteonecrosis of the femoral head in the chip data set GSE123568 (including 13 males and 17 females, aged 23.07 ± 3.01 years, ARCO stage I-II patients accounted for 10 patients, ARCO III patients accounted for 10 cases, ARCO IV patients accounted for 10 cases) and 10 normal serum samples (including 7 males and 3 females, aged 25.02 ± 2.87 years). Besides, there are 4 hip cartilage samples from patients with femoral head necrosis   Table 3. According to the adjust P value < 0.01, the expression change range is greater than or equal to twice (|log2 FC|≥1.0) as the criteria for screening differential genes. A total of 750 differentially expressed genes were screened in the GSE123568 data set, including 454 upregulated genes and 296 downregulated genes. A total of 5211 differentially expressed genes were screened in the GSE74089 data set, including 1594 upregulated genes and 3617 downregulated genes. The DEGs in GSE123568 and GSE74089 were selected based on the adjust P value to screen the top 25 most significant upregulated and downregulated differentially expressed genes, so as to draw heat maps (Fig. 1). The -log10 conversion is performed for the adjust P value in the data processed by GEO2R. Then, according to log2 FC, -log10 (adjust P value) was divided into groups of upregulated genome, downregulated genome, and non-statistically different genome. The processed data was imported into Graph-Pad Prism 7 to draw a volcano map (Fig. 2). Furthermore, a total of 118 intersecting DEGs were obtained between the peripheral blood and cartilage samples, including 40 upregulated genes and 78 downregulated genes (Fig. 3).

Enrichment analysis of KEGG pathway and GO terms
In this study, the DAVID database was used to perform GO and KEGG pathway enrichment analysis on the upregulated and downregulated genes. GO analysis of upregulated DEGs demonstrates that the biological process (BP) is concentrated in inflammation, immune response, and signal transduction. Cellular component (CC) focuses on plasma membrane, endoplasmic reticulum membrane, and lysosome. Molecular functions (MF) are manifested in immune regulation process, IL-8 receptor activity, superoxide-producing NADPH oxidase activity, and MAP kinase activity (Table 4, Fig. 4). Besides, GO analysis of downregulated genes revealed that BP is rich in cell cycle regulation, apoptosis process, translation process regulation, autophagy, and cell response to hypoxia. CC is mainly condensed in cytoplasm, nucleus, mitochondria, autophagosomes, and exosomes. MF is manifested in protein binding, transcription factor binding, and DNA binding. The GO analysis of the top 10 downregulated genes with the smallest adjust P value is illustrated in Table 5 and Fig. 5. Enrichment analysis of the upregulated gene KEGG pathway revealed that it focuses on the signaling pathways related to staphylococcus aureus infection, leishmaniasis, antigen processing and presentation, asthma, and graft-versus-host disease. Downregulated genes are concentrated in the FoxO signaling pathway, AMPK signaling pathway, signaling pathway regulating stem cell pluripotency, and mTOR signaling pathway (Table 6, Fig. 6).

Construction of PPI network and identification of hub genes
Three different algorithms (MNC, MCC, and Degree) in the Cytoscape software's cytoHubba plug-in were used to calculate top 10 upregulated key genes and top 10 downregulated key genes (Table 7 and Table 8, Fig. 7 and Fig. 8) and obtain the intersection. Hub genes for bone necrosis are as follows: (1)

Validation of the hub genes in ONFH patients
Among the identified DEGs, CXCR1, FPR1, MAPK1, FOXO3, FPR2, CXCR2, and TYROBP were selected to verify the integrated result. Peripheral serum samples from 12 steroid-induced ONFH patients and 8 healthy adults were collected to determine mRNA levels. The result revealed that CXCR1, FPR1, and TYROBP were upregulated, and MAPK1 was downregulated in peripheral blood of steroid-induced ONFH patients (P < 0.05).
There was no significant difference in the expression of FOXO3, FPR2, and CXCR2 in peripheral serum samples between the two groups (P > 0.05) (Fig. 9). The expression of CXCR1, FPR1, MAPK1, and TYROBP were consistent with bioinformatics analysis except for FOXO3, FPR2, and CXCR2. Moreover, CXCR1, FPR1, and TYR-OBP were upregulated, and MAPK1 was downregulated in the peripheral blood of steroid-induced ONFH patients.

Discussion
Steroid-induced ONFH is a devastating and difficult joint disease and primarily caused by the interruption of blood supply and the dysfunction of the coagulation system, resulting in the death of bone cells, the decrease of the elastic modulus of the femoral head, and the collapse of the femoral head [21]. Various molecular and genetic studies have explored the etiology and pathogenesis [22,23]. However, the exact pathogenesis remains unclear.
In this paper, the genome-wide method was used to study the differential expression of genes in the peripheral blood serum and hip joint cartilage specimens of steroid-induced ONFH patients and healthy controls.    Through bioinformatics analysis of GSE123568 and GSE74089, potentially important genes in ONFH were identified. In addition, the role of the identified genes in ONFH and their interaction were analyzed. A total of 118 intersecting DEGs were obtained in peripheral blood and cartilage samples, of which 40 were upregulated genes and 78 were downregulated genes. These DEGs may be important biomarkers related to the pathogenesis and progression of ONFH. Finally, qRT-PCR verified the expression patterns of several genes such as CXCR1, FPR1, MAPK1, FOXO3, FPR2, CXCR2, and TYROBP in vitro. However, due to different inclusion criteria and a small number of patients in the validation set, some results are inconsistent with microarray analysis.
The DEGs screened in this study contribute to the identification of new molecules or new pathways related to steroid-induced ONFH. These new molecules or pathways may be essential targets for disease diagnosis and treatment. The upregulated DEGs were generally involved in inflammation, immune response, and signal transduction. Molecular functions were manifested in immune regulation process, IL-8 receptor activity, superoxideproducing NADPH oxidase activity, and MAP kinase activity. The upregulated DEGs were enriched in pathways engaged in staphylococcus aureus infection, leishmaniasis, antigen processing and presentation, asthma, and graftversus-host disease. Particularly, the background of this femoral head necrosis can be understood using corticosteroids to treat related acute or chronic graft-versus-host   disease (aGVHD, cGVHD) [24]. Socie et al. [25] retrospectively analyzed the risk factors of osteonecrosis after allogeneic hematopoietic stem cell transplantation. The study revealed that the odds ratios of aGVHD and cGVHD were 3.73 and 3.52, respectively. The severity of graft-versus-host disease was associated with a higher frequency of osteonecrosis. Additionally, GO analysis of downregulated genes demonstrated that BP focused on cell cycle regulation, apoptosis process, translation process regulation, autophagy, and cell response to hypoxia. MF was manifested in protein binding, transcription factor binding, and DNA binding. The downregulated genes were concentrated in the FoxO signaling pathway, AMPK signaling pathway, signaling pathway regulating stem cell pluripotency, and mTOR signaling pathway. Moreover, AMPK was an essential energy state sensor responding to energy supply and demand through coordinated metabolic pathways [26]. Studies have indicated that AMPK can promote cell survival under various stress conditions [27]. The forced activation of AMPK is a good strategy to protect osteoblasts from Dex [28,29]. The mTOR pathway may be a vital mediator of bone homeostasis [30]. Furthermore, the impaired differentiation of the osteogenic/adipogenic lineage of bone marrow mesenchymal stem cells (BMSCs) was caused by the activation of the mTOR signaling pathway [31]. Meanwhile, blocking the mTOR pathway can prevent the development of phenotypes related to aging and enhance the proliferation ability of BMSCs [32]. Therefore, regulating the mTOR signaling pathway may be one of the crucial methods to treat steroid-induced ONFH. The construction of the PPI network of DEGs was used to study the relationship between GO enrichment and important proteins identified by pathway analysis. The upregulated DEGs included TYROBP, CXCR1, FPR1, LPAR2, FPR2, HLA-DRB5, HLA-DRB1, CXCR2, and MNDA. TYROBP is a protein, which is involved in osteoclast differentiation and function, such as the generation of actin cytoskeleton, which is very important for bone resorption [33]. In this study, both biosynthesis analysis and specimen verification found that the expression of TYROBP was upregulated, indicating that  Fig. 9 Expression of hub genes in ONFH patients. a-c CXCR1, FPR1, and TYROBP were upregulated in peripheral blood of steroid-induced ONFH patients(P < 0.05). d MAPK1 was downregulated in peripheral blood of steroid-induced ONFH patients (P < 0.05); e-g there was no significant difference in the expression of CXCR2, FPR2, and FOXO3 in peripheral serum samples between the two groups (P > 0.05) TYROBP may play an important role in regulating osteoclast differentiation in steroid-induced ONFH. Specifically, CXCR1 (IL 8 receptor), a chemokine, is recognized to be abundant in the synovial fluid of bone marrow mesenchymal stem cells and has been verified to be a potent inducer of bone marrow mesenchymal stem cell migration [33]. Chemokines can induce the migration response of bone marrow mesenchymal stem cells to cartilage defect areas [34]. Hence, CXCR1 may play an important role in the repair of steroid-induced ONFH cartilage. Besides, FPR is a chemotactic receptor, and its expression is often identified in phagocytes (monocytes and neutrophils, etc.). FPR may play a role in promoting stem cell reproliferation because MSCs functionally express FPR and are positive for n-formyl peptides [35]. Matsushita et al. [36] found an important role of fMLP or FRP1 in MSC osteogenic and adipogenic differentiation. Moreover, the important role of fMLP is to promote the osteoblast differentiation of MSCs, and the knockdown of FPR1 inhibits this osteogenic differentiation. Downregulated DEGs identified in the present study included MAPK1, AKT2, MKRN1, FBXL4, GABARAPL2, MAP1LC3B, SIAH2, UBE2H, and FOXO3. MAPK was a type of serine/threonine protein kinase in cells that can regulate cell growth and differentiation. Experiments have demonstrated that activating the MAPK signaling pathway can promote the mineralization of osteoblasts and the expression of ALP and BMP-2, contributing to bone formation. ERK, a downstream factor of MAPK, can inhibit osteoclast formation when it is inactivated [35]. We discovered that MAPK1 was downregulated in steroid-induced ONFH blood samples compared with normal controls, providing another pathogenic role in bone disease. Moreover, AKT2 is a pro-survival protein in the AKT family and is activated through the PI3K pathway. Zhang et al. [37] revealed that AKT2 plays a crucial role in the development of zebrafish bone. Furthermore, inhibiting the PI3k/Akt signaling pathway may be one of the vital mechanisms of glucocorticoid-induced osteoblast apoptosis [38]. In the present study, the decreased expression of AKT2 was observed in steroid-induced ONFH, suggesting a vital role in bone formation and metabolism.
There are limitations to our study. First, the sample size in the qRT-PCR data set was small, and larger numbers of blood samples of ONFH patients are needed for further research. Second, the genes investigated and their pathways were not confirmed through in vitro studies or other functional studies. Therefore, the potential role of DEGs and PPIs identified in ONFH patients requires to be further investigated in vivo and in vitro. For example, the effects of downregulated or upregulated DEGs in cellular or in vivo models could be determined.
In conclusion, the present study would provide novel insight into the genes and associated pathways involved in steroid-induced ONFH. CXCR1, FPR1, TYROBP, and MAPK1 could serve as potential drug targets and biomarkers for the diagnosis and prognosis of steroidinduced ONFH.
Abbreviations ONFH: Osteonecrosis of the femoral head; KEGG: Kyoto Encyclopedia of Genome; GO: Gene Ontology; PPI: Protein-protein interaction; NCBI: National Center for Biotechnology Information