Screening for differentially expressed circRNA between Kashin–Beck disease and osteoarthritis patients based on circRNA chips CURRENT STATUS: POSTED

Objective: The purpose of this research was to explore differentially expressed circRNA between OA and KBD and the potential differential diagnosis biomarkers. Methods: Total RNA was extracted from 5 pairs of KBD and OA knee joint cartilage specimens, and the expression of circRNAs was analyzed by Chip Scanning Analysis. The microarray data was verified by quantitative polymerase chain reaction (qRT-PCR) analysis. GO enrichment analysis and KEGG pathway were used to predict the functions of the differentially expressed circRNAs. A circRNA-miRNA network was constructed to predict targeting microRNAs of differentially expressed circRNA genes. On the basis of microarray, we expanded the sample size, peripheral blood samples from 25 KBD patients and 25 OA patients (five pairs of patients for chondrocyte microarray study included) were collected and qRT-PCR was performed for hsa_circRNA_0020014 verification. Diagnostic value was evaluated by the area under the receiver operator characteristic (ROC) curve. Results: A total of 1627 circRNAs were differentially expressed between OA and KBD (P<0.05; 0.5<fold change>2), 1328 were up-regulated and 299 were down-regulated among them. Five differentially expressed genes associated with bone and joint disease were chosen for further qRT-PCR validation. After obtaining the parental genes of them, functional annotations were performed on the top ten enrichment GO items and KEGG pathways. The difference in expression profile of hsa_circRNA_0020014 was confirmed by qRT-PCR at articular cartilage level, and its circRNA-miRNA regulation network was set up. The ROC curve demonstrated that hsa_circ_0020014_CBC1 in peripheral blood could distinguish patients with KBD and OA (area under the curve: 0.6415, P< 0.01). The hsa_circ_0020014_CBC1 may be a potential biomarker for differential diagnosis between KBD and OA patients Conclusion: Our results suggested that the expression profiles of circRNA were significantly different between OA and KBD. hsa_circRNA_0020014 is a potential biomarker In our study, we found a large number of differentially expressed circRNAs in OA and KBD cartilage specimens. The microarray analysis between KBD and osteoarthritis chondrocytes was conducted to analyse the differentially expressed circRNAs. A total of 1627 circRNAs were differentially expressed, include 1328 upregulated circRNAs and 299 downregulated circRNAs. We use qRT-PCR to detect the expression levels of five selected circRNAs, the results were consistent with the results of the analysis of the microarray analysis data, demonstrating the reliability of the spectrogram. GO and KEGG pathway analyses were used for functional annotation of differentially expressed circRNA genes. We found that these genes are involved in RNA transport, RNA degradation, and the p53 pathway. The most enriched and meaningful disease pathways were related to “vascular diseases, short rib −polydactyly syndrome, distal hereditary motor neuropathies” (from KEGG-disease database). “generalized atrophic benign, chondrodysplasia with joint dislocations (GPAPP type), and Parkinson’s disease” (from OMIM database).

found in plants and confirmed to encode a sub-virus [1], they are characterized by being singlestranded, covalent, and closed loops [2]. Later, circRNA was found to be highly expressed in mammalian brain tissues [3]. It has been demonstrated that the circRNA transcripts of sexdetermining region Y (SRY) and circular RNA CDR1as [4] can be used as miRNA sponges [5]. Studies have shown that circRNA is widely involved in gene expression regulation, and has been proved to play an important role in the occurrence and development of nervous system diseases [6]. Due to their high degree of conservation [7], stablility [8] and accessible characteristics, the circRNAs could serve as biomarkers for diagnostics of clinical disease, such as gastric cancer [9], hepatocellular carcinoma [10] and oesophageal cancer [11]. The hsa_circ_0005836 was found as a biomarker and therapeutic target of active pulmonary tuberculosis [12].
Osteoarthritis (OA) is a joint disease characterized by cartilage degeneration and secondary subchondral bone thickening, osteophyte formation and synovitis, which eventually leads to joint stiffness, pain, swelling and functional loss [13,14]. Due to tissue wear and chondrocyte death, cytokines and growth factors lead to increased catabolism and abnormal differentiation in OA, leading to degradation of extracellular matrix (ECM) [15].
Kashin-Beck disease (KBD) is a unique osteoarthropathy in China, diagnosed based on the national KBD diagnostic criteria (WS/T 207-2010). The disease is mainly distributed in the low-selenium areas of the soil from the northeast to the southwest of China. At present, there are more than 2.5 million current KBD patients and about 30 million residents at high risk [16]. Similar to OA, such as chondrocyte degeneration, ECM degradation and other basic pathological features were displayed [17,18], although specific differences between OA and KBD chondrocyte gene expressions have been revealed [19]. In contrast to OA, which commonly appears in middle-aged and elderly people, the onset of KBD occurs mostly in 3-12 years old children. The main clinical manifestations of KBD were enlarged finger joints, short fingers, curvature of distal fingers and short stature [17]. Until now, the diagnosis of KBD has been based on clinical features and typical imaging characteristics of metaphyseal bone, typical manifestations are irregular sclerosis and the appearance of calcification.
KBD may be classified as early, stage І, II and III, based on the "Diagnosis Criteria of KBD in China (WS/T207-2010). The above methods are not sensitive enough to the early stage diagnosis of the disease, and it is difficult to intervene early in the late performance such as short stature and joint deformation.
At present, the generally accepted pathogenic theories of environmental factors concerning the pathogenesis of KBD include environmental selenium deficiency, high levels of humic acid in drinking water and contamination of grains by mycotoxins [20]. The exact pathologic mechanism of KBD is not yet fully understood. Because of the similar clinical manifestations between KBD and OA, KBD is known as an endemic OA. In recent years, circRNA has been proved to be a new biomarker in human peripheral blood [21] for diagnosis of diseases, such as coronary artery disease, OA and rheumatoid arthritis (RA) [22][23][24]. Until now, there have been no identification and characterization of circRNAs for OA and KBD in bone or cartilage tissue as well as in peripheral blood. Therefore, we focused on articular chondrocytes of patients with KBD and OA, to explore the genetic features that could be associated with early diagnosis.  [25]. Patients with genetic bone disease and RA were excluded. In this research, a total of 10 participants were selected and matched according to gender and age and were divided into two groups (5 patients with KBD and 5 patients with OA). Samples were selected randomly from outpatients and inpatients at The First Affiliated Hospital of Xi'an Jiaotong University, and from the KBD-endemic area Linyou County (Baoji, China) from April 2016 to October 2017. The study was approved by the Human Ethics Committee of Xi'an Jiaotong University, project approval number , and all subjects signed informed consent. KBD and OA chondrocyte samples were taken from the same anatomical region of the femoral condyle of the knee joints. The chondrocyte specimens were obtained and cut part rapidly, and then frozen in liquid nitrogen and immediately stored at -80°C until RNA extraction.

Collection of peripheral blood samples
Total of 3 ml peripheral blood was drawn from each donor before breakfast by peripheral venipuncture and gathered in ethylene diamine tetra acetic acid (EDTA) anticoagulant vacutainers, and kept in the -80℃ degree refrigerator until used for RNA isolation. Then, we expanded the sample size, total RNA was extracted from the peripheral blood samples of twenty-five patients with OA and twenty-five patients with KBD (Five pairs of patients for chondrocyte microarray study included).

CircRNA microarray analysis
The circRNA microarray analysis was summarized, normalized and quality control using GeneSpring software v13.0 (Agilent Technologies). The differentially expressed circRNAs in OA and KBD chondrocytes were identified through fold change, hierarchical clustering was used to display differentially expressed circRNAs between them. The limit values are as follows: fold change for upregulated circRNAs was ≥ 2 and for down-regulated ones ≤ 0.5 and P < 0.05.

qRT-PCR validation of microarray data
While writing this paper, we have consulted a lot of relevant literature. Matta C et al [26] published a review article in 2014 showing that the calcium signaling pathway may be related to OA, while the target gene PDE1C of hsa_circ_0134111_CBC1 participates in the same pathway. Wang L et al [27] reported in 2015 that the PRKCH gene is associated with RA, systemic lupus erythematosus, ankylosing spondylitis and OA, our research found that the target gene of hsa_circ_0032131_CBC1 is PRKCH. Sun HY et al reported in 2017 that the MAPK signaling pathway is related to OA, and the target gene dusp5 of hsa-circ-0020014-cbc1 is involved in the MAPK signaling pathway [28]. In 2018, Hu J et al [29] found that the development of microRNA-17-5p was closely related to the development of OA, and the prediction of microRNA-17-5p could be combined with hsa_circ_0094742_CBC1. Fang Y et al reported in 2019 [30] that OA may be associated with protein digestion and absorption, while col3a1, the target gene of hsa-circ-0057421-CBC1, is involved in protein digestion and absorption.
The above literature is the theoretical basis for us to select the five differentially expressed genes for further qRT-PCR verification. The selected circRNAs were analyzed using the ABI PRISM 7500 Sequence Detection System (Applied Biosystems, Life Technologies, Waltham, MA, USA). The primers used in our study are listed in Table 1.

qRT-PCR verification of related circRNAs in peripheral blood samples
To further validate the results of microarray analysis, the significant differentially expressed gene was chosen as target genes for qRTPCR. Total RNA in peripheral blood samples were reverse transcribed into cDNA using superscript II reverse transcriptase (invitrogen; Thermo Fisher Scientific, Inc.

Data analysis
First, raw data extraction from scanned images by Genepix Pro 6.0 software (Axon; Molecular Devices, LLC, Sunnyvale, CA, USA). Next, the original data is normalized and processed using the R software package. Finally, exclude low-intensity substances. The fold change of each cirRNA was calculated and the difference between OA and KBD patients was compared. The limit values are as follows: P < 0.05and fold change≥ 2 and ≤ 0.5, the differentially expressed circRNAs were arranged according to the P-value and fold change.      and "mRNA capping", "phosphatidylinositol signaling system", "RNA degradation" and "RNA transport" (from the KEGG pathway database). The most enriched and meaningful pathways from the PANTHER database are related to the "PDGF signaling pathway", "inflammation mediated by chemokine and cytokine signaling pathway", "endothelin signaling pathway", and the "p53 pathway".

Related target miRNAs and the miRNA/circRNA interaction network
Considering microarray analysis and verification results, hsa_circRNA_0020014 was expressed significantly differentially in chondrocytes between OA patients and KBD patients and closely related to bone and joint diseases. To explore target miRNAs that may be associated with OA and KBD, the validated hsa_circRNA_0020014 was selected for further analysis and the miRNA/circRNA interaction network was predicted by miRNA targeting prediction software. We ranked the top 10 from low to high according to the tot energy (table 3), and the circRNA/microRNA regulatory network diagram of hsa_circ_0020014_CBC1 was drawn on this basis (Figure 7).

Figure 7. CircRNA-miRNA regulation network constructed and visualized by Target Scan and miRanda.
Network diagram shows the interaction based on the predicted target genes of hsa_circRNA_0020014 targeting miRNAs.

Discussion
Protein-coding genes and their transcription in eukaryotic cells have been extensively studied [32].
The reality is that a large number of sequences in the human genome do not encode proteins, and non-coding RNA accounts for almost 95% of the total RNA transcribed from the eukaryotic genome, including circular RNA [33]. There are two main types of RNA molecules: endogenous RNA molecules formed by nonlinear reverse splicing and intron-derived cyclic RNA molecules [33]. In the 1970s, circular RNA was first discovered in RNA viruses [34,35]. CircRNA was also found in the transcripts of human cells in 1993, but at that time, it was considered to be a low-abundance RNA molecule formed by mis-splicing of exons during transcription [36]. With the rapid development of RNA sequencing and bioinformatics, circRNAs are not a rare phenomenon, but present widely in the cytoplasm of eukaryotic cells and play an important role in regulating gene expression at the post-transcriptional level. Studies have shown that circRNAs acts as a microRNA sponge that competitively inhibits miRNA activity and thereby regulating target genes of miRNA. Research has shown that circRNAs are involved in the initiation and development of many kinds of diseases. For example, due to the deficiency in the "sponge" effect of CIRS-7, miRNA-7 is up-regulated, lead to down-regulating Alzheimer's disease-relevant target genes [6]. Some scholars have reported that circRNAs play an important role in the progression of cardiovascular diseases, such as hsa_circ_0124644, which may be associated with coronary artery disease [22]. And many studies have reported that circRNAs are closely related to a variety of tumors. In one study, Weihai Liu et al [37] noted that hsa_circRNA_103801 and hsa_circRNA_104980 can be used as biomarkers for the diagnosis and treatment of osteosarcoma. Research on circRNAs relating to osteoarticular diseases mainly includes the following results. Some studies have systematically explored the potential functions of non-coding RNAs in osteoclast genesis by microarray technology, and constructed co-expression networks of lncRNA-mRNA and circRNA-miRNA [38].
In another study, circRNA was shown to potentially affect the development and progression of RA and play an important role in the regulation of related gene expression (24). Liu Qiang et al explored the function of human cartilage degeneration mechanical stress-related circRNA (circRNAs-MSR) and concluded that inhibition of circRNAs-MSR can limit the degradation of the ECM of chondrocytes, therefore, circRNAs-MSR may be a potential target for the treatment of osteoarthritis [39]. However, as a prevalent degenerative joint disease, the role of circRNAs in the pathogenesis of OA is not very clearly. A recent study indicated that hsa-circ-0045714 inhibits the progression of OA by promoting aggregation of type II collagen and proteoglycan expression and chondrocyte proliferation [40]. For the special endemic OA (KBD) in China, there has been no report related to circRNA until now.
In our study, we found a large number of differentially expressed circRNAs in OA and KBD cartilage specimens. The microarray analysis between KBD and osteoarthritis chondrocytes was conducted to analyse the differentially expressed circRNAs. A total of 1627 circRNAs were differentially expressed, include 1328 upregulated circRNAs and 299 downregulated circRNAs. We use qRT-PCR to detect the expression levels of five selected circRNAs, the results were consistent with the results of the analysis of the microarray analysis data, demonstrating the reliability of the spectrogram. GO and KEGG pathway analyses were used for functional annotation of differentially expressed circRNA genes. We found that these genes are involved in RNA transport, RNA degradation, and the p53 pathway. The most enriched and meaningful disease pathways were related to "vascular diseases, short rib −polydactyly syndrome, distal hereditary motor neuropathies" (from KEGG-disease database).
In addition, the most significantly differentially expressed and confirmed circRNA (hsa_circRNA_0020014) was further explored and revealed the potential molecular mechanisms of regulation by circRNAs in joint disease pathogenesis. we expanded the number of samples of the subjects, 25 KBD patients and 25 OA patients were selected and divided into two groups according to the age and gender matching principle, peripheral blood samples were collected and RNAs were extracted for qRT-PCR verification. Peripheral blood verification results were consistent with chondrocyte verification results, and hsa_circ_0020014_CBC1 was significantly different between OA and KBD. Moreover, ROC curve was used to evaluate the diagnostic value and reliability, the levels of hsa_circ_0020014_CBC1 in peripheral blood appeared potential diagnostic value between KBD and OA, and a high ROC AUC value (0.6415, P< 0.01) compared with traditional bone metabolism indicators, it has great potential as a diagnostic biomarker. Furthermore, research on the function of hsa_circ_0020014_CBC1 can improve the understanding of the mechanism of these two diseases occurrence and development.
According to the sponge theory, we predicted the potential target miRNAs and constructed a circRNA-miRNA interaction regulation network. The predicted microRNAs were selected according to tot energy using Cytoscape software, and the smaller the energy was, the tighter interaction between circRNAs and miRNAs. Top ten predicted miRNAs of hsa_circRNA_0020014 were identified in the networks. Therefore, they are theoretically related to osteoarthrosis and require further study. Two of the top 10 predicted miRNAs have been reported in a recent study, hsa-circ-0061012 was upregulated in psoriasis lesions compared to normal healthy skin tissue [41]. Some scholars reported that miR-6743-5p acts as an oncomiRNA in glioma by targeting GRIM-19 and STAT3 (39). And miR-6743-5p expression is downregulated in patients with lymph node metastasis in oesophageal cancer [42]. There is no definitive report on OA and KBD in these predicted targeted microRNAs. The gene symbol corresponding to hsa_circRNA_0020014 is dual-specificity phosphatase 5 (DUSP-5), studies have shown that in in DUSP5-overexpressing mice, the severity of arthritis, as shown by clinical arthritis scores and tissue inflammation and cartilage damage was attenuated [43]. Pramanik noted that mutations in dusp-5 and snrk-1 have been found in tissues of patients with vascular disorder [44]. Wickramasekera NT and his team have drawn conclusions in a study showing that DUSP-5 appears to regulate the signalling of pressure-dependent myogenic cerebral arterial constriction, it is essential for maintaining the blood supply to the brain [45].

Conclusion
We used high-throughput information analysis to explore the potential biological function of circRNAs in human primary OA and endemic OA (KBD) in China for the first time. In compliance with WHO requirements for minimally invasive biomarkers with clinical significance we enlarged the sample size and extracted RNA from peripheral blood for further validation of previous research results. From this, we can draw the following conclusion: hsa_circ_0020014_CBC1 can be detected from peripheral blood at a relatively low cost, fast, and minimally invasive. Therefore, hsa_circ_0020014_CBC1 can be used as a biomarker for differential diagnosis between OA and KBD.      Classification of differentially expressed genes according to genome origin.

Figure 3
The consistency of qRT-PCR and chip analysis of 5 differentially expressed circRNA molecules. From left to right represent Hsa_circ_0020014_CBC1, hsa_circ_0032131_CBC1, hsa_circ_0057421_CBC1, hsa_circ_0094742_CBC1, hsa_circ_0134111_CBC1.    CircRNA-miRNA regulation network constructed and visualized by Target Scan and miRanda.
Network diagram shows the interaction based on the predicted target genes of hsa_circRNA_0020014 targeting miRNAs. Table 3: Top 10 most likely miRNA targets of hsa_circRNA_0020014 predicted and selected based on tot energy.

Figure 8
ROC curve of hsa_circ_0020014_CBC1, hsa_circ_0032131_CBC1 in peripheral blood of KBD patients and OA patients (two ROC graphs will be combined into one picture later).