Constructing a competing endogenous RNA network for osteoarthritis
Introduction
The characteristics of the joint disorder of osteoarthritis (OA) include cartilage degeneration, secondary thickening of subchondral bones, and osteophyte formation (1). OA is a common disease worldwide. In 2005, about 26.9 million adults suffered from OA in the United States, which represents a marked increase from 5.9 million patients in 1990 (OA, 2015. Accessed 13 December 2015, at https://www.cdc.gov/arthritis/basics/osteoarthritis.htm). OA has been ranked as the 10th leading contributor to global years lived with disability (YLD) (1). Given that we live in an aging society, the incidence of OA will continue to increase. This is a big challenge for public health management.
Despite extensive research on its pathogenesis, the occurrence and progression mechanisms of OA remain obscure. Currently, several main underlying mechanisms are thought to explain the occurrence process of OA, including the accumulation of micro- and macro-injuries, and the abnormal activation of the repair response due to an injured joint, such as inflammation and innate immunity, which in turn leads to cell stress and extracellular matrix degradation. The mechanisms underlying the pathogenesis and progression of OA need to be further understood to enable the early diagnosis, better management, and assessment of OA.
Due to advances and developments in genetic biological methods, genetic alteration, long non-coding ribonucleic acid (RNA), and microRNA (miRNA) have been proven to be critical to OA occurrence and progression, and thus are diagnostic and therapeutic targets. Wang et al. found that lncRNA FOXD2 Adjacent Opposite Strand RNA 1 (FOXD2-AS1) serves as the protector for OA patients by inducing chondrocyte proliferation by sponging miR-27a-3p (2). Similarly, Cao et al. found that lncRNA FOXD2-AS1 inhibits the OA process by sponging miR-206 to regulate the expression of Cyclin D1 (CCND1) (3). Additionally, several lncRNAs promote OA progression; for example, LncRNA LOC101928134 acts as a promoter of OA, and the downregulation of lncRNA LOC101928134 can block the OA process (4).
mirRNAs also work as significant regulators of OA. MiR-132-3p acts as an inhibitor of OA by regulating the expression of ADAM metallopeptidase with thrombospondin type 1 motif 5 (ADAMTS-5) (5). Hyperbaric oxygen induces the upregulation of mir-107 in human osteoarthritic chondrocytes, which in turn inhibits the high mobility group box 1 (HMGB1)/advanced glycosylation end-product specific receptor (RAGE) signaling pathway. Conversely, miR-483-5p and microRNA-384-5p serve as promoters of OA (6,7). LncRNAs and miRNAs usually serve as significant regulators of disease by regulating gene expression levels. For example, miR-132-3p regulates the expression of ADAMTS-5, microRNA-145 inhibits MKK4, and lncRNA FOXD2-AS1 indirectly regulates the expression of CCND1 (3,5,8). OA progression also was deeply affected via OA patients microenvironments (9,10). And there are several approved that genes, miRNAs, and lncRNAs influence OA process through effecting the immune cells (11,12). Herein, we conducted a bioinformatics analysis to explore the effect of the ceRNA network on OA. Our analysis of the lncRNA/miRNA/mRNA network may inform novel or complementary theories on therapy methods for OA patients.
We present the following article in accordance with the STREGA reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-6711/rc).
Methods
Data collection
We downloaded the mRNA data set, GSE48556, the miRNA data set, GSE105027, and the lncRNA data set, GSE126963, from the Gene Expression Omnibus (GEO) database. GSE48556 contained 33 healthy samples, and 106 OA samples. GSE105027 contained 12 healthy samples and 12 OA samples. GSE126963 contained 3 healthy samples and 3 OA samples. As the data is relatively large, we use log2 logarithm processing, while the logarithm of values less than 1 will produce negative numbers, so a log2 +1 was used to normalize the data sets that had not been standardized
Analysis of DEGs
The healthy and OA samples were grouped, and the differentially expressed genes (DEGs) in the GSE48556, GSE105027, and GSE126963 data sets were analyzed using the limma package of R software. A log2 fold change |logFC| >0.3 and a P value <0.05 were set as the thresholds for DEGs. In the analysis of the miRNAs, the thresholds were set as |logFC| >0.1 and a P value <0.05, as there were fewer DEGs.
Construction of the ceRNA network
We took the differentially expressed miRNAs (DEmiRNAs) as the core of ceRNA network, and predicted the mRNAs and lncRNAs targeted by the DEmiRNAs through the miRcode, miRDB, miRtarbase, and TargetScan databases. Next, the genes of the ceRNA network were obtained by intersecting the targeted mRNAs and the DEmRNAs, and intersecting the targeted lncRNAs and the differentially expressed lncRNAs (DElncRNAs). Finally, Cytoscape software was used to construct the ceRNA network from the DEmiRNAs, differentially expressed mRNAs (DEmRNAs), and DElncRNAs in OA patients.
Functional enrichment analysis of DEmRNAs in OA
We performed a Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis of DEmRNAs in the GSE48556 data set via the clusterProfiler package in R, and used the ggplot2 package to present the results.
Gene expression levels of ceRNA network
Wilcoxon assays were conducted to analyze the gene expression levels of the ceRNA network between healthy and OA samples.
Immune infiltration analysis
A single-sample gene set enrichment analysis was conducted to examine the abundance of immune cells in healthy and OA patients based on the mRNA expression of the matrix, and 28 types of immune cells were identified (13). We compared the infiltration of 28 immune cells between the healthy and OA samples. We also analyzed the relationship between the abundance of immune cells and mRNA expression levels in the ceRNA network.
Statistical analysis
All the statistical analyses were completed using R software (version 4.0.3). The LIMMA package was used for the differential gene analysis, and the thresholds were set to |logFC|>0.3 and a P value <0.05. The Wilcoxon non-parameter test was used for comparisons between 2 groups, and the Spearman method was used for the correlation analysis. A P value <0.05 was considered statistically significant.
Ethical statement
The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Results
Results of the DEmRNA, DEmiRNA, and DElncRNA analyses
We identified 52 upregulated mRNAs and 127 downregulated mRNAs in the GSE48556 data set (Figure 1A,1B), 10 upregulated miRNAs and 25 downregulated miRNAs in the GSE105027 data set (Figure 1C,1D), and 1,239 upregulated lncRNAs and 1,209 downregulated lncRNAs in the GSE126963 data set (Figure 1E,1F).
The ceRNA network of OA
The target lncRNAs of 35 DEmiRNAs were predicted by the miRcode database, and only DElncRNAs with related pairs were retained. In total, 1 miRNA (hsa-mir-425-3p) and 24 lncRNAs were identified. We also predicted 53 target mRNAs of hsa-mir-425-3p via the miRDB, miRtarbase and TargetScan databases. We then took the intersection of 53 target mRNAs and 179 DEmRNAs to identify dual-specificity phosphatase 1 (DUSP1) (Figure 2A). Ultimately, 1 miRNA, 1 mRNA, and 24 lncRNAs were included in the ceRNA network (Figure 2B).
Expression levels of the ceRNA network genes
In the GSE48556 data set, the expression level of DUSP1 of the OA samples was lower than that of the healthy samples (Figure 3A). In the GSE105027 data set, hsa-mir-425-3p was also lowly expressed in the OA samples (Figure 3B). Among the 24 lncRNAs of the ceRNA network, in the OA sample of the GSE126963 data set, 17 lncRNAs were upregulated and 7 lncRNAs were downregulated (Figure 3C).
Functional enrichment analysis of DEmRNAs
The biological processes (BPs) of DEmRNAs were mainly enriched in the regulation of hemopoiesis, myeloid cell differentiation, positive regulation of hemopoiesis, leukocyte migration regulation, and negative anion transmembrane transport regulation. The cell components (CC) were mainly enriched in the nuclear speck, spliceosomal complex, platelet alpha granule, exon-exon junction complex, and cell body membrane. The molecular functions (MFs) were mainly enriched in Nuclear Factor Kappa B (NF-kappa B) binding, mitogen-activated protein kinase binding, protein tyrosine/threonine phosphatase activity, MAP kinase tyrosine/serine/threonine phosphatase activity, and MAP kinase phosphatase activity (Figure 4A,4B).
The KEGG signaling pathways of the DEmRNAs were mainly enriched in the cytokine-cytokine receptor interactions, pathogenic Escherichia coli infections, interleukin (IL)-17 signaling pathway, tumor necrosis factor (TNF) signaling pathway, osteoclast differentiation, apoptosis, measles, legionellosis, B cell receptor signaling pathway, and malaria (Figure 4C,4D).
Abundance of immune cell infiltrates in healthy and OA samples
The volcanic diagram of the differential expression analysis indicated that compared to the healthy samples, T follicular helper cells, type 17 T helper cells, natural killer cells and cluster of differentiation (CD)56 dim natural killer cells were highly infiltrated in OA patients, while plasmacytoid dendritic cells were lowly infiltrated in OA patients (Figure 5).
Correlation analysis between expression levels of DUSP1 and immune cells
The correlation analysis results showed that the high expression of DUSP1 leads to the low infiltration of activated B cells (Figure 6A), but the high infiltration of activated dendritic cells, central memory CD8 T cells, eosinophils, gamma delta T cells, macrophages, mast cells, monocytes, neutrophils, plasmacytoid dendritic cells, and regulatory T cells (Figure 6B-6K).
Discussion
Due to the aging population and the increasing numbers of obesity and chronic sports injuries, OA has become the one of the most common diseases to cause YLD. Currently, most diagnostic OA methods are only based on images and patients’ symptoms; thus, most OA patients are not diagnosed early, and late and terminal stage of patients can only undergo surgery to rehabilitate their daily activity. Molecules (e.g., lnRNAs, circular RNAs, miRNAs, and mRNA) in various diseases, including OA, have been proven to have a significant effect on OA progression (2-5). Several molecules, such as fatty acid-binding protein 4, transglutaminase 2, and circ_101178, have found to have significantly different expression levels between OA and normal samples (14-16). Thus, the alteration of molecules not only affects the progression OA, but can also serve as early diagnostic biomarkers. Further, the molecules also interact with each other. We conducted an analysis of a ceRNA network to identify the DElnRNAs, DEmiRNAs, and DEmRNAs between OA and normal samples, and further examined the interactions among them. From this ceRNA network, we can provide the some useful novel early biomarkers for diagnosing even as the therapy targets for OA patients.
In the present study, we identified 52 upregulated mRNAs and 127 downregulated mRNAs in the GSE48556 data set, 10 upregulated miRNAs and 25 downregulated miRNAs in the GSE105027data set, and 1,239 upregulated lncRNAs and 1,209 downregulated lncRNAs in the GSE126963 data set. To further examine the correlations among the differential expressions of the lncRNAs, miRNAs, and mRNAs, we constructed a ceRNA network. Ultimately, 1 miRNA, 1 mRNA, and 24 lncRNAs were enrolled in the ceRNA network. To date, no study appears to have been conducted in relation to hsa-mir-425-3p in OA patients; however, several studies have investigated its role in depression, Alzheimer’s disease, and aortic dissection (17-19). DUSP1 acts as a critical regulator in several OA processes. Notably, Choi et al. found that DUSP1 regulates RANKL expression and thus affects bone formation in OA patients (20). Another study showed that DUSP1 may act as the protector of OA fibroblast-like synoviocytes in OA patients by inhibiting the MAPK signaling pathway (21). Interestingly, Pest et al. showed that the deletion of DUSP1 did not increase the incidence of OA in mice (22). Thus, DUSP1 may have different biological functions in different conditions.
In our research, we identified 24 lncRNAs, including FAM182A, SNHG12, AL772363.1, AC025278.1, CLLU1, AL137145.1 AC011497.1, AL645728.1, MIR17HG, TTTY5, C7orf65, E2F3-IT1, TMEM72-AS1, FAM66E, SOX21-AS1, SUCLA2-AS1, LMO7-AS1, CHODL-AS1, SNHG7, LINC00472, LINC00427, CECR3, LINC00534, and MEG8. Xie et al. revealed that lncRNA MEG8 was downregulated in OA patients, and affected the progression of OA via the regulation of chondrocyte cells (23). LncRNA SNHG7 has also been shown to be downregulated in OA patients, and its main biological function is as a critical regulator of chondrocyte cell proliferation, apoptosis, and autophagy (24). In another study, Xu et al. found that lncRNA SNHG7 can act as an inhibitor of IL-1β-induced OA by sponging the miR-214-5p-mediated PPARGC1B signaling pathways (25). Thus, OA patients have several DElncRNAs, and these can have strong effects on OA. Regrettably, the roles of multiple lncRNAs have not yet been classified in OA.
To understand the main biological functions of our ceRNA network, we performed a functional enrichment analysis. The results showed that certain lnRNAs, miRNAs, and mRNAs take part in various significant BPs, such as the regulation of leukocyte migration, MAP kinase tyrosine/serine/threonine phosphatase activity, the IL-17 signaling pathway, the TNF signaling pathway, and osteoclast differentiation. Leukocytes are the most critical factors in the inflammation process. Chibber et al. found that OA-Dehydrozingerone (DHZ) is an anti-inflammation agent that inhibits leukocyte migration and thus provides an obstacle to the OA process (26). IL-17 has been shown to be abnormally expressed in OA patients (27,28). Zhou et al. found that increases in IL-17 expression leads to aggravated cartilage injury via PI3K/AKT/mTOR activation and thus regulates the ciRS-7/ miR-7 axis (28). Faust et al. confirmed that inhibiting IL-17 expression can protect OA patients by reducing joint degeneration (27). The TNF signaling pathway has been proven to be an important pathway in the regulation of the OA process. Chen et al. showed that spermidine promoted RIP1 deubiquitination to inhibit the TNF-α-induced NF-κB/p65 pathway, which in turn inhibited the progression of OA (29).
In addition to inflammation, bone formation and degradation are also critical processes of OA. Several studies have shown that osteoclast differentiation is a significant process in OA progression (30,31). Several functions of the ceRNA network and the roles in OA have not yet been investigated. Thus, further research is required to confirm the role of the ceRNA network in OA patients.
The inflammation process occurs throughout the whole OA process. Immune cells are the main cause of inflammation. We also examined differences in immune cell infiltration between OA and healthy samples, and found that several immune cells, such as T follicular helper cells, type 17 T helper cells, natural killer cells, and CD56 dim natural killer cells, were more abundant in OA patients than healthy individuals. We also found that hub mRNA DUSP1 expression levels affect immune cell infiltration. Thus, lncRNAs, miRNAs, and mRNAs not only directly affect OA but also regulate OA by regulating immune cell infiltration. However, shortage of validation experiments is one of limitation of our study.
Conclusions
In our study, we constructed and analyzed a DUSP1-specific ceRNA network. We not only found the potential underlying pathogenesis of OA but also found that DElncRNAs, DEmiRNAs, and DEmRNAs can be used as early diagnostic biomarkers for OA patients. A limitation of our study is that we did not conduct any experiments to validate the DElncRNA, DEmiRNA, and DEmRNA expression levels and their biological roles in OA.
Acknowledgments
Funding: None.
Footnote
Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-21-6711/rc
Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-21-6711/coif). The authors have no conflicts of interest to declare.
Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.
References
- Wang Q, Wu Y, Li F, et al. Association of genetic polymorphisms in immune-related lncRNA with osteoarthritis susceptibility in Chinese Han population. Per Med 2018;15:103-10. [Crossref] [PubMed]
- Wang Y, Cao L, Wang Q, et al. LncRNA FOXD2-AS1 induces chondrocyte proliferation through sponging miR-27a-3p in osteoarthritis. Artif Cells Nanomed Biotechnol 2019;47:1241-7. [Crossref] [PubMed]
- Cao L, Wang Y, Wang Q, et al. LncRNA FOXD2-AS1 regulates chondrocyte proliferation in osteoarthritis by acting as a sponge of miR-206 to modulate CCND1 expression. Biomed Pharmacother 2018;106:1220-6. [Crossref] [PubMed]
- Yang DW, Zhang X, Qian GB, et al. Downregulation of long noncoding RNA LOC101928134 inhibits the synovial hyperplasia and cartilage destruction of osteoarthritis rats through the activation of the Janus kinase/signal transducers and activators of transcription signaling pathway by upregulating IFNA1. J Cell Physiol 2019;234:10523-34. [Crossref] [PubMed]
- Zhou X, Luo D, Sun H, et al. MiR-132-3p regulates ADAMTS-5 expression and promotes chondrogenic differentiation of rat mesenchymal stem cells. J Cell Biochem 2018;119:2579-87. [Crossref] [PubMed]
- Wang H, Zhang H, Sun Q, et al. Intra-articular Delivery of Antago-miR-483-5p Inhibits Osteoarthritis by Modulating Matrilin 3 and Tissue Inhibitor of Metalloproteinase 2. Mol Ther 2017;25:715-27. [Crossref] [PubMed]
- Zhang W, Cheng P, Hu W, et al. Inhibition of microRNA-384-5p alleviates osteoarthritis through its effects on inhibiting apoptosis of cartilage cells via the NF-κB signaling pathway by targeting SOX9. Cancer Gene Ther 2018;25:326-38. Erratum in: Cancer Gene Ther 2020;27:836-7. [Crossref] [PubMed]
- Hu G, Zhao X, Wang C, et al. MicroRNA-145 attenuates TNF-α-driven cartilage matrix degradation in osteoarthritis via direct suppression of MKK4. Cell Death Dis 2017;8:e3140. [Crossref] [PubMed]
- Woodell-May JE, Sommerfeld SD. Role of Inflammation and the Immune System in the Progression of Osteoarthritis. J Orthop Res 2020;38:253-7. [Crossref] [PubMed]
- Lopes EBP, Filiberti A, Husain SA, et al. Immune Contributions to Osteoarthritis. Curr Osteoporos Rep 2017;15:593-600. [Crossref] [PubMed]
- Zhao C. Identifying the hub gene and immune infiltration of osteoarthritis by bioinformatical methods. Clin Rheumatol 2021;40:1027-37. [Crossref] [PubMed]
- Li Z, Yuan B, Pei Z, et al. Circ_0136474 and MMP-13 suppressed cell proliferation by competitive binding to miR-127-5p in osteoarthritis. J Cell Mol Med 2019;23:6554-64. [Crossref] [PubMed]
- Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
- Chen C. Serum hsa_circ_101178 as a Potential Biomarker for Early Prediction of Osteoarthritis. Clin Lab 2020; [Crossref] [PubMed]
- Tarantino U, Ferlosio A, Arcuri G, et al. Transglutaminase 2 as a biomarker of osteoarthritis: an update. Amino Acids 2013;44:199-207. [Crossref] [PubMed]
- Zhang C, Li T, Chiu KY, et al. FABP4 as a biomarker for knee osteoarthritis. Biomark Med 2018;12:107-18. [Crossref] [PubMed]
- Dakterzada F, Targa A, Benítez ID, et al. Identification and validation of endogenous control miRNAs in plasma samples for normalization of qPCR data for Alzheimer's disease. Alzheimers Res Ther 2020;12:163. [Crossref] [PubMed]
- Ji J, Xu Q, He X, et al. MicroRNA microarray analysis to detect biomarkers of aortic dissection from paraffin-embedded tissue samples. Interact Cardiovasc Thorac Surg 2020;31:239-47. [Crossref] [PubMed]
- Serafini G, Trabucco A, Amerio A, et al. Commentary on the study of Roy et al. Amygdala Based Altered mir-128-3p in Conferring Susceptibility to Depression-like Behavior via Wnt Signalling. Int J Neuropsychopharmacol 2020;23:178-80. [Crossref] [PubMed]
- Choi Y, Yoo JH, Lee Y, et al. Calcium-Phosphate Crystals Promote RANKL Expression via the Downregulation of DUSP1. Mol Cells 2019;42:183-8. [PubMed]
- Peng HZ, Yun Z, Wang W, et al. Dual specificity phosphatase 1 has a protective role in osteoarthritis fibroblast-like synoviocytes via inhibition of the MAPK signaling pathway. Mol Med Rep 2017;16:8441-7. [Crossref] [PubMed]
- Pest MA, Pest CA, Bellini MR, et al. Deletion of Dual Specificity Phosphatase 1 Does Not Predispose Mice to Increased Spontaneous Osteoarthritis. PLoS One 2015;10:e0142822. [Crossref] [PubMed]
- Xie W, Jiang L, Huang X, et al. lncRNA MEG8 is downregulated in osteoarthritis and regulates chondrocyte cell proliferation, apoptosis and inflammation. Exp Ther Med 2021;22:1153. [Crossref] [PubMed]
- Tian F, Wang J, Zhang Z, et al. LncRNA SNHG7/miR-34a-5p/SYVN1 axis plays a vital role in proliferation, apoptosis and autophagy in osteoarthritis. Biol Res 2020;53:9. [Crossref] [PubMed]
- Xu J, Pei Y, Lu J, et al. LncRNA SNHG7 alleviates IL-1β-induced osteoarthritis by inhibiting miR-214-5p-mediated PPARGC1B signaling pathways. Int Immunopharmacol 2021;90:107150. [Crossref] [PubMed]
- Chibber P, Kumar C, Singh A, et al. Anti-inflammatory and analgesic potential of OA-DHZ; a novel semisynthetic derivative of dehydrozingerone. Int Immunopharmacol 2020;83:106469. [Crossref] [PubMed]
- Faust HJ, Zhang H, Han J, et al. IL-17 and immunologically induced senescence regulate response to injury in osteoarthritis. J Clin Invest 2020;130:5493-507. [Crossref] [PubMed]
- Zhou X, Li J, Zhou Y, et al. Down-regulated ciRS-7/up-regulated miR-7 axis aggravated cartilage degradation and autophagy defection by PI3K/AKT/mTOR activation mediated by IL-17A in osteoarthritis. Aging (Albany NY) 2020;12:20163-83. [Crossref] [PubMed]
- Chen Z, Lin CX, Song B, et al. Spermidine activates RIP1 deubiquitination to inhibit TNF-α-induced NF-κB/p65 signaling pathway in osteoarthritis. Cell Death Dis 2020;11:503. [Crossref] [PubMed]
- Sabokbar A, Crawford R, Murray DW, et al. Macrophage-osteoclast differentiation and bone resorption in osteoarthrotic subchondral acetabular cysts. Acta Orthop Scand 2000;71:255-61. [Crossref] [PubMed]
- Xiao M, Hu ZH, Jiang HH, et al. Role of osteoclast differentiation in the occurrence of osteoarthritis of temporomandibular joint. Hua Xi Kou Qiang Yi Xue Za Zhi 2021;39:398-404. [PubMed]
(English Language Editor: L. Huleatt)