Circular RNAs in peripheral blood mononuclear cells from ankylosing spondylitis.

BACKGROUND: Circular RNA (circRNA) is a type of closed circular noncoding RNA (ncRNA), mostly formed by back-splicing or alternative splicing of pre-messenger RNA (mRNA). The aim of this study was to explore the expression profile of circRNA in peripheral blood mononuclear cells (PBMCs) of patients with ankylosing spondylitis (AS) and discover potential molecular markers of AS.
METHODS: The circRNA microarray technology was used to detect the expression of circRNAs in the peripheral blood of 6 patients with AS and 6 healthy controls (HC). To screen the differentially expressed circRNAs by fold change (FC) and P value, these differentially expressed circRNAs were analyzed by bioinformatics. In 60 cases of AS and 30 cases of HC, 4 circRNAs were subjected to real-time fluorescence quantitative polymerase chain reaction (RT-qPCR), and their correlation with various clinical indicators was analyzed. Finally, the receiver operating characteristic (ROC) curve was used to analyze their potential as AS diagnostic markers.
RESULTS: The microarray results showed that there were 1369 significantly differently expressed (P < 0.05, FC > 1.5) circRNAs between the AS and HC groups (675 upregulated and 694 downregulated). The results of bioinformatics analysis suggested that they were mainly involved in "enzyme binding," "adenosine ribonucleotide binding," "MAPK signaling pathway", etc. The RT-qPCR results showed that the expressions of hsa_circRNA_001544 (U = 486.5, P < 0.05) and hsa_circRNA_102532 (U = 645, P < 0.05) were significantly different between the AS group and the HC group. The AS group was further divided into two subgroups: active AS (ASA) and stable AS (ASS). After analysis, it was found that compared with the HC group, hsa_circRNA_001544 was significantly increased in both ASA (U = 214, P < 0.05) and ASS groups (U = 273, P < 0.05), while hsa_circRNA_008961 (U = 250, P < 0.05) and hsa_circRNA_102532 (U = 295, P < 0.05) were only significantly increased in the ASA group. Furthermore, hsa_circRNA_012732 was significantly different between the ASA and ASS groups (U = 194, P < 0.05), and there was no statistical significance among the remaining groups. Correlation analysis results showed that hsa_circRNA_012732 was negatively correlated with Bath Ankylosing Spondylitis Disease Activity Index (BASDAI), high-sensitivity C-reactive protein (hsCRP), and globulin (GLOB) and positively correlated with lymphocyte count (LY), mean corpusular volume, and albumin (ALB), and hsa_circRNA_008961 was negatively correlated with platelet (PLT) count. ROC curve analysis showed that hsa_circRNA_001544 (95% CI = 0.610-0.831, P < 0.05) and hsa_circRNA_102532 (95% CI = 0.521-0.762, P < 0.05) were statistically significant, and their area under curve (AUC) values were 0.720 and 0.642, respectively.
CONCLUSIONS: There are differentially expressed circRNAs in PBMCs of AS patients, and they may be involved in the occurrence and development of AS. Among these differentially expressed circRNAs, hsa_circRNA_012732 has the potential to become an indicator of disease activity, and hsa_circRNA_001544 has the potential to become a molecular marker for AS diagnosis.

Introduction

Circular RNA (circRNA) is a type of closed circular noncoding RNA (ncRNA), mostly formed by back-splicing or alternative splicing of pre-messenger RNA (mRNA).[ Its main biological function is to act as a competing endogenous RNA (ceRNA), which combines with microRNA (miRNA) through the miRNA response element (MRE) to release miRNA's inhibitory effect on target genes.[ Also, some intron-derived circRNAs can regulate the expression of their parental genes by positively regulating the transcription of RNA polymerase II.[

Ankylosing spondylitis (AS) is a spondyloarthritis (SpA) with inflammatory back pain and enthesitis as the main clinical manifestations, and it is more common in young males. The pathogenesis of AS is mainly caused by the combination of genetic, environmental, infection, and other factors, but the specific pathogenesis is still unclear. The onset of AS is insidious, and it usually has a delay of 5 to 10 years from onset to diagnosis.[ Patients often cannot get effective treatment in time, resulting in restricted physical activity and even disability.

It has been reported that circRNAs are involved in the occurrence and development of many diseases, including many autoimmune diseases, and can be used as markers for the diagnosis and prognosis of diseases or therapeutic targets.[ A recent study has analyzed the circRNA expression profile of AS spinal ligament tissue by RNA sequencing and found that there are differentially expressed circRNAs between AS and lumbar disc herniation in patients.[ In this regard, we are very curious whether differentially expressed circRNAs can also be found in AS peripheral blood mononuclear cells (PBMCs). Therefore, we used the circRNA microarray technology to detect the expression of circRNAs in the PBMCs of patients with AS, construct its expression profile, and look for circRNAs that have the potential to become AS biomarkers.

Methods

Ethical approval

The study was approved by the Ethics Committee of the Affiliated Hospital of North Sichuan Medical College (No. 2021ER(A)142) and was consistent with the principles of the Declaration of Helsinki. All participants gave informed consent.

Clinical data collection

Sixty AS patients from the Department of Rheumatology and Immunology of the Affiliated Hospital of North Sichuan Medical College from April 2019 to April 2021 were collected as the case group (44 males and 16 females, aged 18.0–55.0 years, with an average of 36.9 ± 10.6 years). Inclusion criteria are as follows: at least 18-year-old patients and those who were diagnosed in accordance with the New York standard revised in 1984[; exclusion criteria included patients with diabetes, cardiovascular diseases, tumors, infectious diseases, other autoimmune diseases, etc. In addition, among the selected AS patients, we divided them into active and stable patients based on the Bath Ankylosing Spondylitis Disease Activity Index (BASDAI) and acute phase reactants. The active AS (ASA) was defined as BASDAI ≥ 6 or 6 > BASDAI > 4, erythrocyte sedimentation rate (ESR) > 22 mm/h or 6 > BASDAI > 4, and high-sensitivity C-reactive protein (hsCRP) level > 9 mg/L; stable AS (ASS) was defined as BASDAI ≤ 4. Self-compiled epidemiological questionnaires were used to collect patient gender, age, course of disease, treatment history, past history, laboratory test results (including blood routine examination, ESR, hsCRP level, albumin [ALB], globulin [GLOB], aspartate aminotransferase [AST], and alanine aminotransferase [ALT]), imaging test results, etc. Meanwhile, the BASDAI, Bath Ankylosing Spondylitis Functional Index (BASFI), and other scales were used to assess the patient's disease activity and functional impairment. At the same time, 30 healthy people (20 males and 10 females, aged 16.0–53.0 years, with an average of 35.8 ± 10.8 years) who were also undergoing physical examination at the Affiliated Hospital of North Sichuan Medical College was collected as the HC group. And the laboratory indicators from the HC group were all normal and included no medical history of cardiovascular and cerebrovascular diseases, diabetes, tumors, or infectious diseases. There was no statistically significant difference in gender and age of all study subjects. The characteristics of AS patients are listed in [Table 1].

Table 1

Characteristics of 60 AS patients.

Items	AS	ASA	ASS
HLA-B27 (+/−)	54/6	29/1	25/5
Uveitis history (+/−)	9/51	6/24	3/27
Family history of AS (+/−)	6/54	4/26	2/28
BASDAI	3.72 ± 1.89	5.20 ± 0.98	2.24 ± 1.31
BASFI	1.97 ± 2.17	3.20 ± 2.34	0.74 ± 0.99
ESR (mm/h)	37.23 ± 30.23	58.83 ± 27.04∗	5.90 ± 5.54
hsCRP (mg/L)	21.25 ± 27.05	36.60 ± 31.14∗	15.63 ± 12.75
WBC (×109/L)	7.48 ± 1.78	7.99 ± 1.62∗	6.98 ± 1.81
GR (×109/L)	4.99 ± 1.57	5.59 ± 1.43∗	4.39 ± 1.50
LY (×109/L)	1.82 ± 0.52	1.78 ± 0.56	1.85 ± 0.47
Mo (×109/L)	0.45 ± 0.14	0.47 ± 0.15	0.42 ± 0.13
RBC (×1012/L)	4.89 ± 0.70	4.86 ± 0.81	4.92 ± 0.58
HGC (g/L)	137.63 ± 18.32	130.90 ± 17.39∗	144.37 ± 16.95
HCT	0.43 ± 0.05	0.42 ± 0.05∗	0.45 ± 0.05
MCV (fL)	88.56 ± 6.63	86.28 ± 8.15∗	90.85 ± 3.50
MCHC (g/L)	319.29 ± 10.58	314.94 ± 10.04∗	323.63 ± 9.39
PLT (×109/L)	253.55 ± 87.63	283.07 ± 105.99∗	224.03 ± 50.88
ALT (U/L)	24.34 ± 17.30	24.66 ± 17.75	24.01 ± 17.13
AST (U/L)	23.43 ± 6.93	22.92 ± 7.24	23.94 ± 6.70
ALB (g/L)	44.19 ± 4.95	42.74 ± 4.83∗	45.65 ± 4.71
GLOB (g/L)	36.87 ± 6.14	40.30 ± 6.08∗	33.44 ± 3.93

Data are expressed as mean ± standard deviation or n.

The difference between the ASA group and ASS group was statistically significant (P < 0.05).

ALT: Alanine aminotransferase; ALB: Albumin; AS: Ankylosing spondylitis; AST: Aspartate aminotransferase; ASA: Active AS; ASS: Stable AS; BASDAI: Bath Ankylosing Spondylitis Disease Activity Index; BASFI: Bath Ankylosing Spondylitis Functional Index; ESR: Erythrocyte sedimentation rate; GR: Neutrophil granule count; GLOB: Globulin; HCT: Hematocrit; HGC: Hemoglobin; hsCRP: High-sensitivity C-reactive protein; LY: Lymphocyte count; MCV: Mean red blood cell volume; MCHC: Mean red blood cell hemoglobin concentration; Mo: Monocytes count; PLT: Platelet; RBC: Erythrocyte red blood cell; WBC: White blood cell count.

RNA extraction and quality testing

Four ml of fasting peripheral blood of all subjects was collected with a heparin anticoagulation tube; human PBMCs were extracted and isolated from peripheral blood; then, total RNAs from PMBCs were extracted according to the RNA extraction procedure provided by Trizol (Takara, Shiga, and Japan), and next, the purity and concentration of RNA was detected by an ultraviolet (UV) spectrophotometer (Thermo Fisher, Sunnyvale, CA, USA). Finally, all the samples were stored at −80°C.

CircRNA microarray

Six AS and 6 HC samples of qualified and age-matched individuals were selected for microarray experiment. The purity and concentration of RNA were detected with NanoDrop ND-1000 (Thermo Fisher), and the integrity of RNA was evaluated by electrophoresis on denatured agarose gel. The circRNAs were amplified and labeled by Arraystar Super RNA labeling kits (Arraystar, Rockville, MD, USA) and then using Arraystar Human circRNA Array v2 (8 × 15K) (Arraystar) for hybridization, and the array was scanned with an Agilent scanner G2505C (Agilent, Santa Clara, CA, USA). The scanned images were analyzed by Agilent feature extraction software (version 11.0.1.1); then, the Limma package in R software was used to perform quantile normalization and subsequent data processing. At last, the volcano plot and hierarchical cluster diagram were used to identify the two groups of differentially expressed circRNAs (circRNA microarray experiment and related data analysis were commissioned by Shanghai Kangcheng Biological Company, Shanghai, China).

Bioinformatics analysis

The differentially expressed circRNAs were analyzed by Gene Ontology (GO) (www.geneontology.org)[ and Kyoto Encyclopedia of Genes and Genomes (KEGG) (www.genome.jp/kegg/)[ pathway analysis to obtain their biological functions and related signal transduction pathways.

Arraystar home-made miRNA target prediction software based on TargetScan[ and miRanda[ were used to predict the circRNA/miRNA interaction.

Real-time fluorescence quantitative polymerase chain reaction (RT-qPCR)

Taking 60 AS patients and 30 HCs as the research subjects, 4 circRNAs with obvious differential expression and target genes predicted to be involved in the pathogenesis of AS were selected for RT-qPCR. Following the instructions of TaKaRa Reverse Transcription Kit (Takara, Japan), 2 μL of total RNA was taken for reverse transcription to synthesize cDNA. Then, RT-qPCR was performed on the Q12K Flex real-time PCR instrument (ABI, USA). The reaction system consisted of 12.5 μL SYBR Green PCR Master Mix, 0.4 μL ROX Dye II, 0.15 μL each of the forward and reverse primers, 2 μL cDNA, and ddH2O as required to make up to 20 μL. The reaction conditions were as follows: 95°C for 30 seconds, 95°C for 5 seconds, and 60°C for 34 seconds, with a total of 40 cycles, and 95°C for 15 seconds, 60°C for 60 seconds, and 95°C for 15 seconds, with a total of one cycle. All samples were provided with multiple holes. The Ct value of the target gene minus the Ct value of the internal reference gene was ΔCt, and the 2−ΔΔCt value represented the relative molecular weight of the target gene to distinguish the difference in target gene expression between different samples. The basic information of four circRNAs is listed in [Table 2], and the primers for RT-qPCR are listed in [Table 3].

Table 2

Basic information of four circRNAs.

CircRNA	Alias	Chrom	Strand	CircRNA type	Best transcript	Gene symbol	P values	FC	Regulation
hsa_circRNA_001544	hsa_circ_0001544	chr5	−	Intronic	ENST00000416954	NR3C1	<0.05	1.60	Up
hsa_circRNA_102532	hsa_circ_0008813	chr19	+	Exonic	NM_005499	UBA2	<0.02	1.58	Up
hsa_circRNA_008961	hsa_circ_0008961	chr2	+	Exonic	NM_005966	NAB1	<0.02	1.53	Up
hsa_circRNA_012732	hsa_circ_0012732	chr1	−	Exonic	NM_001085487	MYSM1	<0.002	1.89	Up

circRNA: Circle RNA; FC: Fold change; MYSM1: Myb-like SWIRM and MPN domain 1; NR3C1: Nuclear receptor subfamily 3 group C member 1.

Table 3

Primers for RT-qPCR.

Gene name	Forward primer	Reverse primer
β-actin	GAGCTACGAGCTGCCTGACG	GTAGTTTCGTGGATGCCACAG
hsa_circRNA_001544	CAGAACAGCAACATTTGAAGGGCC	ACTGGTCGACCTATTGAGGTTTGC
hsa_circRNA_102532	CCGAAGAGGGAGAGACGGAAG	CCTGCTGGTCTTTCAGGCCTAAC
hsa_circRNA_008961	TGCCCTGCTGACAAGAAGAGATG	AGCAACTTCCACTTCCACTTCCATC
hsa_circRNA_012732	TGAGCATCTCTCTGCTGAGGAGTTG	TCACTTGTAAAACAGTGCGGCTTCC

circRNA: Circular RNA; RT-qPCR: Real-time fluorescence quantitative polymerase chain reaction.

Statistical analysis

Statistical analysis and graphing were performed using statistical software SPSS, version 24.0 (IBM Corp, Armonk, NY, USA), and GraphPad Prism, version 8.0. Quantitative data with approximate normal distribution was expressed as mean ± standard deviation (SD), and data with non-normal distribution were expressed as median ± interquartile range (IQR). Mann-Whitney U test was used to measure the difference between each group; Spearman correlation was used to analyze the correlation between variables, and the diagnostic efficacy of biomarkers was evaluated by a receiver operating characteristic (ROC) curve; P < 0.05 was taken as the statistical significance of the difference.

Results

More than 1300 circRNAs (675 circRNAs were upregulated and 693 were downregulated) (FC > 1.5, P < 0.05) differentially expressed between AS and HC groups were found among 12,690 circRNAs. The differences were determined by the volcano plot and hierarchical clustering. The source of differentially expressed circRNAs was analyzed, and the upregulated circRNAs were distributed in all chromosomes except chrY, while the downregulated circRNAs were distributed in all chromosomes; in addition, there was no statistically significant difference in the distribution of these differentially expressed circRNAs among chromosomes (χ2 = 34.73, P > 0.05), while there were statistically significant differences between the upregulated genes and downregulated genes on the 2nd, 4th, 8th, 17th, 19th, 21st, and 22nd chromosomes (χ2 = 8.18; χ2 = 26.81; χ2 = 6.07; χ2 = 7.11; χ2 = 17.02; χ2 = 5.35; χ2 = 5.33, P < 0.05). And most circRNAs were derived from exons. Next, the eight circRNAs screened out were verified by RT-qPCR in the 12 samples originally submitted for inspection, and they were found to be consistent with the microarray chip results [Figure 1].

Figure 1

Differential expression between AS and HC groups. (A) Volcano plot of differentially expressed circRNAs. (B) Hierarchical clustering of differentially expressed circRNAs. (C) Chromosome distribution of differentially expressed circRNAs. (D) Source distribution of differentially expressed circRNAs. (E) Relative expression levels of eight circRNAs in six AS patients and six HCs. (The Y-axis is the ratio of the relative expression level of circRNA in the AS group to the HC group.) AS: Ankylosing spondylitis; circRNA: Circular RNA; HC: Healthy controls.

GO analysis and KEGG pathway analysis results

The GO analyses of differentially expressed circRNAs were carried out from three aspects: molecular function, cell composition, and biological process. The results showed that the annotation of molecular functions in the upregulated circRNAs was mainly related to “enzyme binding” (GO: 0019899), the annotation of cell composition was mainly concentrated in “intracellular” (GO: 0005622), and the annotation of biological processes was mainly in “organelle organization” (GO: 0006996) and “cellular macromolecule metabolic process” (GO: 0044260); the annotations on the molecular functions of downregulated circRNAs were mainly enriched in “adenyl ribonucleotide binding” (GO: 0032559) and “adenyl nucleotide binding” (GO: 0030554), cell composition annotations were mainly concentrated in “intracellular” (GO: 0005622) and “cytoplasm” (GO: 0005737), and biological process annotations were mainly related to “cellular component organization or biogenesis” (GO: 0071840) and “cellular component organization” (GO: 0016043). Furthermore, the upregulated differentially expressed circRNAs in the KEGG pathway were mainly enriched in “MAPK signaling pathway” (hsa04010) and “human T cell leukemia virus 1 infection” (hsa05166), while downregulated circRNAs were mainly enriched in “endometrium cancer” (hsa05213) pathway [Figure 2].

Figure 2

Bioinformatics analysis results of differentially expressed circRNAs. (A) GO analysis results of upregulated circRNAs. (B) GO analysis results of downregulated circRNAs. (C) KEGG pathway analysis results of upregulated circRNAs. (D) KEGG pathway analysis results of downregulated circRNAs. circRNA: Circle RNA; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

RT-qPCR

The RT-qPCR results showed that, compared with the HC group, the expressions of hsa_circRNA_001544 (U = 487, P < 0.05) and hsa_circRNA_102532 (U = 645, P < 0.05) in the AS group was significantly increased, which was consistent with the microarray results. However, the other two circRNAs were not statistically significant between the AS group and the HC group. Next, we divided the AS group into two subgroups of the ASA group (n = 30) and the ASS group (n = 30) and then performed statistical analysis again in three groups. It was found that compared with the HC group, hsa_circRNA_001544 was significantly increased in both ASA (U = 214, P < 0.05) and ASS groups (U = 273, P < 0.05), while hsa_circRNA_008961 (U = 250, P < 0.05) and hsa_circRNA_102532 (U = 295, P < 0.05) were only significantly increased in the ASA group. Furthermore, hsa_circRNA_012732 was significantly different between the ASA and ASS groups (U = 194, P < 0.05), and there was no statistical significance among the remaining groups [Figure 3].

Figure 3

RT-qPCR results of circRNA. (A) Differentially expressed circRNA between the AS group and HC group (hsa_circRNA_001544, hsa_circRNA_102532, hsa_circRNA_008961, and hsa_circRNA_12732). (B) Differentially expressed circRNAs among the ASA group, ASS group, and HC group (hsa_circRNA_001544, hsa_circRNA_102532, hsa_circRNA_008961, and hsa_circRNA_12732). AS: Ankylosing spondylitis; ASA: Active AS; ASS: Stable AS; circRNA: Circular RNA; HC: Healthy controls; RT-qPCR: Real-time fluorescence quantitative polymerase chain reaction. ∗ means P<0.05, † means P<0.001, ns means no significance.

Correlation analysis of circRNAs and AS clinical indicators

Spearman correlation analysis showed that the expression level of hsa_circRNA_012732 was negatively correlated with the BASDAI (r = −0.284, P < 0.05), BASFI (r = −0.279, P < 0.05), hsCRP level (r = −0.334, P < 0.05), and GLOB (r = −0.431, P < 0.05) and positively correlated with lymphocyte count (LY) (r = 0.260, P < 0.05), mean corpusular volume (MCV) (r = 0.367, P < 0.05), and ALB (r = 0.307, P < 0.05) [Figure 4]. And hsa_circRNA_008961 was negatively correlated with platelet (PLT) count (r = −0.334, P < 0.05). Correspondingly, there was no statistical significance among the remaining indicators.

Figure 4

Correlation analysis of circRNA and AS-related indicators. The correlation analysis among hsa_circRNA_012732 and LY, MCV, ALB, BASDAI, BASFI, hsCRP level, and GLOB. The correlation analysis between hsa_circRNA_008961and PLT. ALB: Albumin; AS: Ankylosing spondylitis; BASDAI: Bath Ankylosing Spondylitis Disease Activity Index; BASFI: Bath Ankylosing Spondylitis Functional Index; circRNA: Circle RNA; GLOB: Globulin; hsCRP: High-sensitivity C-reactive protein; LY: Lymphocyte count; MCV: Mean corpusular volume; PLT: Platelet.

ROC curve analysis

ROC curve analysis showed that the expression levels of two circRNAs were statistically different between the AS group and the HC group (P < 0.05); the area under the ROC curve (AUC) (95% CI) values of hsa_circRNA_001544 and hsa_circRNA_102532 were 0.720 (0.609–0.831) and 0.642 (0.521–0.762), respectively [Figure 5].

Figure 5

AUC of hsa_circRNA_001544, hsa_circRNA_102532, hsa_circRNA_008961, and hsa_circRNA_012732. AUC: Area under curve; circRNA: Circle RNA.

Target miRNA analysis

A co-expression network of circRNA and miRNA was constructed, and the five miRNAs with the highest matching value were selected as the predicted target miRNAs of circRNA. Table 4 shows the predicted miRNAs of hsa_circRNA_001544, hsa_circRNA_012732, hsa_circRNA_102532, and hsa_circRNA_008961.

Table 4

The target miRNAs of six circRNAs.

CircRNA	Target miRNAs
hsa_circRNA_001544	has-miR-3681-5p, hsa-miR-3691-5p, hsa-miR-450b-5p, hsa-miR-412-5p, hsa-miR-6855 -5p
hsa_circRNA_012732	hsa-miR-4659a-3p, hsa-miR-153-5p, hsa-miR-3680-3p, hsa-miR-4778-3p, hsa-miR-4659b-3p
hsa_circRNA_102532	hsa-miR-658, hsa-miR-215-3p, hsa-miR-432-5p, hsa-miR-92a-1-5p, hsa-miR-144-5p
hsa_circRNA_008961	hsa-miR-6881-3p, hsa-miR-6783-3p, hsa-miR-6857-5p, hsa-miR-450a-2-3p, hsa-miR-1248

circRNA: Circular RNA; miRNA: microRNA.

Discussion

The pathogenesis of AS is still unclear, but studies have found that genetic factors account for >90%, and HLA-B27 plays a very critical role in it.[ In an epidemiological survey abroad, it was found that the positive rate of HLA-B27 in white patients can reach 90%, but it can also have a positive rate of 6% in the normal population;[ therefore, a positive HLA-B27 is not enough to diagnose AS. Generally, the evidence of imaging changes is required for diagnosis of AS in clinics. However, when the patient's imaging changes are clearly observed, the disease usually has developed into irreversible damage, and the patient may miss the optimal intervention time, thereby affecting their prognosis. Under such conditions, there is an urgent need for new and more specific indicators for the diagnosis of AS. Similarly, the need for disease activity evaluation and treatment targets is equally important.

CircRNA is one of the ncRNA, which is thought to encode almost no proteins.[ Although ncRNA plays a huge role in the process of life, participating in biological growth,[ differentiation,[ metabolism,[ disease occurrence,[ and other processes, its specific mechanism is still an unsolved mystery. Studies have shown that circRNA can participate in the occurrence of many autoimmune diseases and can regulate the course of diseases by regulating the body's immune response[ and immune tolerance.[ In our experiments, we clearly discovered the differential expression of many circRNAs in AS, further revealing that circRNA plays an important role in autoimmune diseases.

In the results of microarray analysis, there was no statistical difference in the distribution of circRNA parental genes among the chromosomes, but there were statistically significant differences between the upregulated genes and downregulated genes on the 2nd, 4th, 8th, 17th, 19th, 21st, and 22nd chromosomes, which suggested that the genes distributed on these chromosomes might have different degrees of mutations in the pathogenesis of AS, and then participated in the regulation of the pathogenesis of AS.

By analyzing the results of bioinformatics analysis, we can predict the potential mechanism of circRNA involved in the occurrence and development of AS to a certain extent. The results of previous microarray analysis suggested that most of the differentially expressed circRNAs were derived from exons, while exon-derived circRNAs are often obtained through lariat-driven circularization and back-splicing and are often located in the cytoplasm.[ This coincides with the finding that the differentially expressed circRNAs were mainly located in the cytoplasm in the GO analysis. In the KEGG pathway analysis, the upregulated genes were mostly enriched in different inflammation or infection-related pathways, including “MAPK signaling pathway,” “TNF signaling pathway,” “Yersinia infection,” “human immunodeficiency virus 1 infection,” etc., which suggested that inflammation activation was closely related to the occurrence and development of AS. The downregulated genes were mostly enriched in tumor-related pathways, including endometrial cancer, gastric cancer, colon cancer, breast cancer, and so on. A study from the United States reviewed the status of solid cancer in elderly patients with AS from 1999 to 2010 and found that the risk of esophageal cancer, gastric cancer, and lung cancer was lower and the risk of melanoma, urinary system cancer, and breast cancer was higher.[ Two other studies have similar findings.[ Combined with the enrichment results of this study, there may be similar pathways but different feedback mechanisms between AS and different tumors, especially colon cancer, gastric cancer, and other gastrointestinal cancers. This may have important implications for the treatment and prognosis of the disease.

In the pathogenesis of AS, inflammation and immune responses are both prominent manifestations, which are fully reflected in the experimental results of this time. We analyzed the RT-qPCR of each group and found that hsa_circRNA_001544 and hsa_circRNA_102532 were significantly different between the AS group and HC group, while there was no statistical difference between the ASA group and ASS group. This indicated that hsa_circRNA_001544 and hsa_circRNA_102532 may not be related to the inflammatory response of AS, and more importantly, they were related to the disease itself. The parent gene of hsa_circRNA_001544 is nuclear receptor subfamily 3 group C member 1 (NR3C1), which encodes the glucocorticoid receptor. At present, studies have found that it may be related to the pathogenesis of rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE), but it may affect the role of hormones in the treatment of diseases.[ Therefore, hsa_circRNA_001544 whether it is related to the expression of glucocorticoid receptors, there is still much to be explored. At the same time, hsa_circRNA_012732 was significantly lower in the ASA group than in the ASS group, which suggested that hsa_circRNA_012732 may be related to inflammation and it may participate in the disease process by participating in the inflammation regulation of AS. In addition, we also found that hsa_circRNA_008961 was significantly higher in the ASA group than in the HC group, and it also suggested that there might be a correlation between hsa_circRNA_008961 and the inflammatory of AS.

We further analyzed the correlation between four circRNAs and the clinical indicators of AS, and hsa_circRNA_012732 was found to be negatively correlated with the BASDAI, BASFI, hsCRP level, and GLOB but positively correlated with LY, MCV, and ALB. Both the BASDAI and hsCRP level in clinical work can be used to reflect the inflammation level of AS and assess whether the patient is in the active stage of disease. GLOB, LY, and ALB play an important role in the pathogenesis of autoimmune diseases.[ BASFI is often used to assess the functional status of the patient's body, which can reflect the prognosis of the patient to a certain extent.[ Correlation analysis demonstrated that hsa_circRNA_012732 was downregulated during inflammation, but the specific mechanism is unclear. Myb-like SWIRM and MPN domain 1 (MYSM1) is the parent gene of hsa_circRNA_012732. Studies have found that MYSM1 was a key inhibitor of innate immunity and autoimmunity. It can inhibit cell proliferation, plasma cell differentiation, and antibody production and can be used as a potential therapeutic agent for inflammatory and autoimmune diseases.[ This suggested that hsa_circRNA_012732 and MYSM1 may have opposite effects on the mechanism of action. The hsa_circRNA_012732 may act as a pro-inflammatory factor in AS, participate in AS disease progression through immune-mediated inflammation, and has the potential to become the indicator of disease activity and prognosis, even become the target of AS treatment. In the correlation analysis, we also found that hsa_circRNA_008961 was negatively correlated with PLT. An elevated PLT count was found in a meta-analysis to assess the disease activity of AS.[ Therefore, hsa_circRNA_008961 may also have the ability to assess the disease activity of AS.

More and more studies have shown that circRNA can participate in the regulation of disease by acting as the miRNA sponge.[ Hsa_circRNA_102532 target miRNAs were hsa-miR-92a-1-5p, miR-144-5p, etc. Experimental verification has shown that miR-92a-1-5p can negatively regulate β-catenin, a positive regulator of osteogenic differentiation.[ There were also studies that the overexpression of miR-144-5p might reduce the viability of THP-1 macrophages by inhibiting the expression of Toll-like receptor 2 (TLR2) and inhibiting the activation of nuclear factor kappa-B (NK-κB) signaling and inhibit the expression of tumor necrosis factor (TNF)-α, interleukin (IL)-6, and IL-8 in the cells. It was also proposed that miR-144-5p had the potential as a new therapeutic target for RA.[ Therefore, hsa_circRNA_102532 can participate in the disease development of AS through many aspects in the process of acting as the miRNA sponge. Similarly, the above suggests that the mutual regulation between circRNA and miRNA may often be a multidirectional network rather than a single axis.

The ROC curve can reflect the sensitivity and specificity of the disease diagnosis index to a certain extent. It is drawn by setting the sensitivity% as the ordinate and 1-specificity% as the abscissa, using AUC indicates the diagnostic efficiency of the diagnostic index. The range of AUC is 0.5–1. The value closer to one indicates the better the diagnostic efficiency. When AUC < 0.5, it means the index has no diagnostic power; when 0.5 < AUC < 0.7, the diagnostic power is low; when 0.7 < AUC < 0.9, it means it has a certain diagnostic power; when AUC > 0.9, the index has higher diagnostic power. We performed ROC curve analysis on these four circRNAs in AS. Among them, the AUC of hsa_circRNA_001544 was >0.7. The AUC (95% CI) of hsa_circRNA_001544 was 0.729 (0.622–0.837). Therefore, hsa_circRNA_001544 is expected to become a molecular marker for the diagnosis of AS, and hsa_circRNA_102532 also has low diagnostic efficiency for AS, which could provide a new method for the diagnosis of AS.

Conclusion

In this study, we found that the differential expression of circRNAs in PBMCs of AS patients might be related to the occurrence and development of AS. Among these differentially expressed circRNAs, hsa_circRNA_102532 and hsa_circRNA_001544 have the potential to become molecular markers for the diagnosis of AS, while hsa_circRNA_012732 has the potential to become an indicator of disease activity and even a therapeutic target.

Funding

This work was supported by grants from the National Natural Science Foundation of China (No. 81974250) and Nanchong City Science and Technology Project (No. 18SXHZ0522).

Conflicts of Interest

None.