Computational and Experimental Biology Reveals Dihydroartemisinin’s Efficacy Against Steroid-Induced Osteonecrosis of the Femoral Head Adjusting Ferroptosis via CCL17-PRDX6

Introduction

Steroid-induced osteonecrosis of the femoral head (SONFH) is a significant complication that arises due to prolonged use of glucocorticoids.¹ This disorder is characterized by the progressive worsening of clinical symptoms. The initial symptoms of this disorder primarily manifest as insidious, deep, dull pain in the hip or groin region.^2,3 In addition, the physical examination yielded positive results for the hip internal rotation and abduction test, the four-sign test, and an abnormal gait pattern.⁴ The principal patient population affected by SONFH is comprised of young adults between the ages of 35 to 45, exhibiting a notable predominance within this age bracket.⁵ Nevertheless, the prevention and treatment of this condition have become a pressing challenge in the field of orthopedics. On one hand, epidemiological data from the Chinese Femoral Head Necrosis Database indicate that patients with a complete history of glucocorticoid use account for 26.84% of osteonecrosis of the femoral head cases.⁶ A more thorough investigation reveals that autoimmune diseases, skin diseases, and respiratory diseases continue to be the primary underlying causes of SONFH.³ K-H Koo et al conducted a study that followed 22 patients with a complete history of glucocorticoid use across four large hospitals.⁷ The study found that 21 of the patients developed osteonecrosis of the femoral head with subsequent collapse within 12 months. Michael A. Mont et al also observed a positive correlation between the use of high-dose glucocorticoids and the occurrence of SONFH.⁸ It has been demonstrated that an escalation of 10 milligrams in the quotidian dosage of metformin results in an augmentation of the incidence rate of osteonecrosis to 3.6%. Conversely, approximately 1% to 3% of the global population receive systemic glucocorticoid therapy, and 0.5% to 1.8% of the population uses oral glucocorticoids long-term.⁹

Epidemiological data reflect the grim situation facing the prevention and treatment of SONFH. Despite the evident correlation between SONFH and glucocorticoid use, as established by existing studies, the insidious nature of early symptoms often results in disease progression for the majority of patients. Furthermore, approximately 60% of patients with asymptomatic osteonecrosis of the femoral head who do not receive timely and effective treatment may miss the optimal window for hip-preserving therapy, ultimately leading to femoral head collapse on imaging studies.¹⁰ Moreover, due to the multifaceted pathogenesis of osteonecrosis of the femoral head, 70% to 80% of patients with clear risk factors, such as glucocorticoid use, trauma, or alcohol abuse, will develop osteonecrosis of the femoral head within 1 to 3 years.¹¹ In light of the limited efficacy of hip-preserving surgical interventions, such as core decompression, Will Jiang et al conducted a comprehensive analysis of 3025 patients who underwent core decompression.¹² Their analysis, which drew upon the PearlDiver M157 administrative database, revealed that within a five-year period, 12.8% of these patients required total hip arthroplasty. Consequently, this may necessitate a greater number of revision surgeries, potentially requiring two to three procedures in the future.¹³ These surgeries could result in the permanent loss of labor capacity and a significant financial and emotional burden on patients. Therefore, the development of efficacious treatment methods that would delay the progression of SONFH is imperative.

The prevailing conservative treatment modalities encompass immobilization therapy and pharmacotherapy with phosphate binders, lipid-lowering agents, and vasodilators. Nevertheless, contemporary conservative treatment modalities demonstrate notable limitations in addressing SONFH. Non-weight-bearing therapy has been demonstrated to alleviate pain in early-stage SONFH patients. Despite this, the majority of studies have found that it has limited efficacy in improving the Harris Hip Score. Statistical analysis indicates that approximately 22.7% of patients can effectively delay disease progression after non-weight-bearing therapy.¹⁴ Phosphate drugs are effective for the treatment of osteoporosis. Nonetheless, due to their excessive inhibition of bone resorption, studies have found that they may induce drug-induced osteonecrosis of the jaw, with an incidence rate ranging from 1% to 9%.¹⁵ It is estimated that statins are one of the most widely used pharmaceuticals globally, and they are frequently prescribed as initial treatment for hypercholesterolemia.¹⁶ Conversely, Massimiliano Ruscica et al have posited that the prolonged utilization of statins may result in the clinical manifestations, including muscle discomfort and elevated transaminases.¹⁷ Kai-yun Chen et al found that when statins were used as a monotherapy to treat SONFH, trabecular bone thinning and increased porosity were observed by 9 weeks.¹⁸ This suggests that the protective effect of statins may weaken over time. Angiotensin-converting enzyme inhibitors have been proven to be effective in promoting vasodilation and inhibiting vasoconstriction.¹⁹ Although these medications are generally well-tolerated, common side effects include angioedema and persistent dry cough. Hence, there is a pressing need to identify medications that exhibit both efficacy and reduced toxicity.

As a novel form of programmed cell death, ferroptosis is characterized by abnormal iron accumulation and elevated lipid reactive oxygen species (ROS).²⁰ Prolonged excessive glucocorticoid use substantially increases intracellular ROS levels, exacerbating oxidative stress imbalance in the bone microenvironment.²¹ This aberrantly activates the ferroptosis pathway, directly inhibiting osteoblast proliferation and differentiation, disrupting bone metabolism, and accelerating the onset of SONFH.²² Notably, Huihui Xu et al confirmed via clinical and animal studies that osteoblast apoptosis in SONFH is closely linked to activation of the NOX/ROS/NF-κB pathway.²³ This finding is particularly significant because excessive ROS accumulation triggers both ferroptosis and inflammatory injury mediated by NF-κB, two intertwined pathological events in SONFH. Concurrently, dihydroartemisinin (DHA) has demonstrated considerable potential for application in bone repair as a semi-synthetic derivative of artemisinin.^24,25 Mechanistically, DHA acts on normal cells by upregulating glutathione peroxidase 4 (GPX4) levels, which reduces ROS and inhibits ferroptosis.²⁶ Lin Zhou et al found that DHA inhibits osteoclast formation in a dose-dependent manner by suppressing NF-κB and activating T-cell nuclear factor.²⁷ Furthermore, the study revealed that the administration of DHA led to a decrease in the mRNA levels of cathepsin K, calcitonin receptor, and tartrate-resistant acid phosphatase in the treated group. The observed decrease in these mRNA levels contributed to a reduction in the progression of osteoporosis, a condition triggered by estrogen deficiency. Consistent with this, Dong Ding et al confirmed via in vitro and in vivo osteoarthritis models that DHA inhibits osteoclast formation through suppressing NF-κB, MAPK, and NFATc1 expression.²⁸ Moreover, the study revealed that while DHA suppresses bone resorption, it does not influence the osteogenic differentiation and proliferation of bone marrow mesenchymal stem cells. Consequently, the selection of DHA for the treatment of SONFH is a judicious choice. However, the current research mechanism remains opaque, and the core genes have yet to be identified. Mendelian randomization (MR) analysis employs genetic variation as an instrumental variable to infer causal relationships, thereby effectively addressing confounding factors in observational studies.^29,30 MR directly utilizes publicly available genome-wide association study (GWAS) summary data, eliminating the need to re-collect samples and thereby significantly reducing research costs and time.³¹ MR has been identified as a pivotal instrument in the validation of drug targets and the decision-making process for precision medicine.^32,33 The confirmation of the causal association between changes in gene exposure and disease risk enables MR to make prospective predictions regarding drug efficacy and side effects, accelerate target screening, and provide clues for the discovery of new biomarkers.³⁴ Therefore, our study will use multi-omics to investigate the core mechanisms through which DHA exerts its therapeutic effects on SONFH by targeting ferroptosis (Figure 1). This figure was created using BioGDP.³⁵

Figure 1 Flowchart of multi-omics study on the mechanism of action of DHA in regulating ferroptosis to treat SONFH.

Materials and Methods

Transcriptome Sequencing and Enrichment Analysis

Experimental Animals

This study used 8-week-old healthy adult male Sprague-Dawley (SD) rats weighing 280±20 g. All animal experimental protocols were approved by the Animal Experimental Ethics Committee of Tianjin Hospital, Tianjin University (Approval No. 2025YLS145). All procedures were performed in strict accordance with the guidelines for laboratory animal ethics and welfare. The rats were housed in a quiet environment with a controlled temperature between 21 and 24 °C and a 12-hour light/dark cycle. The cages and water dispensers were regularly cleaned and disinfected, and fresh feed and drinking water were provided daily. All surgeries and specimen collections were performed under isoflurane anesthesia to minimize pain and stress. The number of animals used during the experiment was reasonably controlled in compliance with animal welfare requirements.

Experimental Groups

The SD rats were randomly divided into control group (n = 3), model group (n = 3), and DHA group (n = 3). To induce SONFH, the DHA and model groups received a single intravenous injection of LPS (0.2 mg/kg, Solarbio) during the first week. Then, they received intraperitoneal injections of methylprednisolone sodium succinate (100 mg/kg, Pfizer) once daily for three days, followed by 40 mg/kg three times weekly for three weeks.³⁶ After successful modeling, the DHA group received intraperitoneal injectionsof DHA (13 mg/kg, Macklin) four times weekly for four weeks. The control group and the model group received equivalent doses of saline. After the treatment period, all surviving rats were euthanized and their tissues collected in accordance with institutional ethical guidelines.

Transcriptome Sequencing

Tissue samples from the femoral head were obtained from the control group, the model group, and the DHA group, and then subjected to transcriptome sequencing. To obtain total RNA, the TRIzol reagent was utilized according to the manufacturer’s guidelines. In order to determine the purity and quantity of RNA, a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, USA) was utilized. The integrity of RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Transcriptome libraries were constructed using the VAHTS Universal V10 RNA-seq Library Prep Kit (Premixed Version) following the manufacturer’s instructions. Subsequently, the samples were loaded onto an Illumina sequencing platform for high-throughput sequencing, which generated 150 base pair paired-end reads. The raw reads in fastq format were processed using fastp software, and low-quality reads were removed to obtain clean reads for subsequent data analysis. The sequencing and analysis of the transcriptomes were carried out by Shanghai Ouyi Biotechnology Co., Ltd. (Shanghai, China).

Differential Gene Expression Analysis and Enrichment Analysis

HISAT2 software was used for reference genome alignment, and gene expression levels were calculated. HTSeq-count was utilized to obtain the read counts for each gene. R (v3.2.0) was selected for the execution of PCA analysis and the graphical representation of genes, with the objective of evaluating biological replicates of samples.

Differentially expressed genes (DEGs) were identified by DESeq2, with the screening criteria set at q-value < 0.05 and absolute fold change (|FC|) ≥ 2.³⁷ Additionally, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for the DEGs to screen for significantly enriched functional terms and signaling pathways. The R (v3.2.0) software was utilized to generate bar charts, chord diagrams, and enrichment analysis circle diagrams, which were employed to visualize significantly enriched functional entries. Gene Set Enrichment Analysis (GSEA) was performed to identify coordinated expression alterations of gene sets associated with biological pathways, based on genome-wide gene expression profiles.

MR Analysis

Study Design

A two-sample MR method was employed to investigate the impact of protein quantitative trait loci (pQTLs) associated with ferroptosis on SONFH. The present study adhered to the fundamental assumptions of MR analysis. Firstly, genetic variation is closely related only to the exposure factor under study. Furthermore, genetic variation is unrelated to any confounding factors other than the exposure factor. Additionally, the results are not directly influenced by any factors other than the exposure factor.

The present study was conducted in three parts. A two-way MR study was conducted to ascertain the causal relationship between pQTLs data from participants in the UK Biobank, serving as the exposure factor, and SONFH, functioning as the outcome, thereby deriving the total effect value. To this end, we employed ferroptosis-related pQTLs data as the exposure factor and SONFH as the outcome, thereby exploring the causal relationship between them and screening for positive mediators. The mediating effect of ferroptosis-positive pQTLs was calculated.

Data Sources

All data utilized in this study were obtained from online databases that institutional ethics committees have approved. The UK Biobank Proteomics Public Project (UKBPPP) has developed a comprehensive pQTLs map for 2923 proteins, which includes plasma proteomics data from 54,219 UK Biobank participants. The generation of these data was facilitated by the antibody-based OlinkExplore3072 platform. The 4907 plasma protein data were obtained from 35,559 European blood samples, which detected 27.2 million sequence variants associated with plasma proteins and identified 18,084 pre-protein quantitative trait locus associations. The core targets of ferroptosis were retrieved using FerrDBV2 (http://www.zhounan.org/ferrdb/current/). The data that satisfied the specified criteria were then intersected with the established key targets of ferroptosis. Outcome data for SONFH were obtained from release 11 of the FinnGen database (finngen_R11_OSTEON_DRUGS). Cases were defined as individuals with glucocorticoid-related osteonecrosis of the femoral head, based on Finnish nationwide inpatient and outpatient healthcare registers. Controls were defined as individuals without relevant ICD-10 codes for osteonecrosis, no history of osteonecrosis, and no long-term glucocorticoid use. The final pooled dataset comprised 218,792 individuals of European ancestry, including 101 SONFH cases and 218,691 controls.

Selection of Instrumental Variables

Instrumental variables were selected based on the following criteria. The genome-wide significance threshold for screening single nucleotide polymorphisms was set at p < 5e-8. The screening threshold for SONFH in the Finnish database was set at p < 5e-8. To reduce data heterogeneity and reflect more realistic causal relationships, linkage disequilibrium in exposure data was removed by setting kb = 10,000 and r² = 0.1. F > 10 was set to obtain instrumental variables strongly associated with exposure while minimizing bias.

Data Analysis

To correct for potential biases in genetic variation and ensure the robustness of the results, we employed inverse variance weighting as the primary testing method. As auxiliary testing methods, MREgger, weighted median, simple mode, and weighted mode were used, with the requirement that the directions of the five test results be consistent. The aforementioned primary method corresponds to p < 0.05 to ensure the significance of the results. The false discovery rate correction parameter was set to 0.2 to avoid false positives due to multiple testing. Cochran’s Q test was used to assess the heterogeneity of instrumental variables, evaluating their consistency across different conditions or samples. Based on the Cochran’s Q test, we failed to reject the null hypothesis of no heterogeneity (p > 0.05), indicating that the observed data did not provide sufficient evidence for the presence of heterogeneity. MR-Egger intercept analysis was used to detect horizontal pleiotropy at the data level. If the analysis results were not zero and p < 0.05, it indicated that the instrumental variables might directly influence the outcome factors through other pathways, resulting in horizontal pleiotropy. MR-Pleiotropy RESidual Sum and Outlier (MR-PRESSO) analysis identified studies with instrumental variables exhibiting abnormal behavior. Subsequently, the anomalous instrumental variables should be removed and reanalyzed. By sequentially eliminating instrumental variables, leave-one-out sensitivity analysis can identify data with excessive influence on causal relationships, thereby enhancing the reliability of the results. Finally, in conjunction with reverse MR studies, observe whether there is a significant causal relationship between the exposure data and the outcome data. If no significant relationship is found in the analysis, it implies a robust causal relationship between them. β1 represents the effect value of UKBPPP-related pQTLs on ferroptosis-related pQTLs. β2 represents the effect value of ferroptosis-related pQTLs on SONFH. β represents the effect value of UKBPPP-related pQTLs on SONFH. The mediating effect is calculated using the established formula: β1 × β2, and the mediated proportion is (β1 × β2)/β.

Computer-Aided Drug Design

The three-dimensional structure models of all core target proteins were retrieved from the AlphaFold Database (AlphaFold DB, https://alphafold.ebi.ac.uk/), with a strict screening threshold of overall predicted local distance difference test score ≥ 70 to ensure the high confidence and reliability of the models. The UniProt database (https://www.uniprot.org/) was only used to verify the consistency between the amino acid sequence derived from the downloaded PDB-formatted structure files and the canonical sequence of the target protein, to ensure the accuracy of the protein structure for subsequent analysis. Standardized protein preprocessing was performed sequentially: solvents, non-specific ligands and heteroatoms were removed via PyMOL (v3.1.1), followed by polar hydrogen addition, Gasteiger charge calculation and PDBQT format conversion using AutoDockTools (v1.5.7). The SDF structure of active ingredient was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The SDF structure was converted to an mol2 structure using OpenBabelGUI software (version 3.1.1). The protein and molecule underwent conversion to pdbqt format using AutoDockTools software (version 1.5.7). A grid parameter file was generated using AutoDockTools software. AutoDock vina software was utilized to perform docking calculations and identify the most stable binding mode between the ligand and receptor.³⁸ A vina score close to −10 kcal/mol or lower is generally indicative of favorable protein-molecule binding. When the vina score approaches −5 kcal/mol, it signifies that the protein binds to the molecule at an average level. The GraphPad Prism software (version 9.5) was utilized to generate a heatmap illustrating the molecular docking scores. Furthermore, the visualization of molecular docking results were visualized using PyMOL software (version 3.1.1).

The subsequent analysis of the docking of DHA with the four targets selected for this study was based on the results of MR analysis and molecular docking. Molecular dynamics simulations were performed using Gromacs 2022. The force field parameters were obtained using the pdb2gmx tool in Gromacs and the AutoFF website. During the simulation, the molecular parameters of the receptor protein were based on the CHARMM36 force field, while those of the ligand were based on the CGenff force field. To this end, a 1 nm TIP3P-type cubic water box was added around the system for solvent simulation. Ions were incorporated into the system through the utilization of the gmx genion tool to attain electrical neutrality. The study employed the Particle Mesh Ewald (PME) method to facilitate the management of long-range electrostatic interactions, establishing a cutoff distance of 1 nm. The SHAKE algorithm was employed to constrain all bonds, and the molecular dynamics simulation was executed using the Verlet leapfrog algorithm with an integration step size of 1 fs. Prior to the molecular dynamics simulation, the system underwent an energy optimization process. The energy minimization process comprised 3000 steps of steepest descent optimization, followed by 2000 steps of conjugate gradient optimization. Initially, the solute was constrained, and the water molecules were minimized. Subsequently, the counterions were constrained, and the system was minimized. Consequently, the entire system was minimized without constraints. The simulation was executed at a temperature of 310 K under constant pressure in an NPT system, with a simulation time of 100 ns.³⁹ During the simulation, the g-RMSD, g-RMSF, g-HBonds, g-Rg, and g-SASA tools were used to calculate the root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), hydrogen bonds (HBonds), radius of gyration (Rg), and solvent-accessible surface area (SASA), respectively. The gmx_mmpbsa online tool was utilized to calculate the binding free energy between proteins and molecules based on the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) method. The binding free energy of the complex was then analyzed.

Results

Transcriptome Sequencing and Enrichment Analysis

Principal component analysis (PCA) was used to evaluate intragroup reproducibility and intergroup variability. As shown in Figure 2a, PC1 accounted for 50.37% and PC2 accounted for 15.53%, indicating comprehensive data coverage. Additionally, samples clustered within each group, indicating good reproducibility. Intergroup dispersion indicated intergroup differences, facilitating subsequent differential gene screening. A total of 125 upregulated genes and 270 downregulated genes were identified between the model group and the control group (Figure 2b). A total of 176 upregulated genes and 86 downregulated genes were identified between the DHA group and the model group (Figure 2c). The results of the GO enrichment analysis were categorized into three categories: biological process (BP), molecular function (MF), and cellular component (CC). The significant enrichment analysis results for differentially expressed genes between the model and control groups included 291 BP, 53 MF, and 42 CC (Figure 2d). The significant enrichment analysis results for differentially expressed genes between the DHA and model groups included 183 BP, 40 MF, and 32 CC (Figure 2e). Notably, in the BP group, the three results of B cell proliferation, B cell receptor signaling pathway, and B cell activation showed a significant downregulation trend in the model vs control group. However, these three results showed a significant upregulation trend in the DHA vs model group. In the MF group, the three results of transmembrane signaling receptor activity, DNA-binding transcription factor activity, and phosphorus-oxygen lyase activity showed a significant downregulation trend in the model vs control group. Conversely, these three results showed a significant upregulation trend in the DHA vs model group. Additionally, in the CC group, the three results of external side of plasma membrane, B cell receptor complex, and plasma membrane showed a significant downregulation trend in the model vs control group. Conversely, these three results showed a significant upregulation trend in the DHA vs model group.

Figure 2 Transcriptome sequencing and enrichment analysis result charts. (a) PCA diagram. (b and c) Volcano plots of DEGs. (d and e) Bar charts of GO enrichment analysis. (f and g) Bubble plots of KEGG enrichment analysis. (h and i) Enrichment score line charts for GSEA.

KEGG enrichment analysis revealed that the significantly enriched pathways in the model vs control group included 31 pathways (Figure 2f). The significantly enriched pathways in the DHA vs model group included 23 pathways (Figure 2g). Analysis of the above results revealed that the expression of the three results—primary immunodeficiency, B cell receptor signaling pathway, and non-homologous end-joining—showed a significant downregulation trend in the model vs control group. However, the expression of these three results showed a significant upregulation trend in the DHA vs model group.

Additionally, to detect activated high expression pathways and avoid false negatives. KEGG was selected as the gene set source for GSEA. The significant enrichment analysis results of differentially expressed genes between the model and control groups included 327 pathways. The significant enrichment analysis results of differentially expressed genes between the DHA and model groups included 326 pathways. Further analysis of the above results revealed that the ferroptosis signaling pathway showed a significant downregulation trend in the model vs control group (Figure 2h), while this pathway showed a significant upregulation trend in the DHA vs model group (Figure 2i). To further verify the regulatory effect of DHA on the ferroptosis pathway in SONFH, we analyzed the differential expression of core ferroptosis regulatory genes between groups (Table S1-2). Compared with the control group, the SONFH model group showed significantly downregulated expression of the ferroptosis suppressor genes GPX4, solute carrier family 7 member 11 (SLC7A11) and ferritin heavy chain 1(FTH1), and significantly upregulated expression of the ferroptosis promoter genes acyl-CoA synthetase long-chain family member 4(ACSL4) and transferrin receptor 1(TFRC) (p < 0.05). After DHA intervention, the expression of GPX4, SLC7A11 and FTH1 was significantly upregulated, while the expression of ACSL4 and TFRC was significantly downregulated (p < 0.05). These genes are well-recognized molecular markers of ferroptosis, and their consistent expression changes, combined with the GSEA results of the ferroptosis signaling pathway, provide solid molecular-level evidence for the regulatory effect of DHA on ferroptosis in SONFH. These results indicate that DHA can effectively reverse the abnormal ferroptosis-related gene expression in SONFH, thereby delaying the progression of the disease.

MR Analysis

Screening Instrumental Variables

A total of 49,251 SNPs related to UKBPPP data were obtained according to the screening criteria. For the final included SNPs associated with UKBPPP plasma protein pQTLs exposure, the F-statistics ranged from 29.715 to 24649.71, with a mean value of 229.54 (Table S3). The core targets of ferroptosis were deduplicated and merged to obtain a total of 564 targets. The pQTLs data of 4907 plasma proteins were intersected with the core targets of ferroptosis to obtain a total of 8218 SNPs. For SNPs associated with ferroptosis-related core genes exposure, the F-statistics ranged from 29.720 to 13642.73, with a mean value of 124.80 (Table S4). All F-statistics were well above the 10 threshold, confirming the absence of weak instrumental variables in the analysis.

Causal Relationship Between UKBPPP and SONFH

According to the analysis method and screening criteria, a total of 110 targets with a strong causal relationship with SONFH were obtained (Table S5), including serine protease inhibitor Kazal-type 4 (SPINK4), serine protease inhibitor Kazal-type 2 (SPINK2), zona pellucida sperm-binding protein 3 (ZP3), apolipoprotein E (APOE), protein cholesin (C7ORF50), and C-C motif chemokine 17 (CCL17), among others (Figure 3a). The detailed SNP-level analysis results for the causal relationship between UKBPPP and SONFH were shown in Table S6. Heterogeneity tests, sensitivity analyses, and horizontal pleiotropy tests were conducted on these targets simultaneously (Table S7-8). The analysis results showed no abnormalities, confirming the robustness of the main MR results. The Cochran’s Q test yielded p > 0.05, indicating no heterogeneity in the data. MR-Egger and MR-PRESSO tests did not identify any horizontal pleiotropy in the analysis results.

Figure 3 Forest diagrams for MR analysis. (a) Causal relationship between UKBPPP key genes and SONFH. (b) Causal relationship between ferroptosis core genes and SONFH. (c) Causal relationship between UKBPPP key genes and ferroptosis core genes. Bold p values indicate significant associations (p < 0.05). Abbreviations: OR, odds ratio; CI, confidence interval; nsnp, number of SNPs.

To eliminate reverse causation bias, we performed a reverse MR analysis, considering SONFH as the exposure and UKBPPP key genes as the outcomes. We rigorously quality-controlled the exposure and outcome data. The instrumental variables had an F-statistic greater than 10. Additionally, we excluded linkage disequilibrium. SNP locus matching revealed no overlap between the SONFH genetic instrumental variable SNPs and the expression-regulating SNPs of the UKBPPP key genes (Table S9). These results suggest that genetic susceptibility to SONFH may not influence the expression of UKBPPP key genes at the genomic level. Thus, a reverse causal relationship is precluded from a genetic mechanistic perspective.

Causal Relationship Between Ferroptosis and SONFH

Based on the analysis methods and screening criteria, two targets with strong causal relationships with SONFH were identified (Table S10), including Ferritin mitochondrial (FTMT) and Peroxiredoxin-6 (PRDX6) (Figure 3b). The detailed SNP-level analysis results for the causal relationship between ferroptosis and SONFH were shown in Table S11. After heterogeneity testing, sensitivity analysis, and horizontal pleiotropy testing, no abnormalities were observed in the above results, confirming the robustness of the main MR results (Table S12-13). The Cochran’s Q test result had p > 0.05, indicating no heterogeneity in the data. The MR-Egger and MR-PRESSO tests did not detect any multi-effect in the analysis results.

Causal Relationship Analysis Between UKBPPP and Ferroptosis

Using the aforementioned analytical methods and screening criteria, we investigated the causal relationships between UKBPPP key genes and ferroptosis core genes (Figure 3c). These were used as exposure data to examine their causal relationship with FTMT and PRDX6 (Table S14). SPINK4 and ZP3 showed strong correlations with FTMT. SPINK2, ZP3, APOE, C7ORF50, and CCL17 showed strong correlations with PRDX6. The detailed SNP-level analysis results for the causal relationship between UKBPPP key genes and ferroptosis core genes were shown in Table S15. To ensure the robustness of the analysis results and avoid false-positive outcomes, we simultaneously conducted heterogeneity tests, sensitivity analyses, and horizontal pleiotropy tests on these genes (Table S16-17). The analysis results showed no abnormalities, confirming the robustness of the main MR results. The Cochran’s Q test yielded p > 0.05, indicating no heterogeneity in the data. The MR-Egger and MR-PRESSO tests did not identify any multi-effect issues in the analysis results.

Mediator Effect Analysis

Combining the results of the above two steps of MR analysis, a total of seven signaling pathways were obtained. Among these seven pathways, the elevation of core targets in UKBPPP leads to a decrease in FTMT and PRDX6, thereby accelerating ferroptosis and aggravating the occurrence of SONFH. Among them, CCL17 regulates ferroptosis and has the most obvious effect in mediating the occurrence of SONFH. With the increase of CCL17, the expression of PRDX6 will show a downward trend. With the decrease of PRDX6, the occurrence of ferroptosis in SONFH patients will be further aggravated, thereby accelerating the disease process (Table 1).

Table 1 Mediated Pathways Regulating Ferroptosis Between UKBPPP and SONFH

Computer-Aided Drug Design

The specific AlphaFold DB structure IDs of the core targets are as follows: CCL17 (AF-Q92583-F1), PRDX6 (AF-P30041-F1), C7ORF50 (AF-Q9BRJ6-F1), APOE (AF-P02649-F1), FTMT (AF-Q8N4E7-F1), SPINK4 (AF-O60575-F1), ZP3 (AF-P21754-F1), SPINK2 (AF-P20155-F1). In our study, all eight core targets exhibited Vina scores below −5 kcal/mol (Figure 4). Among them, C7ORF50 had the best docking result, with a vina score of −6.9 kcal/mol. But the worst docking result was for SPINK2, with a vina score of −5.5 kcal/mol. At the same time, based on the MR analysis results and molecular docking scores from Vina, the docking results of C7ORF50, APOE, PRDX6, and CCL17 with DHA were selected for further analysis.

Figure 4 Schematic diagrams of the core targets and DHA molecule docking. (a) C7ORF50-DHA (−6.9 kcal/mol). (b) APOE-DHA (−6.7 kcal/mol). (c) CCL17-DHA (−5.8 kcal/mol). (d) FTMT-DHA (−6.6 kcal/mol). (e) PRDX6-DHA (−6.6 kcal/mol). (f) SPINK4-DHA (−6.8 kcal/mol). (g) ZP3-DHA (−5.9 kcal/mol). (h) SPINK2-DHA (−5.5 kcal/mol).

Root-mean-square deviation (RMSD) is a good indicator of the conformational stability of proteins and ligands, as well as a measure of the degree of deviation of atomic positions from their initial positions. The smaller the deviation, the better the conformational stability. Therefore, RMSD was used to assess the equilibrium of the simulated system. As shown in Figure 5a, the PRDX6-DHA complex system reached equilibrium after 5 ns and ultimately fluctuated around 2.3 Å. The C7ORF50-DHA complex system reached equilibrium after 30 ns, with final fluctuations around 12.1 Å. The APOE-DHA complex system reached equilibrium after 35 ns, with final fluctuations around 1.6 Å. The CCL17-DHA complex system exhibited no significant fluctuations. Among these, CCL17-DHA complex RMSD values were relatively low. Therefore, the small molecule DHA exhibits high stability when binding to the CCL17 target protein. The radius of gyration (Rg) can be used to describe changes in the overall structure and characterize the compactness of protein structures. A larger Rg change indicates greater system expansion. The CCL17-DHA, PRDX6-DHA, and APOE-DHA complexes exhibit relatively stable fluctuations during motion. This indicates that the small molecule-target protein complexes do not undergo significant expansion or contraction during motion. The C7ORF50-DHA complex system exhibits fluctuating changes during motion. This indicates that the small molecule-target protein complexes undergo conformational changes during motion (Figure 5b). Solvent-accessible surface area (SASA) is an indicator for assessing protein surface area. In this simulation, the solvent-accessible surface area between the target protein and the small molecule was calculated (Figure 5c). The results show that the PRDX6-DHA, APOE-DHA, C7ORF50-DHA, and CCL17-DHA complex systems exhibit slight fluctuations and gradually stabilize. This demonstrates that binding to small molecules affects the binding microenvironment and leads to changes in SASA to some extent. Hydrogen bonds play a crucial role in the binding of ligands to proteins. The number of hydrogen bonds between small molecules and target proteins during the kinetic process is shown in Figure 5d. The number of hydrogen bonds between small molecules and target proteins ranges from 0 to 10, with most complexes having approximately 2 hydrogen bonds. This indicates that the ligand exhibits good hydrogen bonding interactions with the target protein. The root mean square fluctuation (RMSF) can indicate the flexibility of amino acid residues in a protein. As shown in Figure 5e, the RMSF values of the PRDX6-DHA, APOE-DHA, C7ORF50-DHA, and CCL17-DHA complexes are relatively low (mostly below 6 Å), indicating lower flexibility and higher stability.

Figure 5 Analytical diagrams of molecular dynamics simulations. (a) RMSD analysis. (b) Rg analysis. (c) SASA analysis. (d) Hydrogen bond diagram. (e) RMSF analysis.

Subsequently, using the binding conformation of the complexes, the binding free energy between the small molecules and the target proteins was calculated using the MM/PBSA method (Figure 6a). The binding free energies of the C7ORF50-DHA, PRDX6-DHA, APOE-DHA, and CCL17-DHA complexes were −17.270 kcal/mol, −4.8 kcal/mol, −0.01 kcal/mol, and −27.75kcal/mol, respectively. Negative values indicate that the molecule has binding affinity for the target protein, with lower values indicating stronger binding. Based on the above analysis, both target proteins exhibit strong molecular binding capabilities. Among them, CCL17-DHA demonstrates the highest binding affinity and will be the focus of our subsequent research.

Figure 6 Molecular binding free energy calculation and key amino acid residue interaction analysis. (a) Analysis of binding free energy components and total binding free energy for complexes of four core targets with DHA. (b) Binding energy decomposition analysis of key residues in the CCL17-DHA complex. (c) Binding energy decomposition analysis of key residues in the PRDX6-DHA complex. (d) Binding energy decomposition analysis of key residues in the APOE-DHA complex. (e) Binding energy decomposition analysis of key residues in the C7ORF50-DHA complex. Abbreviations: VDWAALS, van der Waals interaction energy; EEL, electrostatic interaction energy; EPB, polar solvation energy; ENPOLAR, non-polar solvation energy; GGAS, gas phase total free energy; GSOLV, solvation free energy; TOTAL, total binding free energy of the complex.

CCL17-DHA had the lowest total binding free energy (Figure 6b), 1.6-fold lower than C7ORF50-DHA, reflecting the strongest spontaneous binding ability. PRDX6-DHA maintained a stable negative binding free energy (Figure 6c), and its binding residues (PHE14, ASN17, ARG22) were located in the peroxidase active domain, ensuring functional relevance. In contrast, APOE-DHA had almost no binding affinity (Figure 6d), and C7ORF50-DHA’s high binding energy was offset by poor conformational stability (Figure 6e). The combination of MD simulation and MM/PBSA results confirmed that CCL17 and PRDX6 maintained stable binding conformations and effective intermolecular interactions throughout the simulation, while C7ORF50’s conformational fluctuation and APOE’s weak binding rendered them non-functional targets. Consequently, CCL17 and PRDX6 may serve as optimal core functional targets for DHA in the treatment of SONFH.

Discussion

Ferroptosis is a novel type of programmed cell death, which is characterized by excessive accumulation of intracellular iron and lipid peroxides.^40,41 In recent years, many studies have found that ferroptosis is closely related to bone metabolism balance, and excessive oxidative stress may lead to osteoblast damage and dysfunction.⁴² Yike Wang et al found through human bone tissue proteomics and in vitro and in vivo experiments that reduced expression of the key protein COPB1 will exacerbate ferroptosis, thereby accelerating bone loss.⁴³ Shuyuan Li et al confirmed that NRF2 can inhibit endogenous ferroptosis in osteoblasts.⁴⁴ The study further experimentally verified that targeting NRF2 and promoting the activation of the NRF2/ARE signaling pathway are conducive to inducing the growth and mineralization of intramembranous bone grafts. Through transcriptome sequencing, we found that the DHA intervention rat femoral heads had low expression and high ranking of ferroptosis signaling pathways compared with the model group. Therefore, we conducted research on the role of DHA in regulating ferroptosis in SONFH (Figure 7). This figure was created using BioGDP.³⁵

Figure 7 Schematic diagram of multi-omics combined investigation of DHA regulation of ferroptosis treatment for SONFH.

To further validate the causal relationship between the ferroptosis-related targets and SONFH, we performed two-sample MR analysis based on large-scale population genetic data. The UKBPPP database provides the largest, most comprehensively annotated, and strictly quality-controlled plasma proteomic pQTLs dataset currently available, which effectively avoids population stratification bias within the same ancestry group and meets the core assumptions of two-sample MR analysis, ensuring the reliability of our causal inference results.

Among the ferroptosis-related targets identified by our MR analysis, PRDX6 emerged as a core mediator in the pathological process of SONFH. According to previous studies, the cysteine residue within the PRDX6 structure serves as the catalytic center responsible for scavenging hydrogen peroxide. This inhibits ferroptosis and protects cells from damage. However, our follow-up MR study found that high expression of CCL17, ZP3, C7ORF50, APOE, and SPINK2 promotes PRDX6 upregulation. Ultimately, this exacerbates the occurrence of SONFH. Combining the latest research reports, PRDX6 overexpression accelerates the progression of cancers such as cervical cancer, gastric cancer, and skin cancer.^45,46 Nevertheless, multiple studies have conducted reverse validation and found that knocking down PRDX6 gene expression can inhibit tumor cell proliferation, migration, and invasion capabilities.^47–49 Jaeyul Kwon et al conducted experiments using multiple cell lines and found that the PLA₂ and peroxidase activities of PRDX6 are essential products supporting NOX1 activity.⁵⁰ The study further demonstrated that overexpression of PRDX6 directly promotes NOX1-mediated superoxide production, while knocking down PRDX6 inhibits NOX1 activity. Jin Ran Chen et al discovered substantial bone resorption in NOX4-deficient mice, and bone formation markers, including bone-specific alkaline phosphatase and procollagen type I N-terminal propeptide, decreased significantly at 3 weeks of age.⁵¹ This finding lends further credence to the notion that PRDX6 may serve as a pivotal therapeutic target for the treatment of SONFH, which aligns with the conclusions of our own research endeavours.

Our MR analysis further identified CCL17 as a key upstream regulator of PRDX6, with a strong causal association with SONFH risk. CCL17 plays a crucial role as a key member of the C-C chemokine family, primarily exerting its effects by binding to the chemokine receptor C-C chemokine receptor type 4 (CCR4) on target cell surfaces.⁵² Th22 cells represent a subset of CD4⁺ T cells characterized by high surface expression of three chemokine receptors, CCR4, CCR6, and CCR10, with CCL17 serving as the specific ligand for CCR4.^53,54 Yusuke Miyazaki et al demonstrated that CCL17, highly expressed in the synovial tissue of rheumatoid arthritis patients, guides Th22 cells from peripheral blood to cross the vascular wall.⁵⁵ These cells then migrate directionally to the site of synovial inflammation and accumulate in large numbers. This process intensifies the production of pro-inflammatory factors, thereby exacerbating bone destruction. Kevin M.-C. Lee et al established arthritis models using CCL17 knockout mice and revealed that pathological changes such as cartilage erosion and synovial hyperplasia were significantly suppressed.⁵⁶ The study further validated the upstream-downstream relationship between CCL17 and tumor necrosis factor(TNF)/cyclooxygenase(COX) using anti-TNF monoclonal antibodies and COX inhibitors. This demonstrated that the functional effects of CCL17 are closely intertwined with those of both molecules. The team also confirmed that granulocyte-macrophage colony-stimulating factor acts as an upstream regulatory factor, inducing CCL17 production and thereby driving inflammatory progression and pain development.⁵⁷ Although CCL17 is widely reported to be involved in inflammatory regulation, an increasing number of recent studies have confirmed its close regulatory relationship with the ferroptosis pathway.⁵⁸ The CCL17/CCR4 axis functions as a pivotal upstream regulator of the NF-κB signaling pathway.⁵⁵ Concurrently, a substantial body of research has validated NF-κB as a pivotal regulatory node of ferroptosis. It has been demonstrated that NF-κB can directly regulate the expression of GPX4, SLC7A11, and ACSL4, which are key ferroptosis-related genes.⁵⁹ Recent experimental research has demonstrated that the activation of CCL17/CCR4 signaling significantly promotes the production of intracellular ROS, a hallmark of ferroptosis.⁶⁰ Study also found that elevated CCL17 levels were accompanied by increased oxidative stress markers such as MDA.⁶¹ Our two-step MR analysis further identified PRDX6 as a key mediator between CCL17 and SONFH, and PRDX6 is a well-documented key regulator of ferroptosis that exerts an antioxidant effect by clearing lipid peroxides.

The remaining two targets, C7ORF50 and APOE, exhibited divergent characteristics. While C7ORF50 demonstrated robust binding affinity for DHA in molecular docking and MM/PBSA analysis, its regulatory impact on the ferroptosis-SONFH pathway was comparatively weaker than that of CCL17 and PRDX6 in MR mediation analysis. Conversely, APOE exhibited an exceptionally low binding affinity for DHA. These observations suggest that these two targets might not function as the primary therapeutic targets of DHA in addressing SONFH.

Conclusion

This study suggested that the CCL17-PRDX6 axis may be the core pathway through which DHA regulates ferroptosis to treat SONFH, and our findings provide a valuable reference for the mechanistic research and clinical management of SONFH.

Several limitations of this work should be acknowledged. The core targets have not been further validated through molecular functional experiments, and the sample size of the animal experiment is relatively small. In addition, large-scale, high-quality proteomic pQTLs data and SONFH-specific GWAS data from the Chinese population are not yet available, which limits the implementation of population-specific MR analysis at this stage. Correspondingly, follow-up studies will combine in vivo and in vitro experiments to confirm the role of the ferroptosis pathway in DHA-mediated regulation of bone metabolism, and focus on collecting multi-omics data from Chinese patients with SONFH to verify our findings and explore the population-specific pathogenesis of SONFH.