Comprehensive analyses of long non-coding RNA expression profiles by RNA sequencing and exploration of their potency as biomarkers in psoriatic arthritis patients

Background The aim of the current study was to investigate the long non-coding RNA (lncRNA) expression profiles in psoriatic arthritis (PSA) patients by RNA sequencing, and to further explore potential biomarkers that were able to predict PSA risk and activity. Methods LncRNA and mRNA expression profiles in peripheral blood mononuclear cells (PBMC) of 4 PSA patients and 4 normal controls (NCs) were detected by RNA sequencing, followed by comprehensive bioinformatic analyses. Subsequently, 3 top upregulated and 2 top downregulated lncRNAs were chosen for further validation in 93 PSA patients and 93 NCs by quantitative polymerase chain reaction (qPCR) assay. Results Totally 76 upregulated and 54 downregulated lncRNAs, as well as 231 upregulated and 102 downregulated mRNAs were discovered in PSA patients compared with NCs. Enrichment analyses revealed that they were mostly associated with nucleosome, extracellular exosome and extracellular matrix, and the top enriched pathways were systemic lupus erythematosus (SLE), alcoholism and viral carcinogenesis. qPCR assay showed that lnc-RP11-701H24.7 and lnc-RNU12 were upregulated in PSA patients compared with NCs, and they could predict PSA risk with high area under curves. Besides, lnc-RP11-701H24.7 was positively associated with ESR, SJC, TJC and pain VAS score while lnc-RNU12 was positively correlated with PASI score, CRP and PGA score, implying that both of them were positively correlated with disease activity. Conclusion Our study facilitates comprehensive understanding of lncRNA expression profiles in PSA pathogenesis, and discovers that lnc-RP11-701H24.7 and lnc-RNU12 might be served as novel biomarkers for PSA risk and activity. Electronic supplementary material The online version of this article (10.1186/s12865-019-0297-9) contains supplementary material, which is available to authorized users.


Background
Psoriatic arthritis (PSA), a complicated and heterogeneous autoimmune disease that causes severe joint damage and reduces quality of life, is characterized by different symptoms such as enthesitis, dactylitis, nail dystrophy, uveitis, osteitis and several comorbidities including cardiovascular disease (CAD), obesity and metabolic syndrome [1,2]. It's reported that approximately 30% of psoriatic patients would occur PSA in their lifetime and there is no "golden standard" diagnostic test for PSA in clinical practices so far, thus the diagnosis that only relies on various pieces of evidence remains a challenge [2][3][4][5][6]. More importantly, the detailed mechanisms underlying pathogenesis of PSA are still poorly illuminated, which greatly obstruct the improvement of long-term outcome of PSA patients [2].
Long non-coding RNA (lncRNA), a class of RNAs that have more than 200 nucleotides and lack open reading frames, is implicated different cancers, CAD and autoimmune diseases [7,8]. In the past few decades, lncRNA expression profiles of autoimmune diseases such as systemic lupus erythematosus (SLE), rheumatic arthritis (RA) and inflammatory bowel disease (IBD) have been investigated via microarray or RNA sequencing, through which the roles of lncRNA expression profiles in etiology of several autoimmune diseases have been disclosed and a few lncRNAs that present with predicting values for disease risk, inflammation level and activity have been discovered [9][10][11]. For instance, an interesting study investigates the lncRNA expression profiles in T cells of SLE patients using microarray, which finds 1935 differentially expressed lncRNAs (DELs) in SLE patients compared with healthy controls (HCs), further analyses reveal that both uc001ykl.1 and ENST00000448942 expressions could differentiate SLE patients from controls and positively correlate with disease activity as well as inflammation level [9]. In another study, lncRNA expression profiles in peripheral blood mononuclear cells (PBMC) of RA patients are investigated, then 5045 DELs are identified in RA patients compared with HCs, enrichment analyses indicate that these DELs are mostly implicated bile secretion, T cell receptor signaling pathways and SLE [10].
Considering that PSA is also a polygenic autoimmune disease, exploring lncRNA expression profiles of PSA might also contribute to unveiling development and progression of PSA and discovering biomarkers for predicting PSA risk and activity. However, just one study explores the lncRNA expression profiles in PBMC of PSA patients by microarray, and the DELs in this study are only validated in 20 PSA patients and 20 controls by quantitative polymerase chain reaction (qPCR) assay [12]. Thus, the in-depth knowledge of lncRNA expression profile in PSA pathogenesis remains to be further investigated, and lncRNAs that are able to predict PSA risk and disease activity remain to be discovered. Therefore, the objective of the current study was to explore lncRNA expression profiles in PBMC of PSA patients by RNA sequencing, and to further explore biomarkers (from DELs of PSA) that were able to predict PSA risk and activity.

Materials and methods
Patients and samples

Study design and ethics approval
The detailed study design was exhibited in Fig. 1

RNA sequencing
Automate quality control and adapter trimming were conducted using Trim Galore, Cutadapt and FastQC. The trimmed reads were mapped to the human genome Hg38 by HISAT2 with default parameters according to the methods in a previous report [13], and mapping quality control was performed using RSeQC referring to the methods in another previous paper [14]. The read counts of lncRNA and mRNA were then calculated using feature-Counts based on the annotation file (Homo_ sapiens.GRCh38.83.gtf) in Ensembl database according to the methods described in a previous report [15]. In this study, lncRNAs and mRNAs discovered in less than 50% samples were wiped off for analysis, only the remaining lncRNAs and mRNAs were used for follow-up analysis, and their raw reads counts were normalized, logarithmic transformation was applied, and DELs and differentially expressed mRNAs (DEMs) were detected using DeSeq2 as method that was previous described [16], and the statistical Fig. 3 Volcano plot and heatmap analysis. Volcano plot showed that 76 lncRNAs were upregulated and 54 lncRNAs were downregulated in PSA patients compared with NCs (a). Heatmap analysis revealed that these DELs were able to separate PSA patients from NCs (b). 231 mRNAs were upregulated while 102 mRNAs were downregulated in PSA patients compared with NCs (c), and these DEMs could also separate PSA patients from NCs (d). lncRNA, long non-coding RNA; PSA, psoriatic arthritis; NCs, normal controls; DELs, differentially expressed lncRNAs; DEMs, differentially expressed mRNAs significance was defined as adjusted P value < 0.05 (Benjamini & Hochberg (BH) adjusted method was applied to reduce the false positive result), and the biological significance was defined as a difference of at least 2.0 folds, which meant abs (log 2 (fold change)) > =1.0.

Bioinformatic analyses
Principal component analysis (PCA) of lncRNA and mRNA expression patterns was performed using Stats package. Heatmap plot analysis of total expression pattern, DELs and DEMs were performed using Pheatmap package. Gene Ontology (GO) and Kyoko Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of DELs were performed using DAVID web servers referring to the method in a previous report [17]. Regulatory networks of DELs by cis targets and trans targets were drawn using igraph package. Cis targets of lncRNAs were determined by the location and Trans targets of lncRNAs were determined by pearson correlation coefficient of expressing pattern. All the bioinformatics analysis was performed using R software (Version 3.3.3).

Data collection
Comprehensive data of PSA patients were collected including age, gender, family history of PSA, disease duration of psoriasis, disease duration of arthritis, psoriasis subtype, arthritis subtype, affected body surface area (BSA), psoriasis area and severity index (PASI) score, Creactive protein (CRP), erythrocyte sedimentation rate (ESR), swollen joint count (SJC), tender joint count (TJC), Pain visual analogue scale (VAS) score and patients' global assessment (PGA) score.

qPCR assay
Total RNAs were extracted from PBMC samples of 93 PSA patients and 93 NCs using Trizol reagent (Invitrogen, USA) according to manufacturer's instructions, Qubit® 2.0 Flurometer (Life Technologies, USA) was applied for examination of RNA concentration and purity. Subsequently, cDNA was synthesized by QuantiNova Reverse Transcription Kit (Qiagen, German) then it was subjected to qPCR with SYBR Green kit (TaKaRa, Japan). The protocol of qPCR included a single cycle at 95°C for 5 mins, followed by 40 cycles at 95°C for 10 s, and annealing at 60°C for 30 s. After that, the expressions of 5 selected lncRNAs were calculated using the 2 -ΔΔCt methods with phosphoglyceraldehyde dehydrogenase (GAPDH) as internal reference. Primers used in qPCR validation were listed in Additional file 1: Table S1.

Statistics
Statistical analysis was conducted using SPSS 22.0 software (IBM Corp, USA) and GraphPad Prism 5.01

PCA plot and heatmap analysis of lncRNA and mRNA profiles
To determine whether these 8 samples of PSA patients and NCs could be grouped, PCA was performed, which showed a clear segregation between 4 PSA patients and 4 NCs by lncRNA profiles (Fig. 2a). The heatmap analysis of lncRNA profiles showed that lncRNA profiles were able to differentiate PSA patients from NCs (Fig. 2b). PCA plot (Fig. 2c) and heatmap analysis (Fig. 2d) of mRNAs profiles disclosed that mRNA profiles also clearly segregated PSA patients from NCs.

Volcano plot and heatmap analysis for DELs and DEMs
A total of 10,079 lncRNAs which were detected in no less than 50% of samples were identified and included into analysis, then 76 remarkably upregulated and 54 remarkably downregulated lncRNAs were found in PSA patients compared with NCs using BH adjusted method (p < 0.05 and fold change≥2) (Fig. 3a). And the top 10 upregulated and top 10 downregulated lncRNAs were  Table 1. Heatmap analysis revealed that these DELs could well separate PSA patients from NCs (Fig. 3b). As to DEMs, totally 13,737 mRNAs were identified, among which 231 mRNAs were obviously upregulated while 102 mRNAs were obviously downregulated in PSA patients compared with NCs (Fig. 3c), and these DEMs were also able to separate PSA patients from NCs (Fig. 3d).

GO and KEGG enrichment analyses for DELs and DEMs
The enrichment analyses are commonly used to illuminate the molecular functions, biological processes, cellular components and signaling pathways that DELs and DEMs might be implicated. In the current study, GO enrichment analysis for DELs showed that DELs were mostly associated with nucleosome, extracellular exosome and extracellular matrix (Fig. 4a), KEGG enrichment analyses for DELs revealed that DELs were mainly enriched in SLE, alcoholism and viral carcinogenesis (Fig. 4b).
As for DEMs, their enrichment analysis results were consistent with that of DELs (Figs. 3c, 4d). These results illustrated the potential pathways and functions that the DELs and DEMs in PSA were involved.

Regulatory networks of DELs
In order to better understand the functions of DELs in pathogenesis of PSA, the regulatory networks of DELs with their target mRNAs were drawn, including cis targets (Fig. 5a) and trans targets (Fig. 5b). Some DELs have no target mRNA (not displayed in the networks), while other DELs directly regulated one or more target mRNAs. The dots represented the DELs, the squares represented mRNAs. Upregulated genes were colored red, downregulated genes were colored blue, and the regulation insignificant genes were colored green.

Characteristics of 93 PSA patients in validation stage
To explore potential values of 5 candidate lncRNAs in predicting disease risk and activity of PSA, 93 PSA patients were recruited for further validation, and the characteristics of these PSA patients were displayed in  Fig. 6a) and lnc-RNU12 (P < 0.001, Fig. 6b) levels were increased in PSA patients compared with NCs, while lnc-SNORD3A (P = 0.215, Fig. 5 Regulatory networks. Regulatory networks of DELs with their target mRNAs were drawn, including cis targets (a) and trans targets (b). The dots stood for the DELs, the squares stood for mRNAs; upregulated genes were colored red, downregulated genes were colored blue, and the regulation insignificant genes were colored green. DELs Differentially expressed lncRNAs Fig. 6c), lnc-TRAV1-2 (P = 0.096, Fig. 6d) and lnc-TRAV1-1 (P = 0.320, Fig. 6e) expressions between two groups were similar.

Discussion
LncRNA, a group of RNAs that is described in 2002 for the first time, presents with several characteristics, including poor conservation, tissue specificity and so on [7,18,19]. LncRNAs exist in almost every branch of life and are involved various important biological activities such as genomic imprinting, cell differentiation and organogenesis [20,21]. LncRNAs regulating gene expressions are mainly through three different ways: transcription regulation, posttranscription regulation or epigenetic regulation [22]. For transcription regulation, lncRNAs are able to regulate different transcription processes such as modulate transcription factor activity, mediate RNA polymerase (RNAP) II activity and regulate the functions of RNA binding protein.
For post-transcription regulation, lncRNAs also modulate various post-transcriptional processes, such as splicing, editing, transport, translation and degradation. Besides, lncRNAs are also involved epigenetic regulations such as imprinting and X-chromosome inactivation by recruiting chromatin remodeling complexes to specific genomic loci [23,24]. In recent years, the role of lncRNAs in a number of diseases have been intensively investigated, especially in complicated diseases such as cancers, CAD and autoimmune diseases [18]. Taken psoriasis as an example, a study explores the lncRNA expression profiles in skin sample of psoriasis patients, which identifies 1080 novel differentially expressed lncRNAs in psoriasis patients compared with controls, further analyses revealed that these dysregulated lncRNAs are mostly enriched in biosynthesis of unsaturated fatty acid and cytokine activity [25]. However, the implications of lncRNA expression profiles in PSA are still obscure. PSA is a systemic inflammatory disease that both genetic and environmental factors play crucial roles in the etiology of this disease [5]. Since lncRNA expression profiles are reported to be paramount in etiopathologies of autoimmune diseases such as SLE, IBD and RA, exploring lncRNA expression profiles of PSA might also be able to provide new insights for PSA pathogenesis and find novel biomarkers for predicting PSA risk and activity [9][10][11]. Whereas relevant study is rarely found until now, with only one study reported [12]. In this previous study, lncRNA expression profiles in PBMC of 10 PSA patients are analyzed by microarray, which discovers 259 DELs in PSA patients compared with HCs, further analyses reveal that these DELs are implicated PSArelated signaling pathways such as immune response, glycolipid metabolism, bone remodeling and type 1 interferon. Subsequently, 7 lncRNAs (TRIM55-1, EPB41L4A-AS1, LINC00657, LINC00909, RP11-539 L10.3, LA16c-360H6.3 and LUCAT1) are selected for validation in 20 PSA patients and 20 HCs, which discloses that LA16c-360H6.3 is upregulated while the rest of 6 lncRNAs are downregulated in PSA patients compared with HCs [12]. However, this study uses microarray to detect DELs rather than RNA sequencing. What's more, the sample size for validation is small, with only 20 samples in each group, which also decreases statistical power. Considering that no more studies are reported to explore the role of lncRNA expression profiles in PSA, and the potential pathways or mechanisms that lncRNA expression profiles might be involved are largely unelucidated, we conducted the current study by RNA sequencing, which discovered 130 DELs in PSA patients compared to NCs, among which 76 lncRNAs were upregulated while 54 lncRNAs were downregulated. The number of DELs in our study was smaller than the previous study, which might be due to the following reasons: (1) our study discarded those lncRNAs that were discovered in less than 50% samples, whereas the previous study does not. (2) Our study utilized RNA sequencing to detect DELs, while the reported study uses microarray to detect DELs, which might cause differences.
(3) DELs in our study must meet an adjusted P value < 0.05 (BH adjusted method) and a fold change≥2.0, while BH adjusted method is not used in their study. Taken together, the DELs in the current study were relatively fewer than that of reported study. In addition, we compared our data-set with that previous study, discovering that 1 DEL and 41 DEMs in our data-set were overlapped with that previous study, indicating that the remaining 129 DELs and 292 DEMS in our study were first identified. The enrichment analyses are prevalent approaches to illuminate the molecular functions, biological processes, cellular components and signaling pathways that differentially expressed genes might be implicated [26,27]. In the current study, enrichment analyses disclosed that DELs and DEMs were mostly associated with nucleosome, extracellular exosome and extracellular matrix, and the top enriched pathways were SLE, alcoholism and viral carcinogenesis. The possible explanations for the results might be that: for nucleosome, it is reported that its loss contributes to macrophage activities, indicating that DELs involved PSA might be through regulating DEMs then affecting the numbers of nucleosome thereby activating macrophages [28]. For extracellular exosome, it is discovered to induce inflammatory immune response in a couple of autoimmune diseases such as SLE and RA, and extracellular exosome in PSA patients is able to stimulate osteoclast differentiation, suggesting that DELs implicated PSA might also be through regulating DEMs then modulating exosome thereby mediated osteoclast differentiation [29,30]. For extracellular matrix, it serves important functions to cell adhesion, cell-to-cell communication and cell differentiation, and the dysregulated degradation of extracellular matrix is associated with a number of diseases such as osteoarthritis and PSA, therefore DELs involved PSA might also be through regulating DEMs then mediating degradation of extracellular matrix [31,32]. For SLE, it closely correlates with autoimmunity and inflammation, which are similar to PSA. As for alcoholism and viral carcinogenesis, their associations with PSA need further investigation.
The disease manifestations of PSA are sometimes varied from patients to patients, which make the diagnosis difficult. To explore if there are DELs that could be served as biomarkers for disease risk, inflammation or activity of PSA, 3 top upregulated and 2 downregulated lncRNAs were selected and their expressions were validated in 93 PSA patients and 93 NCs by qPCR assay, which revealed that both lnc-RP11-701H24.7 and lnc-RNU12 expressions were increased in PSA patients compared with NCs, and they disclosed good predictive values for PSA risk with high AUCs. Besides, lnc-RP11-701H24.7 and lnc-RNU12 were also positively associated with inflammation level and disease activity of PSA. The possible reasons might be that: both lnc-RP11-701H24.7 and lnc-RNU12 might be implicated PSA development and progress through (i) epigenetic regulations such as histone and DNA methylation.
(ii) Mediating gene transcription; (iii) regulating its target genes such as COL1A2, RNU12 and MAP1A. However, the in-depth mechanisms of lnc-RP11-701H24.7 and lnc-RNU12 being implicated PSA pathogenesis need further investigation [33]. In addition, we also observed that lnc-TRAV1-2 was negatively associated with PASI score and PGA score, indicating that lnc-TRAV1-2 might be also correlated with decreased disease activity of PSA. In short, it was the first study to discover specific lncRNAs that might be regarded as biomarkers for PSA risk and activity, which might remarkably promote the early diagnosis and contribute to discovering novel approaches for improving long-term outcomes of PSA patients. There were some limitations in the current study. First, the study did not evaluate the treatment response of PSA patients, so the role of lncRNA expression profiles in treatment response and long-term outcome of PSA patients were not investigated. Second, most patients in this study were from East China, which might bring in selection bias. Third, although we discovered the positive correlations of lnc-RP11-701H24.7 and lnc-RNU12 with disease risk and activity, the detailed mechanisms of these two lncRNAs in development and progression of PSA remained unclear, which needed to be further explored in our future work.

Conclusion
In summary, our study facilitates comprehensive understanding of lncRNA expression profiles in PSA development and progression, and discovers that lnc-RP11-701H24.7 and lnc-RNU12 might be served as novel biomarkers for PSA risk and activity.

Additional file
Additional file 1: Table S1.