Comprehensive Analysis of Aberrantly Expressed CircRNAs in The Ossication of The Posterior Longitudinal Ligament

With producing the severe neurological symptoms and movement disorders, ossication of the posterior longitudinal ligament(OPLL) have attracted attention in mounting numbers. Increasing evidences show that circRNAs act a vital regulatory role in the biological process of human disease development. However, the potential mechanism of circRNAs in OPLL was not fully known. In this study, 3 OPLL patients and 3 controls were selected for RNA-seq analysis and the expression proles of circRNAs were established. TargetScan and miRanda's miRNA target prediction software were used to predict the circRNA-microRNA interaction. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used to identify the biological process and key pathways. we used STRING database to analyze and construct a Protein-Protein Interaction (PPI) network. Several circRNAs were randomly selected for quantitative real-time PCR (qRT-PCR) validation. Microarray data proling indicated that 489 circRNAs (234 up and 255 down) was highly differentially expressed with setting screening criteria (FC ≥ 2.0 and p value ≤ 0.05). GO enrichment and KEGG pathways analysis revealed some processes and pathways might inuence the progress of OPLL such as osteoclast differentiation, vasopressin-regulated water reabsorption, Wnt, MAPK, VEGF signaling pathway. We have obtained some signicant genes such as AKT3, ACVR1, RBPJ, NFATC1, PTK2, SLC8A1. The osteoclast activity regulated by AKTs may compensate for the promotion of bone ectopic hyperplasia on OPLL. Long non-coding RNA MALAT1 functions as miR-1 sponge to regulate Connexin 43-mediated ossication of the posterior longitudinal ligament.

Introduction Ossi cation of the posterior longitudinal ligament (OPLL) is a progressive disease of abnormal calci cation of the spinal ligament [1]. The ossi ed tissue encroached on the volume of the spinal canal resulting in compression of the spinal cord and serious nervous system problems [2,3]. The stiff ligament limits the normal curvature and mobility of the spine, which seriously affects the patient's quality of life [4]. Compared to North America and Europe, East Asians have a higher prevalence of OPLL [5,6].
Research on OPLL in recent decades have shown that it is a multifactorial disease, in which both genetic and environmental factors play a vital role [7]. Cardiovascular disease, diabetes, cerebrovascular disease are signi cantly associated with OPLL [8]. Finding the molecular mediator of posterior longitudinal ligament ossi cation and revealing its mechanism of action will help expand our understanding of the pathogenesis of OPLL and provide a theoretical basis for OPLL to prevent and even reverse the occurrence of OPLL.
Studies have found that multiple mechanisms and signaling pathways are involved in the ossi cation process of posterior longitudinal ligament broblasts. Environmental factors such as stress stimulation can up-regulate the expression of broblast gap junction protein 43 (Cx43), activate extracellular regulatory protein kinase 1/2 (ERK1/2), mitogen protein kinase (MAPK), and JNK pathways to make osteoblast differentiation [9]. The transcription factor osterix promotes the differentiation of broblasts into osteoblasts by activating the bone formation regulator Runx2 while inhibiting the β-catenin signaling pathway [10]. Transforming growth factor-β (TGF-β) and other cytokines are also involved in the process of broblast differentiation into osteoblasts [3]. With the development of molecular biology research, some non-coding RNA has been found to be involved in the occurrence of OPLL.
MicroRNAs (miRNAs) are small non-coding single-stranded RNAs with a length of about 18-22 nt, which bind to the 3'-UTR of the target gene to inhibit the translation of the target gene mRNA or directly degrade it. Abnormal expression of miRNAs is closely related to a variety of bone diseases, including osteoarthritis, osteosarcoma, osteoporosis, and lumbar degenerative disease [11][12][13]. For example, miR-122 inhibits the proliferation/differentiation of osteoblasts in osteoporosis by activating the PCP4mediated JNK pathway [13]. Long non-coding RNAs (LncRNAs), a kind of non-coding RNA with a length greater than 200 nucleotides, plays a complex regulatory function in gene expression [14]. A previous study showed that Long non-coding RNA MALAT1 acts as a miR-1 sponge and regulates the OPLL mediated by Cx43 [15].
Circular RNAs(circRNAs), a class of non-coding RNAs, have continuous closed loop structure formed by covalent bonds that widely found in eukaryotic cells and do not contain 5'caps and 3'polyA tails. In the past, circRNAs was thought to be RNA transcribed from chromosomal DNA, formed by splicing errors mediated by spliceosome, and was a kind of non-functional product or intermediate product. However, it seems that circRNAs are involved in the occurrence of diseases by regulating their downstream pathways by combining with microRNAs. CircRNAs can bind to the corresponding miRNAs, releasing the inhibition of the corresponding target gene, and play the role of competitive endogenous RNA(ceRNA) [16,17]. For example, circRNA LARP4 regulates miRNA-424 through ceRNA mechanism, which can regulate the occurrence and prognosis of osteosarcoma [18].
Up to now, there is no research report on which role circRNAs plays in the mechanism of OPLL. In this study, we performed a circRNAs microarray analysis for the rst time to establish the expression pro le of circRNAs in the posterior longitudinal ligament of OPLL patients. We used bioinformatics prediction methods to annotate the differentially expressed circRNAs. Our ndings provide new evidence for exploring the molecular mechanism of CircRNAs in the pathogenesis of OPLL.

Sample collection
The specimens of the posterior longitudinal ligament were provided by patients from the First A liated Hospital of Harbin Medical University. Every patient underwent computed tomography (CT) examination, and OPLL was diagnosed if the thickness of the ectopic OPLL was more than 2 mm [7]. We collected OPLL tissue (n=3) and normal PLL tissue (n=3) during the operation. The study protocol was approved by the Institutional Review Broad of Harbin Medical University, Harbin, Heilongjiang Province, China. Informed consent was acquired from the all study participants. All research methods are within relevant ethical principles.

Sequencing of circRNAs and Microarray analysis
The purity and concentration of total RNA samples were determined with NanoDrop ND-1000. The total RNA of each sample was treated with RNase R to enrich circular RNA. The enriched circular RNA was ampli ed and transcribed into uorescent cRNA utilizing random primer. After hybridizing, incubating and washing, raw data was extracted by Agilent Feature Extraction software.
A series of data processing including quantile normalization were performed using the R software limma package. The circRNAs that at least 3 out of 6 samples have ags in "P" or "M" (de ned by GeneSpring software) were retained for further differential analyses. The microarray data have been uploaded in the National Center for Biotechnology Information Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [GSE164546: GEO].

Differential expression analysis
Comparing two groups of pro les, the fold change (FC, the ratio of the group average) between each circRNA was calculated . The statistical signi cance of the difference may be conveniently estimated with the t-test. CircRNA with FC ≥ 2.0 and p value ≤ 0.05 was selected as the signi cantly differentially expressed . Microsoft Excel's Data/Sort & Filter functionalities were used to lter the analysis output, and the differentially expressed circRNA was sorted according to the fold change, p value, etc.

Annotation for CircRNAs binding site
TargetScan and miRanda's miRNA target prediction software were used to predict the circRNA-microRNA interaction, and the differentially expressed circRNAs within all the comparisons were annotated in detail with the circRNA-miRNA interaction information. The miRNA response elements (MRE) predicted by the two algorithms were retained. We predicted the top 5 miRNAs of each circRNA based on the pairing scores, and constructed a miRNA-circRNA network through cytoscape_v3.8.0 software (http://cytoscape.org/).

GO annotation and KEGG pathway enrichment analyses of the DEGs
Gene Ontology (GO), including biological processes (BPs), cellular components (CCs), molecular functions (MFs) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses of DEGs were performed using the Database for Annotation, Visualization and Integrated Discovery (P<0.05 as the criteria for enrichment signi cance).

Protein-Protein Interaction (PPI) Network construction
In order to further explore the interaction of different target genes and the molecular mechanism of OPLL, we used STRING (https://string-db.org/) to analyze and construct a PPI network. Then the interaction network and the top 50 hub genes were visualized by cytoHubba in Cytoscape. The nodes in the network were represented as target genes, and the lines between two nodes were denoted interactions.

Reverse transcription and quantitative PCR (RT-qPCR)
Total RNA were extracted from OPLL tissues or cultured cells using Trizol reagent (Invitrogen, CA), and a reverse-transcription kit (Takara, Japan) was used to synthesize total RNA into cDNA. To measure the expression of circular RNAs, linear RNA was degraded with RNase-R. Real-time PCR was performed with SYBR PreMix Ex Taq kit (Takara) and applied biological system (Foster City, CA). Glyceraldehyde-3phosphate dehydrogenase (GAPDH) was used as the internal control of circRNAs. The fold change was calculated using the 2-ΔΔCT method.

Statistical analyses
All statistical analyses were analyzed with SPSS 17.0(Chicago, IL, USA). The t test was used for comparison of samples between test and control groups and the Benjamini Hochberg FDR (the FDR cutoff was 0.05) was used for multiple-testing correction. P<0.05 was considered statistically signi cant.

Microarray data
Comparing the results of OPLL groups and control groups, analyzing and normalizing the results, we found 13,617 circRNAs. Through the setting screening criteria (FC ≥ 2.0 and p value ≤ 0.05), we identi ed 234 up-regulated circRNAs and 255 down-regulated circRNAs. Based on these data, a hierarchical cluster analysis was performed to generate a heat map of differentially expressed circRNAs.
( Figure 1A) The top 10 up-regulated and top 10 down-regulated circRNAs were listed in Table 1. Among the differentially expressed circRNAs, which derived from linear transcript exons accounted for the majority with 414 (84.7%) ( Figure 1B, Table 2). We used a scatter plot to compare the distribution of differentially expressed circRNAs between samples. (Figure 1C) Based on fold change values and p values, we used a volcano plot to show the overall data differential expression distribution of circRNAs in ligament cells of OPPL patients and normal patients. ( Figure 1D) After data standardization, the median of each chip were at the same level, and the probes distribution have been also relatively similar, which means that the chip data quality was reliable ( Figure 1E). In summary, comparing the circRNAs expression data of OPLL and healthy ligament tissue cells, we found signi cant differences in circRNA expression between the two groups.

CircRNAs binding site prediction
The differentially expressed circRNAs were predicted by TargetScan and miRanda targets. we retained the top 5 predicted miRNAs according to their scores. (Table 3) MiRNA-circRNA network was constructed using cytoscape_v3.8.0. (Figure 2) CircRNA/miRNA interaction information provided detailed annotations for all differentially expressed circRNAs in the comparison. (Figure 3) Using circRNA-targeted miRNA-gene network, we chose four differentially expressed cirRNAs to show the relationship between circRNA/miRNA. (Figure 4, Table 4)

Enrichment analyses of the target genes
We performed GO enrichment and KEGG pathway analysis of the 489 differentially expressed circRNA target genes to further clarify their functional characteristics. ( Figure 5) In our research, we found that it is involved in biological processes, cellular component and molecular functions such as ossi cation, tissue morphogenesis, site of polarized growth, protein kinase activity, ubiquitin-like protein binding, and etc. Several signi cant pathways were also discovered, including osteoclast differentiation, Wnt signaling pathway, vasopressin-regulated water reabsorption, MAPK signaling pathway, VEGF signaling pathway, etc. (Table 5) 3

.4 PPI Network Analysis
The STRING online database was used to distinguish the connections between the target genes, and a network composed of 206 nodes and 707 edges was obtained, with a con dence score > 0.4 as signi cant. We used ClusterOne in Cytoscape for gene clustering. The top 50 scoring genes were represented by orange circles, and the color shades were used to correlate signi cance. (Figure 6 a. b)

qRT-PCR validation of sequencing data
Six differently expressed circRNAs were randomly selected for qRT-PCR experiments in OPLL patients and control samples, including three up-regulated circRNAs (hsa_circRNA_101725, hsa_circRNA_000950, hsa_circRNA_405283) and three down-regulated cirRNAs (hsa_circRNA_048764, hsa_circRNA_100811, hsa_circRNA_000675). The results of qRT-PCR were consistent with the results of circRNA-seq data, indicating relative reliability of the sequencing data. (Figure 8)

Discussion
With the increasing incidence of ossi cation of the posterior longitudinal ligament, the problems of di cult treatment and poor prognosis needs to be overcome urgently [19]. Due to regional differences and ethnic characteristics, OPLL has been restricted in epidemiology and pathogenesis research [20]. At present, surgery is widely adopted as the main treatment method for myelopathy caused by OPLL, but it does not fundamentally solve the problem of heterotopic ossi cation [21,22]. The development of proteomics and high-throughput sequencing tools greatly expanded our understanding of the molecular mechanism of OPLL. In previous studies, it was revealed that some non-coding RNAs play an important role in ossi cation [15,23]. Several treatment concepts for dysregulated genes have been proposed in OPLL animal models, and encouraging results have been achieved [22].
Up to now, there is no report on the in uence of circRNAs in OPLL. It was not until the derivation of new sequencing and transcript calculation methods that more functions of circRNAs were revealed [24]. To this end, we collected samples from OPLL patients for high-throughput sequencing and performed bioinformatics analysis. We found that many circRNAs were abnormally expressed during OPLL, which suggests that circRNAs plays an important role in the occurrence and development of OPLL. Further research found that most of the differentially expressed circRNAs are derived from exons, which indicate that circRNAs mostly play a regulatory role in the cytoplasm. However, the speci c shear mode and function retention of circRNAs still needed to be explored [25,26].
This study used bioinformatics methods to explore cell functions and biological pathways related to OPLL. We found that circRNAs are involved in osteoblast differentiation, regulation of ossi cation, tissue morphogenesis, and cellular protein metabolic process, which indicates that circRNAs may regulate the abnormal ossi cation of posterior longitudinal ligament through molecular function and cellular component. Such as the Wnt/β-Catenin signaling pathway is involved in the occurrence of osteoarthritis and myeloma [27]. The MAPK/SPARC/ERK signaling pathway regulates heterotopic ossi cation [28], and the YAP/TAZ signaling pathway plays a role in enhancing bone growth and angiogenesis [29]. The ErbB signaling pathway promotes the development of cartilage and bone [30]. In vitro culture of broblasts under the action of BMP-2 combined with a certain concentration of tumor necrosis factor (TNF-α) can detect the secretion of ALP and OC and the formation of bone nodules [31]. In the process of osteoporosis, inhibiting the expression of miR-185 can activate the BMP/Smad signaling pathway to increase bone formation [32]. Our result of KEGG pathway enrichment also veri ed that the mentioned pathways may also participate in the occurrence and development of OPLL. We have obtained some signi cant genes such as AKT3, ACVR1, RBPJ, NFATC1, PTK2, SLC8A1 and etc, which have also been con rmed to be speci cally expressed in other related diseases. AKT3 can regulate cartilage ossi cation In mouse chondrocytes, and NF-κB/MAPK-mediated osteoarthritis is further promoted by ACVR1 [33][34][35].
Previously, most of the research on OPLL started from the direction of bone ectopic hyperplasia. lncRNA XIST aggravated the ossi cation of ligaments through the miR-17-5p/BMP2 signaling pathway [36], and miR-10a promoted the bone formation of mouse ligament cells through the ID3/RUNX2 pathway [37]. However, our research found that the biological process of osteoclasts also seems to be involved in the development of OPLL. Previous studies have proved that the targeted knockout of AKT1 promoted the differentiation of osteoblasts, but knocking out AKT2 has reversed the differentiation of osteoblasts and restored normal osteoclast production [38]. We speculate that OPLL is the result of an imbalance in the tight coupling between bone formation and osteoclastosis. Osteoclastic activity may compensate for reducing the promotion of bone ectopic hyperplasia on OPLL. It seems that promoting local osteoclasts process can be a new way to treat OPLL.
This study also has several limitations: 1. We select scarce number of samples for microarray analysis, which may affect the validity of the sequencing results. 2. Due to the limitations of research methods, we only predicted the function of differentially expressed circRNAs, but could not clarify the detailed mechanism of action in these circRNAs. Subsequently, the functional mechanism of circRNAs can be explored by RIP-QCR, RNA pull down, Luciferase assay, etc methods [39,40].

Conclusions
It is the rst study for portraying circRNA expression pro les of OPLL. This provides new ideas for the prevention and alleviation of OPLL . The speci c functions and mechanisms of circRNAs in OPLL should be further veri ed to provide more clinical treatment options for OPLL patients.

Declarations
Ethical Approval and Consent to participate The study protocol was approved by the Institutional Review Broad of Harbin Medical University, Harbin, Heilongjiang Province, China. Informed consent was acquired from the all study participants. All research methods are within relevant ethical principles.

Consent for publication
All patients or relatives of the organization donors were fully understood the purpose and plan of the study, and agreed to publish the results of the study. Subsequently, a written informed consent was obtained.
Availability of supporting data The microarray data have been uploaded in the National Center for Biotechnology Information Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [GSE164546: GEO] JAG and QCG conceived the study; ZZZ, QCG, TS, WLT, ZGY, YHH, FCL and HTS obtained the samples and carried out experiments; and QCG, JAG, FCL and ZGY analyzed data. All authors were involved in writing the paper and had nal approval of the submitted and published versions. Figure 4 The network includes 11 miRNAs, 4 circRNAs, 106 target genes, 272 edges. The orange nodes represented miRNA, the red nodes represented circRNAs, and the green nodes represented target genes.