Abstract

Total hip arthroplasty (THA) is a cost-effective treatment for osteoarthritis (OA), and osteolysis is a common complication of THA. This study was aimed at exploring the relevant molecular biomarkers for osteolysis after THA. We performed RNA sequence to identify and characterize expressed mRNAs and lncRNAs in OA and osteolysis. Differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) in OA and osteolysis were acquired, as well as shared DEmRNAs/DElncRNAs in OA and osteolysis and osteolysis-specific DEmRNAs/DElncRNAs. Then, shared and osteolysis-specific DElncRNA-DEmRNA coexpression networks were constructed to further investigate the function of DElncRNAs and DEmRNAs in OA and osteolysis. In total, 343 DEmRNAs and 25 DElncRNAs in OA, 908 DEmRNAs and 107 DElncRNAs in osteolysis, and 406 DEmRNAs and 46 DElncRNAs between OA and osteolysis were acquired. A total of 136 shared DEmRNAs and 9 shared DElncRNAs in OA and osteolysis and 736 osteolysis-specific DEmRNAs and 103 osteolysis-specific DElncRNAs were acquired. Then, 128 shared DElncRNA-DEmRNA coexpression pairs and 522 osteolysis-specific DElncRNA-DEmRNA coexpression pairs were identified. The present study highlighted the roles of four interaction pairs, including two shared lncRNA-mRNA interaction pairs in OA and osteolysis (AC111000.4 and AC016831.6), which may function in the immune process of OA and osteolysis by regulating CD8A and CD8B, respectively, and two osteolysis-specific interaction pairs (AC090607.4-FOXO3 and TAL1-ABALON), which may play an important role in osteoclastogenesis.

1. Introduction

Osteoarthritis (OA) is a leading cause of chronic disability in old people. For the treatment, the total hip arthroplasty (THA) is a cost-effective way, which can reduce joint pain, restore joint function, and increase the quality of life of patients [1]. Despite the improvement of the quality of polyethylene, osteolysis remains a risk for older designs and younger, active patients. Osteolysis is a progressive, active, biologic cascade, a phenomenon due to a foreign body response to particulate wear debris from the prosthetic joint [2]. Increased wear particles activate osteoclast formation, and overweight osteoclasts caused much bone resorption, which eventually led to osteolysis [3]. Osteolysis, as a complication of THA, leads to prosthesis failure and bringing about additional suffering and burden for patients [4].

Long noncoding RNAs (lncRNAs), as a type of noncoding RNA, have been recognized as key regulatory molecules with diverse roles in gene expression, epigenetic modification, and protein activity [5]. Recently, lncRNA has been revealed to be involved in osteolysis. lncRNA TSIX was involved in particle-induced osteolysis by regulating miR-30a-5p to promote osteoblast apoptosis [6]. lncRNA DANCR inhibits osteoblast differentiation in osteolysis after THA through regulating FOXO1 [7]. lncRNA KCNQ1OT1 could ameliorate particle-induced osteolysis, by inhibiting miR-21a-5p to induce macrophage polarization [8]. However, more studies focused on the role of lncRNAs in osteolysis after THA need to be performed.

In this study, we, respectively, investigated the gene expression profiles of lncRNAs and mRNAs in patients with osteolysis after THA and OA, attempting to screen out differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) associated with osteolysis after THA and OA. The objective of this study was to explore the underlying mechanism and the relevant molecular biomarkers for osteolysis after THA.

2. Materials and Methods

2.1. Subjects and Samples

The cohort subjected to RNA-Seq comprised 3 patients with OA, 3 patients with osteolysis after THA, and 3 healthy individuals. Inclusion criteria for patients with osteolysis after THA: (1) with current radiographic evidence of osteolysis and (2) received a THA after failure to improve function and pain after at least 6 months of conservative treatment. Subjects were excluded if they had any history of inflammatory arthropathy or known secondary causes of hip arthritis such as trauma, avascular necrosis, or developmental or childhood hip disease. Subjects were also excluded if they had taken courses of immunosuppressant agents or bisphosphonates for a continuous period of greater than 6 months since THA. OA was diagnosed according to the criteria of the American College of Rheumatology. Healthy individuals with no personal or family history of OA, no symptoms or signs of OA, or any other type of arthritis, or any painful condition of the joints, were included as controls. The participants with history of joint diseases, including inflammatory arthritis (rheumatoid arthritis or any other autoimmune disease), posttraumatic or postseptic arthritis, poliomyelitis, and skeletal dysplasia, were excluded. Table 1 described the characteristics of all these participants. All samples were collected after obtaining written informed consent from every participant. This study was approved by the Ethics Committee of Shandong Provincial Hospital (No. 2020-123) and performed in accordance with the Declaration of Helsinki. Peripheral whole blood (2.5 mL) drawn from each subject was collected in PAXgene® RNA blood tubes and stored at -80°C prior to processing. RNA isolation was performed with PAXgene blood RNA kit. RNA integrity and concentration were evaluated with an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit). Total RNA samples used in subsequent experiments fulfilled the following requirements: RNA integrity and . In brief, total RNA was subjected to ribosomal RNA (rRNA) removal using Ribo-Zero. A total amount of 3 μg RNA was used for library preparation. Libraries for sequencing were constructed according to the manufacturer’s protocol. The quality of the libraries was determined using an Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System. Based on the BGIseq platform, sequencing was performed. The raw sequencing data were submitted to sequencing quality control by FastQC to assess whether they will be used for subsequent data analysis. Reads with low quality (adaptor sequences, sequences with a quality , and sequences with an N base rate of raw %) were removed. Clean reads were aligned with the human reference genome, Ensemble GRCh38.p7.

2.2. Identification of DEmRNAs and Functional Annotation

StringTie software was used to compare the results to the known transcriptome and calculate the transcriptional abundance. Then, Ballgown was used for quantification of gene expression levels, as well as analysis of intergroup differential expression of genes. DESeq2 was applied to identify DEmRNAs in OA vs. control and osteolysis after THA vs. control with a value < 0.05. DEmRNAs between OA and osteolysis after THA were obtained with a value < 0.05 as well. Hierarchical clustering analysis of DEmRNAs was performed with R package “pheatmap.” Then, shared DEmRNAs in OA and osteolysis after THA and osteolysis-specific DEmRNAs (DEmRNAs in osteolysis after THA but no differences in OA) were further identified with Venny 2.1.0. The Database for Annotation, Visualization and Integrated Discovery (DAVID), which is a web-based tool, provides integrated solutions for the annotation and analysis of genome-scale datasets from high-throughput sequencing. DAVID 6.8 was used to perform GO and KEGG enrichment analysis for shared DEmRNAs in OA and osteolysis after THA and osteolysis-specific DEmRNAs. The lower the value, the more significant are the GO term and the pathway. A value of was considered to be represented statistically significant.

2.3. Identification of DElncRNAs

DESeq2 was applied to identify DElncRNAs in OA vs. control and osteolysis after THA vs. control with a value < 0.05 and . DElncRNAs between OA and osteolysis after THA were obtained with a value < 0.05 and as well. Hierarchical clustering analysis of DElncRNAs was performed with R package “pheatmap.” Then, shared DElncRNAs in OA and osteolysis after THA and osteolysis-specific DElncRNAs (DElncRNAs in osteolysis after THA but no differences in OA) were further identified with Venny 2.1.0.

2.4. DElncRNA-DEmRNA Coexpression Network and Functional Annotation

The shared and osteolysis-specific DElncRNA-DEmRNA coexpression networks were constructed to further investigate the potential functions of lncRNAs and mRNAs in OA and osteolysis. Pearson correlation coefficients were calculated between the expression values of DElncRNAs and DEmRNAs. The pairs with and in shared DElncRNA-DEmRNA pairs and pairs with and in osteolysis-specific DElncRNA-DEmRNA pairs were defined as coexpressed DElncRNA-DEmRNA pairs, respectively. Then coexpressed networks were visualized by using Cytoscape. DAVID 6.8 was used to perform GO and KEGG enrichment analysis for DEmRNAs in shared and osteolysis-specific DElncRNA-DEmRNA coexpression network. A value of was considered to be represented statistically significant.

2.5. Statistical Analysis

Statistical analyses were performed using R software v3.5.3 (R Foundation for Statistical Computing, Vienna, Austria). All tests were two-tailed, and values < 0.05 were considered statistically significant. In addition, the GraphPad Prism Software, version 7.0, was used for the statistical analysis of experimental data. The results are expressed as (SDs). In addition, we used Shapiro-Wilk to test the normal distribution of data. However, the statistical analysis of clinical data showed that the data are not normal distribution. Therefore, we choose the Kruskal-Wallis test for statistical analysis. Pearson correlation coefficients were calculated between the expression values of DElncRNAs and DEmRNAs. was considered to indicate a significant difference between the groups.

3. Results

3.1. Identification of DEmRNAs and Functional Annotation

DESeq2 was applied to identify DEmRNAs, and our results showed that a total of 343 DEmRNAs (184 up- and 159 downregulated) in OA vs. control, 908 DEmRNAs (429 up- and 479 downregulated) in osteolysis vs. control, and 406 DEmRNAs (112 up- and 294 downregulated) in OA vs. osteolysis were identified. Of these, HIST2H4A and CD8A, ZNF778 and TMEM56-RWDD3, and EPHB4 and MTURN were the most up- and downregulated DEmRNAs in OA vs. control, osteolysis vs. control, and OA vs. osteolysis, respectively (Table 2). The heatmap of the top 100 up- and downregulated DEmRNAs was shown in Figures 1(a)1(c). The volcano plots of DEmRNAs are shown in Figures S1A–C. A total of 136 shared DEmRNAs (71 up- and 65 downregulated) in OA vs. control and osteolysis vs. control and 736 osteolysis-specific DEmRNAs (381 up- and 355 downregulated) were acquired (Figure 1(d)).

To investigate the functions of shared and osteolysis-specific DEmRNAs, DAVID 6.8 was used to perform GO and KEGG enrichment analysis. For shared DEmRNAs, cell activation (), leukocyte activation (), integral to plasma membrane (), and protein homodimerization activity () were several significantly enriched GO terms, and mTOR signaling pathway () and cell cycle () were significantly enriched KEGG pathways (Figures S2, S4 A and B). For osteolysis-specific DEmRNAs, intracellular signaling cascade (), cell adhesion (), plasma membrane (), and cytoskeletal protein binding () were several significantly enriched GO terms, and hematopoietic cell lineage (), cell adhesion molecules (CAMs) (), porphyrin and chlorophyll metabolism (), and systemic lupus erythematosus () were significantly enriched KEGG pathways (Figures S3, S4 C–F).

3.2. Identification of DElncRNAs

DESeq2 was applied to identify DElncRNAs, and our results showed that a total of 25 DElncRNAs (15 up- and 10 downregulated) in OA vs. control, 107 DElncRNAs (58 up- and 49 downregulated) in osteolysis vs. control, and 46 DElncRNAs (17 up- and 29 downregulated) in OA vs. osteolysis were identified. Of these, NEAT1 and AC005726.2, FLJ42393 and AC123912.4, and AC016737.1 and AC123912.4 were the most up- and downregulated DElncRNAs in OA vs. control, osteolysis vs. control, and OA vs. osteolysis, respectively (Table 3). The heatmap of DElncRNAs is shown in Figures 2(a)2(c). The volcano plots of DEmRNAs are shown in Figure S1D-F. A total of 9 shared DElncRNAs (6 up- and 3 downregulated) in OA vs. control and osteolysis vs. control and 103 osteolysis-specific DElncRNAs (52 up- and 51 downregulated) were acquired (Figure 2(d)).

3.3. DElncRNA-DEmRNA Coexpression Network and Functional Annotation

To further investigate the potential functions of lncRNAs and mRNAs in OA and osteolysis, the shared and osteolysis-specific DElncRNA-DEmRNA coexpression networks were constructed. A total of 128 shared DElncRNA-DEmRNA coexpression pairs including 9 DElncRNAs and 81 DEmRNAs were obtained (Figure 3(a)). Among them, AC018761.2 (), AC090607.4 (), and ABALON () were the top 3 DElncRNAs that covered the most DEmRNAs. Then, to investigate the functions of DEmRNAs in shared and osteolysis-specific DElncRNA-DEmRNA coexpression network, DAVID 6.8 was used to perform GO and KEGG enrichment analysis. For DEmRNAs in shared DElncRNA-DEmRNA coexpression network, T cell activation (), leukocyte differentiation (), leukocyte activation (), immune system development (), T cell differentiation (), and mTOR signaling pathway () were several significantly enriched pathways (Figure 3(b)).

A total of 522 osteolysis-specific DElncRNA-DEmRNA coexpression pairs including 36 DElncRNAs and 194 DEmRNAs were obtained (Figure 4(a)). Among them, AC111000.4 (), OVCH1-AS1 (), and AC016831.6 () were the top 3 DElncRNAs that covered the most DEmRNAs. For DEmRNAs in osteolysis-specific DElncRNA-DEmRNA coexpression network, myeloid cell differentiation (), immune system development (), cell proliferation (), and positive regulation of myeloid cell differentiation () were several significantly enriched pathways (Figure 4(b)).

4. Discussion

Osteolytic lesions may develop after THA from a biologic reaction to particulate debris [9]. It is a major complication of THA, causing additional suffering and burden for patients. Aggravating evidence indicates that lncRNAs regulate gene expression via cis- and/or transregulation mechanisms and participate in various biological processes, including chromatin modification, DNA synthesis, cell proliferation, differentiation, and apoptosis [10]. In this present study, we screened out critical genes and lncRNAs correlated with OA and osteolysis by RNA sequencing and bioinformatics analysis.

Due to the close interaction between immune cells and bone cells, immune disorders can lead to abnormal bone metabolism [11]. T lymphocyte differentiation could be made by the “clusters of differentiation” (CD), and the common receptors are CD4 and CD8. CD8, a cell surface glycoprotein, is predominantly expressed on the surface of cytotoxic T killer cells and implicated in the immune system. The CD8 antigen is a specifically recognized antigen displayed by the class major histocompatibility complex I molecules. Long et al. identified a 16-gene biomarker panel to differentiate RA from OA and indicated the lower expression level of CD8A in OA than in RA [12]. Landgraeber et al. reported that the CD4+/CD8+ ratio was associated with the stage of osteolysis in aseptic loosening [13]. In this study, both CD8A and CD8B were significantly decreased in OA vs. control and osteolysis vs. control. Functional annotation revealed that CD8A and CD8B were enriched in immune-related pathways, such as T cell activation, leukocyte differentiation, leukocyte activation, immune system development, and T cell differentiation. In addition, CD8A was coexpressed with AC111000.4, and CD8B was coexpressed with AC016831.6 in shared the DElncRNA-DEmRNA coexpression network. Taken together, we speculated that AC111000.4 and AC016831.6 may function in the immune process of OA and osteolysis by regulating CD8A and CD8B, respectively.

The FOXO proteins are an evolutionarily conserved family of transcription factors which comprised FOXO1, FOXO3, FOXO4, and FOXO6 in mammals [14]. The four FOXO members share obvious sequence homology and are ubiquitously expressed in various organs, including bone [15]. FOXOs regulate diverse cellular processes, including oxidative stress, metabolism, apoptosis, and inflammation [16]. In addition, FOXOs have been revealed to regulate bone cell functions, osteoblast differentiation, and the maintenance of skeletal homeostasis [17]. Bartell et al. suggested that FOXO3 regulate receptor activator of NF-κB ligand- (RANKL-) induced osteoclast differentiation [18]. Miller et al. reported that FOXO3 play important inhibitory roles in TNF-α-mediated osteoclastogenesis and bone resorption [19]. Ambrogini et al. found that overexpression of FOXO3 in osteoblasts may reduce osteoclast numbers [20]. Increased osteoclastogenesis leads to osteolysis. Here, FOXO3 was determined to be a significantly decreased osteolysis-specific DEmRNA and coexpressed with AC090607.4, which indicated the importance of AC090607.4 in osteolysis by targeting FOXO3.

TAL1, also known as SCL, plays a crucial role in hematopoiesis and is causally connected to T cell acute lymphatic leukemia [21, 22]. The fact that embryonic stem cells of Tal1-knockout mice did not develop into osteoclasts in vitro supported that Tal1 plays a role in osteoclast differentiation [23]. The maintenance of bone homeostasis mainly depends on the balance between bone-forming osteoblasts and bone-resorbing osteoclasts. The transcription factor TAL1 was of great importance in cell cycle progression and proliferation in differentiating murine bone marrow monocyte precursors [24]. Knockdown of Tal1 in osteoclast progenitors leads to larger osteoclasts, containing more nuclei, and altered expression of a large number genes [25]. Wang et al. found that TAL1 was a differentially expressed transcription factor in osteoporosis [26]. In this study, TAL1 was coexpressed with ABALON in the osteolysis-specific DElncRNA-DEmRNA coexpression network. Therefore, the role of TAL1-BALON in osteolysis should be paid attention.

5. Conclusions

In conclusion, we highlighted the roles of four interaction pairs, including two shared lncRNA-mRNA interaction pairs in OA and osteolysis (AC111000.4-CD8A and AC016831.6-CD8B) and two osteolysis-specific interaction pairs (AC090607.4-FOXO3 and TAL1-BALON) in this present work. Taken as a whole, our work made contribution to understanding the pathophysiology of osteolysis. Our study also had a limitation. The sample size for RNA sequencing in this study was small, and further studies with larger sample size are warranted to confirm these results.

Data Availability

The clinical data used to support the findings of this study are restricted by the Ethics Committee of Shandong Provincial Hospital (No. 2020-123) in order to protect patient privacy. Data are available from Guang Yang ([email protected]) for researchers who meet the criteria for access to confidential data.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of Shandong (No. ZR2020QH073) and National Natural Science Foundation of China (No. 81972056).

Supplementary Materials

Figure S1: the volcano plots of DEmRNAs in OA (A), osteolysis (B), OA vs osteolysis (C), DElncRNAs in OA (D), osteolysis (E), and OA vs osteolysis (F). Figure S2: significantly enriched GO terms (A) and KEGG (B) pathways of shared DEmRNAs in OA and osteolysis after THA. Figure S3: significantly enriched GO terms (A) and KEGG (B) pathways of osteolysis-specific DEmRNAs. Figure S4: mTOR signaling pathway and cell cycle were significantly enriched KEGG pathways of shared DEmRNAs (A, B). Hematopoietic cell lineage, cell adhesion molecules (CAMs), porphyrin and chlorophyll metabolism, and systemic lupus erythematosus were significantly enriched KEGG pathways of osteolysis-specific DEmRNAs (C–F). (Supplementary Materials)