Long Non-Coding RNAs Play a Role in the Pathogenesis of Psoriatic Arthritis by Regulating MicroRNAs and Genes Involved in Inflammation and Metabolic Syndrome.

Psoriatic arthritis (PsA) is an inflammatory arthritis, characterized by inflammation of entheses and synovium, leading to joint erosions and new bone formation. It affects 10-30% of patients with psoriasis, and has an estimated prevalence of approximately 1%. PsA is considered to be primarily an autoimmune disease, driven by autoreactive T cells directed against autoantigens present in the skin and in the joints. However, an autoinflammatory origin has recently been proposed. Long noncoding RNAs (lncRNAs) are RNAs more than 200 nucleotides in length that do not encode proteins. LncRNAs play important roles in several biological processes, including chromatin remodeling, transcription control, and post-transcriptional processing. Several studies have shown that lncRNAs are expressed in a stage-specific or lineage-specific manner in immune cells that have a role in the development, activation, and effector functions of immune cells. LncRNAs are thought to play a role in several diseases, including autoimmune disorders. Indeed, a few lncRNAs have been identified in systemic lupus erythematosus, rheumatoid arthritis, and psoriasis. Although several high-throughput studies have been performed to identify lncRNAs, their biological and pathological relevance are still unknown, and most transcriptome studies in autoimmune diseases have only assessed protein-coding transcripts. No data are currently available on lncRNAs in PsA. Therefore, by microarray analysis, we have investigated the expression profiles of more than 50,000 human lncRNAs in blood samples from PsA patients and healthy controls using Human Clariom D Affymetrix chips, suitable to detect rare and low-expressing transcripts otherwise unnoticed by common sequencing methodologies. Network analysis identified lncRNAs targeting highly connected genes in the PsA transcriptome. Such genes are involved in molecular pathways crucial for PsA pathogenesis, including immune response, glycolipid metabolism, bone remodeling, type 1 interferon, wingless related integration site, and tumor necrosis factor signaling. Selected lncRNAs were validated by RT-PCR in an expanded cohort of patients. Moreover, modulated genes belonging to meaningful pathways were validated by RT-PCR in PsA PBMCs and/or by ELISA in PsA sera. The findings indicate that lncRNAs are involved in PsA pathogenesis by regulating both microRNAs and genes and open new avenues for the identification of new biomarkers and therapeutical targets.

Introduction

Psoriatic arthritis (PsA) is a chronic, immune-mediated, asymmetric inflammatory arthritis characterized by inflammation at tendon or ligament insertion sites into bone (enthesitis) and by synovitis, eventually leading to joint erosions and new bone formation (1).

Up to 30% of patients with skin psoriasis may develop PsA and its prevalence is estimated in 1% in the general population. PsA shares genetic and clinical features with other forms of seronegative spondyloarthritis (2, 3). Diagnostic criteria for PsA have not been validated, but the Classification Criteria for PsA (CASPAR criteria), published in 2006, define PsA for the purpose of enrolling patients in clinical trials and provide guidance to clinicians (4, 5). Therefore, the diagnosis of PsA is mainly performed on clinical features after the exclusion of other seronegative arthritides and no diagnostic tests are available so far.

The pathogenesis of PsA is still poorly understood and both autoinflammation and autoimmunity are believed to play a pivotal role in the disease. Synovial tissue in PsA is characterized by T-cell infiltrate, by marked angiogenesis, and by synovial hyperplasia with increased secretion of cytokines and proteases, which may amplify the local inflammatory process eventually leading to joint destruction (6). Tumor necrosis factor-alpha (TNF-α) is a very important inflammatory mediator and has been implicated in the pathogenesis of articular damage in PsA (6). TNF-α inhibitors are currently used in PsA treatment; however, a high percentage of PsA patients does not respond to TNF-α antagonists (1, 7). Therefore, other cytokines have recently become targets of biological agents, such as interleukin-12 (IL-12), interleukin-23 (IL-23), and interleukin-17 (IL-17) (1, 7). Indeed IL-17 plays a fundamental role in disease development and progression (8).

We have reported the findings of the transcriptome analysis in paired synovial tissue and peripheral blood cells of patients with PsA (9). The upregulation of Th-17 cells related genes and of type I interferon (IFN) inducible genes in PsA patients strengthened the hypothesis that PsA has a strong autoimmune origin, since the coactivity of type I IFN and IL-17 pathways is typical of autoimmunity (9). Moreover, we confirmed these findings with a miRNA microarray analysis in PBMCs of PsA patients showing that pathway enrichment analysis on gene targets of deregulated microRNAs (miRNAs) revealed signaling pathways typically implicated in PsA, such as TNF, mitogen-activated protein (MAP) kinase, and wingless related integration site (WNT) cascades (10).

By this study we wanted to provide a more in-depth knowledge on the epigenetic mechanisms that regulate the PsA pathogenesis by analyzing the expression profiles of long non-coding RNAs (lncRNAs) in the same cohort of patients that we studied in our previous work (10). lncRNAs are important molecules that regulate gene expression through multiple mechanisms and are involved in immune and inflammatory pathways (11). As far as we know, no study has yet taken into consideration lncRNAs expression profiles in PsA patients and only a few data have been reported the deregulation of some lncRNAs in psoriasis (12, 13).

Moreover, in this study we offer a sophisticated and integrated analysis of lncRNAs, miRNAs, and gene expression profiles in PsA patients that allows to identify lncRNAs that regulate transcripts effectively modulated in the disease and that are involved in pathogenetically relevant molecular pathways.

Materials and Methods

Patients

We studied a cohort of 10 patients (6 males and 4 females, mean age: 53.5 years) affected by PsA, attending the Unit of Autoimmune Diseases, at the University Hospital of Verona, Italy. All patients fulfilled the CASPAR criteria for the diagnosis of PsA: inflammatory musculoskeletal involvement combined with at least three features: (1) evidence of current psoriasis, personal history of psoriasis, and family history of psoriasis in unaffected patients; (2) affected nails (onycholysis and pitting); (3) dactylitis; (4) negative rheumatoid factor (RF); and (5) radiographic evidence of new juxta-articular bone formation (excluding osteophytes) (4). All the patients underwent clinical examination and laboratory evaluation comprehensive of inflammatory markers, such as C-reactive protein and erythrocytes sedimentation rate; RF and anti-cyclic citrullinated peptide antibody detected by ELISA test; antinuclear antibody detected by indirect immunofluorescence on HeLa-derived HEp-2 cells; and genetic screening for the association with the allele HLA-B27. All patients underwent the following instrumental investigations: ultrasonography with Power Doppler to investigate subclinical enthesopathy and synovitis in asymptomatic patients, conventional radiography, magnetic resonance imaging, and scintigraphy. The radiological features of peripheral PsA included asymmetric distribution, participation of distal interphalangeal joints, periostitis, bone density preservation, bone ankylosis, and pencil-in-cup deformity.

The patients were affected by cutaneous or nails psoriasis and were enrolled in the study at diagnosis of peripheral PsA before starting immunosuppressive treatment.

All the participants to the study signed a written informed consent and the local Ethical Committee of the University Hospital of Verona, Verona, Italy, had approved the study protocol. All the investigations have been performed according to the principles of the Helsinki declaration.

Microarray Analysis

Blood samples were collected in BD Vacutainer K2EDTA tubes using a 21-gauge needle. PBMCs were obtained upon stratification on Lympholyte® cell separation density gradient (Cedarlane, Burlington, ON, Canada). PBMCs composition was similar between patients and controls. Total RNA extraction from PBMCs was performed with miRNeasy mini kit following manufacturer’s protocol (Qiagen GmbH, Hilden, Germany). 500 ng of total RNA were used for sample preparation starting from 5 ml of blood. cRNA preparation, samples hybridization and scanning were performed following the Affymetrix (Affymetrix, Santa Clara, CA, USA) provided protocols, by Cogentech Affymetrix microarray unit (Campus IFOM IEO, Milan, Italy). All samples were hybridized on Human Clariom D (Thermo Fisher Scientific) gene chip and were analyzed using the Transcriptome Analysis Console 4.0 software (Applied Biosystem, Foster City, CA, USA by Thermo Fisher Scientific, Waltham, MA, USA).

The Human Clariom D arrays allow to interrogate more than 540,000 transcripts sourced from the largest number of public databases starting from as little as 100 pg of total RNA.

The Signal Space Transformation-Robust Multi-Array Average algorithm was applied to background-adjust, normalize, and log-transform signals intensity.

Relative gene expression levels of each transcript were validated applying a One-Way analysis of variance (p ≤ 0.01) and multiple testing correction. Coding genes and lncRNAs that displayed an expression level at least 1.5-fold different in the test sample versus control sample (p ≤ 0.01) were carried forward in the analysis.

The targets (including microRNAs and genes) of all the lncRNAs that satisfied the above-mentioned FC and p-value criteria were screened using NPInter v3.0. This database allows the efficient recovery of all lncRNAs interactions experimentally validated by high-throughput experimental technologies (14, 15).

The list of gene targets of miRNAs that were targeted by significantly modulated lncRNAs was obtained using the FunRich database (16).

Protein–Protein Interaction (PPI) Network Construction and Network Clustering

The Search Tool for the Retrieval of Interacting Genes (STRING version 10.5) was used to obtain PPIs pairs that were validated by experimental studies (17) and to construct the PPI networks. Network topological analysis was performed using the Cytoscape software (18).

High-flow areas (highly connected regions) of the network were detected using the MCODE plugin of Cytoscape, based on the thresholds of k-core = 3 and node score cutoff = 0.2.

Gene Functional Classification and Enrichment Analysis

Genes were functionally classified into canonical biological processes (BPs) on the basis of the gene ontology (GO) annotations.

Biological processes and Pathways enrichment analysis was performed employing FunRich (hypergeometric p-value ≤0.05).

Real-Time PCR

lncRNAs Modulated in PsA

For each sample, 500 ng of total RNA was treated with 1 unit of DNase I Amplification Grade (Invitrogen; Carlsbad, CA, USA) according to the manufacturer’s protocol. First-strand cDNA was generated using the SuperScript IV First-Strand Synthesis System (Invitrogen; Carlsbad, CA, USA) with random hexamers, according to the manufacturer’s protocol. Real-time PCR was performed in triplicate with PowerUp™ Sybr® Green reagent (Applied Biosystems; Foster City, CA, USA) in a QuantStudio 6 Flex system (Applied Biosystems; Foster City, CA, USA). Relative expression levels were calculated for each sample after normalization against the geometric mean of the housekeeping genes GAPDH and beta-actin (ACTB) expression. The ΔΔCt method was used for comparing relative fold expression differences. The data are expressed as fold changes with respect to healthy.

Genes Modulated in PsA

First-strand cDNA was generated using the SuperScript III First-Strand Synthesis System for RT-PCR Kit (Invitrogen), with random hexamers, according to the manufacturer’s protocol. PCR was performed in a total volume of 25 µl containing 1× Taqman Universal PCR Master mix, no AmpErase UNG, and 2.5 µl of cDNA; pre-designed, gene-specific primers, and probe sets for each gene were obtained from Assay-on-Demande Gene Expression Products (Applied Biosystems).

Real-time PCR reactions were carried out in a two-tube system and in singleplex. The real-time amplifications included 10 min at 95°C (AmpliTaq Gold activation), followed by 40 cycles at 95°C for 15 s and at 60°C for 1 min. Thermocycling and signal detection were performed with 7500 Sequence Detector (Applied Biosystems). Signals were detected according to the manufacturer’s instructions and the relative expression levels were calculated as it has been previously described (9).

Detection of Soluble Mediators in Sera of PsA Patients

Serum levels of glypican-4, IFN-γ, Wnt-2, mTOR, TNF-α, SPD-1, NFKB p65, NOTCH1, omentin, and adiponectin were detected using commercially available ELISA kits that were supplied by antibodies-online (glypican-4, Wnt-2, and mTOR), LifeSpan BioSciences (TNF-α, NFKB p65, and omentin), IBL International (IFN-γ), R&D (adiponectin), and Ray Biotech (sPD-1).

Plasmid Construction and Luciferase Reporter Assay

The plasmids to knockdown LINC00909 and LINC00657 were purchased from GenePharma (Shanghai, China).

Human HEK (human embryonic kidney) 293 T cells (1.5 × 104) grown in a 96-well plate were co-transfected with 150 ng of empty vector, 50 ng of firefly luciferase reporter comprising the lncRNAs mut vectors, (Promega, Madison, WI, USA) using Lipofectamie 2000 (Invitrogen, USA). Cells were harvested 48 h after transfection and analyzed using the Dual-Luciferase Reporter Assay System (Promega) according to the manufacturer’s protocol.

Results

High-Throughput Long Non-Coding RNA and Gene Expression Profiling in Peripheral Blood Mononuclear Cells of PsA

In order to evaluate the potential role played by lncRNAs in PsA pathogenesis, we performed a gene array study using the Clariom D human gene chip that enabled us to analyze, at the same time, both conventional gene and lncRNA expression profiles.

We compared the lncRNA expression profiles of PBMC samples obtained from 10 PsA patients with 10 PBMC samples obtained from age and sex matched healthy donors and we found that 259 lncRNAs satisfied the Bonferroni-corrected p value criterion (p ≤ 0.01) and the fold change criterion (FC ≥ |1.5|), displaying robust and statistically significant variation between PsA and healthy controls samples (Table S1 in Supplementary Material). The study was implemented by the analysis of conventional gene expression profiles in the same PsA samples and we found that 1,922 differently expressed genes satisfied the above-mentioned criteria. The complete list of modulated genes can be found in Table S2 in Supplementary Material. In both cases, the arrays were validated by real-time PCR. LncRNAs LUCAT1 and TRIM55-1 were validated by real-time PCR in the entire series of patients analyzed. Significantly different expression levels were found for all tested lncRNAs in PsA as compared to healthy controls (see Figure S1 in Supplementary Material).

Moreover, real-time PCR analysis for seven lncRNAs was carried out in an expanded panel of PsA patients (20 patients) and healthy controls (20 subjects). A significant modulation of all these lncRNAs was found in all tested patients thus confirming gene array results (see Figure 1).

Figure 1

Expression of selected long non-coding RNAs in an expanded panel of psoriatic arthritis patients (20 patients) and healthy controls (20 subjects). Bars indicate SD.

To gain meaningful insights on the potential role played by modulated lncRNAs in PsA pathogenesis, the complete list of modulated lncRNA was filtered, extracting only those transcripts for which a bona fide target annotations was present in NPInter. By this method 92 lncRNAs were selected (Table S3 in Supplementary Material) and, simultaneously, the list of all gene and microRNA targets of the selected lncRNAs, experimentally validated by high-throughput technologies, was extracted from the same database.

To corroborate our results we narrowed down our analysis to modulate lncRNA that targeted genes that were significantly modulated in the array and miRNAs that we found deregulated in our previous analysis of PsA PBMCs from the same cohort of patients (10) (Table 1). In particular, we found that 15 of these miRNAs modulated in PsA (hsa-miR-130a-3p, hsa-miR-148a-3p, hsa-miR-151a-5p, hsa-miR-17-5p, hsa-miR-186-5p, hsa-miR-199a-3p, hsa-miR-199a-5p, hsa-miR-28-5p, hsa-miR-3135b, hsa-miR-320c, hsa-miR-320d, hsa-miR-331-3p, hsa-miR-423-5p, hsa-miR-451a, and hsa-miR-92a-3p) were targeted by selected lncRNAs. We then extracted from the FunRich database the annotated gene targets of the above-mentioned miRNAs selecting only transcripts that also resulted when modulated in the Clariom D array (Table S4 in Supplementary Material). Table 2 recapitulates the above-selected lncRNAs and targets.

Table 1

Selected modulated long non-coding RNAs in psoriatic arthritis patients versus healthy controls.

ID	Fold change	p-value	Gene symbol	mRNA accession
TC0500008318.hg.1	−2.37	0.0026	EPB41L4A-AS1	ENST00000413221.2
TC0700007000.hg.1	−2.02	0.0015	HOTAIRM1	ENST00000616712
TC0600008510.hg.1	−2.25	0.0087	KCNQ5-IT1	ENST00000445310
TC0700013567.hg.1	2.05	0.0013	LINC00174	ENST00000416366
TC0700007277.hg.1	2.46	0.0001	LINC00265	ENST00000340510.4
TC1500007707.hg.1	−1.58	0.0089	LINC00593	ENST00000558385.1
TC2000008995.hg.1	−2.27	0.0017	LINC00657	ENST00000565493
TC1800009043.hg.1	−1.84	0.0025	LINC00909	ENST00000577806
TC0200007199.hg.1	2.27	0.0007	LINC00486	ENST00000414054
TC0100009691.hg.1	−1.97	0.0003	RP11-403I13.5	ENST00000443018.1
TC0200010127.hg.1	−3.24	<0.0001	RP11-171I2.4	ENST00000605334.1
TC0200011420.hg.1	−2.68	0.0005	AC133528.2	ENST00000433036.1
TC0400009914.hg.1	−2.37	0.0027	RP11-539L10.3	ENST00000513179.1
TC0500009465.hg.1	−1.74	0.007	RP11-779O18.3	ENST00000523005.1
TC0800007847.hg.1	−14.9	<0.0001	AC084082.3	ENST00000517961.2
TC1100011278.hg.1	1.91	0.0003	RP11-867G23.3	ENST00000501708.1
TC1200006772.hg.1	−1.86	0.0044	RP11-75L1.1	ENST00000541404.1
TC1200010732.hg.1	−3.11	0.0058	RP11-1100L3.8	ENST00000564363.1
TC1400006719.hg.1	−2.54	0.0066	RP11-468E2.5	ENST00000558478.1
TC1400009275.hg.1	−1.61	0.0027	RP11-930O11.2	ENST00000560296.1
TC1600009188.hg.1	2.21	0.0021	LA16c-360H6.3	ENST00000574245.1
TC1700007241.hg.1	2.64	0.0008	RP11-283C24.1	ENST00000578585.1
TC2100007843.hg.1	−1.78	0.0016	AF131217.1	ENST00000430247.1
TC2200008462.hg.1	−3.38	0.0058	RP3-430N8.10	ENST00000602955.1
TC1500010312.hg.1	−2.52	0.0064	RP11-815J21.2	ENST00000561409.1
TC1800007426.hg.1	−3.49	0.0032	RP11-1151B14.4	ENST00000591360.1
TC1900011833.hg.1	1.81	0.0031	CTB-25B13.12	ENST00000588225.1
TC1900007159.hg.1	−1.77	0.0084	CTB-55O6.10	ENST00000590715.1
TC1200008393.hg.1	−1.74	0.0028	RP11-981P6.1	ENST00000552778.1
TC1200008425.hg.1	−2.1	0.0024	RP11-796E2.4	ENST00000499685.2
TC1400009962.hg.1	−2.61	0.0006	RP11-471B22.2	ENST00000555853.1
TC1600006833.hg.1	2.16	0.0024	RP11-77H9.5	ENST00000564919.1
TC1400009667.hg.1	1.88	0.0092	RP4-693M11.3	ENST00000557304.1
TC1000009009.hg.1	−1.6	0.0098	RP11-498B4.5	ENST00000433600.1
TC1400010386.hg.1	1.71	0.0064	CTD-3051D23.4	ENST00000553344.2
TC1200008527.hg.1	−2.91	0.0065	RP11-256L6.3	ENST00000551849.1
TC0200007485.hg.1	1.61	0.0048	AC016722.4	ENST00000429761.1
TC1400007302.hg.1	1.95	0.0041	CTD-2002H8.2	ENST00000557322.1
TC0500008150.hg.1	−2.4	0.0066	CTD-2260A17.1	ENST00000512856.1
TC0600010636.hg.1	−1.89	0.0005	RP3-406P24.3	ENST00000415144.1

Table 2

Long non-coding RNAs (LncRNAs) and their targets modulated in psoriatic arthritis.

LncRNAs	Gene targets	miRNA targets
LINC00174	CPSF7
LINC00265	IGF2BP2
LINC00593	UPF1
LINC00657	CPSF7, FXR1, HNRNPC, NUDT21 UPF1, ZC3H7B	hsa-miR-130a-3p
		hsa-miR-17-5p
		hsa-miR-186-5p
		hsa-miR-199a-3p
		hsa-miR-199a-5p
		hsa-miR-28-5p
		hsa-miR-320c
		hsa-miR-320d
		hsa-miR-331-3p
		hsa-miR-423-5p
		hsa-miR-451a
LINC00909	UPF1	hsa-miR-130a-3p
		hsa-miR-148a-3p
		hsa-miR-28-5p
		hsa-miR-320c
		hsa-miR-320d
LINC00486	FXR1, UPF1
EPB41L4A-AS1	CPSF7, UPF1	hsa-miR-130a-3p
		hsa-miR-17-5p
HOTAIRM1	CPSF7, UPF1
KCNQ5-IT1	CPSF7, UPF1
RP11-403I13.5		hsa-miR-92a-3p
RP11-171I2.4	CPSF7
AC133528.2	CPSF7
RP11-539L10.3		hsa-miR-148a-3p
		hsa-miR-199a-3p
RP11-779O18.3	HNRNPC, UPF1, ZC3H7B
AC084082.3	CPSF7, NUDT21, UPF1	hsa-miR-331-3p
RP11-867G23.3	UPF1, IGF2BP2
RP11-75L1.1	UPF1
RP11-1100L3.8	UPF1	hsa-miR-17-5p
		hsa-miR-423-5p
RP11-468E2.5	UPF1
RP11-930O11.2	UPF1
LA16c-360H6.3	IGF2BP2	hsa-miR-17-5p
RP11-283C24.1	UPF1
AF131217.1	UPF1
RP3-430N8.10		hsa-miR-331-3p
RP11-815J21.2		hsa-miR-186-5p
RP11-1151B14.4	UPF1
CTB-25B13.12	UPF1
CTB-55O6.10	CPSF7, UPF1
RP11-981P6.1	CPSF7, UPF1, ZC3H7B
RP11-796E2.4	UPF1	hsa-miR-3135b
RP11-471B22.2	CPSF7, IGF2BP2, UPF1, ZC3H7B
RP11-77H9.5	CPSF7, HNRNPC, IGF2BP2, UPF1
RP4-693M11.3	UPF1
RP11-498B4.5	UPF1
CTD-3051D23.4	IGF2BP2
RP11-256L6.3	NUDT21, UPF1
AC016722.4	CPSF7, UPF1, ZC3H7B
CTD-2002H8.2	CPSF7, HNRNPC, NUDT21, UPF1, ZC3H7B
CTD-2260A17.1		hsa-miR-151a-5p
RP3-406P24.3	CPSF7

In summary, we selected only those lncRNAs, that targeted miRNAs and genes with evidence of modulation in our PsA samples to trace, with a good confidence, lncRNAs–miRNAs–genes interactions that are expected to be established in the course of PsA.

To validate our results we conducted a LINC00909siRNA and a LINC00657siRNA silencing in human 293 T cells, to explore whether this knockdown altered the expression levels of selected miRNAs targeted by LINC00909 (miR-148a-3p and mir-28-5p) and by LINC00657 (miR-130a-3p and miR-17-5p). We observed that the silencing of the two lncRNAs significantly increased the level of expression of their targeted miRNAs.

The registered percentages of increase were: 70 ± 3.8 and 78 ± 1.4% for miR-148a-3p and mir-28-5p, respectively, and 75 ± 2.3 and 80 ± 1.53% for miR-130a-3p and miR-17-5p, respectively.

PPI Network of Modulated lncRNAs in PsA

The lncRNA expression profiling of PsA PBMCs was integrated with a network analysis. We, therefore, inspected, by a bioinformatic methodology, all the functional and experimentally validated interactions among the protein products of genes targeted by lncRNAs and by the 15 above-mentioned miRNAs (230 genes) that were selected as previously described. Then, taking in account of these interactions, we constructed a PPI network which showed a good PPI enrichment p-value (0.00041). In the obtained network 229 of the above-mentioned genes (nodes), were connected by 195 pairs of interactions (edges). Since 229 out of 230 of lncRNAs and miRNAs gene targets that we used as input for the network analysis, resulted connected in the PPI network, we could observe that the selected lncRNAs may act in an integrated fashion in the disease.

We then imported in Cytoscape the PPI network, adding to its scaffold the lncRNAs-genes, miRNAs-genes, and lncRNAs–miRNAs interactions that we selected as we explained above. Thus, we obtained an implemented network that is showed in Figure 2. A topological analysis of this network was performed to highlight lncRNAs that were highly connected to the PPI network genes (i.e., targeting a large number of PPI-network genes) also considering all the lncRNA–gene interactions mediated by the selected miRNAs (i.e., lncRNA–miRNA–gene interactions). By this approach we could select seven lncRNAs, including EPB41L4A-AS1, LA16c-360H6.3, LINC00657, LINC00909, RP11-1100L3.8, RP11-539L10.3, and RP11-403I13.5 that displayed the highest connectivity in the network. As shown in Figures 3A–G (where genes of the PPI network are ordered around a circle based on their degree of connectivity, i.e., number of edges) the above-mentioned lncRNAs directly or indirectly targeted genes that, in most cases, were distributed in areas of the network characterized by a high degree of connectivity. This observation indicated that, since these lncRNAs targeted highly interconnected genes, they may have a broad effect in the disease that goes beyond the modulation of single target gene(s).

Figure 2

Protein–protein interaction (PPI)-network of differently expressed genes in psoriatic arthritis (PsA) targeted by the selected long non-coding RNAs and microRNAs modulated in PsA.

Figure 3

Protein–protein interaction (PPI) network of differently expressed genes targeted by selected microRNAs and long non-coding RNAs (lncRNAs). Graphical representation of the connections among genes in the PPI network and EPB41L4A-AS1 (A), LA16c-360H6.3 (B), LINC00657 (C), LINC00909 (D), RP11-1100L3.8 (E), RP11-539L10.3 (F), and RP11-403I13.5 (G). In the PPI network all genes are ordered around a circle based on their degree of connectivity (i.e., number of edges). Genes directly or indirectly targeted by the lncRNAs are highlighted in yellow.

Modular Analysis of Genes and lncRNAs Modulated in PsA

In a second part of our analysis we wanted to dissect the global impact of lncRNAs on the entire PsA transcriptome (i.e., on the 1,922 modulated genes in PsA). Since the targeting of genes with more interacting partners may have a broader impact on the global transciptome than the modulation of few isolated genes, we aimed to inspect if the modulated lncRNAs may display connections with highly interacting genes in PsA. To this purpose, we built the PPI network that included the protein products of all modulated genes in PsA that showed experimentally validated interactions obtaining a wide network containing 1,637 nodes and 9,415 edges and showing a very good PPI enrichment p-value (p < 10–16). Then we performed a modular analysis to find areas of the network that included the most highly connected genes (modules). By this analysis we could extract 16 modules (M1–M16, see Table S5 in Supplementary Material). All the 16 modules were imported in Cytoscape along with all the previously selected interactions with miRNAs and modulated lncRNAs that we also considered in the network analysis described in the previous paragraph. The graphical visualization of genes connected in modules and their interactions with lncRNAs showed that the majority of the 16 modules (10/16) were directly or indirectly targeted by modulated lncRNAs, including M1, M2, M4, M5, M8, M10, M11, M12, M14, and M15 (Figure 4). Moreover, we could observe that the vast majority of modules was targeted by the above-mentioned lncRNAs that included LINC00909 (targeting 9/10 modules) LINC00657 (targeting 8/10 modules), EPB41L4A-AS1, RP11-539L10.3 (targeting 7/10 modules), LA16c-360H6.3, RP11-403I13.5, and RP11-1100L3.8 (targeting 5/10 modules) (Table S6 in Supplementary Material).

Figure 4

Gene modules that are targeted by the selected long non-coding RNAs (lncRNAs) modulated in psoriatic arthritis. In each module (M) genes are represented by light blue rectangles whereas targeting microRNAs and lncRNAs are represented by orange and yellow rectangles, respectively.

Functional Analysis of the Targeted Modules and Their Targeting lncRNAs

To dissect the possible role played by the module-targeting-lncRNAs we performed a functional classification of genes included in the 10 targeted modules, based on the GO classification criteria. Figure 5 recapitulates all the most relevant GO BPs in which the genes of targeted modules are classified.

Figure 5

Functional classification of genes included in the targeted modules. Graphical representation of long non-coding RNAs targeting the gene modules and of the most relevant biological processes (BPs) in which genes in each module can be classified. The BPs are represented by colored circles (each color corresponds to a different BP). The size of circles does not correspond to the number of genes that are classified in each BP.

M1 was the most targeted module and 35 lncRNAs showed direct or indirect connections with it (see Table S6 in Supplementary Material). The 88 genes included in this module were mostly related to translation, mRNA splicing, and protein ubiquitination GO BPs. Among genes included in protein ubiquitination, nine genes belonged to the “proteasome-mediated ubiquitin-dependent protein catabolic process” BP, including CUL1, CUL3, FBXW5, FBXW7, FBXL19, HECTD3, NEDD4L, RNF126, and UBE2A. Interestingly, it has been recently observed that ubiquitin–proteasome-mediated protein degradation plays a role in modulating the balance between bone formation and bone resorption, since it is involved in signaling pathways that are crucial for bone homeostasis like RANK/NF-κB pathway and the Wnt/β-catenin pathway (19). In addition CUL3 is also associated with the “Wnt signaling pathway” BP. 12 genes were targeted in the M1, including CUL5, TRA2B, SRRM1, ASB6, HNRNPUL1, CPSF7, HNRNPC, NUDT21, UPF1, and the above-mentioned FBXL19, NEDD4L, and FBXW7. UPF1 was the most targeted gene in M1. It is critical for the activity of Regnase-1 also called MCPIP1 (monocyte chemotactic protein-1-induced protein-1), a crucial molecule that regulates immune response avoiding over-activation of the immune system (20). Interestingly, the gene FBXL19 represents a susceptibility locus for PsA (21). Finally, other two genes of M1, TRIM32, and HERC5 play a role in the innate immune response and, in particular, HERC5 activates IRF3 (22).

M2 was targeted by 6 lncRNAs and included 40 genes that were mainly involved in cell division, inflammatory response (i.e., CCR7, CXCR1, CXCR2, PF4, and PPBP/CXCL7), tumor necrosis factor (TNF)-mediated signaling pathway (i.e., TNFRSF1B, TNFSF12, TNFRSF12A, and PSMA7), and Wnt signaling pathway (i.e., AXIN1, DVL3, GNB2, PPP2R5C, PSMA4, PSMB1, PSMC4, PSMD4, and RNF146) BPs. Interestingly, the expression of CXCL7 is increased in Rheumatoid arthritis (AR) synovia (23). Moreover AXIN2 another gene included in M2 is related to the intramembranous ossification BP. Eight genes were targeted in M2, including ADCY3, GNAI1, PSME3, PSMF1, RELA, TRAF3, and the above-mentioned TNFRSF1B and CCR7. In particular, the expression of CCR7 is stimulated by TNF-α and interleukin 1-beta in osteoclasts (24). Moreover, RELA is involved in peripheral regulatory T cell-induced tolerance (25) and TRAF3 is required for T cell and invariant natural killer T cell effector functions (26). In addition, TRAF3 is also involved in limiting osteoclastogenesis, indeed RANKL-induced degradation of TRAF3 enhances TNF-induced osteoclastogenesis (27).

M4 was targeted by seven lncRNAs and included 57 genes that were mainly involved in cell division, DNA repair, endocytosis, mitochondrial translation, mRNA export from nucleus, regulation of glycolytic process, and transcription. In addition, this module included two genes (i.e., MRE11A and XRCC6) associated with the BP “positive regulation of type I interferon production” and the gene POU2F2/OCT2 that increases the differentiation of antibodies secreting activated B cells (28). In this module the genes GBF1, MYBL2, RAD23A, and WASL were targeted.

M5 was targeted by the RP11-403I13.5 lncRNA and included 19 genes. Protein targeting to ER and translation were the most represented BPs in this module.

M8 was targeted by four lncRNAs (EPB41L4A-AS1, LINC00657, LINC00909, and RP11-539L10.3) and in this module seven genes were connected, including BRPF3, TAF12, TAF2, TAF7, TBPL1, TRIM24, and BPTF almost associated with the “translation” BP. Among these, BPTF is important for the homeostasis of T cells and is crucial for the maintenance and function of regulatory T (Treg) cells (29).

M10 was targeted by 8 lncRNAs and included 49 genes. These genes were associated with many BPs, including cell proliferation, signal transduction, transcription, Wnt signaling pathway (i.e., BCL9L, FZD6, LRP6, ATP6V0C, BCL9, and PYGO2), vasculogenesis (i.e., ENG and TGFBR2), bone resorption (i.e., TCIRG1), negative regulation of bone mineralization (i.e., SRGN), osteoblast differentiation (i.e., LGR4), and osteoblast development (i.e., LRP5). Several genes were associated with BPs related to glucose homeostasis like, for example, PPARGC1A (cellular glucose homeostasis), GSK3A (cellular response to insulin), IGF1R (insulin receptor signaling pathway), RAF1 (insulin secretion), and RAP1A (positive regulation of glucose import). In addition other genes were classified into BPs related to adaptive immune response like T cell differentiation (i.e., RIT1 and TCF7), and T cell cytokine production (KDELR1) and B cell homeostasis (SOS1). In M10 12 genes were targeted, including SHOC2, CCNL1, CDK13, MED16, MMS19, and the above-mentioned FZD6, PYGO2, KDELR1, PPARGC1A, RAP1A, SOS1, and TGFBR2.

M11 was targeted by 2 lncRNAs and included 28 genes. A large number of these genes were assigned to “lipid metabolic process” (i.e., PLCD3, PPAP2B, SLC27A1, and TBL1X), “phospholipid metabolic process” (i.e., AGPAT2, LPCAT3, LPCAT4, and PPAP2C), and to “phosphatidylinositol biosynthetic process” (CDIPT, MTMR14, PI4K2B, and PIK3C2B) BPs. Interestingly, two genes namely CARM1 and FASN were associated with “endochondral bone morphogenesis” and “osteoblast differentiation” BPs, respectively. In this module, TERF2, YPEL3, and the above-mentioned MTMR14 were targeted.

M12 was targeted by 6 lncRNAs and included 30 genes. We observed that several genes of this module classified in BPs related to the immune response, including SPI1/PU.1 (lymphoid progenitor cell differentiation), CD244 and KLRF1 (innate immune response), CD226 (positive regulation of natural killer cell cytokine production), CD7 (T cell activation), CD247 (T cell receptor signaling pathway), GBP1 (regulation of T cell receptor signaling pathway), CD79A (B cell activation), GBP1, GBP2, and CIITA (type I interferon signaling pathway), HCST, LILRB2, and KLRB1 (regulation of immune response), GZMA (immune response) and SPI1/PU.1 (lymphoid progenitor cell differentiation), and CIITA (positive regulation of MHC class I and of MHC class II biosynthetic process). In particular, the transcription factor SPI1/PU.1 is crucial for the development of interleukin-9-producing helper T cells (Th9) (30, 31) that have been recently involved in the PsA pathogenesis (32).

Moreover five genes, including DAG1, COL14A1, LAMA2, LAMA3, and LAMC1 were associated with the “extracellular matrix organization” BP. LAMA3 and FAM168B were targeted in this module.

M14 was targeted by three lncRNAs and included 14 genes that were mainly associated with protein targeting to mitochondrion (TOMM20 and TOMM40), protein import into mitochondrial inner membrane (TIMM22), ATP biosynthetic process (ATP5I), apoptosis (BCL2L11 and CASP4), and in Wnt signaling pathway (TLE1 and TLE3) BPs. In this module the genes BCL2L11 and VDAC1 were targeted.

Finally, M15 was targeted by 6 lncRNAs and included 44 genes. In this module, we observed the presence of genes classified in several BPS, including cell division (TUBA1C, TFDP1, CEP192, CEP63, and MAPRE1) glycolytic process (PFKL, PFKM, and PKM), Wnt signaling pathway (MARK2, CTNNBIP1, EIF2C1, and TNRC6C), osteoblast differentiation (ASF1A and MAPK11), and regulation of type I interferon-mediated signaling pathway (CDC37). Genes targeted in this module were: AHNAK, CALB1, CCT6A, CEP63, CEP97, GAB1, MAPRE1, PFKM, and TNRC6C.

CALB1 plays a role in differentiation and matrix formation in osteoblasts (33), GAB1 is important for normal postnatal bone homeostasis (34) and MAPRE1 regulates cell–cell adhesion-induced osteoblast differentiation (35).

Since disease can be viewed as the result of perturbations of complex signaling networks, a pathway enrichment analysis was performed, in each module, using FUnrich. This analysis showed that the most represented enriched pathways were related to immune response (including signaling of both adaptive and innate immune response, type I interferon and gamma interferon signaling), bone homeostasis, metabolism, and cell adhesion.

The results obtained are shown in Figure 6, which shows a schematic representation of the most relevant enriched signaling pathways (p ≤ 0.05).

Figure 6

Biological processes and pathways enrichment analysis of the targeted modules. Selection of the different enriched pathways (p < 0.05) in the 10 targeted modules (M). Pathways that are enriched in more than one module are labeled by the same color in the different modules.

To provide experimental evidence that pathways identified by the throughput analysis of modulated genes and lncRNAs in PsA are active in the disease, we performed a two step validation of the presented results. Indeed, first we validated by real-time PCR the modulation of genes ascribed to several signaling pathways, including SH3GL1, EPN1, and NUP214 (mTOR signaling), TRAF2, TNFRSF1B, and RELA (TNF/NF-κB signaling), MAPK7 and MAPK11 (toll-like receptors signaling), TCF7, BCL9, and LRP6 (Wnt signaling), MED15, ENG, and ACTN4 (glypican signaling), GSK3A, IGF1R, and PARGC1A (insulin signaling), and CIITA IFI27 and MYD88 (type I interferon signaling). All the tested transcripts were overexpressed in PsA patients when compared to healthy subjects (see Figure S2 in Supplementary Material).

Second, in the sera of patients with PsA, we evaluated the presence of several molecules associated with the pathways modulated in PsA. Indeed, we decided to quantify the serum levels of glypican-4, IFN-γ, Wnt-2, mTOR, TNF-α, sPD-1, NFKB p65, NOTCH1, omentin, and adiponectin. Figure 7 shows the concentration of these molecules in the sera of the PsA patients. The serum levels of all the tested molecules were significantly increased in PsA patients when compared to healthy donors with the exception of adiponectin and omentin that were significantly decreased in PsA sera.

Figure 7

Serum levels of selected molecules in psoriatic arthritis (PsA) patients and in normal subjects. The histograms represent the mean of the results obtained in 20 healthy donors and in 20 PsA patients.

Discussion

Psoriatic arthritis is a chronic inflammatory arthritis that affects 10–30% of patients with skin psoriasis. IL-17, TNF-α, and type I IFNs play a fundamental role in the pathogenesis of the disease and monoclonal antibodies targeting IL-17 and TNF-α along with IL-12/IL-23 are the main biotechnological treatments that have a dramatic good response in PsA and in skin psoriasis. As we have previously reported (9), the involvement of the mentioned cytokines is typical of an autoimmune process and support the hypothesis of the autoimmune origin of the disease. In particular, IL-17 is able to synergize with TNF-α, IL-22, and other cytokines, including IL-6 and IL-8, in sustaining the inflammatory process at different sites and in favoring the development of comorbidities that are typical of psoriasis, such as PsA, metabolic syndrome, obesity, and cardiovascular disease.

Long non-coding RNAs control gene expression at multiple levels, including epigenetic regulation, chromatin remodeling, and post-transcriptional gene regulation. Accumulating evidences indicate that lncRNAs can be involved in different types of immune-mediated human diseases, including autoimmune diseases (11). The role played by lncRNAs in PSA pathogenesis has not been elucidated yet and a comprehensive analysis of lncRNA expression profiles in PsA has not been performed. We, therefore, aimed to identify lncRNA expression signatures associated with PsA, analyzing the expression of a vast array of lncRNA in PBMC cells from PsA patients by microarray analysis.

The applied criteria for the selection of modulated lncRNAs included a multiple step process that combined the lncRNAs expression analysis to the conventional gene expression and miRNAs profiling in the same cohort of patients.

By this approach we were able to select deregulated lncRNAs that were as much as possible correlated to the PsA transcriptome.

Through a complex network analysis of the whole set of molecular interactions among modulated genes and lncRNAs we further wanted to select those lncRNAs that targeted modules of most highly connected genes in the PsA interactome, since they may have a major impact on PsA gene modulation. That, for example, has led us to favor lncRNAs with high connectivity rather than poorly connected lncRNAs (even) with higher fold change. The functional analysis of the above-mentioned highly connected modules revealed that the modulated lncRNAs targeted genes involved in BPs that play a pivotal role in the PsA pathogenesis such as, for example, immune response (including B and T cell activation) inflammatory response, TNF, Wnt and type I interferon-mediated signaling, bone resorption, bone mineralization, and glyco-lipid metabolic process (see Results), thus indicating that the selected lncRNAs regulated the three major aspects of the disease, including skin involvement, osteo-articular features, and metabolic syndrome.

As it has recently acquired in molecular biology, diseases are now viewed in the context of signaling pathways perturbation. We, therefore, performed a pathway enrichment analysis of all targeted gene modules and, by this approach we reached a non-ambiguous identification of important pathways that can be modulated by the selected lncRNAs. Our analysis confirmed that these lncRNAs play a role in the regulation of signaling pathways associated with both adaptive and immune response and with autoimmune diseases (i.e., type I interferon signaling pathway). In addition, they regulate important signaling involved in bone homeostasis including glypican, syndecan, and Wnt pathways.

Moreover, these results revealed the presence in PsA of lncRNAs that modulate the signaling pathway associated with influenza, thus highlighting the possible role played by influenza viruses in triggering PsA and suggesting the possible involvement of an epigenetic control. The selected lncRNAs were involved in the modulation of lipid metabolism pathway (see Results, module M11), and it is well known that dyslipidemia is frequently associated with PsA (36) and, in particular, patients with PsA have been shown to have elevated cholesterol levels (37) and increase risk of cardiovascular disease.

Another deregulated signaling pathway is the S6K1 (S6 kinase 1) and mTOR (mammalian target of rapamycin) pathway (see Results, modules M1, M4, and M5), a nutrient sensing system that respond to nutrient overload leading to insulin resistance, obesity, and diabetes (38, 39). These observations suggest that lncRNAs may also be involved in the metabolic syndrome frequently associated with PsA. Hyperuricemia is another feature of metabolic syndrome and is also known to be more frequent in psoriatic subjects than in the normal population. In humans, purines are metabolized into uric acid, which is a strong stimulator of innate immunity. In this regard, hyperuricemia may increase uric acid crystallization in and around joints, thereby inducing PsA in psoriatic subjects, thus representing an independent risk factor for PsA. Activation of pathways modulation of genes involved in related to nucleobase, nucleoside, nucleotide, and nucleic acid metabolism (modules M1, M4, and M8) may reflect this particular aspect of PsA.

Finally, lncRNAs control several pathways involved in cell adhesion (see Results, modules M4, M10, M11, M12, M14, and M15). The modulation of these pathways may contribute to the increased mucosal permeability of the gut that has been observed in PsA, even in the absence of bowel pathology (40). Indeed, bowel disease is six times more frequent in patients with PsA. Since an altered intestinal barrier is a feature of PsA, the presence of gastrointestinal involvement in PsA should be carefully investigated during the diagnostic work up.

The detection of significantly different serum levels in PsA patients compared to healthy subjects of several soluble mediators, ascribed to the main pathways targeted by the selected lncRNAs, strongly suggests that they modulate signaling networks that are active in the disease. In particular, the increased serum levels of glypican-4 and the decrease of adiponectin and omentin highlight the typical metabolic aspects of PsA patients.

Indeed glypican-4 has been considered a novel adipokine that boosts insulin signaling and, interestingly glypican-4 levels are significantly increased in obese patients with insulin resistance (41). Moreover, both adiponectin and omentin have been showed decreased in sera from psoriatic patients (42, 43) and decreased levels of both molecules have been associated with, and involved in, metabolic syndrome (44, 45).

In conclusion, this study is the first report on lncRNAs that may exert an epigenetic control on PsA pathogenesis. The original approach that integrates genes and lncRNAs expression profiles presented in the study, allows to identify deregulated lncRNAs that modulate crucial molecular signaling associated with the typical features of the disease. These findings further support that PsA is an autoimmune disease with systemic inflammation associated with obesity and metabolic syndrome leading to cardiovascular disease.

Finally, we may suggest that these lncRNAs may be useful to design novel therapeutic strategies in PsA.

Ethics Statement

All the participants to the study signed a written informed consent and the local Ethical Committee of the University Hospital of Verona, Verona, Italy, had approved the study protocol. All the investigations have been performed according to the principles of the Helsinki declaration.

Author Contributions

APU, CL, and MD conceived and designed the experiments. MD, APE, PF, GP, and ET performed the experiments. MD and APE analyzed the data. GP and ET contributed reagents and collected the patients’ samples. MD wrote the paper with input from APU and CL.

Conflict of Interest Statement

The research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.