Comprehensive Analysis of Differentially Expressed Circular RNAs in Patients with Senile Osteoporotic Vertebral Compression Fracture

Aim Circular RNAs (circRNAs) have been found to contribute to the regulation of many diseases and are abundantly expressed in various organisms. The present study is aimed at systematically characterizing the circRNA expression profiles in patients with senile osteoporotic vertebral compression fracture (OVCF) and predicting the potential functions of the regulatory networks correlated with these differentially expressed circRNAs. Methods The circRNA expression profile in patients with senile OVCF was explored by using RNA sequencing. The prediction of the enriched signaling pathways and circRNA-miRNA networks was conducted by bioinformatics analysis. Real-time quantitative PCR was used to validate the selected differentially expressed circRNAs from 20 patients with senile OVCF relative to 20 matched healthy controls. Results A total of 884 differentially expressed circRNAs were identified, of which 554 were upregulated and 330 were downregulated. The top 15 signaling pathways associated with these differentially expressed circRNAs were predicted. The result of qRT-PCR of the selected circRNAs was consistent with RNA sequencing. Conclusions CircRNAs are differentially expressed in patients with senile OVCF, which might contribute to the pathophysiological mechanism of senile osteoporosis.


Introduction
Osteoporosis has been the most common age-related bone disorder characterized by low bone mineral density (BMD) and deterioration of microarchitecture with increased bone fragility and subsequent susceptibility to fractures [1,2]. Senior population with osteoporosis often suffers from a fragility fracture, which will remarkably influence their life quality and increase family economic burden [3]. With the aging of the population, this burden has been increasing drastically. In addition, every low-energy fracture is linked to a raised risk of a subsequent fracture, which indicates a high risk of poor outcomes, and this risk is higher in men than women in the elderly [4]. In spite of its severe harm, osteoporosis is normally asymptomatic clinically at the early stage, and it is always ignored until the first fragility fracture happens. Clinically, primary osteoporosis can be divided into postmenopausal osteoporosis and senile osteoporosis. Senile osteoporosis, a common age-related degenerative disease, is characterized by impaired bone formation principally with a low bone turnover state [5,6], whereas postmenopausal osteoporosis is related to excessive bone absorption owing to estrogen deficiency [7]. Despite the advances in therapeutic strategies, the diagnosis of senile osteoporosis at an early stage and the poor prognosis remains challenging. Therefore, there is an urgent need for a better understanding of the pathogenesis of senile osteoporosis in order to identify novel biomarkers and to develop new treatment strategies for the prevention and treatment of senile osteoporosis.
Emerging evidences demonstrate that microRNAs (miRNAs) are essential for bone remodeling and regeneration through the posttranscriptional regulation of gene expression [8][9][10]. With the development of next-generation sequencing technology, more noncoding RNAs have been discovered. As a large class of noncoding RNAs, circular RNAs (circRNAs) have multiple biological properties, such as high conservation, widespread expression, specific expression, and high stability [11][12][13]. Furthermore, some studies have reported that cir-cRNAs were involved in age-related disease and played important roles in many biological processes as miRNA sponges [14]. Nevertheless, information about circRNAs on senile osteoporosis is very limited.
In this study, we used RNA sequencing, a powerful and unbiased approach to identify and analyze the differences in expression profiles of circRNAs and their regulatory interaction networks in patients with senile osteoporotic vertebral compression fracture (OVCF). These findings may help elucidate the mechanism underlying senile osteoporosis development and uncover novel diagnostic biomarkers for senile osteoporosis.

Materials and Methods
2.1. Patients and Samples. This study was approved by the Ethics Committee of our hospital. Fresh peripheral blood and written informed consent were attained from patients with newly diagnosed senile OVCF and healthy volunteers receiving a regular physical examination at the same hospital. The diagnosis for senile OVCF was confirmed through magnetic resonance imaging, dual-energy X-ray absorptiometry, and the history of patients. To decrease sources of circRNA expression heterogeneity not connected with disease status, selected subjects were all male, with similar age and body mass index. 2 ml blood was collected in ethylenediaminetetraacetic acid (EDTA) tube and mixed with three volumes of RNASaferTM LS Reagent (Magen, Guangzhou, China), then stored at -80C for batch processing.

RNA Extraction and High-Throughput Sequencing
Analysis. Total RNA was extracted by Hipure PX Blood RNA Mini Kit (Magen) according to the manufacturer's protocol. The concentration and integrity of the isolated RNA were evaluated with Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, California) and Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA). High-throughput sequencing was performed by Illumina Hiseq Xten to determine cir-cRNAs in patients with senile OVCF (n = 3) and the controls (n = 3).

Data Analysis.
We used Bowtie2 version 2.1.0 to map the reads to the latest UCSC transcript set and RSEM v1.2.15 to estimate the gene expression levels. We performed TMM to normalize the gene expression. We used edgeR program to perform quantile normalization and subsequent data processing. We used STAR to map the reads to the genome. We used DCC to identify the circRNAs and estimate their expression. CircRNAs were considered differentially expressed according to P < 0:01 and absolute log fold change >2.

Bioinformatics Analysis. Gene ontology (GO), Kyoto
Encyclopedia of Genes and Genomes (KEGG), and Reactome pathway analyses of the differentially expressed cir-cRNAs were performed by clusterProfiler R package. The related target miRNAs of circRNAs were identified by miRanda (v3.3). The circRNA-miRNA-mRNA network was constructed by Cytoscape software (v3.7.1).

Expression Profiles of Differentially Expressed circRNAs.
RNA sequencing was performed to identify circRNAs in peripheral blood in three patients with senile OVCF and three healthy controls. We discovered 884 differentially expressed circRNAs (554 upregulated and 330 downregulated). The variations in circRNA expression were demonstrated by scatter maps and volcano plots (Figures 2(a) and 2(b)). Hierarchical clustering heat map analysis presented a recognizable circRNA expression profile in these two groups ( Figure 2(c)).

GO, KEGG, Reactome
Analysis, and circRNA-miRNA Network Analysis. To explore the role of differentially expressed circRNAs, GO, KEGG, and Reactome analyses and circRNA-miRNA network analysis were performed. The top three GO terms in the biological process (BP) were "DNA replication, nuclear chromosome segregation, mitotic nuclear division", and "peptidyl−serine phosphorylation, protein methylation, histone methylation" for upregulated and downregulated circRNAs, respectively. The top three GO terms in the cellular component (CC) were "chromosome, centromeric region, microtubule organizing center part, condensed chromosome", and "ribonucleoprotein granule, nuclear transcription factor complex, cytoplasmic ribonucleoprotein granule" for upregulated and downregulated circRNAs, respectively. The top three GO terms in the molecular function (MF) were "histone binding, protein C−terminus binding, modification−dependent protein binding", and "Ras guanyl−nucleotide exchange factor activity, phosphatidylinositol binding, histone methyltransferase activity" for upregulated and downregulated circRNAs, respectively.
The top 10 most enriched GO terms of the upregulated and downregulated circRNAs in BP, CC, and MF are shown in Figures 3(a) and 4(a), respectively. The top three KEGG pathways were "RNA transport, cell cycle, and Hippo signaling pathway", and "MAPK signaling pathway, human T−cell leukemia virus 1 infection, and T cell receptor signaling pathway" for upregulated and downregulated circRNAs, respectively. The top three Reactome pathways were "M Phase,    BioMed Research International transcriptional regulation by TP53, and DNA repair", and "chromatin organization, chromatin modifying enzymes, and regulation of TP53 activity" for upregulated and downregulated circRNAs, respectively. The top 15 KEGG and Reactome pathways related to the upregulated and downregulated circRNAs are shown in Figures 3(b), 3(c), 4(b) and 4(c), respectively. The circRNA-miRNA coexpression network analysis of these differentially expressed circRNAs is shown in Figure 5. 3.3. qRT-PCR Validation and circRNA-miRNA-mRNA Network of the Selected circRNAs. The result of qRT-PCR was consistent with the sequencing data and showed significant differences in expression for circ_0079449 and circ_ 0122913 (Figure 1(b)). Next, we showed the five miRNAs that most potentially bind to circ_0079449 or circ_0122913, and the top five target genes to each miRNA in the illustration (Figure 1(c)).

Discussion
In the present study, after strict filtering of approximately 180 million reads, we detected 76547 circRNAs, including 884 circRNAs that were differentially expressed between patients  9 BioMed Research International miRNA is one kind of most studied regulators in the bone remodeling process and likely to be the diagnostic biomarker for osteoporosis [15]. Previous studies revealed miRNAs might play critical roles in the development of osteoporosis. For example, Mäkitie et al. [16] found circulating miRNAs in serum reflected status of WNT1 mutation, which leads to skeletal pathologies. Zhang et al. [17] found that miRNA-9-5p was lowly expressed in peripheral blood of osteoporosis patients and promoted the occurrence and progression of osteoporosis through inhibiting osteogenesis and promoting adipogenesis via targeting Wnt3a. In addition, miRNA-based therapies have shown great potential to become novel approaches against osteoporosis and associated fractures [8]. As miRNA sponges, circRNAs can adsorb

10
BioMed Research International miRNAs and eliminate the inhibitory effects of miRNAs on their target genes, leading to upregulation of target genes [18]. Besides, circRNAs are more stable than linear RNA due to their closed-loop structure [19]. CircRNAs have been considered potentially promising blood biomarkers in some musculoskeletal-related disease [20][21][22]. However, information about circRNAs on osteoporosis is limited. Zhao et al. [23] used chip microarray analysis to investigate circRNA expression in postmenopausal osteoporosis and healthy menopausal women. Their data identified 203 upregulated and 178 downregulated circRNAs and revealed circ_0001275 was a potential diagnostic biomarker according to qRT-PCR results. Huang et al. [24] performed microarray analysis and PCR validation, and their result showed that circ_0006873 and circ_0002060 were linked to the low BMD state. In our study, rather than investigating by microarray which could identify only thousands of circRNAs, we chose to sequence all circRNAs that were attained in blood samples from patients with newly diagnosed senile OVCF and healthy controls.
Even though we identified the differentially expressed cir-cRNAs in patients with senile OVCF, the underlying mechanism is still poorly understood. The GO terms provide computable knowledge regarding the activities of gene products and have been good predictors of them [25]. KEGG is an integrated database for biological interpretation of genome sequences and widely used in the enrichment analysis [26]. Reactome is a knowledgebase of biomolecular pathways and can be combined with other databases [27,28]. Therefore, we systematically predicted the circRNA-related functions and corresponding pathways in senile osteoporosis with GO, KEGG, and Reactome analyses. In the present study, we found several pathways in the enrichment analysis were associated with osteoporosis. These pathways include the mitogen-activated protein kinase (MAPK) signaling pathway, DNA damage, DNA repair, and Th17 cell differentiation. Previous studies reported that the MAPK signaling pathway plays important roles in bone metabolism [29][30][31]. Other pathways were also involved in the development of osteoporosis, such as DNA damage, DNA repair, and Th17
There are some limitations in this study. First, our patients were carefully selected. We chose male and age, BMI-matched subjects from a single center. Although this decreased heterogeneity for this initial study, it would limit the generalizability of our findings. Second, the sample of this study is small. Therefore, a larger number of clinical samples are needed to confirm our results. Furthermore, the mechanism and function of the differentially expressed circRNAs should be further verified by strict molecular biological experiments research and studied deeply in future work.
In summary, we reported the circRNA expression profile in patients with senile OVCF by RNA sequencing analyses. Their regulatory interaction networks were constructed based on the sequencing data. These original discoveries might provide some clues for the biological functions of circRNAs in the pathophysiological mechanism of senile osteoporosis.

Data Availability
The RNA-seq data have been uploaded to NCBI SRA database (Number: PRJNA562388).

Ethical Approval
The Institutional Review Board of The Xiaoshan First People's Hospital approved this study (Protocol Number: 2019-XS-001).