Long Non-Coding RNAs Target Pathogenetically Relevant Genes and Pathways in Rheumatoid Arthritis

Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease driven by genetic, environmental and epigenetic factors. Long non-coding RNAs (LncRNAs) are a key component of the epigenetic mechanisms and are known to be involved in the development of autoimmune diseases. In this work we aimed to identify significantly differentially expressed LncRNAs (DE-LncRNAs) that are functionally connected to modulated genes strictly associated with RA. In total, 542,500 transcripts have been profiled in peripheral blood mononuclear cells (PBMCs) from four patients with early onset RA prior any treatment and four healthy donors using Clariom D arrays. Results were confirmed by real-time PCR in 20 patients and 20 controls. Six DE-LncRNAs target experimentally validated miRNAs able to regulate differentially expressed genes (DEGs) in RA; among them, only FTX, HNRNPU-AS1 and RP11-498C9.15 targeted a large number of DEGs. Most importantly, RP11-498C9.15 targeted the largest number of signalling pathways that were found to be enriched by the global amount of RA-DEGs and that have already been associated with RA and RA–synoviocytes. Moreover, RP11-498C9.15 targeted the most highly connected genes in the RA interactome, thus suggesting its involvement in crucial gene regulation. These results indicate that, by modulating both microRNAs and gene expression, RP11-498C9.15 may play a pivotal role in RA pathogenesis.


Introduction
Rheumatoid arthritis (RA) is an autoimmune disease characterized by chronic inflammation of the joints with severe pain and swelling, joint damage and disability, which ultimately leads to joint destruction and loss of function [1]. Several genome-wide association studies have identified genetic variants that confer RA risk. However, these variants can explain less than 20% of susceptibility in RA. Several factors have been shown to contribute to the onset of the disease, such as genetic susceptibility, environmental factors including smoking, and epigenetic mechanisms [2].
LncRNAs are epigenetic regulators of gene expression and are involved in immune and inflammatory molecular networks. Moreover, it has been demonstrated that they also play a role in several autoimmune diseases [3].
The modulation of some lncRNAs has been described in RA [4], but the interplay between lncRNAs and gene modulation in RA has only been partially explored. Indeed, only a few microarray studies described the gene expression profiling analysis of a congruous number of lncRNAs and mRNA transcripts in RA PBMCs [5,6], but in these works, only target prediction methods were used to inspect all the possible interactions among lncRNAs and coding genes. Moreover, only lncRNA-mRNA pairs v2.0 (http://starbase.sysu.edu.cn/starbase2/index.php), where lncRNAs interactions, experimentally validated by high-throughput experimental technologies, are registered [8].

Protein-Protein Interaction (PPI) Network Construction and Network Clustering
The PPI network was constructed upon the experimentally validated protein-protein interactions using STRING (Search Tool for the Retrieval of Interacting Genes) version 10.5 (http://string-db.org/) [10]. Network topological analysis was performed using Cytoscape software (http://www.cytoscape.org/) [11]. High-flow areas (highly connected regions) of the network (modules) were detected using the MCODE plugin of Cytoscape (k-core = 4 and node score cutoff = 0.2).
Pathway classification and enrichment (Bonferroni corrected p-value ≤ 0.05) analysis were achieved with FunRich.

Real-Time PCR of LncRNA
First, 500 ng of total RNA were treated with one unit of DNase I Amplification Grade (Invitrogen, Carlsbad, CA, USA). First-strand cDNA was generated using the SuperScript IV First-Strand Synthesis System (Invitrogen) with random hexamers, according to the manufacturer's protocol. Real-time PCR was performed in triplicate with a PowerUp™ Sybr ® Green reagent (Applied Biosystems) in a QuantStudio 6 Flex system (Applied Biosystems). Transcripts' relative expression levels were obtained 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. Results are expressed as fold changes with respect to healthy patients.

Real-Time PCR of Genes Modulated in RA Patients
First-strand cDNA was obtained using the SuperScript III First-Strand Synthesis System for RT-PCR Kit (Invitrogen), with random hexamers, following 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-Demand Gene Expression Products service (Applied Biosystems).
Real-time PCR reactions were carried out in a two-tube system and in singleplex. The real-time amplifications encompassed 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 a 7500 Sequence Detector (Applied Biosystems). Signals were detected by following the manufacturer's instructions. This methodology allows for the identification of the cycling point where the PCR product is detectable by means of fluorescence emission (threshold cycle or Ct value). The Ct value correlates to the quantity of target mRNA. Relative expression levels were calculated for each sample after normalization against the housekeeping genes GAPDH, beta-actin and 18s ribosomal RNA (rRNA), using the ∆∆Ct method for comparing relative fold expression differences. Ct values for each reaction were determined using TaqMan SDS analysis software (Applied Biosystems). For each amount of RNA tested, triplicate Ct values were averaged. Since Ct values vary linearly with the logarithm of the amount of RNA, this average represents a geometric mean.

Real-Time PCR of MicroRNA
miRNA expression was evaluated by TaqMan ® Advanced miRNA assays chemistry (Applied Biosystems). Briefly, 10 ng of total RNA was reverse transcribed and pre-amplified with TaqMan ® Advanced miRNA cDNA synthesis kit according to the manufacturer's instructions (Applied Biosystems). Pre-amplified cDNA was diluted 1/10 in nuclease-free water and 5 µL of diluted cDNA for each replicate were loaded in PCR. 20 µL PCR reactions were composed by 2× Fast Advanced Master Mix and TaqMan ® Advanced miRNA assay for miR-520e. The mean of Ct for hsa-miR-16-5p and hsa-miR-26a-5p expression was used to normalize miRNA expression. Real-time PCR was carried out in triplicate on a QuantStudio 6 Flex instrument (Applied Biosystems). Expression values were reported as fold change with respect to healthy controls by the ∆∆Ct method, employing QuantStudio Real-Time PCR system software v. 1.3.

Statistical Analysis
Statistical testing was performed using SPSS Statistics 2 software (IBM, Armonk, NY, USA). Data obtained from RT-PCR analysis of RA samples and healthy controls were analysed using the Mann-Whitney Test.

High-Throughput Gene and Long Non-Coding RNA Expression Profiling in Peripheral Blood Mononuclear Cells of RA Patients
We simultaneously profiled the expression of more than 540,000 human transcripts, including those ascribed to more than 50,000 long non-coding RNAs (lncRNAs), in four PBMC samples from patients with clinically diagnosed RA symptoms, with the purpose of identifying lncRNAs potentially involved in RA pathogenesis. RA-associated transcriptional profiles were compared to those obtained from four age-and sex-matched healthy subjects and, 97 lncRNAs and 942 coding genes were selected applying a robust filtering approach (FDR-corrected p-value ≤ 0.05 and fold change ≥ |1.5|) (Tables S1 and S2).
The functional classification by Gene Ontology (http://www.geneontology.org/) of the 942 differentially expressed genes (DEGs) highlighted the modulation of transcripts that play a role in biological processes (BPs) strictly associated with RA, including apoptosis, cell proliferation, cell migration, inflammatory response, immune response, angiogenesis, extracellular matrix degradation and bone resorption. In particular, we observed that the bone resorption BP included upregulated genes that negatively regulate osteoblast functions, such as HES1, and genes that are involved in osteoclast development, like ATP6AP1, TCTA and SBNO2 [14][15][16]. In this regard, we also have to mention the upregulation of CSF1/MG-CSF that is one of the most important soluble factors responsible for osteoclast maturation and survival [17]. Interestingly, PLCB1, a positive regulator of osteoblast differentiation [18], was downregulated in RA samples.
Several DEGs were involved in well-known pathways including Wnt, TNF, type I interferon, p38 MAP kinase, NF-kB, Toll-like receptors, Jak-Stat, PI3K and mTOR signalling that have already been associated with RA pathogenesis. A selection of genes involved in the abovementioned functional classes is given in Tables 1 and 2.   All the differentially expressed transcripts were submitted to a pathway enrichment analysis that highlighted other meaningful signalling networks in which modulated genes were involved. These pathways included, for example, signalling that operates in vascular biology (i.e., PAR, uPA/uPAR, PDGFR, endothelins and VEGF signalling), interferon-gamma, EGF-receptor, Arf6, IL-5, IL-3 and S1P1 signalling pathways. All the enriched pathways (Bonferroni corrected p-value ≤ 0.05) are listed in Table 3.

Selected Long Non-Coding RNAs Modulated in RA Patients Have the Potential to Regulate Genes Differentially Expressed in the Disease
To strengthen the significance of our analysis, we interrogated the StarBase database to select only modulated lncRNAs for which experimentally validated microRNA (miRNA) targets had already been annotated and, by this criterion, six out of 97 lncRNAs were filtered: FTX, HNRNPU-AS1, MIATNB, RP11-498C9.15, RP4-714D9.5 and RP11-73E17.2. All these lncRNAs were downregulated except RP11-498C9.15 and RP11-73E17.2, which were overexpressed in RA samples (Table 4).
To find all the possible interactions among modulated genes and the selected lncRNAs, we verified if they could target miRNAs able to regulate RA-DEGs. We therefore analysed the complete list of genes regulated by the miRNA targets of the six lncRNA that were validated by high-throughput technologies, and selected only those microRNAs that could modulate differentially expressed genes in RA patients. We thus observed that all the selected lncRNAs, via their miRNA targets, could control genes differentially expressed in RA patients, but only FTX, HNRNPU-AS1 and RP11-498C9. 15 targeted quite a large number of DEGs (Tables 4 and S3). To gain insight on the potential role played by the selected lncRNAs in regulating gene clusters that were most probably associated with the disease pathogenesis, we performed a pathways enrichment analysis of all their targeted DEGs. This approach led us to observe that genes targeted by RP11-498C9.15 significantly enriched (Bonferroni p-value < 0.05) the largest number of signalling pathways and, interestingly, these pathways were almost the same as those that were globally enriched by the 942 modulated genes (Table 5). On the contrary, FTX and HNRNPU-AS1 targeted a small number of enriched pathways that were also targeted by RP11-498C9.15, whereas MIATNB RP4-714D9.5 and RP11-73E17.2 did not target any enriched pathways (Table S4). These results led us to suppose that RP11-498C9.15 may exert a major impact on gene modulation associated with RA so, to examine its possible involvement in the disease pathogenesis, we focused our analysis on genes targeted by this lncRNA.  We observed that several DEGs that were targeted by RP11-498C9.15 were involved in meaningful biological settings, such as the immune and inflammatory response, bone metabolism, apoptosis regulation, etc. (Figure 1). Moreover, the modulation of several of them has already been associated with RA. Upregulated targeted genes were, for example, implicated in T cell survival (DUSP5) [19] or activation (i.e., MAL, LAT and CD81) [20][21][22], whereas others were involved in B cell development (i.e., PRDM1/BLIMP-1, MEF2D, LRRC8A and ZBTB7A) [23][24][25]. Notably, CD81 has been found to be upregulated in RA synoviocytes [26], while single-nucleotide polymorphisms of the gene BLIMP-1 have already been described in RA [23]. A fair number of targeted genes played a role in the inflammatory response, like DDIT4/REDD-1, COTL1/CLP, SPATA2, PTGER4, PTGES2, NR4A2, LTB and KLF2. Among these, the pro-inflammatory gene DDIT4 has a role in NF-kB activation [27], which has been shown to modulate the production of inflammatory cytokines implicated in RA joint pathology [28]. SPATA2 is a component of the TNF-alpha signalling [29], while COTL1 is able to regulate the production of leukotriene A4 [30] and NR4A2 and LTB are highly expressed in inflamed RA synovial tissues [31,32]. The gene product of PTGES2 converts prostaglandin H2 to prostaglandin E2, whereas PTGER4 encode for the prostaglandin E2 receptor and, interestingly, this molecule is involved in the differentiation and expansion of T helper lymphocytes, a process that is involved in RA onset [33]. In addition, this molecule exerts an inhibitory action on human bone marrow stromal cells-mediated bone matrix mineralization [34].
Notably, other upregulated targeted genes included transcripts that are involved in bone erosion, like CRTC2, a negative regulator of BMP2-induced osteogenic cell differentiation [35], and JUNB, involved in osteoclast development [36] and strongly expressed in RA fibroblast-like cells [37]. In addition, it has been demonstrated that JUNB promotes Th17 cells' development, but inhibits T regulatory cells' fate during chronic autoimmunity [38].
We also observed that RP11-498C9.15 targeted upregulated genes involved in the negative control of the apoptotic process (i.e., SH3BGRL3 and WEE1) and in autophagy (i.e., ATG4D), a mechanism that seems to be implicated in apoptosis resistance in RA as well as in the development of several autoimmune diseases including RA [39]. Interestingly, RP11-498C9.15 targeted FOXJ3 and gene polymorphism of this transcript have been associated with RA [40].

LncRNA RP11-498C9.15 Targets RA-Associated Meaningful Signalling Pathways
Since, as mentioned above, DEGs targeted by RP11-498C9.15 enriched signalling pathways that were also over-represented in the whole RA dataset, we analysed the entire list of these molecular networks ( Table 5) and found that almost all are involved in pathogenetic mechanisms that may play a role in the development of RA and, interestingly, that the activation of several of them has already been associated with RA pathogenesis. Indeed, the involvement of phosphatidylinositide-3-kinase (PI-3K), AKT, mTOR and sphingosine-1 phosphate (S1P1) signalling in RA pathogenesis has been extensively documented [41,42].
The most enriched pathway was the beta-1 integrin signalling, a cluster of molecules that are well represented at the synovial lining layer, where synovial cells adhere to cartilage. It has been observed that the RA-associated pro-inflammatory milieu promotes the synthesis of beta-1 integrins [43] and indeed, fibroblast, macrophages and endothelial cells of RA synovial tissue express high levels of these molecules. Notably, the beta-1 integrin binding to laminin at the synovial lining, strongly increases the expression of metalloproteases [43] and signalling by beta-1 integrins leads to immune cell activation, stimulates cell migration, cytokine production and new vessels formation. Besides the integrins pathway, interferon-gamma was the second most enriched signalling-not surprising since this cytokine is rapidly accumulated in RA synovial tissue by activated T cells and its level of expression correlates with the RA radiographic severity [44].
A high level of thrombin activity has been found in RA patients, and it has been suggested that this molecule may have a strong mitogenic effect on synovial fibroblast-like cells, thus possibly playing a significant role in RA pathogenesis. Moreover, it has been observed that the high levels of thrombin that have been detected in RA synovium are associated with increased expression of platelet derived growth factor beta (PDGF-b) [45] and, interestingly, its pathway was also enriched. Remarkably, PDGF receptor activation promotes the RA synoviocytes' pro-destructive behaviour [46].
RA synovial fibroblasts were reported to express high levels of plasminogen activator and urokinase (uPA), and a pro-inflammatory role has been hypothesized for this molecule, since anti-inflammatory glucocorticoids strongly suppress uPA gene expression [47]. Finally, endogenous endothelins are involved in articular inflammation by modulating inflammatory pain, oedema formation and leukocyte migration. In addition, an elevated plasma level of endothelin-1 has been observed in RA, which may be associated with the symptoms of vascular dysregulation frequently observed in RA patients [48]. RP11-498C9.15 also targeted members of the proteoglycan glypican and syndecan pathways, which are structural molecules thought to govern cell migration, tissue invasion and angiogenesis [49], thus possibly playing a crucial role in fibroblast-like synoviocytes' behaviour [50]. Other enriched pathways involved in angiogenesis were the signalling of focal adhesion kinases (FAKs), epidermal growth factor (EGF) receptor and VEGF/VEGF-receptor. The former has been implicated in both RA inflammatory angiogenesis and in RA synovial fibroblasts' pro-invasive activity [51,52], whereas the second is particularly involved in synovial fibroblast proliferation. Indeed, EGF and EGFR serum concentrations are increased in RA patients compared to healthy subjects [53]. VEGF also has pro-inflammatory and bone-destructive skills and its level correlates with the RA disease activity [54].
It has been described that RA osteoclasts exhibit increased activity and, interestingly, pathways involved in bone erosion were also significantly enriched, including Arf6, as well as hepatocyte growth factor (HGF) signalling. The Arf6 pathway plays an important role in osteoclast maturation [55], while HGF promotes osteoclast activation and inhibits osteoblast differentiation, thus favouring bone-erosive processes. Interestingly, plasma levels of this cytokine can predict joint damage in RA [56].
RP11-498C9.15 also targeted enriched pathways involved in glucose metabolism like "insulin" and "insulin growth factor (IGF1)" signalling. The activation of these pathways may reflect the abnormal glucose metabolism that characterizes a good percentage of RA patients [57] and has been correlated with the degree of systemic inflammation [58]. In addition, it has been observed that the pro-inflammatory cytokine TNF may promote insulin resistance by phosphorylation of the insulin receptor [59].
Enriched targeted pathways also included the granulocyte-macrophage colony-stimulating factor (GM-CSF) signalling, a cytokine highly expressed in both RA synovial fluid and tissue and on circulating mononuclear cells from RA patients. Given its prominent role in macrophage differentiation and activation, it has been suggested that GM-CSF inhibition may represent a favourable therapeutic approach to treat RA. Indeed, early phases of clinical trials evaluating anti-GM-CSF therapy demonstrated its potential clinical benefit in RA patients [60].
Finally, other enriched targeted pathways included interleukin-5 (IL-5) and interleukin-3 (IL-3) signalling. IL-5 is known to stimulate B cells' growth, increasing immunoglobulin secretion, whereas IL-3 is a potent inducer of RANKL expression in human basophils and an important role for this molecule in the early phase of collagen-induced arthritis has been demonstrated [61].

LncRNA RP11-498C9.15 Targets Highly Connected Genes in the RA Transcriptome
Since it is well known that the targeting of highly connected genes can have a more pronounced impact on the development of a disease than the modulation of transcripts that show no functional interactions, we wanted to verify that RP11-498C9.15 could control highly interacting genes in the RA transcriptome.
With this purpose in mind, we first built a protein-protein interaction (PPI) network that included all the experimentally validated functional interactions among the protein products of the 942 modulated genes in RA; thereafter, we performed a modular analysis to find the areas of the network in which the most highly connected genes were clustered. The obtained network included 755 nodes (genes) and 2496 edges (pairs of interactions) and exhibited a good enrichment p-value (p < 1.0e −16 ) (Figure 2). Interestingly, the topological analysis of the PPI network revealed that, via their miRNAs targets, RP11-498C9.15 was connected to genes with a high degree of connectivity (Figure 3), whereas the modular analysis highlighted the presence of six modules (Table S5) that included genes regulated by RP11-498C9.15 (Table 6).  In module M1, CUL3, HNRNPA0 and SOCS3 were targeted and, notably, the latter has been found to be upregulated in RA fibroblast-like synoviocytes and RA PBMCs [62]. In module M2, RP11-498C9. 15 targeted ANO6, which is involved in bone mineralization, whereas in module M3, the lncRNA targeted ATP5G2, TERF2 and the abovementioned PTGES2, which, as previously mentioned, operates the conversion of prostaglandin H2 to prostaglandin E2. Module M4 included seven targeted genes, ATM, NACA, PCGF5, SEL1L, SURF4, THBS1 and VPS45, while in module M5, AKT2, JUND, RELA and SRC were targeted. In this regard, we have to mention that AKT2 induces the proliferation and migration of RA fibroblast-like synoviocytes [63], JUND can contribute to bone erosion [64], RELA is a key component of the NF-kB transcription factor crucially involved in RA [65], and SRC has bone resorptive properties.
In module M6, STX12, which is involved in autophagy, was targeted.
The level of expression of RP11-498C9.15, selected miRNA and gene targets were validated by RT-PCR ( Figure S1). Statistically significant differences between patients and healthy subjects were found in the expression levels of all the tested transcripts.

Discussion
RA is an inflammatory chronic autoimmune disease and its pathogenesis is influenced by genetic, environmental and epigenetic factors [1]. LncRNAs are key components of the epigenetic machinery that can regulate chromatin remodelling and gene expression by interacting with other epigenetic factors and with genes. Interestingly, lncRNAs seem to be involved in the development of several autoimmune diseases [3].
In the present work, we simultaneously profiled a large number of coding and non-coding transcripts in the same cohort of early-phase-RA patients in order to identify the lncRNAs that most probably shape the outcomes of crucial biological processes strictly associated with RA pathogenesis.
The criteria adopted to select deregulated lncRNAs in RA rely on the analysis of experimentally validated functional interaction among coding and non-coding partners of the RA transcriptome. Using this approach, six lncRNAs were filtered; such lncRNAs, via their miRNA targets, can modulate differentially expressed genes in RA patients.
Since it is now a common notion that disease can be explained in terms of molecular pathways' perturbation, we performed a pathway enrichment analysis of all modulated genes that were targeted by the six selected lncRNAs. Through this analysis we observed that transcripts targeted by RP11-498C9.15, unlike transcripts targeted by the other five lncRNAs, showed a good coverage of pathways significantly modulated in the disease. Indeed, RP11-498C9.15 was able to target almost all the signalling pathways in which genes modulated in RA samples are involved, thus showing the best correlation to the RA transcriptome.
Notably, RP11-498C9.15 targeted meaningful signalling pathways that have been already associated with RA and RA-fibroblast-like synoviocytes' behaviours, while the other selected lncRNAs targeted only a few pathways that were also modulated by RP11-498C9. 15.
By analysing the whole pattern of functional interactions among the modulated genes, we could observe that RP11-498C9.15 targeted modules of the most highly connected genes in the RA interactome, which are believed to be principally involved in the disease onset.
Also, for this reason, although we cannot exclude the involvement of the other selected lncRNAs, we believe that RP11-498C9.15 may play a crucial role in the pathogenesis of RA.
We are aware that a limitation of this work is the small number of samples analysed, but this is mainly due to the difficulty of recruiting RA patients in the early phase of the disease and in the absence of any treatment.
We believe that, especially if our results are confirmed on a larger cohort of patients, the lncRNA RP11-498C9.15 deserves to be identified as a candidate in the design of novel therapeutic strategies in RA.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4409/8/8/816/s1, Table S1: Genes modulated in RA patients versus healthy subjects, Table S2: LncRNAs modulated in RA patients versus healthy subjects, Table S3: Modulated genes in RA patients that are targeted by miRNA targets of the selected lncRNAs, Table S4: Enriched signalling pathways involving modulated genes that are targeted by the selected lncRNAs, Table S5: Modulated genes in RA patients that are included in the six modules, Figure S1: Expression of selected genes, RP11-498C9.15 and miR-520e in RA patients versus healthy subjects. Bars indicate SD. The histograms represent the mean of the results obtained in 20 healthy donors and in 20 RA patients. The Mann-Whitney test was used to compare the two groups' means.

Conflicts of Interest:
The authors declare no conflict of interest.