A landscape of circulating long non-coding RNA (lncRNA) expression profile and the predictive value of candidate lncRNAs for disease risk of knee osteoarthritis.

BACKGROUND: This study aimed to investigate the plasma long non-coding RNA (lncRNA) expression profile in knee osteoarthritis (KOA) patients and the value of candidate lncRNAs for predicting KOA risk.
METHODS: Plasma was obtained for RNA sequencing (RNA-seq) in eight KOA patients and eight healthy controls (Ctrls). Ten candidate lncRNAs were then selected from the differentially expressed (DE) lncRNAs according to the rank of absolute value of Log2 (fold change). Afterward, RT-qPCR was used to examine 10 candidate lncRNAs expressions in plasma of 100 KOA patients and 100 Ctrls.
RESULTS: In eight KOA patients and eight Ctrls, principal component analysis and heatmap plots disclosed that lncRNA and mRNA expression profile could distinguish KOA patients from Ctrls. Then Volcano plot identified 418 upregulated lncRNAs, 347 downregulated lncRNAs, 521 upregulated mRNAs, and 333 downregulated mRNAs in KOA patients compared to Ctrls. Next, enrichment analyses revealed that DE lncRNAs were mainly enriched in biological processes, molecular functions, and signaling pathways related to inflammation and bone formation. In 100 KOA patients and 100 Ctrls, eight candidate lncRNAs were dysregulated in KOA patients compared to Ctrls, including lncRNA ABCF2P2, lncRNA RP13-16H11.7, lncRNA CTC-340A15.2, lncRNA RP4-735C1.6, lncRNA RP11-293G6-B.8, lncRNA RP11-1246C19.1, lncRNA RP11-303E16.6, and lncRNA RP5-882C2.2. Receiver operating characteristic curve analysis revealed that these eight candidate lncRNAs presented with values for predicting KOA risk. Furthermore, multivariate logistic regression elucidated that six candidate lncRNAs could independently predict KOA risk.
CONCLUSION: We disclosed a landscape of circulating lncRNA expression profile in KOA patients, and discovered several specific lncRNAs which could assist in KOA management.

INTRODUCTION

Knee osteoarthritis (KOA) is one of the worldwide major causes of functional disability in the joints attacking approximately 30% of the older population aged over 60. , , KOA is predominantly a disease of the knee joint degeneration, which is pathogenetically presented as cartilage bone damage accompanied with hyperostosis, cystic change of the cartilage, edema in the synovial tissue and so on. , , , The patients' life quality is seriously injured by KOA with consistent pain from morning to night, knee pain in daily activities, and the most severe condition‐disability, which result in a situation that requires early diagnosis and prompt intervention of KOA. However, the definite diagnosis of KOA demands comprehensive data including clinical symptoms and imaging examinations, besides, no clinical manifestations or biomarkers with satisfactory sensitivity or specificity have been identified until now.

Long non‐coding RNA (lncRNA) is a class of RNAs with almost no protein‐coding ability and is composed of more than 200 nucleotides, which is unique compared to other non‐coding RNAs, such as microRNAs (miRNAs). The number of discovered lncRNAs is growing consistently, and in KOA, there are also reports indicating the roles of lncRNAs. For instance, lnc‐HOTAIR is reported to aggravate KOA by enhancing cell apoptosis of chondrocytes through sponging miR‐130a‐3p. However, there are only limited studies that investigate the lncRNA expression profile in KOA, in which the lncRNA expression profile is assessed in tissue samples (synovial tissue and cartilage tissue) that is not as convenient as the circulating samples regarding clinical utility. ,

Thus, in this study, we aimed to investigate the plasma lncRNA expression profile in KOA patients, and the value of candidate lncRNAs for predicting KOA risk.

METHODS

Participants

From July 2018 to June 2019, a total of 100 KOA patients were consecutively recruited in our hospital. All patients were newly diagnosed as KOA according to the European league against rheumatism (EULAR) evidence‐based recommendations for the diagnosis of knee osteoarthritis (2010). The screening criteria were: (a) age ≥40 years old; (b) not received non‐steroidal anti‐inflammatory drugs (NSAIDS), hyaluronic acid (HA) or immunosuppressor in the past 3 months; (c) without autoimmune disease, cardio‐cerebrovascular disease, hematological malignancies or solid tumors; (d) no pregnancy or lactation. In addition, 100 healthy subjects whose age, gender, and body mass index (BMI) matched with KOA patients were enrolled as control (Ctrl) cohort, between January 2019 and October 2019. This study was approved by the Institutional Review Board of our hospital, and all procedures were conducted according to the Declaration of Helsinki. All participants provided informed consents before enrollment.

Data and sample collection

After enrollment, the demographic characteristics (age, gender, and BMI) of participants were recorded. Peripheral blood (PB) samples of KOA patients were collected before treatment, and PB samples of Ctrl were collected after enrollment. All PB samples were centrifuged at 3000 g for 20 minutes (4°C) to separate plasma. After the plasma samples were isolated, total RNA was extracted from plasma samples using Trizol reagent (Thermo Fisher Scientific) according to the manufacturer's protocol and then stored at −80°C until detection.

Library preparation and RNA sequencing (RNA‐seq)

Eight RNA samples from 100 KOA patients and eight RNA samples from 100 matched Ctrl subjects were selected for RNA‐seq analysis. In brief, the RNA concentration and purity were measured by NanoDrop ND‐1000 spectrophotometer (Thermo), and the RNA integrity was assessed by Agilent 2100 Bioanalyzer (Agilent). After that, the mRNA and non‐coding RNAs were enriched by removing ribosomal RNA from the total RNA using Epicentre Ribo‐zero™ rRNA Removal Kit (Epicentre), then the cDNA library was constructed by NEBNext® Ultra™ Directional RNA Library Prep Kit (NEB) according to the instructions provided by the manufacturer. After the amplification by polymerase chain reaction (PCR) and the quantification by Bioanalyzer 2100 system (Agilent Technologies), the cDNA library of mRNA and lncRNA was sequenced on Illumina Hiseq X10 platform (Illumina), and 150 bp paired‐end reads were obtained. Besides, the RNA‐seq was performed in one replicate.

Workflow of lncRNAs and mRNAs sequencing data

Quality assessment (pre‐trimming), adapter removal, trimming, and quality assessment (post‐trimming) of raw reads were performed as described in a previous study. Afterward, the trimmed reads were aligned to the reference genome to the human reference genome (hg38) using Tophat (http://tophat.cbcb.umd.edu/) with the default parameters. The reads count of lncRNA and mRNA was then calculated using featureCounts based on the annotation file (Homo_sapiens.GRCh38.83.gtf) in Ensembl database. The counts normalization and differential expression analysis were performed using DeSeq2 (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html). The DE lncRNAs and mRNAs were defined as the lncRNAs or mRNAs with a Log2 (fold change [FC]) ≥1.0 and an adjusted P value (Benjamini and Hochberg's BH multiple test correction) <.05.

Bioinformatic analysis

The bioinformatics analyses were mainly performed using R software (Version 3.3.3). Principal component analysis (PCA) and heatmap plot were created to obtain an overview of the expression profile of lncRNA and mRNA using factoextra and pheatmap packages. Volcano plots were constructed to display DE lncRNAs and DE mRNAs using ggplot2 package. The Fisher's exact test was used to assess statistical significance for DE lncRNAs and DE mRNAs enrichment analysis, including GO (Gene ontology: http://www.geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomes: http://www.genome.jp/kegg/). The target genes of DE lncRNAs were predicted by calculating the Pearson correlation coefficients and P values among multiple genes, and the |correlation| ≥ 0.9 were used to filter the transcripts. Regulatory network of top 20 DE lncRNAs including top 10 upregulated lncRNAs and top 10 downregulated lncRNAs by trans targets (mRNAs) was plotted using igraph package.

Reverse transcription–quantitative PCR (RT‐qPCR) validation

For RT‐qPCR validation, 10 lncRNAs including top five upregulated lncRNAs and top five downregulated lncRNAs were screened out according to the rank of Log2 fold change (FC) absolute value from DE lncRNAs identified in the RNA‐seq analysis as candidate lncRNAs. And the expression of these 10 candidate lncRNAs was further verified in 100 KOA patients and 100 Ctrl subjects by RT‐qPCR. The RT‐qPCR was conducted as follows: first, we used TRIzol™ Reagent (Thermo Fisher Scientific) to extract total RNA from plasma, then the total RNA was reversely transcribed by iScript™ cDNA Synthesis Kit (Bio‐Rad); second, THUNDERBIRD® SYBR® qPCR Mix (Toyobo) was applied for PCR procedure; third, the relative expression of 10 lncRNAs was measured using GADPH as internal reference and was calculated by 2−ΔΔct formula. The primers were shown in Table S1. Besides, the RT‐qPCR was performed in triplicates.

Statistical analysis

Statistical analyses were performed using IBM SPSS version 22.0 (IBM), and figures were plotted using GraphPad Prism version 7.00 (GraphPad Software). Comparison of demographic characteristics between KOA patients and Ctrl subjects was determined by Student's t test or chi‐square test. Comparison of 10 candidate lncRNAs between KOA patients and Ctrl subjects was determined by Wilcoxon rank‐sum test. The performances of 10 candidate lncRNAs in predicting KOA risk were illuminated by receiver operating characteristic (ROC) curve and area under the curve (AUC) with 95% confidence interval (CI). Predicting values of 10 candidate lncRNAs for KOA risk were analyzed by univariate and back backward stepwise multivariate logistic regression model. P value < .05 was considered significant.

RESULTS

Demographic characteristics of KOA patients and Ctrls in RNA‐seq

In eight KOA patients and eight Ctrls for the RNA‐seq and bioinformatics analyses, the mean age of KOA patients was 50.8 ± 6.3 and 51.9 ± 6.6 years in Ctrls (P = .732) (Table 1). The number of females and males was both 4 (50.0%) in KOA patients and in Ctrls, respectively (P = 1.000). In addition, the mean value of BMI was 24.4 ± 1.3 and 24.1 ± 1.0 kg/m2 in KOA patients and Ctrls, respectively (P = .660).

TABLE 1

Demographic characteristics

Items	RNA‐Seq			RT‐qPCR
	KOA (N = 8)	Ctrl (N = 8)	P value	KOA (N = 100)	Ctrl (N = 100)	P value
Age (y), mean ± SD	50.8 ± 6.3	51.9 ± 6.6	.732	52.9 ± 6.2	51.6 ± 8.4	.228
Gender, No. (%)
Female	4 (50.0)	4 (50.0)	1.000	73 (73.0)	67 (67.0)	.355
Male	4 (50.0)	4 (50.0)		27 (27.0)	33 (33.0)
BMI (kg/m2), mean ± SD	24.4 ± 1.3	24.1 ± 1.0	.660	23.2 ± 2.3	22.9 ± 2.5	.291

Comparison was determined by Student's t test or chi‐square test.

Abbreviations: BMI, body mass index; KOA, knee osteoarthritis; RNA‐Seq, RNA sequencing; RT‐qPCR, reverse transcription–quantitative polymerase chain reaction; SD, standard deviation.

LncRNA and mRNA expression pattern assessed by PCA and heatmap plot

In eight KOA patients and eight Ctrls, according to PCA pot analysis, lncRNA expression profile (Figure 1A) and mRNA expression profile (Figure 1B) in plasma both clearly distinguished KOA patients from Ctrls. Furthermore, the heatmap plot also disclosed that KOA patients and Ctrls could be differentiated by lncRNA expression profile (Figure 2A) as well as mRNA expression profile (Figure 2B).

FIGURE 1

PCA plot. The lncRNA (A) and mRNA (B) expression profile by PCA plot. Ctrl, control; KOA, knee osteoarthritis; lncRNA, long non‐coding RNA; PCA, principal component analysis

FIGURE 2

Heatmap plot. The lncRNA and mRNA expression profile by heatmap plot. Ctrl, control; KOA, knee osteoarthritis; lncRNA, long non‐coding RNA

DE lncRNAs and DE mRNAs assessed by Volcano plot

Furthermore, the volcano plot showed that there were 418 upregulated lncRNAs and 347 downregulated lncRNAs in KOA patients compared with Ctrls (Figure 3A). As for DE mRNAs, the volcano plot identified 521 upregulated mRNAs and 333 downregulated mRNAs in KOA patients compared to Ctrls (Figure 3B).

FIGURE 3

Volcano plot. The upregulated/unchanged/downregulated lncRNAs (A), and the upregulated/unchanged/downregulated mRNAs by Volcano plot (B). DEG, differentially expressed gene; DOWN, downregulated lncRNAs; lncRNA, long non‐coding RNA; UP, upregulated lncRNAs

Enrichment analysis of DE lncRNAs and DE mRNAs

The GO enrichment analysis displayed that the DE lncRNAs were mainly enriched in biological processes including positive regulation of NF‐kappa B transcription, inflammatory response, and signal transduction; enriched in molecular functions consisting of identical protein binding, signal transducer activity, and protease binding; in addition, the DE lncRNAs were mostly located at cell surface, plasma membrane, and extracellular exosomes. (Figure 4A). The KEGG enrichment analysis illuminated that the DE lncRNAs were predominantly enriched in inflammatory and bone formation‐related pathways such as the NF‐kappa B signaling pathway, Wnt signaling pathway, and osteoclast differentiation. (Figure 4B). In terms of DE mRNAs, the GO enrichment elucidated that they were notably enriched in the biological processes of positive regulation of NF‐kappaB transcription, inflammatory response, and ossification; enriched in the molecular functions of identical protein binding, calcium ion binding, and protease binding; and they were mainly located in cell components of plasma membrane, cell surface, and extracellular exosome. (Figure 5A). Besides, the KEGG enrichment analysis showed that the DE mRNAs were markedly enriched in NF‐kappaB signaling pathway, pathways in cancer and Wnt signaling pathway, etc (Figure 5B). Furthermore, the top 20 DE lncRNAs including 10 upregulated lncRNAs and 10 downregulated lncRNAs were selected according to the rank of absolute value of Log2FC (Table 2), and the regulatory network between the top 20 DE lncRNAs and DE mRNAs was presented in Figure 6.

FIGURE 4

Enrichment analysis of DE lncRNAs. GO enrichment (A) and KEGG enrichment (B) of DE lncRNAs. BP, biological processes; CC, cell component; DE, differentially expressed; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long non‐coding RNA; MF, molecular function

FIGURE 5

Enrichment analysis of DE mRNAs. GO enrichment (A) and KEGG enrichment (B) of DE mRNAs. BP, biological processes; CC, cell component; DE, differentially expressed; GO, Gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function

TABLE 2

Top 20 DE lncRNAs (10 upregulated and 10 downregulated) by RNA‐Seq

Gene‐ID	Symbol	Log2FC	P value	Adjusted P value	Trend
ENSG00000228769	LncRNA ABCF2P2	5.799762	.001545	.021914	Up
ENSG00000234296	LncRNA RP13‐16H11.7	5.731294	.001379	.020256	Up
ENSG00000241956	LncRNA CTC‐340A15.2	5.594905	.001511	.021592	Up
ENSG00000228420	LncRNA RP4‐735C1.6	5.507139	.001803	.02441	Up
ENSG00000233651	LncRNA RP11‐111F5.8	5.385424	8.88E‐09	2.47E‐06	Up
ENSG00000229950	LncRNA TFAP2A‐AS1	5.276863	8.73E‐05	.002719	Up
ENSG00000228613	LncRNA AC144450.1	5.055785	7.38E‐06	.000418	Up
ENSG00000260725	LncRNA AC005307.1	4.956591	.000805	.013847	Up
ENSG00000240411	LncRNA RPL5P16	4.89303	4.86E‐06	.00031	Up
ENSG00000231291	LncRNA RP11‐445O3.1	4.830801	3.93E‐06	.000264	Up
ENSG00000236993	LncRNA GAPDHP21	−5.77175	2.62E‐06	.000189	Down
ENSG00000270710	LncRNA RP11‐293G6‐B.8	−5.58435	9.51E‐08	1.48E‐05	Down
ENSG00000273230	LncRNA RP11‐1246C19.1	−5.58141	6.22E‐06	.000369	Down
ENSG00000261838	LncRNA RP11‐303E16.6	−5.21083	5.36E‐07	5.57E‐05	Down
ENSG00000260793	LncRNA RP5‐882C2.2	−5.14379	2.36E‐07	3.03E‐05	Down
ENSG00000224116	LncRNA INHBA‐AS1	−5.14119	.00054	.010446	Down
ENSG00000232407	LncRNA RP11‐641C17.2	−4.94134	1.20E‐06	.000105	Down
ENSG00000272367	LncRNA CTC‐428H11.2	−4.93645	1.67E‐05	.000775	Down
ENSG00000224934	LncRNA RP11‐441O15.3	−4.83866	.00084	.014233	Down
ENSG00000237803	LncRNA LINC00211	−4.79562	.000111	.003276	Down

Top 10 upregulated and 10 downregulated lncRNAs were selected by the rank of absolute value of Log2FC.

Abbreviations: DE lncRNAs, differentially expressed long non‐coding RNAs; FC, fold change; RNA‐Seq, RNA sequencing.

FIGURE 6

Regulatory network showing the interactions of the top 20 DE lncRNAs with DE mRNAs. DE, differentially expressed; lncRNA, long non‐coding RNA

Demographic characteristics of KOA patients and Ctrls in RT‐qPCR

In the total 100 KOA patients and 100 Ctrls, the mean age was 52.9 ± 6.2 and 51.6 ± 8.4 years, respectively (P = .228) (Table 1). There were 73 (73.0%) females and 27 (27.0%) males in KOA patients, 67 (67.0%) females and 33 (33.0%) males in Ctrls (P = .355). In addition, a mean BMI of 23.2 ± 2.3 kg/m2 was found in KOA patients, and the Ctrls had a mean BMI of 22.9 ± 2.5 kg/m2 (P = .291).

Comparison of candidate lncRNAs between KOA and Ctrls

RT‐qPCR validation of 10 candidate lncRNAs, which included top five upregulated lncRNAs and top five downregulated lncRNAs, was conducted in a large sample population containing 100 KOA patients and 100 Ctrls. And the results of RT‐qPCR validation showed that lncRNA ABCF2P2 (P < .001) (Figure 7A), lncRNA RP13‐16H11.7 (P < .001) (Figure 7B), lncRNA CTC‐340A15.2 (P = .001) (Figure 7C) and lncRNA RP4‐735C1.6 (P < .001) (Figure 7D) were upregulated, while lncRNA RP11‐293G6‐B.8 (P < .001) (Figure 7G), lncRNA RP11‐1246C19.1 (P < .001) (Figure 7H), lncRNA RP11‐303E16.6 (P < .001) (Figure 7I) and lncRNA RP5‐882C2.2 (P = .026) (Figure 7J) were downregulated in KOA patients compared to Ctrls. However, lncRNA RP11‐111F5.8 (P = .882) (Figure 7E) and lncRNA GAPDHP21 (P = .157) (Figure 7F) were of no difference between KOA patients and Ctrls.

FIGURE 7

Ten candidate lncRNAs expressions in plasma by RT‐qPCR conducted in large sample size populations. The expressions of lncRNA ABCF2P2 (A), lncRNA RP13‐16H11.7 (B), lncRNA CTC‐340A15.2 (C), lncRNA RP4‐735C1.6 (D), lncRNA RP11‐111F5.8 (E), lncRNA GAPDHP21 (F), lncRNA RP11‐293G6‐B.8 (G), lncRNA RP11‐1246C19.1 (H), lncRNA RP11‐303E16.6 (I) and lncRNA RP5‐882C2.2 (J) in KOA patients and Ctrls. Ctrl, control; KOA, knee osteoarthritis; lncRNAs, long non‐coding RNAs; RT‐qPCR, reverse transcription–quantitative polymerase chain reaction

ROC curve analyses of candidate lncRNAs for KOA risk

According to ROC curve analyses, lncRNA ABCF2P2, lncRNA RP13‐16H11.7, lncRNA RP4‐735C1.6 and lnc RNA CTC‐34‐A15.2 could predict increased KOA risk, and their AUC values were 0.847 (95% CI: 0.794‐0.900), 0.766 (95% CI: 0.699‐0.834), 0.819 (95% CI: 0.762‐0.876) as well as 0.628 (95% CI: 0.551‐0.705), respectively (Figure 8A). Whereas lncRNA RP11‐293G6‐B.8, lncRNA RP11‐1246C19.1, lncRNA RP11‐303E16.6 and lncRNA RP5‐882C2.2 exhibited values for predicting decreased KOA risk, the AUCs were 0.825 (95% CI: 0.768‐0.881), 0.801 (95% CI: 0.741‐0.862), 0.754 (95% CI: 0.687‐0.820) and 0.600 (95% CI: 0.521‐0.678), respectively (Figure 8B). However, lncRNA RP11‐111F5.8 (Figure 8A) and lncRNA GAPDHP21 (Figure 8B) exhibited no value for predicting KOA risk, and their AUC values were 0.533; 95% CI: 0.452‐0.613 as well as 0.548; 95% CI: 0.468‐0.628, respectively.

FIGURE 8

ROC curve analysis of 10 candidate lncRNAs. The ROC curve analysis of 5 upregulated candidate lncRNAs (A) and 5 downregulated candidate lncRNAs (B). lncRNA, long non‐coding RNA; ROC, receiver operating characteristic

Values of candidate lncRNAs for predicting KOA risk

The univariate logistic regression model illuminated that lncRNA ABCF2P2 (P < .001, OR = 4.270), lncRNA RP13‐16H11.7 (P < .001, OR = 3.253), lncRNA CTC‐340A15.2 (P < .001, OR = 1.995) and lncRNA RP4‐735C1.6 (P < .001, OR = 4.100) correlated with higher risk of KOA, while lncRNA RP11‐293G6‐B.8 (P < .001, OR = 0.050), lncRNA RP11‐1246C19.1 (P < .001, OR = 0.081), lncRNA RP11‐303E16.6 (P < .001, OR = 0.257) and LncRNA RP5‐882C2.2 (P = .045, OR = 0.652) associated with lower risk of KOA (Table 3). However, lncRNA RP11‐111F5.8 (P = .072, OR = 1.368), lncRNA GAPDHP21 (P = .311, OR = 0.795), age (P = .227, OR = 1.024), male gender (P = .355, OR = 0.751) and BMI (P = .290, OR = 1.065) were not correlated with KOA risk. Subsequently, multivariate logistic regression using backward stepwise method demonstrated that lncRNA ABCF2P2 (P < .001, OR = 39.098), lncRNA RP13‐16H11.7 (P < .001, OR = 19.095) and lncRNA CTC‐340A15.2 (P = .019, OR = 17.057) were independently correlated with increased risk of KOA, nonetheless, lncRNA RP11‐293G6‐B.8 (P < .001, OR = 0.001), lncRNA RP11‐1246C19.1 (P = .003, OR = 0.004) and lncRNA RP11‐303E16.6 (P = .001, OR = 0.003) were independently associated with decreased risk of KOA. Furthermore, ROC curve analysis revealed that the combination of these 6 independent lncRNAs presented with great value in predicting KOA risk with AUC 0.982, 95% CI: 0.964‐1.000 (Figure 9).

TABLE 3

Predicting values of 10 candidate lncRNAs for KOA risk

Items	Logistic regression model
	P value	OR	95% CI
			Lower	Higher
Univariate logistic regression
LncRNA ABCF2P2	<.001	4.270	2.696	6.764
LncRNA RP13‐16H11.7	<.001	3.253	2.214	4.780
LncRNA CTC‐340A15.2	<.001	1.995	1.362	2.923
LncRNA RP4‐735C1.6	<.001	4.100	2.593	6.485
LncRNA RP11‐111F5.8	.072	1.368	0.973	1.925
LncRNA GAPDHP21	.311	0.795	0.511	1.239
LncRNA RP11‐293G6‐B.8	<.001	0.050	0.021	0.122
LncRNA RP11‐1246C19.1	<.001	0.081	0.037	0.178
LncRNA RP11‐303E16.6	<.001	0.257	0.150	0.441
LncRNA RP5‐882C2.2	.045	0.652	0.430	0.990
Age	.227	1.024	0.985	1.063
Male	.355	0.751	0.409	1.378
BMI	.290	1.065	0.948	1.198
Backward stepwise multivariate logistic regression
LncRNA ABCF2P2	<.001	39.098	5.171	295.637
LncRNA RP13‐16H11.7	<.001	19.095	4.192	86.973
LncRNA CTC‐340A15.2	.019	17.057	1.582	183.898
LncRNA RP11‐293G6‐B.8	<.001	0.001	0.000	0.046
LncRNA RP11‐1246C19.1	.003	0.004	0.000	0.150
LncRNA RP11‐303E16.6	.001	0.003	0.000	0.082

Predicting values of 10 candidate lncRNAs for KOA risk were analyzed by univariate and backward stepwise multivariate logistic regression models.

Abbreviations: BMI, body mass index; CI, confidence interval; KOA, knee osteoarthritis; lncRNA, long non‐coding RNA; OR, odds ratio; RT‐qPCR, reverse transcription–quantitative polymerase chain reaction.

FIGURE 9

ROC curve analysis of predictive model based on independent predictive lncRNAs on KOA risk. KOA, knee osteoarthritis; lncRNA, long non‐coding RNA; ROC, receiver operating characteristic

DISCUSSION

In this study, we discovered that (a) the lncRNA expression profile and mRNA expression profile both distinguished KOA patients from Ctrls; (b) the DE lncRNAs were mainly located in cell surface, membrane, and exosomes, and were predominantly enriched in biological processes related to inflammation, molecular functions of protein binding and signaling pathways of inflammation and bone formation; (c) RT‐qPCR validation in large sample size revealed that among 10 candidate lncRNAs, 8 candidate lncRNAs were markedly dysregulated in KOA patients, and they also predicted KOA risk; in addition, six of them were independent predictive factors for KOA risk.

Some recent reports have presented findings of lncRNA expression profile in OA patients. However, little effort is done aiming to investigate lncRNA expression profile in KOA. With regard to OA, a previous study elucidates that there are 7082 differentially expressed RNAs (including lncRNAs) in mice model of OA, then the OA‐related core RNAs (including lncRNAs) are identified, which are mainly enriched in signaling pathways related to inflammation and bone formation. Another study reveals 2042 differentially expressed lncRNAs in OA cartilage (collected from OA patients) compared with non‐OA cartilage (collected from patients who undergo invasive procedures), including 1137 upregulated lncRNAs and 905 downregulated lncRNAs. As for KOA, a prior study reveals that in cartilage, 3007 upregulated lncRNAs and 1707 downregulated lncRNAs are in KOA patients compared with Ctrls. In addition, a recent study reports that there are 1870 abnormally expressed lncRNA‐mRNA interactions, which consist of 476 gain interactions and 1394 loss interactions involving 131 lncRNAs and 1251 mRNAs in KOA patients with mild and severe pain. Besides, another study discovers two differentially expressed lncRNAs in synovial tissue of KOA patients with high pain (visual analogue scale score 7‐10) compared to KOA patients with low pain (visual analogue scale score 0‐3) using the hypothesis‐free next‐generation RNA‐sequation These findings indicate the promising role of lncRNA expression profile in KOA. In our study, we found that the lncRNA expression profile as well as mRNA expression profile in plasma could differentiate KOA patients from Ctrls. Furthermore, we also identified 418 upregulated lncRNAs and 347 downregulated lncRNAs in KOA patients, which were mostly located in surface, plasma membrane and extracellular exosome of cells, and were mainly enriched in biological processes of inflammation, molecular functions of protein binding as well as signaling pathways of inflammation and bone formation. Compared with the previous studies, the number of DE lncRNAs found in our study was quite different, either larger or smaller, which might be caused by the discrepancies in sample size, type of samples, or RNA‐seq techniques. In addition, our study was conducted in circulating samples (plasma) but not tissues (cartilage tissue or synovial tissue). This might provide more valuable information about the utility of lncRNA expression profile in KOA management compared with the previous similar studies, since that circulating sample was easier to be obtained from patients in the clinical practice.

The value of lncRNAs for KOA management is a promising area which needs further exploration. There are only several studies uncovering the diagnostic values of lncRNAs for KOA. For instance, a study reveals that the putative functional variants of lncRNA are correlated with risk of KOA, including H19 rs2067051 T allele and MEG3 rs4378559 T allele (the former associates with lower KOA risk while the latter correlates with higher KOA risk). Another study reports that lnc‐MFI2‐AS1 level is elevated in cartilage tissue of KOA patients compared to normal cartilage tissue (obtained from undamaged areas of KOA patients). In our study, the RT‐qPCR validation in large sample size revealed that among 10 candidate lncRNAs, there were four upregulated and four downregulated candidate lncRNAs in KOA patients compared to Ctrls, and ROC curve disclosed that these eight candidate lncRNAs could predict KOA risk. Further multivariate logistic regression elucidated that six of these eight candidate lncRNAs independently predicted the KOA risk. These results indicated that the six lncRNAs might serve as biomarkers assisting the KOA diagnosis in clinical practice, which, however, needed to be verified in a larger cohort of KOA patients. As for the probable explanation to these results, it might be derived from the following fact: the target mRNAs of these six lncRNAs were involved in the regulation of inflammation and bone formation as reported previously. For instance, the target mRNA of lncRNA ABCF2P2, angiopoietin‐like 2 (ANGPTL2), the reduction of which could repress osteoclast generation via regulating nuclear factor‐kappaB (NF‐κB)/mitogen‐activated protein kinases (MAPKs)/cyclin signaling pathway, and the target mRNA of LncRNA RP11‐303E16.6, receptor‐interacting kinase 1 (RIPK1), is a crucial modulator of inflammation. , , , , , , , , , , , , Therefore, these six lncRNAs might regulate inflammation and bone formation in KOA through interacting with their target mRNAs, which contributed to their predictive values for KOA risk. Besides, since age, gender and BMI were already matched between KOA patients and Ctrls, it might not be necessary to include them in the multivariate analysis. While as the recommendation by the statistical expert of our hospital, the addition of them in the multivariate would further diminish the influence of age, gender, and BMI on the KOA risk. In addition, since the original aim of our study was to investigate the dysregulated lncRNA expression profile in KOA patients vs Ctrls, we did not document the disease assessment score or disease course when these KOA patients were enrolled, thus the correlation of candidate lncRNAs with disease severity and course of KOA was not analyzed in this study.

Collectively, this study revealed an aberrant plasma lncRNA expression profile in KOA patients and identified several specific lncRNAs which could assist in KOA management in clinical practice.

Table S1

Click here for additional data file.