A Long Noncoding RNA (lncRNA)-Associated Competing Endogenous RNA (ceRNA) Network Identifies Eight lncRNA Biomarkers in Patients with Osteoarthritis of the Knee.

BACKGROUND Osteoarthritis (OA) of the knee is a common disease that is associated with chronic pain. This study aimed to identify and investigate the functional role of biomarkers associated with long noncoding RNA (lncRNA) in the progression of OA of the knee by lncRNA-associated competing endogenous RNA (ceRNA) integrated network analysis. MATERIAL AND METHODS High-quality microRNA (miRNA)-lncRNA and miRNA-mRNA interactions and lncRNA and mRNA expression profiles for patients with OA of the knee with mild and severe pain were obtained from the Gene Expression Omnibus (GEO) database (GSE99662). A three-step computational method was used to construct the lncRNA-associated ceRNA interaction network in OA by integrating miRNA-lncRNA/mRNA interactions and lncRNA/mRNA expression profiles in patients with OA with mild and severe pain. RESULTS A total of 1,870 dysregulated lncRNA-mRNA interactions were obtained in the lncRNA-associated ceRNA network in OA, including 476 gain and 1,394 loss interactions, covering 131 lncRNAs and 1,251 mRNAs. Characterization of the lncRNA-associated ceRNA network in OA indicated that lncRNAs had roles in the network. Further differential expression analysis identified eight lncRNA biomarkers, which could distinguish between patients with OA with mild pain and severe pain. These lncRNA-associated interactions showed significantly different co-expression patterns in samples from patients with OA of the knee associated with mild pain. CONCLUSIONS Integrated network analysis of lncRNA-associated ceRNA identified eight lncRNA molecular biomarkers associated with the progression of OA of the knee.

Background

Primary osteoarthritis (OA) of the knee is a common chronic disease with increasing incidence and prevalence [1-3]. Effective therapies remain to be discovered to prevent the progression of OA and reduce its severity, and OA is predicted to become an increasing economic burden as the population ages [1,4,5]. Chronic pain from OA is a common public health problem that has a detrimental impact on patient health and function [6]. Chronic pain also reduces the quality of life in patients with OA. Therefore, an accurate assessment of the severity of OA is clinically important [7,8].

Currently, several studies have been published on the association between the expression of long noncoding RNA (lncRNA) and the progression of OA. The lncRNA has a length of more than 200 nucleotides and has widespread biological functions, including regulation of transcriptional patterns and protein activity [9,10]. Recent studies have shown that some lncRNAs are involved in the pathogenesis of OA [11-16]. Dudek et al. reported that lncRNA H19 was overexpressed in chondrocytes and was regulated by SOX9 [17]. Large-scale RNA sequencing studies in patients with OA of the knee with severe pain and mild pain showed that several lncRNAs were differentially expressed, which indicates that these dysregulated lncRNAs may be involved in the progression of OA [1,18]. However, the molecular mechanisms of the role of dysregulated lncRNAs in the progression of OA remain poorly understood and further functional studies are needed to determine these roles.

In 2011, the hypothesis was proposed that competing endogenous RNAs (ceRNAs) could act as miRNA sponges through their miRNA response elements (MREs) and regulate their expression [19]. Several studies have shown that lncRNAs can act as ceRNAs and have important roles in malignancy and chronic diseases, including OA [20-24]. Shen et al. found that the lncRNA-SNHG5 could function as a ceRNA to competitively sponge miR-26a to regulate the expression of mRNA SOX2, and had a role in OA [23]. Also, Li et al. identified the modulating effect of lncRNA X-inactive specific transcript (XIST) in OA [24]. These authors showed that lncRNA-XIST served as a ceRNA for miR-211 to counteract miR-211-mediated CXCR4 suppression and modulate chondrocyte proliferation and apoptosis through downstream mitogen-activated protein kinase (MAPK) signaling [24]. These studies have demonstrated the importance of lncRNAs as ceRNAs in OA. However, the identification of genome-wide expression of lncRNA associated with ceRNA-ceRNA interactions in OA remains unknown.

Therefore, this study aimed to identify and investigate the functional role of biomarkers associated with lncRNA in the progression of OA of the knee by lncRNA-associated ceRNA integrated network analysis.

Material and Methods

Data collection

Human microRNA- long noncoding RNA (miRNA-lncRNA) and miRNA-mRNA target data were collected from StarBase [25,26]. High-quality miRNA-lncRNA and miRNA-mRNA interactions were validated by experiments performed in triplicate that included crosslinking immunoprecipitation (CLIP) and high-throughput sequencing. The miRNA-mRNA interactions were predicted by more than two computing methods as the final dataset. In total, 18,482 and 69,7637 interactions were identified for miRNA-lncRNA and miRNA-mRNA, respectively. Also, lncRNA and mRNA expression profiles for patients with osteoarthritis (OA) of the knee and with mild and severe pain were obtained from the Gene Expression Omnibus (GEO) database (accession GSE99662) [1].

Construction of the integrated network analysis of lncRNA-associated ceRNA

A three-step computational method was used to integrate miRNA-lncRNA/mRNA interactions and lncRNA/mRNA expression profiles in patients with OA with mild and severe pain.

The lncRNA and mRNA pairs were identified that were regulated by same miRNAs. The hypergeometric test was used to determine whether a given lncRNA and mRNA pair significantly shared miRNAs. The Bonferroni correction was performed to control the multiple hypotheses. The lncRNA and mRNA pairs with a false discovery rate (FDR) <0.01 were identified as candidate ceRNA-ceRNA interactions.

Pearson’s correlation coefficients between each lncRNA-mRNA pair in the context of mild pain and severe pain due to OA were identified. Then, the Pearson’s correlation coefficients (PCCs) were defined to measure the lncRNA-mRNA correlation difference between patients with OA of the knee associated with mild pain and severe pain as dPCC=PCCsevere-pain–PCCmild-pain. The lncRNA-mRNA pairs with dPCC >0.5 and significantly positive PCCsevere-pain were considered as gain ceRNA interactions, while lncRNA-mRNA pairs with dPCC <−0.5 and significantly positive PCCmild-pain were considered as loss ceRNA interactions. Assembling all the gain and loss of the ceRNA-ceRNA interactions, a lncRNA-associated ceRNA network in OA was constructed (Supplementary Table 1, available on request from authors). Each node represented a lncRNA or a mRNA, and two nodes were connected if they were dysregulated in OA associated with severe pain.

Network visualization

The networks were visualized by Cytoscape version 3.3.0 [27], including the lncRNA-associated ceRNA network in OA and the differential expression associated subnetwork.

Identification of differential expression of mRNA and lncRNA

Student’s t-test was used to evaluate the significance of differential expression of RNAs between samples of OA associated with mild pain and severe pain. Considering the small sample sizes, the RNAs were selected with P <0.2 as the differentially expressed RNAs [28,29]. Also, the edges were extracted that were related to differentially expressed RNAs and the differential expression subnetwork was constructed of the lncRNA-associated ceRNA network in OA (Supplementary Table 2, available on request from authors).

Results

Construction and characterization of the long noncoding RNA (lncRNA)-associated competing endogenous RNA (ceRNA) network in osteoarthritis (OA)

From the matched lncRNA and mRNA expression profiles for patients with osteoarthritis (OA) and the experimentally validated miRNA-lncRNA/mRNA target interactions, the lncRNA-associated ceRNA network was constructed in OA. The steps taken to construct the network were as described in the Methods. Following the hypergeometric test and Bonferroni correction test to identify the co-regulated lncRNA-mRNA pairs, 40,298 target relationships between 278 lncRNAs and 4,138 mRNAs were identified.

Each pair of lncRNA-mRNAs was found to share at least four miRNAs. The Pearsons’ correlation coefficients (PCCs) were calculated and the significance of each lncRNA-mRNA pair in patients with OA of the knee with mild pain and severe pain was determined. The lncRNA-mRNA pairs with R>0 and P<0.05 were considered as the lncRNA-associated ceRNA network in OA of the knee. The dPCC was defined to measure the different correlation between each lncRNA-mRNA pair in patients with OA of the knee associated with mild pain and severe pain. Combining the dPCC and the p-value using Pearson’s correlation of significance, there were 1,870 dysregulated lncRNA-mRNA interactions identified and the network was constructed of 131 lncRNAs and 1,251 mRNAs (Figure 1A, 1B).

Figure 1

Dysregulated long noncoding RNA (lncRNA)-associated competing endogenous RNA (ceRNA) network in osteoarthritis (OA) of the knee. (A) Global view of the lncRNA-associated ceRNA network in osteoarthritis. The purple and orange nodes represent lncRNA and mRNA respectively. Red and green edges represent gain and loss interactions in samples from patients with severe osteoarthritis (OA) of the knee compared with samples from patients with mild OA of the knee. (B) Statistics of the nodes and edges of the lncRNA-associated ceRNA network in OA of the knee. (C) The degree of distribution of the lncRNA-associated ceRNA network in OA of the knee.

The lncRNA-associated ceRNA network in OA of the knee included a quarter of the gain ceRNA-ceRNA interactions and three-quarters of the loss of ceRNA-ceRNA interactions (Figure 1B). This result suggested that there were relatively more ceRNA-ceRNA that lost their consistent expression patterns in samples from patients with severe pain associated with OA of the knee. To explore the structure and features of the ceRNA-ceRNA interaction network in OA, the degree of distribution analysis of the nodes in the network was performed. As a result, most nodes had relatively few interactions while a small proportion of nodes had several interactions, which fits with a power law distribution, indicating that the network was scale-free and different from randomly generated networks (Figure 1C).

The lncRNA had central roles in the integrated network of lncRNA-associated ceRNA in OA

To characterize the roles of lncRNAs and mRNAs in the ceRNA-ceRNA interaction network in OA of the knee, the network degree of each lncRNA and mRNA were calculated. The lncRNAs were found to have more dysregulated neighbors than mRNAs (P<2.2e-16) (Figure 2A). One lncRNA was linked to a maximum of 94 mRNAs while one mRNA was dysregulated with a maximum of seven lncRNAs. This result indicated that lncRNAs tended to be hub nodes and lncRNAs had centralized roles in the integrated network of lncRNA-associated ceRNA in OA. The categories of lncRNAs in the dysregulated network were analyzed and the results showed that a variety of lncRNAs tended to be dysregulated in OA, especially antisense lncRNA and lncRNA (Figure 2B).

Figure 2

Characteristics of the long noncoding RNA (lncRNA) and mRNAs in the competing endogenous RNA (ceRNA) network interaction in osteoarthritis (OA). (A) The degree of lncRNAs and mRNAs. (B) The category of lncRNAs in the ceRNA-ceRNA network interaction in osteoarthritis (OA). (C) The proportion of gain and loss interactions for each lncRNA in the ceRNA-ceRNA network interaction in OA. (D) The proportion of gain and loss interactions for each mRNA in the ceRNA-ceRNA network interaction in OA.

The gain and loss interactions for each lncRNA were analyzed. For each lncRNA, their proportion of gain and loss interactions were analyzed in all interactions with mRNAs. There were 60% lncRNAs that lost their interactions and 18% lncRNAs had mostly gain interactions with associated mRNAs in patients with OA of the knee that was associated with severe pain (Figure 2C). Also, the gain and loss interactions for each mRNA were calculated. As a result, 69% and 24% of mRNAs were found to be related to loss and gain interactions with lncRNAs in patients with OA of the knee and severe pain (Figure 2D). Only 22% of lncRNAs and 7% of mRNAs showed nondirectional interactions with their neighbors (30–70% gain and 70–30% loss (Figures 2C, 2D).

The differential expression pattern of lncRNAs in the dysregulated ceRNA network

Previous studies have shown that dysregulated ceRNAs contributed to the abnormal transcriptome [30]. Therefore, the expression pattern of lncRNAs and mRNAs in the ceRNA-ceRNA interaction network in OA of the knee were investigated. The upregulated and down-regulated lncRNAs and mRNAs were extracted in samples from patients with OA of the knee associated with severe pain and compared with samples from patients with OA associated with mild pain and identified 16 and 180 differential expression lncRNAs and mRNAs, respectively. Also, the edges were extracted that linked the differentially expressed RNAs and we constructed the differential expression subnetwork with 38 edges, covering eight lncRNAs and 35 mRNAs (Figure 3A). Next, we analyzed the differential expression direction of dysregulated lncRNA-associated ceRNA pairs. As a result, 89% of the gain of ceRNA pairs showed the same dysregulated expression direction, while only 31% of the loss of ceRNA pairs had the same dysregulated expression direction (Figure 3B). This finding was consistent with previous studies [31], and suggested that ceRNA pairs with enhanced co-regulation in samples from patients with OA of the knee associated with severe pain had a consistent expression pattern.

Figure 3

The eight long noncoding RNAs (lncRNAs) identified as possible biomarkers of osteoarthritis (OA). (A) The differential expression subnetwork of the competing endogenous RNA (ceRNA) network in osteoarthritis. The purple and orange nodes represent lncRNAs and mRNAs, respectively. Red and green edges represent gain and loss interactions. (B) The differential expression pattern of ceRNA pairs in the subnetwork. (C) The unsupervised hierarchical clustering heatmap of samples from patients with OA of the knee associated with mild pain and severe pain based on the expression profile of eight lncRNA biomarkers.

To test whether the eight lncRNAs with a dysregulated expression pattern could distinguish between samples from patients with OA associated with severe pain from samples from patients with OA associated with mild pain, the unsupervised hierarchical clustering for five OA with mild pain and five OA with severe pain with the expression profiles of the eight lncRNAs were performed. The result showed that samples from patients with OA associated with mild pain and samples from patients with OA associated with severe pain were separately classified into two clusters (Figure 3C), suggesting that the eight dysregulated lncRNAs could be biomarkers in prediction of pain intensity in patients with OA of the knee.

SNHG16 and LINC00473 affected OA by dysregulated ceRNA interactions

In the differential expression subnetwork, eight dysregulated lncRNAs were identified as candidate biomarkers for OA of the knee, including AC124045.1, AL391244.1, AC106038.1, LINC00473, AL357140.2, SNHG16, AF111167.2, and AC116667.1. Zhu et al. previously showed that SNHG16 was highly expressed in osteosarcoma [32]. In this study, SNHG16 was upregulated in samples from patients with OA of the knee associated with severe pain, suggesting that SNHG16 could promote the severity of OA. Also, the findings suggested that not only the abnormal expression pattern of SNHG16 but also the loss of associated ceRNA interactions (SNHG16 and POLD3, SNHG16 and DNAJA2) could accelerate the severity of OA of the knee (Figure 4A, 4B). However, Zhang et al. found that ZBTB7A could inhibit cisplatin-included apoptosis by repressing LINC00473 expression [33]. We found that LINC00473 was upregulated in patients with severe pain associated with OA and the associated ceRNA interactions were enhanced (Figure 4C, 4D), which showed that LINC00473 and the associated dysregulated ceRNAs might be potential biomarkers for OA.

Figure 4

The differences in co-expression patterns of long noncoding RNAs (lncRNAs) as biomarkers associated with competing endogenous RNA (ceRNA) pairs. The co-expression patterns of SNHG16 and POLD3/DNAJA2 in patients with osteoarthritis (OA) of the knee associated with mild pain (A) and severe pain (B), and the co-expression patterns of LINC00473 and PPRC1/DIEXF in patients with osteoarthritis (OA) of the knee associated with mild pain (C) and severe pain (D).

Discussion

In this study, a long noncoding RNA (lncRNA)-associated competing endogenous RNA (ceRNA) integrated network analysis was performed in patients with osteoarthritis (OA) of the knee. Gain and loss interactions were obtained. Gain interactions represented the enhanced ceRNA-ceRNA interactions in samples from patients with OA associated with severe pain compared with samples from patients with OA associated with mild pain, while the loss of interactions represented the weakened ceRNA-ceRNA interactions. Analysis of the dysregulated ceRNA network showed that lncRNAs played central roles in the network and most of the lncRNAs or mRNAs were related to one type of dysregulated interaction, gain or loss. Integrating the differential expression analysis of lncRNAs identified eight lncRNA biomarkers. Further analysis showed the efficient prediction of the pain level in OA patients based on the expression profiles of eight lncRNAs with an abnormal expression pattern, indicating that these lncRNAs might have potential as biomarkers for OA. The findings of this study showed novel candidate biomarkers associated with the progression of OA of the knee and identified their specific functional mechanism in promoting OA as dysregulated ceRNAs.

Most previous molecular studies have focused on the aberrant expression pattern of single RNA molecules [1,34]. In the present study, we investigated the difference between samples from patients with OA of the knee associated with mild pain and severe pain by evaluating the competitive interactions between ceRNAs instead of single lncRNAs or mRNAs. Only 12% of lncRNAs and 14% of mRNAs in the ceRNA-ceRNA interaction network in OA of the knee showed differential expression. Also, only 2% of the edges in the ceRNA-ceRNA interaction network in OA contained two differentially expressed RNAs. These results suggested that not only abnormal expression but also the abnormal interaction with other RNAs could influence the progression of OA [31]. Calculating the percentage of gain and loss interactions for each lncRNA and mRNA showed that most lncRNAs or mRNAs tended to enhance or weaken the interactions with their neighbors. The dysregulated interactions between the lncRNA and mRNA could be caused by the abnormal expression of the related miRNAs, the abnormal regulation of miRNAs to the lncRNA or mRNA, and the abnormal length of the 3′UTR of lncRNA or mRNA [35,36]. Further studies are needed to explain these findings and to understand the molecular mechanism of the dysregulation of lncRNA-associated ceRNA interactions in OA of the knee.

Conclusions

By integrating the microRNA (miRNA), long noncoding RNA (lncRNA), and mRNA target interactions and expression of lncRNAs and mRNAs in patients with osteoarthritis (OA) of the knee associated with mild pain and severe pain, it was possible to construct a lncRNA-associated competing endogenous RNA (ceRNA) integrated network. Analysis of the network characteristics identified eight lncRNA biomarkers. The findings support the need for further studies on the potential role of lncRNA biomarkers in OA.