Expression of long non‑coding RNAs in chondrocytes from proximal interphalangeal joints.

Osteoarthritis (OA) of hand is a common disease, resulting in disability of the hands. The pathogenesis of hand (H) OA remains to be elucidated, and findings from knee and hip joints cannot be simply applied to HOA. To improve knowledge on the specific biology and pathobiology of HOA, the present study performed bioinformatics analyses to analyze the long non‑coding (lnc) RNA expression profile in human chondrocytes of proximal interphalangeal (PIP) finger joints and knee joints. Gene expression data were downloaded from the Gene Expression Omnibus database, and PIP and knee chondrocytes were analyzed (n=3/group). Probes of the Affymetrix Human Gene 2.0 ST Microarray were annotated to obtain information about lncRNA expression profile. Compared with chondrocytes from knee joints, chondrocytes derived from PIP joints had significantly different lncRNA expression profiles, and 1,172 lncRNAs were differentially expressed. Compared with chondrocyte from knee joints, 534 lncRNAs were upregulated and 638 lncRNAs were downregulated in chondrocytes from PIP joints. A co‑expression network was constructed to analyze the correlation between lncRNAs and protein‑coding genes. Function annotation analyses suggested that protein‑coding genes that are co‑expressed with lncRNAs are enriched in the biological processes of bone morphogenesis, bone development and cartilage development. In conclusion, the present study demonstrated that chondrocytes derived from PIP joints exhibit a significant difference in lncRNA expression compared with chondrocytes derived from knee joints.

Introduction

Osteoarthritis (OA) is a comment disease that affects many major joints including the hip, spine, knee and the hand (1,2). OA leads to severe functional and emotional burdens to patients and causes heavy economic burdens (3). Hand osteoarthritis (HOA) is one of the most common kinds of OA and HOA is a leading cause of disability of the hands. The incidence of HOA is ~20% in people aged 65 or older (4).

HOA is considered as a primarily degenerative process with inflammation (5), and degradation of articular cartilage is the hallmark of HOA. The cartilage degradation process of the knee and hip joints has been well investigated using animal models or human chondrocyte culture. However, the molecular mechanism of HOA is largely unknown. Limited by different anatomical shapes and mechanical requirements, findings from knee or hip joints cannot be simply applied to other joints. For HOA, Stradner et al (6) successfully established a chondrocyte culture of proximal interphalangeal (PIP) cartilage, and demonstrated that chondrocytes from PIP and knee joints have distinct expression patterns of protein-coding genes.

Previous evidence has indicated that non-coding RNAs serve important roles in various biological processes including carcinogenesis, cell differentiation, metabolism, and immunity responses (7–12). Long non-coding RNA (lncRNA) is a kind of non-coding RNA transcript, which is >200 nt long and does not encode proteins (7). Researchers have demonstrated that lncRNAs are involved in regulating the expression levels of various factors involved in the pathological process of OA, including maternally expressed 3 (13), imprinted maternally expressed transcript (14) and HOX transcript antisense RNA (15). However, the lncRNA expression profile from hip and knee joints cannot be simply applied to HOA, and the expression pattern of lncRNA in HOA remains unknown. Therefore, the present study performed a bioinformatics analysis to systematically analyze lncRNA expression profile in the human cartilage of the PIP finger joints by annotation of microarray data to improve the understanding about the pathology of HOA.

Materials and methods

Dataset selection

The gene expression data used in the current study was obtained from the publicly available Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). This GSE68038 dataset was identified and analyzed, which contains 3 samples of knee chondrocytes and 3 samples of PIP chondrocytes. The GSE68038 dataset was based on the Affymetrix Human Gene 2.0 ST Array [transcript (gene) version] (Affymetrix; Thermo Fisher Scientific, Inc.). The raw CEL files were downloaded and quantile normalized and background adjusted using Robust Multichip Average (Windows version) software v1.20.0 (16). After normalization, the expression value of each probe was obtained. The normalized data were then analyzed with Significant Analysis of Microarray software v4.01 (17). Unpaired t-test was used to calculate differentially expressed genes and P<0.05 was considered to indicate a statistically significant difference. Differentially expressed lncRNAs were clustered and visualized by Cluster 3.0 software (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm).

Annotation of lncRNA probes

To identify probe sets that target lncRNA transcripts, the BioMart data portal (www.biomart.org/) was used to download annotations of lncRNA probes. Following this, the probe sets were filtered according to Ensemble transcript type, and transcripts with ‘antisense’, ‘processed_transcripts’, ‘lincRNA’, ‘non_sense_mediated_decay’, ‘sense intronic’ and ‘sense overlapping’ were selected as lncRNA (18,19). After creation of annotation file, the differentially probe sets were matched and annotated.

Co-expression network construction

The co-expression network of lncRNA protein-coding genes was built to identify the interactions between coding genes and lncRNA (20). Differentially expressed lncRNAs and protein-coding genes with P<0.01 were used to construct the co-expression network, and the normalized expression values were retrieved. For each pair of coding gene-lncRNA, coding-coding genes or lncRNA-lncRNA, the P value of Pearson correlation and significant correlation pairs (P<0.05) were used to construct the network (21). Cytoscape software v3.3.0 was utilized to visualize the co-expression network (www.cytoscape.org/).

Functional enrichment analysis

The gene functional enrichment analysis was performed by the Database for Annotation, Visualization and Integrated Discovery online tools (http://david.ncifcrf.gov/) (22). P<0.05 and false discovery rate<0.05 was considered as significant enrichment.

Results

Gene expression profiles of chondrocytes were analyzed with Affymetrix Human Gene 2.0 ST Arrays. After annotation, ~18,000 probes mapped to lncRNAs were identified. According to the statistical threshold, 1,172 differentially expressed lncRNAs were identified. Compared with chondrocytes from knee joints, 534 lncRNAs were upregulated and 638 lncRNAs were downregulated in chondrocytes from PIP joints. The expression patterns of lncRNAs were different between chondrocytes from knee and PIP joints, and the top 100 differentially expressed lncRNAs are presented as a heat map (Fig. 1). As presented, a set of lncRNAs were significantly differentially expressed between PIP and knee joints. The detailed information of 20 most upregulated and downregulated lncRNAs are presented in Table I. Compared with knee joints chondrocyte, ENST00000451530 was mostly upregulated in chondrocytes from PIP, and ENST00000462445 was mostly downregulated.

Figure 1.

Heat map of differentially expressed lncRNAs between chondrocytes from PIP and knee joints. Each row represents an lncRNA and each column represents a sample, the left 3 columns are chondrocytes from knee joints and the right 3 columns are chondrocytes from PIP joints; red, upregulated lncRNAs, green, downregulated lncRNAs. lncRNA, long non-coding RNA; PIP, proximal interphalangeal.

Table I.

Differentially expressed lncRNAs between chondrocytes from PIP and knee joints.

Ensembl gene ID	Ensembl transcript ID	Regulation (PIP vs. knee)	P-value	logFC	FC	Strand	Chromosome name	Transcript start (bp)	Transcript end (bp)	Transcript type
ENSG00000238258	ENST00000451530	Up	0.000105	3.55224	11.73089	1	10	33211277	33213804	Antisense
ENSG00000253357	ENST00000517346	Up	0.00013	8.75114	430.8793	1	5	1.68E+08	1.68E+08	Sense intronic
ENSG00000260758	ENST00000562495	Up	0.000197	1.14294	2.208306	−1	15	86078743	86079792	lincRNA
ENSG00000226677	ENST00000554337	Up	0.000212	1.05432	2.076739	1	14	34939324	34940332	Processed_pseudogene
ENSG00000125430	ENST00000466596	Up	0.000262	8.60867	390.3623	1	17	14301083	14349144	Nonsense_mediated_decay
ENSG00000264425	ENST00000585107	Up	0.000281	2.49343	5.631152	1	7	1.01E+08	1.01E+08	miRNA
ENSG00000231019	ENST00000430214	Up	0.000307	0.59894	1.514603	1	13	88142867	88236082	lincRNA
ENSG00000222276	ENST00000410344	Up	0.000541	2.09069	4.259517	1	14	96384624	96384815	snRNA
ENSG00000232677	ENST00000590657	Up	0.000544	0.82319	1.769314	−1	19	36313067	36315737	lincRNA
ENSG00000173267	ENST00000465679	Up	0.000715	1.04375	2.061579	1	10	86958618	86962873	Processed_transcript
ENSG00000257175	ENST00000548656	Up	0.000781	1.4772	2.784079	−1	14	18634955	18637208	Processed_pseudogene
ENSG00000256879	ENST00000535755	Up	0.000794	2.73779	6.670477	−1	12	20361732	20370262	Antisense
ENSG00000246876	ENST00000509105	Up	0.000804	0.8389	1.788686	−1	4	1.3E+08	1.3E+08	lincRNA
ENSG00000258285	ENST00000547006	Up	0.000828	1.19266	2.285738	1	12	1.17E+08	1.17E+08	lincRNA
ENSG00000233639	ENST00000413121	Up	0.00085	0.68007	1.602217	−1	2	1.05E+08	1.05E+08	lincRNA
ENSG00000183760	ENST00000601575	Up	0.001074	0.6835	1.606031	1	19	39083913	39111493	Nonsense_mediated_decay
ENSG00000226277	ENST00000449119	Up	0.001085	1.02464	2.034452	1	2	421057	422156	lincRNA
ENSG00000104964	ENST00000592414	Up	0.001142	0.65003	1.569201	−1	19	3052910	3056768	Retained_intron
ENSG00000249115	ENST00000587439	Up	0.001171	0.52637	1.440301	1	19	35612744	35623609	Nonsense_mediated_decay
ENSG00000129038	ENST00000566011	Up	0.001176	2.30828	4.952922	1	15	73925989	73951919	Nonsense_mediated_decay
ENSG00000122862	ENST00000462445	Down	3.26E-05	−3.2672	0.103866	1	10	69088106	69104541	Processed_transcript
ENSG00000207036	ENST00000384309	Down	0.000144	−2.3052	0.202333	1	8	70357347	70357448	misc_RNA
ENSG00000212951	ENST00000473402	Down	0.000315	−0.82787	0.56336	−1	9	62532397	62532853	Unprocessed_pseudogene
ENSG00000265520	ENST00000584165	Down	0.000349	−1.70144	0.307479	−1	8	17681578	17681657	miRNA
ENSG00000215304	ENST00000562108	Down	0.000371	−2.93247	0.13099	−1	15	32398956	32435233	Processed_transcript
ENSG00000223929	ENST00000441598	Down	0.000492	−0.70198	0.614728	−1	2	60359720	60383036	lincRNA
ENSG00000255366	ENST00000528139	Down	0.000607	−3.02043	0.123242	1	8	47190772	47193262	lincRNA
ENSG00000178162	ENST00000424873	Down	0.000674	−5.04991	0.030187	−1	2	1.3E+08	1.3E+08	Processed_transcript
ENSG00000236166	ENST00000430428	Down	0.000747	−1.67661	0.312817	1	6	1.32E+08	1.32E+08	lincRNA
ENSG00000226853	ENST00000444196	Down	0.000786	−1.35524	0.39087	1	2	1.74E+08	1.74E+08	lincRNA
ENSG00000207009	ENST00000384282	Down	0.000931	−1.31169	0.402849	−1	4	1683420	1683529	Misc_RNA
ENSG00000188257	ENST00000469162	Down	0.001124	−8.2513	0.003282	−1	1	19975431	19979607	Processed_transcript
ENSG00000262117	ENST00000571259	Down	0.001136	−0.75445	0.592772	−1	16	11819829	11828828	lincRNA
ENSG00000273259	ENST00000553947	Down	0.001138	−7.74996	0.004645	1	14	94592058	94624646	Nonsense_mediated_decay
ENSG00000255480	ENST00000529078	Down	0.001249	−3.04775	0.12093	1	11	30425552	30429268	Antisense
ENSG00000243478	ENST00000487742	Down	0.001249	−0.88687	0.540786	1	2	2.01E+08	2.01E+08	Unitary_pseudogene
ENSG00000258827	ENST00000554595	Down	0.001278	−0.43912	0.737584	−1	14	37097062	37098563	Sense_intronic
ENSG00000273172	ENST00000576171	Down	0.001315	−0.49803	0.708073	−1	17	183824	191587	lincRNA
ENSG00000211836	ENST00000390484	Down	0.001331	−0.4889	0.712568	1	14	22482287	22482346	TR_J_gene
ENSG00000275890	ENST00000616723	Down	0.001356	−5.29061	0.025549	1	GL000205.2	9151	9226	miRNA

lincRNA, long intervening non-coding RNA; miRNA, microRNA; PIP, proximal interphalangeal; lncRNA, long non-coding RNA; FC, fold change; snRNA, small nuclear RNA.

Co-expression networks of lncRNA and protein-coding genes may provide information to infer the biological function of lncRNA, as genes in the same signaling pathway or with the same function may have similar expression patterns. Additionally, co-expression networks have been widely used to predict the potential biological function of lncRNAs. Therefore, a co-expression network was built between lncRNA and protein-coding genes. As presented in Fig. 2A, lncRNAs and protein-coding genes are closely connected. In the co-expression network, RP11-6E9.4 (ENST00000508955) and RP11-713P17.3 (ENST00000529070) were 2 lncRNAs that co-expressed with >30 protein-coding genes. In addition, many inflammation-associated genes were co-expressed with these lncRNAs, including interleukin (IL)-7R, IL-19 and chemokine (C-C) motif ligand 1. On the other hand, several Wnt pathway genes, including Wnt5A, were also co-expressed with these lncRNAs.

Figure 2.

(A) Co-expression network of lncRNAs and protein-coding genes. Yellow dots, lncRNAs; green dots, protein-coding genes; light red dots, centers of the sub-networks. (B) Enrichment plot of gene ontology results. Biological processes associated with bone and cartilage are highlighted with red, and grey dots are terms not associated with bone or cartilage. (C) Column chart of gene ontology results, and bone and cartilage-associated gene ontology terms are highlighted with red boxes. lncRNA, long non-coding RNA.

Co-expression networks are a potential method to predict lncRNA function, as genes involved in the same biological processes are usually co-expressed. Therefore, a functional enrichment analyses for protein-coding genes in the co-expression network was performed to further predict their biological functions. As presented in Fig. 2C and D, gene ontology analyses demonstrated that these protein-coding genes were significantly associated with bone morphogenesis, bone development, skeletal system development and cartilage development, suggesting these protein-coding genes and lncRNAs may be involved in OA pathogenesis.

Discussion

Articular cartilage degradation is a hallmark of OA. Human chondrocyte culture has been extensively used to investigate the cartilage degradation of knee and hip joints. However, little is known is about the cartilage degradation in HOA, and findings from hip and knee joints cannot be simply applied to HOA due to different anatomical shapes and mechanical requirements between joints (23,24).

The present study compared lncRNA expression profiles between chondrocytes from PIP and knee joints. The results demonstrated that chondrocytes from PIP and knee joints have different lncRNA expression patterns. Compared with chondrocytes from knee joints, 534 lncRNAs were upregulated and 638 lncRNAs were downregulated in chondrocytes from PIP joints. The different expression profile of lncRNA between chondrocyte from PIP and knee joints supported previous reports that lncRNA expression is highly temporally and specially specific (7). Therefore, these differentially expressed lncRNAs in PIP joints may be involved in the biological function of chondrocytes, and the pathological process of HOA.

Previous research has investigated the lncRNA expression profile in patients with rheumatoid arthritis (25,26), and it was concluded that lncRNA may contribute to the pathogenesis of rheumatoid arthritis (26). Compared with these studies in rheumatoid arthritis, lncRNA expression in HOA was analyzed with bioinformatics analyses, and the potential association between lncRNAs and coding genes was inferred. However, HOA and rheumatoid arthritis have different pathologies and clinical features; therefore, findings between the two cannot be easily compared.

Co-expression networks between lncRNA and protein-coding genes were constructed to infer the potential function of these lncRNAs. GO enrichment analyses demonstrated that protein-coding genes in the co-expression network are associated with many skeletal system-specific items, including bone, cartilage and skeletal system development. The results suggested that these lncRNAs were highly likely involved in these biological processes, and these lncRNAs may also serve important roles in the pathological process of HOA.

Subsequently, the present study demonstrated that RP11-6E9.4 (ENST00000508955) and RP11-713P17.3 (ENST00000529070) were co-expressed with the most protein-coding genes. Compared with chondrocytes from knee joints, RP11-6E9.4 and RP11-713P17.3 were significantly downregulated in PIP joints. In addition, we constructed a co-expression network of coding genes and lncRNAs and found RP11-6E9.4 and RP11-713P17.3 were co-expressed with several inflammation and cartilage specific genes, including collagen type I alpha I and transcription factor SOX9 (27). Therefore, RP11-6E9.4 and RP11-713P17.3 are potentially involved in the pathogenesis of OA; however, further experiments are required to validate the functional roles of these lncRNAs.

The biological function and molecular mechanism have been widely investigated in various biological and pathological processes. To the best of our knowledge, there is currently no literature on the expression profile of HOA in lncRNA. However, the present study has several limitations. Firstly, this is only a bioinformatics analysis; thus, these findings should be further validated by laboratory experiments. Secondly, the sample sizes analyzed were small, which might limit the application of these findings.

In conclusion, the present study demonstrated that chondrocytes from PIP joints have different lncRNA expression profiles. These findings may improve knowledge on the biological roles of lncRNA and the pathology of HOA.