Integrated Analysis of a circRNA-miRNA-mRNA Axis in Intervertebral Disc Degeneration


 Background: Low back pain (LBP) is a common symptom in daily life and one of the primary causes is intervertebral disc degeneration (IDD). Growing studies have indicated that circular RNAs (circRNAs) are intimately associated with IDD; however, the underlying mechanism has not yet been elucidated. We aimed to explore how circRNAs regulate IDD in an effort to provide novel insight for clinical diagnosis and treatment. Methods: The sequencing data of circRNAs, microRNAs (miRNAs), and mRNA were acquired from Gene Expression Omnibus (GEO) datasets. By analyzing the dataset consisting of a control group and degenerated group, differentially expressed circRNAs, miRNAs, and mRNAs were collected, and then the intersection of circRNAs, miRNAs, and mRNAs was screened. According to these intersectional RNAs, we constructed an integrally circRNA-miRNA-mRNA network. Finally, using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, we further clarified functions of the intersectional mRNA in IDD. Results: we obtained 620 differentially expressed circRNAs(DEcircRNAs), 13 miRNA (DEmiRNA), 273 mRNAs(DEmRNAs), 12 intersectional miRNAs, and 47 intersectional mRNAs. Finally, based on interactional 8 circRNA, 5 miRNAs and 15 mRNAs, an integrally circRNA-miRNA-mRNA network was constructed. Eight circRNAs, contained hsa_circ_0032254, hsa_circ_0003183, hsa_circ_0032253, hsa_circ_0001293, hsa_circ_0004565, hsa_circ_0091570, hsa_circ_0077526, and hsa_circ_0057552, may regulate IDD onset and progression by acting as competing endogenous RNAs. The results of GO and KEGG analyses implied that the targeted genes might significantly correlate to IDD.Conclusion: our findings improved a better understanding of the circRNA-related ceRNA regulatory mechanism in IDD and offered possible targets for IDD treatment.


Introduction
Low back pain is the primary cause of disability in developed countries and a growing number of people are affected worldwide [1]. Although the reason for low back pain remains unclear, some research has demonstrated that low back pain is correlated with intervertebral disc degeneration (IDD) [2][3][4][5][6][7].
Anatomically, the intervertebral disc consists of three components: I) an inner cord termed the nucleus pulposus (NP), II) annulus brous surrounding the NP, and III) the cartilaginous endplate [8,9]. Nucleus pulposus cells (NPCs) in the NP generate several components of the extracellular matrix (ECM) including collagen 2 and proteoglycans in physiological conditions [9,10]. These proteins are vital in maintaining disc structure and resistance to various mechanical compression forms from external factors [11]. However, multiple abnormal stimuli, such as neutrophil proteases, can upregulate the expression of in ammatory cytokines (e.g., IL-1β, TNF-α), and then breakdown the balance of NPCs through several factors including matrix metalloproteinases (MMPs), collagen II, metalloproteinases with thrombospondin motifs (ADAMTS), and aggrecan. These imbalances trigger or accelerate the degeneration of the intervertebral disc [12,13]. However, the deeper pathological mechanism underlying IDD is still unknown. Thus, it is essential to seek a more effective way to moderate the in ammatory response and reverse the NP microenvironment imbalance. Circular RNAs (circRNAs) are a type of single-stranded loop RNA that are primarily considered to act as competing endogenous RNAs (ceRNAs), binding microRNAs (miRNAs) to further regulate the expressed level of targeted mRNAs [14][15][16][17]. For instance, upregulation of circ-4099 has been shown to inhibit IDD development by modulating the circ-4099-miR-616-5p-Sox9 network [18]. An increasing number of reviews demonstrated that circRNAs have a crucial role in the pathogenesis of IDD [19,20]. In order to better know the underlying mechanism of circRNAs in IDD, we generated a circRNA-miRNA-mRNA network by sequencing data from the GEO database. This study may provide a new target in the pathogenesis and therapeutic strategies of IDD.

Data acquisition
Three expression pro les of circRNAs (GSE67566), miRNAs (GSE63492), and mRNAs (GSE124272) were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). And then we analyzed these data by limma R package. Gene symbols corresponding to the probes were obtained using the annotation information in the platform. Identi cation of the differential expression of circRNAs, miRNAs and mRNAs The differentially expressed pro les of circRNAs, miRNAs, and mRNAs between normal and degenerated NP tissues were acquired using the limma R package. It is necessary and statistically signi cant thresholds that Log fold change (FC) > 1.0 and adjusted P-value < 0.05 were needed.

Construction of the ceRNA network
The circBase website (http://www.circbase.org/) was used to obtain information regarding differentially expressed circRNAs(DEcircRNAs), including the chromosomal position of each DEcircRNA and precise localization of DEcircRNAs at the gene sequence level. Then, we determined the miRNAs (MREs) that combined with DEcircRNAs by the cancer-speci c circRNA database (CSCD) (https://gb.whu.edu.cn/CSCD/). We rst overlapped MREs and differentially expressed miRNAs (DEmiRNAs) for potential targeted miRNAs; then, the potential target mRNA was obtained using the TargetScan and miRDB databases. At last, the targeted genes were identi ed by the interaction of potential target mRNAs and DEmRNAs. Therefore, the regulated circRNA-miRNA-mRNA network was constructed by screening out intersectional circRNAs, miRNAs, and mRNAs. The regulatory network is presented using Cytoscape 3.7.1.

Functional and pathway enrichment analyses form databases of GO and KEGG
To ascertain the function of targeted genes in intervertebral disc degeneration, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) analyses were performed using the online tool Database of DAVID (http://david.abcc.ncifcrf.gov, version 6.8). An adjusted P-value < 0.05 was regarded as statistically signi cant. Table 1 listed the expressional basic information of the circRNAs, miRNAs and mRNAs in normal and degenerated NP tissues, which come from three sequencing data. Based on Log [FC] > 1.0 and adjusted Pvalue < 0.05, we identi ed 620 DEcircRNAs including 337 upregulated and 283 downregulated DEcircRNAs (Fig. 1). The basic information of the 620 DEcircRNAs was downloaded from circBase (http://www.circbase.org/cgi-bin/listsearch.cgi). Using the CSCD database, we found a total of 2,390 miRNAs combined with DEcircRNAs. In the differentially expressed dataset of miRNAs, eight DEmiRNAs were upregulated and ve DEmiRNAs were downregulated (Fig. 2). Twelve overlapped miRNAs were collected using the limma R package (Fig. 3). Using the TargetScan and miRDB databases, a total of 5,203 targeted miRNAs were predicted. We nd 182 mRNAs upregulated and 91 mRNAs downregulated in the IDD groups (Fig. 4). Finally, 47 intersected genes were identi ed (Fig. 5). Table 1 Basic information of the 3 microarray datasets from GEO Data source Platform Series Sample size (T/N)

Differential expression analysis
Construction of the ceRNA network To better understand the effect of circRNAs on IDD, we constructed an integrated circRNA-miRNA-mRNA network by linking circRNA-miRNA and miRNA-mRNA. The network contained eight circRNAs (two upregulated and six downregulated), ve miRNAs (two upregulated and three downregulated), and 15 mRNAs (eight upregulated and seven downregulated), which are shown using Cytoscape 3.7.1 (Fig. 6).

Functional and pathway enrichment analysis
The mainly enriched genes based on GO analysis are shown. The rst ve gathered terms were male sex determination, positive regulation of myelination, sex determination, regulation of myelination, positive regulation of nervous system process in the biological process (BP) category (Fig. 7a); speci c granule, tertiary granule, host cell cytoplasm part, host cell cytoplasm, and host intracellular part in the cellular component (CC) category (Fig. 7b); and nuclear receptor activity, ligand-activated transcription factor activity, steroid hormone receptor activity, RAGE receptor binding, and S100 protein binding in the molecular function (MF) category (Fig. 7c). The p53 signaling pathway and cortisol synthesis and secretion were strongly enriched in KEGG pathway analysis (Fig. 7d).

Discussion
The circRNAs were originally discovered in plant viroid and yeast mitochondrial RNAs [7,21] and have been rarely tested by transcriptomic polyadenylated RNA pro ling owing to the lack of polyadenylation at their 3' ends [22][23][24][25]. With the development of novel technology including calculation procedure and high-throughput sequencing [22], non-polyadenylated human RNA sequencing analysis and the signals of RNA sequencing can be obtained to assess back-spliced exons [26] or stabilized introns [27]. With their covalently closed structure, circRNAs are highly stable non-coding RNAs that play a critical role in cell and tissue development and several disease processes, including cancer, IDD, and neurodegeneration [19,[28][29][30].
Although the molecular mechanisms of circRNAs in mammalian cells remain mostly unclear, circRNA regulation of the level of mRNAs by a circRNA-miRNA-mRNA network may be an important mechanism in the initiation and development of IDD. Wang et al. showed that the circular RNA SEMA4B could weaken the impact of IL-1β on NPC senescence, ECM and aggrecan degradation in IDD via inhibition of miR-431 (a targeted miRNA of circRNA SEMA4B) [20]. Moreover, circRNA_104670 suppressed the expression level of collagen II and the proliferation of NPCs, and promoted NPCs' apoptosis by circRNA_104670/miR-17-3p/MMP-2 [19]. The circRNA VMA21, regarded as a sponge of mir-200c, plays a role in regulating the level of in ammatory cytokines by targeting mir200c and XIAP [31]. To research the role of circRNAs in IDD, we identi ed 620 DEcircRNAs by comparing normal and degenerated groups' expression pro les. Similarly, we discovered 13 DEmiRNAs and 273 DEmRNAs. To further understand the underlying mechanisms of these RNA molecules, we constructed a circRNA-miRNA-mRNA network for IDD. Our research provided several novel targets and molecular biomarkers for IDD. The hsa_circ_0091570, a sponge of miR-1307 that acted to regulate ISM1 expression in hepatocellular cancer, was screened from the ceRNA network constructed by 8 circRNAs, 13 miRNAs and 15 genes [32]. Three additional genes, CHEK1, TSGA10, and GSTCD, have been reported to downregulate genes that inhibit cellular proliferation or reduce cell numbers [33][34][35]. what's more, NPCs' proliferation and apoptosis play a key role in IDD [36]. Meanwhile, we con rmed that the level of these genes was downregulated in IDD. Thus, Hsa_circ_0091570 may suppress NPC proliferation by downregulating the expression of CHEK1, TSGA10, and GSTCD. The regulated network of hsa_circ_0091570-hsa-miR-508-5p-CHEK1/TSGA10/GSTCD was constructed and the RNAs (circRNAs and miRNAs) and genes in the network may all be novel targets for IDD. In parallel, Li et al. reported that miR-183-3p regulates autophagy in gastric cancer cells by the PI3K/AKT/mTOR pathway [37]. According to previous research, the PI3K/AKT/mTOR pathway has a signi cant effect on cell growth, survival, metabolism, proliferation, and angiogenesis [38]. Moreover, overexpressed PLG can promote the proliferation of keratinocytes [39]. We constructed a network with downregulated hsa_circ_0004565, upregulated hsa-miR-183-3p, and downregulated PLG. However, whether hsa_circ_0004565 depresses the proliferation of NPCs via downregulating the expression of PLG which further regulates PI3K/AKT/mTOR signaling is unclear yet. Some studies have shown that S100B regulates in ammation by increasing the expression of TNF-α and IL-1β during osteoarthritis [40].
Overexpression of TNF-α and IL-1β is also signi cant in the breaking of intervertebral discs [13]. However, we found that the level of S100B expression was decreased in IDD. The S100B-related ceRNA network (hsa_circ_0003183/hsa_circ_0032253/hsa_circ_0001293-hsa-miR-4534-S100B) was constructed. We discover that circRNA (hsa_circ_0003183/hsa_circ_0032253/hsa_circ_0001293) may act as an independent protective factor to inhibit IDD by decreasing the expression of S100B. Many miRNAs are involved in different diseases, such as hsa-miR-1827 and hsa-miR-508-5p [41,42]. In our research, we found that hsa-miR-1827 regulated the seven upregulated target genes, including GBP6, CEACAM8, C7orf34, KCNV2, OLFM4, NR5A1, and ZNF488. Studies have demonstrated that hsa-miR-1827 and its target(s) are closely associated with breast cancer, cervical cancer, and tongue squamous cell carcinoma [42][43][44][45][46][47][48]; we also explored whether these genes or miRNAs are involved in IDD, and then we further verify the connection of circRNA and IDD by constructing an integrated circRNA-miRNA-mRNA network. For example, the high expression of target genes CEACAM8 and OLFM4 plays a crucial role in the in ammatory process by increasing the level of in ammatory cytokines [49,50]. A great number of in ammatory cytokines promote extracellular matrix degradation and disorder of NPCs' inner or outer environment, which plays a signi cant role in IDD [51]. We also constructed the hsa_circ_0057552-hsa-miR-1827-GBP6/CEACAM8/C7orf34/KCNV2/OLFM4/NR5A1/ZNF488 network. We nd that hsa_circ_0057552 may take part in the in ammatory process and degeneration in IDD by enhancing the expression of CEACAM8 and OLFM4. Additionally, it was reported that ZNF488-overexpressing cells could increase the expression of collagen IV [52] which obviously decreased in IDD [53]. Hence, hsa_circ_0057552 may be a protective factor to resist IDD by increasing the expression of ZNF488. Although different miRNAs or genes involved in the network have been reported in different studies, this ceRNA-regulated network has not yet been reported in IDD. Moreover, the molecular circadian function of the hsa_circ_0057552-hsa-miR-1827-GBP6/CEACAM8/C7orf34/KCNV2/OLFM4/NR5A1/ZNF488 network remains unclear and will require further research.
KEGG analysis showed that the majority of genes were enriched in the p53 signaling pathway. Although p53 is a nuclear transcription factor that plays a vital role in regulating cell cycle progression, senescence, and cell death [54], its role in IDD will require further study. GO analysis revealed genes mostly focused on the regulation of the nervous system process, myelination, and steroid hormone receptor activity.

Limitation
This study has some restricted conditions. At rst, the number of samples is not enough. We need an additional cohort to further validate the expressed level of circRNAs, miRNAs and mRNAs in IDD. Secondly, based on our preliminary screening study, the interactions of identi ed ceRNA should be con rmed in future studies.  Figure 1 Heatmap of the 620 differentially expressed circRNAs in the 10 microarray datasets. By analysing the circRNA dataset, 620 differentially expressed circRNAs were identi ed (|log 2 (fold change (FC))| > 1 and an adjusted P value < 0.05 were considered Page 14/19

Figure 2
A total of 13 differentially expressed miRNAs were obtained by analysing the mRNA dataset (|log 2 (fold change (FC)) | > 1 and an adjusted P value < 0.05 were considered)  47 mRNAs were obtained from the intersection of mRNAs predicted and the differentially expressed mRNAs Figure 6 CeRNA network of circRNA-miRNA-mRNA interactions in IDD. Red indicates mRNAs, green indicates miRNAs, and blue indicates circRNAs. The deep color and light color represent up-regulation and down-regulation, respectively