Constructing a competing endogenous RNA network for osteoarthritis.

Background: Osteoarthritis (OA) is one of the most common diseases in elderly people; however, the correlation between molecular alterations and the occurrence and progression of OA are still not well understood. We conducted this study to investigate the molecular changes in OA via the competing endogenous ribonucleic acid (ceRNA) network.
Methods: We downloaded the messenger RNA (mRNA) data set, GSE48556, the microRNA (miRNA) data set, GSE105027, and the long non-coding (lncRNA) data set, GSE126963 from the Gene Expression Omnibus (GEO) database, and examined the differentially expressed genes (DEGs) in these data sets. Further, we constructed a ceRNA network of the differentially expressed miRNAs, mRNAs, and lncRNAs. To determine the biological functions of the ceRNA network, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses. Finally, we conducted an immune cell infiltration analysisusing single-sample gene set enrichment analysis to examine the abundance of immune cells in healthy and OA patients, and 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.
Results: Ultimately hsa-mir-425-3p, dual specificity phosphatase 1, and 24 lncRNAs were identified in the ceRNA network. The functional enrichment analyses showed that these lncRNAs, miRNAs, and mRNAs are involved in various significant biological process, such as the regulation of leukocyte migration, Mitogen-Activated Protein (MAP) kinase tyrosine/serine/threonine phosphatase activity, the interleukin-17 signaling pathway, the tumor necrosis factor signaling pathway, and osteoclast differentiation, and can also have a strong effect on immune cell infiltration. Conclusions: The dual-specificity phosphatase 1-specific ceRNA network can be used as a diagnostic tool to assess the progression of OA patients. 2022 Annals of Translational Medicine. All rights reserved.

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 (), 10 upregulated miRNAs and 25 downregulated miRNAs in the GSE105027 data set (), and 1,239 upregulated lncRNAs and 1,209 downregulated lncRNAs in the GSE126963 data set ().

Figure 1

Results of DEmRNA, DEmiRNA, and DElncRNA analysis of 3 data sets. (A) Heatmap of top 50 rank up/downregulated mRNAs in OA; (B) volcanic plot of DEmRNAs in OA; (C) heatmap of all up/downregulated miRNAs in OA; (D) volcanic plot of DEmiRNAs in OA; (E) heatmap of top 50 rank up/downregulated lncRNAs in OA; (F) volcanic plot of DElncRNAs in OA. DEmRNAs, differentially expressed mRNAs; DEmiRNAs, differentially expressed miRNAs; DElncRNAs, differentially expressed lncRNAs.

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) (). Ultimately, 1 miRNA, 1 mRNA, and 24 lncRNAs were included in the ceRNA network ().

Figure 2

The competing endogenous RNA network of OA. (A) Venn diagram of the intersection of DEmRNAs and the target mRNA of hsa-mir-425-3p in the GSE48556 data set; (B) the ceRNA network of OA; red represents miRNAs, blue represents mRNAs, and green represents lncRNAs. DEmRNAs, differentially expressed mRNAs; OA, osteoarthritis.

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 (). In the GSE105027 data set, hsa-mir-425-3p was also lowly expressed in the OA samples (). 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 3

Expression level of ceRNA network genes of OA. (A) Expression level between healthy and OA samples of DUSP1 in the GSE48556 data set; (B) expression level between healthy and OA samples of hsa-mir-425-3p in the GSE105027 data set; (C) expression level between healthy and OA samples of 24 lncRNAs in the GSE48556 data set. ceRNA, competing endogenous RNA; OA, osteoarthritis; DUSP1, dual specificity phosphatase 1; lncRNAs, long non-coding RNAs.

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 4

Functional enrichment of DEmRNAs. (A) Bubble diagram of GO functional enrichment analysis of DEmRNAs; (B) circle diagram of the DEmRNAs specifically enriched in biological processes; (C) bubble diagram of the top 10 KEGG signaling pathways; (D) circle diagram of DEmRNAs specifically enriched in the signaling pathways. The size of the bubble represents the number of enrichment DEmRNAs, and the color represents the q value and P value. DEmRNAs, differentially expressed RNAs; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

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 ().

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

Abundance of 28 types of immune cell infiltrates in healthy samples and OA patients. OA, osteoarthritis.

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 (), 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 6

Correlation analysis between the expression level of DUSP1 and immune cells. (A) Activated B cells; (B) activated dendritic cells; (C) central memory CD8 T cells; (D) eosinophils; (E) gamma delta T cells; (F) macrophages; (G) mast cells; (H) monocytes; (I) neutrophils; (J) plasmacytoid dendritic cells; (K) regulatory T cells. DUSP1, dual specificity phosphatase 1.

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.

The article’s supplementary files as