RNA Sequencing for Gene Expression Profiles in Peripheral Blood Mononuclear Cells with Ankylosing Spondylitis RNA

Several previous studies have attempted to investigate the regulatory mechanisms underlying gene expression in ankylosing spondylitis (AS). However, the specific molecular pathways underlying this condition remain unclear. Previous research used next-generation RNA sequencing to identify a series of differentially expressed genes (DEGs) in peripheral blood mononuclear cells (PBMCs) when compared between patients with AS and healthy controls, thus implying that these DEGs may be related to AS. Furthermore, by screening these DEGS, it may be possible to facilitate clinical diagnosis and optimize treatment strategies. In order to test this hypothesis, we recruited 15 patients with AS and 15 healthy controls. We randomly selected five subjects from each group of patients for RNA sequencing analysis. Sequence reads were generated by an Illumina HiSeq2500 platform and mapped on to the human reference genome using HISAT2. We successfully identified 973 significant DEGs (p < 0.05) in PBMCs. When compared with controls, 644 of these genes were upregulated (with a fold change (FC) > 2) in AS patients and 329 were downregulated (FC < 0.5). Our analysis identified numerous genes related to immune response. Gene Ontology (GO) analysis indicated that these DEGs were significantly related to the positive regulation of epidermal growth factor-activated receptor activity, the positive regulation of the ERBB (erb-b2 receptor tyrosine kinase) signaling pathway, the differentiation of trophoblast giant cells, oxygen transport, immune-related pathways, and inflammation-related pathways. The DEGs were also closely related to the TNF and NF-κB signaling pathways. Six DEGs were verified by quantitative real-time polymerase chain reaction (qRT-PCR). Receiver operating characteristic (ROC) curve analysis indicated that IL6 may represent a useful biomarker for diagnosing AS. The development of new biomarkers may help us to elucidate the specific mechanisms involved in the development and progression of AS.


Introduction
Ankylosing spondylitis (AS) is an immune-mediated chronic inflammatory form of arthritis and is characterized by chronic nonspecific inflammation and pathological bone formation, the latter representing a common clinical form of spondyloarthritis (SpA) [1]. The incidence of AS in China is approximately 0.3% and predominantly affects adults (mean age: 25 years; range: 15-35 years) [2]. Chronic inflammation of the spinal joints can lead to severe chronic pain and stiffness, ultimately leading to bone stiffness in the spine; this can also exert impact on several other systems [3]. AS is characterized by chronic progressive and refractory characteristics and can cause irreversible damage to the central axis of the spine; this results in the spine fusing with the sacroiliac joint, thus resulting in reduced spinal activity [4]. AS exerts serious effects on a patient's quality of life and is associated with a significant economic burden to society and families. Previous research suggested that this disease is highly correlated with the MHC (major histocompatibility complex) class I gene, HLA-B27 [5]. However, the specific cause of AS remains unclear, and therapeutic options for the treatment of AS remain inadequate. Over the last decade, new biological agents have been developed that have had a profound effect on the success rates of AS treatment. However, approximately 30% of patients fail to tolerate these drugs or experience differing degrees of adverse reactions [6]. Therefore, there is an urgent need to identify new biomarkers that may act as diagnostic or prognostic indicators for AS. The discovery of such biomarkers is likely to prove invaluable in the prevention, treatment, and control of this disease.
Previous researches involving the identification of molecular mechanisms and novel biomarkers associated with cancer, stroke, and diabetes have involved mRNA expression profiling performed by microarray analysis or high-throughput RNA sequencing [7][8][9][10]. Other studies have described the application of high-throughput technology for autoimmune diseases, such as systemic lupus erythematosus, autoimmune thyroid, and rheumatoid arthritis [11][12][13][14]. High-throughput methodology has already been used to study AS [15]; however, this previous study only focused on synovial tissue [15]. In the present study, we used RNA sequencing to construct a protein-coding gene regulation network in peripheral blood mononuclear cells (PBMCs) isolated from AS patients and healthy controls.

Patients and Controls.
Fifteen patients with AS were recruited from the Department of Rheumatology of the First Affiliated Hospital of Anhui University of Chinese Medicine. These patients were diagnosed by visiting staff according to the American College of Rheumatology (ACR) modified New York criteria [16,17]. In addition, we also recruited 15 age-and sex-matched healthy subjects as controls. None of the patients or controls had any previous history of cardiovascular disease, diabetes, hepatitis, malignancy, or other autoimmune and inflammatory illnesses. The

Functional Enrichment
Analysis. Next, we aimed to gain a better understanding of the functionality of the DEGs identified. To do this, we used the R package (v 3.5.1) [27] and clusterProfiler to perform GO term enrichment analysis [28,29] and KEGG [30] pathway analysis.
2.5. Protein-Protein-Interaction (PPI) Network Construction and Module Analysis. Next, we used the STRING tool [31] to map PPIs for all DEGs with a composite interaction score ≥ 0:4.   Figure 1: Volcano plot (a, blue and red indicate >twofold decreased and increased expression in AS, respectively) and scatter plot (b, blue and red indicate >twofold decreased and increased expression in AS, respectively) of DEGs between the AS and control groups. Gray indicates no significant difference.

BioMed Research International
Cytoscape [32] was used to visualize the PPI network, and modules were filtered by the Molecular Complex Detection (MCOD) plug-in [33] using the following parameters: degree cut − off = 2, k − core = 2, node score cut − off = 0:2, and max depth = 100. Functional enrichment within each module was considered if the MCODE score was ≥4 and the node frequency was ≥10. GO and KEGG enrichment analysis for the DEGs within the four modules was performed in clusterProfiler.
2.6. Validation of DEGs. Quantitative real-time polymerase chain reaction (qRT-PCR) was used to verify the RNA sequencing data using β-actin as the internal control. Relative mRNA expression was calculated by the 2 -ΔΔCT method [34]. In total, the expression levels of six genes were quantified (TNFAIP3, IL1β, IL6, GPR55, CCR2, and CXCL5). The primers used for qRT-PCR are shown in Table 1.

Statistical
Analysis. Data are presented as the mean ± standard error of the mean. Data were analyzed using Stu-dent's t-test. p < 0:05 was considered to indicate a statistically significant difference. Student's t-test and ROC were analyzed by GraphPad Prism (version 8) and presented as the mean ± standard error of the mean.  Tables 2 and 3. 3.2. GO and Pathway Enrichment Analysis. Next, we used clusterProfiler to perform GO enrichment analysis for the  Figure 3 shows the top 30 subclasses. GO analysis revealed that the identified DEGs were predominantly associated with a range of biological processes, including positive regulation of the epidermal growth factor-activated receptor activity, positive regulation of the ERBB signaling pathway, the differentiation of trophoblast giant cells, and oxygen transport. In terms of cellular   (Figure 3(a)).

Results
Pathway analysis showed that the pathways most closely related to the DEGs were the TNF signaling pathway, the PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and the NF-κB signaling pathway (Figure 3(b)).

PPI Network
Analysis. The DEGs identified in our analysis were then used to build a PPI network based on string databases (Figure 4(a)). The PPI network consisted of 768 nodes and 3824 edges. Network analysis showed that hundreds of genes were able to interact with 10 more other genes. In particular, IL6, EGF, and CDH1 were shown to interact with more than 100 genes; node distribution is shown in Figure 4(b). In total, 29 modules were identified by MCODE, operating with default criteria. Table 4 lists these modules in descending order according to MCODE scores > 2. We selected five modules (modules 1, 2, 3, 4, and 5) with an MCODE score ≥ 3 and ≥10 nodes for module network visualization ( Figure 5).

Verification of DEGs by qRT-PCR.
Compared with the control group, qRT-PCR detected significantly higher levels of expression for TNFAIP3, IL1β, and IL6, in the group of AS patients ( Figure 6). AS patients also had lower expression levels of GPR55, CCR2, and CXCL5; this was consistent with the results derived from RNA sequencing, thus indicating that the qRT-PCR verification was reliable.
3.5. Receiver Operating Characteristic (ROC) Curve Analysis of Confirmed mRNAs in PBMCs. ROC curve analysis was used to evaluate the potential diagnostic value of differentially expressed mRNAs that showed statistical significance. Our analysis showed that the levels of TNFAIP3, IL1β, IL6, GPR55, and CXCL5 could discriminate between AS patients and controls. Analysis showed that IL6 had the highest area . Therefore, our analyses suggested that IL6 (p < 0:0001) may be more valuable than the other four mRNAs as a biomarker for AS diagnosis (Figure 7).

Discussion
In the present study, we used high-throughput RNA sequencing to analyze protein-coding mRNA expression profiles in PBMCs isolated from AS patients and controls. We successfully identified 973 DEGs and then performed a range of analyses (GO, KEGG pathway, PPI, and PPI module) in an attempt to identify novel mechanisms for AS.

Conclusions
Collectively, our findings provide clinically useful information relating to the mRNA profile of PBMCs in patients with AS. In addition, bioinformatic methodology was used to predict the potential functional roles of DEGs and explore their possible roles in the pathogenesis of AS. This information is likely to provide us with a better understanding of the pathogenic processes leading to AS. Future research should be aimed at investigating the specific biological functions and molecular mechanisms underlying the roles of these DEGs in the pathogenesis of AS.

Data Availability
The RNA-seq data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.