Comprehensive expression profile of long non-coding RNAs in Peripheral blood mononuclear cells from patients with neuropsychiatric systemic lupus erythematosus

Introduction

Systemic lupus erythematosus (SLE) is a prototypic systemic autoimmune disease characterized by autoantibody production and immune complex deposition, resulting in multiple organ damages including kidneys, heart, lungs, and nervous system (1). Neuropsychiatric SLE (NPSLE) is one of the most severe manifestations of SLE and encompasses a variety of neurological and psychiatric disorders, leading to substantial morbidity and mortality (2). Complex interrelated pathogenic mechanisms have been implicated in NPSLE pathogenesis, including genetic factors, vasculopathy, neuroendocrine-immune imbalance, neuronal damage mediated by autoantibodies, inflammatory mediators and blood brain barrier dysfunction (2). Recent studies found that epigenetic alterations, such as altered DNA methylation, histone modifications, and one class of non-coding RNAs (ncRNAs) named microRNAs mediated posttranscriptional gene silencing may affect the initiation and progression of NPSLE (3). However, the potential role of other types of ncRNAs such as long non-coding RNAs (lncRNAs) in etiology and pathogenesis of NPSLE still largely unknown.

LncRNAs are generally defined as non-protein-coding transcripts larger than 200 nucleotides, located within the intergenic stretches or overlapping antisense transcripts of protein coding genes (4). Although without protein coding capability, accumulating evidence suggested that lncRNAs participate in diverse biological processes, including genomic imprinting, chromosome modification, transcriptional activation, protein stability, and miRNAs inhibition. LncRNAs were reported to be involved in pathogenesis of multiple diseases such as neurodegenerative diseases, diabetic mellitus and cancer (5-7). Recently, the roles of lncRNAs in immune responses have aroused increasing attention (8,9). LncRNAs have been found to be differentially expressed (DE) in peripheral blood mononuclear cells (PBMCs), immortalized B cells, and kidney biopsy specimens from SLE patients (10-13), and play crucial roles in autoimmune diseases including SLE (14,15). However, the precise role of lncRNAs in NPSLE pathogenesis remains unclear.

In the present study, we performed high-throughput RNA sequencing assays on PBMCs from NPSLE patients. Outstanding lncRNAs functions were annotated based on co-expression genes and a gene ontology (GO) biological analysis process. The correlations of lncRNAs with disease activity indicators, system and visceral involvement indicators including the nervous system of NPSLE patients were revealed through correlation analyses.

Methods

Study subjects

Totally 57 active SLE patients and 30 age-matched healthy subjects with no history of autoimmune diseases were included. All patients fulfilled the 1997 American College of Rheumatology (ACR) revised criteria for SLE (16). SLE disease activity was recorded by Systemic Lupus Erythematosus Disease Activity Index 2000 (SLEDAI) score at the time of blood draw. Patients with SLEDAI ≥8 were categorized as active SLE patients. Characteristically, active SLE patients can be classified into two groups, NPSLE group are 26 patients with acute symptoms of NPSLE, and Non-NPSLE group are 31 SLE patients who never experienced neuropsychiatric manifestations (17). NPSLE was diagnosed if the patent had at least one of 19 neuropsychiatric syndromes according to ACR classification and no other explanation related with SLE. Demographic and clinical characteristics of these patient samples are shown in the Table S1. All subjects were given informed consents for the collection of peripheral blood. This study was approved by the Ethics Committee of the Affiliated Drum Tower Hospital of Nanjing University Medical School (ID: SC201700201) and was conducted in accordance with the principle set forth under the 1989 Declaration of Helsinki.

Cell isolation and RNA extraction

10 ml EDTA-anticoagulated blood were collected from each subject. PBMCs were prepared by Ficoll-Hypaque density gradient centrifugation. Total RNA was extracted from PBMCs by using Trizol Reagent (Invitrogen, Carlsbad, USA) according to the manufacture’s protocol. The purity and concentration of RNA was evaluated with the NanoDrop ND-1000 (Thermo Fisher Scientific, Wilmington, DE).

High-throughput RNA sequencing and analysis

Briefly, total RNA was extracted using RNeasy mini kit (Qiagen, Germany) according to the manufacturer’s instructions, quantified using NanoDrop ND-2000, and checked for RNA integrity by an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). RNA libraries were generated using TruSeq Stranded RNA Sample Prep Kit (Illumina) following manufacturer’s protocol, quantified with Qubit 2.0 Fluorometer (Life Technologies, USA), validated by Agilent 2100 bioanalyzer (Agilent Technologies, USA) to confirm the insert size and calculate the mole concentration, and subjected to sequencing using Illumina HiSeq 2500 (Illumina, USA). Acquired data were processed to obtain raw reads, which were then extracted for genome mapping with TopHat (version 2.0.9). After genome mapping, Cufflinks v2.1.1 was run with a reference annotation to generate FPKM values for known gene models. DE genes were identified using Cuffdiff. The P value significance threshold in multiple tests was set by the false discovery rate (FDR). The fold changes (FC) were also estimated according to the FPKM in each sample. The DE genes were selected using the following filter criteria: FDR ≤0.01 and FC ≥2. Genes with differential expression levels were identified using edgeR, with GO enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The library construction and sequencing were performed at Shanghai Biotechnology Corporation.

Screening of DE lncRNAs for validation

We screened DE lncRNAs for further validation using the following approaches to validate the results of sequencing experiments using the same cohort together with another independent cohort and investigate the correlations between gene expression levels and clinical characteristics.

First, lncRNAs were evaluated based on the data revealed by sequencing experiments using the following criteria: (I) the FC of DE lncRNAs (P<0.01) in NPSLE group must be greater than two-fold compared to Non-NPSLE group; (II) lncRNA levels may have no significant difference between Non-NPSLE group and HC group; or if the difference is significant, the change tendencies of abnormally expressed lncRNAs between these two group should be same as that between NPSLE-Non-NPSLE-HC group; (III) the contexts of lncRNA-mRNA network are associated with immune system; and (IV) genes with lncRNA-mRNA repeated sequences and without information in databases were excluded.

Second, the degree of DE genes between NPSLE patients and active SLE patients was compared based on the co-expression network results. The top 10 genes with high different degrees and the qualified requirements for signal values were selected.

Quantitative real-time PCR (qRT-PCR) validation

Extracted total RNA of 87 participants was reverse-transcribed to complementary DNA (cDNA) with random primer using a RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific). The qRT-PCR reaction was performed in duplicate with a total volume of 10 µL, including 5 µL Maxima SYBR Green/ROX qPCR master mix (Thermo Scientific), 0.3 µL PCR Forward Primer (10 µM), 0.3 µL PCR Reverse Primer (10 µM), 2 µL cDNA, and 2.4 µL double-distilled water. ViiA 7 Real-Time PCR System (Applied Biosystems) was used and the samples were subjected to the following thermal cycling: 95 °C for 10 min; 40 cycles of 95 °C for 15 s and annealing/extension for 1 min at 60 °C, followed by a melting curve from 65 °C to 95 °C in 0.5 °C increments. The qRT-PCR primer sequences were displayed in Table S2. Beta-actin was used as internal control and the relative gene expression was calculated using 2^−ΔCt comparative threshold cycle method.

Bioinformatic analysis

Pathway analysis and GO analysis were applied to explore the potential roles that the DE mRNAs played in biological pathways or GO terms, including the following three categories: biological process, cellular component, and molecular function. Fisher’s exact test was used to find if there was more overlap between the DE list and the GO annotation list than would be expected by chance. The P value denoted the significance of GO terms enrichment in the DE genes. The KEGG (http://www.genome.jp/kegg) database was used to determine the pathways associated with DE gene. The enrichment and statistics calculation of KEGG pathway were similar to the GO analysis. A co-expression analysis was performed by associating the expression profiles of DE-lncRNAs with DE-mRNAs. This association was based on the information of the lncRNA-associated coding regions, which was supplied by the sequencing analysis.

Construction of coding-non-coding gene (CNC) coexpression network

Taking the role of lncRNAs in modulating transcription into account, we analyzed the relationship between DE lncRNAs and mRNAs in hope of uncovering possible mRNA targets for differential lncRNAs. Given the normalized signal intensity (>6), FC (≥2), FDR (<0.01), and the relevance of mRNA in SLE, 36 DE lncRNAs and 151 DE mRNAs [these mRNAs were broadly selected for several reasons: (I) their abnormal expressions in SLE are validated in other studies; (II) their corresponding genes are susceptible genes in SLE; (III) their novel functions in immune response are established by others.] were selected for construction of CNC gene coexpression network. The network construction was based on the correlation analysis between these DE lncRNAs and mRNAs. Pearson correlation coefficients over 0.90 were selected for the construction with Cytoscape software (The Cytoscape Consortium, San Diego, USA) and shown in Supplementary material.

Olfactory function assessment by computerized testing

The odor threshold test measures a person’s ability to detect an odor. Thirteen concentrations of an industry standard odorant are used. The participant smells two samples, one blank and one with an odorant. They are then asked to determine which of the two samples is the stronger. The concentrations are decreased or increased depending on whether the stimulus is identified correctly.

The odor memory test assesses a person’s ability to remember odors that were previously presented. The test begins with the presentation of 10 odors followed by a 10-minute break. Next, 20 odors containing the 10 original and 10 new odors are presented. After each odor presentation, the participant is asked to identify the odor and to determine if the odor is one of the original odors or a new odor. This provides measures of both semantic and episodic odor memory.

All the tests were done according to the manual of OLFACT^TM as described in the website: http://www.osmicenterprises.com/index.html.

Statistical analysis

All data represent mean ± SEM. Differences between two groups were determined by unpaired Student’s t-test if the variance was normally distributed. Comparisons among three or more groups were conducted using one-way ANOVA. Randomization and blinding strategy were used whenever possible. Data were analyzed and visualized with GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA) and a two-tailed P value <0.05 was considered statistically significant.

Results

LncRNAs and mRNAs expression profiles in NPSLE patients

To determine whether the expression levels of lncRNAs were associated with SLE pathogenesis, we firstly performed whole-transcriptome sequencing (RNA-Seq) on the peripheral blood mononuclear cells (PBMCs) of 11 SLE patients and 6 healthy controls. A total of 1,515 lncRNAs and 2,098 mRNAs were DE between these two groups using scatter plot and hierarchical clustering analyses (Figure 1, http://cdn.amegroups.cn/static/other/b191243af9a0d8a7f8fa5ee14a9a4620/atm.2020.03.25-1.xls, http://cdn.amegroups.cn/static/other/2284aef2c7af028a615bd36c215970a9/atm.2020.03.25-2.xlsx; FC ≥2, P<0.01). To reveal DE lncRNA profiles in NPSLE patients, these 11 SLE patients were classified into NPSLE and Non-NPSLE group, and DE lncRNAs in Figure 1 were further analyzed between the two groups. As shown in Figure 2 and http://cdn.amegroups.cn/static/other/ea6fab857009712577755c3a449f940a/atm.2020.03.25-3.xlsx, http://cdn.amegroups.cn/static/other/ae9e70ffaa1ad61a20e131857d1e8749/atm.2020.03.25-4.xlsx, 451 lncRNAs and 272 mRNAs were DE between NPSLE and Non-NPSLE group, among which, 29 lncRNAs and 43 mRNAs were upregulated (http://cdn.amegroups.cn/static/other/a690c7fc04662153fadadc94c29ae7f1/atm.2020.03.25-5.xls, http://cdn.amegroups.cn/static/other/03bf6ea47a10ad0b983eaadb34f49b0b/atm.2020.03.25-6.xls), and 19 lncRNAs and 5 mRNAs were downregulated (http://cdn.amegroups.cn/static/other/a690c7fc04662153fadadc94c29ae7f1/atm.2020.03.25-5.xls, http://cdn.amegroups.cn/static/other/03bf6ea47a10ad0b983eaadb34f49b0b/atm.2020.03.25-6.xls) both in Figure 1 (http://cdn.amegroups.cn/static/other/b191243af9a0d8a7f8fa5ee14a9a4620/atm.2020.03.25-1.xls, http://cdn.amegroups.cn/static/other/2284aef2c7af028a615bd36c215970a9/atm.2020.03.25-2.xlsx) and Figure 2 (http://cdn.amegroups.cn/static/other/ea6fab857009712577755c3a449f940a/atm.2020.03.25-3.xlsx, http://cdn.amegroups.cn/static/other/ae9e70ffaa1ad61a20e131857d1e8749/atm.2020.03.25-4.xlsx), suggesting that these lncRNAs and mRNAs might play roles in NPSLE pathogenesis.

Figure 1 Expression profiles of lncRNAs and mRNAs in SLE patients and HCs. (A,C) Heat map showing hierarchical clustering of lncRNAs and mRNAs with expression changes greater than two-fold and P value <0.01. SLE patient group: SLE1-11; healthy control group: HC1-6. Red and green colors represent up- and downregulated genes, respectively; (B,D) Volcano plot of differentially expressed lncRNAs and mRNAs. The horizontal white line represents a P value of 0.01, and vertical white lines represent 2.0-fold changes up and down. X axes are the fold change values (log₂ scaled), and Y axes are the P values (log₂ scaled). Red and blue plots represent up- and downregulated genes, respectively.

Figure 2 Expression profiles of lncRNAs and mRNAs in NPSLE patients and Non-NPSLE patients. (A,C) Hierarchical clustering of lncRNAs and mRNAs with expression changes greater than two-fold and P value <0.01. NPSLE patient group: NPSLE1-6; Non-NPSLE patient group: Non-NPSLE1-5. Red and green colors represent up- and downregulated genes, respectively; (B,D) Volcano plot of differentially expressed lncRNAs and mRNAs. The horizontal white line represents a P value of 0.01, and vertical white lines represent 2.0-fold changes up and down. X axes are the fold change values (log2 scaled), and Y axes are the P values (log2 scaled). Red and blue plots represent up- and downregulated genes, respectively.

The top 20 DE lncRNAs/mRNAs between NPSLE and Non-NPSLE groups are listed in Table 1 and Table 2, respectively. The RNA-seq data have been submitted to NCBI Gene Expression Omnibus and are accessible with the GEO Series accession number GSE90524 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE90524).

Table 1 Top 20 DE lncRNAs between NPSLE and Non-NPSLE patients.
Full table

Table 2 Top 20 DE mRNAs between NPSLE and Non-NPSLE patients
Full table

Verification of lncRNA expression profiles and their correlations with clinical characteristics

To independently validate our data, DE lncRNAs between NPSLE and Non-NPSLE group were explored using more stringent criteria (as shown in Methods of “Screening of DE lncRNAs for validation”). In this way, 10 aberrantly expressed lncRNAs were selected for qRT-PCR validation in PBMCs from 26 NPSLE patients, 31 Non-NPSLE patients and 30 healthy subjects (Figure 3).

Figure 3 Relative expression levels of ten selected lncRNAs in PBMCs of NPSLE patients, Non-NPSLE patients and HCs identified using qRT-PCR. (A,D,E,F,G,I,J) the expression levels of seven lncRNAs (NONHSAT208182.1, NONHSAT182114.1, NONHSAT106801.2, NONHSAT039491.2, ENST00000356215, NONHSAT087499.2 and NONHSAT207026.1) were significantly upregulated NPSLE patients compared to Non-NPSLE patients; (B,C) the expression levels of two lncRNAs (NONHSAT001281.2 and NONHSAT024353.2) were significantly downregulated in NPSLE patients compared to Non-NPSLE patients; (H) the difference in lncRNA ENST00000414046 was not significant. n(NPSLE)= 26, n(Non-NPSLE) =31, n (HC) =30. *, P<0.05; **, P<0.01; ***, P<0.001.

To determine the specific roles of lncRNAs in NPSLE, correlations of these DE lncRNA expression levels with experimental characteristics of NPSLE patients were analyzed. As shown in Figure 4A,B,C,D,E,F, the expression levels of NONHSAT039491.2, NONHSAT106801.2, NONHSAT024353.2 and NONHSAT207026.1 were negatively associated with SLEDAI scores, frequency of CD19⁺ B cells, serum IgG levels, and glomerular filtration rate (GFR), respectively; whereas NONHSAT207026.1 and NONHSAT208182.1 levels positively correlated with serum uric acid (UA) and urinary red blood cells (RBC) counts, respectively; no significant correlation existed between these DE lncRNAs and characteristics including age, disease duration, serum immunoglobulin (Ig) E, IgA, complement 3, complement 4, C-reactive protein, 24-hour proteinuria, albumin and IgG in cerebrospinal fluid (CSF, the detailed correlation analysis results was shown in Table S3). Further, the role of these DE genes in visceral involvement were analyzed, and found the expression levels of NONHSAT208182.1, NONHSAT024353.2, NONHSAT106801.2 (r=−0.05, P=0.049) and ENST00000356215 were significantly association with pulmonary arterial pressure (PAH) levels of these NPSLE patients (Figure 4G,H,I). Olfactory dysfunction has been confirmed in various central nervous system (CNS) diseases, as well as in NPSLE (18). To explore the roles of lncRNAs in NPSLE, the relationship of lncRNAs with olfactory functions of NPSLE patients was studied, and found negative association of ENST00000356215 with olfactory threshold (Figure 4J), NONHSAT039491.2 and NONHSAT001281.2 with olfactory memory (Figure 4K,L), respectively.

Figure 4 Correlation between DE lncRNAs and experimental characteristics in NPSLE patients. (A,B,C,D) The expression levels of NONHSAT039491.2, NONHSAT106801.2, NONHSAT024353.2 and NONHSAT207026.1 were negatively associated with systemic lupus erythematosus disease activity index (SLEDAI) scores, frequency of CD19+ B cells, serum IgG levels, and glomerular filtration rate (GFR), respectively. (E,F) NONHSAT207026.1 and NONHSAT208182.1 levels were positively correlated with serum uric acid and urinary red blood cells counts, respectively; (G,H,I) NONHSAT208182.1, NONHSAT024353.2 and ENST00000356215 showed significant association with pulmonary arterial pressure levels of NPSLE patients; (J,K,L) the negative association of ENST00000356215 with olfactory threshold, NONHSAT039491.2 and NONHSAT001281.2 with olfactory memory.

Given the multitude of clinical manifestations related to NPSLE, the correlations of expression levels of DE lncRNAs with clinical manifestations were analyzed. Compared with patients without such symptoms, expression levels of NONHSAT106801.2 and NONHSAT208182.1, NONHSAT207026.1, and NONHSAT039491.2 were significantly downregulated in those with fever (Figure 5A,B), arthralgia (Figure 5C), dizziness and headache (Figure 5D), respectively; while expression levels of ENST00000356215, NONHSAT024353.2 and NONHSAT208182.1 were upregulated in patients with oral ulcer (Figure 5E), consciousness disorder (Figure 5F) and unstable walking (Figure 5G), respectively.

Figure 5 Correlation between DE lncRNAs and clinical characteristics in NPSLE patients. (A,B) NONHSAT106801.2 and NONHSAT208182.1 expression levels in NPSLE patients with and without fever; (C) NONHSAT207026.1 expression level in NPSLE patients with and without arthralgia; (D) NONHSAT039491.2 expression level in NPSLE patients with and without dizziness and headache symptom; (E) ENST00000356215 expression level in NPSLE patients with and without oral ulcer; (F) NONHSAT024353.2 expression level in NPSLE patients with and without consciousness disorder; (G) NONHSAT208182.1 expression level in NPSLE patients with and without unstable walking symptom; (H,I) NONHSAT106801.2 and NONHSAT024353.2 expression levels in anti-SSA+ and anti-SSA- NPSLE patients; (J,K) NONHSAT087499.2 and NONHSAT039491.2 expression levels in anti-RNP+ and anti-RNP- NPSLE patients; (L) NONHSAT039491.2 expression level in anti-dsDNA+ and anti-dsDNA- NPSLE patients. *, P<0.05; **, P<0.01; ***, P<0.001.

Available studies on the significance of different autoantibodies in NPSLE have shown controversial results (19). In our study, NONHSAT106801.2 was significantly upregulated while NONHSAT024353.2 was downregulated in anti-SSA⁺ patients compared to anti-SSA^- patients (Figure 5H,I); NONHSAT087499.2 was significantly upregulated in anti-ribonucleoprotein antibody (anti-RNP)⁺ patients compared to anti-RNP^- patients while NONHSAT039491.2 was down-regulated both in anti-RNP⁺ and anti-dsDNA⁺ patients compared to anti-RNP^- and anti-dsDNA^- patients, respectively (Figure 5J,K,L); no difference of lncRNA levels was found between groups positive and negative for anti-cardiolipin antibody, anti-nuclear antibody, anti-smith antibody, anti-histone antibody, anti-ribosomal P antibody and anti-nucleosome antibody, respectively.

Considering that clinical treatment including corticosteroids and immunosuppressants may affect the expression of lncRNAs, we analyzed the relationship between the expression level of lncRNA and the therapeutic dosage of glucocorticoid (GLU) and immunosuppressive agents including hydroxychloroquine (HCQ), cyclophosphamide (CYC), and mycophenolate (MMF). It was found that only NONHSAT024353.2 was highly expressed in CTX-treated patients compared with patients without CTX treatment (0.232±0.07019, n=9 vs. 0.5652±0.06684, n=17, P=0.0042), but no relationship was found in the expression of other lncRNAs between immunosuppressive and non-immunosuppressive group, as well as different doses of corticosteroids treatment group (Table S4), these above findings indicating that these abnormally expressed lncRNAs between NPSLE and Non-NPSLE group might not be influenced by clinical therapy.

GO and KEGG pathway analyses of DE mRNAs

Next, to further examine the potential mechanism of NPSLE, GO enrichment analysis was performed to investigate the over-representation of biological processes, cellular components, and specific molecular function associating protein-coding mRNAs. As shown in Figure 6, the top 30 GO from the DE mRNAs were listed. Our data demonstrated that the most significantly enriched cellular components of DE mRNAs in PBMCs of NPSLE were the external side of the axon, intrinsic and integral component of plasma membrane (Figure 6A, http://cdn.amegroups.cn/static/other/5984c3b0471eb93f988c76e7878328c7/atm.2020.03.25-7.xls). The most significantly enriched molecular functions of DE mRNAs were carbon hydrate binding, protein homodimerization activity and substrate-specific transmembrane transporter activity (Figure 6A, http://cdn.amegroups.cn/static/other/5984c3b0471eb93f988c76e7878328c7/atm.2020.03.25-7.xls); The most significantly enriched biological processes of DE mRNAs were leukocyte chemotaxis and migration, and cell chemotaxis (Figure 6A, http://cdn.amegroups.cn/static/other/5984c3b0471eb93f988c76e7878328c7/atm.2020.03.25-7.xls); these together reveal the roles of significantly dysregulated mRNAs in NPSLE.

Figure 6 Enrichment analysis of GO terms and pathways for DE mRNAs between NPSLE and Non-NPSLE group. (A) The top 30 GO analyses which consisted of significant molecular function, cellular component and biological process of DE mRNAs; (B) significantly enriched pathways of DE mRNAs. P <0.01, FC >2.0 or <−2.0.

KEGG pathway enrichment analysis for significantly dysregulated mRNAs is useful to reveal related pathways and molecular interactions among genes. The results revealed that the DE mRNAs in PBMCs of NPSLE were significantly involved in hematopoietic cell lineage, staphylococcus aureus infection, and osteoclast differentiation (Figure 6B, http://cdn.amegroups.cn/static/other/a916f2ba1dfe9c79658b5d5b12c5539b/atm.2020.03.25-8.xls). Our data showed that the DE mRNAs were each individually associated with 11 pathways, among them, the metabolic pathway, neurotrophin signaling pathway, and chemokine signaling pathway were three of the most enriched pathways within the set of DE mRNAs (Figure 6B, http://cdn.amegroups.cn/static/other/a916f2ba1dfe9c79658b5d5b12c5539b/atm.2020.03.25-8.xls), suggesting that these pathways may be involved in the development of NPSLE.

lncRNA-mRNA co-expression network

Co-expression network analysis was performed between the 193 DE lncRNAs and the 49 DE mRNAs based on the Max-relevance Min-redundancy (mRMR) method. In total, 170 lncRNAs and 46 mRNAs were included in the co-expression network (Figure 7, http://cdn.amegroups.cn/static/other/a7791bdae56436daf1f955a9b6288064/atm.2020.03.25-9.xls). In addition, our data showed that the co-expression network was composed of 216 network nodes and 251 connections. To further address DE mRNAs as possible predicted target of these DE lncRNAs, the expression levels of top 10 DE mRNAs that have characterized proteins were performed using RT-PCR. Compared with that of Non-NPSLE patients, the relative expression levels of CFH, P2RY14, VNN1, IDO1, BCL2A1and KCNJ15 were significantly increased in PBMCs from NPSLE patients (Figure S1). Considering of the co-expression network of DE lncRNAs with DE mRNAs (http://cdn.amegroups.cn/static/other/a7791bdae56436daf1f955a9b6288064/atm.2020.03.25-9.xls), and the negative correlated role of IDO1 and BCL2A1 with NONHSAT001281.2 (Table S3), our data support that differential lncRNAs might operate via modulating expressions of their correlated, relevant mRNAs in NPSLE.

Figure 7 Network graph of DE lncRNA-related mRNAs between NPSLE and Non-NPSLE group. Octagon and circles represent DE lncRNAs and their related DE genes, respectively. Red and green represent upregulated and downregulated lncRNAs or genes, respectively. P<0.01, FC >2.0 or <−2.0.

Discussion

Recent discoveries have revealed that lncRNAs are aberrantly expressed in a broad spectrum of autoimmune diseases (AD) including SLE (14,15,20), and thus, play key roles in promoting and maintaining AD initiation and progression, demonstrating their clinical potential as biomarkers and therapeutic targets. Studies have shown that serum lnc-gas5, lnc0597, lnc-DC could be used as diagnostic markers for SLE (12). LncRNA NEAT1 acted as an inflammatory regulator to activate the MAPK pathway involved in SLE pathogenesis (21). LncRNA MALAT1 could affect the expression level of IL-21 in monocytes by regulating SIRT1 signaling and participate in SLE development (13). lncRNA-Gas5 binds to the DNA binding region of the glucocorticoid receptor, preventing the receptor from interacting with the glucocorticoid response element of genomic DNA, thereby affecting the patient's response to glucocorticoids (22). These above findings supported the pathogenic role of lncRNAs in SLE pathogenesis (23), but the roles of lncRNAs in NPSLE remain largely unknown.

Growing evidence showed that blood–brain barrier (BBB) dysfunction may be essential to the development of NPSLE (24,25), allowing the passive diffusion of auto-reactive antibodies and cytokines in CSF of SLE patients. Considering of the specific link between PBMCs, serum and CSF (26), CSF are rare, difficult to obtain and traumatic, while peripheral blood is easy to obtain, with more clinical utility. In our study, compared with that of Non-NPSLE patients, hundreds of DE lncRNAs were found by RNA-seq and 10 of them were validated by qRT-PCR in PBMCs from NPSLE patients (Figure 3, Figure S1), and the correlation of expression levels of DE lncRNAs with experimental indicators in CSF from NPSLE patients were explored. Unfortunately, there was no significant correlation of expression levels of lncRNA with protein or IgG levels in CSF of NPSLE patients (Table S3). Since protein or IgG levels in CSF positively associated with SLEDAI of NPSLE patients, while in our study, none of these lncRNAs displayed significantly association with SLEDAI except NONHSAT039491.2, thus the link between these lncRNAs with NPSLE activity still needs further investigation.

Multitude of clinical manifestations were related to NPSLE, each syndrome may present various different patterns, which increases the complexity of clinical outcomes (27). In our study, the expression levels of NONHSAT024353.2, NONHSAT039491.2, NONHSAT208182.1 were associated with clinical index of consciousness, dizziness and headache, as well as unstable walking, respectively (Figure 5D,F,G). To predict their possible target gene, except for lncRNA-mRNA co-expression network analysis (Figure 7), the correlations of DE lncRNAs with DE mRNAs confirmed by RT-PCR were analyzed, and found negative correlation of expression levels IDO1 and BCL2A1 with that of NONHSAT208182.1 in PBMCs from NPSLE patients (Table S3), together with the positive association of NONHSAT208182.1 levels with unstable walking (Figure 5G), it is likely that NONHSAT208182.1 may be involved in the development of unstable walking related-NPSLE through targeting IDO1 or BCL2A1, further RNA-pulldown and CHIP experiments will be performed to support the possible target of NONHSAT001281.2 in our next experiments.

Olfactory disturbance was one of the manifestations of NPSLE (28). Recently, the correlation between decreased olfactory abilities and disease activity, especially CNS manifestations of NPSLE patients is supported by increasing studies, indicating that smell decrement might be a useful and easy tool for the physician in early diagnosis of CNS involvement in autoimmune diseases (18). Interestingly, negative association of ENST00000356215 with olfactory threshold, NONHSAT039491.2 and NONHSAT001281.2 with olfactory memory were found in our NPSLE cohort. Since accumulating studies found that decreased sense of smell was associated with history of neuropsychiatric manifestations and higher disease activity, ENST00000356215, NONHSAT039491.2 and NONHSAT001281.2 may be used as predictors of diffuse disorders in NPSLE patients but still need further investigation.

As one of the major causes of death in SLE, SLE related PAH could also be postulated as vasculitis led by immune-mediated vasculopathy (29), a potential link between PAH and NPSLE manifestations, however, have never been studied. In our study, NONHSAT208182.1, NONHSAT024353.2, NONHSAT106801.2 and ENST00000356215 all showed significant association with PAH levels of these NPSLE patients. Furthermore, NONHSAT106801.2 and NONHSAT024353.2 levels were also associated with the presence of anti-SSA, which had been reported to be an independent variable with PAH in SLE (30), supporting that some specific lncRNAs together with autoantibodies could be recommended in early diagnosis and precise management of NPSLE-associated PAH patients.

Auto-antibody mediated vasculopathy is one of the mechanisms for NPSLE (19). To date, more than 116 auto-antibodies have been reported in SLE, and at least 20 of them, including 11 brain-specific and 9 systemic auto-antibodies, have been controversially associated with NPSLE. Among these auto-antibodies, ACA and anti-RNP were significantly associated with specific manifestations of neuropsychiatric disease attributed to SLE (19). In our study, our data support important correlations of NONHSAT106801.2, NONHSAT024353.2, NONHSAT039491.2 and NONHSAT087499.2 in NPSLE patients with some specific autoantibodies including anti-SSA, anti-RNP, and anti-dsDNA, respectively (Figure 5H,I,J,K,L). Although the detection of auto-antibodies in serum or cerebrospinal fluid has not been included as a diagnostic marker within the ACR classification (19,27), our data indicate that lncRNAs together with their specific autoantibodies might be used for precise management of NPSLE. Various inflammatory cytokines and chemokines, including IL-6, IL-1, IL-8, IL-10, tumor necrosis factor-α (TNF-α), interferon-α (IFN-α), IP-10/CXCL10, MCP-1/CCL2, normal T-cell expressed and secreted (RANTES)/CCL5 as well as CX3CR1, have been reported to increase the permeability of the BBB, providing access for circulating autoantibodies and leukocytes into the CNS, and thus, involved in the pathogenesis of NPSLE (31-33). It is worth noting that, several GO terms were enriched in gene lists through a bioinformatics analysis. Since lncRNAs have been reported to play important roles in immune cell homeostasis (34), and chemokines are crucial for NPSLE development, our results suggest important effects of these lncRNAs together with chemokines in the diagnosis of NPSLE.

Plenty of studies have shown that lncRNAs are expressed in a highly lineage-specific manner. Cao et al. elucidated for the first time that DCs-specific lncRNA lnc-DC can inhibit SHP1-mediated de-phosphorylation by directly binding to the structural region near the STAT3 phosphorylation site Tyr705, thereby enhancing STAT3 signaling and participating in DCs differentiation and functional maintenance (35). Chen et al. found that PU.1 regulated lnc-MC can regulate mononuclear macrophage differentiation by "sponge" adsorption of miRNA-199a-5p (36). In our study, compared with that of Non-NPSLE patients, NONHSAT208182.1 and NONHSAT087499.2 were specifically highly expressed while NONHSAT001281.2 were reduced in CD19⁺ B cells from NPSLE patients, NONHSAT182.114.1 and ENST00000356215 were specifically highly expressed in CD4⁺ T cells, while NONHSAT024353.2 were significantly reduced in CD19⁺ B cells and NK cells (data not shown), supporting that DE lncRNAs may play a role in the development of NPSLE by affecting different subsets of PBMCs from NPSLE patients.

In summary, to the best of our knowledge, this is the first study that describes the expression profiles of human lncRNAs in NPSLE using RNA-seq. Altered lncRNAs may play potential roles in the pathogenesis and/or development of “the most clinical challenging ‘visceral’ involvement of SLE”. To accurately and comprehensively elucidate the role of lncRNAs in NPSLE, more comprehensive laboratory and clinical researches will be needed to confirm whether these lncRNAs play essential roles in the pathogenesis, treatment, and prognosis of NPSLE.

Conclusions

We have demonstrated a comprehensive lncRNA expression profile in PBMCs of NPSLE patients, whose deregulation could be related to the pathophysiology of the disease. More specifically, this study has highlighted some pathways that could be crucial in the onset or the maintenance of the disease. Our finding shed light on the understanding of the molecular mechanisms of lncRNAs in the pathogenesis of NPSLE.

Figure S1 The expression levels of DE mRNAs in PBMCs from NPSLE, Non-NPSLE, and healthy subjects. All data are means ± SEM. n (NPSLE) =26, n (Non-NPSLE) =31, n (HC) =30. *, P<0.05; **, P<0.01; ***, P<0.001.

Table S1 Demographic and clinical characteristics of SLE patients
Full table

Table S2 Primers sequences for qRT-PCR
Full table

Table S3 Correlation between expression levels of 9 abnormally expressed lncRNAs and clinical characteristics of NPSLE patients
Full table

Table S4 Correlation between expression levels of 9 abnormally expressed lncRNAs and clinical treatments of NPSLE patients
Full table

Acknowledgments

The authors thank the patients who were involved in this study and the individuals from the rheumatology departments who took part in the care of these patients.

Funding: The work was supported by the Major International (Regional) Joint Research Project (No. 81720108020), National Natural Science Foundation of China (No. 81871283, 81571586, 81501347 and 81601418), Research Project of Jiangsu Province Health Committee (H2019060), National Natural Science Foundation of Jiangsu (BK20150098), Jiangsu Provincial Special Program of Medical Science (BE2015602), Medical Science and Technology Development Key Project of Nanjing (ZKX15018), and Jiangsu Province Six Talent Peaks Project (No. 2016-WSN-148).

Footnote

Conflicts of Interest: The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was approved by the Ethics Committee of the Affiliated Drum Tower Hospital of Nanjing University Medical School (SC201700201) and was conducted in accordance with the principle set forth under the 1989 Declaration of Helsinki.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Tsokos GC, Lo MS, Costa Reis P, et al. New insights into the immunopathogenesis of systemic lupus erythematosus. Nat Rev Rheumatol 2016;12:716-30. [Crossref] [PubMed]
2. Govoni M, Bortoluzzi A, Padovan M, et al. The diagnosis and clinical management of the neuropsychiatric manifestations of lupus. J Autoimmun 2016;74:41-72. [Crossref] [PubMed]
3. Sharaf-Eldin WE, Kishk NA, Gad YZ, et al. Extracellular miR-145, miR-223 and miR-326 expression signature allow for differential diagnosis of immune-mediated neuroinflammatory diseases. J Neurol Sci 2017;383:188-98. [Crossref] [PubMed]
4. Ransohoff JD, Wei Y, Khavari PA. The functions and unique features of long intergenic non-coding RNA. Nat Rev Mol Cell Biol 2018;19:143-57. [Crossref] [PubMed]
5. Arun G, Diermeier SD, Spector DL. Therapeutic Targeting of Long Non-Coding RNAs in Cancer. Trends Mol Med 2018;24:257-77. [Crossref] [PubMed]
6. Wang DQ, Fu P, Yao C, et al. Long Non-coding RNAs, Novel Culprits, or Bodyguards in Neurodegenerative Diseases. Mol Ther Nucleic Acids 2018;10:269-76. [Crossref] [PubMed]
7. De Rosa S, Arcidiacono B, Chiefari E, et al. Type 2 Diabetes Mellitus and Cardiovascular Disease: Genetic and Epigenetic Links. Front Endocrinol (Lausanne) 2018;9:2. [Crossref] [PubMed]
8. Chen YG, Satpathy AT, Chang HY. Gene regulation in the immune system by long noncoding RNAs. Nat Immunol 2017;18:962-72. [Crossref] [PubMed]
9. Ranzani V, Rossetti G, Panzeri I, et al. The long intergenic noncoding RNA landscape of human lymphocytes highlights the regulation of T cell differentiation by linc-MAF-4. Nat Immunol 2015;16:318-25. [Crossref] [PubMed]
10. Li LJ, Zhao W, Tao SS, et al. Comprehensive long non-coding RNA expression profiling reveals their potential roles in systemic lupus erythematosus. Cell Immunol 2017;319:17-27. [Crossref] [PubMed]
11. Luo Q, Li X, Xu C, et al. Integrative analysis of long non-coding RNAs and messenger RNA expression profiles in systemic lupus erythematosus. Mol Med Rep 2018;17:3489-96. [PubMed]
12. Wu GC, Li J, Leng RX, et al. Identification of long non-coding RNAs GAS5, linc0597 and lnc-DC in plasma as novel biomarkers for systemic lupus erythematosus. Oncotarget 2017;8:23650-63. [PubMed]
13. Yang H, Liang N, Wang M, et al. Long noncoding RNA MALAT-1 is a novel inflammatory regulator in human systemic lupus erythematosus. Oncotarget 2017;8:77400-6. [PubMed]
14. Tang Y, Zhou T, Yu X, et al. The role of long non-coding RNAs in rheumatic diseases. Nat Rev Rheumatol 2017;13:657-69. [Crossref] [PubMed]
15. Wu GC, Pan HF, Leng RX, et al. Emerging role of long noncoding RNAs in autoimmune diseases. Autoimmun Rev 2015;14:798-805. [Crossref] [PubMed]
16. Hochberg MC. Updating the American College of Rheumatology revised criteria for the classification of systemic lupus erythematosus. Arthritis Rheum 1997;40:1725. [Crossref] [PubMed]
17. The American College of Rheumatology nomenclature and case definitions for neuropsychiatric lupus syndromes. Arthritis Rheum 1999;42:599-608. [Crossref] [PubMed]
18. Shoenfeld N, Agmon-Levin N, Flitman-Katzevman I, et al. The sense of smell in systemic lupus erythematosus. Arthritis Rheum 2009;60:1484-7. [Crossref] [PubMed]
19. Sciascia S, Bertolaccini ML, Roccatello D, et al. Autoantibodies involved in neuropsychiatric manifestations associated with systemic lupus erythematosus: a systematic review. J Neurol 2014;261:1706-14. [Crossref] [PubMed]
20. Li J, Wu GC, Zhang TP, et al. Association of long noncoding RNAs expression levels and their gene polymorphisms with systemic lupus erythematosus. Sci Rep 2017;7:15119. [Crossref] [PubMed]
21. Zhang F, Wu L, Qian J, et al. Identification of the long noncoding RNA NEAT1 as a novel inflammatory regulator acting through MAPK pathway in human lupus. J Autoimmun 2016;75:96-104. [Crossref] [PubMed]
22. Mayama T, Marr AK, Kino T. Differential Expression of Glucocorticoid Receptor Noncoding RNA Repressor Gas5 in Autoimmune and Inflammatory Diseases. Horm Metab Res 2016;48:550-7. [Crossref] [PubMed]
23. Long H, Yin H, Wang L, et al. The critical role of epigenetics in systemic lupus erythematosus and autoimmunity. J Autoimmun 2016;74:118-38. [Crossref] [PubMed]
24. Gelb S, Stock AD, Anzi S, et al. Mechanisms of neuropsychiatric lupus: The relative roles of the blood-cerebrospinal fluid barrier versus blood-brain barrier. J Autoimmun 2018;91:34-44. [Crossref] [PubMed]
25. Hirohata S, Sakuma Y, Matsueda Y, et al. Role of serum autoantibodies in blood brain barrier damages in neuropsychiatric systemic lupus erythematosus. Clin Exp Rheumatol 2018;36:1003-7. [PubMed]
26. Yoshio T, Okamoto H, Kurasawa K, et al. IL-6, IL-8, IP-10, MCP-1 and G-CSF are significantly increased in cerebrospinal fluid but not in sera of patients with central neuropsychiatric lupus erythematosus. Lupus 2016;25:997-1003. [Crossref] [PubMed]
27. Faria R, Goncalves J, Dias R. Neuropsychiatric Systemic Lupus Erythematosus Involvement: Towards a Tailored Approach to Our Patients? Rambam Maimonides Med J 2017; [Crossref] [PubMed]
28. Bombini MF, Peres FA, Lapa AT, et al. Olfactory function in systemic lupus erythematosus and systemic sclerosis. A longitudinal study and review of the literature. Autoimmun Rev 2018;17:405-12. [Crossref] [PubMed]
29. Tselios K, Gladman DD, Urowitz MB. Systemic lupus erythematosus and pulmonary arterial hypertension: links, risks, and management strategies. Open Access Rheumatol 2016;9:1-9. [Crossref] [PubMed]
30. Huang C, Li M, Liu Y, et al. Baseline Characteristics and Risk Factors of Pulmonary Arterial Hypertension in Systemic Lupus Erythematosus Patients. Medicine (Baltimore) 2016;95:e2761. [Crossref] [PubMed]
31. Fragoso-Loyo H, Richaud-Patin Y, Orozco-Narvaez A, et al. Interleukin-6 and chemokines in the neuropsychiatric manifestations of systemic lupus erythematosus. Arthritis Rheum 2007;56:1242-50. [Crossref] [PubMed]
32. Guo L, Lu X, Wang Y, et al. Elevated levels of soluble fractalkine and increased expression of CX3CR1 in neuropsychiatric systemic lupus erythematosus. Exp Ther Med 2017;14:3153-8. [Crossref] [PubMed]
33. Pranzatelli MR. Advances in Biomarker-Guided Therapy for Pediatric- and Adult-Onset Neuroinflammatory Disorders: Targeting Chemokines/Cytokines. Front Immunol 2018;9:557. [Crossref] [PubMed]
34. Mowel WK, Kotzin JJ, McCright SJ, et al. Control of Immune Cell Homeostasis and Function by lncRNAs. Trends Immunol 2018;39:55-69. [Crossref] [PubMed]
35. Wang P, Xue Y, Han Y, et al. The STAT3-binding long noncoding RNA lnc-DC controls human dendritic cell differentiation. Science 2014;344:310-3. [Crossref] [PubMed]