Transcriptomic analysis reveals genetic factors regulating early steroid-induced osteonecrosis of the femoral head.

The present study aimed to explore the signaling pathways involved in development of early steroid-induced osteonecrosis of the femoral head (SONFH) and identify diagnostic biomarkers regulating peripheral blood in SONFH patients. We downloaded transcriptome data and identified differentially expressed genes (DEGs) using the R software. We used ClusterProfiler to perform enrichment analysis of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes, and analyzed protein-protein interactions using the STRING database. Network X was used to visualize the networks in Python. A total of 584 DEGs were identified, of which 294 and 290 were upregulated and downregulated, respectively. Enrichment analysis showed that the DEGs were mainly involved in red blood cell differentiation, cell protein catabolism, gas transportation, activation of myeloid leukocytes, phagocytosis, and inflammatory response. Pathway analysis revealed that these DEGs were involved in regulation of mitophagy-animal, human T-cell leukemia virus-1 infection, Forkhead box O, phagocytosis, osteoclast differentiation, and cytokine-cytokine receptor interaction. Quantitative real-time polymerase chain reaction results were consistent with findings from protein-protein interaction network analysis. Several genes, including peroxiredoxin 2, haptoglobin, matrix metallopeptidase 8, formyl peptide receptor 2, and integrin subunit alpha X, promote SONFH occurrence by regulating the redox, inflammatory response, and osteoblast and osteoclast structure and function pathways. They may be important targets for designing approaches for early diagnosis and treatment of SONFH.

1. Introduction

Osteonecrosis of the femoral head (ONFH) is a common disease mainly characterized by apoptosis of bone cells that leads to collapse of the femoral head and cartilage damage. Notably, the round surface of the head and acetabulum joint flattens and eventually leads to secondary osteoarthritis.[ Studies have shown that the disease is mainly prevalent in patients with long-term use of hormonal drugs, long-term alcohol abuse, and hip trauma.[ Steroid-induced ONFH (SONFH) is a common type of femoral head necrosis. SONFH patients are frequently misdiagnosed as hip synovitis or osteoarthritis owing to lack of obvious signs coupled with imaging abnormalities during early stages of disease development.[ Although transcriptomic analyses have contributed to development of clinical treatment strategies for SONFH, the precise underlying mechanism of SONFH remains unclear.[

Bioinformatics analyses have significantly improved the search and identification of key genes regulating the occurrence of various diseases, thereby generating new insights to guide development of treatment strategies. SONFH has a complex pathophysiological mechanism. Genotyping of SONFH in adults revealed that gene coding protein polymorphism was associated with SONFH.[ It has been hypothesized that progression of SONFH to the end stage in adult patients may be due to differential expression of protein-coding genes involved in thrombosis, angiogenesis inhibition, programmed cell death, lipid biosynthesis, and bone homeostasis. In the present study, we aimed to identify differentially expressed genes in peripheral blood of patients with early SONFH and further explored the pathogenesis of this disease using transcriptome data from the Gene Expression Omnibus database.

2. Materials and Methods

2.1. Data acquisition

The Gene Expression Omnibus database not only stores many gene chips and sequence data but also provides tools for users to query and download gene expression profiles. In the present study, we searched for relevant datasets using “steroid-induced osteonecrosis of femoral head” as the keyword, with focus on SONFH-related datasets. Consequently, we downloaded transcriptome data (GSE123568) belonging to 30 SONFH and 10 non-SONFH (following steroid administration) subjects.

2.2. Sample correlation analysis

Overall expression profiles of the transcriptome data were clustered using ClusterProfiler packages implemented in R software (version 4.0.2). Only dimension features containing the vast majority of variances were retained, with those containing zero variances discarded based on a principal component analysis (PCA) plot. Using weighted gene correlation network analysis (WGCNA) package in R language to perform WGCNA on data. WGCNA adds phenotypic weight parameters in the process of constructing gene co-expression network and optimizes classification using scale-free clustering and dynamic tree cut algorithm to achieve accurate and efficient data analysis.

2.3. Functional annotation of DEGs

Functional annotation of key genes was performed using the ClusterProfiler package in R and followed by comprehensive functional correlation of candidate genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to evaluate related functional categories. Significant differentially expressed genes (DEGs) were filtered, and bar charts and bubble charts were drawn according to their functions and pathways, respectively.

2.4. Construction of protein–protein interaction networks

Protein–protein interaction (PPI) networks of the DEGs were generated using the STRING online tool and Cytoscape software. Summarily, DEGs were entered into the STRING online tool to screen for interacting proteins with a combined score of > 0.9. The obtained PPI results were imported into Cytoscape software, then prominent modules in PPI were screened to construct the network through molecular complex detection. NetworkX package was then used to visualize the network in Python, with node sizes taken to be proportional to node degree.

2.5. Quantitative real-time polymerase chain reaction

Peripheral blood samples were collected from in- or out-patients at the Funan County Hospital of Traditional Chinese Medicine. The study group comprised 4 non-SONFH (following steroid administration for at least 3 months and 1 female and 3 males) and 8 SONFH (5 females and 3 males) patients, with samples collected between January 2020 and July 2021. Efforts were made to explain the study objectives to the subjects, who voluntarily signed a written informed consent provided by the Ethics Committee of Funan County Hospital of Traditional Chinese Medicine before sample collection. RNA was extracted from by mixing 1 mL of TRIzol containing RNase inhibitor with 100 μL of blood sample. RNA purity was measured and then reverse transcribed to complementary DNA. Primers were designed using the Primer-blast tool targeting some genes randomly selected in the PPI networks for both up- and downregulated genes. The complementary DNA was subjected to quantitative real-time polymerase chain reaction (qRT-PCR) on the Applied Biosystems 7500 RT-PCR System. Relative gene expression was determined using the 2-△△Ct method. Primers used in the study were synthesized by Shenzhen Zhicheng Medical Technology Co., Ltd and are outlined in Table 1.

Table 1

The primer sequences used for qRT-PCR validation.

Gene	Forward primers	Reverse primers
ITGAX	5′-GGGATGCCGCCAAAATTCTC-3′	5′-ATTGCATAGCGGATGATGCCT-3′
LILRB2	5′-GCATCTTGGATTACACGGATACG-3′	5′-CTGACAGCCATATCGCCCTG-3′
PPP3CB	5′-CCCCAACACATCGCTTGACAT-3′	5′-GGCAGCACCCTCATTGATAATTC-3′
MMP8	5′-TGCTCTTACTCCATGTGCAGA-3′	5′-TCCAGGTAGTCCTGAACAGTTT-3′

qRT-PCR = quantitative real-time polymerase chain reaction.

2.6. Statistical analysis

Statistical analyses were performed in GraphPad 8.0 software, and data presented as means ± standard deviations. Student t test and Mann–Whitney U were applied to compare data conforming to normal and nonnormal distribution, respectively. At least 3 independent experiments were used in the study, and data followed by P ≤ .05 were considered statistically significant.

3. Results

3.1. DEGs in SONFH

We analyzed transcriptome data for 40 samples (30 SONFH alongside 10 non-SONFH subjects) and then performed consistency analysis. Transcriptome data annotations and PCAs are shown in Figure 1A and B. The correlation within non-SONFH and SONFH groups was high, and higher than that between groups. PCA of the total samples under positive ion mode demonstrated that the transcriptome expression profile of the SONFH group changed by >49% compared to that of the non-SONFH group. The difference between transcriptomes in both groups was >9%. These data indicated high consistency between non-SONFH and SONFH groups, although differences were observed between groups. Cluster analysis results of the samples and the eigengene adjacency heatmap for the top 14 DEGs, respectively are shown in Figure 1C and D. Next, we analyzed differences in gene expression between the 2 groups, using a fold change >2 and a P value <.05 as the threshold. Results revealed a total of 584 DEGs, of which 294 and 290 were upregulated and downregulated, respectively. We also generated a volcano plot of these DEGs to show distribution of the genes between the 2 groups (Fig. 1E). The top 10 up- and downregulated genes in non-SONFH and SONFH groups and their functional annotation are listed in Tables 2 and 3, respectively.

Figure 1.

Identification of DEGs. (A, B) Two sets of transcriptome data annotations and PCAs. (C) The cluster dendrogram of the sample. (D) Heatmap showing clustering of the top 14 DEGs. (E) Distribution of statistically significant DEGs between the 2 groups. DEG = differentially expressed gene, PCA = principal components analysis.

Table 2

Top 10 genes with significantly different downregulated gene expression in the 2 groups.

Name	Log2 (fold change)	P value	Correction P value	Function
EIF1AY	–3.59866737	.00142	.00845	Stabilize the binding of Met-tRNA to the 40S ribosomal subunit
CA1	–3.5690905	2.52E–09	6.95E–07	Catalytic reversible hydration of carbon dioxide
BPGM	–3.3879619	4.35E–12	8.09E–09	Synthesis and degradation of 2,3-diphosphoglycerate
IGF2/ INS-IGF2	–3.3839209	8.15E–08	8.00E–06	Participate in development and growth
GYPA	–3.30603157	3.87E–12	7.95E–09	Coding for major sialoglycoproteins on the human erythrocyte membrane
RAP1GAP	–3.25195777	7.56E–11	5.32E–08	Down regulate the activity of the Ras-related RAP1 protein
RHCE/ RHD	–3.24127427	1.99E–09	5.74E–07	Determine Rh blood group
SNCA	–3.0973601	1.36E–09	4.53E–07	For integrating presynaptic signaling and membrane trafficking
IFIT1B	–2.9970138	5.50E–10	2.44E–07	May be involved in the defense response against bacteria and viruses
GYPB	–2.98665627	1.02E–08	2.00E–06	Coding for major sialoglycoproteins on the human erythrocyte membrane

Table 3

Top 10 genes with significantly different upregulated gene expression in the 2 groups.

Name	Log2 (fold change)	P value	Correction P value	Function
CCR3	2.26577837	1.48E–06	6.23E–05	Bind and respond to a variety of chemokines
CHI3L1	1.99439617	8.30E–05	.00109	Participate in inflammation and tissue remodeling
KRT23	1.95157603	3.49E–06	.000113	Responsible for the integrity of epithelial cell structure
LGALS2	1.89826127	7.99E–06	.000202	Encodes and regulates homodimer levels
FAM198B	1.83099127	8.41E–06	.00021	Tumor proliferation, invasion, metastasis, angiogenesis and receptor 13 activation
SHTN1	1.80011493	6.34E–06	.000172	Involved in positive regulation of neuron migration
CFD	1.76660107	1.20E–06	5.46E–05	Cleavage of catalytic factor B
HAL	1.7645862	4.52E–07	2.75E–05	Catalyzing the first reaction in histidine catabolism
SDPR	1.75829207	2.36E–05	.000441	May regulate the formation of caveolae and act as tumor suppressor
WLS	1.75482067	5.28E–08	6.13E–06	Involved in positive regulation of cell communication and protein transport

3.2. Functional annotation and pathway enrichment of the DEGs

GO and KEGG analyses revealed the biological function and enriched pathways associated with the DEGs in pathogenesis of SONFH. GO terms revealed that downregulated genes were mainly involved in red blood cell differentiation, cell protein catabolism, and gas transportation (Fig. 2A), whereas upregulated ones were involved in activation of myeloid leukocytes, phagocytosis, and inflammatory response (Fig. 2B). KEGG pathway analysis indicated that downregulated DEGs mainly regulated mitophagy-animal, human T-cell leukemia virus-1 (HTLV-1) infection, and Forkhead box O(FOXO) signaling pathways (Fig. 2C). On the other hand, upregulated DEGs were mainly involved in phagocytosis, osteoclast differentiation, and cytokine–cytokine receptor interaction pathways (Fig. 2D).

Figure 2.

Function enrichment analysis and construction of a PPI network. (A) Downregulated genes were enriched in red blood cell differentiation, cell protein catabolism, and gas transportation. (B) Upregulated genes were enriched in activation of myeloid leukocytes, phagocytosis, and inflammatory response. (C) Downregulated genes were involved in regulation of mitophagy-animal, HTLV-1 infection, and FOXO signaling pathways. (D) Upregulated genes regulated phagocytosis, osteoclast differentiation, and cytokine–cytokine receptor interaction pathways. (E) The top 2 MCODE components of downregulated PPIs. (F) The top 2 of MCODE components of upregulated PPIs. FOXO = Forkhead Box O, HTLV-1 = human T-cell leukemia virus-1, PPI = protein–protein interaction.

3.3. PPI network

We used the STRING online tool and Cytoscape software to construct a PPI network of the 584 DEGs. Next, we extracted the top 4 MCODE components from the PPIs and fund that the top 2 MCODE components for downregulated PPIs included hexokinase 1, protein phosphatase 3 catalytic subunit beta (PPP3CB), FK506 binding proteins 1B, peroxiredoxin 2 (PRDX2), bisphosphoglycerate mutase, DnaJ heat shock protein family (Hsp40) member A4, haptoglobin (HP), olfactomedin 4, and matrix metallopeptidase 8 (MMP8; Fig. 2E). The top 2 MCODE components for upregulated PPIs included formyl peptide receptor 2 (FPR2), integrin subunit alpha X (ITGAX), signal regulatory protein alpha, leukocyte immunoglobulin-like receptor B2 (LILRB2), and platelet-activating factor receptor (Fig. 2F).

3.4. qRT-PCR results

We randomly selected some genes, including ITGAX, LILRB2, PPP3CB, and MMP8, from the PPI network and validated their levels of expression using qRT-PCR. Results were consistent with profiles observed during differential gene expression analysis. As shown in Fig. 3, particularly, ITGAX and LILRB2 were significantly upregulated, whereas PPP3CB and MMP8 were significantly downregulated in peripheral blood of SONFH patients (P < .05; Fig. 3).

Figure 3.

qRT-PCR results showing significant upregulation of ITGAX and LILRB2 and significant downregulation of PPP3CB and MMP8 in SONFH patients., ITGAX = integrin subunit alpha X, LILRB2 = leukocyte immunoglobulin like receptor B2, MMP8 = matrix metallopeptidase 8, PPP3CB = protein phosphatase 3 catalytic subunit beta, qRT-PCR = quantitative real-time polymerase chain reaction, SONFH = steroid-induced osteonecrosis of the femoral head.

4. Discussion

More than 30 years ago, researchers discovered a close relationship between use of glucocorticoids with the occurrence of femoral head necrosis.[ Results from a retrospective analysis showed that hormone-related nontraumatic ONFH is dominant.[ Investigating the genetic factors regulating SONFH characteristics is an important strategy, and scholars have employed bioinformatics to explore sequence data. However, these studies also present some defects, such as inability to accurately interpret high-throughput sequence data to ensure that a certain deviation exists between PCR verification data and the data significance.[ In the present study, we explored the signaling pathways regulating development of early SONFH and predicted diagnostic biomarkers in peripheral blood from this group of patients.

Our results revealed that hundreds of coding genes, including tumor necrosis factor-α (TNF-α), phosphoinositide 3-kinase Akt (PI3K-Akt), hypoxia-inducible factor-1 (HIF-1), adherens junction, and mitogen-activated protein kinase signaling pathways, were significantly upregulated in peripheral blood of patients with early SONFH. Notably, these DEGs participate in regulation of redox, inflammatory response, angiogenesis, and cell metabolism by secreting inflammatory chemokine signaling proteins. We constructed a PPI network of these DEGs and found that PRDX2, HP, MMP8, FPR2, and ITGAX may be associated with the SONFH process. Among them, FPR2 is a low-affinity receptor of N-formyl methionine peptides and a powerful neutrophil chemokine.[ Family with sequence similarity 19 member A5 is a neural factor or brain-specific chemokine that was shown to remarkably downregulate expression of osteoclast-related genes, such as receptor activator of nuclear factor kappa and TNF receptor-associated factor 6, and mediate upregulation of osteoclast-producing negative regulators, such as musculoaponeurotic fibrosarcoma oncogene homolog B and interferon regulatory factor-8.[ Interestingly, FPR2 antagonist WRW4 or FPR2 deficiency were shown to markedly reverse inhibition of Family with sequence similarity 19 member A5-inducedosteoclastogenesis, indicating that FPR2 plays an important role in regulation of osteoclastogenesis as well as occurrence and development of SONFH disease.[

Integrin alpha-X/beta-2, which is encoded by ITGAX, is a receptor that recognizes the G protein-coupled receptor sequence in fibrinogen, thereby mediating cell–cell interactions during inflammatory responses. Functionally, it plays an important role in monocyte adhesion and chemotaxis.[ In the present study, we found that PRDX2 was a thiol-specific peroxidase in downregulated genes in SONFH relative to non-SONFH group. Functionally, this gene catalyzes reduction of hydrogen peroxide and organic peroxide to produce water and alcohol, thereby causing inactivation of toxic substances.[ Studies have shown that PRDX2, a hydrogen peroxide-mediated signal transduction sensor that prevents oxidative stress, not only plays a role in cell protection but may also be participating in signal cascade reaction of growth factor and TNF-α by regulating the intracellular concentration of hydrogen peroxide.[

HP-related diseases, which include haptoglobinemia and Plasmodium falciparum malaria,[ are regulated by the innate immune system and vesicle-mediated transport. HP accumulates in the kidney and is secreted in the urine due to hemolysis. Moreover, HP reportedly captures and binds to free plasma HP, thereby preventing kidney damage. HP also acts as an antioxidant, presents antibacterial activity, and plays a role in many acute phase regulations.[ Therefore, the observed downregulation of PRDX2 and HP may suggest that both genes promote development of SONFH disease through the redox pathway.

Matrix metalloproteinases (MMPs) are considered the key regulator in maintenance of bone homeostasis.[ Among them, MMP8, which is encoded by the MMP8 gene, also known as neutrophil or PMNL collagenase, is a collagen lyase that generally exists in mammalian connective tissues. The reduction of collagenase is the premise of effective osteoclast fusion. Inhibition of endogenous collagenase expression or activity was associated with a significant increase in osteoclast fusion, while fusion of trap-positive mononuclear osteoclast precursors is essential in the formation of fully functional mature osteoclasts.[ Moreover, MMP8 deficiency was associated with increased occurrence of arthritis and bone loss in the serum a metastatic arthritis model.[ Results from in situ hybridization of rat hindlimb tissue sections revealed that profiles of MMP8 expression were upregulated in osteoblast precursor cells, osteoblasts, and chondrocytes during embryonic development of the bone and cartilage.[ These findings suggest that MMP8 expression may be playing a role in occurrence and progression of SONFH by affecting the structure and function of osteoblasts and osteoclasts.

This study had several potential shortcomings. First, only the important abnormal gene groups were randomly selected and verified, due to limited resources. It is possible that the overall results may be some deviation from the actual results. Second, our sample size for qRT-PCR validation was small, as evidenced by only 12 volunteers, which may not provide strong statistical support. Finally, the authors believe that in subsequent studies, clear and effective interventions for SONFH can be used to observe the changes of peripheral blood biomarkers in patients with SONFH, so as to achieve indirect verification and make the research translate to clinical practice. Unfortunately, this study does not refine these works.

5. Conclusions

Transcriptomic analysis revealed upregulation of key genes that regulate important signaling pathways during SONFH development, including TNF-α, PI3K-Akt, HIF-1, adherens junction, and mitogen-activated protein kinase signaling pathways. Overall, these genes promote SONFH occurrence by regulating redox pathways, inflammatory reactions, and structure and function of osteoblasts and osteoclasts, suggesting that they may be important targets for designing approaches for early diagnosis and treatment of SONFH.

Authors’ contributions

Wenhui Shao and Honghai Zhou collected and analyzed data. Cong Tian wrote the manuscript. All authors read and approved the final manuscript.