Identify CRNDE and LINC00152 as the key lncRNAs in age-related degeneration of articular cartilage through comprehensive and integrative analysis

Introduction

Osteoarthritis (OA) is a common chronic joint disease in elderly individuals, with weight-bearing joints such as the knees, hips, and the spine being the most commonly affected. Between 2011 and 2012, the prevalence of symptomatic knee OA in China was 8.1% (Tang et al., 2016). OA progresses slowly over time, eventually resulting in joint pain, stiffness, and dysfunction (Hunter & Felson, 2006). With the growing number of patients suffering from OA, the treatment costs of OA are expected to present a huge economic burden to both individuals and health systems worldwide. While the exact mechanisms underlying OA remain poorly understood, numerous risk factors have been identified, including age, obesity, sex, genetic predisposition, and joint injury (Loeser, 2017; Silverwood et al., 2015). Among these factors, age is one of the most important critical influences on the degeneration of articular cartilage. Previous studies have demonstrated that the prevalence of radiographic knee OA increases with age (Shane Anderson & Loeser, 2010). Aging-related changes in joint tissues include accumulation of senescent secretory phenotype cells, increased expressions of catabolic factors, a gradual loss of extracellular cartilage matrix, and oxidative stress (Campisi & d’Addadi Fagagna, 2007; Ding et al., 2005; Henrotin, Bruckner & Pujol, 2003; Wu et al., 2002). Moreover, age-related increases in mesenchymal stromal cells in the bone marrow are also thought to regulate of osteoblast formation and bone remodeling (Ganguly et al., 2017).

Long non-coding RNAs (lncRNAs), defined as non-coding capacity RNAs >200 bp in length, play an important role in transcriptional regulation (Ponting, Oliver & Reik, 2009). Recent studies have identified several lncRNA-associated functions involved in the development of human diseases, including several forms of cancer, cardiovascular disease, diabetes, and osteoporosis (Cao et al., 2013; Hao et al., 2015; Huang, 2018; Leti & DiStefano, 2017; Wang et al., 2018b), as well as OA. Using DNA microarrays and bioinformatic approaches, Ming et al. identified a number of lncRNAs that are differentially expressed in OA cartilage compared to normal samples (Fu et al., 2015). Similarly, Xiang et al. found that synovium samples obtained from OA patients showed differential expression of numerous mRNAs and lncRNAs (Xiang et al., 2018). In addition, an increasing number of lncRNAs have been implicated in chondrocytes proliferation and cartilage metabolism. For example, increased expression of lncRNA-ZFAS1 was shown to induce the proliferation and migration of chondrocytes (Ye et al., 2018), while increased expression of lncRNA-HOTAIR contributes to chondrocytes apoptosis and matrix degradation (Hu et al., 2018). Finally, lncRNA-CIR acts as a sponge for miRNA-27b, resulting in increased extracellular matrix degradation via its effects on MMP13 expression (Li et al., 2017c). Taken together, these data suggest that lncRNAs play key roles in the pathogenesis of OA.

Because OA is an age-related chronic disease, morphological and molecular biological characters may appear differently in different stages. Here, to explore differential expression of age-related lncRNAs in each stage of OA, we examined two previously published OA gene datasets available in the Gene Expression Omnibus (GEO) database. Differential expression of lncRNAs across different developmental stages was examined by comprehensive bioinformatic analyses. Then the differentially expressed lncRNAs were correlated with other DEGs and potential diseases based on data obtained from the RAID and MNDR databases. Moreover, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were also performed, and a protein–protein interaction (PPI) network was used to screen for crucial genes and lncRNAs.

Materials and Methods

Expression profile dataset

Gene expression profiles for age-related degeneration in rat knee articular cartilage at different development stages (GSE66554) were obtained from the GEO database and analyzed alongside a second gene expression dataset (GSE113825) which investigated lncRNA and mRNA expression profiles of human knee OA. A detailed workflow of the data analyses is shown in Fig. 1. The quality of gene expression data was analyzed and visualized using the ggplot2 package of R software for each sample.

Results

Data distribution analyses and DEG screening

In GSE113825, a total of 95,722 genes were detected in each sample. The expression values (ranging from 0 to 20 log₂ [fragments per kilobase of transcript, per million fragments sequenced, FPKM]) and the distributions were similar between normal group and OA group (Fig. 2A). Principal component and h-cluster analyses showed that samples were easily grouped into different groups. As for GSE66554, a total of 36685 genes were detected in each sample, with the values and distributions of these genes similar across the five groups (Fig. 2B). Based on these assessments, further bioinformatics analyses could be performed based on the available data.

Following data processing using limma, we identified 4470 DEGs (1658 upregulated, 2812 downregulated) in human normal vs. OA, including 2891 DEGs (1596 upregulated, 1295 downregulated) in rat T0 vs. rat T1, 4354 DEGs (2451 upregulated, 1903 downregulated) in rat T0 vs. rat T2, 5606 DEGs (2703 upregulated, 2903 downregulated) in rat T0 vs. rat T3, and 14063 DEGs (6432 upregulated, 7631 downregulated) in rat T0 vs. rat T4 (Fig. 2C).

Intersection between the predicted targets of lncRNAs and DEGs

Using the GSE66554 rat OA dataset, we identified a total of 1355 DEGs common across each of the different stages (Fig. 3A). Hierarchical clustering of the identified DEGs is displayed as a heatmap in Fig. 3B. Next, homologous human genes were identified for each of the rat DEGs and compared against DEGs identified in the human OA dataset (Fig. 3C), resulting in an overlap of 174 age-related OA genes, including 33 upregulated and 141 downregulated genes (Table S1). Alignment of these to the human genome is shown in Fig. 3D, as intuitively, green means go and red means stop/downregulated.

GO annotation and KEGG pathway enrichment analyses

GO annotation and KEGG pathway enrichment analyses were performed to better understand the functional significance of the DEGs. Biological processes (BP) were significantly enriched for cell adhesion, endodermal cell differentiation, collagen fibril organization, regulation of cell migration, and cell–matrix adhesion (Fig. 4A). Significant cellular component (CC) terms revealed enrichment in extracellular exosome, extracellular space, proteinaceous extracellular matrix, focal adhesion, and sarcolemma. Finally, the top five molecular functions (MF) identified were calcium ion binding, extracellular matrix structural constituent, semaphorin receptor binding, chemorepellent activity, and spectrin binding. All significantly enriched entries are shown in a histogram in Fig. 4B. KEGG analyses of the DEGs revealed enrichment for 18 terms (Fig. 4C). The top five significantly enriched KEGG pathways were the PI3K–Akt signaling pathway, focal adhesion, ECM–receptor interaction, protein digestion and absorption, and pathways in cancer. Detailed information regarding each of the DEGs involved in the BP analyses and KEGG pathway are listed in Fig. 4.

PPI network construction and hub lncRNAs identification

We constructed a putative PPI network map for the overlapping DEGs using the STRING database, which was visualized with Cytoscape. Excluding the DEGs distributed on the edge of the PPI network, the remaining 110 central DEGs were evaluated using the centiscape plugin. After degree calculation, a total of 34 hub genes were colored red in the PPI network (Fig. 5A and Table 2). Using this same approach, we also constructed a lncRNA–gene co-expression network. Each hub lncRNA was determined by calculating its downstream hub genes whose count should be ranked in the top three (minimum of two genes). A total of 23 lncRNAs were screened and are presented in Table 3. From this analysis, the lncRNAs CRNDE (target DEGs: COL4A2, COL6A2, COL6A1, COL5A1, AGRN, POSTN and SGCG) and LINC00152 (target DEG: CDKN1A) stood out. We also identified an important hub gene, PPARGC1A, which exhibited tight connections with other lncRNAs (Fig. 5C).

qRT-PCR validation

To confirm the results of our bioinformatic analysis, we examined the expression of CRNDE and LNC00152 by qRT-PCR in 16 cartilage samples. According to Fig. 6, the transcriptional levels of CRNDE and LNC00152 were significantly decreased in OA group compared with the normal group. These results were consistent with our previous integrative analysis listed in Table 3 and showed the same trends of these hub lncRNAs.

Table 2:

Hub-genes identified by centiscape (degree ≥ 5) among 174 DEGs.

Name	logFC	AveExpr	P.Value	Transcript_id	Degree
ITGB1	−1.75	16.20	5.68E–04	NM_002211	24
COL4A1	−2.08	6.32	1.06E–05	NM_001845	15
DMD	−2.54	5.21	1.16E–04	NM_000109	15
PIK3CG	−2.76	5.00	2.46E–03	NM_002649	15
ACTA2	−3.98	8.71	2.63E–05	NM_001141945	13
COL3A1	−2.17	14.51	9.08E–03	NM_000090	13
COL4A2	−2.39	5.86	1.27E–06	NM_001846	13
COL6A2	−2.16	7.79	5.55E–03	NM_001849	13
COL6A1	−3.24	12.98	2.16E–03	NM_001848	12
ITGB2	−1.57	5.28	2.96E–02	NM_000211	12
COL5A1	−3.98	12.58	4.02E–03	NM_000093	11
LAMB1	−4.11	7.02	6.24E–05	NM_002291	11
NCAM1	−2.70	9.34	8.57E–03	NM_000615	11
AGRN	−1.29	5.82	3.74E–05	NM_198576	10
CD34	−2.22	8.14	8.45E–03	NM_001025109	10
FBN1	−2.32	12.40	2.56E–05	NM_000138	10
COL12A1	−2.64	14.08	4.18E–05	NM_004370	9
COL15A1	−3.08	8.66	7.14E–03	NM_001855	9
IGF1	−2.17	8.53	2.33E–04	NM_000618	9
LAMC1	−1.14	5.81	4.06E–03	NM_002293	9
PLEK	−3.59	7.30	2.46E–04	NM_002664	9
COL14A1	−4.13	12.30	7.96E–04	NM_021110	8
LAMA4	−1.82	7.39	8.45E–05	NM_001105206	8
NID1	−4.19	9.38	2.22E–08	NM_002508	8
RRAD	−1.19	4.96	3.31E–03	NM_001128850	7
VCAN	−4.09	10.96	2.68E–04	NM_001126336	7
CAP2	−1.17	6.67	2.07E–03	NM_006366	6
H2AFX	−1.19	6.96	2.48E–04	NM_002105	6
LOX	−2.88	11.46	1.92E–03	NM_001178102	6
PLEC	−1.05	8.23	4.22E–03	NM_000445	6
CDKN1A	−1.73	7.04	6.11E–03	NM_000389	5
NID2	−3.56	9.41	3.20E–05	NM_007361	5
POSTN	−7.29	10.68	1.97E–03	NM_001135934	5
SGCG	−1.26	4.39	3.62E–02	NM_000231	5

DOI: 10.7717/peerj.7024/table-2

Table 3:

All differentially expressed lncRNAs upstream of DEGs.

Gene symbol	logFC	AveExpr	P.Value	change
ACOX2	−1.29	5.88	1.83E–05	DOWN
ANKAR	1.44	9.79	6.63E–05	UP
BST2	−1.26	6.01	8.68E–04	DOWN
CALD1	−2.03	11.76	8.24E–05	DOWN
CRNDE	−1.62	6.96	5.43E–06	DOWN
DDIT4	−1.08	6.48	6.00E–03	DOWN
DNAJA1	−1.18	10.99	1.41E–04	DOWN
ENO1	−1.24	12.23	3.24E–02	DOWN
EPS8	1.39	10.73	2.24E–04	UP
ERBB3	1.35	6.31	1.90E–03	UP
FASTKD1	−1.27	10.04	2.31E–03	DOWN
HSPA8	−1.20	14.94	7.94E–05	DOWN
ITGAV	−1.35	9.69	6.15E–04	DOWN
KHDRBS3	−2.58	7.15	4.01E–04	DOWN
LINC00152	−2.20	8.75	2.19E–04	DOWN
MYO7A	1.09	6.56	4.67E–03	UP
N4BP2L2-IT2	1.27	11.17	1.64E–04	UP
NDRG2	1.83	9.05	9.92E–03	UP
RAPGEF5	1.88	5.55	2.01E–05	UP
SLC20A1	1.23	11.19	1.77E–02	UP
TMEM161B-AS1	1.02	7.91	1.70E–05	UP
TRIP12	−1.67	10.53	1.84E–04	DOWN
XPNPEP3	1.22	7.23	1.04E–03	UP

DOI: 10.7717/peerj.7024/table-3

Discussion

LncRNAs have emerged as critical modulators of transcriptional, post-transcriptional, and epigenetic gene regulation, with emerging evidence suggesting a relationship between lncRNAs and OA. With the development of high-throughput sequencing technology, an increasing number of RNA-sequencing projects have been performed, identifying key functional mRNAs, miRNAs, lncRNAs, and circRNAs in the development of OA. Xing et al. (2014) reported six lncRNAs, including HOTAIR, GAS5, PMS2L2, RP11–445H22.4, H19, and CTD-2574D22.4, which were differently expressed in OA samples in microarray analyses (Xing et al., 2014). Subsequent studies have since described lncRNAs expression patterns in human osteoarthritic cartilage using a combination of DNA microarray and bioinformatic analyses (Fu et al., 2015). Additional lncRNAs including PACER, CILinc01, and CILinc02, have also been reported in human hip OA chondrocytes, and were significantly associated with the OA inflammatory response (Pearson et al., 2016). Similarly, when comparing the expression of lncRNAs in OA cartilage of variable severity, four hub lncRNAs, SNHG5, ZFAS1, GAS5, and DANCR, were identified as key functional mediators in OA pathogenesis (Xiao et al., 2018).

Although many bioinformatics analyses have revealed differences in gene expression patterns in OA samples, few have examined what, if any, age-related lncRNAs are involved in this process. Therefore, a detailed database of lncRNA observed in rat knee articular cartilage at different developmental stages was used in our study. By comparing differential gene expression in five stages, newborn (T0), youth (T1), adult (T2), early-stage elderly (T3), and later-stage elderly (T4), a set of age-related DEGs were identified. Next, we identified the homologous human genes for each of our age-related rat DEGs, which then were compared against a separate set of DEGs extracted from a human OA lncRNAs expression database. Our results broadly confirmed dysregulated gene in age-related cartilage, including ITGB1, COL4A1, DMD, PIK3CG, ACTA2, COL3A1, COL4A2, COL6A2, and COL6A1. We also identified two hub lncRNAs, CRNDE and LINC00152, which were not identified in previous studies.

LINC00152 is an 828 bp lncRNA located on chromosome 2p11.2. Previous studies have found that LINC00152 plays an oncogenic role in the development of a wide range of tumor types (Chen et al., 2018a; Chen et al., 2018c). Enhanced LINC00152 expression has also been found to be a potential prognostic biomarker in patients with lung and colorectal cancers (Chen et al., 2018b; Li et al., 2017b). LINC00152 also acts as a sponge for various miRNAs, as shown in a recent report describing the interaction between LINC00152 and miR-139 in colorectal cancer cells (Bian et al., 2017). This observation is important, as emerging evidence has shown increased expression of miR-139 in OA cartilage, with increased miR-139 expression resulting in higher expression of IL-6 and chondrocyte apoptosis (Hu et al., 2016; Makki & Haqqi, 2015). Furthermore, KEGG analyses of LINC00152 (Fig. 7A) revealed strong associations between osteoclast differentiation, cell cycle, and Wnt signaling pathways and LINC00152 in OA. Given these observations, the LINC00152/miRNA axis will likely become a major focus of future OA studies.

Another hub lncRNA identified in our study, CRNDE, is a 1910 nt cancer-secreted lncRNA transcribed by human chromosome 16 (16q12.2). This lncRNA has been observed in a variety of cancers, including osteosarcoma, colorectal, cervical, and gastric cancers (Ding et al., 2017; Hu et al., 2017; Li et al., 2018; Meng et al., 2017). Huan et al. (2017) reported an important role for CRNDE in the pathophysiology of human breast cancer. CRNDE was shown to modulate the Wnt/β-catenin signaling pathway by repressing miR-136 expression via the miRNA sponge mechanism, resulting in significant difference in proliferation, migration, and invasion of breast cancer cells. Coincidentally, another group demonstrated that miR-136 plays an important role in the regulation of chondrogenesis in human adipose-derived stem cells (Zhang et al., 2012). Further research revealed that miR-136 bound to the 3′-UTR of MMP13, a key catabolic enzyme involved in the degradation of ECM (Li et al., 2017a). By competitively binding with miR-136, circRNA-CER regulates MMP13 expression, further promoting chondrocyte ECM degradation (Liu et al., 2016). Beyond miR136, decades of research have shown that CRNDE modulates miR-384 activity in a number of tumor cell types (Chen et al., 2016; Zheng et al., 2016). Recently, miR-384–5p was shown to induce OA by regulating the expression of SOX9 and downregulating activity of the NF-κB signaling pathway (Zhang et al., 2018). Taken together, these data suggest that CRNDE can act as a miRNA sponge, competing for miRNA binding with protein-coding transcripts, although the exact mechanism of lncRNA CRNDE-mediated regulation of OA pathogenesis remains to be investigated. Other potential miRNAs related to these two hub lncRNAs have been curated from databases (RAID 2.0, starBASE v2.0, miRTarBase), and are listed in Table 4. Targets suggested by KEGG analyses suggest that the PI3K–Akt signaling pathway, proteoglycan metabolism, and the Ras signaling pathway are promising targets (Fig. 7B).

Outside of lncRNAs, the PPI network identified PPARGC1A as a novel mediator of OA pathogenesis due to its close relationship with various age-related lncRNAs. PPARGC1A, also known as human-accelerated region 20, encodes the peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PGC-1α). PGC-1 α is widely viewed as a master regulator of mitochondrial biogenesis (Liang & Ward, 2006). Reduced expression of PGC-1α has been observed in knee cartilage of aged mice, where it may attenuate oxidative stress and cartilage erosion in OA (Zhao et al., 2014). Activation of PGC-1α expression significantly enhances mitochondrial biogenesis and inhibits oxidative phosphorylation in OA chondrocytes. Similarly, the AMP-activated protein kinase/PGC-1 signaling pathway is an effective target for OA treatment in a rat model (Wang et al., 2018a).

Conclusions

The current study identified a series of DEGs and lncRNAs in each developmental stage of articular cartilage. Using a series of bioinformatics analyses, two pivotal lncRNAs, CRNDE and LINC00152, were identified, which are strongly associated with age-related cartilage degradation. Further functional analyses suggest a potential mechanism through which these hub lncRNAs mediate OA pathogenesis. Our results highlight the important role of lncRNAs in the pathogenesis of OA, although more research is necessary to confirm our findings.

Supplemental Information