Abstract

   Primary Sjögren’s syndrome (pSjS) is a chronic systemic autoimmune
   disorder, primarily affecting exocrine glands; its pathogenesis is
   still unclear. Long non-coding RNAs (lncRNAs) are thought to play a
   role in the pathogenesis of autoimmune diseases and a comprehensive
   analysis of lncRNAs expression in pSjS is still lacking. To this aim,
   the expression of more than 540,000 human transcripts, including those
   ascribed to more than 50,000 lncRNAs is profiled at the same time, in a
   cohort of 16 peripheral blood mononuclear cells PBMCs samples (eight
   pSjS and eight healthy subjects). A complex network analysis is carried
   out on the global set of molecular interactions among modulated genes
   and lncRNAs, leading to the identification of reliable
   lncRNA-miRNA-gene functional interactions. Taking this approach, a few
   lncRNAs are identified as targeting highly connected genes in the pSjS
   transcriptome, since they have a major impact on gene modulation in the
   disease. Such genes are involved in biological processes and molecular
   pathways crucial in the pathogenesis of pSjS, including immune
   response, B cell development and function, inflammation, apoptosis,
   type I and gamma interferon, epithelial cell adhesion and polarization.
   The identification of deregulated lncRNAs that modulate genes involved
   in the typical features of the disease provides insight in disease
   pathogenesis and opens avenues for the design of novel therapeutic
   strategies.

   Keywords: primary Sjögren’s syndrome, long non-coding RNA, signaling
   pathway, protein–protein (PPI) network, gene module

1. Introduction

   Sjögren’s syndrome (SjS) is a chronic systemic autoimmune disease of
   still unknown origin, that may be present either alone (defined as
   primary SjS: pSjS) or associated with other autoimmune diseases
   (defined as secondary SjS). The prevalence of pSjS in the general
   population is between 0.01 and 0.1% with a higher prevalence of the
   disease in women [[34]1,[35]2]. The clinical symptoms of pSjS include
   exocrinopathy, resulting in dry eyes and mouth and extraglandular
   systemic manifestations such as arthralgias, fatigue, vasculitis,
   pulmonary fibrosis and pulmonary hypertension, interstitial nephritis
   and central and peripheral involvement of the nervous system [[36]3].

   Although little is known on the pathogenesis of the disease, several
   factors have been shown to contribute to its onset, such as genetic
   background, environmental factors including viral infections [[37]4]
   and epigenetic mechanisms, such as microRNAs (miRNAs).

   Our attention is focused on the epigenetic mechanisms that may be
   involved in pSjS pathogenesis, through the analysis of the modulation
   of long non-coding RNAs (lncRNAs) expression in pSjS patients.

   The identification of lncRNAs modulated in the disease and an
   integrated analysis of lncRNAs, miRNAs and gene expression profiles in
   patients affected by pSjS are reported here as an historical first, to
   the best of the authors’ knowledge. Interestingly, the identified
   lncRNAs are able to modulate pathogenetically relevant molecular
   pathways of the disease.

2. Materials and Methods

2.1. Patients

   Eight female patients (mean age 52 ± 15 years old) affected by pSjS,
   attending the Unit of Autoimmune Diseases at the University Hospital of
   Verona, Northern Italy, and 8 sex and age matched healthy controls were
   enrolled. Both patients and controls were subjects of Caucasian origin
   from Northern Italy. All patients enrolled in this study were diagnosed
   according to the American College of Rheumatology (ACR)–European League
   against Rheumatism (EULAR) criteria [[38]5]. They underwent labial
   salivary gland biopsy and disease activity was evaluated at the moment
   of enrollment in the study.

   A written informed consent was obtained from all the participants of
   the study and the study protocol was approved by the Ethical Committee
   of the Azienda Ospedaliera Universitaria Integrata di Verona
   (identification code 1538, version 3, date of approval 12 March 2012).
   All the investigations have been performed according to the principles
   contained in the Helsinki declaration.

2.2. Microarray Analysis

   Blood sample collection was carried out using BD Vacutainer K2EDTA
   tubes (Becton Dickinson, Franklin Lakes, NJ, USA) and 21-gauge needles.

   PBMCs isolation was performed by Ficoll–HyPaque (Pharmacia Biotech,
   Piscataway, NJ, USA) gradient centrifugation. Patients and controls had
   similar PBMCs distribution. Total RNA was extracted from PBMCs (10^7
   cells) using an miRNeasy mini kit (Qiagen GmbH, Hilden, Germany). cRNA
   preparation, samples hybridization and scanning were performed
   following the Affymetrix (Affymetrix, Santa Clara, CA, USA) provided
   protocols using a Cogentech Affymetrix microarray unit (Campus IFOM
   IEO, Milan, Italy). All samples were hybridized on a Human Clariom D
   (Thermo Fisher Scientific, Waltham, MA, USA) gene chip. Signal
   intensities were analysed with the Transcriptome Analysis Console (TAC)
   4.0 software (Applied Biosystem, Foster City, CA, USA by Thermo Fisher
   Scientific, Waltham, MA, USA).

   Using the Human Clariom D arrays, more than 540,000 human transcripts
   were interrogated, starting from as little as 100 pg of total RNA.
   Signal 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.

   Target annotations of long non-coding RNAs were retrieved using
   starBase v2.0 ([39]http://starbase.sysu.edu.cn/) where lncRNAs
   interactions, experimentally validated by high-throughput experimental
   technologies, are registered [[40]6].

   The list of gene targets of microRNAs (miRNAs) that were targeted by
   lncRNAs were inferred interrogating the FunRich database
   ([41]www.funrich.org/) [[42]7].

2.3. 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
   ([43]http://string-db.org/) [[44]8].

   Network topological analysis was performed using the Cytoscape software
   (version 3.7, ([45]www.cytoscape.org) [[46]9].

   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).

2.4. Gene Functional Classification and Enrichment Analysis

   Genes were classified functionally into Biological Processes (BPs)
   according to the Gene Ontology (GO) project annotations
   ([47]www.geneontology.org) by the Panther expression analysis tools
   ([48]http://pantherdb.org/) [[49]10].

   Pathway classification and enrichment (hypergeometric p-value ≤ 0.05)
   analysis were achieved with FunRich.

2.5. Real Time PCR of lncRNAs

   Five hundred 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). Transcripts of relative expression levels were obtained after
   normalization against the geometric mean of the housekeeping genes
   Glyceraldehyde 3-phosphate dehydrogenase GAPDH and beta-actin (ACTB)
   expression. The ΔΔCt method of relative quanrtification was used for
   comparing relative fold expression differences. Results are expressed
   as fold changes with respect to healthy subjects.

2.6. Real Time PCR of Genes Modulated in pSjS 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-Demande Gene
   Expression Products service (Applied Biosystems).

   Real-Time PCR reactions were carried out in a two-tube system and in
   singleplex. The Real- Timeamplifications encompassed 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 a 7500 Sequence Detector (Applied Biosystems). Signals
   were detected following the manufacturer’s instructions. This
   methodology allowed the identification of the cycling point where the
   PCR product was detectable by means of fluorescence emission (Threshold
   cycle or Ct value). The Ct value correlated to the quantity of target
   mRNA. Relative expression levels were calculated for each sample after
   normalization against the housekeeping genes GAPDH, beta-actin and 18 s
   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. Each amount of RNA tested had
   triplicate Ct values averaged. Since Ct values varied linearly with the
   logarithm of the amount of RNA, this average represented a geometric
   mean.

2.7. Real Time PCR of microRNAs

   miRNAs expression was evaluated 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 a TaqMan^®
   Advanced miRNA cDNA synthesis kit according to the 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 was loaded in PCR. Twenty µL PCR reactions were
   composed by 2× Fast Advanced Master Mix and TaqMan^® Advanced miRNA
   assays for miR-30e-5p, miR-32-5p, miR-155-5p, miR-195-5p, miR-378a-3p
   and miR-30b-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 was
   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 employing
   QuantStudio Real-Time PCR system software v. 1.3.

3. Results

3.1. Patients Characteristics

   All the patients enrolled in the study had antinuclear antibodies
   >1:640, nuclear dots or homogeneous patterns, and were positive for
   anti-ENA antibodies SSA/Ro; four of them also were positive for
   anti-SSB antibodies. Labial salivary gland histopathology showed a
   focal lymphocytic sialadenitis with a focus score >1 [[50]11]. The
   median Eular Sjogren’s syndrome disease activity index (ESSDAI) score
   was 6. The median Eular Sjogren’s syndrome Patient Reported Index
   (ESSPRI) score was 7 [[51]12]. Three patients presented germinal
   centre-like structures in their labial salivary gland biopsy and one
   patient developed a B cell lymphoma soon after enrollment in the study.

3.2. High-Throughput Gene and Long Non-Coding RNA Expression Profiling in
Peripheral Blood Mononuclear Cells of Patients with pSjS

   Intending to identify lncRNAs potentially involved in pSjS
   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 16 PBMCs samples (8 pSjS and 8 healthy
   subjects). Transcriptional profiles of pSjS patients and healthy
   subjects were compared and, after a robust filtering procedure
   (Bonferroni-corrected p-value ≤ 0.01 and fold change ≥ |1.5|), 2503
   coding-genes were modulated significantly ([52]Table S1).

   The functional classification by Gene Ontology
   ([53]http://www.geneontology.org/) of the 2503 differentially expressed
   genes (DEGs) revealed that they were involved in meaningful biological
   processes (BPs) known to play a role in the disease, including
   apoptosis, cell adhesion, immune response, type I Interferon signaling,
   Interferon-gamma signaling, extracellular matrix (ECM) organization and
   morphogenesis of a branching epithelium. A selection of genes that play
   a role in the above mentioned BPs is shown in [54]Table 1.

Table 1.

   Selection of significant differentially expressed genes in pSjS
   patients versus healthy subjects, grouped according to the biological
   processes to which they are ascribed.
   Gene Symbol Description Fold Change p-Value mRNA Accession
   Apoptosis
   CCAR2 cell cycle and apoptosis regulator 2 1.83 0.002 [55]NM_021174
   CASP10 caspase 10 2.26 0.004 [56]NM_001206524
   BCL2L13 BCL2-like 13 (apoptosis facilitator) 2.49 <0.001
   [57]NM_001270726
   BCL2L12 BCL2-like 12 (proline rich) 2.44 0.002 [58]NM_001040668
   PYCARD PYD and CARD domain containing 2.52 0.004 [59]NM_013258
   DAPK3 death-associated protein kinase 3 2.09 <0.001 [60]NM_001348
   PDCD6 programmed cell death 6 2.93 <0.001 [61]NM_001267556
   CASP9 caspase 9 2.01 <0.001 [62]NM_001229
   DEDD2 death effector domain containing 2 3.82 <0.001 [63]NM_001270614
   BMF Bcl2 modifying factor 1.88 0.003 [64]NM_001003940
   PDCD7 programmed cell death 7 1.98 0.001 [65]NM_005707
   FADD Fas (TNFRSF6)-associated via death domain 2.06 <0.001
   [66]NM_003824
   ANP32B acidic nuclear phosphoprotein 32 family member B 2.34 <0.001
   [67]NM_006401
   BRI3 brain protein I3 1.87 0.003 [68]NM_001159491
   CLPTM1L CLPTM1-like 2.8 <0.001 [69]NM_030782
   ING2 inhibitor of growth family member 2 1.97 <0.001 [70]NM_001291959
   BFAR bifunctional apoptosis regulator 1.98 0.003 [71]NM_016561
   DAP death-associated protein 2.42 0.004 [72]NM_001291963
   SORT1 sortilin 1 2.51 0.008 [73]NM_001205228
   RRAGA Ras-related GTP binding A 1.76 0.001 [74]NM_006570
   Cell adhesion
   LIN7C lin-7 homolog C (C. elegans) −1.76 0.001 [75]NM_018362
   ITGB3BP integrin beta 3 binding protein (beta3-endonexin) −2.87 0.001
   [76]NM_001206739
   PRKCE protein kinase C, epsilon −1.62 0.005 [77]NM_005400
   FER fer (fps/fes related) tyrosine kinase −2.3 0.001 [78]NM_001308028
   ADAM8 ADAM metallopeptidase domain 8 3.13 <0.001 [79]NM_001109
   ADAM15 ADAM metallopeptidase domain 15 3.59 <0.001 [80]NM_001261464
   CLDN5 claudin 5 1.69 0.007 [81]NM_001130861
   EMP2 epithelial membrane protein 2 1.71 0.005 [82]NM_001424
   ICAM3 intercellular adhesion molecule 3 2.86 <0.001 [83]NM_002162
   HN1 hematological and neurological expressed 1 2.71 <0.001
   [84]NM_001002032
   ITGA5 integrin alpha 5 2.22 <0.001 [85]NM_002205
   MISP mitotic spindle positioning 1.86 0.003 [86]NM_173481
   PTPRJ protein tyrosine phosphatase, receptor type, J 1.71 0.006
   [87]NM_001098503
   RHOC ras homolog family member C 2.5 0.004 [88]NM_001042678
   ZEB1 zinc finger E-box binding homeobox 1 −4.18 <0.001 [89]NM_001128128
   ZFYVE21 zinc finger, FYVE domain containing 21 3.07 0.001
   [90]NM_001198953
   Immune response
   LAMTOR2 late endosomal/lysosomal adaptor, MAPK and MTOR activator 2
   3.08 0.002 [91]NM_001145264
   IL6R interleukin 6 receptor 1.79 0.001 [92]NM_000565
   IGKV3D-7 immunoglobulin kappa variable 3D-7 4.17 0.004
   OTTHUMT00000476805
   IGLV3-19 immunoglobulin lambda variable 3-19 6.58 0.002
   OTTHUMT00000321830
   HLA-G major histocompatibility complex, class I, G 2.98 <0.001
   [93]NM_002127
   NCF1 neutrophil cytosolic factor 1 2.17 0.003 [94]NM_000265
   CD68 CD68 molecule 2.18 0.002 [95]NM_001040059
   DEF6 DEF6 guanine nucleotide exchange factor 1.58 0.003 [96]NM_022047
   LILRA6 leukocyte immunoglobulin-like receptor, subfamily A, member 6
   2.37 0.004 [97]NM_024318
   TRBV24-1 T cell receptor beta variable 24-1 3.57 0.001
   OTTHUMT00000352499
   IL5RA interleukin 5 receptor, alpha 2.14 <0.001 [98]NM_000564
   CD6 CD6 molecule 3.71 <0.001 [99]NM_001254750
   CD7 CD7 molecule 1.93 <0.001 [100]NM_006137
   BTK Bruton agammaglobulinemia tyrosine kinase 1.79 0.006 [101]NM_000061
   BAX BCL2-associated X protein 3.14 <0.001 [102]NM_001291428
   BAK1 BCL2-antagonist/killer 1 2.46 0.005 [103]NM_001188
   GATA3 GATA binding protein 3 1.74 0.008 [104]NM_001002295
   TAP1 antigen peptide transporter 1 2.58 <0.001 [105]NM_001292022
   TAP2 antigen peptide transporter 2 2.23 0.002 [106]NM_018833
   IL17RA interleukin 17 receptor A 2.52 0.001 [107]NM_001289905
   IL23A interleukin 23, alpha subunit p19 −2.9 0.009 [108]NM_016584
   IL4R interleukin 4 receptor 2.37 <0.001 [109]NM_000418
   CD33 CD33 molecule 3.39 0.002 [110]NM_001082618
   IL2RA interleukin 2 receptor, alpha 1.89 0.006 [111]NM_000417
   LAT linker for activation of T-cells 3.25 <0.001 [112]NM_001014987
   C1RL complement component 1, r subcomponent-like 1.52 0.010
   [113]NM_001297640
   FCGR1A Fc fragment of IgG, high affinity Ia, receptor (CD64) 3.2 <0.001
   [114]NM_000566
   CD81 CD81 molecule 2.04 0.008 [115]NM_001297649
   KLHL20 kelch-like family member 20 −1.94 0.004 [116]NM_014458
   KLHL6 kelch-like family member 6 2.18 0.005 [117]NM_130446
   ORAI1 ORAI calcium release-activated calcium modulator 1 2.97 <0.001
   [118]NM_032790
   CBX8 chromobox homolog 8 1.62 0.008 [119]NM_020649
   ENO1 enolase 1, (alpha) 1.91 0.004 [120]NM_001201483
   GNAI2 guanine nucleotide binding protein, alpha inhibiting activity
   polypeptide 2 1.67 0.007 [121]NM_001166425
   SPI1 Spi-1 proto-oncogene 1.73 0.007 [122]NM_001080547
   TRAF3 TNF receptor-associated factor 3 −1.63 0.009 [123]NM_001199427
   HNRNPL heterogeneous nuclear ribonucleoprotein L 2.16 0.001
   [124]NM_001005335
   FOSL1 FOS-like antigen 1 −2.14 0.008 [125]NM_001300844
   Type I interferon signaling
   IRF5 interferon regulatory factor 5 2.24 0.003 [126]NM_001098627
   IRF7 interferon regulatory factor 7 3.14 0.002 [127]NM_001572
   IRF9 interferon regulatory factor 9 2.68 <0.001 [128]NM_006084
   MYD88 myeloid differentiation primary response 88 1.86 0.005
   [129]NM_001172566
   HLA-H major histocompatibility complex, class I, H (pseudogene) 2.18
   0.004 [130]NR_001434
   OAS1 2-5-oligoadenylate synthetase 1 3.58 0.002 [131]NM_001032409
   IFI35 interferon-induced protein 35 1.78 0.007 [132]NM_005533
   IFITM3 interferon induced transmembrane protein 3 2.54 0.006
   [133]NM_021034
   IFNA10 interferon, alpha 10 1.99 0.002 [134]NM_002171
   ADAR adenosine deaminase, RNA-specific 3.15 0.005 [135]NM_001025107
   OAS3 2-5-oligoadenylate synthetase 3 4.27 0.008 [136]NM_006187
   STAT1 signal transducer and activator of transcription 1 1.75 0.008
   [137]NM_007315
   KLHL20 kelch-like family member 20 −1.94 0.004 [138]NM_014458
   Interferon-gamma signaling
   HLA-G major histocompatibility complex, class I, G 2.98 <0.001
   [139]NM_002127
   FCGR1B Fc fragment of IgG, high affinity Ib, receptor (CD64) 2.25 0.002
   [140]NM_001004340
   PML promyelocytic leukemia 3.18 0.002 [141]NM_002675
   STAT1 signal transducer and activator of transcription 1 1.75 0.008
   [142]NM_007315
   TRIM21 tripartite motif containing 21 2.13 0.002 [143]NM_003141
   FCGR1A Fc fragment of IgG, high affinity Ia, receptor (CD64) 3.2 <0.001
   [144]NM_000566
   Inflammatory response
   CCR4 chemokine (C-C motif) receptor 4 3.73 0.005 [145]NM_005508
   CCR8 chemokine (C-C motif) receptor 8 2.24 0.004 [146]NM_005201
   IL6R interleukin 6 receptor 1.79 0.001 [147]NM_000565
   MYD88 myeloid differentiation primary response 88 1.86 0.005
   [148]NM_001172566
   CSF1 colony stimulating factor 1 (macrophage) 2.23 0.005 [149]NM_000757
   MIF macrophage migration inhibitory factor (glycosylation-inhibiting
   factor) 2.48 <0.001 [150]NM_002415
   TNFRSF1A tumor necrosis factor receptor superfamily, member 1A 4.39
   <0.001 [151]NM_001065
   TGFB1 transforming growth factor beta 1 1.71 0.003 [152]NM_000660
   LTB4R leukotriene B4 receptor 2.46 <0.001 [153]NM_001143919
   ALOX5 arachidonate 5-lipoxygenase 2.57 0.007 [154]NM_000698
   IL23A interleukin 23, alpha subunit p19 −2.9 0.009 [155]NM_016584
   CXCR3 chemokine (C-X-C motif) receptor 3 2.25 0.002 [156]NM_001142797
   MIF macrophage migration inhibitory factor (glycosylation-inhibiting
   factor) 2.48 <0.001 [157]NM_002415
   CEBPD CCAAT/enhancer binding protein (C/EBP), delta 1.94 0.007
   [158]NM_005195
   MAP3K12 mitogen-activated protein kinase kinase kinase 12 1.88 0.005
   [159]NM_001193511
   TRIB2 tribbles pseudokinase 2 2.19 0.005 [160]NM_021643
   ECM organization
   MMP9 matrix metallopeptidase 9 3.31 0.002 [161]NM_004994
   EFEMP2 EGF containing fibulin-like extracellular matrix protein 2 2.43
   0.001 [162]NM_016938
   ADAM15 ADAM metallopeptidase domain 15 3.59 <0.001 [163]NM_001261464
   ADAM8 ADAM metallopeptidase domain 8 3.13 <0.001 [164]NM_001109
   BSG basigin (Ok blood group) 2.07 0.001 [165]NM_001728
   DAG1 dystroglycan 1 (dystrophin-associated glycoprotein 1) 1.57 0.003
   [166]NM_001165928
   CTGF connective tissue growth factor 1.54 0.006 [167]NM_001901
   TGFB1 transforming growth factor beta 1 1.71 0.003 [168]NM_000660
   LOXL3 lysyl oxidase-like 3 2.02 0.002 [169]NM_001289164
   SPOCK2 sparc/osteonectin, cwcv and kazal-like domains proteoglycan
   (testican) 2 2.2 0.006 [170]NM_001134434
   Morphogenesis of a branching epithelium
   CTSZ cathepsin Z 2.99 0.002 [171]NM_001336
   ILK integrin linked kinase 1.82 0.001 [172]NM_001014794
   DAG1 dystroglycan 1 (dystrophin-associated glycoprotein 1) 1.57 0.003
   [173]NM_001165928
   EDN1 endothelin 1 −2.17 0.002 [174]NM_001168319
   TGFB1 transforming growth factor beta 1 1.71 0.003 [175]NM_000660
   GREM1 gremlin 1 1.53 0.009 ENST00000633992
   ENG endoglin 2.42 <0.001 [176]NM_000118
   Epithelial cells polarization
   ARF6 ADP-ribosylation factor 6 1.53 0.008 [177]NM_001663
   FRMD4A FERM domain containing 4A 1.85 0.003 [178]NM_018027
   RHOQ ras homolog family member Q 1.79 0.005 [179]NM_012249
   [180]Open in a new tab

   LncRNA profiling showed a statistically significant variation
   (Bonferroni-corrected p-value ≤  0.01 and fold change  ≥ |1.5|) for 199
   long non-coding RNAs ([181]Table S2).

   Since modulated genes are well representative of BPs and pathways
   strictly connected to the disease, we decided to verify whether the
   modulated lncRNAs could be connected functionally to the pSjS
   transcriptome, thus playing a role in pSjS pathophysiology.

   The identified 199 lncRNAs were then analysed, retaining only those
   transcripts for which experimentally validated microRNA (miRNA) targets
   had already been annotated in starBase and, among these, only lncRNAs
   that targeted at least 10 miRNAs were selected.

   Using this approach, we obtained 6 lncRNAs, including CTD-2020K17.1,
   LINC00657, RP11-169K16.9, LINC00511, RP11-372K14.2 and RP11-214O1.2.
   ([182]Table 2)

Table 2.

   Selected long non-coding RNAs modulated in pSjS patients versus healthy
   subjects.
   Gene Symbol Description Fold Change p-Value Public Gene IDs
   LINC00657 long intergenic non-protein coding RNA 657 1.8 0.001
   [183]NR_027451
   LINC00511 long intergenic non-protein coding RNA 511 −2.0 0.008
   [184]NR_036488
   CTD-2020K17.1 novel transcript, antisense to FMNL1 2.7 0.000
   ENST00000585471.1
   RP11-169K16.9 uncharacterized LOC729614 1.7 0.008 [185]NR_024279
   RP11-214O1.2 uncharacterized protein MGC12916 −2.1 0.010 [186]NR_026880
   RP11-372K14.2 novel transcript, antisense to SH3D19 −2.1 0.004
   ENST00000603472.1
   Gene Symbol miRNA Targets Total Number of Targeted Modulated Genes
   Targeted Modules Total Number of Targeted Module-Associated Genes
   LINC00657 67 313 7 75
   LINC00511 11 194 7 41
   CTD-2020K17.1 11 120 7 25
   RP11-169K16.9 13 90 5 15
   RP11-214O1.2 12 86 4 13
   RP11-372K14.2 12 75 4 12
   [187]Open in a new tab

   To strengthen the significance of our analysis we wanted to verify the
   ability of the selected lncRNAs to regulate genes that are
   differentially expressed in pSjS samples through the interaction with
   their miRNA targets. We therefore evaluated the entire lists of their
   miRNA targets that were validated experimentally by high-throughput
   technologies and selected only those miRNAs that targeted genes
   modulated in pSjS patients to bona fide outline authentic interactions
   that are well established in pSjS.

   Shown in [188]Table 2, LINC00657, LINC00511 and CTD-2020K17.1 targeted
   the highest number of modulated genes ([189]Table S3). The expression
   data were confirmed by RT-PCR ([190]Figure S1).

3.3. PPI Network and Modular Analysis of Genes and lncRNAs Modulated in pSjS

   Since the modulation of highly connected genes can have a more
   pronounced effect in the disease pathogenesis than the modulation of
   genes that are not functionally connected, we evaluated whether the 6
   selected lncRNAs could target highly interacting genes in pSjS.

   Considering this, to prioritize transcripts that may have a role in the
   pSjS pathogenesis, we performed a network analysis dissecting all the
   differentially expressed genes in pSjS that were connected
   functionally. We therefore constructed a protein–protein interaction
   (PPI) network that included all the experimentally validated functional
   interactions among the protein products of the 2503 modulated genes in
   pSjS ([191]Figure 1A). The obtained network included 2500 modulated
   genes (nodes) that were connected by 14,169 pairs of interactions
   (edges) and showed a good PPI enrichment p-value (<10^−16).

Figure 1.

   [192]Figure 1
   [193]Open in a new tab

   Network analysis of differentially expressed genes (DEGs) in pSjS
   patients. (A) Protein–protein interaction (PPI) network of DEGs; (B)
   Modules originated from the interaction network.

   Moreover, since genes that are strictly connected to the disease
   phenotype display a strong tendency to cluster together in few network
   regions [[194]13], we performed a modular analysis to find areas in
   which the most highly connected genes were clustered. Using this
   approach, we could identify 7 gene modules that are most likely to be
   involved in the disease pathogenesis ([195]Figure 1B and [196]Table
   S4).

   All 7 modules were imported in Cytoscape adding to their scaffolds
   miRNAs-genes, and lncRNAs—miRNAs interactions.

   The topological analysis of such implemented modules highlighted, for
   the 6 lncRNAs, the lncRNA–gene interactions mediated by their targeted
   miRNAs. We observed that LINC00657, LINC00511 and CTD-2020K17.1
   targeted highly connected genes and the highest number of
   module-associated genes that were distributed in all the 7 modules
   ([197]Table 2 and [198]Table S5).

   We therefore narrowed our analysis to CTD-2020K17.1, LINC00511, and
   LINC00657 lncRNAs, since they most probably had a major impact on pSjS
   transcriptome (as also suggested by the network analysis), ([199]Figure
   2).

Figure 2.

   [200]Figure 2
   [201]Open in a new tab

   Functional interactions among the three selected lncRNAs and genes
   modulated in pSjS patients. Degree sorted circle layouts of the
   protein–protein interaction (PPI) network of differentially expressed
   genes in pSjS patients are shown. Genes (nodes) are ordered around a
   circle based on their degree of connectivity (number of edges).

   Noteworthy, targeted genes included several important transcripts
   involved in T cell development (GATA3), in the response to type I
   interferon (IRF5, IRF9 and KLHL20), in inflammatory response (IL6R and
   CEBPD) and in B-cell physiology and malignancy (BAK1, BAX, CBX8, ENO1,
   GNAI2, HNRNPL, LTBR, TRAF3 and WDR5).

   When we analysed the list of miRNAs regulated by the 3 lncRNAs
   ([202]Table S6), we found that 6 miRNAs played a role in B cell
   development, for example mir-17-5p, mir-20b-5p, mir-34a-5p, mir-34c-5p,
   mir-155-5p and mir-93-5p [[203]14,[204]15], and that 51 miRNAs have
   been previously associated with several types of human B cell
   lymphomas. Noteworthy, 15 of the above mentioned miRNAs were already
   associated to pSjS (see [205]Table 3).

Table 3.

   miRNA targets of the selected lncRNAs that have been associated to
   lymphoma and/or pSjS.
   MiRNAs Previously Associated to Lymphoma
   lncRNA miRNA Target References