Identification of lncRNAs associated with the pathogenesis of ankylosing spondylitis

Abstract

Background

Ankylosing spondylitis (AS) is a chronic autoimmune disease affecting the sacroiliac joint. To date, few studies have examined the association between long non-coding RNAs (lncRNAs) and AS pathogenesis. As such, we herein sought to characterize patterns of AS-related lncRNA expression and to evaluate the potential role played by these lncRNAs in this complex autoimmune context.

Methods

We conducted a RNA-seq analysis of peripheral blood mononuclear cell (PBMC) samples isolated from five AS patients and corresponding controls. These data were then leveraged to characterize AS-related lncRNA expression patterns. We further conducted GO and KEGG enrichment analyses of the parental genes encoding these lncRNAs, and we confirmed the validity of our RNA-seq data by assessing the expression of six lncRNAs via qRT-PCR in 15 AS and control patient samples. Pearson correlation analyses were additionally employed to examine the associations between the expression levels of these six lncRNAs and patient clinical index values.

Results

We detected 56,575 total lncRNAs in AS and control patient samples during our initial RNA-seq analysis, of which 200 and 70 were found to be up- and down-regulated (FC > 2 or < 0.05; P < 0.05), respectively, in AS samples relative to controls. In qRT-PCR validation assays, we confirmed the significant upregulation of NONHSAT118801.2, ENST00000444046, and NONHSAT183847.1 and the significant downregulation of NONHSAT205110.1, NONHSAT105444.2, and NONHSAT051856.2 in AS patient samples. We further found the expression of NONHSAT118801.2 and NONHSAT183847.1 to be positively correlated with disease severity.

Conclusion

Overall, our findings highlight several lncRNAs that are specifically expressed in PBMCs of AS patients, indicating that they may play key functions in the pathogenesis of this autoimmune disease. Specifically, we determined that NONHSAT118801.2 and NONHSAT183847.1 may influence the occurrence and development of AS.

Background

Ankylosing spondylitis (AS) is a chronic autoimmune inflammatory disease [1]. While AS is often asymptomatic during the early stages when it is most amenable to treatment, during the later stages of disease irreversible axial spine damage can develop, resulting in loss of mobility and debilitating pain in affected patients [2]. AS development is believed to be driven by numerous environmental and genetic factors [3], including macrophage activation status, infections with particular bacteria, and the expression of specific genes [4,5,6]. The exact etiology of this disease, however, remains incompletely characterized, and further in-depth analyses of the molecular basis of AS are thus warranted.

Long non-coding RNAs (lncRNAs) are RNAs > 200 nucleotides long that largely lack coding potential, and yet are critically important in physiological and pathological contexts [7]. Indeed, cell- and disease-specific lncRNA expression patterns have been characterized to date [8]. From a functional perspective, lncRNAs can control gene and protein expression at the epigenetic, transcriptional, and post-transcriptional levels via controlling genomic imprinting, chromatin remodeling, histone modification, splicing, transcription, and cell cycle progression [9, 10]. There is also evidence that lncRNAs are related to the incidence of autoimmune diseases including systemic lupus erythematosus [11,12,13], multiple sclerosis [14, 15], rheumatoid arthritis, and psoriasis [16,17,18]. The functional relevance of lncRNAs expressed in peripheral blood mononuclear cells (PBMCs) from AS patients, however, remain to be clarified.

Herein, we assessed AS-related lncRNA expression profiles via the RNA-seq analysis of PBMC samples from five AS patients and controls. We then used qRT-PCR to validate these RNA-seq findings. Together, our data offer a new theoretical framework for understanding the pathological basis of AS, highlighting novel disease-related lncRNAs for future research and therapeutic targeting efforts.

Materials and methods

Sample collection

We enrolled 15 total patients from the Department of Rheumatology at the First Affiliated Hospital of Anhui University of Chinese Medicine that had been diagnosed with AS based on the modified American College of Rheumatology (ACR) New York criteria (1984) [19]. We also enrolled 15 age- and sex-matched healthy controls. No enrolled patients had any history of prior diabetes, hepatitis, malignancies, cardiovascular disease, or other autoimmune/inflammatory diseases. All enrolled AS patients completed Bath ankylosing spondylitis diseases active index (BASDAI) [20], Bath ankylosing spondylitis functional index (BASFI) [21], and visual analog scale (VAS) assessments. In addition, we determined key clinical and laboratory parameters for each of these patients including erythrocyte sedimentation rate (ESR), C-reactive protein (CRP) levels, and immunoglobulin G (IgG) levels. We randomly selected five AS patients and five healthy controls for RNA sequencing analysis.

PBMC preparation and total RNA extraction

Lymphoprep (Stemcell, USA) was utilized to separate PBMCs from the blood of AS patients and controls at room temperature via density gradient centrifugation. A MiRNeasy Mini Kit (Qiagen, Germany) was then used to extract total sample RNA, after which RNA quantity was evaluated using a NanoDrop 2000 instrument (Thermo Fisher Scientific, USA), and an Agilent 2100 bioanalyzer (Agilent Technologies, USA) was utilized to establish RNA quality [22].

RNA-sequencing

A Truseq® chain total RNA sample preparation kit (Illumina, USA) was used to prepare RNA libraries based upon provided instructions, after which quantification was conducted using a Qubit 2.0 fluorometer, and quality was assessed via an Agilent 2100 Bioanalyzer. CBOT was then utilized to dilute the library to 10 pm in order to generate clusters, and an Illumina Hiseq 2500 instrument (Illumina) was utilized for sequencing. OG Biotech Inc. (AO Ji biotech, Shanghai, China) conducted all library preparation and sequencing.

Data analysis

FastQC1 (v. 0.11.3) was initially utilized for quality control (QC) assessment of RNA-seq reads, after which rRNA reads, adapter sequences, and low-quality reads were trimmed with the seqtk2 software. BWA-MEM (v2.0.4) was then used to map the trimmed reads to the human reference genome (hg38), after which CIRI was used to predict lncRNAs in these RNA-seq data, and SRPBM was used to quantify lncRNA counts therein. Next, lncRNAs that were differentially expressed between AS and control patients were detected using the EdgeR software based on the following selection criteria: FC ≥ 2 or < 0.5; P < 0.05. The potential functional roles of these lncRNAs were assessed via GO and KEGG enrichment analyses of the parental genes harboring these differentially expressed lncRNAs. In addition, Pearson correlation analyses were used to detect associations between lncRNA and mRNA expression profiles, with a co-expression network being constructed by incorporating interactions with a Pearson’s R > 0.9 and P < 0.05. Cytoscape was then used to visualize this lncRNA-mRNA interaction network.

qRT-PCR

To validate our RNA-seq data, we randomly selected three up-regulated (NONHSAT118801.2, ENST00000444046, and NONHSAT183847.1) and three down-regulated (NONHSAT205110.1, NONHSAT105444.2, and NONHSAT051856.2) lncRNAs from 270 differentially expressed lncRNAs via qRT-PCR in 15 AS and control samples with a Universal SYBR Green Master mix. OG-Biotech Inc. synthesized all primers used for this analysis (Table 1). The 2 − ΔΔCt method was used to assess relative lncRNA expression levels.

Table 1 Primers used in the present study

Statistical analysis

GraphPad Prism 6.0 was used for the assessment of relative lncRNA expression and figure construction. Data are given as means ± standard deviation (SD), and were compared via Student’s t-tests. Relationships between different parameters were compared via Spearman and Pearson correlation analyses. P < 0.05 was the significance threshold for this study.

Results

Patient characteristics

We detected no significant differences in age or sex when comparing the AS and control patient cohorts in the present study (Table 2).

Table 2 Study population clinical characteristics

Identification of AS-related differentially expressed lncRNAs

We next identified lncRNAs that were differentially expressed between AS and control patient PBMCs using EdgeR (FC ≥ 2 or ≤ 0.5; P < 0.05). In total, we identified 270 lncRNAs that were differentially expressed between these two patient groups, of which 200 and 70 were up- and down-regulated in AS patient samples, respectively (Fig. 1a-c). The top 10 up- and down-regulated lncRNAs identified through these analyses are compiled in Tables 3 and 4.

Fig. 1

Scatter plot (a), Volcano plot (b), and Hierarchical clustering (c) analyses highlighting patterns of differential lncRNA expression in PBMC samples from AS patients and controls. Green and blue are used to indicate a > 2-fold reduction in the expression of the indicated lncRNAs in the AS group relative to controls in respective hierarchical clustering analyses and volcano plots (P < 0.05). Red corresponds to a > 2-fold increase in the expression of the indicated lncRNA in the AS group relative to the control group (P < 0.05)

Table 3 The top 10 known differentially expressed lncRNAs that were upregulated in AS patient PMBCs relative to control samples

Table 4 The top 10 known differentially expressed lncRNAs that were downregulated in AS patient PMBCs relative to control samples

Pathway enrichment analyses of differentially expressed lncRNAs

In an effort to more fully explore the potential functional roles of these differentially expressed lncRNAs, we next conducted GO and KEGG enrichment analyses of the parental genes harboring these lncRNAs using the ClusterProfler tool.

GO analyses enabled us to identify key biological processes (BPs), molecular functions (MFs), and cellular components (CCs) enriched among these lncRNAs (Fig. 2a). Significantly enriched BPs included the regulation of NF-κB nuclear import into the nucleus, responses to electrical stimuli, negative regulation of I-κB kinase/NF-κB signaling, and positive regulation of transcription factor import into the nucleus, while significantly enriched CC terms included basement membrane, and significantly enriched MFs included transmembrane receptor protein tyrosine kinase activity.

Fig. 2

The results of GO and KEGG pathway enrichment analyses of parental genes harboring AS-related differentially expressed lncRNAs. a The top 30 enriched GO enrichment terms. b The top 30 classes of KEGG pathway enrichment terms

In a subsequent KEGG analysis, 134 pathways were found to be significantly enriched for these lncRNAs, with the top 30 being shown in Fig. 2b. The most significantly enriched of these pathways included the TNF, PI3K-Akt, HIF-1, Rap1, Hippo, ECM-receptor interaction, and arrhythmogenic right ventricular cardiomyopathy (ARVC) signaling pathways.

Preparation of a lncRNA-mRNA co-expression network

We next prepared a lncRNA-mRNA co-expression network using our RNA-seq data that incorporated pairs of lncRNAs and mRNAs that yielded Pearson correlation coefficient values ≥0.9, with Cytoscape being used to visualize this network. In total, 948 differentially expressed mRNA were obtained, including 628 upregulated expressed and 320 downregulated expressed. Co-expression networks generated using control and AS patient samples different substantially with respect to the number of incorporated nodes and connections (Fig. 3).

Fig. 3

AS-related lncRNA-mRNA interaction network. An interaction network was constructed incorporating differentially expressed genes, with up- and down-regulated genes being marked in red and green, respectively

RNA-seq result validation

We next validated our RNA-seq data by assessing the expression of 6 randomly selected differentially expressed lncRNAs in PBMC samples from 15 AS and 15 control patients via qRT-PCR. This analysis enabled us to confirm that NONHSAT118801.2, ENST00000444046, and NONHSAT183847.1 were significantly upregulated in AS patient PMBCs, whereas NONHSAT205110.1, NONHSAT205110.1, and NONHSAT051856.2 were significantly downregulated in these samples (Fig. 4). The qRT-PCR results were consistent with the RNA-seq results.

Fig. 4

qRT-PCR-based validation of the expression of six differentially expressed lncRNAs in control and AS patient PBMCs

Correlations between lncRNA expression patterns and AS patient clinical findings

Lastly, we evaluated the association between the expression of the identified differentially expressed lncRNAs in AS patients and patient disease status. Specifically, we conducted Pearson correlation analyses comparing the expression of the six validated lncRNAs (NONHSAT118801.2, ENST00000444046, NONHSAT183847.1, NONHSAT205110.1, NONHSAT205110.1, and NONHSAT051856.2) and AS clinical features (ESR, hs-CRP, IgG, BASDAI, BASFAI, and VAS). We determined that NONHSAT118801.2 expression levels were positively correlated with ESR, BASDAI, and BASFAI levels in AS patients (Fig. 5a-c), while NONHSAT183847.1 levels were positively associated with ESR, BASDAI, and CRP (Fig. 5d-f).

Fig. 5

Correlations between lncRNA expression profiles and AS clinical disease activity

Discussion

While structurally similar to mRNAs, lncRNAs generally lack the ability to encode peptides. Although they were originally considered to be irrelevant and a form of transcriptional noise, lncRNAs are now recognized to be key regulators of gene and protein expression at the epigenetic, transcriptional, and post-transcriptional levels. These lncRNAs are thus key determinants of metabolic activity, development, evolution, and disease pathophysiology [22].

A previous microarray study identified five lncRNAs, NDRG1-AS6, CSNK1D-AS8, CD46-AS9, SMYD5-AS2, and NR_045553 that are involved in hip joint ligament tissues in patients with AS [23]. In addition, MSTRG.8559 and LINC00987 were identified as two hub differentially expressed lncRNAs in AS patients [24]. In one recent study, peripheral blood lncRNA-AK001085 levels were found to be reduced in AS patients and to be negatively correlated with CRP and ESR levels in these individuals, suggesting a potential protective role for this lncRNA in this pathological context [25]. Although the role of LncRNAs in the pathogenesis of AS has been preliminarily understood, the regulatory mechanism of most LncRNA in AS remains unclear due to its large system, complex mechanism and diverse functions. Our data provide novel insights into lncRNA expression profiles in AS patient PBMCs, as we were able to identify 270 total differentially expressed lncRNAs, including 200 and 70 that were up- and down-regulated in AS patient samples, respectively. We validated the differential expression of three up- and down-regulated lncRNAs via qRT-PCR to confirm the reliability of our RNA-seq data, The qRT-PCR results were consistent with the RNA-seq resultsand selected these lncRNAs for further analysis. Together, our data have the potential to serve as a framework for future studies of the role of these lncRNAs in the context of AS disease processes.

Our GO analysis results indicate potential roles for these differentially expressed lncRNAs as important regulators of NF-κB, which is a key transcription factor and regulator of inflammatory and oxidative stress-related processes in cells. These lncRNAs were also associated with signaling pathways including the TNF, PI3K-Akt, Rap1, Hippo, and HIF-1 pathways, which are primarily involved in inflammatory immune regulation and osteogenic differentiation. In some prior studies, the PI3K-Akt signaling pathway has been suggested to promote fibroblastic ossification and inflammation in the context of AS [26]. Some of the identified pathways, including the Rap1, Hippo, and HIF-1 pathways, have not previously been studied in the context of AS, and thus warrant future investigation.

To confirm our RNA-seq findings, we utilized a qRT-PCR approach to validate the differential expression of three up-regulated lncRNAs (NONHSAT118801.2, ENST00000444046, NONHSAT183847.1) and three down-regulated lncRNAs (NONHSAT205110.1, NONHSAT205110.1, NONHSAT051856.2) in 15 AS patient samples relative to 15 control samples. This approach confirmed that all six of these lncRNAs were differentially expressed in the expected directions (all P < 0.01), thereby reaffirming the reliability of our RNA-seq data. We also assessed correlations between the expression levels of these lncRNAs and AS patient disease indices, leading us to determine that NONHSAT118801.2 expression was positively correlated with ESR, BASDAI, and BASFAI levels in AS (Fig. 5a-c), while NONHSAT183847.1 expression was positively correlated with ESR, BASDAI, and CRP. These data suggest that NONHSAT118801.2 and NONHSAT183847.1 may be key regulators of the pathogenesis of AS, although additional research will be necessary in order to understand the mechanisms whereby these candidate lncRNAs impact disease progression.

Together, our data offer a new mechanistic basis for understanding the pathogenesis of AS. However, we recognized that the present study has some limitations. First, the sample size was small, necessitating the collection of a larger sample from different regions and involving subjects belonging to different races. Second, we did not assess the ability of the markers to effectively differentiate AS from other rheumatic diseases such as Rheumatoid Arthritis (RA), systemic lupus erythematosus (SLE), and multiple sclerosis (MS). These other rheumatic diseases should be considered in similar future studies so as to strengthen the argument supporting that differentially expressed lncRNAs are useful diagnostic biomarkers for AS.

Conclusion

In summary, we utilized an RNA-seq approach to characterize AS-related lncRNA expression profiles in patient PBMCs. The identified differentially expressed lncRNAs may offer new insight into the molecular etiology of this complex and debilitating disease. However, future mechanistic analyses will be required in order to validate our findings and to fully explore the clinical and therapeutic relevance of these lncRNAs in this disease.