In Systemic Sclerosis, a Unique Long Non Coding RNA Regulates Genes and Pathways Involved in the Three Main Features of the Disease (Vasculopathy, Fibrosis and Autoimmunity) and in Carcinogenesis

Systemic sclerosis (SSc) is an autoimmune disease characterized by three main features: vasculopathy, immune system dysregulation and fibrosis. Long non-coding RNAs (lncRNAs) may play a role in the pathogenesis of autoimmune diseases and a comprehensive analysis of lncRNAs expression in SSc is still lacking. We profiled 542,500 transcripts in peripheral blood mononuclear cells (PBMCs) from 20 SSc patients and 20 healthy donors using Clariom D arrays, confirming the results by Reverse Transcription Polymerase-chain reaction (RT-PCR). A total of 837 coding-genes were modulated in SSc patients, whereas only one lncRNA, heterogeneous nuclear ribonucleoprotein U processed transcript (ncRNA00201), was significantly downregulated. This transcript regulates tumor proliferation and its gene target hnRNPC (Heterogeneous nuclear ribonucleoproteins C) encodes for a SSc-associated auto-antigen. NcRNA00201 targeted micro RNAs (miRNAs) regulating the most highly connected genes in the Protein-Protein interaction (PPI) network of the SSc transcriptome. A total of 26 of these miRNAs targeted genes involved in pathways connected to the three main features of SSc and to cancer development including Epidermal growth factor (EGF) receptor, ErbB1 downstream, Sphingosine 1 phosphate receptor 1 (S1P1), Activin receptor-like kinase 1 (ALK1), Endothelins, Ras homolog family member A (RhoA), Class I Phosphoinositide 3-kinase (PI3K), mammalian target of rapamycin (mTOR), p38 mitogen-activated protein kinase (MAPK), Ras-related C3 botulinum toxin substrate 1 (RAC1), Transforming growth factor (TGF)-beta receptor, Myeloid differentiation primary response 88 (MyD88) and Toll-like receptors (TLRs) pathways. In SSc, the identification of a unique deregulated lncRNA that regulates genes involved in the three main features of the disease and in tumor-associated pathways, provides insight in disease pathogenesis and opens avenues for the design of novel therapeutic strategies.


Introduction
Systemic Sclerosis (SSc) is a rare autoimmune disease characterized by endothelial disfunction, dysregulation of the immune system and increased extracellular matrix deposition leading to extended involvement and fibrosis of skin and internal organs and resulting in a remarkable heterogeneity in clinical features and in disease course [1].
Several factors have been shown to contribute to the onset of the disease, such as genetic susceptibility, environmental factors including viral infections [2,3] and epigenetic mechanisms, such as microRNAs (miRNAs) [4,5].
We have already reported the gene expression profiling of peripheral blood mononuclear cells (PBMCs) obtained from patients affected by SSc [6] showing the presence of cancer-related gene signatures.
Indeed, an increased frequency of different types of cancer, including breast, lung, and hematologic malignancies in SSc has been already described [7], and has been associated with the presence of particular autoantibodies [8]. In addition, SSc may also be a paraneoplastic disease [9,10] suggesting a strong link with cancer.
In this study we focused our attention on the epigenetic mechanisms that may be involved in SSc pathogenesis by analyzing the expression profiles of long non-coding RNAs (lncRNAs) in SSc patients.
For the first time, we propose an integrated analysis of lncRNAs, miRNAs and gene expression profiles in SSc patients leading to the identification of lncRNAs modulated in the disease. We identified a lncRNA involved both in pathogenetically relevant molecular pathways of the disease and in tumor-associated pathways, strongly associated to the disease.

Patients
Twenty patients affected by SSc, attending the Unit of Autoimmune Diseases at the University Hospital of Verona, Northern Italy, and 20 sex and age matched healthy controls were enrolled. Both patients and controls were subjects of Caucasian origin from Northern Italy. All patients were diagnosed accordingly to the American College of Rheumatology (ACR)/European League Against Rheumatism (Eular) classification criteria for SSc [11]. The distinction between limited (lSSc) and diffuse cutaneous SSc (dSSc) was performed following the criteria proposed by Le Roy et al. [12]. Ten patients were affected by lSSc and 10 by dSSc. The clinical and demographic features of the patients enrolled in the study are reported in Table 1.
Samples obtained from 10 patients with lSSc, 10 patients with dSSc and 20 healthy controls were used for the gene and lncRNAs expression analysis. Blood samples were collected from patients with active disease (European Scleroderma Study Group activity index > 3) [13] within 24 months after diagnosis and in the absence of immunosuppressive therapies. Samples were processed immediately after the blood draw.
Twenty more patients and 20 healthy controls have been analyzed as validation groups (Supplementary Table S1).
A written informed consent was obtained by all the participants to the study and the study protocol was approved by the Ethical Committee of the Azienda Ospedaliera Universitaria Integrata di Verona. All the investigations have been performed according to the principles contained in the Helsinki declaration. Table 1. Demographic and clinical and features of patients and healthy subjects enrolled in the study.

Microarray Analysis
Blood samples collection was performed using BD (Becton Dickinson)Vacutainer K2EDTA tubes and 21-gauge needles.
PBMCs were isolated by Ficoll-HyPaque (Pharmacia Biotech, Quebec, Canada) gradient centrifugation. The PBMCs distribution was similar between patients and controls. Total RNA extraction from PBMCs (10 7 cells) was achieved with miRNeasy mini kit (Qiagen GmbH, Hilden, Germany). cRNA preparation, samples hybridization and scanning were carried out 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, Waltham, MA, USA) gene chip. Signal intensities were analyzed with the Transcriptome Analysis Console (TAC) 4.0 software (Applied Biosystem, Foster City, CA USA by Thermo Fisher Scientific).
The Human Clariom D arrays interrogates more than 540,000 human transcripts starting from as little as 100 pg of total RNA. Signals intensity was background-adjusted, normalized, and log-transformed using the Signal Space Transformation (SST)-Robust Multi-Array Average algorithm (RMA).
Differentially expressed genes that showed an expression level at least 1.5 fold different in the test sample versus control sample at a significant level (p ≤ 0.01) were chosen for final consideration. p-values were calculated applying the unpaired t-test and Bonferroni multiple testing correction.
The list of Gene targets of microRNAs (miRNAs) that are targeted by ncRNA00201 were obtained from the FunRich database (www.funrich.org/) [15].

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/) [16].
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).

Real Time PCR of ncRNA00201
A total of 500 ng of total RNA was treated with 1 unit of DNase I Amplification Grade (Invitrogen; Carlsbad, CA, USA.). 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. The primers used for the detection were: GCA GGA GAA TCG CTT GAA C (forward) and ACC CTA CCA TCC AAC TTC AC (reverse).

Real Time PCR of Genes Modulated in SSc Patients
First-strand cDNA was generated using the SuperScript III First-Strand Synthesis System for Reverse transcription Polymerase-chain Reaction (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 one minute. Thermocycling and signal detection were performed with 7500 Sequence Detector (Applied Biosystems). Signals were detected according to the manufacturer's instructions. This technique allows the identification of the cycling point where 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 Software and Digital Store (SDS) analysis software (version 2.4.0). For each amount of RNA tested triplicate Ct values were averaged. Because Ct values vary linearly with the logarithm of the amount of RNA, this average represents a geometric mean.

Real Time PCR of microRNAs Targeted by ncRNA00201
miRNAs expression was assayed by TaqMan ® Advanced miRNA assays chemistry (Applied Biosystems, Foster City, CA, USA). Briefly, 10 ng of total RNA was reverse transcribed and pre-amplified with TaqMan ® Advanced miRNA cDNA synthesis kit following manufacturer's instructions (Applied Biosystems, Foster City, CA, USA). 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 2X Fast Advanced Master Mix and TaqMan ® Advanced miRNA assays for miR-30b-5p, miR-30e-5p, miR-302a-3p, miR-520d-3p, miR-613 and miR-9-5p. The mean of Ct for hsa-miR-16-5p and hsa-miR-26a-5p expression was used to normalize miRNA expression. Real time PCR were carried out in triplicate on a QuantStudio 6 Flex instrument (Applied Biosystems, Foster City, CA, USA). Expression values were reported as fold change with respect to healthy controls by ∆∆Ct method using QuantStudio Real-Time PCR system software v. 1.3.

Results
With the aim of identifying lncRNAs potentially involved in SSc pathogenesis, the expression of more than 540,000 human transcripts, including those ascribed to more than 50,000 lncRNAs was profiled at the same time, in a cohort of 40 PBMCs samples (20 SSc and 20 healthy subjects). Transcriptional profiles of SSc patients and healthy subjects were compared and, after a robust filtering procedure (Bonferroni-corrected p-value ≤ 0.01 and fold change ≥ |1.5|), 837 coding-genes were significantly modulated (Supplementary Table S2). More than one hundred lncRNAs were detectable, but only one, long non-coding heterogeneous nuclear ribonucleoprotein U antisense RNA 1 (HNRNPU-AS1, also known as ncRNA00201) satisfied the above mentioned criteria (statistical and fold change) showing a significant modulation ( Table 2). ncRNA00201 has been showed to be overexpressed in pancreatic ductal adenocarcinoma (PDAC) tissues and cell lines, and plays a role in the regulation of cell proliferation, invasion, and migration. Indeed, the knockdown of this lncRNA reduces proliferation and in vitro migration and invasion of PDAC cell lines [19].
The functional classification by Gene Ontology (http://www.geneontology.org/) of the 837 differentially expressed genes (DEGs) revealed the presence of a large number of differentially expressed coding genes involved in biological processes (BPs) and signaling pathways that are strictly connected to SSc pathogenesis, including, immune response BP, interferon alpha/beta signaling, IL-6 mediated signaling events BP, inflammatory response BP, apoptosis BP, angiogenesis BP, Vascular Endothelial Growth Factor (VEGF) and VEGFR signaling, cell adhesion BP, blood coagulation BP, endothelin signaling, PDGFR-beta signaling, EGF receptor (ErbB1) signaling, TGF-beta receptor signaling, extracellular matrix (ECM) component and organization, proteoglycan syndecan mediated signaling events, and positive regulation of fibroblast proliferation BP.    Table 2 shows a selection of DEGs involved in these meaningful biological processes and pathways.
The modulation of genes involved in immune response highlighted the immune system alterations that characterize SSc and other autoimmune diseases. Not surprising, in SSc samples we found upregulation of genes belonging to the type alpha and beta interferon signaling pathways, that is a common trait of autoimmune diseases and, this finding has also been previously overserved in SSc PBMCs [20]. In addition, we found a strong upregulation of interleukin-6 (IL-6) (F.C. 12.69) a cytokine involved in inflammatory and autoimmune diseases including SSc and, besides that, we found overexpression of several members of its signaling pathway (see Table 2).
The presence of modulated genes involved in the inflammatory process reflected the chronic inflammation that accompanies SSc. Moreover the upregulation of genes involved in apoptosis is not surprising, since apoptosis is the first event of vascular damage in early SSc [2]. The modulation of genes involved in angiogenesis gives evidence of the angiogenic response to tissue ischemia and vascular damage that occurs in the initial stages of SSc. We also have to mention that, members of the proangiogenetic VEGF and VEGFR signaling pathway were modulated in SSc (see Table 2).
The vascular response to endothelial damage consists of endothelial activation with increased expression of adhesion molecules and, indeed, several adhesion molecules were overexpressed in SSc samples.
As it is well known, the spectrum of vascular abnormalities that characterizes SSc also involves disorders of the blood coagulation-fibrinolysis system [21] and, noteworthy, we found upregulation of several genes involved in blood coagulation.
The earliest vascular changes in SSc also include a dysregulated control of vascular tone and indeed we have to mention that endothelins are among the strongest vasoconstrictors known, and are involved in vascular diseases of several organs. Interestingly in SSc samples several components of the endothelins pathway were overexpressed. Moreover, we found modulation of members of the platelet derived growth factor receptor-beta (PDGFR-B). PDGFR-B has been implicated in SSc [22] and it plays a role in activation of the epidermal growth factor receptor (EGFR) that is also involved in the disease. In particular, in scleroderma fibroblasts, aberrant activation of EGF-mediated signaling pathways, can upregulate the receptor II of transforming growth factor beta (TGFBRII) [23], the growth factor primary playing a role in the development of connective tissue fibrosis that is typically involved in SSc pathogenesis and, noteworthy, members of the EGFR pathway were modulated in SSc. Moreover the deregulation of both PDGF and EGF pathways has been already described in a previous study on gene expression profiles in SSc peripheral blood [24].
Consistently with these findings, other modulated genes were involved in the signaling of the profibrotic cytokine TGF-beta.
Fibrosis is supported by extracellular matrix (ECM) remodeling and, consistently, we found overexpression of ECM components and of molecules involved in ECM organization. Moreover, we found an increased expression of molecules involved in proteoglycan syndecan-mediated signaling events (see Table 2) and of genes involved in the positive regulation of fibroblasts proliferation.
The expression levels of ncRNA00201, selected miRNAs and coding genes were validated by RT-PCR in the entire series of the analyzed patients and in an expanded panel of SSc patients (20) and healthy controls (20 subjects) (Supplementary Table S1). Significantly different expression levels were found for all tested transcripts as compared to healthy controls ( Supplementary Figures S1 and S2).
Gene expression results were also confirmed by the detection of several soluble mediators (including interleukin-6, IL-6; macrophage migration inhibitory factor, MIF; selectin P, SELP; chemokine C-X-C motif ligand 10, CXCL10 and basigin, BSG) in the sera of SSc patients. As shown in Supplementary Figure S3 the serum levels of all the tested molecules were significantly different in SSc patients when compared to healthy donors.
Noteworthy, the pathways enrichment analysis confirmed that the important signaling networks described in Table 2 were all included in the list of pathways significantly (p < 0.05) enriched in the 837 genes, along with others that were likewise connected to autoimmune and inflammatory response, vascular damage (i.e., apoptosis, platelet activation and aggregation, Urokinase-type plasminogen activator (uPA) and uPAR-mediated, Thrombin/protease-activated receptor (PAR) signaling etc.) and fibrosis (i.e., glypican, syndecan-3 mediated, regulation of nuclear Small Mother Against Decapentaplegic (SMAD) 2/3 signaling etc.) (Supplementary Table S3).
Since modulated genes were well representative of the main features of the disease, we decided to verify if the only modulated lncRNA could be functionally connected to SSc transcriptome, thus playing a role in SSc pathogenesis. To this purpose we extracted the complete list of experimentally validated genes and microRNAs (miRNAs) targets of ncRNA00201, and we found that 56 miRNAs and 31 genes were annotated as target of ncRNA00201 (Supplementary Table S4). One of the gene targets, named hnRNPC, belonged to a subfamily of ubiquitously expressed heterogeneous nuclear ribonucleoproteins (hnRNPs) and interestingly encoded for a known antigen identified in SSc [25].
To dissect all the possible connections between ncRNA00201 and the SSc transcriptome we analyzed the lists of genes targeted by each of the 56 miRNAs (3759 genes) selecting only transcripts that also resulted significantly modulated on the array. Among the 3759 gene targets, 138 were modulated in SSc patients and were further analyzed along with their targeting miRNAs (47/56) (Supplementary Table S5). We therefore selected only those miRNAs that targeted genes with evidence of modulation in SSc patients to bona fide outline authentic interactions that are most probably established in SSc.
The analysis of signaling pathways, in which the 138 DEGs may be involved, showed that these genes were present in the 83% (147/176) of pathways significantly enriched in the SSc transcriptome (Supplementary Table S6).
A protein-protein interaction (PPI) network including all the protein products of the 837 modulated genes that showed experimentally validated interactions was constructed and, the obtained network included 693 nodes (interacting partners) and 2226 edges (interactions), showing a good PPI enrichment p-value (<10 −16 ) indicating that the network contained significantly more interactions than expected. In other words, the connected proteins had more interactions among themselves than what would be expected for a random set of proteins of similar size, drawn from the genome. Such an enrichment indicates that these proteins are biologically connected, as a group.
The topological analysis of the PPI network showed the presence of a large number of targeted genes in areas with high density of connections (Figure 1). When a modular analysis was performed, we could identify six areas of the PPI network in which the most highly connected genes were present (modules).
Since the modulation of highly connected genes is expected to have a wider impact on the transcriptome than the targeting of several isolated genes, we checked if ncRNA00201 could target genes included in the modules and, noteworthy, we observed that genes modulated by miRNA targets of ncRNA00201 were present in each of the 6 modules (Supplementary Table S7).
To dissect the most relevant interactions through which ncRNA00201 could influence the pathogenesis of the disease, the miRNAs targeted by ncRNA00201 were filtered, prioritizing those that targeted at least 7 genes in the whole SSc transcriptome (i.e., the 837 modulated genes) and at least two modules. By these criteria, 26 miRNAs target of ncRNA00201 were selected. Table 3 shows the 26 selected miRNAs and summarizes their targeted genes in the SSc transcriptome and in the modules. When a modular analysis was performed, we could identify six areas of the PPI network in which the most highly connected genes were present (modules).
Since the modulation of highly connected genes is expected to have a wider impact on the transcriptome than the targeting of several isolated genes, we checked if ncRNA00201 could target genes included in the modules and, noteworthy, we observed that genes modulated by miRNA targets of ncRNA00201 were present in each of the 6 modules (Supplementary Table S7).
To dissect the most relevant interactions through which ncRNA00201 could influence the pathogenesis of the disease, the miRNAs targeted by ncRNA00201 were filtered, prioritizing those that targeted at least 7 genes in the whole SSc transcriptome (i.e., the 837 modulated genes) and at least two modules. By these criteria, 26 miRNAs target of ncRNA00201 were selected. Table 3 shows the 26 selected miRNAs and summarizes their targeted genes in the SSc transcriptome and in the modules. Interestingly, among the selected miRNAs, miR-30b-5p and miR-31-5p were previously described as deregulated in SSc whereas, all the 26 miRNAs are modulated in several types of human cancer ( Table 3).
The selected miRNAs targeted module-associated-genes with important biological significance, indeed these genes were involved in inflammation, in the modulation of immune response regulating signaling, in maintaining vascular homeostasis, in tissue remodeling, and also in metastasis development in many cancers.
Indeed, Module 1 (M1) was targeted by 20 miRNAs (see Table 3) that globally modulated SOCS3, UBE2F, UBE2J1, IRF9, and CBLB. SOCS3, a major regulator of inflammation, is rapidly induced by JAK-STAT (Janus kinase-signal transducer and activator of transcription) signaling during the inflammatory response. UBE2F and UBE2J1 are ubiquitin conjugating enzymes that are involved in inflammation and in metastasis development in many cancers [54]. IRF9 is strongly induced by type I interferons and CBLB is involved in the regulation of immune response modulating signaling by both T-cell and B-cell receptors [55].
M3 was targeted by miR-31-5p that modulated KDELR2. This molecule is involved in extracellular matrix degradation and in the signaling of the Src family proteins that is activated during cell invasion leading to metastasis development [58].
In M4 the gene APIS3 was targeted by members of the miR-181s family (including miR-181-a/b/c/d-5p). Lacking of AP1S3 gene has been associated to the development of skin inflammation [59] and interestingly, AP1S3 was down-modulated in SSc samples (see Supplementary  Table S1).
M5 was the most targeted module since it was targeted by 21 out of the 26 selected miRNAs (see Table 3). In this module 5 genes were targeted including GCH1, GUCY1A3, PAX5, SPARC and BCL2. GCH1, also named GTPCH1 is implicated in the synthesis of the three nitric oxide synthases (NOS1-3) [60] and endothelial cell GCH1-dependent endothelial nitric oxide (eNOS) regulation plays an important role in maintaining vascular homeostasis [61]. GUCY1A3 encodes for the α1 subunit of nitric oxide (NO) receptor and is involved in the control of smooth muscle cells relaxation [62].
PAX5 is indispensable for the B-cell lineage specification and maintenance [63] and has been implicated in human B cell malignancies, including acute lymphoblastic leukemias and non-Hodgkin lymphomas [64].
SPARC (also known as osteonectin) is a glycoprotein of the extracellular matrix expressed in concomitance with tissue remodeling and wound repair. Interestingly SPARC upregulation has been described in cultured dermal fibroblasts from patients with SSc [65] and it has been observed that SPARC silencing can reduce collagens overproduction in SSc skin fibroblasts [66].
The BCL2 gene is involved in apoptosis and its co-expression with the oncogene MYC is associated with a poor prognosis in B-cell lymphoma patients [67]. Interestingly both BCL2 and MYC were upregulated in SSc patients (see Supplementary Table S1).
In module M6, PSME1 was targeted by miR-204-5p and miR-211-5p. PSME1, also named PA28, is a proteasome component that influences antigen processing. Proteasomal system plays crucial roles in important biological processes including cell differentiation, proliferation, apoptosis, transcriptional activation and angiogenesis, and is a pivotal target for treatment of several diseases including autoimmune diseases and cancer [68].
We also observed that a large proportion of the pathways enriched (p < 0.05) in modulated genes targeted by the 26 selected miRNAs were linked to the three main features of the disease (i.e., immune and inflammatory response, vasculitis and fibrosis) thus confirming their possible involvement in SSc pathogenesis. Figure 2 shows a selection of the above mentioned enriched signaling pathways in which are involved modulated genes targeted by 11 of the 26 selected miRNAs. Interestingly, several of the aforementioned pathways are also deregulated in many forms of cancer. These pathways included for example EGF receptor, ErbB1 downstream, Sphingosine-1phosphate receptor 1 (S1P1), Arf6 downstream, ALK1, Endothelins, RhoA, Class I PI3K (Phosphatidylinositol 3-kinase), mTOR, p38 MAPK, RAC1, TGF-beta receptor, MyD88 and Toll-like receptors signaling pathways (Figure 3). Interestingly, several of the aforementioned pathways are also deregulated in many forms of cancer. These pathways included for example EGF receptor, ErbB1 downstream, Sphingosine-1-phosphate receptor 1 (S1P1), Arf6 downstream, ALK1, Endothelins, RhoA, Class I PI3K (Phosphatidylinositol 3-kinase), mTOR, p38 MAPK, RAC1, TGF-beta receptor, MyD88 and Toll-like receptors signaling pathways (Figure 3). We thereafter aimed at identifying modulated genes targeted by the selected miRNAs whose deregulation could have the widest impact on the global disease-associated gene modulation. We therefore narrowed our analysis to modulated genes targeted by the 26 miRNAs that were included in the six modules and we identified five transcripts, namely IRF9, GUCY1A3, SOCS3, BCL2 and GNAI2 that were targeted by the highest number of selected miRNAs (see Table 2).
Interestingly, we observed that they were associated to a large number of pathways (136) and that 82 (60%) of such pathways are significantly enriched (p < 0.05) in the whole SSc transcriptome such as for example interferon alpha/beta, apoptosis, endothelins, PDGF receptor beta, TGF-beta receptor and regulation of SMAD (Supplementary Table S8).

Discussion
Though it has been widely recognized that lncRNAs play a pivotal role in the regulation of autoimmune diseases [69], an extensive analysis of lncRNAs in SSc is still lacking. We therefore performed a comprehensive analysis assessing the expression profiles of a very large number of lncRNAs in SSc patients, and we could find that, differently from what we had already observed in other autoimmune diseases [70], a single lncRNA, namely ncRNA00201, appeared to be significantly modulated in SSc.
Among the gene targets of ncRNA00201, there is the hnRNPC, which encodes for a known autoantigen in SSc [25].
Interestingly, ncRNA00201 has been showed to be involved in cancer proliferation [19] and this finding may reinforce the hypothesis of a link between SSc and tumor development, as we have recently suggested [6].
By a combined analysis of the expression profiles of a vast number of coding and non-coding transcripts, we have been able to observe that ncRNA00201 alone controls biological processes and pathways closely related to the three main features of SSc, including immune and inflammatory response, vasculitis and fibrosis.
We moreover identified 26 miRNAs by which ncRNA00201 could modulate the expression of several genes crucial to the disease, thus highlighting the most effective lncRNA-miRNA-gene interactions in SSc. Interestingly, all the 26 miRNAs are modulated in several types of human cancer.
Indeed, changes in EGFR expression and autoantibodies to this receptor have been observed in patients with SSc [71], nevertheless there is evidence for the involvement of EGFR signaling in the development of different carcinoma types [72].
The activation of S1P signaling pathway can induce many of the alterations observed in SSc patients [73], moreover, S1P receptors regulate cell proliferation in various tumors [74].
Arf6 is involved in angiogenesis [75] and is implicated in tumor metastasis [76]. ALK1 inhibition reduces pro-fibrotic genes expression in SSc fibroblasts [77] moreover, this gene is an emerging target We thereafter aimed at identifying modulated genes targeted by the selected miRNAs whose deregulation could have the widest impact on the global disease-associated gene modulation. We therefore narrowed our analysis to modulated genes targeted by the 26 miRNAs that were included in the six modules and we identified five transcripts, namely IRF9, GUCY1A3, SOCS3, BCL2 and GNAI2 that were targeted by the highest number of selected miRNAs (see Table 2).
Interestingly, we observed that they were associated to a large number of pathways (136) and that 82 (60%) of such pathways are significantly enriched (p < 0.05) in the whole SSc transcriptome such as for example interferon alpha/beta, apoptosis, endothelins, PDGF receptor beta, TGF-beta receptor and regulation of SMAD (Supplementary Table S8).

Discussion
Though it has been widely recognized that lncRNAs play a pivotal role in the regulation of autoimmune diseases [69], an extensive analysis of lncRNAs in SSc is still lacking. We therefore performed a comprehensive analysis assessing the expression profiles of a very large number of lncRNAs in SSc patients, and we could find that, differently from what we had already observed in other autoimmune diseases [70], a single lncRNA, namely ncRNA00201, appeared to be significantly modulated in SSc.
Among the gene targets of ncRNA00201, there is the hnRNPC, which encodes for a known autoantigen in SSc [25].
Interestingly, ncRNA00201 has been showed to be involved in cancer proliferation [19] and this finding may reinforce the hypothesis of a link between SSc and tumor development, as we have recently suggested [6].
By a combined analysis of the expression profiles of a vast number of coding and non-coding transcripts, we have been able to observe that ncRNA00201 alone controls biological processes and pathways closely related to the three main features of SSc, including immune and inflammatory response, vasculitis and fibrosis.
We moreover identified 26 miRNAs by which ncRNA00201 could modulate the expression of several genes crucial to the disease, thus highlighting the most effective lncRNA-miRNA-gene interactions in SSc. Interestingly, all the 26 miRNAs are modulated in several types of human cancer.
Indeed, changes in EGFR expression and autoantibodies to this receptor have been observed in patients with SSc [71], nevertheless there is evidence for the involvement of EGFR signaling in the development of different carcinoma types [72].
The activation of S1P signaling pathway can induce many of the alterations observed in SSc patients [73], moreover, S1P receptors regulate cell proliferation in various tumors [74].
Arf6 is involved in angiogenesis [75] and is implicated in tumor metastasis [76]. ALK1 inhibition reduces pro-fibrotic genes expression in SSc fibroblasts [77] moreover, this gene is an emerging target for antiangiogenic therapy in cancer [78]. Endothelins (ET) play important roles in the in the fibrotic process associated to SSc and have a documented role in the development of cancer [79]. RhoA signaling, is frequently altered in many types of human tumors [80] and the expression of RhoA is elevated in SSc endothelial cells [81].
PI3K signaling is involved in collagen gene expression and in the antiapoptotic phenotype that characterizes SSc fibroblasts [82,83]. This pathway acts in synergy with the mTOR signaling in the development of many cancers [84]. Interestingly, mTOR blockade can reduce the pro-fibrotic action of TGF-beta and its therapeutic implications for SSc have been suggested [85].
The p38 MAPK and the Rac1 signaling were also associated to SSc fibrosis [86,87] and to cancer progression [88,89].
In conclusion, this study is the first report of a unique lncRNA that controls gene networks that may contribute to SSc pathogenesis and that may play a predisposing role to the well documented development of malignancies in SSc.
Finally, we may propose that ncRNA00201 may be a potential target in the development of novel therapeutic strategies.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-0383/8/3/320/s1, Table S1. Demographic and Clinical features of SSc patients and healthy controls that were enrolled as validation cohort. Table S2. Genes modulated in SSc patients versus healthy subjects. Table S3. Signaling pathways significantly enriched (p < 0.05) in genes modulated in SSc patients. Table S4. MicroRNAs and genes that are annotated as targets of ncRNA00201 in StarBase. Table S5. Modulated genes in SSc that are regulated by miRNA targets of ncRNA00201. Table S6. Pathways in which are involved the 138 modulated genes that are targeted by miRNAs modulated by ncRNA00201. Pathways significantly enriched in the 837 genes that were modulated in SSc patients are highlighted in orange. Table S7. Modulated genes in SSc that are included in modules. Genes targeted in the six modules are written in bold characters. Table S8. Signaling pathways in which are involved the five selected genes modulated in SSc. Figure S1. Expression of ncRNA00201 and of selected miRNAs targeted by ncRNA00201 in the 20 SSc patients and healthy controls analyzed in the study and in an expanded panel of SSc patients (20 patients) and healthy controls (20 subjects). Bars indicate SD. Figure S2. Expression of selected genes modulated in SSc patients and regulated by miRNAs targeted by ncRNA00201. Bars indicate SD. Figure S3. Serum levels of selected molecules in SSc patients and in normal subjects. The histograms represent the mean of the results obtained in 20 healthy donors and in 20 SSc patients. Mann-Whitney Test was used to compare the two groups means. Writing-review and editing, A.P. and C.L.

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