Identification and validation of ferroptosis key genes in bone mesenchymal stromal cells of primary osteoporosis based on bioinformatics analysis

Primary osteoporosis has long been underdiagnosed and undertreated. Currently, ferroptosis may be a promising research direction in the prevention and treatment of primary osteoporosis. However, the specific mechanism of ferroptosis in primary osteoporosis remains a mystery. Differentially expressed genes (DEGs) were identified in bone mesenchymal stromal cells (BMSCs) of primary osteoporosis and heathy patients from the GEO databases with the help of bioinformatics analysis. Then, we intersected these DEGs with the ferroptosis dataset and obtained 80 Ferr-DEGs. Several bioinformatics algorithms (PCA, RLE, Limma, BC, MCC, etc.) were adopted to integrate the results. Additionally, we explored the potential functional roles of the Ferr-DEGs via GO and KEGG. Protein–protein interactions (PPI) were used to predict potential interactive networks. Finally, 80 Ferr-DEGs and 5 key Ferr-DEGs were calculated. The 5 key Ferr-DEGs were further verified in the OVX mouse model. In conclusion, through a variety of bioinformatics methods, our research successfully identified 5 key Ferr-DEGs associated with primary osteoporosis and ferroptosis, namely, sirtuin 1(SIRT1), heat shock protein family A (Hsp70) member 5 (HSPA5), mechanistic target of rapamycin kinase (MTOR), hypoxia inducible factor 1 subunit alpha (HIF1A) and beclin 1 (BECN1), which were verified in an animal model.


Introduction
As a complex bone metabolism disorder mainly characterized by bone loss and bone microstructure changes, primary osteoporosis has long been underdiagnosed and undertreated (1)(2)(3)(4). Approximately 40% of postmenopausal white women are affected by osteoporosis. The aging population will exacerbate the medical and socioeconomic effects of osteoporosis. A patient with osteoporosis has a 40% lifetime fracture risk, and most of them occur in the spine, hip or wrist (5). Osteoporotic fractures affect 8.9 million people around the world per year. It is expected that osteoporotic hip fractures among elderly men will increase by 310% and 240% in women by 2050 (6,7). Osteoporosis and its complications have become a worldwide health crisis.
As a recently identified regulated cell death (RCD), ferroptosis differs from apoptosis (8), necrosis (9), autophagy (10), pyroptosis (11), etc. The main characteristics of ferroptosis are iron-dependent lipid peroxidation and reactive oxygen species (ROS) accumulation (12,13). It involves multiple signaling pathways and their regulators (14) and has been extensively studied in the treatment of cancers, such as renal carcinoma and leukemia (15,16) and hepatocellular carcinoma (17). Ferroptosis has been observed to be associated with various noncancer diseases, such as neurological diseases (18, 19), viral infection, ischemia and reperfusion injury, and atherosclerosis (20). Interestingly, a number of studies have suggested that ferroptosis may be a promising research direction in the prevention and treatment of osteoporosis (21)(22)(23).
Bone mesenchymal stromal cells (BMSCs) are pluripotent mesenchymal cells that can differentiate into different lineages, such as osteoblasts and adipocytes (24). BMSCs play an important role in maintaining bone homeostasis. As a result, osteoporosis is always accompanied by the low osteogenic potential of circulating BMSCs (25), the two form a vicious cycle, exacerbating the progression of osteoporosis. The delicate balance between osteogenic and lipogenic differentiation of BMSCs will be disrupted during the onset of osteoporosis (26). Under certain conditions, ferroptosis can also occur during the process of directional differentiation of BMSCs (27). At present, some researchers have also studied the mechanisms related to ferroptosis and several signaling pathways had been confirmed to play a role in the process of ferroptosis related osteoporosis (28). It is the existing basic research that provides inspiration and guidance for this study. Despite this, the specific mechanisms and related signaling pathways of ferroptosis remain mysterious in osteoporosis.
Bioinformatics analysis of ferroptosis genes in primary osteoporosis has not been reported to date. Thus, we identified DEGs in BMSCs of primary osteoporosis and healthy patients from the GEO databases with the help of bioinformatics analysis. By intersecting these DEGs with the ferroptosis dataset, ferroptosis DEGs (Ferr-DEGs) were obtained. After that, protein-protein interactions (PPI) were used to predict potential interactive networks. Then, the key Ferr-DEGs we obtained were further verified by an experimental OVX model. We have discovered the potential functional role of ferroptosis in BMSCs for the first time and have laid a possible foundation for understanding the pathological mechanism of primary osteoporosis.

Microarray data collection
Information on primary osteoporosis patients was downloaded from the Gene Expression Omnibus (GEO; http:// www.ncbi.nlm.nih.gov/geo), a public database that collects gene sequencing results (29). We downloaded the raw data of GSE35958 uploaded by Benisch from GEO (30). Five samples of primary osteoporosis (GSM878104, GSM878105, GSM878106, GSM878107, GSM878108; op group) and four controls (GSM878100, GSM878101, GSM878102, GSM878103; ctrl group) in GSE35958 were collected. A detailed table of patients' parameters can be seen in Supplementary Table S1.

Datasets analysis
Bioconductor was used to download the annotated R 4.0.3 package, and R was used to convert the microarray probes to symbols. When multiple probes corresponded to one gene ID, the average value was taken for analysis. Principal component analysis (PCA) was used to determine the significant difference dimensions with a P value < 0.05 (31). Relative log expression (RLE) was used for quality control (32). We used Limma to identify the DEGs. P-Value <0.05 and |log2-fold change (FC)| >1 were considered statistically significant (33). In addition, the Ferroptosis Database (FerrDb; http://www.zhounan.org/ferrdb) provided us with a dataset containing 259 genes (34). Then, we intersected it with the DEGs of GSE35958 to obtain ferroptosis DEGs (Ferr-DEGs).

Functional enrichment analysis and visualization
Metascape (http://metascape.org) was used to determine gene functions (35). Biological process (BP), cellular component (CC) and molecular function (MF) are the three major components of gene ontology (GO), which are used to describe gene functions and interactions (36).. The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database is an extensive database for mapping pathway annotation results (37,38). GO and KEGG analyses were performed to explore the functional roles of the Ferr-DEGs via Metascape, and Cytoscape (V3.7.2) was used for visualization (39,40). We then chose molecular complex detection (MCODE) to identify the interconnected central genes. Different colors are used to represent each MCODE network and the close interactions between molecules. During the interaction enrichment analysis, the minimum and maximum sizes of the selection network were 3 and 500, respectively.

Protein-protein interaction network analysis
Then, Ferr-DEGs were uploaded to STRING (https://stringdb.org/) to predict and construct PPI networks with a confidence level of >0.4 (41, 42). The node degree of proteins analyzed by STRING was used for ranking. Cytoscape was used to visualize the network. Nodes represent genes, and edges represent connections between genes. Genes were ranked by size according to the betweenness centrality (BC) score analyzed by the CytoNCA plugin. The cytoHubba plugin was used to rank genes by depth of the color corresponding to the weighted score based upon maximal clique centrality (MCC) algorithms (43,44).

Identification of key Ferr-DEGs
Protein and gene expression profiles at the tissue level (high, medium, low, NA) were obtained from the Human Protein Atlas database (HPA, https://www.proteinatlas.org). After calculation, we regarded the top ranked genes of each algorithm as hub genes and then took the trend of tissue expression and intersected most of those algorithms as the 5 most likely key Ferr-DEGs.

OVX model construction and cell culture
We selected twelve eight-week-old female C57BL/6 mice (provided by the Experimental Animal Center of Soochow University) for osteoporosis induction experiments (45). Mice were randomly divided into two groups: the ovariectomy (OVX) and sham groups. Eight weeks after surgery, all the specimens of the left femurs were removed for microCT and histological experiments. The right femurs were removed for BMSCs collection. Blood was collected by eyeballs removing, and the serum was separated by centrifugation at 2,000 r/m for 20 min and frozen at -80°C for future analysis. All procedures and experiments were approved by the Animal Ethics Committee of the First Affiliated Hospital of Soochow University.
BMSCs were obtained from the right femurs according to the protocol (46). Cells were seeded in a 10-cm 2 dish and incubated in a 37°C incubator with 5% CO 2 . BMSCs after three passages were used for subsequent experiments.

Micro-CT analysis
Left femur samples were scanned with a SkyScan 1176 micro-CT (Aartselaar, Belgium) (n = 6/group) after fixation in 10% buffered formalin for 48 h. Three-dimensional (3D) histomorphometric images were constructed by using NRecon software (SkyScan micro-CT, Aartselaar, Belgium). Bone mineral density (BMD), bone volume (BV), trabecular number (Tb. N), bone volume per tissue volume (BV/TV) and trabecular separation (Tb.sp) were used to evaluate bone mass. The region of interest (ROI) began with 100 pieces below the femur growth plate, and 200 pieces (6 µm each) were read per sample.

Lipid peroxidation assay
A Lipid Peroxidation Malondialdehyde (MDA) Assay Kit (Beyotime, China) was used to measure lipid peroxidation levels. Serum can be used directly for MDA detection. Cells were homogenized on ice in cell lysis buffer (Beyotime, China). After homogenization followed by centrifugation at 12,000×g for 10 minutes, the supernatant was collected for MDA assays. A 100 ml sample was added to each vial containing 200 ml of MDA detection working solution. The samples were mixed and incubated at 100°C for 15 mins. The samples were cooled to room temperature in a water bath and centrifuged at 1000×g for 10 minutes at room temperature. Then, 200 ml supernatant from each sample was pipetted into a 96-well plate, and the absorbance was measured at 532 nm with a microplate reader. MDA content was determined according to the standard curve.

Iron assay
Serum and cellular iron levels were measured by an iron Content Assay Kit (Solarbio, China). Serum can be used directly for iron assay. 100ml serum was mixed with 250 ml of working solution and incubated at 100°C for 5 mins. The samples were cooled to room temperature in a water bath, mixed with 62 ml chloroform, and centrifuged at 10000×g for 10 minutes at room temperature. Then, 210 ml supernatant from each sample was pipetted into a 96-well plate, and the absorbance was measured at 520 nm. BMSCs were homogenized on ice in cell lysis buffer (Beyotime, China). After homogenization followed by centrifugation at 12,000×g for 10 minutes, 20 ml supernatant was collected and mixed with 180 ml working buffer at 25°C for 10 min, and then the absorbance was measured at 510 nm. A standard curve was simultaneously generated according to the manufacturer's instructions.

Immunofluorescence staining
We used immunofluorescence (IF) staining to detect the expression of Ferr-DEGs between the two groups. Antigen repair of the sections was performed, and the sections were blocked with horse serum. Finally, the primary antibody was incubated at 4°C overnight, and fluorescent-labeled secondary antibodies (ab150079, abcam, UK) were incubated at room temperature for 1 hour. The localization and protein expression level of Becn1 (A7353, ABclonal, China) were observed by IF staining under an AxioCam HRC microscope (Carl Zeiss, Germany).

Statistical analysis
The quantitative data are presented as the mean ± SD. Student's t test was used to compare the differences between two groups. Each assay condition was performed in triplicate for all quantitative assays. A P value <0.05 was defined as statistically significant.

Data quality control and identification of Ferr-DEGs
A schematic diagram of this study and the main findings is shown in Figure 1. Samples from GEO were divided into two groups [primary osteoporosis (OP) group versus control (Ctrl) group]. PCA and RLE analysis were performed to determine available dimensions and screen correlated samples according to quality control standards and the normalization of raw data. PCA demonstrated significant differences between the op group and the ctrl group (Supplementary Figure S1A). The RLE plots showed that the normalization was acceptable (Supplementary Figure S1B).
We first analyzed the DEGs of GSE35958 by using the "Limma" package and presented it in a volcano plot  Figure S2). Then, we obtained a gene set including 259 genes from FerrDb and intersected them with the DEGs of GSE35958 to identify Ferr-DEGs. Eighty Ferr-DEGs were found and are shown in a heatmap and Venn diagram (Figures 2A, C). The Ferr-DEGs were further divided into driver, suppressor and marker according to FerrDb (Table 1).

Enrichment analysis of Ferr-DEGs
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database were used to analyze the biological classification of Ferr-DEGs by using Meyascape. The respective results of BP, CC, MF and KEGG are shown in detail (Supplementary Figure S3). The results showed that the significantly enriched genes were involved in the cellular response to stress, response to starvation, ferroptosis, autophagy, response to iron ion, etc. To further understand the relationships between different terms, we visualized it as a network plot (Figures 3A, B). Details of the top 20 terms are listed in the table ( Figure 3C).
The MCODE algorithm was applied to explore those highly connected modules, and three of them were found ( Figures 4A-D), The genes in Module 1 were associated with nongenomic actions of 1,25 dihydroxyvitamin D3, transcriptional activation by NFE2 like bZIP transcription factor 2(NFE2L2) in response to phytochemicals and regulation of protein transport. The genes in Module 2 were associated with autophagy and the response to starvation. The genes in Module 3 were related to protein processing and cellular responses to stress ( Figure 4E).

PPI network analysis and construction of Ferr-DEGs
Because one of the 80 genes was not related to other molecules and did not form a molecular network, we obtained a PPI network containing 79 nodes and 348 edges. Node degree

Selection and analysis of key Ferr-DEGs
We took the top gene in each algorithm as hub genes of different algorithms. After calculation, 10 Ferr-DEGs were obtained for further research. According to the literature, BMSCs tend to differentiate into adipogenic cells rather than osteogenic cells in osteoporosis patients. Therefore, protein and gene expression profiles at the tissue level (high, medium, low, NA) were obtained from HPA. We used the differences in their expression levels in adipose and bone marrow tissues as a preliminary reference for key Ferr-DEGs (Supplementary Figure S4). At the protein level, the expression of SIRT1, MAPK8(Mitogen-activated protein kinase 8), NFE2L2, HIF1A and BECN1 in bone marrow was higher than that in adipose tissues. In contrast, the expression of ATG5 was higher in adipose tissues, and the expression level of HSPA5 was the same, while XBP1(X-box binding protein 1) and JUN were not detected in either tissue. At the gene level, the expression of SIRT1, MAPK8, XBP1, HSPA5, HIF1A and BECN1 in bone marrow was higher than that in adipose tissues, and in contrast, the expression of ATG5, NFE2L2, MTOR and JUN was higher in adipose tissues. Finally, we considered the trends of protein level and gene expression in tissues and intersected most of those algorithms as the 5 most likely key Ferr-DEGs, including SIRT1, HSPA5, MTOR, HIF1A and BECN1 (Figures 6A, B).

Validation of the key Ferr-DEGs in vivo
To verify the key Ferr-DEGs in vivo, an OVX mouse model was used for validation. The OVX model was successfully established in C57BL/6 mice 8 weeks after surgery. The microCT results showed that the BMD, BV, BV/TV and Tb. N of the OVX group were decreased, while Tb.sp was increased (Figure 7). Meanwhile, the mRNAs of BMSCs isolated from the femurs were used to detect the expression of those 5 key Ferr-DEGs. The results showed that the relative mRNA expression levels of Sirt1, Hspa5, Mtor, Hif1a and Becn1 in the OVX group were significantly lower than those in the sham group (Figures 8A-E). To further validate the expression of these five differentially expressed genes in vivo, we chose Becn1 as the representative target, and IF staining was used for evaluation. The results showed that BECN1 was significantly more highly expressed in the sham group than in the OVX group (Figures 8F, G).
To further demonstrate the involvement of ferroptosis in the development of osteoporosis, lipid peroxidation levels and iron concentrations in serum and BMSCs from OVX mice were measured. The results showed that MDA levels and iron concentrations were significantly increased both in serum and BMSCs from OVX group (Figures 9A-D). Lipid peroxidation and iron overload are the characteristic manifestations of ferroptosis.
In summary, these in vivo results are consistent with the bioinformatics analysis, demonstrating that there is a high correlation between ferroptosis genes and the progression of osteoporosis in a mouse model. However, the exact mechanism needs to be further elucidated.

Discussion
Different types of osteoporosis, such as diabetic osteoporosis, glucocorticoid-induced osteoporosis, and postmenopausal osteoporosis, were confirmed to be associated with ferroptosis (22, 47-50). However, the specific mechanisms and signaling pathways of ferroptosis involved in osteoporosis remain a mystery. We first intersected microarrays on primary osteoporosis patients in GEO with FerrDb, especially the detection and analysis of BMSCs, which has not been reported in previous studies. We aimed to identify the key genes of ferroptosis and explore its mechanisms in primary osteoporosis. Therefore, we compared the gene expression of primary osteoporosis patients with controls. To be more statistically persuasive, several bioinformatics algorithms (PCA, RLE, Limma, BC, MCC, etc.) were adopted to integrate the results. Additionally, PPI network construction, GO and KEGG were performed to explore the potential functional roles of the Ferr-DEGs. Finally, 80 Ferr-DEGs and 5 key Ferr-DEGs were calculated for further study.
As a newly discovered form of RCD, ferroptosis is accompanied by iron overload, lipid peroxidation and oxidative stress and is regulated by manifold genes. There appears to be a synergistic or antagonistic effect between ferroptosis and other RCDs because they have cross regulators in some pathways (51-53). Interestingly, according to our findings in this research, autophagy was regarded as the most crucial biological process after ferroptosis. Autophagy is an intracellular degradation system that maintains the stability of the intracellular environment (54, 55). Including three subtypes: microautophagy, macroautophagy and chaperonemediated autophagy (CMA) (56)(57)(58). Several stress-related proteins were reported to be important regulators of autophagydependent ferroptosis in cancer therapy (59). New evidence suggests that abnormal levels of autophagy and mitochondrial autophagy disrupt the balance of bone metabolism. Osteoporotic BMSCs exhibit impaired osteogenic and increased adipogenic differentiation due to decreased autophagy (60, 61). Our results indicate that ferroptosis is involved in the process of osteoporosis and that autophagy is also involved. This conclusion was further supported by the MCODE results, which suggested that the Ferr-DEGs in different modules were connected with nongenomic actions of 1,25-dihydroxyvitamin D3, autophagy, and response to starvation and stress. We verified the potential connection between ferroptosis and autophagy in osteoporosis via various bioinformatics methods. The interaction between ferroptosis and autophagy deserves further study.
As a group of cells with the potential for multidirectional differentiation, BMSCs can differentiate into osteoblasts, chondrocytes and adipocytes. There is an opposite relationship between the differentiation of BMSCs into osteoblasts and adipocytes (62, 63). When osteoporosis occurs, BMSCs tend to differentiate into adipocytes rather than osteoblasts (64, 65). To test and verify our hypothesis, based on the HPA database, we compared the protein and gene expression levels of 10 hub genes between human bone marrow and adipose tissue. We considered the trends of protein and gene expression in tissues and intersected most of those algorithms as the 5 most likely key Ferr-DEGs, including SIRT1, HSPA5, MTOR, HIF1A and BECN1. Ultimately, these 5 key Ferr-DEGs were considered to be a preliminary reference for determining whether they are differentially expressed when osteoporosis occurs.
SIRT1 is an important regulator of bone homeostasis. The results of HPA showed that the expression level of SIRT1 in adipose tissue is NA, which is not entirely consistent with what is known thus far as SIRT1 is a potential regulator in adipose tissue inflammation and metabolism. The conflicting results may be due A B FIGURE 6 The expression analysis of hub genes in human tissues according to Human Protein Atlas data base. to delays in updating the site. However, the relative expression levels of SIRT1 in adipose tissue and bone marrow tissue were consistent, which did not affect the conclusion of this study. Sirt1 KO mice presented with osteoporosis characterized by decreased osteogenesis and increased adipogenesis in BMSCs (66). In contrast, BMSCs overexpressing Sirt1 repressed the increased expression of superoxide dismutase 2 (Sod2) and forkhead box O3A (Foxo3a), which promoted the differentiation of BMSCs into osteoblasts and delayed senescence (67). A study in head and neck cancer cells showed that silencing SIRT1 expression can inhibit epithelial-mesenchymal transition and decrease ferroptosis, while SIRT1 agonists can promote ferroptosis (68). SIRT1 was also reported to be associated with mitophagy (69). HSPA5 is a chaperone protein mainly expressed in the Endoplasmic reticulum (ER) that maintains the stability of glutathione peroxidase 4 (GPX4) by forming the HSPA5-GPX4 complex, thereby causing resistance to ferroptosis (70, 71). There seems to be no current research on the potential role of HSPA5 in osteoporosis, which would be an interesting direction. HIF1A degradation, circadian rhythm and lipid peroxidation play a role in the regulation of ferroptosis, and HIF1A is a key factor (72).   cell line was significantly reduced after knockdown of autophagy related 7 (Atg7) and Becn1 (77). Our previous study also found that tet methylcytosine dioxygenase 2(Tet2) promoted bone loss in OVX mice. It positively regulates BECN1-dependent autophagy by inhibiting BCL2 expression and promoting osteoclast differentiation (78). However, whether BECN1 is involved in ferroptosis during BMSCs differentiation is still unclear and will be further researched. Acts as a central regulator of the cellular response to growth stimulation, MTOR can inhibit autophagydependent ferroptosis. In cancer cells, interactions between MTOR and GPX4 can regulate autophagy-dependent ferroptosis by inhibiting CMA (79), and MTOR inhibitors may promote GPX4 degradation by activating the CMA pathway. MTOR also has a protective effect on ferroptosis in BMSCs, and autophagy mediated by the MTOR pathway can regulate the regeneration function of BMSCs, thereby influencing the occurrence and development of osteoporosis in postmenopausal women (61). Collectively, there are few studies about the specific functions and related mechanisms of these 5 key Ferr-DEGs in ferroptosis during osteoporosis progression. Our results provide a very meaningful direction for future research.
To further verify our conclusions, an osteoporotic model was successfully established in OVX mice. We found that the expression of Sirt1, Hspa5, Mtor, Hif1a and Becn1 in BMSCs from the OVX group was lower than that in BMSCs from the control group, which is consistent with our bioinformatic analysis results. The results of immunofluorescence staining were also consistent with our expectations. Detection of lipid peroxidation and iron levels at serum and cellular levels also confirmed that ferroptosis may be involved in the development of primary osteoporosis. These outcomes suggest that ferroptosis of BMSCs may be involved in the pathological process of primary osteoporosis. Finally, we believe that our results will offer new insight into the role of ferroptosis in primary osteoporosis and identify the 5 key Ferr-DEGs as potential biomarkers for primary osteoporosis diagnosis and treatment.

Conclusion
In summary, through a variety of bioinformatics methods, our research successfully identified 5 key Ferr-DEGs associated with primary osteoporosis and ferroptosis, namely, SIRT1, HSPA5, MTOR, HIF1A and BECN1. In addition, autophagy may also be involved in ferroptosis-related primary osteoporosis. In the future, we intend to collect enough bone tissues from patients with primary osteoporosis for further research. We aimed to clarify the regulatory roles of these genes to explore the differentiation mechanisms of BMSCs.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement
The animal study was reviewed and approved by Animal Ethics Committee of the First Affiliated Hospital of Soochow University.