BioMed Research International
Volume 2020 (2020), Article ID 7165893, 12 pages
https://doi.org/10.1155/2020/7165893

Research Article

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

Jianqiang Kou,¹ Guoming Liu,¹ Xiangyun Liu,¹ Tianmi Li,² Ying Wei,² Yuanliang Sun,¹ Ting Wang,¹ Yingzhen Wang,¹ and Xiujun Zheng¹

¹Department of Orthopedics, The Affiliated Hospital of Qingdao University, Qingdao, Shandong 266000, China
²Department of Operating Room, The Affiliated Hospital of Qingdao University, Qingdao, Shandong 266000, China

Correspondence should be addressed to Xiujun Zheng; qyfykjq@126.com

Received 3 March 2020; Revised 8 May 2020; Accepted 27 May 2020; Published 15 June 2020

Guest Editor: Jiangning Song

Copyright © 2020 Jianqiang Kou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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

1. 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 tendon-bone interface, leading to bone proliferation [2–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–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 circRNAs 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.

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

Table 1: Personal characteristics of participants in the sequencing.

Items
Control group
Experimental group

C1
C2
C3
AS1
AS2
AS3

Age (years)
58
69
64
64
58
76

Sex
Female
Male
Male
Male
Male
Male

Clinical symptoms
Lumbar disc herniation
Lumbar disc herniation
Lumbar disc herniation
AS for 40 years, lumbar kyphosis deformity
AS for 8 years, thoracic fracture
AS for 30 years, thoracic fracture and kyphosis deformity

Procedures
Posterior lumbar decompression and fusion
Posterior lumbar decompression and fusion
Posterior lumbar decompression and fusion
Posterior lumbar decompression and fusion
Posterior thoracic decompression and fusion
Posterior thoracic decompression and fusion

HLA B27⁺
No
No
No
Yes
Yes
Yes

Bilateral damage of the sacroiliac joint
No
No
No
Yes
Yes
Yes

Spinal and sacroiliac joint fusion
No
No
No
Yes
Yes
Yes

Non-steroidal anti-inflammatory drugs
No
No
No
No
No
No

Other type of autoimmune disease.
No
No
No
No
No
No

Complications
No
No
No
No
No
No

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 first-strand 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.

2.3. 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 values, ). 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 circRNAs [26]. As previously described [27], we controlled the false discovery rate to yield the value (adjusted value) to obtain the values and determine the threshold of the value. In addition, the spliced reads per billion mapping (SRPBM) value was used to estimate the fold change of circRNA expression in each sample. CircRNAs meeting the condition of and fold 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 circRNAs, 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 and overlap gene [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 deviation. Comparisons between the two groups were tested using two-tailed Student’s test. Differences were considered statistically significant at .

3. Results

3.1. 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 circRNAs were transcribed from all other chromosomes except for chr 22, which only had one downregulated circRNA (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 and value < 0.05 were considered to be significant differentially expressed between the two groups, and among these differentially expressed circRNAs, 57 circRNAs 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 circRNAs, 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Differential expression profiles of circular RNAs (circRNAs) in spinal ligament tissue from patients with AS and controls by RNA-seq. (a) Length distribution of circRNAs was identified in this study. Orange represents differentially expressed circRNAs with statistical significance ( value < 0.05 and fold ), and blue represents the circRNAs with no any changes in expression between the two groups (). (b) Distribution of identified circRNAs on human chromosomes. The length of the lines represents the relative size of fold change, and red and blue indicate the upregulated and downregulated circRNAs, respectively. (c) The scatter plot of circRNA expression between the two groups. The dots in the figure represents circRNAs. Compared with the controls, the red dots and green dots, respectively, indicate the upregulated and downregulated circRNA expressions in the spinal ligament tissues of AS patients ( value < 0.05 and fold ), and the blue dots indicate no change in circRNA expression. (d) The volcano plot depicting circRNA expression profiles in the spinal ligament tissues of the two groups. The two vertical blue lines refer to a 1.5-fold (log2 fold change) upregulation and downregulation, respectively. The horizontal blue line corresponds to a value of 0.05 (−log10 value). Compared with the controls, the red and green dots, respectively, indicate upregulated and downregulated circRNA expressions in the spinal ligament tissues of patients with AS ( value < 0.05 and fold ). (e) Hierarchical clustering analysis of circRNAs that were significantly differentially expressed ( value < 0.05 and fold ) in the spinal ligament tissue between the two groups. (f) Categories of differentially expressed circRNAs based on their genomic origin.

Table 2: The top 10 circRNAs with the significant upregulation or downregulation.

Position
circBankID
circBaseID
Gene symbol
Length
Type
log2FC
value

Regulation

chr11_92085261_92088570_+
hsa_circFAT3_006
hsa_circ_0000348
FAT3
3309
Exon-exon
3.7819
0.0000
0.0058
Up

chr13_76195898_76335174_+
NA
NA
LMO7
NA
Exon-exon
3.1527
0.0005
0.0639
Up

chrX_147743428_147744289_+
hsa_circAFF2_005
hsa_circ_0001947
AFF2
861
Exon-exon
2.9490
0.0000
0.0000
Up

chr8_13356557_13357705_-
hsa_circDLC1_002
hsa_circ_0135780
DLC1
1148
Exon-exon
2.3865
0.0005
0.0639
Up

chr2_200233327_200298237_-
hsa_circSATB2_015
hsa_circ_0003915
SATB2
531
Exon-exon
2.2778
0.0013
0.1144
Up

chr9_4286037_4286523_-
hsa_circGLIS3_003
hsa_circ_0002874
GLIS3
486
Exon-exon
2.1186
0.0000
0.0047
Up

chr12_122825299_122826244_-
hsa_circCLIP1_019
hsa_circ_0029069
CLIP1
945
Exon-exon
2.0273
0.0001
0.0123
Up

chr3_18419661_18462483_-
hsa_circSATB1_004
hsa_circ_0064555
SATB1
1599
Exon-exon
1.9818
0.0084
0.2278
Up

chr21_36206706_36231875_-
hsa_circRUNX1_005
hsa_circ_0002360
RUNX1
297
Exon-exon
1.8979
0.0066
0.2154
Up

chr9_37424841_37426651_+
hsa_circGRHPR_002
hsa_circ_0001861
GRHPR
321
Exon-exon
1.8159
0.0008
0.0828
Up

chr6_144858717_144864006_+
hsa_circUTRN_065
hsa_circ_0130908
UTRN
362
Exon-exon
-2.1886
0.0082
0.2278
Down

chr19_4511522_4511918_-
NA
NA
PLIN4
NA
Exon-exon
-2.1979
0.0022
0.1357
Down

chr8_17570722_17581342_-
hsa_circMTUS1_004
hsa_circ_0083443
MTUS1
336
Exon-exon
-2.3338
0.0032
0.1500
Down

chr5_102432244_102433485_-
hsa_circGIN1_005
hsa_circ_0006499
GIN1
655
Exon-exon
-2.6637
0.0019
0.1295
Down

chr1_154207066_154207767_+
hsa_circUBAP2L_011
hsa_circ_0110845
UBAP2L
298
Exon-exon
-3.1140
0.0027
0.1426
Down

chr9_14120438_14179779_-
hsa_circNFIB_012
hsa_circ_0138300
NFIB
683
Exon-exon
-3.1464
0.0002
0.0326
Down

chr10_34558584_34573173_-
hsa_circPARD3_036
hsa_circ_0018168
PARD3
354
Exon-exon
-3.1686
0.0020
0.1295
Down

chr3_100565215_100581229_-
hsa_circABI3BP_001
hsa_circ_0121334
ABI3BP
513
Exon-exon
-3.4042
0.0000
0.0000
Down

chr3_100565215_100570784_-
NA
NA
ABI3BP
NA
Exon-exon
-4.1847
0.0000
0.0000
Down

chr17_58012553_58018304_+
NA
NA
RPS6KB1
NA
Exon-exon
-5.5909
0.0000
0.0079
Down

Position: chromosome location of circRNA; log2FC: log2Fold change; : adjusted value; NA: not available.

3.2. 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 (). 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 circRNAs 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.

(a)

(b)

Figure 2: GO and KEGG pathway analyses of differentially expressed circRNA. (a) GO analysis of differentially expressed circRNAs ( value < 0.05 and overlap gene ). The horizontal axis is the enrichment score for each GO term, and the vertical axis is the GO term. The enrichment score was calculated as -log10 ( value). The number on the histogram columns represents the number of overlapping genes. (b) The KEGG enrichment scatter plot of differentially expressed circRNAs with the 10 highest enrichment scores ( value < 0.05 and overlap gene ). The horizontal axis is the enrichment score for each KEGG term, and the vertical axis is the KEGG pathway names.

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

3.3. 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 circRNAs’ 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 circRNAs and their target miRNAs (Figure 3). The networks 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.

Figure 3: circRNA/miRNA interaction networks for differentially expressed circRNAs. The interaction networks of differentially expressed circRNAs and their complementary binding miRNAs were visualized using Cytoscape 3.7.1 software. The red and green triangles represent upregulated and downregulated circRNAs, respectively, and the blue spots indicate microRNAs. The networks included 60 differentially expressed circRNAs and 221 miRNAs.

4. 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 circRNAs 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 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 circRNAs 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 circRNAs, 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 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 ( value < 0.05 and ; 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 3-pinases (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 differentiation 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.

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

Conflicts of Interest

The authors have no conflicts of interest to declare.

Authors’ Contributions

YW and XZ conceived and designed the experiments. All authors performed the experiments. JK, GL, and XZ analyzed the data. JK and GL were significant contributors in the manuscript. All authors revised the manuscript and have agreed to the publication of this manuscript.

Acknowledgments

The authors would like to thank Guangzhou Geneseed Biotech Co., Ltd. (Guangzhou, China) for their assistance during high-throughput sequencing and data analysis and also thank LetPub (http://www.letpub.com) for its linguistic assistance during the preparation of this manuscript. This research was supported by the Special Fund for Youth of Applied Foundational Research Program of Qingdao (no. 18-2-2-44-jch).

Supplementary Materials

Sequencing identified a total of 1,172 circRNAs, and the details on total circRNA expression profiling were shown in supplementary data. (Supplementary Materials)

References

W. P. Maksymowych, “Disease modification in ankylosing spondylitis,” Nature Reviews Rheumatology, vol. 6, no. 2, pp. 75–81, 2010.

R. Landewe, M. Dougados, H. Mielants, H. van der Tempel, and D. van der Heijde, “Physical function in ankylosing spondylitis is independently determined by both disease activity and radiographic damage of the spine,” Annals of the Rheumatic Diseases, vol. 68, no. 6, pp. 863–867, 2009.

D. Simone, M. H. Al Mossawi, and P. Bowness, “Progress in our understanding of the pathogenesis of ankylosing spondylitis,” Rheumatology, vol. 57, Supplement 6, pp. vi4–vi9, 2018.

A. El Maghraoui, “Extra-articular manifestations of ankylosing spondylitis: prevalence, characteristics and therapeutic implications,” European Journal of Internal Medicine, vol. 22, no. 6, pp. 554–560, 2011.

M. A. Brown and B. P. Wordsworth, “Genetics in ankylosing spondylitis - current state of the art and translation into clinical outcomes,” Best Practice & Research Clinical Rheumatology, vol. 31, no. 6, pp. 763–776, 2017.

R. Ashwal-Fluss, M. Meyer, N. R. Pamudurti et al., “circRNA biogenesis competes with pre-mRNA splicing,” Molecular Cell, vol. 56, no. 1, pp. 55–66, 2014.

M. T. Hsu and M. Coca-Prados, “Electron microscopic evidence for the circular form of RNA in the cytoplasm of eukaryotic cells,” Nature, vol. 280, no. 5720, pp. 339-340, 1979.

S. Memczak, M. Jens, A. Elefsinioti et al., “Circular RNAs are a large class of animal RNAs with regulatory potency,” Nature, vol. 495, no. 7441, pp. 333–338, 2013.

T. B. Hansen, T. I. Jensen, B. H. Clausen et al., “Natural RNA circles function as efficient microRNA sponges,” Nature, vol. 495, no. 7441, pp. 384–388, 2013.

S. Qu, X. Yang, X. Li et al., “Circular RNA: a new star of noncoding RNAs,” Cancer Letters, vol. 365, no. 2, pp. 141–148, 2015.

W. Chen and E. Schuman, “Circular RNAs in brain and other tissues: a functional enigma,” Trends in Neurosciences, vol. 39, no. 9, pp. 597–604, 2016.

L. Santer, C. Bar, and T. Thum, “Circular RNAs: a novel class of functional RNA molecules with a therapeutic perspective,” Molecular Therapy, vol. 27, no. 8, pp. 1350–1363, 2019.

Y. Wang, J. Liu, J. Ma et al., “Exosomal circRNAs: biogenesis, effect and application in human diseases,” Molecular Cancer, vol. 18, no. 1, p. 116, 2019.

W. Yang, X. Yan, Q. Xia et al., “Predisposition of six well-characterized microRNAs to syndesmophytes among Chinese patients with ankylosing spondylitis,” Modern Rheumatology, vol. 29, no. 1, pp. 173–180, 2018.

B. P. Qian, M. L. Ji, Y. Qiu et al., “Identification of serum miR-146a and miR-155 as novel noninvasive complementary biomarkers for ankylosing spondylitis,” Spine, vol. 41, no. 9, pp. 735–742, 2016.

Z. Xu, X. Zhou, H. Li, Q. Chen, and G. Chen, “Identification of the key genes and long non‑coding RNAs in ankylosing spondylitis using RNA sequencing,” International Journal of Molecular Medicine, vol. 43, no. 3, pp. 1179–1192, 2019.

S. van der Linden, H. A. Valkenburg, and A. Cats, “Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria,” Arthritis and Rheumatism, vol. 27, no. 4, pp. 361–368, 1984.

H. Xu, C. Wang, H. Song, Y. Xu, and G. Ji, “RNA-seq profiling of circular RNAs in human colorectal cancer liver metastasis and the potential biomarkers,” Molecular Cancer, vol. 18, no. 1, p. 8, 2019.

L. P. Wang, X. Y. Peng, X. Q. Lv et al., “High throughput circRNAs sequencing profile of follicle fluid exosomes of polycystic ovary syndrome patients,” Journal of Cellular Physiology, vol. 234, no. 9, pp. 15537–15547, 2019.

M. Zhao, J. Xu, S. Zhong et al., “Expression profiles and potential functions of circular RNAs in extracellular vesicles isolated from radioresistant glioma cells,” Oncology Reports, vol. 41, no. 3, pp. 1893–1900, 2019.

K. P. Zhu, C. L. Zhang, X. L. Ma, J. P. Hu, T. Cai, and L. Zhang, “Analyzing the interactions of mRNAs and ncRNAs to predict competing endogenous RNA networks in osteosarcoma chemo-resistance,” Molecular Therapy, vol. 27, no. 3, pp. 518–530, 2019.

B. Langmead and S. L. Salzberg, “Fast gapped-read alignment with Bowtie 2,” Nature Methods, vol. 9, no. 4, pp. 357–359, 2012.

B. Li and C. N. Dewey, “RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome,” BMC Bioinformatics, vol. 12, no. 1, p. 323, 2011.

A. Dobin, C. A. Davis, F. Schlesinger et al., “STAR: ultrafast universal RNA-seq aligner,” Bioinformatics, vol. 29, no. 1, pp. 15–21, 2013.

J. Cheng, F. Metge, and C. Dieterich, “Specific identification and quantification of circular RNAs from sequencing data,” Bioinformatics, vol. 32, no. 7, pp. 1094–1096, 2016.

M. D. Robinson, D. J. McCarthy, and G. K. Smyth, “edgeR: a bioconductor package for differential expression analysis of digital gene expression data,” Bioinformatics, vol. 26, no. 1, pp. 139-140, 2009.

Z. Zhang, N. Song, Y. Wang et al., “Analysis of differentially expressed circular RNAs for the identification of a coexpression RNA network and signature in colorectal cancer,” Journal of Cellular Biochemistry, vol. 120, no. 4, pp. 6409–6419, 2019.

C. Zhang, J. Liu, M. Lai et al., “Circular RNA expression profiling of granulosa cells in women of reproductive age with polycystic ovary syndrome,” Archives of Gynecology and Obstetrics, vol. 300, no. 2, pp. 431–440, 2019.

Y. Ou, M. Liu, L. Zhu et al., “The expression profile of circRNA and its potential regulatory targets in the placentas of severe pre-eclampsia,” Taiwanese Journal of Obstetrics & Gynecology, vol. 58, no. 6, pp. 769–777, 2019.

A. E. Pasquinelli, “MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship,” Nature Reviews Genetics, vol. 13, no. 4, pp. 271–282, 2012.

X. Gu, M. Li, Y. Jin, D. Liu, and F. Wei, “Identification and integrated analysis of differentially expressed lncRNAs and circRNAs reveal the potential ceRNA networks during PDLSC osteogenic differentiation,” BMC Genetics, vol. 18, no. 1, p. 100, 2017.

L. Li, J. Guo, Y. Chen, C. Chang, and C. Xu, “Comprehensive CircRNA expression profile and selection of key CircRNAs during priming phase of rat liver regeneration,” BMC Genomics, vol. 18, no. 1, p. 80, 2017.

Y. Yan, R. Zhang, X. Zhang, A. Zhang, Y. Zhang, and X. Bu, “RNA-Seq profiling of circular RNAs and potential function of hsa_circ_0002360 in human lung adenocarcinom,” American Journal of Translational Research, vol. 11, no. 1, pp. 160–175, 2019.

D. M. Evans, C. C. Spencer, J. J. Pointon et al., “Interaction between ERAP1 and HLA-B27 in ankylosing spondylitis implicates peptide handling in the mechanism for HLA-B27 in disease susceptibility,” Nature Genetics, vol. 43, no. 8, pp. 761–767, 2011.

C. Jandus, G. Bioley, J. P. Rivals, J. Dudler, D. Speiser, and P. Romero, “Increased numbers of circulating polyfunctional Th17 memory cells in patients with seronegative spondylarthritides,” Arthritis and Rheumatism, vol. 58, no. 8, pp. 2307–2317, 2008.

Z. Xie, J. Li, P. Wang et al., “Differential expression profiles of long noncoding RNA and mRNA of osteogenically differentiated mesenchymal stem cells in ankylosing spondylitis,” The Journal of Rheumatology, vol. 43, no. 8, pp. 1523–1531, 2016.

X. Li, W. Chai, G. Zhang et al., “Down-regulation of lncRNA-AK001085 and its influences on the diagnosis of ankylosing spondylitis,” Medical Science Monitor, vol. 23, pp. 11–16, 2017.

Z. B. Zhou, D. Du, G. X. Huang, A. Chen, and L. Zhu, “Circular RNA Atp9b, a competing endogenous RNA, regulates the progression of osteoarthritis by targeting miR-138-5p,” Gene, vol. 646, pp. 203–209, 2018.

Y. Wu, Y. Zhang, Y. Zhang, and J. J. Wang, “CircRNA hsa_circ_0005105 upregulates NAMPT expression and promotes chondrocyte extracellular matrix degradation by sponging miR-26a,” Cell Biology International, vol. 41, no. 12, pp. 1283–1289, 2017.

Q. Liu, X. Zhang, X. Hu et al., “Emerging roles of circRNA related to the mechanical stress in human cartilage degradation of osteoarthritis,” Mol Ther Nucleic Acids, vol. 7, pp. 223–230, 2017.

M. Benjamin and D. McGonagle, “Entheses: tendon and ligament attachment sites,” Scandinavian Journal of Medicine & Science in Sports, vol. 19, no. 4, pp. 520–527, 2009.

T. J. Kim, T. H. Kim, J. B. Jun, K. B. Joo, and W. S. Uhm, “Prevalence of ossification of posterior longitudinal ligament in patients with ankylosing spondylitis,” The Journal of Rheumatology, vol. 34, no. 12, pp. 2460–2462, 2007.

J. Salzman, C. Gawad, P. L. Wang, N. Lacayo, and P. O. Brown, “Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types,” PLoS One, vol. 7, no. 2, article e30733, 2012.

B. Han, J. Chao, and H. Yao, “Circular RNA and its mechanisms in disease: from the bench to the clinic,” Pharmacology & Therapeutics, vol. 187, pp. 31–44, 2018.

W. W. Du, C. Zhang, W. Yang, T. Yong, F. M. Awan, and B. B. Yang, “Identifying and characterizing circRNA-protein interaction,” Theranostics, vol. 7, no. 17, pp. 4183–4191, 2017.

W. W. Du, W. Yang, E. Liu, Z. Yang, P. Dhaliwal, and B. B. Yang, “Foxo3 circular RNA retards cell cycle progression via forming ternary complexes with p21 and CDK2,” Nucleic Acids Research, vol. 44, no. 6, pp. 2846–2858, 2016.

Y. Yang, X. Gao, M. Zhang et al., “Novel role of FBXW7 circular RNA in repressing glioma tumorigenesis,” JNCI: Journal of the National Cancer Institute, vol. 110, no. 3, pp. 304–315, 2018.

Y. Yang, X. Fan, M. Mao et al., “Extensive translation of circular RNAs driven by N6-methyladenosine,” Cell Research, vol. 27, no. 5, pp. 626–641, 2017.

P. F. Hu, F. F. Sun, and J. Qian, “Leonurine exerts anti-catabolic and anti-apoptotic effects via nuclear factor kappa B (NF-κB) and mitogen-activated protein kinase (MAPK) signaling pathways in chondrocytes,” Medical Science Monitor, vol. 25, pp. 6271–6280, 2019.

N. A. Guadagno and C. Progida, “Rab GTPases: switching to human diseases,” Cells, vol. 8, no. 8, p. 909, 2019.

A. Ridley, H. Hatano, I. Wong-Baeza et al., “Activation-induced killer cell immunoglobulin-like receptor 3DL2 binding to HLA-B27 licenses pathogenic T cell differentiation in spondyloarthritis,” Arthritis & Rhematology, vol. 68, no. 4, pp. 901–914, 2016.

N. Vanaki, S. Aslani, A. Jamshidi, and M. Mahmoudi, “Role of innate immune system in the pathogenesis of ankylosing spondylitis,” Biomedicine & Pharmacotherapy, vol. 105, pp. 130–143, 2018.

B. Xu, T. Yang, Z. Wang, Y. Zhang, S. Liu, and M. Shen, “CircRNA CDR1as/miR-7 signals promote tumor growth of osteosarcoma with a potential therapeutic and diagnostic value,” Cancer Management and Research, vol. Volume 10, pp. 4871–4880, 2018.

M. Piwecka, P. Glažar, L. R. Hernandez-Miranda et al., “Loss of a mammalian circular RNA locus causes miRNA deregulation and affects brain function,” Science, vol. 357, no. 6357, article eaam8526, 2017.