Transcriptome Sequencing Identifies Abnormal lncRNAs and mRNAs and Reveals Potentially Hub Immune-Related mRNA in Osteoporosis with Vertebral Fracture

Introduction

Osteoporosis is a systemic bone disease characterized by reduced bone mass and degradation of bone microstructure, which leads to increased brittleness of bone and is prone to fracture.^1,2 Osteoporosis could occur at any age, but it is more common in the elderly,³ and postmenopausal women.⁴ With the increase in world population and life expectancy and the arrival of aging, the threat of this disease will become greater and greater. The current situation of osteoporosis prevention and treatment in China is not optimistic. In a recent national prevalence analysis of osteoporosis,⁵ the prevalence of osteoporosis is 5.0% in males aged 40 years or older and was 20.6% among women. In addition, the prevalence of osteoporosis increased with age. A meta-analysis indicated the prevalence of osteoporosis in the elderly aged ≥ 60 years in China was up to 37.7%.⁶ The formation and absorption of bone is a process of dynamic balance. In childhood and adolescence, bones grow, strength, and mineral content, and bone formation exceeds bone absorption.⁷ However, the osteoblast function gradually decreases after reaching its peak, the reduced mineral and organic matrix of bone cause the decline of bone quality, which is the main cause of osteoporosis.^8,9 The etiology of osteoporosis mainly includes endocrine disorder,¹⁰ metabolic disorder,¹¹ and mechanical factors.¹²

Recent studies have put forward the viewpoint of “bone immunology”, which holds that the immune system and immune factors play an important regulatory role in the occurrence and development of osteoporosis.^13,14 Bone immunology specializes in the interaction between the immune system, hematopoietic system, and bone.¹⁵ Take T cells in the immune system as an example, the T lymphoid cells secrete inflammatory factors and Wnt ligands, which promote bone formation and absorption. In addition, T lymphocytes regulate the dynamic balance of metabolism between bone stromal cells and osteoblasts through CD40 ligands and costimulatory molecules.¹³

Fracture is a serious complication of osteoporosis, which affects people’s normal life. Genomics research has continuously confirmed that diseases are caused by external risk factors and internal gene changes, including skeletal system diseases.¹⁶ In this study, transcriptome sequencing analysis was performed in blood samples from healthy control (HC), patients with osteoporosis without vertebral fracture (OWF), and patients with osteoporosis with vertebral fracture (OVF). Subsequently, differential expression mRNAs (DEmRNAs) and differential expression lncRNAs (DElncRNAs) were screened out. The focus is on screening fracture-related immune mRNAs and lncRNA and provides a useful reference for identifying patients with high-risk fractures.

Methods

Patients and Samples

The peripheral blood samples were obtained from 3 groups of female subjects with age over 55 years: healthy control (HC), osteoporosis patients without vertebral fracture (OWF), and osteoporosis patients combined with vertebral fracture (OVF). The ethnicity of all individuals is Han. The HC, OWF and OVF groups included 4, 6, and 6 subjects, respectively. In this study, dual-energy X-ray absorptiometry (DXA) was used to evaluate the bone mineral density (BMD)¹⁷ in recruited subjects. BMD analysis was based on femur. World Health Organization (WHO) criteria for osteoporosis using BMD.¹⁸ Classification according to T-Score value: ① Normal: greater than −1.0; ② Osteopenia: −1.0 to −2.5 ③ Osteoporosis: less than −2.5. Clinical manifestations such as pain, swelling, and abnormal activity, and combined general manifestations of osteoporosis. The diagnostic criteria of patients with OVF include: the decrease of BMD meets the diagnostic criteria of osteoporosis and may present clinical manifestations such as chest and waist pain, shortened height, and limited activity; the fracture was scanned and determined by X-ray, computed tomography (CT), magnetic resonance imaging (MRI) or whole-body bone-scan by emission CT (ECT). The vertebral fracture time range in the OVF group was from November 1, 2022 to December 30, 2022. Inclusion criteria for OWF: except for no fracture, other diagnostic criteria meet the diagnostic criteria for osteoporosis. Subjects with metastatic bone tumor, thoracolumbar tuberculosis, multiple myeloma, endocrine diseases such as hyperparathyroidism, immune diseases such as rheumatoid arthritis, taking drugs affecting bone metabolism, or various congenital or acquired abnormal bone metabolism diseases were excluded. The individuals in the HC group were gender and age matched with the case group and had no disease before sampling.

RNA Sequencing Library Preparation

The total RNA in blood was extracted by Trizol reagent. The DNA fragments in samples were digested by DNase I, purified and recovered by magnetic beads. The rRNA was removed with the commercial kit and the purification of RNA samples was verified by Agilent 2100 Bioanalyzer. After the thermal interruption, a strand of cDNA was synthesized according to the corresponding procedure. Then, through the synthesis of reverse transcription double-strand and terminal repair, the cDNA was double linked. The double-stranded DNA was digested by the uracil-N-glycosylase (UNG) digestion reaction system and the product was purified and recovered with magnetic beads, and then amplified by polymerase chain reaction (PCR). The quality of the RNA library was measured by Agilent 2100 Bioanalyzer and ABI StepOnePlus real-time PCR system after purification.

RNA Sequencing

The qualified library was denatured into a single chain by NaOH and diluted to optimal concentration on the FlowCell. The RNA fragment was hybridized with the connector on FlowCell. Hybridized molecules were amplified by bridge PCR on cBot and finally sequenced using the BGIseq platform with library strategy: PE100.

Raw Data Processing

The bases with mass<20 and N greater than 10% of the 5′ and 3′ segments of reads were trimmed using fastp. The obtained clean reads were aligned to the human reference genome (GRCh38) in the Ensemble database (gene annotation: Ensemble 92) using HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts, https://ccb.jhu.edu/software/hisat2/index.shtml) v 2.1.0. The expressions level of mRNA and lncRNA were quantified and output using the “standard” approach by Stringtie v1.3.3b (http://ccb.jhu.edu/software/stringtie/).

Differential Expression Analysis

Difference analysis mainly included three steps: we firstly normalized the original read count mainly by correcting the sequencing depth. Then, the hypothesis test probability (P-value, P) was calculated through the statistical model. The corrected P value (false discovery rate, FDR) was obtained by conducting a multiple hypothesis tests (Benjamini-Hochberg) correction. The Deseq2 was used to compare the differential expression differences of lncRNA and mRNA between different groups. Parameters with P < 0.01 and | log2 fold change | (| log2 FC |) > 1 were screened.

Identification and Co-Expression Analysis of Fracture-Related DEmRNAs and DElncRNAs

The common DEmRNAs and DElncRNAs of the OWF vs HC and OVF vs OWF groups were obtained by venn package, but they were not differentially expressed in OWF vs HC group. These common DEmRNAs and DElncRNAs were considered as fracture-related DEmRNAs and DElncRNAs. The Pearson correlation coefficient method was used to analyze the correlation between lncRNA and distant mRNA among samples. The screening criteria were |correlation coefficient| (|r|) ≥ 0.8 and P <0.001. Subsequently, lncRNA-mRNA pairs associated with fracture were screened based on |r| ≥ 0.8 and P <0.001 and co-expression network were mapped.

Functional Enrichment Analysis

The Gene Ontology (GO) consists of three main categories, namely Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). The genecodis 4.0 database (https://genecodis.genyo.es/) was used for Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) functional enrichment analysis of DEmRNAs in co-expression network. The screening criteria was pval_adj <0.05.

Identification of Immune-Related mRNAs in Fracture-Related DEmRNAs

The immune-related mRNAs were downloaded from the ImmPort database (https://www.immport.org/shared/home). Then, the immune-related mRNAs were intersected with the fracture-related DEmRNAs to obtain the fracture-related immune DEmRNAs. Then, the protein-protein interaction (PPI) network was constructed by STRING (https://cn.string-db.org/) to explore the interaction between fracture-related immune DEmRNAs.

Identification of Hub DEmRNAs

To screen hub DEmRNAs, DEmRNAs obtained from STRING were imported into the Cytoscape software (http://www.cytoscape.org/) and the Cytohubba plug-in was used to screen the hub DEmRNAs. Four algorithms, including density of maximum neighborhood component (DMNC), Maximal Clique Centrality (MCC), Degree, and Edge Percolated Component (EPC), were adopted. The top 10 DEmRNAs of each algorithm were taken to intersect, and the intersection were regarded as the hub DEmRNAs.

Construction of Hub DEmRNAs Related lncRNA-miRNA-mRNA Network and Transcription Factors (TFs) Network

The TarBase database of the NetworkAnalyst platform (https://www.networkanalyst.ca/) was used to predict miRNAs associated with fracture-related DEmRNAs interactions to obtain mRNA-miRNA pairs. The LncBase v.2 database (https://dianalab.e-ce.uth.gr/html/diana/web/index.php?r=lncbasev2) was used to predict miRNAs associated with fracture-related DElncRNAs interactions to obtain lncRNA-miRNA pairs (score >0.7). Then, the mRNA-miRNA pairs and lncRNA-miRNA pairs were fused to construct the lncRNA-miRNA-mRNA network. Subsequently, the lncRNA-miRNA-mRNA sub-network related to hub DEmRNAs was isolated for visualization. In addition, the JASPAR database of the NetworkAnalyst platform was used to predict TFs related to hub DEmRNAs and construct a TFs network.

Validation of Real-Time PCR

The real-time PCR analysis sample was inconsistent with the preparation of the library sample. The individual samples included in real-time PCR analysis were consistent with the sample screening criteria included in sequencing. The vertebral fracture time range in the OVF group was from March 8, 2023 to April 7, 2023. A total of 32 (including 11 HC, 9 OWF and 6 OVF) peripheral blood samples were selected for real-time PCR verification. Total RNA was extracted using RNAliquid hyperspeed whole blood (liquid sample) total RNA extraction kit (Beijing Hutian Oriental Technology Co., LTD.). The kit used for RNA reverse transcription into DNA was the FastKing cDNA first strand synthesis kit (TIANGEN). The reaction system and amplification procedure of real-time PCR were carried out according to the instruction manual of SuperReal PreMix Plus (SYBR Green). GAPDH and ACTB were internal reference genes. The 2^{−ΔΔ CT} method was used for relative quantitative analysis of data.

Results

Identification of DEmRNAs

DEmRNAs were screened with criteria of P< 0.01 and |log2FC|> 1. Compared with the HC group, 95 DEmRNAs (32 up-regulated and 63 down-regulated) were obtained in the OWF group. The Volcano and heat maps of the DEmRNAs are shown in Figure 1A and B, respectively. Compared with the HC group, 464 DEmRNAs (163 up-regulated and 301 down-regulated) were obtained in the OVF group. The Volcano and heat maps of the DEmRNAs are shown in Figure 1C and D, respectively. Compared with the OWF group, 585 DEmRNAs (147 up-regulated and 438 down-regulated) were obtained in the OVF group. The Volcano and heat maps of the DEmRNAs are shown in Figure 1E and F, respectively. The top 20 DEmRNAs in OWF vs HC, OVF vs HC and OVF vs OWF groups were listed in Tables S1–S3.

Figure 1 Volcano maps (left) and heat maps (right) of DEmRNAs. (A and B): Volcano map (left) and heat map (right) of DEmRNAs in the OWF vs HC group; (C and D): Volcano map (left) and heat map (right) of DEmRNAs in the OVF vs HC group; (E and F), Volcano map (left) and heat map (right) of DEmRNAs in the OVF vs OWF group. HC, OWF and OVF represent the healthy control, osteoporosis patients without vertebral fracture and osteoporosis patients combined with vertebral fracture, respectively.

Identification of DElncRNAs

In this study, differential expression analysis was performed in OWF vs HC, OVF vs OWF, and OVF vs HC groups, respectively. The screening criteria were | log2FC | > 1 and P < 0.05. Compared with the HC group, 26 DElncRNAs (9 up-regulated and 17 down-regulated) were obtained in the OWF group. The Volcano and heat maps of the DElncRNAs are shown in Figure 2A and B, respectively. Compared with the HC group, 162 DElncRNAs (51 up-regulated and 111 down-regulated) were obtained in the OVF group. The Volcano and heat maps of the DElncRNAs are shown in Figure 2C and D, respectively. Compared with the OWF group, 170 DElncRNAs (45 up-regulated and 125 down-regulated) were obtained in the OVF group. The Volcano and heat maps of the DElncRNAs are shown in Figure 2E and F, respectively. The top 20 DElncRNAs in OWF vs HC, OVF vs HC and OVF vs OWF were listed in Tables S4–S6.

Figure 2 The volcano map (left) and heat map (right) of DElncRNAs. (A and B): Volcano map (left) and heat map (right) of DElncRNAs in the OWF vs HC group; (C and D): Volcano map (left) and heat map (right) of DElncRNAs in the OVF vs HC group; (E and F): Volcano map (left) and heat map (right) of DElncRNAs in the OVF vs OWF group. HC, OWF and OVF represent the healthy control, osteoporosis patients without vertebral fracture and osteoporosis patients combined with vertebral fracture, respectively.

Analysis of Fracture-Related DEmRNAs and DElncRNAs

Venn diagram results showed that there were 13 common DEmRNAs were found in OWF vs HC and OVF vs HC groups, and 33 common DEmRNAs were found in OWF vs HC and OVF vs OWF groups; there were 217 common DEmRNAs were found in OVF vs HC and OVF vs OWF groups (Figure 3A). The 217 DEmRNAs shared by the OVF vs HC and OVF vs OWF groups were considered fracture-related DEmRNAs. In addition, 58 DElncRNAs (Figure 3B) among OVF vs HC and OVF vs OWF while having no differential expression between OWF vs HC were selected and classified as fracture-related DElncRNAs.

Figure 3 The Venn diagram of DEmRNAs and DElncRNAs in OWF vs HC, OVF vs HC and OVF vs OWF groups. (A): Venn diagram of DEmRNAs; (B), Venn diagram of DElncRNAs. HC, OWF and OVF represent the healthy control, osteoporosis patients without vertebral fracture and osteoporosis patients combined with vertebral fracture, respectively.

Fracture-Related lncRNAs and mRNAs Co-Expression

The target genes were predicted by analyzing the correlation between the expression of lncRNA and protein-coding genes among samples. A total of 3378 lncRNA-mRNA pairs were obtained, and their co-expression network diagram was listed in Figure 4. GO and KEGG functional enrichment analysis was performed on the GeneCodis4.0 database using 193 co-expressed DEmRNAs with DElncRNAs. Screening criteria: pval_adj < 0.05, the enrichment results of top 15 GO terms and KEGG terms are shown in Figure 5A–D. These results demonstrate that lncRNA co-expressed mRNAs were mainly enriched in immune-related pathways, especially in GO-BP analysis (Figure 5A).

Figure 4 The co-expression network diagram of obtained mRNA and lncRNA. Circle, V-shape, red and blue represents DEmRNA, DElncRNA, up-regulation and down-regulation, respectively.

Figure 5 Bubble plot of GO and KEGG pathway analysis of co-expressed DEmRNAs with DElncRNAs. (A): Top 15 significantly enriched biological processes (BP); (B): Top 15 significantly enriched cell component (CC); (C): Top 15 significantly enriched molecular function (MF); (D): Top 15 significantly enriched KEGG pathway.

Immune-Related Co-Expressed Fracture-Related DElncRNAs and DEmRNAs

A total of immune-related 1793 mRNAs were obtained from the ImmPort database. After the intersection with fracture-related DEmRNAs, we obtained 33 immune-related DEmRNAs, 25 of which were co-expressed with DElncRNAs (Figure 6A). The co-expression network was listed in Figure 6B. The PPI network was established (Figure 6C) to further explore the interaction between 33 immune-related DEmRNAs in the fracture. A total of 8 hub DEmRNAs were obtained after the intersection of the top 10 DEmRNAs of each algorithm (Figure 6D), including IL18R1, IL18RAP, SLC11A1, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1.

Figure 6 Immune-related co-expressed fracture-related lncRNAs and mRNAs. (A): Venn diagram of 33 immune-related DEmRNAs; (B): The co-expression network of 25 immune-related DEmRNAs; (C): the PPI network of 33 immune-related DEmRNAs; (D): 8 core mRNAs (IL18R1, IL18RAP, SLC11A1, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1) screened by 4 algorithms.

Regulatory Networks Associated with Hub DEmRNAs

The lncRNA-miRNA-mRNA sub-network consists of 50 nodes and 87 edges (Figure 7A). The 50 nodes included 32 miRNAs, 7 hub mRNAs (IL18R1, IL18RAP, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1) and 11 lncRNAs. The TFs network consisted of 58 mRNA-TF pairs, including 8 hub mRNAs (IL18R1, IL18RAP, SLC11A1, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1) and 37 TFs (Figure 7B). In addition, for the 8 hub DEmRNAs obtained above, we also constructed the hub mRNA-lncRNA co-expression sub-network. The results showed that 7 hub DEmRNAs (IL18R1, IL18RAP, SLC11A1, CSF2RA, IL-1R2, PGLYRP1, and IL-1R1) have lncRNA-mRNA co-expression relationship (Figure 8).

Figure 7 The lncRNA-miRNA-mRNA network and TFs network associated with hub DEmRNAs. (A): The lncRNA-miRNA-mRNA network network associated with hub DEmRNAs. Circles, squares and V-shapes represent DEmRNA, miRNA and DElncRNA, respectively. (B): The TFs network associated with hub DEmRNAs. Circles and squares represent DEmRNA and TFs, respectively.

Figure 8 The hub mRNA-lncRNA co-expression sub-network. Seven genes (IL18R1, IL18RAP, SLC11A1, CSF2RA, IL1R2, PGLYRP1, and IL1R1) were included. Circle, V-shape, red and blue represents DEmRNA, DElncRNA, up-regulation and down-regulation, respectively.

Expression Validation of AC138035.2, AL844892.2, CSF2RA and IL1R2 by Real-Time PCR

AC138035.2, AL844892.2, CSF2RA and IL1R2 were randomly selected for real-time PCR validation. The primers used for real-time PCR were shown in Table 1. Compared with the HC, the expression of AC138035.2, CSF2RA and IL1R2 in OVF group were down-regulated trend, while the expression of AL844892.2 was up-regulated (Figure 9A). Similarly, compared with the OWF, the expression of AC138035.2, CSF2RA and IL1R2 in OVF group were also down-regulated trend, while the expression of AL844892.2 was also up-regulated trend (Figure 9B). The expression trend of real-time PCR results was consistent with the previous sequencing results. However, although the expression trend was consistent, most genes lacked significance, which may be caused by sample heterogeneity, RNA degradation caused by improper operation process and small sample size. Therefore, further research will be conducted by expanding the sample size in the future.

Table 1 All Primer Sequences Used for Real-Time PCR Validation

Figure 9 Expression validation of AC138035.2, AL844892.2, CSF2RA and IL1R2 by real-time PCR. (A): Relative expression levels of AC138035.2, AL844892.2, CSF2RA and IL1R2 in OVF vs HC group; (B): Relative expression levels of AC138035.2, AL844892.2, CSF2RA and IL1R2 in OVF vs OWF group. *P <0.05. HC, OWF and OVF represent the healthy control, osteoporosis patients without vertebral fracture and osteoporosis patients combined with vertebral fracture, respectively.

Discussions

Osteoporosis is a systemic and metabolic bone system disease. Its pathological characteristics are decreased bone mass, destruction of bone microstructure, increased bone fragility, decreased bone strength, and prone to fracture.¹⁹ Bone remodeling is regulated by the balance between osteoblasts/osteoclasts and the interaction between bone and the immune system. Osteoblasts originate from bone marrow mesenchymal stem cells. Their formation and function are regulated by bone morphogenetic protein, transforming growth factor-β, Wnt, and other molecules.²⁰ The key factors of osteoblast differentiation and proliferation are Runt-related transcription factor 2, bone matrix protein, and other molecules. Moreover, osteoblasts produce positive and negative regulators of osteoclastogenesis, which are receptor activator NFkB ligand (RANKL) and osteoprotegerin (OPG), respectively.²¹ Osteoclasts originate from hematopoietic stem cells and their differentiation and activation are induced by RANKL²² and macrophage colony-stimulating factor.²³ Immune cells are involved in the regulation of osteoblasts and osteoclasts to varying degrees. T cells are differentiated from lymphoid stem cells in the thymus and are the main sources of RANKL and TNF- α.²⁴ After the inflammatory factors are combined with T cells, they release factors such as IL-1 and IL-17, which have a significant impact on bone metabolism.^25–27 The B cells in immune cells are formed by the development and differentiation of hematopoietic stem cells (HSCs) in the bone marrow. Its development depends on the expression of RANKL, CXCL12, and IL-7 in HSCs.²⁸

Osteoporosis increases fracture risk, and high inflammatory status is an important influencing factor of osteoporosis and fracture.²⁹ In this study, we measured differential expression profiles of mRNA and lncRNA in OVF vs HC, OWF vs HC, and OVF vs OWF. The GO-BP analysis results indicated innate immune response, inflammatory response, and immune response are the main functioned process of mRNAs co-expressed with lncRNAs. Subsequently, 33 immune-related DEmRNAs that may be involved in fracture were obtained by intersecting immune-related mRNAs with fracture-related DEmRNAs. CytoHubba plug-in was used to screen hub immune-related DEmRNAs based on 33 immune-related DEmRNAs. Finally, 8 hub immune-related DEmRNAs were identified (IL18R1, IL18RAP, SLC11A1, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1). These mRNAs are involved in the process of bone immunomodulation.

The interleukin-1 receptor (IL1R) is a member of the immunoglobulin superfamily. IL1R1 and IL1R2 are two subtypes of IL1R, and they bind different chains of IL-1. IL1R1 is the predominant form found on T cells and fibroblasts, whereas IL1R2 is the predominant form found on B cells, monocytes, neutrophils, and bone marrow cells.³⁰ In pathologically activated osteoclasts, preferentially expressed of IL1R1 and inhibition expressed of IL1R2 can cause severe bone destruction.³¹ IL1R2 acts as a decoy receptor that inhibits IL-1-mediated inflammatory responses.³² While stimulation of IL1R1 activates the proinflammatory signaling cascade that induces TNF induction.³³ Our GO enrichment analysis, especially the biological process also demonstrated that IL1R1 participates in positive regulation of interferon-gamma production, while IL1R2 is negative regulating interleukin-1 alpha production and participates in the pathogenesis of osteoporosis.³⁴ Research has indicated that interferon-gamma promote osteoblast differentiation.³⁵ Further, IL1R1 knockout mice have significantly reduced osteoclast numbers in the chondro-osseous junction, trabecular bone, and cortical bone.³⁶ In the co-expression network, we found that IL1R1 and IL1R2 are down-regulated. This indicated the anti-inflammation effect of IL1R2 is inhibited. However, no relevant explanation has been found for the down-regulation of IL1R1. This may need further exploration in follow-up research.

Interleukin 18 (IL-18) is a member of the IL-1 family. It activates NFkB and induces inflammatory mediators via the same signaling pathway as IL-1.³⁷ Interleukin 18 receptor 1 (IL18R1), also known as IL18Ra, is a decoy receptor or binding protein (IL18BP) of IL18. The report indicated that IL18R1 is decreased in osteoporotic women,³⁸ which is consistent with our study. Interleukin 18 receptor accessory protein (IL18RAP) combined with IL18R1 to dampen IL-18 activities as a negative feedback loop.³⁹ Down-regulation of IL18RAP means that the negative regulation of IL-18 mediated immune and inflammatory response in patients with an osteoporotic fracture is suppressed. Our GO analysis revealed that the most enrichment pathways of IL18R1 and IL18RAP were coincident and mainly involved in inflammatory and immune-related biological processes. The 2 genes regulate the pathogenesis of osteoporosis through IL-1 signaling pathway.⁴⁰ This result provides a research direction for further research on osteoporosis.

Colony Stimulating Factor 2 Receptor Subunit Alpha (CSF2RA) is a subunit receptor of Granulocyte macrophage-colony Stimulating Factor (GM-CSF).⁴¹ In a transcription analysis, Anam and Davis found that CSF2RA is specifically up-regulated in hematopoietic stem/progenitor cells (HSPCs) relative to bone marrow mesenchymal stromal cells (BMMSCs).⁴² The overexpression of CSF2RA has been reported in multiple tumors.⁴³ In a profile analysis of osteoblast, Kalajzic et al⁴⁴ reported CSF2RA has higher expression levels in unsorted differentiating cultures of bone cell populations while showing significantly lower expression in sorted mature osteoblasts. This indicates that CSF2RA could promote osteoblast differentiation, and its expression is impaired in the process of osteoporosis. The same as our GO enrichment results of biological process and molecular function, that CSF2RA participates in cytokine-mediated signaling pathway and regulating protein binding and cytokine receptor activity, which is related to the osteoblast function.⁴⁵

Peptidoglycan recognition protein 1 (PGLYRP1) plays a role in innate immunity and has bacteriostatic activity. Yang et al⁴⁶ identified it work as a hub mRNA of the osteoarthritis subchondral bone. Our research shows that it is also a hub mRNA regulated by lncRNA in osteoporosis. The biological process of GO enrichment of our study indicated PGLYRP1 is involved in innate immune response and negative regulation of interferon-gamma production.³⁵ Decreased expression of solute carrier family 11 member 1 (SLC11A1) has been reported to differentiate osteoclasts.⁴⁷ We noticed SLC11A1 is positive regulation of interferon-gamma production³⁵ and involved in cellular cadmium ion homeostasis⁴⁸ in our GO analysis. These results suggest that PGLYRP1 and SLC11A1 may be involved in the progression of osteoporosis by regulating multiple biological processes. However, its specific regulatory process in osteoporotic fractures remains to be revealed. Dysregulation of calcium plays an important role in osteoporosis.⁴⁸ Through binding a variety of chemokines, C-C type chemokine (CCR3) subsequently transduces a signal by increasing the intracellular calcium ions level. Among the hub DEmRNA, CCR3 is the only mRNA upregulated, which is consistent with a previous report,⁴⁹ that CCR3 expression is highly upregulated during osteoclastogenesis. Guerard et al⁵⁰ indicated that CCR3 promotes homing of prostate cancer cells to bone, which suggests that CCR3 might also involve the process of osteoclast migration. Our GO and KEGG analysis demonstrated that CCR3 was not enriched in any pathway yet, follow up research is needed to further clarify its regulation pathway in osteoporosis.

However, this experiment has certain limitations. The specific molecular mechanism of the identified mRNAs and lncRNAs in osteoporosis is still unclear. This means that the regulatory function of these molecules needs to be further clarified in the follow-up study. In summary, we found 8 hub immune-related mRNAs (including IL18R1, IL18RAP, SLC11A1, CSF2RA, CCR3, IL1R2, PGLYRP1, and IL1R1). These mRNAs are related to the occurrence and development of an osteoporotic fracture. In order to further understand the potential molecular mechanism of hub immune-related mRNAs, lncRNA-miRNA-mRNA sub-network, hub mRNA-lncRNA co-expression sub-network and TFs network were constructed. Our analysis should provide some valuable information for future studies of osteoporotic fracture.