Integrating Bioinformatic Strategies with Real-World Data to Infer Distinctive Immunocyte Infiltration Landscape and Immunologically Relevant Transcriptome Fingerprints in Ossification of Ligamentum Flavum.

PURPOSE: Ossification of the ligamentum flavum (OLF) is a multifactorial disease characterized by an insidious and debilitating process of abnormal bone formation in ligamentum tissues. However, its definite pathogenesis has not been fully elucidated. Potential links between the immune system and various forms of heterotopic ossification have been discussed for many years, whereas no research investigated the immune effects on the initiation and development of OLF. Therefore, we attempt to shed light on this issue.
METHODS: A series of bioinformatic algorithms were integrated to evaluate the immune score and the immunocyte infiltration patterns between OLF and normal samples, screen OLF-related and immune-related differentially expressed genes (OIDEGs), and analyze their biological functions. Correlation analysis inferred OIDEGs-related differentially expressed lncRNAs (OIDELs) and infiltrating immune cells (OIICs) to construct an immunoregulatory network.
RESULTS: Differential immune score and immune cell infiltration were determined between two groups, and 10 OIDEGs with diverse biological function annotations were identified and verified. A lncRNA-gene-immunocyte regulatory network further revealed 10 OIDEGs, 41 OIDELs and 7 OIICs that were highly correlated. Among them, CD1E and STAT3 were predicted as hub genes whether at the expression level or interaction level. cDCs emerged as having the most prominent differences and the highest degree of connectivity. FO393414.3, AC096734.1, LINC01137 and DLX6-AS1 with the greatest number of OIDEGs were thought to be more likely to participate in immunoregulation of OLF.
CONCLUSION: This is the first research to preliminarily elucidate OLF-related immunocyte infiltration landscape and immune-associated transcriptome signatures based on bioinformatic strategies and real-world data, which may provide compelling insights into the pathogenesis and therapeutic targets of OLF.

Introduction

Ossification of ligamentum flavum (OLF) is characterized by a progressive pathological process of heterotopic bone formation from ligamentous tissues that causes spinal canal stenosis and severe myeloradiculopathy without response to conservative managements.1,2 OLF predominately affects East Asian population, with an estimated prevalence of 6.1% in Chinese, 12% in Japanese and 21.8% in Korean populations.3–5 The onset of OLF is so insidious that early detection is difficult. At the point when OLF is diagnosed, patients usually develop notable neurological symptoms, and have to require challenging surgical interventions with multiple complications.6 Therefore, it is essential to explore potential diagnostic biomarkers and therapeutic targets of OLF based on scientific research of its molecular mechanisms.

It is generally believed that genetic and environmental factors can contribute to the occurrence and development of OLF.7 However, its definite pathomechanism remains inadequately understood. Histologically, its pathological nature is considered as endochondral osteogenesis involving the differentiation of ligament fibroblasts into osteoblasts.7 Therefore, abnormally expressed osteogenic-related genes in ligamentum flavum (LF) cells are indispensable contributors to the pathogenesis of OLF. Multiple osteogenic markers (ALP, RUNX2, OCN, OPN) and signaling pathway molecules (BMP2, Wnt/β-catenin, Notch) have been identified to be significantly up-regulated in OLF.8,9 In addition, potential causative factors including mechanical stressors, metabolic elements, inflammatory markers and angiogenic effects can influence the osteogenic differentiation of LF cells by regulating these target genes.10–13 Transcription factors and epigenetic modification can also be involved in the development of OLF by endogenously altering the expression of osteogenic targets.14,15 With the development of second-generation sequencing techniques, emerging RNA sequencing and bioinformatics analysis have revealed differential OLF-related transcriptome expression profiles and functional characteristics, including identification of mRNAs, miRNAs, lncRNAs, circRNAs and construction of ceRNA network.16–18 However, there is lack of specialized investigation on the immune-associated transcriptome signatures in OLF.

Osteoimmunology, immunoregulatory effects on bone homeostasis and bone metabolism, has gained substantial interest in recent years. Accumulated researches have revealed immune cells and their inflammatory factors secreted (eg, IL-3, IL-6, TNF) participating in multiple forms of heterotopic ossification (HO) disorders, such as ankylosing spondylitis (AS), fibrodysplasia ossificans progressiva (FOP), neurogenic heterotopic ossification (NHO) and vascular calcification.19–21 Moreover, a growing body of evidence has uncovered how the inflammatory response functions in the formation of HO, and several specific cell populations in the innate and adaptive immunity, such as neutrophils, macrophages, mast cells, B cells and T cells, have been particularly implicated in the pathophysiology of HO via interactions with osteoprogenitor cells and the release of osteogenic growth factors.20,22 Historically, differential responses of peripheral lymphocytes were observed among patients with the continuous-type ossification of posterior longitudinal ligament (OPLL), those with the segmental-type OPLL and healthy volunteers.23 Moreover, Saito et al demonstrated that infiltrating macrophages could contribute to the progression of LF hypertrophy by stimulating collagen production in fibroblasts.24 Inflammatory response plays an important role in the immune system, and is also involved in the occurrence and development of OLF.25 TNF-α, a pro-inflammatory cytokine, is mainly produced by activated macrophages, which has been supposed to make for new bone formation.26 Consistently, our previous iTRAQ proteomics study found that the expression level of TNF-α was significantly elevated in OLF samples compared to normal controls, and TNF-α facilitated osteogenic differentiation of LF fibroblasts through activating Osx expression.12 These evidences suggested immune response might be associated with the pathomechanism of OLF. However, no investigation has been performed to shed light on this issue.

To systematically investigate the possible role of immune regulation in OLF, we applied multiple bioinformatic algorithms to evaluate difference of immune score and immune cell infiltration between OLF and normal tissues for the first time. Subsequently, we identified and validated the hub differentially expressed OLF-related and immune-related genes and illustrated their correlative biological functions and signaling pathways based on the real-world data. Finally, we conducted a joint correlation analysis of hub genes, lncRNAs and infiltrating immune cells to predict a potential immunoregulatory network in OLF. These findings preliminarily inferred the potential contribution of an immune-directed mechanism in the occurrence and development of OLF, providing compelling insights into the detailed pathogenesis and novel therapeutic targets of OLF as well as available data resources for future research.

Materials and Methods

Data Collection and Preprocessing

The workflow diagram of this study is described in Figure 1. Gene expression profile data were downloaded from Gene Expression Omnibus (GEO) database (). Eligible GEO datasets were selected according to the following inclusion criteria: 1) organism: Homo sapiens; 2) expression profiling by microarray; 3) samples: OLF ligament tissues and normal ligament tissues. Datasets from other organisms or expression profiling by RT-PCR or genome variation profiling were excluded. Ultimately, the lncRNA/mRNA expression dataset GES106253 was obtained, and clinical data related to these samples was also collected from the published study.16 Ossified ligament samples and normal samples were acquired from 4 TOLF patients and 4 spinal trauma patients during surgery, respectively. Gene expression analysis was based on GPL21827 platform (Agilent-079487 Arraystar Human LncRNA microarray V4), and expression values were normalized with robust multi-array average. The immune-related gene list was obtained from ImmPort database (). The macrophage-related gene list was retrieved from Genecards database ().

Figure 1

The overview of the analysis procedure.

Determination of Immune and Macrophage Score

To investigate whether immune response is implicated in the pathomechanism of OLF, an immune score and a macrophage score for each sample were calculated based on the immune-related gene list and macrophage-related gene list by applying the single sample gene set enrichment analysis (ssGSEA) algorithm.

Evaluation of Immune Infiltration Landscape

Principal component analysis (PCA) was performed to determine whether there was a difference in immune cell infiltration between OLF tissues and normal tissues. ssGSEA was introduced to calculate the relative infiltration scores of different immune cell types based on the expression of reference gene within the gene sets from RNA-seq data. The enrichment score in ssGSEA represented the relative abundance of each immune cell type. xCell, a novel gene signature-based method reliably portraying the cellular heterogeneity landscape of tissue expression profiles, was utilized to infer the estimated proportion of 64 types of immune cells. The cut-off values for the cell analyses were p < 0.05. Heat maps and box diagrams were plotted to visualize the immune infiltration landscape.

Differential Analysis of Infiltrating Immunocyte Types

The significant differential infiltrating immune cells were screened according to p value < 0.05, and the percentage of each type of immune cell in the samples was calculated. Pearson correlation coefficient was used to further understand the relationship between these OLF-related infiltrating immune cells (OIICs).

Identification of Immune-Related Genes and Correlated LncRNAs

To determine differentially expressed genes (DEGs) and differentially expressed lncRNAs (DELs), GEO2R (), an interactive online tool, was engaged. Based on the predetermined statistical threshold of |fold change| >1.5 and p-values < 0.05, DEGs and DELs were screened out independently and utilized for further analysis. Furthermore, OLF-related DEGs and immune-related genes were intersected to obtain the OLF-related and immune-related DEGs (OIDEGs). Subsequently, Pearson correlation analysis (r > 0.9 or < −0.9, p < 0.01) between OIDEGs and DELs was conducted to confirm OIDEG-related DELs (OIDELs). Two-way hierarchical clustering was performed using pheatmap R package to render a heatmap of DEGs, DELs, OIDEGs and OIDELs.

Functional Enrichment Analyses of OIDEGs

Gene ontology (GO) analysis was performed to annotate OIDEGs, and to illustrate their functions in biology process, cell component, and molecular function. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was conducted to investigate the related signaling pathways, and to understand high-level and biological functions from large-scale molecular datasets. The clusterProfiler package was used for the GO and KEGG analysis. The cutoff values for the GO and KEGG analyses were set at p < 0.05.

Correlation Analysis Between OIICs and OIDEGs

In order to further reveal the potential mechanism for the immune dysfunction in OLF, Pearson correlation analysis was used to analyze the relationship between OIDEGs and OIICs in R software. P-value < 0.05 was used as the cut-off standard.

Construction of LncRNA-Gene-Immunocyte Regulatory Networks

Based on above correlation analysis results, the OIDEG-OIDEL correlation pairs with the value of r > 0.95 or < −0.95 together with p < 0.001 and the significant OIDEG-OIIC correlation pairs with the value of r > 0.90 or < −0.90 together with p < 0.001 were integrated to construct a potential regulatory network, and to screen out more valuable regulators.

Validation of Candidate OIDEGs by Quantitative Real-Time PCR

Real-world data were collected to validate the expression changes of above candidate OIDEGs. The study protocol was approved by the Ethics Committee for Human Subjects of the Peking University Third Hospital in accordance with the Declaration of Helsinki (PUTH-REC-SOP-06-3.0-A27, #2014003). Since the homeostasis and function of immune system can be influenced by multiple factors, such as age and other comorbid diseases. Therefore, the rigorous inclusion and exclusion criteria were set to ensure the reliability of the results as far as possible. Experimental group included 10 adult patients diagnosed with TOLF who underwent surgical decompression, and control group included 10 adult patients diagnosed with thoracic disc herniation or spine fracture from June 2020 to December 2020. Patients with ankylosing spondylitis (AS), diffuse idiopathic skeletal hyperostosis (DISH), rheumatoid arthritis, spinal tumors, spinal infections and other systemic autoimmune diseases, were excluded. In addition, patients with hypertension, diabetes and arteriosclerosis were also eliminated. The detailed patients’ clinical data are described in Table 1, which showed there were no significant difference in age, gender and BMI between two groups. Their ossified and normal ligament tissue specimens were obtained during surgery with informed consent. Total RNA from pretreated ligament tissues was extracted with TRIzol reagent (German DBI) for qRT-PCR. Reverse transcription of RNA into complementary DNA was performed using the Bestar qPCR RT Kit (DBI, Germany) following the manufacturer’s instructions. GAPDH was used as the internal control, and the expression of candidate genes was calculated by the 2-ΔΔCt method.

Table 1

Basic Clinical Data of Included Patients from OLF Group and Normal Group

OLF Group					Normal Group
Patient Number	Diagnosis (Location)	Gender	Age (Years)	BMI (kg/m2)	Patient Number	Diagnosis (Location)	Gender	Age (Years)	BMI (kg/m2)
1	TOLF (T10–11)	Male	53	23.88	1	TDH (T11–12)	Female	65	26.35
2	TOLF (T11–12)	Male	64	26.87	2	TF (T6–9)	Male	70	29.41
3	TOLF (T10–12)	Female	60	28.52	3	TDH (T10–12)	Male	62	31.70
4	TOLF (T11–12)	Male	73	24.77	4	TF (T12)	Male	48	19.53
5	TOLF (T8–12)	Male	62	23.88	5	TDH (T10–12)	Male	76	25.97
6	TOLF (T1–3, T6–9)	Male	63	31.02	6	TDH (T9–11)	Male	53	28.41
7	TOLF (T1–5)	Female	53	25.25	7	TDH (T11–12)	Female	61	30.44
8	TOLF (T10–12)	Male	64	33.67	8	TF (T10, T12)	Male	54	24.49
9	TOLF (T11–12)	Male	60	32.04	9	TDH (T11–L1)	Female	58	29.30
10	TOLF (T8–10)	Female	53	27.70	10	TDH (T8–10)	Male	61	25.26

Abbreviations: TOLF, thoracic ossification of ligamentum flavum; TDH, thoracic disc herniation; TF, thoracic fracture; BMI, body mass index.

Statistical Analysis and Visualization

All statistical analyses and visualization were performed using R software 3.6.3 and GraphPad Prism 7 (GraphPad Software Inc.). Differences of immune score, immunocyte fractions and gene expression between two groups were analyzed by Wilcoxon matched pairs test. Correlation analysis was conducted using the Pearson test. The results of qPCR were statistically analyzed by Student’s t-test with SPSS 22.0 (IBM Analytics, United States). All data were expressed as mean ± standard deviation. P < 0.05 was considered as statistically significant.

Results

Confirmation of Differential Immune and Macrophage Score Between OLF and Normal Samples

Based on data of the immune-related gene list and macrophage-related gene list, we first calculated and compared the immune score and macrophage score between OLF and normal tissues by ssGSEA algorithm. As a result, we found that the OLF samples were significantly associated with lower immune score (P=0.019) and macrophage score (P=0.029) than normal samples, which indicated that immune dysregulation might be implicated in the development of OLF (Figure 2).

Figure 2

The comparison of immune score and macrophage score between normal samples and OLF samples based on ssGSEA algorithm. (A) Immune score. (B) Macrophage score.

Distinction of Immune Cell Infiltration Between OLF and Normal Samples

The two-dimensional PCA depicted a significant difference in immune cell infiltration between OLF and control samples (Figure 3A and B). To comprehensively evaluate the landscape of immune cell infiltration between two groups, the percentage of each of the 23 types of immune cells and 64 types of immune cells in each sample was identified by ssGSEA and xCell, respectively (Figures 3C, D and 4A). Fourteen types of OLF-infiltrating immune cells (OIICs) were significantly altered, including natural killer T (NKT) cells, macrophages M2, NK CD56 bright cells, common myeloid progenitor (CMP) cells, basophils, conventional dendritic cells (cDCs), plasmacytoid dendritic cells (pDCs), pro B-cells, memory B cells, B cells, T helper 1 (Th1) cells, mesenchymal stem cells (MSC), smooth muscle cells and preadipocytes (Figure 4B). Correlation heatmap of these 14 types of immune cells revealed that pro B-cells had the strongest positive correlation with Basophils (r = 0.97), whereas NK CD56bright cells had the strongest negative correlation with cDC cells (r = −0.93) (Figure 4C). Among them, macrophages M2, NKT, NK CD56 bright cells, Th1 cells and smooth muscle cells were downregulated, whereas the other cells were upregulated in OLF samples (Figures 5 and 6).

Figure 3

The landscape of immune infiltration between OLF and normal tissues. (A) PCA cluster visualizing 64 types of differentially infiltrated immune cells based on xCell algorithm. (B) PCA visualizing 23 types of differentially infiltrated immune cells based on ssGSEA algorithm. (C) Heatmap of 23 infiltrating immune cells between normal and OLF tissues. (D) Heatmap of 64 infiltrating immune cells between normal and OLF tissues.

Figure 4

The situation of differential immune infiltration between normal and OLF samples. (A) The relative percentage of 14 differential immune cells between normal and OLF samples. Horizontal axes represents 4 × normal samples and 4 × OLF samples, and vertical axes represents the percentage of different immune cells in each sample. Color bars on right of the map represent immune cell types. (B) Radar chart depicting 14 types of infiltrating immune cells with different proportions. *, ** and *** indicate the significance level of 0.05, 0.01 and 0.001, respectively. (C) Correlation heatmap of 14 differential immune cells.

Figure 5

Up-regulated infiltrating immune cells in OLF. Red represents the expression of immune cells in normal LF tissues, and blue represents the expression of immune cells in OLF tissues. (A) Basophils. (B) B-cells. (C) cDC. (D) CMP. (E) Memory B-cells. (F) MSC. (G) pDC. (H) Pro-B cells.

Figure 6

Down-regulated infiltrating immune cells in OLF. Red represents the expression of immune cells in normal LF tissues, and blue represents the expression of immune cells in OLF tissues. (A) Macrophages M2. (B) NKT. (C) Preadipocytes. (D) Smooth muscle cells. (E) NK CD56bright cells. (F) Th1 cells.

Identification of Differentially Expressed Immune-Related Genes Between OLF and Normal Samples

A three-step computational framework was performed to explore hub OIDEGs and OIDELs. Firstly, 121 DEGs (51 downregulated and 70 upregulated) and 131 DELs (32 downregulated and 99 upregulated genes) in OLF were preliminarily screened out according to the predetermined screening cut-off criteria (Figure 7). Secondly, we constructed the intersection between 121 OLF-related DEGs and immune-related genes from ImmPort database to obtain 14 OIDEGs, including 7 downregulated genes (FAM3B, STAT3, SOCS3, IKBKE, LGR5, GDF15, TACR1) and 7 upregulated genes (CRH, IFNB1, CD40LG, CCL28, CD1E, GRP, MUC4) (Figure 8A–C). Table 2 showed the detailed information on these genes. CD1E (p=0.00034) and STAT3 (p=0.00039) were the most significantly up- and downregulated genes, respectively. Correlation analysis showed that STAT3 had the highest positive correlation with SOCS3 (r = 0.98, p = 0.01), but had the highest negative correlation with CD1E (r = −0.98, p = 0.03) (Figure 8D). Thirdly, Pearson correlation analysis inferred 46 featured lncRNAs that were highly correlated with these 14 OIDEGs (r > 0.9, P < 0.01), namely OIDELs, including 10 downregulated genes (eg, BMS1P10, HOPX, GAS2L1P2, SORD2P, ARHGAP22) and 36 upregulated genes (eg, ZNF195, SILC1, ZMAT4, DNAAF5, NUP210P1) (Figure 8E, Table 3). AC106779.1 (p=2.53E-07) and BMS1P10 (p=0.0005) were the most significantly up-regulated and down-regulated OIDELs, respectively.

Figure 7

Differentially expressed genes and lncRNAs between OLF samples and normal samples. (A) Volcano plot of differentially expressed genes. (B) Volcano plot of differentially expressed lncRNAs. (C) Heatmap of differentially expressed genes. (D) Heatmap of differentially expressed lncRNAs. Green dots indicates down-regulated, red indicates up-regulated and gray indicates no differential expression.

Figure 8

Immune-related differentially expressed genes and lncRNAs between OLF samples and normal samples. (A) Venn diagram showing the overlapping between immune-related genes and OLF-related differentially expressed genes. (B) Heatmap of 14 OLF-related and immune-related differentially expressed genes (OIDEGs). (C) Box plot showing expression level of 14 OIDEGs between normal and OLF samples. (D) Correlation heat map of 14 OIDEGs. (E) Heatmap of 46 OIDEGs-related differentially expressed lncRNAs.

Table 2

Detailed Information on 14 Hub OIDEGs

Gene Symbols	Full Names	Gene Function	P value	Regulation
CD1E	CD1e Molecule	It is a T-cell surface glycoprotein that soluble binds diacetylated lipids	3.46E-04	Up
STAT3	Signal Transducer And Activator Of Transcription 3	Mediates the expression of a variety of genes in response to cell stimuli, and thus plays a key role in many cellular processes such as cell growth and apoptosis	3.89E-04	Down
CRH	Corticotropin Releasing Hormone	A major regulator of homeostasis, mediating the autonomic, behavioral and neuroendocrine responses to stress and also play a role in cell growth and survival	6.20E-04	Up
FAM3B	FAM3 Metabolism Regulating Signaling Molecule B	Induces apoptosis of alpha and beta cells in a dose- and time-dependent manner	6.22E-04	Down
GRP	Gastrin Releasing Peptide	Regulates numerous functions of the gastrointestinal and central nervous systems, including release of gastrointestinal hormones, smooth muscle cell contraction, and epithelial cell proliferation	7.70E-04	Up
SOCS3	Suppressor of Cytokine Signaling 3	Cytokine-inducible negative regulators of cytokine signaling and the expression of this gene is induced by various cytokines, including IL6, IL10, and interferon (IFN)-gamma	1.80E-03	Down
IFNB1	Interferon Beta 1	The type I class of interferons involved in cell differentiation, antiviral, antibacterial, anticancer activities, inducing transcription of genes encoding inflammatory cytokines and chemokines.	2.37E-03	Up
CCL28	C-C Motif Chemokine Ligand 28	Chemotactic activity for resting CD4, CD8 T-cells and eosinophils and binds to CCR3 and CCR10 and induces calcium mobilization in a dose-dependent manner.	2.64E-03	Up
CD40LG	CD40 Ligand	Acts as a ligand for CD40 receptor for activation of CD40-CD40LG signaling, which have cell-type dependent effects, such as B-cell activation, NF-kappa-B signaling and anti-apoptotic signaling	3.86E-03	Up
GDF15	Growth Differentiation Factor 15	Regulates food intake, energy expenditure and body weight in response to metabolic and toxin-induced stresses	6.18E-03	Down
MUC4	Mucin 4	Alter cellular behavior through both anti-adhesive effects on cell-cell and cell-extracellular matrix interactions and in its ability to act as an intramembrane ligand for ERBB2	8.05E-03	Up
LGR5	Leucine Rich Repeat Containing G Protein-Coupled Receptor 5	Potentiates the canonical Wnt signaling pathway and acts as a stem cell marker of the intestinal epithelium and the hair follicle	1.06E-02	Down
IKBKE	Inhibitor Of Nuclear Factor Kappa B Kinase Subunit Epsilon	Plays an essential role in regulating inflammatory responses to viral infection, through the activation of the type I IFN, NF-kappa-B and STAT signaling	1.80E-02	Down
TACR1	Tachykinin Receptor 1	It is a receptor for the tachykinin neuropeptide substance P.	1.93E-02	Down

Abbreviation: OIDEGs, OLF-related and immune-related differentially expressed genes.

Table 3

OIDEGs-Related Differentially Expressed LncRNAs Between OLF and Normal Groups

OIDEG-Related DELs	Names	Number
Up-regulated	AC106779.1, NPRL3, CR383656.10, ZMAT4, NUP210P1, AC079089.1, FO393414.3, AC079943.2, MED21, LINC01137, RHOXF1P1, AC022816.1, SILC1, LINC01227, AC015983.2, PTPRM, AC097522.2, DNAAF5, MON2, DLX6-AS1, ADAMTS3, LINC02662, AC007922.1, AC016717.2, TRIM61, AC007953.2, SLC22A2, RFPL3S, LINC01700, AC188617.2, SLC22A16, POLR3A, AL365181.1, ZNF195, AL589987.2, AFAP1L2	36
Down-regulated	BMS1P10, HOPX, SORD2P, ARHGAP22, GAS2L1P2, LINC00601, AC096734.1, AL356495.1, LINC00682, LINC01558	10

Note: OIDEGs-related DELs were set as r > 0.9 and p value < 0.01.

Abbreviations: OIDEGs, OLF-related and immune-related differentially expressed genes; DELs, differentially expressed lncRNAs.

GO Function and KEGG Pathway Enrichment Analysis of OIDEGs

GO analysis showed that OIDEGs were mainly involved in receptor ligand activity, cytokine receptor binding, cytokine activity, G protein-coupled receptor binding, neuropeptide hormone activity, chemokine receptor binding, G protein-coupled peptide receptor activity, hormone activity, protein phosphatase binding, phosphatase binding, hormone receptor activity (Figure 9A and B). KEGG analysis demonstrated that OIDEGs were primarily involved in cytokine–cytokine receptor interaction, JAK-STAT signaling pathway, adipocytokine signaling pathway, RIG-I-like receptor signaling pathway, prolactin signaling pathway, Toll-like receptor signaling pathway, TNF signaling pathway, neuroactive ligand–receptor interaction, insulin resistance, growth hormone synthesis, secretion and action (Figure 9C and D). Among the pathways, “JAK-STAT signaling pathway”, “Toll-like receptor signaling pathway”, “TNF signaling pathway” were all well-known pathways involved in the process of osteogenesis and inflammation. In addition, many of the other pathways suggested metabolism-related involvement in OLF, such as “adipocytokine signaling pathway”, “prolactin signaling pathway” and “insulin resistance”, which suggested the interactions of immune-metabolism mechanism in the development of OLF.

Figure 9

Functional enrichment analysis of OIDEGs. (A) The top 20 GO terms with significantly differential expression from the GO enrichment analyses. (B) The most significant functions of hub enriched OIDEGs. (C) The top 20 KEGG pathways with significantly differential expression from KEGG enrichment analysis. (D) The most significant pathways of hub enriched OIDEGs.

Correlation Between OIDEGs and Infiltrating Immune Cells

According to r > 0.90 and p < 0.001, 21 OIDEG-OIICs correlation pairs were screened, and cDCs were significantly associated with the most OIDEGs (Figure 10A). For example, STAT3, SOCS3 and FAM3B were negatively correlated with cDCs (r = −0.99, p = 5.2E-06; r = −0.96, p = 0.00013; r = −0.95, p = 0.00022) (Figure 10B–D); CD1E and CD40LG were positively correlated with cDCs (r = 0.98, p = 3.3E-08; r = 0.96, p = 0.00022) (Figure 10E and F). In addition, STAT3 was positively correlated with NK CD56 bright cells (r = 0.94, p = 0.00043) (Figure 10G); LGR5 was positively correlated with preadipocytes (r = 0.96, p = 0.00014) (Figure 10H); IFNB1 was negatively correlated with smooth muscle cells (r = –0.94, p = 0.00041) (Figure 10I).

Figure 10

Correlation analysis between OIICs and OIDEGs. (A) Correlation coefficient matrix diagram between 14 OIICs and 14 OIDEGs. (B) STAT3-cDC correlation. (C) SOCS3-cDC correlation. (D) FAM3B-cDC correlation. (E) CD1E-cDC correlation. (F) CD40LG-cDC correlation. (G) STAT3-NK CD56 bright cells correlation. (H) LGR5-preadipocytes correlation. (I) IFNB1-smooth muscle cells correlation.

Prediction of LncRNA-Gene-Immunocyte Regulated Network in OLF

Based on the criteria of r > 0.95 or < −0.95 and p < 0.001, 73 OIDEG-OIDEL co-expression pairs were further selected, indicating that each lncRNA or mRNA could correlate with several lncRNAs or mRNAs (Figure 11A). CD1E and STAT3 were DEmRNAs with the highest degree of connectivity. In addition, SILC1 and ZMAT4 were positively correlated with CD1E, whereas BMS1P10 was negatively correlated with CD1E (Figure 11B–D). DLX6-AS1 was positively correlated with STAT3, while SILC1 was negatively correlated with STAT3 (Figure 11E and F). PTPRM was positively correlated with GRP (Figure 11G); DLX6-AS1 was positively correlated with SOCS3 (Figure 11H); DNAAF5 was negatively correlated with CRH (Figure 11I). Importantly, FO393414.3, AC096734.1, LINC01137 and DLX6-AS1 were OIDELs with the greatest number of OIDEGs, which were more likely to participate in this regulating process. Furthermore, we integrated these correlation results with the above correlation results of OIDEGs and OIICs to generate a specific lncRNA-gene-immunocyte regulatory network map (Figure 12). Targets included 10 OIDEGs, 41 OIDELs and 7 OIICs, which indicated that these regulators might participate in the control of OLF development and deserved further concern.

Figure 11

Correlation analysis between OIDELs and OIDEGs. (A) Correlation heat map between 46 OIDELs and 14 OIDEGs. (B) SILC1-CD1E correlation. (C) ZMAT4-CD1E correlation. (D) BMS1P10-CD1E correlation. (E) DLX6-AS1-STAT3 correlation. (F) SILC1-STAT3 correlation. (G) PTPRM-GRP correlation. (H) DLX6-AS1-SOCS3 correlation. (I) DNAAF5-CRH correlation.

Figure 12

Alluvial plot depicting predicted lncRNA-gene-immunocyte regulated network in OLF, including 10 OIDEGs, 41 OIDELs and 7 OIICs.

RT-PCR Validation of the Candidate OIDEGs Between OLF and Normal Samples

To confirm the accuracy of the prediction results, qRT-PCR was used to detect the expression of the above 10 hub OIDEGs, including CD1E, STAT3, CRH, FAM3B, GRP, SOCS3, IFNB1, CD40LG, MUC4 and LGR5. Eventually, changes in the expression levels of these 10 OIDEGs were paralleled in microarray and qRT-PCR experiments, thereby confirming the reliability and validity of our findings (Figure 13A–J). Among them, CRH and STAT3 were up-regulated and down-regulated genes with the highest fold change (Figure 13K).

Figure 13

Validation of the expression of 10 hub OIDEGs between OLF and normal controls by qRT-PCR. (A) CD1E. (B) STAT3. (C) CRH. (D) FAM3B. (E) GRP. (F) SOCS3. (G) IFNB1. (H) CD40LG. (I) MUC4. (J) LGR5. (K) Fold changes of the mean expression value of OIDEGs. Error bars represented mean ± standard deviation for each group. **P < 0.01, ***P < 0.001.

Discussion

OLF is characterized by heterotopic bone formation of ligamentum flavum tissues through endochondral ossification. In spite of a multitude of investigations, the exact pathogenesis of OLF has not been fully elucidated. Considering potential links between the immune system and abnormal bone formation have been frequently reported, we speculated that immune dysregulation might be implicated in osteogenic differentiation of ligament fibroblasts during OLF. More importantly, it was firstly established that immune score and macrophage score were significantly differential between OLF and normal samples. This aroused our interest in investigating whether immune-related functions and mechanisms in the development of OLF. To the best of our knowledge, this was the first tentative exploration for key factors that have potentially pathogenic effects, such as immune-related genes, molecular pathways and infiltrating immune cells, and their biological and functional correlations in OLF based on integrating bioinformatic analysis with real-world data, which provided novel insights into the pathogenesis and therapeutic targets of OLF.

Different expressional cell patterns of immune infiltration were identified to reveal that dysregulation of initial innate immunity and subsequent adaptive immunity might be associated with OLF. For example, infiltrating M2 macrophages, which are important component of the innate immunity, were lower in OLF samples than normal samples. Considering that M1 macrophages were generally pro-inflammatory, while M2 macrophages were anti-inflammatory,27 lower proportions of M2 macrophages in OLF tissues might contribute to more severe inflammatory and osteogenic processes. On the contrary, recovery of M2 macrophages might antagonize M1 macrophage responses to inhibit inflammation. Additionally, higher proportions of pro B-cells were detected in OLF tissues, representing that dysregulated cell population of adaptive immunity in OLF might be involved in disease development. Admittedly, mature B-cells could inhibit osteoblastic differentiation, thus we speculated that alterations in molecules or transcription factors during OLF precluded normal B-cell differentiation and resulted in sustaining osteogenic differentiation in unpredicted ways.28 Interestingly, cDC cells, which connected the innate and adaptive arms of the immune system, emerged as having the most prominent differential expression. Moreover, cDCs had an extraordinary capacity to produce bioactive molecules that induced and sustained inflammation as well as neovascular formation, which might indicate a crucial regulatory role in new bone formation.29 However, its concrete effects and mechanisms in the initiation and progression of OLF need to be studied by further research.

Through a series of bioinformatics algorithms and experimental verification, CD1E, CD40LG, MUC4, IFNB1, CRH, GRP, FAM3B, STAT3, SOCS3 and LGR5 were identified as hub immune-related genes implicated in OLF. Among them, STAT3 acts as a well-known OLF-associated gene, which plays a significant role in many cellular processes related to human osteogenesis, such as angiogenesis, inflammation, and autophagy.30–32 In addition, STAT3-related immunoregulatory mechanisms during new bone formation were also reported. Nicolaidou et al33 found monocytes/macrophages could potently induce the osteogenic differentiation of human MSCs via inducing the local activation of STAT3 in ankylosing spondylitis. SOCS3, encoding a member of the STAT-induced STAT inhibitor, has been predicted to be enriched in “adaptive immune response” and “chemokine signaling pathway” processes in OLF.17 Furthermore, a recent study revealed the SOCS3-STAT3 pathway, activated by the CUE domain-containing 2 (CUEDC2), accelerated osteoblast differentiation and bone formation.34 IFNB1, also known as Interferon Beta 1, belongs to the Type I interferons, which can induce transcription of genes encoding inflammatory cytokines and chemokines. Emerging studies highlighted the importance of the interaction between IFNB and bone metabolisms.35–37 For instance, osteoblasts treated with IFNβ exhibited a decreased mineralization, showing that IFNβ could be involved in low bone mass or osteoporosis.37 CD40LG, a gene regulating B cell function by engaging CD40, has been considered to regulate bone metabolism. Under appropriate conditions, CD40L suppression could also promote Wnt-10b production, bone formation and accretion of bone mass in the axial skeleton.38 LGR5 is the receptor for R-spondins that potentiates the canonical Wnt signaling pathway. Recently, LGR5 was found to be implicated in the cellular processes of osteogenic differentiation of MSCs through Wnt signaling pathway in fracture healing.39 Importantly, Wnt signaling pathway was highly correlated with the development of OLF; thus, the promising role of LGR5 during OLF should be further elucidated.40 GRP, a member of bombesin-like peptides, has been involved in several physiological and pathological processes. A latest study revealed the role of GRP inhibition on the amelioration of phosphate-induced vascular calcification by inhibiting osteogenic differentiation of vascular smooth muscle cells (VSMCs) in vitro and in vivo.41 As we know, the remaining four genes were the first time to be reported to hold potential relationship with osteogenic processes, which provided new directions in related fields.

Based on GO and KEGG analyses, OIDEGs-related biological processes and pathways were mainly related to inflammation and immunity. KEGG analysis demonstrated that JAK-STAT, Toll-like receptor and TNF signaling pathways might be involved in the development of OLF. A previous study showed that the STAT3 signaling pathway was activated by leptin-induced osteogenic differentiation of OLF cells.42 In addition, the activation of the JAK2-STAT3 pathway also had osteogenic potentials on OPLL cells.43 An iTRAQ proteomics study revealed that the TNF signaling pathway promoted the osteogenic differentiation of LF cells by activating Osx expression.12 Besides, numerous researches demonstrated that the TNF signaling enhanced the osteogenic differentiation ability of periodontal ligament stem cells (PDLSCs) in periodontitis.44–46 However, no studies were conducted to uncover the role of the toll-like receptor signaling pathway in OLF. He et al inferred that toll-like receptor signaling pathway was involved in the development of OPLL by mining microarray gene profiling.47 Moreover, the regulation of toll-like receptor signaling in human PDLSCs under inflammatory conditions has been fully elucidated.48,49 Our results revealed the significant involvement of the toll-like receptor signaling pathway in OLF, and further studies need to explore its underlying mechanisms in mediating the progression of OLF.

Establishment of the potential interactions among OIDEGs, OIDELs and OIICs could provide momentous clues for probing the role of immunoregulatory mechanisms in OLF. LncRNAs can affect the gene expression in various manners, including serving as competing endogenous RNAs (ceRNA), regulating DNA methylation of genes, modulating the promoter of genes.50 AC106779.1 and BMS1P10 were the most significantly up- and down-regulated OIDELs in OLF tissues in the present study, respectively; however, their biological functions remained so obscure that further studies were necessary to validate their potential as diagnostic biomarkers for OLF. Furthermore, FO393414.3, AC096734.1, LINC01137 and DLX6-AS1 were predicted as pivotal lncRNAs of the co-expressed OIDELs-OIDEGs pairs, which modulated the majority of the OIDEGs in OLF. It was hypothesized that these OIDELs might serve crucial roles by regulating the expression of STAT3, CD1E, SOCS3, and their associated signaling pathways. Regarding the immunocyte–gene interaction, not only did cDCs serve as most prominently up-regulated immune cell in OLF tissues, but it also held the highest degree of connectivity with hub OIDEGs, including up-regulated CD40LG and CD1E and down-regulated STAT3 and FAM3B. In addition, the most significantly down-regulated NK CD56bright cells were another valuable object which also concatenated many hub genes. These crucial immune-related regulators in this network deserved special mention.

Several inevitable limitations of this study need to be taken into account when interpreting the findings. First, only a few RNA-sequencing datasets of OLF were available on the GEO database. Moreover, we had no access to clinical data such as disease phenotypes and radiological data from the published study and public database. The small number of samples and incomplete clinical data made subgroup analyses impossible, which may lead to biased interpretation of the results in the case of group heterogeneity and the influence of other factors. Secondly, all of these analyses were obtained by data mining based on a series of bioinformatic algorithms, but these algorithms were commonly used and scientific, and we also preliminarily performed rigorous experimental demonstration to validate the feasibility of these algorithms to some extent. Thirdly, the exact underlying molecular mechanisms of candidate immune-related genes and immune cells in OLF need to be further investigated using in-depth in vitro and in vivo scientific experiments. Finally, the nature of retrospective research limited the clinical value of this work; thus, multicentre or prospective studies are imperative to elucidate the relationship between immune response and the development of OLF.

Conclusion

Taken together, through analyzing an OLF microarray, we found that OLF tissues and normal controls could be clearly distinguished via the overall immune score and featured infiltrating immune cells, and the hub differential immune response-related genes (eg, FAM3B, STAT3, SOCS3, CD40LG, CD1E) were identified. By combining a series of reliable bioinformatic algorithms with real-world data verification, we considered that these genes might participate in the immune response and play important roles in the development of OLF. Furthermore, highly correlated immunocyte-gene-lncRNA network was also preliminarily constructed. These findings may not only provide compelling insights into the underlying immunoregulatory machinery in OLF but also shed lights on promising targets for treatment strategies.