Profiling and Bioinformatics Analysis of Differentially Expressed circRNAs in Spinal Ligament Tissues of Patients with Ankylosing Spondylitis

Recent studies have reported that circular RNAs (circRNAs) play a crucial regulatory role in a variety of human diseases. However, the roles of circRNAs in ankylosing spondylitis (AS) remain unclear. In this study, we conducted circRNA expression profiling of the spinal ligament tissues of patients with AS by RNA sequencing (RNA-seq) and analyzed the potential functions of differentially expressed circRNA by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to investigate the potential mechanisms associated with AS. The results showed that a total of 1,172 circRNAs were detected in the spinal ligament tissue samples, of which 123 circRNAs were significantly differentially expressed by a fold change ≥ 1.5 and p value < 0.05. Among these, 57 circRNAs were upregulated, and 66 were downregulated. GO and KEGG analyses demonstrated that the differentially expressed circRNAs were mainly involved in the regulation of biological processes of peptidyl-serine phosphorylation and human immune system that may be related to AS. In addition, the circRNA/miRNA interaction networks were established to predict the potential roles of differentially expressed circRNAs by bioinformatics analysis. Taken together, these results revealed the expression profiles of circRNAs and the potential functions of the differentially expressed circRNAs in the spinal ligament tissue of patients with AS, which may provide new clues for understanding the mechanisms associated with AS, and proceed to identify novel potential molecular targets for the diagnoses and treatment of AS.


Introduction
Ankylosing spondylitis (AS), which most commonly affects the sacroiliac joint and the axial joint of the spine, is an autoimmune disorder with a global incidence of about 2%-5%, and the patients with AS are often young and middle-aged male [1]. The typical clinicopathological features of AS are inflammation and new bone formation in sacroiliac joint, spine, and peripheral joints (especially hip joint), which finally result in ankylosis [2]. A lot of results have indicated that the inflammation in AS initially occurs at the tendonbone interface, leading to bone proliferation [2][3][4]. Progress in the early diagnosis and treatment of AS has been achieved, while the effect of clinical treatment is not as well as people expected. However, many studies suggested that AS has a high heritability [5]; the pathogenesis of AS has been more likely multifactorial and poorly understood to date. Therefore, to elucidate the pathogenesis of AS would be of great value in theory and clinical applications.
Circular RNAs (circRNAs) are a novel class of endogenous noncoding RNAs with a covalently closed circular structure [6][7][8][9]. Unlike the linear RNAs, circRNAs have no free 3′-end polyA tail and 5′-end cap, which prevents them from being digested by nucleic acid exonuclease [10]. Therefore, the closed circular structure of circRNAs makes them incredibly stable, and that these may be potentially utilized as molecular markers. In addition, circRNAs possess tissue specificity and are often highly conserved among many species [11]. Recent studies have indicated that circRNAs have diverse biological functions and play crucial regulatory roles in many human diseases [12,13]. It has been confirmed that circRNAs regulate the expression of miRNA target genes by acting as microRNA (miRNA) sponges [8,9]. Extensive studies have revealed that numerous miRNAs may be associated with AS, such as miR-29a, miR-335-5p, miR-27a, let-7i, miR-146a, miR-29a, and miR-155 [14,15]. Besides, a previous study indicated that the aberrant expression of a variety of lncRNAs has also been observed in the peripheral blood of patients with AS [16]. However, there is no report on cir-cRNAs in patients with AS.
Therefore, in this study, we recruited three AS patients as the experimental group and three patients with lumbar disc herniation as the control group and collected their spinal ligament tissues. We performed RNA sequencing (RNA-seq), and then, the differentially expressed circRNAs in the spinal ligament tissue between the two groups were analyzed by bioinformatics analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were used to predict the biological functions of the key differentially expressed circRNAs. In addition, the interaction between circRNAs and miRNAs and the networks of circRNA-miRNA was, respectively, predicted and constructed by bioinformatics analysis. This study conducted the circRNA expression profiling of spinal ligament tissues of patients with AS and demonstrated the potential functions of differentially expressed circRNAs, which may provide new clues for studying the mechanisms and potential molecular targets for the treatment of AS.

Materials and Methods
2.1. Patients and Specimens. This study was approved by the Ethics Committee of the Affiliated Hospital of Qingdao University and conformed to the Ethical Guidelines of the Declaration of Helsinki. All of the participants signed informed consent forms upon disclosure of the study details. In this study, the experimental and control groups received posterior lumbar or thoracic decompression and fusion surgery at the Department of Spinal Surgery of TThe Affiliated Hospital of Qingdao University from January 2017 to December 2018. Patients with AS were 64, 58, and 76 years old, and all of them were HLA B27 + . According to kyphosis deformity, the patients with AS have exhibited symptoms of bilateral damage of the sacroiliac joint in the computed tomography results and spinal and sacroiliac joint fusion in the X-ray results, which completely meet the revised New York AS standard [17]. None of the patients with AS were treated with nonsteroidal anti-inflammatory drugs or biological agents and suffered any complication. The participants in control group were patients with lumbar disc herniation, aged 58, 69, and 64 years. All of the participants had no other type of autoimmune diseases, and the characteristics of all of the participants are shown in Table 1. The spinal ligament tissue samples, including supraspinous and interspinous ligaments, were cut and removed from six participants. The tissue samples were snap frozen in liquid nitrogen and stored at -80°C until analysis.
2.2. RNA Sequencing (RNA-seq). We did sampling once from each patient in each group; that is, three biological repeats/ RNA-sequencing in each group and RNA sequencing were used for circRNA expression profiling in the spinal ligament tissue samples from three AS patients and three patients with lumbar disc herniation under the same conditions at a time [18,19]. According to the manufacturer's instructions, total RNA was isolated from the two group samples using a Magen Hipure Total RNA Mini Kit (Magen, Guangzhou, China). Subsequently, the concentration and integrity of the isolated RNA were determined using the Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and the Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA), respectively. RNA-seq libraries were prepared as previously [20]. Briefly, rRNA was removed from total RNA using a KAPA RNA HyperPrep Kit with RiboErase (HMR) for Illumina® (Kapa Biosystems, Inc., Woburn, MA, USA). The RNA samples were fragmented and reversely transcribed to the firststrand cDNA, and then, the directional second strand was synthesized. After cDNA synthesis, a tail and adapter were ligated onto the purified cDNA, and then, the cDNA was amplified. Subsequently, the cDNA library quality and concentration were evaluated using the Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA). For sequencing applications, the qPCR-based KAPA Biosystems Library Quantification kit (Kapa Biosystems, Inc.) was used for quantification of the cDNA library. Finally, the RNA-seq libraries were sequenced on the HiSeq X10 system (Illumina, Inc., San Diego, CA, USA), and 150 bp paired-end (PE150) sequencing was performed on all samples.

Identification of circRNAs.
Raw reads with adaptors and low-quality tags were removed, and the remaining clean reads were used in the subsequent analyses [21]. First, using Bowtie2 version 2.1.0, the clean paired-end reads were mapped to the latest UCSC transcript set [22], and using RSEM v1.2.15, gene expression levels were estimated [23]. CircRNA expression was normalized using TMM (trimmed mean of M values, CPM > 5). For circRNA expression analysis, the reads were firstly mapped to the genome using STAR [24], and then, the circRNA expression was identified and estimated using DCC [25]. The edgeR program (Bioconductor V3.0; Fred Hutchinson Cancer Research Center, Seattle, WA, USA) was used to identify differentially expressed cir-cRNAs [26]. As previously described [27], we controlled the false discovery rate to yield the q value (adjusted p value) to obtain the p values and determine the threshold of the p value. In addition, the spliced reads per billion mapping (SRPBM) value was used to estimate the fold change of cir-cRNA expression in each sample. CircRNAs meeting the condition of p < 0:05 and fold change ≥ 1:5 were considered to be differentially expressed [18].
2.4. GO and KEGG Pathway Analyses. As previously reported [8], circRNAs are alternatively transcribed from their parental genes, and according to the location information of cir-cRNAs, parental genes that could be regulated by circRNAs were obtained. The potential functions of the parental genes corresponding to the differentially expressed circRNAs were predicted by GO (http://geneontology.org/), and the related pathways were analyzed based on the latest KEGG (https:// www.genome.jp/kegg/pathway.html) database. Significant correlations between the parental genes of differentially expressed circRNAs and their potential functions and pathways were determined based on a threshold of p < 0:05 and overlap gene count ≥ 2 [28,29].
2.5. Bioinformatics Analysis of circRNA-miRNA Regulatory Networks. Recently, circRNAs were considered miRNA sponges that regulate gene expression [8,9]. Therefore, we used the miRbase (http://www.mirbase.org) and miRanda (http://www.microrna.org/microrna/home.do) databases to predict potential miRNAs that are associated with the differentially expressed circRNAs in AS [30,31]. Based on the interactions between the differentially expressed circRNAs and their target miRNAs, circRNA-miRNA interaction networks were constructed and visualized with Cytoscape version 3.7.1.
2.6. Statistical Analysis. Statistical analysis was performed using SPSS 19.0 (SPSS Inc., Chicago, IL, USA). Results were expressed as the mean ± standard deviation. Comparisons between the two groups were tested using two-tailed Student's t test. Differences were considered statistically significant at p < 0:05.

Identification of Differentially Expressed circRNAs in AS.
To identify differentially expressed circRNAs in patients with AS, RNA sequencing was used for circRNA expression profiling in the spinal ligament tissue samples from three AS patients and three patients with lumbar disc herniation. The main characteristics of all of the patients were presented in Table 1. Sequencing identified a total of 1,172 circRNAs, and the details on total circRNA expression profiling were shown in supplementary data (available here). Approximately 98% of the identified circRNAs were <2,000 nucleotides (nt) in length, and the length of most differentially expressed circRNAs with statistical significance between the two groups was also <2,000 nt ( Figure 1(a)). The distribution of identified circRNAs on the human chromosomes is shown in Figure 1(b). Briefly, no differentially expressed circRNAs were observed on chr15 and chr20, and only the upregulated circRNAs from chr12, chr18, chr21, and chrY were transcribed, whereas the upregulated and downregulated cir-cRNAs were transcribed from all other chromosomes except for chr 22, which only had one downregulated cir-cRNA ( Figure 1(b)). The scatter plot in Figure 1(c) depicts the variations in circRNA expression levels between the two groups. As shown in the volcano plot, the circRNAs with fold changes ≥ 1:5 and p value < 0.05 were considered to be significant differentially expressed between the two groups, and among these differentially expressed circRNAs, 57 cir-cRNAs were upregulated and 66 circRNAs were downregulated ( Figure 1(d)). The details of the top 10 circRNAs showing the significant upregulation or downregulation are shown in Table 2. Besides, hierarchical clustering revealed that circRNA expression levels were distinguishable in Figure 1(e). In addition, the vast majority of differentially expressed circRNAs in AS are located in the exonic region ( Figure 1(f)). Based on the sequence structure origin of the circRNAs, 114 of the differentially expressed circRNAs in AS were exon-exon circRNAs, 3 were exon-intergenic cir-cRNAs, 4 were intergenic-exon circRNAs, and 2 were intergenic-intergenic circRNAs. In summary, RNA-seq analysis suggests that the expression levels of some circRNAs in the spinal ligament tissue samples had significantly changed between patients with AS and the controls.

GO and KEGG Pathway Analyses of the Differentially Expressed circRNAs in AS.
To explore the functions of differentially expressed circRNAs based on their parental genes in AS, GO and KEGG pathway analyses were conducted to predict the possible functions of circRNAs (p < 0:05). GO analysis comprised three major domains: "biological process" (BP), "cell component" (CC), and "molecular function" (MF). The top five enriched GO terms of differentially expressed cir-cRNAs in each domain were "regulation of GTPase activity", "regulation of cell morphogenesis", "axon development", "peptidyl-serine phosphorylation" and "peptidyl-serine modification" in BP; "axon part", "cell cortex", "distal axon", "growth cone" and "site of polarized growth" in CC; "protein serine/threonine kinase activity", "nucleoside-triphosphatase regulator activity", "Ras GTPase binding", "GTPase regulator activity", and "phosphatidylinositol-3,4,5-trisphosphate binding" in MF (Figure 2(a)). In addition, the detailed information of the enriched GO terms of differentially expressed circRNAs was shown in supplementary data.
KEGG pathway analysis revealed that the differentially expressed circRNAs in AS might mainly be involved in "Yersinia infection", "human immunodeficiency virus 1 infection", and "Human cytomegalovirus infection" (Figure 2(b)). Forty-nine pathways related to the functions of 123 differentially expressed circRNAs were identified by KEGG analysis (supplementary data).

Construction of circRNA/miRNA Interaction Networks in AS.
Growing evidence shows that circRNAs act as miRNA sponges via miRNA response elements (MREs) that regulate the functions of their target mRNAs. To further evaluate the potential functions of the differentially expressed circRNAs in AS, the circRNA-miRNA networks were predicted using the miRBase and miRanda databases and visualized with Cytoscape version 3.7.1. Studies suggest that the more MREs of one miRNA on one circRNA sequence, the microRNAs are more likely regulated by the circRNAs [8,9]. In addition, there are many different MREs on one circRNA sequence, and the MRE of one miRNA also exists in a variety of cir-cRNAs' sequence. Given that the networks consisting of the top 300 circRNA-miRNAs were constructed in previous studies [32,33], the top 300 networks of circRNA/miRNA were selected in this study, and we found that 60 of 123 significantly differentially expressed circRNAs were predicted to interact with 221 miRNAs. The circRNA/miRNA networks were established to clarify the interactions between cir-cRNAs and their target miRNAs (Figure 3). The networks Number of circRNA  show that hsa_circMNT_002 had five MREs of hsa-miR-6722-3p, and hsa-miR-1972, hsa-miR-4706, hsa-miR-6756-5p, and hsa-miR-6812-5p were found to be regulated by four circRNAs, which is higher than that of other miRNAs. Furthermore, hsa_circRUSC2002, hsa_circMNT002, and hsa_ circNFATC1001 were predicted to have complementary binding sites for 45, 24, and 22 miRNAs, respectively, which suggest that the differentially expressed circRNAs play important roles in AS.

Discussion
In this study, we identified numerous circRNAs that are significant differentially expressed in the spinal ligament tissues of AS patients. RNA-seq analysis showed that 123 circRNAs were differentially expressed in the spinal ligament tissues from AS patients compared with the controls. Among these, 57 circRNAs were upregulated and 66 were downregulated. In addition, we performed GO and KEGG pathway analyses and established circRNA/miRNA interaction networks to predict the potential functions of differentially expressed cir-cRNAs in AS.
AS is a chronic autoimmune disease that is characterized by inflammation and pathological osteogenesis and is associated with HLB-27 and T cells [34,35]. Although many studies have significantly improved our understanding of the pathophysiology of AS, its etiology and pathogenesis remain unclear. Recent studies have indicated that a numerous molecular and biochemical changes are involved in the cellular mechanisms of AS [3], and increasing evidence suggests   BioMed Research International  The KEGG enrichment scatter plot of differentially expressed circRNAs with the 10 highest enrichment scores (p value < 0.05 and overlap gene count ≥ 2). The horizontal axis is the enrichment score for each KEGG term, and the vertical axis is the KEGG pathway names. 8 BioMed Research International that ncRNAs, such as lncRNAs and miRNAs, play important roles in the pathogenesis of AS. Some studies have shown that significant changes in the expression of lncRNAs and miRNAs in AS, and the dysfunction of certain lncRNAs and miRNAs has been associated with the pathophysiology of AS. For example, it has been shown that four lncRNAs, namely, lnc-ZNF354A-1, lnc-LIN54-1, lnc-FRG2C-3, and lnc-USP50-2 are involved in osteogenic differentiation of mesenchymal stem cells (MSCS) in patients with AS [36]. A recent study has revealed that lncRNA-AK001085 is downregulated in patients with AS and considered to be a potential diagnostic molecular marker [37]. In addition, the expression of miR-29a, miR-335-5p, miR-27a, and let-7i is upregulated in peripheral blood monocytes of patients with AS [14], and serum miR-146a, miR-29a, and miR-155 levels in AS patients have also significantly increased [15]. In recent years, more studies have focused on circRNAs. However, to date, information on the function of circRNAs in the spinal ligament tissues of patients with AS is limited. To the best of our understanding of the regulatory mechanism of cir-cRNAs in AS, this is the first study that has investigated changes in circRNA expression profiles in the spinal ligament tissues of patients with AS by RNA-seq and bioinformatics analyses, which could be an important step in elucidating the underlying molecular mechanisms associated with AS. circRNAs are highly capable of resisting damage caused by RNase and, thus, make them more stable than linear RNAs. Moreover, circRNAs have highly conserved tissue and cell specificity [11]. Recent studies have shown that cir-cRNAs, such as circular RNA Atp9b [38] and circRNA hsa_circ_0005105 [39], are involved in osteoarthritis [40]. Ligaments connect the ends of bones and allow stability and mobility of most joints. Previous studies have revealed that ligaments may be the main target tissue for inflammation   BioMed Research International and ossification in AS and are responsible for many symptoms of AS patients [17,41,42]. Based on the results of this study, we hypothesize that the differential expression of circRNAs in the ligament tissues of AS patients is more specific and sensitive than those in peripheral blood, and thus, we collected spinal ligament tissues of AS patients for circRNA profiling by RNA-seq. We detected a total of 1,172 circRNAs that are located in 46 chromosomes in the spinal ligament tissue samples, and approximately 98% of the identified circRNAs were <2,000 nt in length (Figure 1(a)). Among these differentially expressed circRNAs, 57 circRNAs were upregulated and 66 circRNAs were downregulated (p value < 0.05 and fold change ≥ 1:5; Figure 1(d)), and the length of most differentially expressed circRNAs was also <2,000 nt (Figure 1(a)). Recent studies have shown that circRNAs are alternatively transcribed from exons, introns, or other regions of their parental genes [43]. The results of sequencing have shown that the vast majority of differentially expressed circRNAs in AS are located in exonic region (Figure 1(f)). Previous studies have also demonstrated that various types of circRNAs have different functions [44]. The canonical pathway of circRNA biologic regulation involves acting as sponges of miRNAs [8,9]. Some circRNAs can affect protein function by directly binding to them [6,45,46] and translated into proteins [47,48].
Because circRNAs have been associated with various human diseases [44], these are considered potential clinical diagnostic and therapeutic molecules. To further investigate the circRNA functions of differentially expressed circRNAs in AS, GO and KEGG analyses were performed to predict their possible biological functions and underlying mechanisms. The results of GO analyses have shown that the differentially expressed circRNAs in AS are enriched in BP and MF, such as "peptidyl-serine phosphorylation", "protein serine/threonine kinase activity", and "peptidyl-threonine phosphorylation" (Figure 2(a)), which may be involved in the physiological functions of mitogen-activated protein kinase (MAPK) signaling pathways and phosphatidylinositol 3pinases (PI3K)/Akt signaling pathways that are associated with inflammation-induced apoptosis in chondrocytes, such as IL-1β [49]. We found that "regulation of GTPase activity", "GTPase regulator activity", and "Ras GTPase binding" are also enriched in BP and MF (Figure 2(a)), while there are few studies about the relationship between GTPase and AS. However, previous studies have shown that GTPase are associated with human diseases [50], which suggested that GTPase may be involved in AS. Moreover, KEGG analysis demonstrated that differentially expressed circRNAs might be related to "Yersinia infection", "Human immunodeficiency virus 1 infection", and "Human cytomegalovirus infection", which are associated with the immune system (Figure 2(b)). In addition, "T cell receptor signaling pathway" and "Th17 cell differentiation" was 11th and 12th in the results of KEGG analysis, respectively (supplementary data). However, previous studies have shown that AS is an autoimmune disorder that is associated with Th17 cells and immune pathways, including the interleukin-(IL-) 17/IL-23 pathway [51,52], and Th17 cells derived from CD4 + T cells can release various kinds of cytokines, such as IL-17A and IL-22, which result in bone erosion/proliferation in AS [3], whereas the differentia-tion of Th17 cells is influenced by multiple inflammatory cytokines, including IL-1β, IL-6, TGF-β, and IL-23 [52].
Current studies have confirmed that circRNAs in the cytoplasm may play crucial roles in posttranscriptional gene regulation by sponging miRNAs as competing endogenous RNAs or RNA-binding proteins [6,8,9,46]. For example, ciRS-7 has been reported as a miR-7 sponge that is involved in osteosarcoma [53] and affected brain function [54]. The results of the present identified numerous differentially expressed circRNAs in the cytoplasm (Figure 2(a)), and the circRNA/miRNA interaction network analysis was established to further demonstrate the potential functions of the differentially expressed circRNAs. We have also predicted the interactions of the differentially expressed circRNAs in AS with miRNAs by bioinformatics analysis. The top 300 networks of circRNA-miRNA were selected, and we predicted that 60 of 123 significantly differentially expressed circRNA interaction with 221 miRNAs (Figure 3). The networks show that hsa_circMNT_002 has five MREs of hsa-miR-6722-3p, and hsa-miR-1972, hsa-miR-4706, hsa-miR-6756-5p, and hsa-miR-6812-5p are regulated by four circRNAs, which is higher than that of other miRNAs. Furthermore, hsa_circRUSC2002, hsa_circMNT002, and hsa_ circNFATC1001 were predicted to have complementary binding sites for 45, 24, and 22 miRNAs, respectively (Figure 3), which suggests that these circRNAs might be associated with the pathogenesis of AS by regulating miRNAs.

Conclusion
This study has provided the first evidence of circRNA expression profiles in the spinal ligament tissue of patients with AS using RNA-seq, and the results indicated that numerous differentially expressed circRNAs may be associated with AS. Furthermore, the potential functions of these circRNAs were investigated by bioinformatics analysis and circRNA/miRNA interaction networks were constructed. These findings may provide novel clues for understanding the mechanisms and have identified novel potential molecular targets for the diagnoses and treatment of AS. However, to better understand the role of circRNAs in AS, additional investigations are warranted.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.