Characterization of Long Non-Coding RNAs in Systemic Sclerosis Monocytes: A Potential Role for PSMB8-AS1 in Altered Cytokine Secretion

Systemic sclerosis (SSc) is a chronic autoimmune disease mainly affecting the connective tissue. In SSc patients, monocytes are increased in circulation, infiltrate affected tissues, and show a pro-inflammatory activation status, including the so-called interferon (IFN) signature. We previously demonstrated that the dysregulation of the IFN response in SSc monocytes is sustained by altered epigenetic factors as well as by upregulation of the long non-coding RNA (lncRNA) NRIR. Considering the enormously diverse molecular functions of lncRNAs in immune regulation, the present study investigated the genome-wide profile of lncRNAs in SSc monocytes, with the aim to further unravel their possible role in monocyte dysregulation and disease pathogenesis. Transcriptomic data from two independent cohorts of SSc patients identified 886 lncRNAs with an altered expression in SSc monocytes. Differentially expressed lncRNAs were correlated with neighboring protein coding genes implicated in the regulation of IFN responses and apoptotic signaling in SSc monocytes. In parallel, gene co-expression network analysis identified the lncRNA PSMB8-AS1 as a top-ranking hub gene in co-expression modules implicated in cell activation and response to viral and external stimuli. Functional characterization of PSMB8-AS1 in monocytes demonstrated that this lncRNA is involved in the secretion of IL-6 and TNFα, two pivotal pro-inflammatory cytokines altered in the circulation of SSc patients and associated with fibrosis and disease severity. Collectively, our data showed that lncRNAs are linked to monocyte dysregulation in SSc, and highlight their potential contribution to disease pathogenesis.

(PCGs) [30], in cis correlation analysis was performed to identify putative target genes of differentially expressed lncRNAs in the Definite cohort. To this end, the expression levels of differentially expressed lncRNAs were correlated with PCGs located 5 kb upstream or downstream of each lncRNA gene ( Figure 1A). Out of 886 differentially expressed, 278 lncRNAs were significantly correlated with PCGs localized in cis (Spearman's rho >0.4 or <−0.4, and p-value ≤ 0.05), allowing for the identification of 332 lncRNA-PCG pairs. Functional enrichment analysis of the correlated PCGs identified pathways associated with IFN response, negative regulation of apoptosis, and inflammatory cell apoptotic processes ( Figure 1B), indicating that lncRNAs altered in SSc monocytes potentially regulate genes involved in these pathways. Figure 1. Differentially expressed lncRNAs are correlated with cis localized protein coding genes relevant for SSc pathogenesis. (A) Schematic overview of the cis correlation approach to identify neighboring long non-coding RNAs (lncRNA, red) and protein coding genes (PCGs, green). Arrows indicate the direction of transcription. (B) GO-term enrichment analysis results of PCGs identified in the cis correlation analysis in the Definite SSc cohort. GO terms for significantly enriched biological processes are given on the y-axis. Bars depict the number of genes identified within the enriched pathway (N genes, bottom x-axis), dashed line indicates B&H corrected p-value of the enrichment (p-value, top x-axis). (C) Correlation coefficients (Spearman's Rho) between lncRNAs and PCGs in the Definite (x-axis) and Non-cutaneous (y-axis) cohorts. Each dot represents a lncRNA-PCG pair. Grey dots represent pairs significantly correlated in the Definite cohort only, while red dots represent pairs significantly correlated in both cohorts (Spearman's rho >0.4 or <−0.4, and p-value ≤ 0.05). (D) GO-term enrichment analysis results of protein coding genes from the in cis correlation analysis replicated the Non-cutaneous cohort.
In order to substantiate these results and identify putative lncRNA-PCG pairs involved already at the early stages of SSc development, we repeated the cis correlation analysis in an additional cohort comprising ncSSc patients and SSc patients at the early Figure 1. Differentially expressed lncRNAs are correlated with cis localized protein coding genes relevant for SSc pathogenesis. (A) Schematic overview of the cis correlation approach to identify neighboring long non-coding RNAs (lncRNA, red) and protein coding genes (PCGs, green). Arrows indicate the direction of transcription. (B) GO-term enrichment analysis results of PCGs identified in the cis correlation analysis in the Definite SSc cohort. GO terms for significantly enriched biological processes are given on the y-axis. Bars depict the number of genes identified within the enriched pathway (N genes, bottom x-axis), dashed line indicates B&H corrected p-value of the enrichment (p-value, top x-axis). (C) Correlation coefficients (Spearman's Rho) between lncRNAs and PCGs in the Definite (x-axis) and Non-cutaneous (y-axis) cohorts. Each dot represents a lncRNA-PCG pair. Grey dots represent pairs significantly correlated in the Definite cohort only, while red dots represent pairs significantly correlated in both cohorts (Spearman's rho >0.4 or <−0.4, and p-value ≤ 0.05). (D) GO-term enrichment analysis results of protein coding genes from the in cis correlation analysis replicated the Non-cutaneous cohort.
In order to substantiate these results and identify putative lncRNA-PCG pairs involved already at the early stages of SSc development, we repeated the cis correlation analysis in an additional cohort comprising ncSSc patients and SSc patients at the early disease stage (eaSSc), as well as individuals with Raynaud's Phenomenon (RP) ("Non-cutaneous cohort", Table 1). This analysis highlighted that 143 out of 332 correlated lncRNA-PCG pairs identified in the Definite cohort were reproduced in the Non-cutaneous cohort (Spearman's rho >0.4 or <−0.4, and p-value ≤ 0.05, Figure 1C). GO-term enrichment analysis of the PCGs identified in the replicated pairs showed a significant enrichment for GO terms related to type I IFN and apoptosis ( Figure 1D), suggesting that lncRNAs altered during SSc development are implicated in these processes. The 15 lncRNA-PCG pairs that were annotated in these biological pathways are given in Table 2.

Weighted Gene Co-Expression Network Analysis Identifies Clusters of Tightly Correlated RNAs Associated to SSc Clinical Features and Relevant Biological Processes
Next to cis regulatory lncRNAs, trans-acting lncRNAs are also emerging as important regulators of gene expression, especially at the post-transcriptional level [31]. To explore the regulatory potential of trans-acting lncRNAs in SSc monocytes, genome-wide co-expression network analysis was performed in parallel to the cis correlation analysis. Weighted gene co-expression network analysis (WGCNA) generated 18 distinct co-expression modules of highly correlated genes in the Definite cohort (Table S2). Correlation between module eigengenes (MEs) and clinical traits for SSc identified 10 co-expression modules that were significantly correlated with clinical parameters associated with SSc ( Figure 2A, Pearson correlation, p-value < 0.05). Comparison of the global ME expression across SSc patients and healthy controls showed that the honeydew1, brown4, darkturquoise, yellowgreen, darkorange2, and paleviolet modules were lower in SSc patients, while the ME expression of the darkgreen module was higher in SSc patients ( Figure 2B). Since the remaining violet, white, and blue modules did not show a distinct ME expression pattern in SSc patients, they were not considered for subsequent analysis ( Figure 2B). Next, GO-term enrichment analysis was performed to annotate the seven modules with a distinct ME expression pattern in SSc patients. The darkgreen module was linked to vesicle transport and the regulation of autophagy, while the darkorange and yellowgreen modules were linked to IkappaB-NF-kappaB signaling, myeloid cell differentiation, and protein modifications ( Figure 2C). No significant enrichment was identified for the remaining modules that were therefore not considered for further investigation ( Figure 2D). loid cell differentiation, and protein modifications ( Figure 2C). No significant enrichment was identified for the remaining modules that were therefore not considered for further investigation ( Figure 2D).  Table indicating the characteristics for selected modules (Column 1), considering correlations to clinical traits (Column 2), distinct ME expression pattern versus heathy controls (Column 3, up = higher ME in SSc, down = lower ME in SSc, and ND = ME not distinct from healthy subjects), and functional enrichment (Column 4). Modules that were selected for subsequent analysis are highlighted in bold.  Table indicating the characteristics for selected modules (Column 1), considering correlations to clinical traits (Column 2), distinct ME expression pattern versus heathy controls (Column 3, up = higher ME in SSc, down = lower ME in SSc, and ND = ME not distinct from healthy subjects), and functional enrichment (Column 4). Modules that were selected for subsequent analysis are highlighted in bold.

Identification of the lncRNA PSMB8-AS1 as a Reproducible Hub Gene Relevant for SSc and Monocyte Biology
The co-expression network analysis was repeated in the Non-cutaneous cohort, and the extent of overlap between the modules from the Definite and Non-cutaneous cohorts was assessed to identify reproducible co-expression modules ( Figure S2, Table S3). Thirteen modules from the Non-cutaneous cohort showed a significant overlap with the 3 selected modules from the Definite cohort (i.e., darkgreen, darkorange, and yellowgreen, Fisher's exact test, p-value < 0.05, Figure 3A), demonstrating that clinically relevant modules are reproducible across different cohorts of SSc monocytes, including pre-clinical SSc stages. the extent of overlap between the modules from the Definite and Non-cutaneous cohorts was assessed to identify reproducible co-expression modules ( Figure S2, Table S3). Thirteen modules from the Non-cutaneous cohort showed a significant overlap with the 3 selected modules from the Definite cohort (i.e., darkgreen, darkorange, and yellowgreen, Fisher's exact test, p-value < 0.05, Figure 3A), demonstrating that clinically relevant modules are reproducible across different cohorts of SSc monocytes, including pre-clinical SSc stages.
Next, by comparing the intramodular connectivity of shared genes across the overlapping modules ( Figure S3), we identified replicated hub genes with high connectivity within modules across both cohorts. While several lncRNAs were present in the replicated modules, PSMB8-AS1 was the only lncRNA identified among the top 25% most highly connected genes ( Figure 3B). Specifically, PSMB8-AS1 was a replicated hub gene in the darkgreen and darkmagenta modules from the Definite and Non-cutaneous networks, respectively, of which overlapping genes are enriched in genes related to immune cell activation and to response to virus and external stimulus ( Figure 3C). Overall, these results pointed to PSMB8-AS1 as a potential central player in the regulation of these molecular pathways that are also relevant processes for monocyte activation in SSc.  Next, by comparing the intramodular connectivity of shared genes across the overlapping modules ( Figure S3), we identified replicated hub genes with high connectivity within modules across both cohorts. While several lncRNAs were present in the replicated modules, PSMB8-AS1 was the only lncRNA identified among the top 25% most highly connected genes ( Figure 3B). Specifically, PSMB8-AS1 was a replicated hub gene in the darkgreen and darkmagenta modules from the Definite and Non-cutaneous networks, respectively, of which overlapping genes are enriched in genes related to immune cell activation and to response to virus and external stimulus ( Figure 3C). Overall, these results pointed to PSMB8-AS1 as a potential central player in the regulation of these molecular pathways that are also relevant processes for monocyte activation in SSc.

Characterization of PSMB8-AS1 Expression in SSc and Healthy Monocytes
Both the cis correlation and co-expression network analysis highlighted PSMB8-AS1 as a putative key regulator of biological processes relevant for SSc pathogenesis and monocyte activation. PSMB8-AS1 expression was significantly upregulated in ncSSc and dcSSc patients in the Definite cohort ( Figure 4A). A trend for PSMB8-AS1 upregulation was also observed in the lcSSc patients in the Definite cohort and in eaSSc and ncSSc patients of the Non-cutaneous cohort ( Figure 4A,B). To confirm the altered expression of PSMB8-AS1 across SSc patients, its expression was evaluated by a target specific RT-qPCR in an additional, independent cohort of SSc patients (Replication cohort). A statistically significant altered expression of PSMB8-AS1 was confirmed in SSc patients with the earliest symptoms (eaSSc) and the most severe phenotype (dcSSc patients) (p-value < 0.05, Figure 4C). In the other SSc groups, only some individuals displayed the upregulation of this lncRNA, but this, however, was not significantly changed when considering the whole group.

Characterization of PSMB8-AS1 Expression in SSc and Healthy Monocytes
Both the cis correlation and co-expression network analysis highlighted PSMB8-AS1 as a putative key regulator of biological processes relevant for SSc pathogenesis and monocyte activation. PSMB8-AS1 expression was significantly upregulated in ncSSc and dcSSc patients in the Definite cohort ( Figure 4A). A trend for PSMB8-AS1 upregulation was also observed in the lcSSc patients in the Definite cohort and in eaSSc and ncSSc patients of the Non-cutaneous cohort ( Figure 4A,B). To confirm the altered expression of PSMB8-AS1 across SSc patients, its expression was evaluated by a target specific RT-qPCR in an additional, independent cohort of SSc patients (Replication cohort). A statistically significant altered expression of PSMB8-AS1 was confirmed in SSc patients with the earliest symptoms (eaSSc) and the most severe phenotype (dcSSc patients) (p-value < 0.05, Figure 4C). In the other SSc groups, only some individuals displayed the upregulation of this lncRNA, but this, however, was not significantly changed when considering the whole group.  To identify the molecular pathways leading to PSMB8-AS1 induction in SSc, its expression was assessed in healthy human monocytes cultured for 2, 5, and 18 h in the presence of LPS (TLR4 ligand), R848 (TLR7/8 ligand), IFNα, and TGFβ, all stimuli linked to monocyte activation and fibrosis in SSc [32][33][34][35]. PSMB8-AS1 expression was strongly induced by LPS, R848, and IFNα, especially after 18 h of stimulation ( Figure 4D, p-value < 0.05 compared to untreated cells). In contrast, PSMB8-AS1 expression was not altered by TGFβ treatment (Figure 4E; control for TGFβ stimulation is given in Figure S4). All together, these results demonstrated that PSMB8-AS1 expression can be induced by pro-inflammatory stimuli relevant for SSc and possibly indicate its implication in the regulation of the downstream pathways, as previously demonstrated for other lncRNAs induced by pro-inflammatory stimuli [20].
Next, as the intracellular localization of lncRNAs partially governs their function [36], the specific cellular location of PSMB8-AS1 was determined by subcellular fractionation of healthy monocytes ( Figure 4F). The percentage of enrichment of IL-8 Primary Transcript and IL-8 mRNA were used as positive controls to verify the appropriate separation of different compartments, as they are expected to be found in the nucleus/chromatin or cytoplasm, respectively. In agreement with previous findings in epithelial, glioma, and pancreatic cells [37][38][39], PSMB8-AS1 was found to be enriched in monocyte cytoplasm ( Figure 4F, third panel), suggesting that this lncRNA may regulate expression of its target(s) at the post-transcriptional level [29].

Characterization of PSMB8-AS1 Function in Monocytes
To further investigate the biological relevance of PSMB8-AS1, its expression was efficiently silenced in resting and R848 stimulated monocytes using a specific small interfering RNA (siRNA) ( Figure 5A, p < 0.001). The impact of PSMB8-AS1 silencing on monocyte apoptosis was first investigated, as this lncRNA has previously been reported as a negative regulator of apoptotic signaling [37][38][39], and apoptotic signaling pathways were also identified by our GO-term enrichment analysis in the cis correlation analysis. However, FACS analysis demonstrated that PSMB8-AS1 silencing in monocytes did not affect apoptosis, neither in resting nor stimulated conditions (p-value > 0.05, Figure S5A,B). Moreover, PSMB8-AS1 silencing also did not affect the expression of its cis correlated PCG PSMB8 ( Figure S5C), showing that this lncRNA does not function via cis regulatory mechanisms, at least in the conditions tested here.
As the co-expression network analysis identified PSMB8-AS1 as a potential trans-acting hub lncRNA implicated in the regulation of immune cell activation and vesicle-related transport, its role in the secretion of pro-inflammatory cytokines was subsequently investigated. To this end, the impact of PSMB8-AS1 silencing was assessed on IL-6, IL-8, and TNFα secretion, three cytokines released by activated monocytes, and at higher levels by stimulated SSc monocytes [15,16]. ELISA analysis showed a significant decrease of IL-6 and TNFα protein concentrations in the cell-free supernatant of si-PSMB8-AS1 monocytes stimulated with R848 for 18 h as compared to monocytes transfected with scramble siRNA ( Figure 5B,C), whereas IL-8 protein levels were not affected ( Figure 5D). Overall, these results demonstrated that PSMB8-AS1 regulates the secretion of specific cytokines in TLR-activated monocytes, and can thereby contribute to monocyte activation in SSc.

Discussion
lncRNAs have been described as pivotal biological regulators of immune responses [19], and previous studies have suggested a role for lncRNAs in SSc [20][21][22][23][24]. However, a thorough characterization of the lncRNA profile in specific cell subsets obtained from SSc patients is still lacking. The aim of this study was to investigate the potential role of lncRNAs in the dysregulation of SSc monocytes, an important cell type implicated in SSc pathogenesis. Cis correlation and co-expression network analysis demonstrated that numerous lncRNAs are altered in SSc monocytes and are potentially implicated in the regulation of various biological processes relevant for SSc pathogenesis.
We have previously shown that the lncRNA NRIR plays an important role in IFN responses of SSc monocytes [20]. The results reported in the present study reinforce the concept that, besides NRIR, other lncRNAs might be strongly implicated in the regulation of IFN responses in these cells. Indeed, we identified a correlation between RP11-326I11.3, a lncRNA highly expressed in both the Definite and Non-cutaneous cohort, and IRF2, an important regulator of IFN responses [40]. Interestingly, a correlation between these two genes has previously been described in brain tissue, where they were linked to IFN signaling in central nervous system homeostasis [41]. Currently available clinical and molecular data suggest that type I IFN dysregulation is a major contributor to SSc pathogenesis [34], and the augmented expression of Type I IFN inducible genes in SSc monocytes has been linked to inflammation and the development of fibrosis [12]. Since our results showed that multiple lncRNAs are potentially involved in the amplified IFN signaling in SSc monocytes, further investigations should address the possible implication of these molecules in the disease pathogenesis.
Next to IFN signaling, cis correlated lncRNA-PCG pairs also linked this class of regulators to the negative regulation of apoptosis, a relevant process that, if inhibited, can potentially contribute to the increased numbers of circulating monocytes observed in SSc patients. RP11-356I2.4, for example, was strongly correlated with TNFAIP3, an apoptosis regulator that is induced through NF-κB signaling [42]. A correlation between

Discussion
lncRNAs have been described as pivotal biological regulators of immune responses [19], and previous studies have suggested a role for lncRNAs in SSc [20][21][22][23][24]. However, a thorough characterization of the lncRNA profile in specific cell subsets obtained from SSc patients is still lacking. The aim of this study was to investigate the potential role of lncRNAs in the dysregulation of SSc monocytes, an important cell type implicated in SSc pathogenesis. Cis correlation and co-expression network analysis demonstrated that numerous lncRNAs are altered in SSc monocytes and are potentially implicated in the regulation of various biological processes relevant for SSc pathogenesis.
We have previously shown that the lncRNA NRIR plays an important role in IFN responses of SSc monocytes [20]. The results reported in the present study reinforce the concept that, besides NRIR, other lncRNAs might be strongly implicated in the regulation of IFN responses in these cells. Indeed, we identified a correlation between RP11-326I11.3, a lncRNA highly expressed in both the Definite and Non-cutaneous cohort, and IRF2, an important regulator of IFN responses [40]. Interestingly, a correlation between these two genes has previously been described in brain tissue, where they were linked to IFN signaling in central nervous system homeostasis [41]. Currently available clinical and molecular data suggest that type I IFN dysregulation is a major contributor to SSc pathogenesis [34], and the augmented expression of Type I IFN inducible genes in SSc monocytes has been linked to inflammation and the development of fibrosis [12]. Since our results showed that multiple lncRNAs are potentially involved in the amplified IFN signaling in SSc monocytes, further investigations should address the possible implication of these molecules in the disease pathogenesis.
Next to IFN signaling, cis correlated lncRNA-PCG pairs also linked this class of regulators to the negative regulation of apoptosis, a relevant process that, if inhibited, can potentially contribute to the increased numbers of circulating monocytes observed in SSc patients. RP11-356I2.4, for example, was strongly correlated with TNFAIP3, an apoptosis regulator that is induced through NF-κB signaling [42]. A correlation between RP11-356I2.4 and TNFAIP3 has also previously been observed in the inflammatory skin disorder chronic actinic dermatitis, where both genes were downregulated in comparison to healthy controls [43]. As actinic dermatitis, like SSc, is also characterized by inflamed and thickened skin, lncRNAs involved in monocyte apoptotic processes could potentially contribute to fibrotic processes in the skin.
Both the cis correlation analysis and co-expression network analysis pointed at the lncRNA PSMB8-AS1 as a key regulator of altered gene-expression in SSc monocytes. Although this lncRNA has previously been described in the context of influenza infection [44] and cancer [37][38][39], our results link PSMB8-AS1 dysregulation with autoimmunity for the first time. The cis correlation analysis predicted PSMB8, the protein coding gene located antisense to PSMB8-AS1, as a potential target of this lncRNA; however, the two genes were not present in the same co-expression modules, and PSMB8 mRNA levels were not altered upon PSMB8-AS1 silencing. These two genes are therefore not functionally related within SSc monocytes, while their expression correlation is probably the result of a shared transcriptional program given that their promoters share a binding site for the IFN inducible transcription factor IRF1 [45]. Consistently, the protein PSMB8 can be directly regulated by IRF1 [46], and treatment of monocytes with IFNα and R848 (a TLR7/8 ligand activating IFN signaling) induces PSMB8-AS1 expression, an event likely mediated by IRF1. These findings demonstrated that IFN-mediated activation is, at least partially, responsible for the upregulation of PSMB8-AS1 in SSc monocytes, in agreement with previous studies showing that this lncRNA can be induced by influenza virus infections (triggering IFN signaling) in other human cells [44].
In line with previous observations of three other studies in pancreatic epithelial cell lines and glioma [37][38][39], our subcellular localization analysis demonstrated that PSMB8-AS1 is present in the cytoplasm of CD14+ monocytes. Given that cytoplasmic lncRNAs can regulate secretory and extracellular vesicles [47], and PSMB8-AS1 was present in coexpression modules annotated to vesicle related transport, it is possible that this lncRNA is involved in the regulation of cytoplasmic vesicles. A more detailed characterization by RNA-FISH of the precise subcellular compartments involved in intracellular vesicle transport (for example, the endoplasmic reticulum or Golgi) [36,48] could provide better insights into the exact molecular processes related to PSMB8-AS1 activity. In addition, since different stimuli can modify the subcellular compartmentalization of lncRNAs [49], it would be interesting to verify whether factors relevant for SSc pathogenesis (e.g., TLRagonists, IFNα, and CXCL4) can influence PSMB8-AS1 localization. The evidence that IL-6 and TNFα protein secretion was repressed by PSMB8-AS1 silencing in monocytes demonstrated that this lncRNA is involved in the positive regulation of these cytokines and possibly links its cytoplasmic localization to the control of cytokine secretion. In contrast, IL-8 release was unaffected, suggesting a possible cytokine-specific action of PSMB8-AS1. This is not surprising, given that specific cytokines can be secreted through distinct pathways in macrophages, adapted to suit specific stimulatory conditions [50]. Alternatively, the different impact of PSMB8-AS1 silencing on IL-6/TNFα and IL-8 secretion could be explained by distinct timeframes in the synthesis/expression of these factors or could be dependent on the specific stimuli used in our experiments. More detailed molecular investigations, exploiting a wide range of pro-inflammatory mediators, as well as stimulation times, are required to unravel the exact role of PSMB8-AS1 in cytokine release.
The positive regulation of IL-6 and TNFα secretion by PMSB8-AS1 directly link this lncRNA to SSc pathogenesis. Elevated levels of IL-6 and TNFα have indeed been reported in the serum of SSc patients in several studies [51,52], and both cytokines are associated with fibrotic processes [53,54], disease progression, and the occurrence of interstitial lung disease [55,56], especially in dcSSc patients. Most importantly, different studies proposed monocytes as the potential source of these cytokines in SSc [15,16,57]. Interestingly, the highest upregulation of PSMB8-AS1 was observed in monocytes from dcSSc patients, where the increase of IL-6 and TNFα levels and the extent of skin fibrosis are the most severe. Even if the upregulation of PSMB8-AS1 was fluctuating in other disease subsets, possibly due to the extreme clinical heterogeneity of these patients, the upregulation of this lncRNA was also confirmed in early and ncSSc patients. Such evidence indicates that PSMB8-AS1 might be involved in the regulation of inflammatory processes from early stages of the disease onward, affecting the secretion of cytokines that eventually contribute to the development and perpetuation of fibrosis.
In conclusion, in-depth bioinformatics analysis unraveled numerous lncRNAs dysregulated in SSc monocytes and highlighted an important regulatory potential of these molecules both in immune activation and disease pathogenesis. Among these, we specifically discovered that the upregulation of PSMB8-AS1 can modulate the secretion of proinflammatory cytokines by monocytes, thereby potentially contributing to the increased activation of these cells in SSc. A more detailed understanding of lncRNAs and their contribution to disease pathogenesis could provide steppingstones for the identification of novel molecular targets for manipulating monocyte activity in SSc, in order to contrast disease onset and progression.

Patient Demographics
Transcriptomic data from SSc patients and matched healthy controls from the Definite and Non-cutaneous cohorts were obtained from the gene expression omnibus (GSE124075). For the replication cohort, peripheral blood samples from SSc patients and age/sex matched healthy controls were obtained from the University Medical Center Utrecht. All participants enrolled in the study signed an informed consent form approved by the local institutional review boards prior to inclusion in this study (METC no. 12-466C, approved 2 October 2012), adherent of the Declaration of Helsinki Principles. Samples and clinical information were treated anonymously immediately after collection. SSc patients fulfilled the ACR/EULAR classification criteria [2] and were classified according to the extent of their skin fibrosis as lcSSc or dcSSc patients. Patients that fulfilled the classification criteria but did not present skin fibrosis are referred to as ncSSc patients throughout the manuscript. Finally, we also included eaSSc patients with Raynaud's Phenomenon and positivity for SSc-specific autoantibodies and/or typical nailfold capillaroscopy patterns, as defined by LeRoy et al. [3]. The demographics and clinical characteristics of the subjects enrolled in these cohorts are reported in Table 1. Ongoing treatment regimens are reported in Table S4.

Purification and Culture of CD14+ Monocytes from Healthy Control Blood and Buffy Coats
PBMCs were isolated from whole heparinized blood samples from SSc patients and healthy controls or from the buffy coats of healthy controls, by density gradient centrifugation using Ficoll-Paque Plus (GE Healthcare, Chicago, IL, USA). CD14+ monocytes were purified from PBMCs using the MACS Human Monocyte Isolation Kit II (Miltenyi Biotec, Bergisch Gladbach, Germany) on the autoMACs Pro Separator (Miltenyi Biotec) according to the manufacturer's instructions. For subsequent analysis, only cell preparations with more than 95% purity (measured by FACS analysis) for CD14+ cells were used.

RNA Purification
Total RNA was purified using the DNA/RNA/miRNA Universal kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. DNAse treatment (RNAse Free DNase I set, Qiagen) on column was performed. RNA was quantified with the Qubit ® RNA Assay Kit (Life Technologies, Carlsbad, CA, USA) on the Qubit ® Fluorometer (Invitrogen, Carlsbad, CA, USA) or on the Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA).

RNA-Sequencing Analysis
Raw sequencing data was obtained from GSE124075 [17]. Sequencing reads were aligned to the GrCh38 reference human genome (Genome Reference consortium) and the H. sapiens transcriptome (Ensembl, version 77) using TopHat [58]. Summed exon read counts per gene were estimated using the HTSeq-count function provided in the HTSeq python package [59] (v. 0.6.1p1). Differential expression analysis was performed using the negative binomial distribution-based method implemented in DESeq2 [60] (v. 1.6.3), and pair wise comparisons between SSc patients and HC groups were tested using the Wald test. Genes with a log 2 (FC) ≥0.58 or ≤−0.58 and a p-value ≤ 0.05 were considered significantly modulated. Normalized gene expression levels were expressed as variance stabilized data (VSD), calculated according to DESeq2 instructions. Gene types were annotated according to the Ensembl 77 database.

In Cis Correlation Analysis
Protein-coding genes (PCGs) that were localized within a region of 5 kb upstream or 5 kb downstream, regardless of the sense of transcription, for each differentially expressed lncRNA were recovered using the Biomart tool available on the Ensembl website. A Spearman rank-order correlation analysis of the expression of the lncRNAs and associated PCGs was then performed. Correlations with Spearman's Rho ≥0.4 or ≤−0.4, and p-value ≤ 0.05 were considered significant. Correlation analysis was performed using the rcorr function implemented in the Hmisc package in R [61].

GO-Term Enrichment Analysis
Gene ontology (GO) enrichment analysis was performed using ToppFun [62]. Enrichment of biological process (BP) associated GO terms was tested using the probability density function. p-value was adjusted according to Benjamini-Hochberg/FDR correction. BP terms significantly enriched (B&H corrected p-value ≤ 0.05) were considered.

Weighted Gene Co-Expression Network Analysis (WGCNA)
Weighted gene co-expression networks were constructed for the Definite and Noncutaneous cohorts separately using the R package WGCNA [63]. The VSD data of all genes with at least 1 count in all samples were used as input. An unsigned network with a scalefree topology was constructed, using a soft threshold power β = 13 for the Definite cohort and a soft threshold power β = 4 for the Non-cutaneous network. Modules were identified using the cutreeDynamic function with a minimum module size of 30. Closely related modules were merged using the mergeCloseModules function (cutHeight = 0.25). Gene expression profiles across the modules were summarized into module eigengene (ME) values based on the first principal component of each module. Fisher's exact test was used to calculate the extent of module overlap between the Definite and Non-cutaneous networks, as previously described [64]. Intramodular connectivity (i.e., the connectivity of a gene to genes nodes within the same module) was obtained using the intramodularConnectivity function from the WGCNA package.

Reverse Transcription Quantitative Real-Time PCR (RT-qPCR)
Purified RNA (200-1000 ng) was reverse transcribed using the SuperScript ® III Reverse Transcriptase kit (Invitrogen), according to the manufacturer's instructions. Gene expression was quantified, in duplicate, by RT-qPCR using 9 ng cDNA with SYBR Select Master Mix (Applied Biosystems, Foster City, CA, USA), in the presence of 400 nM genespecific primers (Table 3), on the ViiA™ 7 Real-Time PCR System (Applied Biosystems). Relative expression of each gene was determined according to the comparative CT (∆∆CT) method using RPL32 as an endogenous control (where the ∆CT equals the CT of the mRNA of interest-the CT of RPL32) [65]. Table 3. Gene specific primer pairs used for RT-qPCR.

Assessment of Cytokine Levels Using ELISA
Concentrations of IL-6, TNFα, and IL-8 in cell-free supernatants from cultured CD14+ monocytes were measured by the sandwich enzyme linked immunosorbent assay (ELISA). IL-6 and IL-8 were quantified using the PeliKine compact human IL-6 and IL-8 ELISA kits (Sanquin Reagents, Amsterdam, The Netherlands), and TNFα was quantified using the Human TNF-α ELISA Set (Diaclone, Besançon, France), according to the manufacturer's instructions.

Statistical Analysis
Data are expressed as mean ± SEM unless otherwise indicated. Unless indicated otherwise, analysis of differences was performed using the Mann Whitney test. For multiple group comparisons, the one-or two-way analysis of variance (ANOVA) was used. p-values < 0.05 were considered statistically significant. Figures were produced using the R package ggplot2 [66] or the GraphPad Prism software (v 8.3, www.graphpad.com, GraphPad Software, Inc., San Diego, CA, USA).

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/ijms22094365/s1. Figure S1: Expression of lncRNAs is altered in SSc monocytes. Figure S2: Module overlap between gene co-expression networks identified in the Definite and Non-Cutaneous cohorts. Figure S3: Connectivity overlap of selected reproduced co-expression modules. Figure S4: MMP2 expression is induced by TGFβ signaling in monocytes. Figure S5: Silencing of PSMB8-AS1 does not affect apoptosis or PSMB8 expression in healthy monocytes. Table S1: Differentially expressed lncRNA identified in SSc monocytes. Table S2: Co-expression network analysis in the Definite cohort. Table S3: Module overlap analysis between the Definite and Non-cutaneous cohorts.  Institutional Review Board Statement: The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the medical ethics committee of the University Medical Center Utrecht (METC no. 12-466C, approved 2 October 2012). All patients gave their written informed consent in accordance with the declaration of Helsinki.
Informed Consent Statement: Informed consent was obtained from all subjects included in the study.
Data Availability Statement: RNA-sequencing data was obtained from the gene expression omnibus (GSE124075). All other relevant raw data from the data presented in the manuscript or the supplementary figures and tables are available by the authors of the study upon request.