Investigation of the Shared Biomarkers in Heterotopic Ossification Between Ossification of the Ligamentum Flavum and Ankylosing Spondylitis

Abstract

Study Design

Bioinformatics analysis of Gene Expression Omnibus (GEO).

Objective

Ossification of the ligamentum flavum (OLF) and ankylosing spondylitis (AS) represent intricate conditions marked by the gradual progression of endochondral ossification. This investigation endeavors to unveil common biomarkers associated with heterotopic ossification and explore the potential molecular regulatory mechanisms.

Methods

Microarray and RNA-sequencing datasets retrieved from the Gene Expression Omnibus (GEO) repository were harnessed to discern differentially expressed genes (DEGs) within the OLF and AS datasets. Subsequently, Weighted Gene Co-expression Network Analysis (WGCNA) was implemented to pinpoint co-expression modules linked to OLF and AS. Common genes were further subjected to an examination of functional pathway enrichment. Moreover, hub intersection genes were identified using the Least Absolute Shrinkage and Selection Operator (LASSO) regression, followed by an evaluation of diagnostic performance in external OLF and AS cohorts. Lastly, an analysis of immune cell infiltration was conducted to scrutinize the correlation of immune cell presence with shared biomarkers in OLF and AS.

Results

A total of 1353 and 91 Differentially Expressed Genes (DEGs) were identified in OLF and AS, respectively. Using the Weighted Gene Co-expression Network Analysis (WGCNA), 2 modules were found to be notably significant for OLF and AS. The integrative bioinformatic analysis revealed 3 hub genes (MAB21L2, MEGF10, ISLR) as shared risk biomarkers, with MAB21L2 being the central focus. Receiver Operating Characteristic (ROC) analysis exhibited a strong diagnostic potential for these hub genes. Gene Ontology (GO) analysis indicated their involvement in the positive regulation of myoblast proliferation. Notably, MAB21L2 was singled out as the optimal common biomarker for OLF and AS. Furthermore, an analysis of immune infiltration demonstrated a correlation between MAB21L2 expression and changes in immune cells. Activated CD8 T cells were identified as shared differential immune infiltrating cells significantly linked to MAB21L2 in both OLF and AS.

Conclusion

This study represents the first instance of identifying MAB21L2 as a prospective diagnostic marker for patients contending with OLF associated with AS. The research results indicate that the ECM-receptor interaction and the cell-cell adhesion may play a role in both disease processes. This newfound knowledge not only enhances our understanding of the pathogenesis behind spinal ligament ossification but also uncovers potential targets for therapeutic interventions.

Keywords

ossification of the ligamentum flavum ankylosing spondylitis endochondral ossification bioinformatics molecular mechanism inflammation immune

Introduction

Ankylosing spondylitis (AS) is a persistent inflammatory ailment that predominantly impacts the axial skeleton, particularly the sacroiliac joints and spine, giving rise to distinctive inflammatory dorsalgia.^1,2 It affects .55% of populations of European ancestry (herein termed Europeans)³ and .23% of Chinese,⁴ but is uncommon in Africans and Japanese. In addition to inflammatory back pain and arthritic degeneration, pathological new bone formation and spinal ankylosis stand as primary hallmarks of AS. When osteophytes span the entirety of a joint, the implicated joint becomes immobilized, culminating in axial joint rigidity, spinal ankyloses, and the potential for enduring disability. Although current treatments effectively quell inflammation in active AS, including nonsteroidal anti-inflammatory drugs (NSAIDs), Tumor Necrosis Factor Alpha (TNF-α) antibodies, and interleukin-17 inhibitors, these therapeutic interventions do little to mitigate the radiographic advancement of AS.^5-7

Ossification of the ligamentum flavum (OLF) is a prevalent spinal disorder characterized by the ectopic ossification of the ligamentum flavum, often serving as a leading etiology of thoracic spinal stenosis.⁸ It is most commonly found in middle-aged and elderly individuals, with a higher frequency of reports in Asia.^9,10 It predominantly occurs in the lower thoracic vertebrae.⁹ The pathological process of OLF involves endochondral ossification. The pathogenesis of OLF is influenced by a combination of systemic and local factors, including genetic predisposition, metabolic disturbances, abnormal secretion of sex hormones, and mechanical stress. However, the precise etiology and mechanisms underlying the transition from normal ligaments to ossified structures in OLF remain elusive. Investigating the molecular targets and cellular mechanisms of ectopic ossification in OLF could potentially provide more effective therapeutic options for this condition.

Ossification of the ligamentum flavum and AS both manifest in spinal ligament ossification. Although the pathogenesis of AS is not yet fully understood, inflammation and new bone formation have been recognized as the primary pathological processes in AS.¹¹ The supraspinous ligaments and ligamentum flava are the main targets of ossification in AS,^12,13 resulting in severe pathological changes. Prior research endeavors have underscored the significance of genetic predisposition and the role of inflammatory responses in OLF and AS, with an exploration of susceptibility loci in each of these maladies.¹⁴ Activation of the inflammatory response likely plays a contributory role in promoting the osteoblastic differentiation of ligament fibroblasts in OLF.¹⁵ MRI examinations of patients with AS indicate that the emergence of new bone formation in vertebral syndesmophytes is more inclined to transpire in regions with prior inflammation rather than in areas devoid of inflammation.¹⁶ Moreover, inflammation promoted the expression of osteoinductive proteins and indirectly facilitated the pathological formation of new bone at the enthesis.¹⁷ Based on the premise that inflammation triggers the reparative process preceding new bone formation, the inefficacy of anti-inflammatory treatments can be attributed to their administration occurring too late to preempt ectopic new bone development in AS patients who conform to diagnostic criteria involving enthesitis and unequivocal radiographic sacroiliitis. Furthermore, the preponderance of new vertebral syndesmophytes emerges at locations devoid of prior evidence of inflammation,^18,19 implying that new bone formation in AS is not solely induced by inflammation.

The pathological progression of OLF and AS entails the process of endochondral ossification.²⁰ However, the precise mechanism underlying pathological new bone formation in OLF and AS remains elusive. Osteophytes arise through the process of endochondral ossification, involving the differentiation of periosteal cells into chondroblasts and subsequently into osteoblasts.¹ Notably, recent observations have revealed chondrocyte differentiation within the ligaments of early-stage AS patients, with subsequent cartilage formation followed by calcification.²¹ Nonetheless, there is a paucity of research dedicated to unraveling the osteogenic mechanisms concerning transcriptional regulation and signaling pathways implicated in syndesmophyte formation in OLF and AS. Consequently, a broader spectrum of perspectives is imperative to delineate shared biomarkers and potential pathways associated with the advancement of ossification in spinal ligaments.

In light of the swift advancement in microarray and high-throughput sequencing technologies, bioinformatics methodologies have become increasingly prominent in investigating the interplay between diseases. The rapid evolution of bioinformatics approaches affords us a more profound comprehension of disease pathobiology at the genetic level.²² Utilizing high-throughput methods and bioinformatics approaches could illuminate the fundamental mechanisms driving ectopic ossification, thereby elucidating the genetic basis for their analogous manifestations in spinal ligament ossification. In this investigation, publicly available gene expression data from the GEO database were scrutinized. A systems biology methodology was subsequently employed to delineate shared gene pathways and diagnostic markers between OLF and AS. Subsequent to this, an assessment of the interaction between these diagnostic markers and infiltrating immune cells was conducted, aiming to attain a more profound comprehension of the pathophysiological mechanisms potentially uniting OLF and AS. The findings offered preliminary insights into shared biomarkers and pathways implicated in the pathogenesis of OLF and AS, thereby unveiling novel therapeutic targets. The research flowchart of this research was shown in Figure 1.

Figure 1.

Research design flow chart.

Materials and Methods

Data Collection and Preprocessing

The datasets GSE106253, GSE134290, GSE113212, GSE181716, GSE153310, and GSE73754 were acquired from the Gene Expression Omnibus (GEO) database at https://www.ncbi.nlm.nih.gov/gds/. For external validation, GSE113212, GSE181716, GSE153310, and GSE73754 were employed. The GSE106253 dataset, hosted on the GPL21827 platform, comprised 8 samples, including 4 OLF patients and 4 healthy controls. All individuals were Chinese. The GSE134290 dataset was based on the GPL570 platform. This dataset consisted of MSCs from 6 samples. MSCs from the enthesis (interspinous ligaments) involved in spinal ankylosis from 3 AS patients who had undergone spinal wedge osteotomy (AS MSCs) were isolated. MSCs from 3 non-AS individuals who underwent traumatic surgery at similar sites (control MSCs) were used as controls. All individuals were Taiwanese. The GSE113212 dataset encompassed 5 cases of OLF patients in the experimental group and 5 normal samples in the control group. All individuals were Japanese. The GSE153310 dataset, utilizing the GPL23227 platform, comprised MSCs from 3 AS patients in the experimental group and MSCs from 3 normal samples in the control group. All individuals were Chinese. The GSE73754 dataset included 52 cases of AS patients in the experimental group and 20 normal samples in the control group. All individuals were Canadian.

Differential Genetic Screening

Differentially expressed genes (DEGs) were identified within the GSE106253 and GSE134290 datasets. The limma R package was applied for DEG screening, utilizing a cutoff criterion of an adjusted P value <.05 and |log FC| ≥ 1.

Construction and Module Analysis of Weighted Gene Co-Expression Network

Weighted Gene Co-expression Network Analysis (WGCNA) was conducted using the R.4.0.3 software package.²³ The WGCNA R package was utilized to build co-expression networks, incorporating relevant clinical characteristics for the DEGs associated with OLF and AS. Before commencing the analysis, hierarchical cluster analysis was executed using the hclust function in the R language to identify and eliminate any outlier samples. Subsequently, in accordance with the criteria for constructing a scale-free network, the “pickSoftThreshold” function in the WGCNA software package was employed to determine the appropriate soft power (ranging from 1 to 20) for automated network generation. The outcomes were subject to topological overlap matrix analysis, resulting in the assignment of modules denoted by distinctive colors and their respective module features (ME). Furthermore, Pearson’s correlation test was applied to assess the relationship between ME and clinical features.

Enrichment Analysis

Gene Ontology (GO) enrichment analysis, a widely employed bioinformatics approach, was employed to extract comprehensive insights from extensive gene data. Furthermore, KEGG pathway enrichment analysis, recognized for illuminating biological mechanisms and functions, played a pivotal role. To enhance the visualization of the GO and KEGG pathways, the GO plot program package was enlisted. The GO enrichment analysis encompassed Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). An adjusted P-value of <.05 was considered as the threshold for statistical significance in the enrichment items. Similarly, in the KEGG analysis, an item was deemed statistically significant if it held a P-value of <.05. Additionally, the gene sets and gene expression matrices for OLF and AS underwent analysis via Gene Set Enrichment Analysis (GSEA) to shed light on potential regulatory pathways. GO enrichment analysis, KEGG pathway analysis and GSEA were also performed in R using the “clusterProfiler” package²⁴ to further investigate the biological roles of DEGs or hub gene.

Candidate Biomarker Expression Levels and Diagnostic Value

In the realm of data analysis, the R software ggplot2 package was harnessed for crafting boxplots that gauged the expression levels of the hub genes (significance at P < .05). To evaluate the efficacy of prospective biomarkers, the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) was employed. This was carried out on datasets GSE113212, GSE181716, GSE153310, and GSE73754, facilitated by the pROC R package.²⁵

Identification and Validation of Potential Shared Diagnostic Markers

The Lasso method, renowned for its effectiveness in regression analysis, capitalizes on the ℓ1 penalty to attain a sparse solution.²⁶ In the context of the intersection genes between DEGs and WGCNA, as mentioned earlier, the Lasso method played a pivotal role in selecting the core markers.

Analysis of Immune Infiltrating Cells

The single sample gene set enrichment analysis (ssGSEA) algorithm was employed to estimate the relative composition of various immune infiltrating cells and their functions based on mRNA expression data within immune gene sets. The immunologic features of all samples were evaluated using the ssGSEA technique using the “GSVA” R package.²⁷ Subsequently, the Wilcoxon test was applied to assess the proportion of immune infiltrating cells, and the outcomes were visually presented using boxplots. Furthermore, Spearman correlation analysis was utilized to evaluate the relationship with MAB21L2.

Results

Identification of DEGs

A total of 1353 differentially expressed genes (DEGs) were discerned between OLF patients and normal controls, comprising 904 upregulated and 449 downregulated DEGs, as illustrated in volcano plots (Figure 2A). The top 40 DEGs among these were visualized in heatmaps (Figure 2B). Furthermore, 91 DEGs were identified when comparing AS and normal controls, with 49 upregulated genes and 42 downregulated genes (Figure 2C). The top 40 DEGs for AS were depicted in heatmaps (Figure 2D). To gain deeper insights into the potential biological alterations associated with these DEGs, KEGG and GO enrichment analyses were conducted. GO results revealed that the DEGs in OLF were predominantly enriched in processes related to ossification and bone morphogenesis (Figure 2E). GO annotation for AS demonstrated that the DEGs were primarily enriched in functions associated with cell adhesion and cell adhesion mediator activity (Figure 2F). Intriguingly, GO and KEGG analysis uncovered a commonality in which both OLF and AS exhibited enrichment in the ECM-receptor interaction and cell adhesion (Figure 2G and H).

Figure 2.

(A) Volcano plots illustrating DEGs in GSE106253. (B) Heatmap plots visualizing DEGs in GSE134290. (C) Volcano plots displaying DEGs in GSE106253. (D) Heatmap plots depicting DEGs in GSE134290. (E) GO term analysis of DEGs in GSE106253. (F) GO term analysis of DEGs in GSE134290. (G) KEGG analysis of DEGs in GSE106253. (H) KEGG analysis of DEGs in GSE134290.

WGCNA Network Construction and Module Identification

WGCNA was utilized to pinpoint clusters of co-expressed genes displaying differential expression patterns in OLF and AS. The correlation between these combined modules and disease characterization was meticulously computed. For the OLF model, a soft threshold power of b = 26 was determined, while in the AS model, the chosen soft threshold power was 7 (Figure 3A and B). Upon amalgamating similar gene modules, eleven distinct modules were discerned in the OLF model, and an additional 13 modules were identified in the AS model (Figure 3C and D). It’s worth noting that the turquoise module exhibited the most robust correlation with OLF incidence (r = .87) (Figure 3C), whereas the brown module exhibited the strongest correlation with AS occurrence (r = −.88) (Figure 3D).

Figure 3.

(A) Cluster dendrogram representing co-expression in OLF. (B) Cluster dendrogram depicting co-expression in AS. (C) Correlation analysis between modules and clinical traits in OLF. (D) Correlation analysis between modules and clinical traits in AS.

Identification of Shared Genes and Shared Pathways

Three central genes (MAB21L2, MEGF10, ISLR) emerged as a nexus between the genes identified through WGCNA and DEGs, representing potential crosstalk genes linking both maladies (Figure 4A). Furthermore, we conducted GO Biological Processes (GO BPs) enrichment analyses on these 3 core genes to unravel shared regulatory pathways. As depicted in Figure 4B, these central genes are primarily associated with several GO BPs, including positive regulation of myoblast proliferation, myoblast migration, skeletal muscle satellite cell proliferation, skeletal muscle cell proliferation, regulation of myoblast proliferation, and myoblast proliferation. These biological processes are intrinsically linked to muscle development and regeneration. Research has indicated that myogenic stem cells, such as myoblasts, may serve as potential progenitor cells for osteoblasts during ectopic bone differentiation.²⁸ Therefore, we hypothesize that this group of hub genes may be involved in heterotopic ossification in both OLF and AS.

Figure 4.

(A) Venn diagram illustrating the convergence and overlap of 3 core genes selected through the intersection of genes identified by WGCNA and DEGs. (B) GO term analysis, with a focus on biological processes, of the 3 core genes.

Identification of hub Genes Associated by Lasso Regression

Subsequently, a LASSO regression algorithm was deployed to pinpoint the potential shared diagnostic genes. In GSE106253, the LASSO analysis revealed 3 of the 3 core cross-genes at the most appropriate λ = .000655 (Figure 5A and B). Similarly, in GSE134290, the LASSO analysis identified 3 of the 3 core cross-genes at the most suitable λ = .000674 (Figure 5C and D). Ultimately, this analysis unveiled 3 intersecting genes (MAB21L2, MEGF10, ISLR) as the optimal shared diagnostic biomarkers for OLF and AS (Figure 5E).

Figure 5.

(A) 10-fold cross-validation utilized to determine the optimal tuning parameter log(lambda) within the GSE106253 database. (B) LASSO coefficient profiles representing the diagnostic genes within the GSE106253 database. (C) 10-fold cross-validation used to identify the most suitable tuning parameter log(lambda) in the GSE134290 database. (D) LASSO coefficient profiles displaying the diagnostic genes in the GSE134290 database. (E) Venn diagram demonstrating the optimal diagnostic biomarkers.

The Validation of the Diagnostic Marker Genes

We assessed the sensitivity and specificity of the candidate biomarkers. Further external validation was conducted to evaluate the diagnostic efficacy of MAB21L2, MEGF10, and ISLR in multiple datasets, including GSE113212, GSE181716, GSE73754, and GSE153310. In the GSE113212 dataset (Figure 6A), these 3 biomarkers exhibited varying diagnostic values: ISLR (AUC = .594), MAB21L2 (AUC = .875), and MEGF10 (AUC = 1.000). In the independent OLF validation set GSE181716, the AUC values for ISLR, MAB21L2, and MEGF10 were .720, .640, and .800, respectively (Figure 6B). In the AS dataset GSE73754 (Figure 6C), MAB21L2 demonstrated a good diagnostic value with an AUC of .633. ISLR and MEGF10, on the other hand, exhibited lower diagnostic values with AUCs of .568 and .579, respectively. Another ROC analysis was performed for GSE153310, where the AUC values for ISLR, MAB21L2, and MEGF10 were .556, 1.000, and .556, respectively (Figure 6D). Ultimately, MAB21L2 emerged as the most promising diagnostic biomarker with stable and robust diagnostic value for both OLF and AS. Interestingly, differential expression analysis revealed that MAB21L2 was significantly overexpressed in AS (Figure 6F) but significantly decreased in OLF (Figure 6E). To further understand the molecular processes affected by MAB21L2, GSEA analysis was conducted. In OLF (Figure 6G and H), the diagnostic gene MAB21L2 was found to be significantly enriched in terms associated with OXIDATIVE PHOSPHORYLATION (NES = −2.04, q-value = .000178). Similarly, in AS, MAB21L2 was enriched in the DNA REPLICATION pathway (NES = 2.632, q-value = .000).

Figure 6.

(A) ROC curve illustrating the performance of the shared diagnostic genes in GSE113212. (B) ROC curve displaying the performance of the shared diagnostic genes in GSE181716. (C) ROC curve demonstrating the performance of the shared diagnostic genes in GSE73754. (D) ROC curve showcasing the performance of the shared diagnostic genes in GSE153310. (E) Expression profile of MAB21L2 in GSE106253. (F) Expression profile of MAB21L2 in GSE134290. (G) Results of single-gene GSEA analyses for DEGs in GSE106253. (H) Results of single-gene GSEA analyses for DEGs in GSE134290.

Immune Cell Infiltration and its correlation with Candidate Biomarker

We delved deeper into the variation in immune cell infiltration among different samples. The landscape of immune cell infiltration in OLF and AS was assessed through the ssGSEA algorithm. The Wilcoxon test unveiled significant alterations in 6 immune subsets between OLF and normal samples. These included activated B cells, central memory CD4 T cells, memory B cells, immature dendritic cells, natural killer cells, and plasmacytoid dendritic cells (Figure 7A and B). Additionally, between AS and normal samples, 3 immune subtypes showed significant differences: activated CD8 T cells, T follicular helper cells, and natural killer T cells (Figure 7C and D). Furthermore, our analysis revealed specific correlations between immune subsets and MAB21L2. In OLF samples, 4 immune subsets, namely immature dendritic cells, natural killer cells, plasmacytoid dendritic cells, and central memory CD4 T cells, displayed positive correlations with MAB21L2, while activated CD8 T cells, memory B cells, and activated B cells exhibited negative correlations with MAB21L2 (Figure 7E). In AS samples, MAB21L2 exhibited positive associations with activated CD8 T cells, T follicular helper cells, immature dendritic cells, and natural killer T cells but a negative association with central memory CD8 T cells (Figure 7F). These findings suggest that MAB21L2 may play a role in B cell and T cell immune responses in both OLF and AS.

Figure 7.

(A) Bar plot and (B) violin plot displaying the distribution of 28 immune cells in OLF. (C) Bar plot and (D) violin plot illustrating the distribution of 28 immune cells in AS. (E) Correlation analysis showing the relationship between immune infiltrating cells and MAB21L2 in OLF. (F) Correlation analysis demonstrating the relationship between immune infiltrating cells and MAB21L2 in AS.

Discussion

Our study represents an integration of transcriptomes from OLF and AS, employing WGCNA for the first time to uncover a shared mechanism between the 2 conditions. This analysis unveiled potential crosstalk genes, common pathways, and associated immune cell interactions. Notably, GO and KEGG analysis uncovered a commonality in which DEGs in both OLF and AS predominantly centered around the ECM-receptor interaction and cell adhesion. The crosstalk between cells and their surrounding ECM and cell-cell adhesion activity play an important role during cartilage development and the establishment of the skeletal framework during endochondral ossification.²⁹ Specific chemical and mechanical signals are transmitted within the ECM and influence many aspects of cell behavior such as matrix synthesis, migration, proliferation, survival, and differentiation.³⁰ In addition, mesenchymal cells migrate toward and condense at the sites of future skeletal elements during endochondral ossification. This condensation results in an increased cell density and increased number of cellcell contacts accompanied by the expression of molecules involved in cell-cell adhesion such as N-cadherin.^31,32 Through the intersection of key genes identified in the WGCNA modules and DEGs, we pinpointed MAB21L2, MEGF10, and ISLR as the pivotal genes mediating the crosstalk between OLF and AS. Further GO BP (biological process) annotation analysis revealed that these central genes were predominantly enriched in processes related to myoblast proliferation, myoblast migration, skeletal muscle satellite cell proliferation, skeletal muscle cell proliferation. Consequently, we hypothesize that the convergent ectopic ossification pathways observed in OLF and AS could be intricately associated with pathological bone formation in the ligamentum flavum and adjacent muscular tissues at its sites of attachment. Heterotopic ossification is a diverse pathologic process, defined as the formation of extraskeletal bone in muscle and soft tissues. Research has shown that Bone morphogenetic protein (BMP) signaling is intricately involved in the processes leading to heterotopic ossification.³³ The patterns of BMPs and BMP receptors expression in individuals with OLF differ from those with nonossified ligamentum flavum (LF).³⁴ Notably, in human LF cells, BMP-2 exhibits a significant capacity to elevate the expression of osteogenic traits and incite the formation of osseous nodules.³⁵ This suggests that BMP might play a role in fostering endochondral ossification at heterotopic ossification sites in OLF.³⁶ Histologic analysis suggests that active BMP signaling is found in the entheses of patients with peripheral SpA and enthesophyte formation.³⁷ Furthermore, recent investigations have uncovered that aberrant chondrocytes promoted ectopic new bone formation in ankylosing spondylitis through BMP6/pSmad1/5 signaling.⁷ To conclude, it is evident that the BMP signaling pathway is integral to the ectopic ossification processes in OLF and AS. Myoblasts (proliferating satellite cells), a population of muscle stem cells, proliferate, differentiate, and fuse into multinuclear myofibers after myofibers are damaged, allowing damaged muscle to be repaired,²⁸ and myoblasts are known to originate from muscle satellite cells.³⁸ While the molecular mechanisms governing their commitment to differentiation remain poorly understood, there is a significant body of research focused on the role of resident stem cells within skeletal muscle in the generation of ectopic bone tissue. Research has indicated that BMP2 appears to have the capability to divert committed myogenic cells away from the myogenic lineage’s differentiation pathway, inducing them to differentiate into the osteogenic lineage.³⁹ Additionally, studies have revealed that sustained expression of Foxc2 in myogenic cells inhibits muscle generation and alters their lineage commitment towards osteogenic differentiation by inducing the Wnt4 and Bmp4 signaling pathways.³⁸ Therefore, we hypothesize that these central genes may be associated with the ectopic ossification of myoblasts in OLF and AS.

In our further study, ROC analysis revealed that MAB21L2 exhibited an AUC exceeding .600 for predicting OLF and AS. As such, MAB21L2 emerged as a diagnostic marker with robust predictive capabilities. Remarkably, the results of the differential expression analysis unveiled a significant downregulation in OLF, whereas it exhibited a notable upregulation of MAB21L2 in cases of AS. At present, investigations into the MAB21L2 gene remain sparse in the scholarly literature. However, through a review of existing literature, we have found evidence that MAB21L2 immunoprecipitates in vivo with the BMP4 effector SMAD1. In vitro, it binds to SMAD1 and the SMAD1-SMAD4 complex. This suggests that MAB21L2 may serve as a repressor of BMP-induced transcription.⁴⁰ BMPs are a group of growth factors known for their ability to induce bone and cartilage development and are a member of the transforming growth factor-beta (TGF-β) superfamily. Activation of BMP receptors results in the initiation of transcription via either the MAP kinase signaling pathway or the phosphorylation of SMAD1/5/8 proteins.^41,42 BMP signaling pathways intricately govern various facets of endochondral bone formation. Studies have provided evidence that inhibitory Smads possess the capability to curtail BMP signaling, thereby modulating BMP-induced impacts on chondrocytes. BMPs are versatile growth factors that play a constructive role in stimulating the proliferation and differentiation of osteoblasts.⁴³ Furthermore, BMP such as BMP2, BMP4, and BMP9 have been shown to be the primary inducers of heterotopic ossification when injected into the muscle tissue of mice^44-46 and are highly expressed in human heterotopic ossification lesions.^47,48 Previous studies have found that BMP is strongly expressed in chondrocytes within the calcified areas in most cases of OLF.^49-51 Therefore, we speculate that the decreased expression level of MAB21L2 in OLF is associated with the ectopic ossification process in OLF. The pathogenesis of AS is not entirely known, and research suggests that skeletal changes in AS extend beyond the characteristic new bone formation leading to ankylosis, including bone loss.⁵² Additionally, studies have confirmed the presence of 3 different forms of ectopic ossification in AS,⁵³ including endochondral ossification, direct bone formation through osteoblast differentiation, and cartilaginous ossification. Furthermore, different BMP signaling pathways and BMPs may be involved at various stages of AS disease development,³⁷ potentially leading to differential expression levels of BMP4 during different processes. Therefore, we speculate that the increased expression level of MAB21L2 in AS may be related to these factors.

Further single-gene GSEA shed light on the potential functions of MAB21L2 in OLF and AS. Notably, these functions encompassed lysosome-related signaling and oxidative phosphorylation. These fundamental cellular processes are intricately connected to the phenomenon of ectopic ossification.^54-56 Additionally, prior research has indicated the presence of immune responses and inflammatory activities in both OLF and AS, involving a variety of immune cell types.⁵⁷ Moreover, local inflammation is believed to play a pivotal role in various phenotypes of ectopic ossification and the bone remodeling process. In this study, we identified differentially expressed immune cells, encompassing both the initial innate immunity and subsequent adaptive immunity, in order to unravel the pathogenesis of OLF and AS. The results of immune infiltration unveiled a higher proportion of T cells in the pathogenesis of both OLF and AS. Prior investigations have documented the emergence and critical contribution of CD4⁺ and CD8⁺ T cells in the inflammatory response at human tendon attachment sites.⁵⁸ Additional research has highlighted the critical importance of the Notch signaling pathway in both the pathological ossification processes within ligamentum flavum tissues and cells,⁵⁹ and in the differentiation and functional modulation of T cells.⁶⁰ Furthermore, it is well-established that the immune system exerts strong paracrine effects that predominantly influence cellular activities associated with ossification. This intricate interplay between the skeletal system and the immune system has ushered in a new avenue of research, referred to as osteoimmunology.⁶¹ Nevertheless, the interaction between adaptive immunity and MAB21L2 necessitates further exploration. We found that MAB21L2 had a significant correlation with activated CD8 T cells in both OLF samples and AS samples. BMP4 has been reported to be able to regulate T cell development in the thymus. In vitro experiments have shown that treatment of thymus explants with exogenous BMP4 negatively regulated differentiation of early thymocyte progenitors and the transition from CD4−CD8− cells to CD4⁺CD8⁺ cells.⁶² As previously mentioned, MAB21L2 acts as a transcriptional inhibitor of BMP4⁴⁰; hence, we speculate that BMP4 may serve as a pivotal link between MAB21L2 and T cell activities. This will highlight the common pathophysiology at the immune level and may be a key factor in understanding the relationship between OLF and AS. However, the interaction between adaptive immunity and MAB21L2 needs to be further explored.

It’s important to acknowledge certain limitations in this research. Firstly, despite our inclusion of all accessible datasets from the GEO database, the volume of data remains restricted. Further investigation should encompass additional multicenter external datasets and involve extensive experimental studies to delve into the precise molecular mechanisms of the hub genes. Secondly, since MAB21L2 does not possess typical receptor characteristics but is primarily involved in cellular development processes,⁶³ there are currently no reports in the literature of specific receptor antagonists or agonists for MAB21L2. However, with ongoing research into the role of this gene in ectopic ossification, future studies may explore key downstream effectors or signaling pathways associated with MAB21L2, aiming to develop targeted small molecule inhibitors or agonists. Such approaches could indirectly modulate the biological functions of MAB21L2, thereby offering new therapeutic strategies for the clinical treatment of ectopic ossification. Lastly, the hypothesis articulated here necessitates substantial empirical backing. Future studies should encompass a variety of experimental techniques, including western blotting, immunohistochemistry assays, and quantitative reverse transcription-PCR, to validate our findings. In light of these constraints, this study aspires to offer an initial perspective on shared biomarkers associated with ectopic ossification in spinal ligaments. In forthcoming research, we plan to delve deeper into the alterations in MAB21L2 expression, utilizing relevant in vitro and in vivo experiments to enhance our understanding.

Conclusions

The current research represents the inaugural attempt to employ bioinformatic analysis to scrutinize the shared pathways and genetic diagnostic markers linked to OLF and AS. The findings suggest a conceivable association between the ECM-receptor interaction, cell-cell adhesion, and the pathogenesis of both OLF and AS. Furthermore, MAB21L2 emerges as a potential diagnostic marker for individuals facing the challenges of OLF concurrent with AS. The correlation analysis of immune infiltration implies a close connection between the pathogenesis of OLF and AS and T cells. In essence, this study offers a fresh perspective for delving into the possible mechanisms of AS secondary to OLF.

Footnotes

Authors’ Contributions

Conceptualization, methodology, Y.L. and Y.L.; software, Y.L. and Y.L.; validation, Y.L., Z.G. and J.Z.; data curation, C.W. and Z.G.; writing—original draft preparation, Y.L. and Y.Q.; writing—review and editing, J.Y. and X.L.; visualization, J.Y.; supervision, X.L. All authors have read and agreed to the published version of the manuscript.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Yishan Liu

Data Availability Statement

The data used in this study were sourced from the publicly accessible GEO database at .

References

Reveille

Weisman

. The epidemiology of back pain, axial spondyloarthritis and HLA-B27 in the United States. Am J Med Sci. 2013;345:431-436. doi:10.1097/maj.0b013e318294457f.

Sieper

Poddubnyy

. Axial spondyloarthritis. Lancet. 2017;390:73-84. doi:10.1016/s0140-6736(16)31591-4.

Akkoc

Khan

. Overestimation of the prevalence of ankylosing spondylitis in the Berlin study: comment on the article by Braun et al. Arthritis Rheum. 2005;52:4048-4049, Author Reply 4049-4050. doi:10.1002/art.21492.

Liao

, et al. Epidemiology of spondyloarthritis in the People's Republic of China: review of the literature and commentary. Semin Arthritis Rheum. 2007;37:39-47. doi:10.1016/j.semarthrit.2007.01.003.

Ward

Deodhar

Akl

, et al. American college of rheumatology/spondylitis association of America/spondyloarthritis research and treatment network 2015 Recommendations for the treatment of ankylosing spondylitis and nonradiographic axial spondyloarthritis. Arthritis Rheumatol. 2016;68:282-298. doi:10.1002/art.39298.

Baeten

Sieper

Braun

, et al. Secukinumab, an interleukin-17A inhibitor, in ankylosing spondylitis. N Engl J Med. 2015;373:2534-2548. doi:10.1056/NEJMoa1505066.

Shao

Liu

Zhu

, et al. Targeting chondrocytes for arresting bony fusion in ankylosing spondylitis. Nat Commun. 2021;12:6540. doi:10.1038/s41467-021-26750-6.

Kawaguchi

Nakano

Yasuda

, et al. Characteristics of ossification of the spinal ligament; incidence of ossification of the ligamentum flavum in patients with cervical ossification of the posterior longitudinal ligament - analysis of the whole spine using multidetector CT. J Orthop Sci. 2016;21:439-445. doi:10.1016/j.jos.2016.04.009.

Guo

Luk

Karppinen

, et al. Prevalence, distribution, and morphology of ossification of the ligamentum flavum: a population study of one thousand seven hundred thirty-six magnetic resonance imaging scans. Spine (Phila Pa 1976). 2010;35:51-56. doi:10.1097/BRS.0b013e3181b3f779.

10.

Aizawa

Sato

Sasaki

, et al. Thoracic myelopathy caused by ossification of the ligamentum flavum: clinical features and surgical results in the Japanese population. J Neurosurg Spine. 2006;5:514-519. doi:10.3171/spi.2006.5.6.514.

11.

Tam

. Pathogenesis of ankylosing spondylitis. Nat Rev Rheumatol. 2010;6:399-405. doi:10.1038/nrrheum.2010.79.

12.

Hermann

Bollow

. Magnetic resonance imaging of the axial skeleton in rheumatoid disease. Best Pract Res Clin Rheumatol. 2004;18:881-907. doi:10.1016/j.berh.2004.06.005.

13.

Ono

Yonenobu

Miyamoto

, et al. Pathology of ossification of the posterior longitudinal ligament and ligamentum flavum. Clin Orthop Relat Res. 1999:18-26. doi:10.1097/00003086-199902000-00003.

14.

Fan

Meng

Sun

, et al. Genome-wide DNA methylation profile analysis in thoracic ossification of the ligamentum flavum. J Cell Mol Med. 2020;24:8753-8762. doi:10.1111/jcmm.15509.

15.

Chen

, et al. Role of Cx43-Mediated NFкB signaling pathway in ossification of posterior longitudinal ligament: an in vivo and in vitro study. Spine (Phila Pa 1976). 2017;42:E1334-e1341. doi:10.1097/brs.0000000000002165.

16.

Baraliakos

Listing

Rudwaleit

, et al. The relationship between inflammation and new bone formation in patients with ankylosing spondylitis. Arthritis Res Ther. 2008;10:R104. doi:10.1186/ar2496.

17.

Wang

Zhan

, et al. Inflammation intensity-dependent expression of osteoinductive Wnt proteins is critical for ectopic new bone formation in ankylosing spondylitis. Arthritis Rheumatol. 2018;70:1056-1070. doi:10.1002/art.40468.

18.

van der Heijde

Machado

Braun

, et al. MRI inflammation at the vertebral unit only marginally predicts new syndesmophyte formation: a multilevel analysis in patients with ankylosing spondylitis. Ann Rheum Dis. 2012;71:369-373. doi:10.1136/annrheumdis-2011-200208.

19.

Chiowchanwisawakit

Lambert

Conner-Spady

, et al. Focal fat lesions at vertebral corners on magnetic resonance imaging predict the development of new syndesmophytes in ankylosing spondylitis. Arthritis Rheum. 2011;63:2215-2225. doi:10.1002/art.30393.

20.

Hijioka

Suzuki

Nakamura

, et al. Light and electron microscopy of hydroxyapatite depositions in the ligamentum flavum. Spine (Phila Pa 1976). 1994;19:2626-2631.

21.

Zhang

Zhu

, et al. Chondrogenesis mediates progression of ankylosing spondylitis through heterotopic ossification. Bone Res. 2021;9:19. doi:10.1038/s41413-021-00140-6.

22.

Yao

Zhang

Gao

, et al. Exploration of the shared gene Signatures and molecular mechanisms between systemic lupus erythematosus and pulmonary arterial hypertension: evidence from transcriptome data. Front Immunol. 2021;12:658341. doi:10.3389/fimmu.2021.658341.

23.

Zhang

Horvath

. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4:17. doi:10.2202/1544-6115.1128.

24.

Wang

Han

, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16:284-287. doi:10.1089/omi.2011.0118.

25.

Robin

Turck

Hainard

, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinf. 2011;12:77. doi:10.1186/1471-2105-12-77.

26.

Friedman

Hastie

Tibshirani

. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33:1-22.

27.

Hänzelmann

Castelo

Guinney

. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;14:7. doi:10.1186/1471-2105-14-7.

28.

, et al. BMP9 functions as a negative regulator in the myogenic differentiation of primary mouse myoblasts. Biosci Biotechnol Biochem. 2023;87:1255-1264. doi:10.1093/bbb/zbad104.

29.

Prein

Beier

. ECM signaling in cartilage development and endochondral ossification. Curr Top Dev Biol. 2019;133:25-47. doi:10.1016/bs.ctdb.2018.11.003.

30.

Geiger

Bershadsky

Pankov

, et al. Transmembrane crosstalk between the extracellular matrix--cytoskeleton crosstalk. Nat Rev Mol Cell Biol. 2001;2:793-805. doi:10.1038/35099066.

31.

Delise

Tuan

. Analysis of N-cadherin function in limb mesenchymal chondrogenesis in vitro. Dev Dyn. 2002;225:195-204. doi:10.1002/dvdy.10151.

32.

Olsen

Reginato

Wang

. Bone development. Annu Rev Cell Dev Biol. 2000;16:191-220. doi:10.1146/annurev.cellbio.16.1.191.

33.

Agarwal

Loder

Breuler

, et al. Strategic targeting of multiple BMP receptors prevents trauma-induced heterotopic ossification. Mol Ther. 2017;25:1974-1987. doi:10.1016/j.ymthe.2017.01.008.

34.

Hayashi

Ishidou

Yonemori

, et al. Expression and localization of bone morphogenetic proteins (BMPs) and BMP receptors in ossification of the ligamentum flavum. Bone 1997: pages 23-30. doi: 10.1016/s8756-3282(97)00080-x

35.

Moon

Park

Kim

, et al. Biologic modification of ligamentum flavum cells by marker gene transfer and recombinant human bone morphogenetic protein-2. Spine (Phila Pa 1976). 2004;29:960-965. doi:10.1097/00007632-200405010-00003.

36.

Chen

Fan

, et al. Two novel BMP-2 variants identified in patients with thoracic ossification of the ligamentum flavum. Eur J Hum Genet. 2017;25:565-571. doi:10.1038/ejhg.2017.2.

37.

Lories

Derese

Luyten

. Modulation of bone morphogenetic protein signaling inhibits the onset and progression of ankylosing enthesitis. J Clin Invest. 2005;115:1571-1579. doi:10.1172/jci23738.

38.

Gozo

Aspuria

Cheon

, et al. Foxc2 induces Wnt4 and Bmp4 expression during muscle regeneration and osteogenesis. Cell Death Differ. 2013;20:1031-1042. doi:10.1038/cdd.2013.34.

39.

Katagiri

Yamaguchi

Komaki

, et al. Bone morphogenetic protein-2 converts the differentiation pathway of C2C12 myoblasts into the osteoblast lineage. J Cell Biol. 1994;127:1755-1766. doi:10.1083/jcb.127.6.1755.

40.

Baldessari

Badaloni

Longhi

, et al. MAB21L2, a vertebrate member of the Male-abnormal 21 family, modulates BMP signaling and interacts with SMAD1. BMC Cell Biol. 2004;5:48. doi:10.1186/1471-2121-5-48.

41.

Rawadi

Vayssière

Dunn

, et al. BMP-2 controls alkaline phosphatase expression and osteoblast mineralization by a Wnt autocrine loop. J Bone Miner Res. 2003;18:1842-1853. doi:10.1359/jbmr.2003.18.10.1842.

42.

Canalis

. Growth factor control of bone mass. J Cell Biochem. 2009;108:769-777. doi:10.1002/jcb.22322.

43.

Matsumoto

Toyama

Chikuda

, et al. Outcomes of fusion surgery for ossification of the posterior longitudinal ligament of the thoracic spine: a multicenter retrospective survey: clinical article. J Neurosurg Spine. 2011;15:380-385. doi:10.3171/2011.6.Spine10816.

44.

Lounev

Ramachandran

Wosczyna

, et al. Identification of progenitor cells that contribute to heterotopic skeletogenesis. J Bone Joint Surg Am. 2009;91:652-663. doi:10.2106/jbjs.H.01177.

45.

Medici

Shore

Lounev

, et al. Conversion of vascular endothelial cells into multipotent stem-like cells. Nat Med. 2010;16:1400-1406. doi:10.1038/nm.2252.

46.

Leblanc

Trensz

Haroun

, et al. BMP-9-induced muscle heterotopic ossification requires changes to the skeletal muscle microenvironment. J Bone Miner Res. 2011;26:1166-1177. doi:10.1002/jbmr.311.

47.

Gannon

Kaplan

Olmsted

, et al. Bone morphogenetic protein 2/4 in early fibromatous lesions of fibrodysplasia ossificans progressiva. Hum Pathol. 1997;28:339-343. doi:10.1016/s0046-8177(97)90133-7.

48.

Grenier

Leblanc

Faucheux

, et al. BMP-9 expression in human traumatic heterotopic ossification: a case report. Skelet Muscle. 2013;3:29. doi:10.1186/2044-5040-3-29.

49.

Miyamoto

Takaoka

Yonenobu

, et al. Ossification of the ligamentum flavum induced by bone morphogenetic protein. An experimental study in mice. J Bone Joint Surg Br. 1992;74:279-283. doi:10.1302/0301-620x.74b2.1544970.

50.

Hoshi

Amizuka

Sakou

, et al. Fibroblasts of spinal ligaments pathologically differentiate into chondrocytes induced by recombinant human bone morphogenetic protein-2: morphological examinations for ossification of spinal ligaments. Bone. 1997;21:155-162. doi:10.1016/s8756-3282(97)00106-3.

51.

Mimatsu

Kishi

Hashizume

. Experimental chronic compression on the spinal cord of the rabbit by ectopic bone formation in the ligamentum flavum with bone morphogenetic protein. Spinal Cord. 1997;35:740-746. doi:10.1038/sj.sc.3100533.

52.

Carter

Lories

. Osteoporosis: a paradox in ankylosing spondylitis. Curr Osteoporos Rep. 2011;9:112-115. doi:10.1007/s11914-011-0058-z.

53.

Benjamin

Toumi

Suzuki

, et al. Evidence for a distinctive pattern of bone formation in enthesophytes. Ann Rheum Dis. 2009;68:1003-1010. doi:10.1136/ard.2008.091074.

54.

Lawrence

Zoncu

. The lysosome as a cellular centre for signalling, metabolism and quality control. Nat Cell Biol. 2019;21:133-142. doi:10.1038/s41556-018-0244-7.

55.

Xiao

Wang

, et al. The interaction between intracellular energy metabolism and signaling pathways during osteogenesis. Front Mol Biosci. 2021;8:807487. doi:10.3389/fmolb.2021.807487.

56.

Salhotra

Shah

Levi

, et al. Mechanisms of bone development and repair. Nat Rev Mol Cell Biol. 2020;21:696-711. doi:10.1038/s41580-020-00279-w.

57.

Lee

Byun

Madhurakkat Perikamana

, et al. Current advances in immunomodulatory biomaterials for bone regeneration. Adv Healthc Mater. 2019;8:e1801106. doi:10.1002/adhm.201801106.

58.

Watad

Rowe

Russell

, et al. Normal human enthesis harbours conventional CD4+ and CD8+ T cells with regulatory features and inducible IL-17A and TNF expression. Ann Rheum Dis. 2020;79:1044-1054. doi:10.1136/annrheumdis-2020-217309.

59.

Chen

Fan

, et al. Notch signaling pathways in human thoracic ossification of the ligamentum flavum. J Orthop Res. 2016;34:1481-1491. doi:10.1002/jor.23303.

60.

Shah

Zúñiga-Pflücker

. Notch receptor-ligand interactions during T cell development, a ligand endocytosis-driven mechanism. Curr Top Microbiol Immunol. 2012;360:19-46. doi:10.1007/82_2012_225.

61.

Schlundt

Schell

Goodman

, et al. Immune modulation as a therapeutic strategy in bone regeneration. J Exp Orthop. 2015;2:1. doi:10.1186/s40634-014-0017-6.

62.

Hager-Theodorides

Ross

Sahni

, et al. Direct BMP2/4 signaling through BMP receptor IA regulates fetal thymocyte progenitor homeostasis and differentiation to CD4+CD8+ double-positive cell. Cell Cycle. 2014;13:324-333. doi:10.4161/cc.27118.

63.

Huang

Xiang

Zhou

, et al. The male abnormal gene family 21 (Mab21) members regulate Eye development. Curr Mol Med. 2016;16:660-667. doi:10.2174/1566524016666160824150729.