Abstract Background Ossification of ligamentum flavum (OLF) is a hidden, indolent disease condition with variable unexplained etiology and pathology. Growing evidences show a correlation between senile osteoporosis (SOP) and OLF, but the fundamental relationship between SOP and OLF remains unclear. Therefore, the purpose of this work is to investigate unique SOP-related genes and their potential functions in OLF. Methods Gene Expression Omnibus (GEO) database was utilized to gather the mRNA expression data ([49]GSE106253) and then analyzed by R software. A variety of methods, including ssGSEA, machine learning (LASSO and SVM-RFE), GO and KEGG enrichment, PPI network, transcription factor enrichment analysis (TFEA), GSEA and xCells were employed to verified the critical genes and signaling pathways. Furthermore, ligamentum flavum cells were cultured and used in vitro to identify the expression of the core genes. Results The preliminary identification of 236 SODEGs revealed their involvement in BP pathways associated with ossification, inflammation, and immune response, including the TNF signaling pathway, PI3K/AKT signaling pathway and osteoclast differentiation. Four down-regulated genes (SERPINE1, SOCS3, AKT1, CCL2) and one up-regulated gene (IFNB1) were among the five hub SODEGs that were validated. Additionally, they were performed by ssGSEA and xCell to show the relationship of immune cells infiltrating in OLF. The most fundamental gene, IFNB1, which was only found in the classical ossification- and inflammation-related pathways, suggested that it may affect OLF via regulating the inflammatory response. In vitro experiment, we found that IFNB1 expression was dramatically higher in cells cocultured with osteogenic induction than in controls. Conclusion As far as we are concerned, this is the first observation using transcriptome data mining to reveal distinct SOP-related gene profiles between OLF and normal controls. Five hub SODEGs were ultimately found using bioinformatics algorithms and experimental verification. These genes may mediate intricate inflammatory/immune responses or signaling pathways in the pathogenesis of OLF, according to the thorough functional annotations. Since IFNB1 was discovered to be a key gene and was connected to numerous immune infiltrates in OLF, it is possible that IFNB1 expression has a substantial impact on the pathogenesis of OLF. Our research will give rise to new possibilities for potential therapeutics that target SOP reverent genes and immune-associated pathways in OLF. Keywords: Ossification of ligamentum flavum, Immune infiltration, Senile osteoporosis, IFBN1, Inflammatory/immune response 1. Introduction Ossification of the ligamentum flavum (OLF) is a rare ectopic ossification disease with unelucidated pathogenesis. Serious spinal stenosis and thoracic myelopathy might develop when the ligamentum flavum, a crucial linking structure in the posterior column of the spine, is replaced by ectopic bone and compresses the posterior column [[50]1]. Due to the insidious and indolent symptoms in the beginning stage, it usually causes the diagnosis to be hysteresis. Once patients experience obvious symptoms such as muscle weakness or numbness, various degrees of difficulty walking, and even paralysis due to the large ectopic ossification, it usually means that patients need surgical interventions [[51]2]. Till now, congenital (genetic factors) and non-genetic factors such as obesity, age, inflammatory responses, biomechanics, etc. have been thought to be implicated in the development of OLF [[52][3], [53][4], [54][5]]; however, its precise etiology remains vague. Epidemiological research had shown that the OLF were closely related to aging and that older individuals were more likely to experience OLF [[55]6,[56]7]. Additionally, the relationship between a variety of bone metabolic diseases and OLF has been studied by a number of investigators. For instance, concrete evidence has shown that individuals with diffuse idiopathic skeletal hyperostosis (DISH) and ossification of the posterior longitudinal ligament (OPLL) experienced lower rates of osteopenia and osteoporosis than normal individuals, and individuals with OLF experienced the same tendency in the rates of osteopenia and osteoporosis, but with no statistically significant difference due to the size of the study cohort [[57][8], [58][9], [59][10]]. These studies suggest that the pathogenesis of OLF may be a potential treatment to reverse bone loss in SOP, but the underlying mechanisms continue to elude researchers. In order to uncover the molecular mechanisms behind SOP-related OLF, advanced bioinformatics techniques for high-throughput mRNA sequencing data were applied to gather SOP-related deferentially expressed genes (SODEGs) in ossification samples and controls and assess their further functions. With regard to this work, we have consistently used a variety of analytical techniques, databases, and experimental methods to validate the involvement of immune infiltration and SOP-related genes in the pathogenesis of OLF. Advanced bioinformatics algorithms were implemented to screen and validate the expression of core genes and enriched pathways, along with the PPI network, edge-betweenness and random walk methods, function enrichment, machine learning, co-expression analysis, transcription factor enrichment analysis (TFEA), immune cell infiltration and correlation analysis. Lastly, the relationship between IFNB1 expression and the immune cell was investigated via qRT-PCR following a series of in vitro experiments. These observations would provide us with a theoretical foundation for understanding the molecular mechanisms that underpin SOP and OLF. 2. Materials and methods 2.1. Collection of OLF-related mRNA expression data and human SOP-related genes For further analysis, dataset [60]GSE106253 was gathered from the GEO database ([61]http://www.ncbi.nlm.nih.gov/gds/) which was derived from [62]GPL21827 platform of Agilent-079,487 Arraystar Human LncRNA microarray V4, it included a total of eight mRNA and LncRNA expression profile samples. Specifically, four of these samples were ligamentum flavum tissue samples from four different OLF patients, while the other four were from four different healthy populations. The Online Mendelian Inheritance in Man (OMIM) database ([63]https://omim.org/), the OsteoporosAtlas database [[64]11] and the GeneCards database ([65]https://www.genecards.org/) were used to gather all human SOP-related genes (SORGs). 2.2. Identification of SOP-related DEGs in OLF The raw data was utilized to perform data normalization by R software. Box diagrams and principal component analysis (PCA) were used to depict the normalized raw data. Differential expansion genes (DEGs) between OLF and normal were filtered out and displayed on volcano maps using the cutoff criteria |fold change| greater than one and an adjusted P less than 0.05. Then, the DEGs of OLF were crossed with SORGs to obtain OLF differential genes associated with SOP (SODEGs). A clustering heatmap was created to illustrate how SODEGs were expressed. 2.3. Functional and pathway enrichment analysis of SODEGs We included a range of methods to perform mutual verification by employing numerous enrichment analysis tools, including Metascape, David, Profiler, ClueGO, WebGestalt, and R software, in order to thoroughly evaluate the GO (gene ontology) and biological pathways of SODEGs. GO analysis was performed to illustrate their functions in BP (biology process), CC (cell component), and MF (molecular function). Enrichment pathway analysis was conducted in three well-known databases: KEGG (Kyoto Encyclopedia of Genes and Genomes), REAC (Reactome), and WP (WikiPathways). 2.4. Construction of protein-protein interactions (PPI) network and select candidate sub-network To forecast protein interactions and build a PPI network, SODEGs were uploaded to the STRING database ([66]http://string-db.org). After taking into account the previous literature and our deliberations [[67]3,[68]12,[69]13], an interaction score greater than 0.4 was finally chose to balance between including a sufficient number of potential interactions for downstream analysis and minimizing the inclusion of false positive interactions. The edge-betweenness and random walk methods were subsequently employed to emphasize sub-networks or neighborhoods of the gene connectivity data from the STRING database, and the divided clusters were then depicted by further function enrichment in R software. According to the degree of the topological attribute of genes, sub-networks were extracted and visualized from the PPI network using Cytoscape. Then these candidate genes in the extracted sub-network were uploaded to the Genemania website ([70]http://genemania.org/search/) to further depict the interrelationships and functions. Then the function analysis of these candidate genes was performed and visualized by ClueGO. Additionally, Metascape ([71]https://metascape.org/gp/index.html) was also utilized to further clarify the potential biological function of these candidate genes. A p value less than 0.05 was considered statistically significant. 2.5. Transcription factors enrichment analysis (TFEA) and the acquisition of hub genes In pursuit of identifying the genes associated with ossification, the OLF-DEGs were meticulously uploaded to the ChEA3 website ([72]https://maayanlab.cloud/chea3/) for performing TFEA. The resulting ossification-related genes were then thoroughly scrutinized by integrating them with the topologically filtered candidate genes and genes in center cluster (Subnetwork) identified through the edge-betweenness and random walk methods to identify the hub genes. In a quest to unravel the core gene regulating OLF via the immune pathway, sophisticated machine learning techniques, including the least absolute shrinkage and selection operator analysis (LASSO analysis) and the support vector machine-recursive feature elimination analysis (SVM-RFE analysis), were employed to unearth the feature genes, which were subsequently merged with the genes enriched in inflammation- and osteogenesis-related pathways. The interplay and relationship between the core gene set were scrutinized using the Spearman's correlation analysis, leaving no stone unturned in the pursuit of gaining a comprehensive understanding of the complex mechanisms underpinning ossification. 2.6. Identified immune infiltrating cells in OLF and their relationship with IFNB1 With a P value less than 0.05, the single sample gene set enrichment analysis (ssGSEA) and xCell algorithms were applied to visualize and depict immune infiltration and function based on mRNA expression data. Boxplots were conducted to depict the immune cells between OLF samples and normal ones. Spearman's correlation analysis was utilized to access the relationship between IFNB1 and immune cells, as well as the immune cell markers, to determine IFNB1's potential role in immunity. Memory B cells, CD8^+ T cells, T-helper 2 (Th 2) cells, T-helper 17 (Th 17) cells, neutrophils, mesenchymal stem cells (MSC), plasmacytoid dendritic cells, and classical dendritic cells (cDC) whose markers were chosen from R&D Systems' website ([73]www.rndsystems.com/cn/resources/cell-markers/immune-cells). 2.7. Clinical specimens from patients and primary ligamentum flavum (LF) cell culture All research processes were carried out under the agreement of the ethical committee of the First Affiliated Hospital of Guangzhou University of Chinese Medicine (Num: JY-2022-223). With the consent of the patients and their families, samples of the LF or OLF were taken from 6 patients (age range, 81–84 years; 2 males, 4 females) during surgical treatment between April 2022 and September 2022 in the first affiliated hospital of Guangzhou University of Chinese Medicine. The comprehensive clinical information of the patients is described in [74]Supplementary Table S1. The isolated surgical specimens were promptly brought to the Lingnan Medical Research Center of the Guangzhou University of Chinese Medicine for subsequent processing. The adherent tissue surrounding the ligament was carefully removed after washing the specimens three times with PBS solution. The cleaned LF specimens were uniformly cut into tissue blocks of approximately 3 mm^2 and digested by 0.25% type I collagenase (Biosharp, Guangzhou, China) for one and a half hours at 37 °C, then the blocks were washed three times again. The cleaned tissue blocks were evenly transplanted into culture dishes containing adequate medium (α-MEM, 10% fetal bovine serum, all from Gibco, Grand Island, NY, USA), and 1% penicillin/streptomycin (Cyagen, Guangzhou, China), which was changed every five days. The culture dishes were displayed in a standard cell culture chamber. At 21 days after transplanting, the fibroblast-like cells were migrating around the tissue; they were passaged using a 0.25% trypsin-containing EDTA solution (Cyagen, Guangzhou, China). Only 3 to 4 passage cells were used in our experiments. All the above operations are performed in a sterile environment. 2.8. Osteogenic differentiation and alkaline phosphatase assay LF cells were seeded into 24-well cell culture plates, after the cell proliferation density reached to 70% the medium was changed to osteogenic induction (OI) medium [[75]14] (Cyagen, Guangzhou, China). Cells were separated into three groups: Group A acquired α- MEM/FBS medium (all from Gibco, Grand Island, NY, USA); Group B received osteogenic induction medium for 7 days; and Group C received osteogenic induction medium for 14 days. Cells were maintained for 14 days by replacing the osteogenic induction medium and basal medium every three days. After 14 days of osteogenic induction, LF cells were treated with 4% paraformaldehyde and assayed for ALP activity using the BCIP/NBT Alkaline Phosphatase kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). 2.9. RNA extraction and real-time RT-qPCR RT-qPCR was used to assess the expression levels of the core genes. The LF or OLF specimens were frozen in liquid nitrogen and crushed into 3 mm^2 tissue blocks before extracting Total RNA using the TRIzol Reagent (Thermo, USA). To reverse transcribe RNA to cDNA, the PrimeScriptTM RT Master Mix Kit (Takara, Shiga, Japan) was utilized. The SYBR Premix Ex TaqTM II Kit was used for quantitative analysis on a Bio-Rad qPCR equipment (Takara, Japan). 2^−ΔΔCt method was used to calculate the relative expressions. The primer sequences were shown in [76]Supplementary Table S2. 2.9.1. Statistical analysis and visualization R software (V4.0.0) was used to conduct statistical analyses and visualization. Nonparametric test and parametric test were performed where appropriate. The data was processed in mean ± SD form. Statistically significant values were identified when P values were less than 0.05. 3. Results 3.1. Flow chart of the entire study and acquisition of SODEGs [77]Fig. 1 depicted the comprehensive flow chart for this research. Specifically, we first conducted differential expression analysis on [78]GSE106253 to obtain DEGs of OLF, and these DEGs were then intersected with SORG to acquire significant SODEGs. Functional enrichment analysis in different tools was conducted to seek the potential biological function of OLF, and then the PPI network and sub-network were acquired. Further function enrichment analysis in GeneMANIA, clueGO, and Metascape were performed to analyze the 25 genes in the sub-network, and then the final five core genes were obtained, among them IFNB1 was thought to may have the capacity to influence OLF through regulating the inflammatory response. The immune infiltration algorithm revealed a correlation of the core genes with immune cells, and qRT-PCR was used to identify the expression of the five core genes. Lastly, we cultured primary human ligament cells and validated the expression of relevant genes after osteogenesis. Fig. 1. [79]Fig. 1 [80]Open in a new tab Workflow of this study. Our work was divided into two parts: bioinformatics and in vitro experiments. In the bioinformatics part, we conducted ample analysis and finally found that IFNB1 was the most significant gene; In the in vitro experiments part, we cultured primary human ligament cells and validated the expression of relevant genes after osteogenesis. The analysis of the express data of OLF samples and normal samples was performed through R software. The PCA diagram showed a significant difference in OLF samples and controls after data normalization, and the median lines of the data expression box plot were at the same level, which allowed further analysis ([81]Fig. 2A and B). R software was then used to identify DEGs in OLF samples and controls, with a cutoff threshold of | FDR | > 1 and an adjusted P value less than 0.05; 577 genes showed higher expression and 391 genes showed lower expression, which were identified ([82]Fig. 2C). After taking intersections with 4964 SOP-Related genes, a total of 236 significant SODEGs were obtained, which were displayed with a heatmap that clearly distinguished OLF samples from controls ([83]Fig. 2D and E). Fig. 2. [84]Fig. 2 [85]Open in a new tab Row Data Normalization and SODEGs Acquisition. (A–B) PCA diagram and data expression box between OLF samples and controls. (C) Volcano map showed 577 up expression genes and 391 down expression genes were identified. (D–F) 236 SOP-related OLF differential genes (SODEGs) were obtained, which overlapped by OLF DEGs and SOP-retaliated genes from OMIN database, GeneCard database and OsteoporosAtlas database. 3.2. GO enrichment analysis revealing potential biological functions To further investigate the potential functions of the SODEGs in OLF, GO functional enrichment analysis was performed by g:profiler, R software, and WebGestalt. GO analysis in MF was primarily enriched in collagen binding, protein binding, and nuclear receptor activity, in addition to 300 notable pathways in BP, 7 key pathways in KEGG, 13 major pathways in REAC, and 17 critical pathways in WP, which were all obtained by g: Profiler ([86]Fig. 3A and B). [87]Supplementary Figs. 1A and B showed specific results of the BP analysis, which suggested that the pathways annotated by WebGestalt—including system development, ossification, regulation of the inflammatory response, and response to oxidative stress—were largely involved in the developmental process and the inflammatory and immune response. Additionally, these SODEGs were significantly enriched in cellular developmental processes such as extracellular matrix organization and extracellular structure organization pathways using the clusterProfiler package in the R software ([88]Fig. 3C and D). Furthermore, multiple significant enrichment pathways, including ossification, leukocyte differential, and negative regulation of immune system processes, were identified by GSEA and shown in [89]Fig. 3E and F. The findings above suggest that these SODEGs may affect OLF formation by regulating inflammatory and immune responses as well as bone metabolism pathways. Fig. 3. [90]Fig. 3 [91]Open in a new tab Functional Analyses of SODEGs. (A) Functional analyses by g: Profiler, included 19 significant pathways for MF, 300 significant pathways for BP, 7 significant pathways for KEGG, 13 significant pathways for REAC, and 17 significant pathways for WP. (B) MF analyses by WebGestalt. (C–D) Scatter Plots about GO terms analysis in BP and specific information on the top ten categories. (E–F) Tree diagram about GSEA analysis and the detailed information about two significant enrichment pathways. Additionally, R software was unitized to assess the functions of up- and down-regulated genes independently. Down-regulated SODEGs were enriched in response to stimulus pathways, such as response to peptide hormone, response to lipopolysaccharide, and response to molecule of bacterial origin, while up-regulated SODEGs were strikingly associated with bone development pathways, including ossification, biomineralization, and extracellular matrix organization ([92]Supplementary Figs. 1C and D). Besides that, MF displays that down-regulated genes were primarily enriched in collagen binding, nuclear receptor activity, and ligand-activated transcription factor activity, and that up-regulated genes were primarily associated with energy metabolism, such as ATPase-coupled transmembrane transporter activity, ATPase-coupled cation transmembrane transporter activity, and ATPase-coupled ion transmembrane transporter activity ([93]Supplementary Figs. 1E and F). These genes could be related to multiple biological pathways orchestrating OLF pathogenesis. 3.3. Pathways enrichment analysis revealed highly overlapping pathways All SODEGs were uploaded to three different datasets to further investigate the potential roles of these genes, which included KEGG, REAC, and WP. Notably, the three datasets show striking similarities in the inflammatory/immune response and ossification-related pathways. Inflammatory/immune response pathways include TNF signaling pathway, JAK-STAT signaling pathway, and T cell receptor signaling pathway in the KEGG dataset from R software ([94]Fig. 4B), interleukin-4 and interleukin-13 signaling, and interleukin-10 signaling in the REAC dataset, and nuclear receptors in the WP dataset from G: Profiler ([95]Fig. 4A). Ossification-related pathways such as endochondral ossification and endochondral ossification with skeletal dysplasias in the WP dataset and ECM-receptor interaction, collagen formation, and collagen degradation in the REAC dataset from the G: Profiler ([96]Fig. 4A). Osteoclast differentiation, TGF-beta signaling pathway and PI3K-Akt signaling pathway in the KEGG dataset ([97]Fig. 4B and C) from R software were also significantly enriched. GSEA analysis showed a tendency for SODEG to be similarly enriched in immune response pathways ([98]Fig. 4F). Fig. 4. [99]Fig. 4 [100]Open in a new tab Pathways Analyses of SODEGs. (A) Three different algorithms for pathway analysis of genes on the website g: Profiler, included KEGG, REAC and WP. (B–C) KEGG pathways analyses by R software and petal chart about 6 significant enrichment pathways. (D) KEGG analyses about up-regulated genes showed 5 pathways were enriched, however no significant pathways were observed in up-regulated genes. (E) Detailed enrichment pathways of down-regulated genes and up-regulated genes through DIVID analysis. (F) Histogram and scatter slots showed the detailed information about significant enrichment pathways through GSEA analysis in RECA dataset. We then employed DAVID, WebGestalt, and R software to evaluate these up and down genes, respectively, to better identify the specific activities of SODEGs. In the WebGestalt RECA dateset, down-regulated genes were most in Interleukin-4 and Interleukin-13 signaling and extracellular matrix organization pathways, but no significant pathways were found in the WP dateset ([101]Supplemental Table S3). Endochondral ossification and extracellular matrix organization pathways were revealed in up-regulated genes from WebGestalt ([102]Supplementary Table S3). Pathways analyzed by DAVID and R software demonstrated that up-regulated genes were both enriched in osteoclast differentiation, ECM-receptor interaction, and TGF-beta signaling pathways, as well as in classical inflammation pathways, including the TNF signaling pathway. Besides, down-regulated genes were mainly associated with the longevity-regulating pathway and the TNF signaling pathway, which were observed in DAVID; however, no significant pathways were observed in R software ([103]Fig. 4D and E). When all database results were combined, the highlight pathways were the TGF-beta signaling pathway, the TNF signaling pathway, and two other remarkable crosstalk pathways, osteoclast differentiation and the PI3K-Akt signaling pathway; the detailed information was shown in [104]Supplementary Fig. 2. 3.4. PPI network construction and functional analysis of sub-network With an interaction score greater than 0.4, we used bioinformatic analysis to forecast the extracted subnetworks or neighborhoods in the PPI network based on functional annotations, which involved the STRING database and R software. The pathway enrichment analysis revealed that the PPI network was primarily divided into five clusters, with the cytokine-cytokine receptor interaction, TNF signaling pathway, osteoclast differentiation, JAK-STAT signaling pathway, and extracellular matrix (ECM) receptor interactions being the most prominent in each of these clusters, followed by the TGF-beta signaling pathway, neuroactive ligand-receptor interaction, focal adhesion, and PI3K-Akt signaling pathway ([105]Fig. 5A). Furthermore, the significantly connected subnetwork from the PPI network was extracted using Cytoscape's degree topology property analysis. It consists of 25 nodes, including AKT1, MMP9, HIF1A, CCL2, IL10, TNFSF11, VCAM1, NFKBIA, ACAN, BMP4, TFRC, ELN, COL2A1, SERPINE1, IFNB1, LOX, SOCS3, NCAM1, BMP7, NES, CHRD, CTSB, COL6A1, STAT6, and CD40LG, and 139 interactions ([106]Fig. 5B, [107]Supplementary Table S5). Interestingly, most of them also exist in the center and largest cluster. Information for these clusters was shown in [108]Supplementary Table S4. Fig. 5. [109]Fig. 5 [110]Open in a new tab PPI network and Further Functional Analysis. (A) PPI network performed by STRING database and R software used edge-betweenness and random walk methods, 5 mainly clusters and their KEGG enrichment analysis were obtained. Subnetworks (Neighborhoods) are colored and annotated with enriched functional categories. Gray lines, connections within a neighborhood; red lines, connections between neighborhoods; dark yellow (#e69f00), cluster1 and 5; light yellow (#f0e442), cluster 2; blue (#0072b2), cluster 3; pink (#cc79a7), cluster 4. (B) Candidate subnetwork selected by degree topology property of genes. (C) Genes in Subnetwork and their co-expression genes were analyzed using GeneMANIA. (D–E) Functional enrichment analysis through ClueGO obtained a number of significant immune response pathways. The results show that these genes in Subnetwork were mainly involved in Interleukin-4 and Interleukin-13 signaling (68.58%) and B cell activation involved in immune (16.35%) pathways. *p < 0.05, **p < 0.01. (F–G) GeneMANIA was utilized to conduct a co-expression analysis to further seek the function of these genes in the potential subnetwork. The results revealed a complex subnetwork whose functions were primarily linked to the control of inflammatory responses and immunological processes ([111]Fig. 5C). We next performed a function network to further analyze the subnetwork, which suggested that inflammatory and immune-related pathways were substantially implicated. It incorporates Interleukin-4 and Interleukin-13 signaling (68.58%), as well as B cell activation involved in immunity (16.35%) ([112]Fig. 5D and E). These genes in the subnetwork were primarily associated with inflammatory and immune-related responses, including immune system cytokine signaling, TNF signaling pathways, IL-18 signaling pathways, and inflammation disease, based on the Metascape database ([113]Fig. 5F and G). Therefore, these genes were the main influences on OLF through inflammation and immune-related pathways. Detailed enrichment pathways of down-regulated genes and up-regulated genes through DIVID analysis. (F) Metascape functional annotation for genes in Subnetwork also pointed to Inflammation/Immunity response. 3.5. Identified ossification-related genes in sub-network and validation core genes by qRT-PCR Five core genes (IFNB1, AKT1, CCL2, SOCS3, and SERPINE1) were obtained by the intersection of genes analyzed by Cytoscape software through degree topology property, cluster 3 in the PPI network, and TFEA through ChEA3 in order to further screen the ossification-related genes ([114]Fig. 6A and B). The violin chart suggests that the expression of all five core genes was statistically significant ([115]Fig. 6C), and the detailed information of genes in the TFEA module were displayed in [116]Supplementary Table S6. Additionally, IFNB1 was identified as the most essential one because it was found in the aforementioned inflammation- and ossification-related pathways and machine learning geneset ([117]Fig. 6B), which implied that it had the capacity to influence OLF through regulating the inflammatory response. The results of machine learning were showed in [118]Supplementary Fig. S3. We then performed a correlation analysis between IFNB1 expression and another four core genes to investigate their interrelationships ([119]Fig. 6D). With a Spearman's correlation value of 0.9, IFNB1 and AKT1 had the largest negative correlation among them, whereas SOCS3 and IFNB1 showed no statistically significant correlation. To further identify the correctness of above predictions, qRT-PCR was used to assess the five core genes' expression described above in clinical specimens ([120]Fig. 6E). Finally, row data and qRT-PCR showed comparable changes in the expression levels of CCL2, AKT1, and IFNB1, although SOCS3 and SERPINE1 showed a different trend ([121]Fig. 6F). Fig. 6. [122]Fig. 6 [123]Open in a new tab TFEA and qRT-PCR Validation of Five Core Genes. (A) TFEA suggested 65 genes significantly associated with ossification. (B) Venn diagram showed 5 core genes from the intersection between ossification-related genes by TFEA, genes selected by degree topology property and PPI central cluster, and the most core genes (IFNB1 and AKT1) were obtained by the overlapped classical ossification-related pathways and inflammation-related pathways. (C) Violin chart for all five core genes. (D) Correlation analysis of 5 core genes. (E) Imaging characteristics of one OLF patient. (F) Validation of the expression of five genes by qRT-PCR. The expression of Akt1, IFNB1 and CCL2 were paralleled in row data, but there was an opposite trend in SOCS3 and SERPINE1 expression. *p < 0.05, **p < 0.01, ***p < 0.001. 3.6. Correlation analysis between IFNB1 expression and immune infiltrating cells as well as immune markers The ssGESA and xCELL algorithms depict the infiltration correlation over core genes and immune cells ([124]Fig. 7A and B). Combining these two algorithms, a total of 15 distinct immune cells showed various infiltrations and were intimately connected with OLF ([125]Fig. 7C). Then we further performed a correlation analysis between IFNB1 expression and these 15 different immune cells and their immune markers to investigate the interrelationships ([126]Fig. 7C). IFNB1 expression demonstrated a statistically significant association with 9 of 15 various immune cells, as indicated in the lollipop chart ([127]Fig. 7D). The degree of pro B-cell infiltration was shown to have a considerable connection (R = 0.95), showing that IFNB1 may influence pro B-cell in OLF. In addition, CMP, Basophils, Type 17 T helper cell, Memory B cell, Neurons, cDC, MSC also showed the positive connection and numerous immune cells were also discovered to be negatively correlated, including the amount of plasmacytoid dendritic cell infiltration (R = −0.95), smooth muscle infiltration (R = −0.9). Fig. 7. [128]Fig. 7 [129]Open in a new tab Immune Infiltration and Correlation Analysis. (A–B) xCcll and ssGSEA algorithms proved the correlation between immue cells and core genes. (C) Box plots showed statistically significant expression of 15 different immune cells between OLF and normal noes. (C) Box plots for the expression of all five core genes analyzed by Wilcox text. (D) Correlation analysis of IFNB1 between 15 different immune cells. (E)The expression of IFNB1was positively correlated with the Th2 cell, Th17 cell and Memory B cell marker (CD19) and MSC markers (CD19 and KIT). (H) The expression of IFNB1 was negatively correlated with the Th17 cell markers (IL1R1, STAT3 and (TGFBR2), Th2 cell marker (STAT6), Plasmacytoid dendritic cell (ILR3A and IRF8), MSC marker (ALCAM), CD8^+ T cell (IL2RA),. *p < 0.05. In addition, we then explored the association between IFNB1 and immunological markers of diverse immune cells ([130]Supplementary Table S7). Here we found that the expression of CD19, which was the immunological marker of Th2 cells, Th17 cells, the memory B cell marker, and MSC, was highly associated with IFNB1 ([131]Fig. 7E, H). Furthermore, the expression of IFNB1 was shown to be substantially linked with KIT in MSC (R = 0.76) ([132]Fig. 7E). Plasmacytoid dendritic cell levels and INFB1 expression had a poor connection (IL3RA, R = −0.81 and IRF8, R = −0.85). Interestingly, we discovered that the majority of T cell markers in OLF, including IL1R1 (R = −0.73), TGFBR2 (R = −0.76), and STAT3 (R = −0.73) for Th17 cells, STAT6 (R = −0.76) for Th2 cells, and IL2RA (R = −0.73) for CD8^+ T cells, were strongly connected with the expression level of IFNB1 ([133]Fig. 7F). These results indicated that IFNB1 could take part in the T-cell immune response in the pathology of OLF. 3.7. IFNB1 shows elevated expression during osteogenic induction We generated normal LF cells from tissue samples taken from SOP patients in order to further illuminate the properties of IFNB1. The 3rd–4th generation LF cells were ultimately employed for osteogenic induction, ALP staining, and an IFNB1 expression assay after the primary culture had been in place for the full 21 days. The primary culture process is depicted in [134]Fig. 8A. Then, we used ALP staining ([135]Fig. 8B) and qRT-PCR to assess IFNB1 expression during osteogenic induction. We found that cells cocultured with osteogenic induction displayed a dramatically increase in IFNB1 expression compared to controls, and the expression of osteogenesis-related genes (BMP4, OCN, and ALP) also increased with increasing induction time ([136]Fig. 8C). Overall, we confirmed IFNB1 can be elevated through osteogenic induction in LF cells. Fig. 8. [137]Fig. 8 [138]Open in a new tab IFNB1 Elevated Dramatically During Osteogenic Induction. (A) Primary LF cell culture process and cell morphology. (B) ALP staining showed positive after osteogenic induction. (C) q RT-PCR revealed the expression of IFNB1 in osteogenic induced LF cells displayed a 6-7-fold increased compared to controls and the expression of osteogenesis-related genes (ALP, OCN, BMP4) showed the same tendency. 4. Discussion OLF was one of the heterotopic ossification diseases where the ligamentum flavum (LF) becomes progressively ossified; it was also one of the commonest causes of low thoracic spinal stenosis and spinal cord compression, often leading to symptomatic spinal canal stenosis. The pathological process of OLF was similar to that of endochondral ossification and was characterized by the formation of heterotopic mature bone, which involved several stages including chondrocyte proliferation, fibroblast hyperplasia, osteogenesis, and maturation [[139]15,[140]16]. Since there was a selectively high prevalence of the condition in certain geographic areas, the underlying mechanisms of OLF were initially considered to be primarily linked to congenital (genetic) factors; however, as research progressed, non-genetic (environmental) factors such as age, inflammation, biomechanics, and so on began to emerge, causing OLF to gradually be considered a multifactorial disease [[141][3], [142][4], [143][5]]. These risk factors largely overlapped with SOP, creating a complex relationship between OLF and SOP. Seil et al. found a highly positive relationship between lumbar spine bone mineral density (BMD) and the number of spine levels exhibiting OLF (P = 0.03) and that OLF patients have a lower prevalence of osteoporosis (P = 0.41), both of which suggested that understanding the pathogenesis of OLF may provide potential treatments for osteoporosis. The underlying mechanisms connecting SOP and OLF, however, have rarely been studied in this domain. On this basis, we seek to identify crosstalk genes between SOP and OLF and discover their probable biological functions thoroughly with a variety of bioinformatics algorithms and in vitro experiments. We identify and validate 5 key SODEGs that have the ability to regulate ossification using qRT-PCR, with the up-regulated IFNB1 being considered to be a critical gene. Additionally, we revealed that IFNB1 was mostly connected with inflammatory signaling pathways (TNF signaling pathway) and numerous immunological infiltrate cells (B cells, T cells, plasmacytoid dendritic cells, smooth muscle, etc.). The current study revealed IFNB1 to be associated with immune infiltrates in then pathology of OLF for the first time, which could also provide a new insight to the therapeutic approaches of SOP-related OLF. Finally, this research revealed the functional signatures of SOP-relevant genes in OLF. The results of GO function analysis showed that all SODEGs were primary associated with bone metabolic pathways and inflammatory/immune response, including ossification, biomineralization, extracellular matrix organization (ECM), regulation of inflammatory response, and response to oxidative stress, which implicated the link between OLF and SOP in bone remodeling and inflammatory response. Research has proven that bone mass, mineralization, and matrix composition are vital factors for bone remodeling, which is greatly associated with the development of osteoporosis [[144]17]. By regulating cell adhesion and mineralization, the ECM was critical for maintaining bone homeostasis in bone tissue, and anomalies in the ECM may be at the root of a number of disorders characterized by an unbalanced state of bone homeostasis [[145]18,[146]19]. Numerous studies have shown that MMPs could indeed carry out a variety of biological functions in bone homeostasis by combining or degrading ECM. These functions include increasing bone resorption (MMP-9) [[147]20], resulting endochondral ossification and bone formation (MMP-14, MMP-13) [[148]21], and directly facilitating bone cell growth and proliferation and skeletal development in the early stages (MMP-2) [[149]22]. MMPs and ECM may represent a potential molecular pathway bridging SOP and OLF, as evidenced by the fact that these MMPs also mediate the chondrogenesis, osteogenesis, and heterotopic ossification processes in numerous heterotopic ossification models [[150][23], [151][24], [152][25]]. Moreover, the inflammatory/immune response pathway—TNF signaling pathway—was also mostly enriched, according to the pathways analysis. Overexpression of TNF-α may influence LF cell proliferation and osteoblast differentiation in OLF, demonstrating the intrinsic relationship between SOP and OLF [[153]4,[154]26,[155]27]. PPI subnetwork containing 25 hub genes was generated by using a novel algorithm paired with the degree topology property method and then subjected to further function analysis [[156]28]. The pathophysiology of OLF and SOP may be closely associated with a variety of cytokines or chemokines involved with inflammation and immunology, including NFKBIA, TNFRSF11, IL-10, and CCL2. CCL2 has been strongly linked to traumatic brain injury, where it may contribute to neurogenic heterotopic ossification, and it can also be secreted by macrophages, where it may play an important role in osteoblast development in heterotopic ossification (HO) [[157]29]. Additionally, a number of studies have hypothesized that CCL2 was an osteoporosis and osteopenia risk factor [[158]30,[159]31]. Bone morphogenetic proteins (BMP4, BMP7) and collagen genes (COL2A1, COL6A1) were also obtained. The latter were proposed to be OLF susceptibility genes in the Chinese Han population [[160]32]. The function of the former were extensively documented, suggesting it may induce ectopic bone production in vivo [[161]33]. Notably, AKT1 could modulate multiple processes, such as metabolism, cell proliferation, growth, and angiogenesis, it has also been demonstrated to have an impact on regulating cartilage calcification during the pathological process of OLF [[162]34]. VCAM-1, a pro-angiogenic gene, has been recognized as a key cell-surface marker correlated with MSCs, which have been implicated in bone creation and absorption in humans [[163]35]. Liu Z et al. found that the activation of VCAM-1 could repress the muscle ring-finger protein-1 (MURF1)-mediated ubiquitylation of PPAR2, resulting in the reduction of osteoblastogenesis from MSCs [[164]36], in addition, study had suggest that the endochondral ossification process of the ossification of spinal may initiate from the MSCs differentiating into chondrocytes [[165]37]. Together, these findings sparked certain innovative speculation on whether VCAM-1 could be a possible treatment target for OLF. Taking account of the results of GeneMANIA and the Metascape database, concrete evidence indicated that biological processes as well as inflammatory and immune responses were highly associated with hub genes, which included Interleukin-4 and Interleukin-13 signaling and B cell activation involved in immunity pathways. Altogether, we proved that these hub genes were essentially in the immune response during the development of SOP and OLF. Then, through qRT-PCR, we demonstrated that IFNB1, SOCS3, SERPINE1, AKT1, and CCL2 were potential DEGs associated with SOP and OLF pathogenesis. In this study, we found IFNB1 was the only immune-related gene that regulates ossification, which was simultaneously enriched in the TNF signaling pathway, osteoclast differentiation, the JAK-STAT signaling pathway, and the PI3K-AKT signaling pathway. As a member of type I IFN, IFNB1 exerts powerful anti-proliferative and pro-apoptotic functions in the immune response by combining multiple immune cells [[166]38]. AAlthough interferons were thought to be important immune system regulators, there is currently little evidence for OLF. Thus, the relationship between IFNB1 and the 15 different immune cell types in OLF was evaluated to investigate the potential immunoregulatory functions of IFNB1 in OLF, which may be the underlying mechanism for the association between OLF and SOP. We observed a substantial correlation between the expression of IFNB1 and various types of T and B cells. Evidence has revealed that B cells and T cells were essential for maintaining bone homeostasis; in particular, T cells played a crucial role in regulating OPG produced by B cells to maintain bone mass. The specific mechanism may be related to CD40 ligand (CD40L) and CD40 costimulation [[167]39], in addition, activated T cells and B cells could also secrete pro-osteoclastogenic factors including RANKL and TNF-α promoting bone loss in inflammatory states [[168]40]. This suggested that IFNB1 could be involved in the regulation of osteoblastic and osteoclastic factors through regulating B cells and T cells during the pathology of OLF. Research also showed that IFNB1 could promote M1 macrophage polarization and inflammation [[169]41]. Our cellular experiments prove a dramatic increase in IFNB1 expression compared to controls during osteogenic induction, which indicates that IFNB1 may be involved in osteogenic differentiation of LF cells; however, more in-depth mechanistic studies are still insufficient. These results suggested that IFNB1 expression may regulate and engage immune infiltrating cells, ultimately inducing inflammation related to OLF. To thoroughly understand the connection between IFNB1 and particular immunological responses in OLF, however, robust experimental investigation is required. This study still has certain restrictions. Unfortunately, since patient both suffer OLF and SOP are extremely rare for which it is challenging to gather adequate clinical samples, we are unable to confirm the underlying relationships and molecular mechanisms of these hub genes at this time. Since the situation may differ between patients with SOP and OLF, individual differences should also be taken into consideration. On the basis of this article, more extensive experimental research will be conducted in the future. Additionally, more research must be done to determine the precise mechanisms of the immunological responses that IFNB1 causes. All in all, we have used actual validation and advanced bioinformatics approaches to comprehend these and provide fresh perspectives that may be informative for subsequent mechanism research. 5. Conclusion As far as we are concerned, this is the first observation using transcriptome data mining to reveal distinct SOP-related gene profiles between OLF and normal controls. Five hub SODEGs were ultimately found using bioinformatics algorithms and experimental verification. These genes may mediate intricate inflammatory/immune responses or signaling pathways in the pathogenesis of OLF, according to the thorough functional annotations. Since IFNB1 was discovered to be a key gene and was connected to numerous immune infiltrates in OLF, it is possible that IFNB1 expression has a substantial impact on the pathogenesis of OLF. Our research will give rise to new possibilities for potential therapeutics that target SOP reverent genes andciated pathways in OLF. Author contribution statement You zhang: conceived and designed the experiments; Performed the experiments; analyzed and interpreted the data and wrote the paper. Hongwei Huang: conceived and designed the experiments. Honglin Chen: conceived and designed the experiments. Peng Zhang: Performed the experiments. Yu Liu: Performed the experiments. Yanchi Gan: Performed the experiments. Xianwei Yan: Performed the experiments. Bin Xie: Performed the experiments. Hao Liu: Performed the experiments. Bowen He: Performed the experiments. Jingjing Tang: Performed the experiments. Gengyang Shen: contributed reagents, materials, analysis tools or data. Xiaobing Jiang: contributed reagents, materials, analysis tools or data. Hui Ren: contributed reagents, materials, analysis tools or data. Data availability statement Data associated with this study has been deposited at NCBI; [170]GSE106253. Additional information Supplementary content related to this article has been published online at [URL]. Source of funding This study was supported by Innovative Team Project of the Department of Education of Guangdong Province (2021KCXTD017), High-Level University Collaborative Innovation Team of GZUCM (2021xk57)， Iternational Program for Postergraduates, Guangzhou University of Chinese Medcine. Declaration of competing interest All authors declare that they have no conflicts of interest. Acknowledgments