Identification and comparison of novel circular RNAs with associated co-expression and competing endogenous RNA networks in postmenopausal osteoporosis

Background Circular RNAs (circRNAs) are emerging as crucial regulators in various human diseases. So far, the expression profile and regulatory mechanism of circRNAs in postmenopausal osteoporosis (PMOP) are less studied and should be deciphered urgently. Herein, we aimed to reveal key circRNAs affecting PMOP and clarify their compounding regulatory actions. Methods To reveal key circRNAs affecting PMOP and clarify their compounding regulatory actions, whole transcriptome sequencing and bioinformatics analysis were performed to identify differentially expressed circRNAs (DECs). The expression pattern and regulatory networks of DECs in peripheral blood mononuclear cells (PBMCs) were unearthed. Results A total of 373 DECs comprising 123 intronic, 100 antisense, 70 exonic, 55 intergenic, and 25 sense-overlapping circRNAs were identified. Among these, 73 circRNAs were upregulated and 300 were downregulated. These DECs exerted pivotal functions in the pathogenesis of PMOP as demonstrated by Gene Ontology (GO) annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The circRNA-miRNA-mRNA co-expression network comprising 28 DECs, 145 miRNAs, and 175 differentially expressed mRNAs predicted the possible mechanism of the pathogenesis and progression of PMOP. Conclusion The results of the present study provided a further comprehension of circRNA-associated competing endogenous RNA regulatory mechanism in PMOP. The steadily expressed and disease-specific DECs may serve as promising diagnostic and prognostic biomarkers for PMOP. Supplementary Information The online version contains supplementary material available at 10.1186/s13018-021-02604-1.


Introduction
Osteoporosis, a pervasive public health concern worldwide, manifests as the depletion of bone mineral with structure deterioration of bone tissue and results in a predisposition to fragility fractures, especially in women [1]. National registry statistics show that osteoporosis affects approximately 40.1% of all postmenopausal women in China and the resulting fragility fractures are associated with significant morbidity, mortality, and financial implications [2]. To alleviate the public health burden of postmenopausal osteoporosis (PMOP), it is a prerequisite to screen out the population at high risk for fracture, which entails the identification of more rapid, specific, and sensitive bone turnover markers (BTMs) [3].
Circular RNAs (circRNAs), a relatively new class of non-coding RNAs (ncRNAs), form a continuous cycle of covalent closures and are highly expressed in the eukaryotic transcriptome [4]. Accumulating evidence has revealed the regulatory roles of circRNAs, such as functioning as micro-RNA (miRNA) sponges, serving as scaffolds in the assembly of protein complexes, regulating the expression of parental genes, and modulating alternate splicing, as well as RNAprotein interactions [5][6][7][8][9]. Moreover, circRNAs are expressed more steadily and abundantly than the standard linear transcription of homologous genes due to the closedloop structure, which prevents degradation by RNA exonuclease [8]. In this way, disease-specific and tissue-specific circRNAs could act as biomarkers for the early diagnosis and prediction of certain diseases theoretically [10]. So far, circRNAs are demonstrated to be expressed in a tissue specific manner and regulate various pathophysiological events, such as organogenesis, tumorigenesis, and organ development [11][12][13][14]. However, specific and sensitive cir-cRNA biomarkers for PMOP have not been fully established. Herein, the present study aimed to uncover accurate circRNA biomarkers and provide clues for exploring the underlying mechanism of PMOP. With the aid of whole transcriptome sequencing, differentially expressed cir-cRNAs (DECs) between PMOP patients and normal controls were identified. Functional annotation and protein-protein interaction network constructions were also performed to explore the biological functions of targeted DECs. These findings will help elucidate the mechanism by which circRNA regulates the balance between osteogenesis and osteoclastogenesis. Moreover, the present study is foundational for subsequent studies on circRNAs in PMOP and may offer insight into prevention and new treatment targets for PMOP.

Ethics statement
This study was approved by the Medical Ethics Committee of Local Institution. Written informed consent was obtained from each participant before enrollment.

Subjects
Postmenopausal patients who received percutaneous kyphoplasty (PKP) surgery or bone mineral density (BMD) examination in the clinic department of our institution between December 2019 and January 2020 were evaluated. Participants were enrolled in this study using the following inclusion and exclusion criteria. The inclusion criteria included (1) age between 55 and 65 years old; (2) at least 1 year after natural menopause; (3) T-scores < −2.5 standard deviation SD (PMOP group) or T-scores > 0 SD (control group) at their lumbar vertebrates [3]; (4) received lumbar PKP surgery (PMOP group) or without PMOP-associated fractures (control group). The exclusion criteria included (1) secondary osteoporosis due to metabolic, blood, thyroid, tumor, drug, and nutritional disorders; (2) premature menopause less than 45 years old; (3) patients who received ovariectomy; (4) had taken calcium, vitamin D, bisphosphonates, and estrogen in the past 3 months. Finally, three paired fresh venous blood samples were collected for whole transcriptome sequencing and 30 paired samples were used for quantitative evaluation using realtime quantitative polymerase chain reaction (qRT-PCR). The demographic and clinical characteristics of the six participants are summarized in Table 1.

Isolation of peripheral blood mononuclear cells (PBMCs)
Early in the morning, fresh venous whole blood was collected from each participant in a 2.5 ml PAXgene tube (BD, Franklin Lakes, NJ, USA), followed by the isolation of PBMCs within 6 h of collection. PBMCs were isolated at room temperature (18-20°C) using Ficoll-Paque PLUS reagent (GE Healthcare, Piscataway, NJ, USA) according to the manufacturer's instruction.

RNA isolation
Total RNA was extracted from the PBMCs and enriched using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instruction. RNA quality and quantity were evaluated on a Nanodrop spectrophotometer (ND-1000, Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity and genomic DNA (gDNA) contamination were determined by gel electrophoresis.

RNA sequence analysis
Paired-end reads were harvested from Illumina HiSeq 4000 sequencer and quality controlled by Q30 (P < 0.001). After 3′ adaptor trimming and low-quality reads removal using the cutadapt software (version 1.9.3), high-quality trimmed reads were aligned to the reference genome/transcriptome (UCSC hg19) guided by the Ensembl Gff gene annotation file with the HISAT2 software (version 2.0.4) and the STAR software (version 2.5.1b). Detection and identification of circRNAs were performed with the DCC software (version 0.4.4). The edgeR software (v3.16.5) was used to normalized the data and obtain the expression profiles of circRNAs in terms of the fragments per kilobase of transcript per million fragments mapped (FPKM) and the cuffdiff software (version 2.2.1) for the expression profiles of mRNAs. Subsequently, the fold-change (FC) and P value were calculated based on FPKM. DECs were compared between the PMOP and the control groups following the criteria of |log(FC)| ≥ 2 and P value < 0.05. Differentially expressed mRNAs (DEMs) were compared between the PMOP and the control groups following the criteria of |log 2 (FC)| ≥ 2 and P value < 0.05.

qRT-PCR
After isolating and enriching the RNA as aforementioned, 8 circRNAs were selected and validated using qRT-PCR. Total RNA was reverse transcribed into cDNA using SuperScript III reverse transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.). Then, a 2X PCR master mix (CloudSeq Biotech, Inc., Shanghai, China) was utilized to proceed with the qRT-PCR with selected RNAs and internal reference, β-actin (ACTB) of the samples. Primers for each circRNA were designed with the CircPrimer software (version 1.2) and listed in Supplementary Table 1.
The program was initiated by denaturation at 95°C for 1 min, and thermal cycling conditions were set to 40 cycles at 95°C for 15 s and 60°C for 30 s. Data were analyzed using the 2 −ΔΔCt method.

Functional group analysis
Database for Annotation, Visualization and Integrated Discovery (DAVID), a bioinformatics web server (https://david.ncifcrf.gov/), was referred to explore the potential functions of the liner transcripts [15]. Gene Ontology (GO) analysis constituted with three domains, namely, the "biological process," "cellular component," and "molecular function," was performed. The pathwayenrichment analysis was also performed according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

circRNA-miRNA-mRNA network prediction
The overlapping miRNAs between the predicted cir-cRNAs target miRNAs and the predicted miRNAs targeting the DEMs were adopted for further analysis. miRNA binding sites were predicted with miRcode (http://www. mircode.org/). Finally, the circRNA-miRNA-mRNA regulatory network was constructed using a combination of circRNA-miRNA pairs and miRNA-mRNA pairs and visualized using the Cytoscape 3.6.0 software. Values were compared between the PMOF and control groups, and the result was statistically significant (P < 0.01)

Statistical analysis
GraphPad Prism version 5.0 (GraphPad Software, San Diego, CA, USA) was used to perform a two-tailed independent t test or analysis of variance (ANOVA) for data processing. Results with P values < 0.05 were considered statistically significant.

Identification of DECs between PMOP and control individuals
In total, 8459 meaningful circRNAs were identified after mapping the sequencing reads to the human genome and abnormal expression was not observed in all six samples (Fig. 1a, Table 2). The sequence length, expression pattern, and distribution mapping of these cir-cRNAs were assessed. The length of these circRNAs ranged from 106 to 94088 bp, with a frequent distribution length of < 200 bp (75%) and 200-300 bp (11%) (Fig. 1b). These circRNAs were unevenly located in all chromosomes, especially on chr1, chr2, chr3, and chr5 ( Fig. 1c). They were categorized into five types based on their functions, including intronic, antisense, intergenic, exonic, and sense-overlapping circRNAs. Specifically, intronic cirRNAs constituted the largest portion (39%), followed by 25% of antisense, 17% of exonic, 13% of intergenic, and only 6% of sense overlapping (Fig. 1d). A hierarchical clustering approach was applied to determine the consistency of all samples. Volcano plot and heat-map analysis were used to identify DECs (Fig. 2a, b). Strikingly, 373 potential DECs comprising 123 intronic, 100 antisense, 70 exonic, 55 intergenic, and 25 senseoverlapping circRNAs were identified (Fig. 2c). Among these, 73 circRNAs were significantly upregulated and 300 were significantly downregulated in PMOP patients compared to the control individuals. The top 15 upregulated and downregulated circRNAs are listed in Table 3.

Functional annotation of DECs
To explore the potential functions of these DECs, GO annotation and KEGG pathway analyses were separately performed with target coding genes of significantly upregulated and downregulated circRNAs. GO analysis comprised three elements, including the "cellular component," "biological process," and "molecular function." The top 10 dysregulated GO processes of each subgroup were analyzed based on the dysregulated, enriched cir-cRNAs derived from gene annotation. The upregulated circRNAs were found to be mostly enriched as follows: G2 DNA damage checkpoint, transcription from RNA polymerase II promoter, and DNA damage checkpoint in the "biological process" subgroup; intracellular part, intracellular and membrane-bounded organelle in the "cellular component" subgroup; and nicotinamide adenine dinucleotide (NAD+) binding, protein binding, and NAD binding in the "molecular function" subgroup ( Fig. 3a, Supplementary Table 2). The significantly downregulated circRNAs were found to be mostly enriched as follows: presynaptic membrane organization, nervous system development, and positive regulation of filopodium assembly in the "biological process" subgroup; filopodium, growth cone, and site of polarized growth in the "cellular component" subgroup; and glutamate receptor activity, ion-gated channel activity, and protein tyrosine phosphatase activity in the "molecular function" subgroup ( Fig. 3b, Supplementary Table 3).

circRNA-miRNA co-expression analysis and miRNA-mRNA prediction
A predicted circRNA-miRNA network comprising 171 nodes and 150 edges was constructed based on the top 15 DECs. It revealed that a single circRNA could target multiple miRNAs and vice versa (Supplementary Figure  2). Certain miRNAs might contribute to the abnormal expression of mRNAs in PMOP patients (Fig. 4a-c). Predicted miRNA-mRNA interactions were also visualized ( Fig. 4d and e).

Construction of competing endogenous RNA (ceRNA) network
We integrated circRNA-miRNA and miRNA-mRNA interactions to construct a circRNA-miRNA-mRNA network, which provided preliminary insight into the association between the top 28 dysregulated circRNAs, 145 intermediate miRNAs, and 175 DEMs (Fig. 5).

Discussion
Due to the absence of 5′ caps and 3′ polyadenylated tails, circRNAs have been ignored in polyadenylated transcriptome studies for quite a long time. With the advent of high-throughput sequencing, computational biology, and biochemical methods, accumulating evidence has demonstrated the critical roles of circRNAs during the pathological process in various diseases [8,11,12]. Benefiting from covalently closed-loop structures, cir-cRNAs could serve as promising diagnostic and prognostic BTMs ascribed to their relative tolerance to exonucleases [6,11]. However, the role of circRNAs in PMOP and the mechanism by which circRNA regulates the balance between osteogenesis and osteoclastogenesis have not been well elucidated. Herein, we explored the DECs in PBMCs between PMOP patients and control individuals using whole transcriptome sequencing and further confirmed the expression level of these DECs by combining them with the public dataset. These findings are foundational for subsequent studies on circRNAs in PMOP and may offer insight into prevention and new treatment targets for PMOP. A consensus has been reached that ideal biomarkers for the diagnosis of PMOP and therapeutic drug monitoring should have the characteristics of abundance, easy-access, and less-invasiveness. Firstly, in the clinical setting, a blood sample is typically the starting point for biomarker search and discovery. The extraction of blood has minimal risks when performed by a trained technician. PBMCs are a great source of DNA for genetic analysis, in vitro culturing, and functional assays or for subsequent isolation of lymphocyte or monocyte subtypes, provided appropriate cryopreservation [16]. Secondly, there is growing evidence that morphological and functional changes in PBMCs, in most lymphocytes and monocytes, may reflect the severity of osteoporosis in postmenopausal women [17,18]. PBMCs are becoming a fairly common subject of research in the field of osteoporosis, and the ease of obtaining the material vs bone biopsy is an indisputable advantage [19,20]. In the present study, 373 DECs, including 73 upregulated and 300 downregulated circRNAs, were identified between PMOP patients and paired healthy individuals. Strikingly, 8 circRNAs (hsa_circ_0000471, hsa_circ_ 0008139, hsa_circ_0007385, hsa_circ_0008631, hsa_circ_ 0001824, hsa_circ_0027464, hsa_circ_0007507, hsa_circ_ 0000008) were significantly upregulated (|FC|>5) and 5 circRNAs (all novel) were significantly downregulated (|FC|>9), which implied that these DECs might exert crucial functions in regulating the pathogenesis and progression of PMOP. These DECs were mainly intronic and antisense circRNAs distributed among all chromosomes, with lengths mainly less than 300 bp. Then, GO annotation and KEGG pathway analyses were performed to elucidate the functions of DECs between PMOP and healthy individuals. The most significant GO items were G2 DNA damage checkpoint, transcription from RNA polymerase II promoter, NAD+ or protein binding, presynaptic membrane organization, and glutamate receptor activity, indicating that the coding genes contributed to the development of PMOP. In KEGG pathway analysis, several important pathways, including hsa04931 (insulin resistance), hsa04062 (chemokine signaling pathway), and has04530 (tight junction), were identified to take pivotal parts in the pathogenesis of PMOP. For instance, insulin resistance might affect osteoclast differentiation, activation, and survival via the tumor necrosis factorrelated cytokine receptor activator of nuclear factor kappa B ligand (RANKL)-induced pathway, which aggravated the osteoporosis process in postmenopausal women [21]. Moreover, activation of C-C chemokine receptor-2 (CCR2) is reported to be crucially involved in the signaling of nuclear factor-kappa B (NF-κB) and extracellular signal-related kinase 1 and 2 (ERK1/2), which contribute to the susceptibility to RANKL-induced osteoclastogenesis [22].
As highly conserved endogenous non-coding RNAs, numerous circRNAs harbor numerous miRNA binding sites, suggesting that they could sponge certain miRNAs and function as ceRNAs to regulate gene expression in several diseases [23][24][25][26]. circNT5E, a cir-cRNA formed from the NT5E genome and regulated by the RNA-editing enzyme, was reported to exert vital regulatory function by sponging glioblastoma suppressor miR-422a [27]. circNHSL1 was also reported to intervene in gastric cancer progression by sponging miR-1306-3p and could serve as a novel biomarker for early diagnosis of gastric cancer [28]. Similarly, highly expressed circ-SFMBT2 was observed in gastric cancer tissues and proved to participate in the pathogenesis of gastric cancer via miR-182-5p sponging [29]. To date, only a few pathogenic cir-cRNAs have been reported in osteoporosis and RNA sequencing tools and bioinformatics analysis have not been employed in these studies [30][31][32]. Herein, we constructed ceRNA networks with DECs and DEMs from whole transcriptome sequencing data, as well as intermediate miRNAs predicted through bioinformatics analysis. These data verified that circRNAs might exert regulatory functions in PMOP through miRNA response elements (MREs), which facilitated our understanding of the pathogenesis of PMOP.
An unavoidable limitation of this study was the relatively small sample size, which was ascribed to the strict inclusion and exclusion criteria. However, the demographic data of all samples were well balanced, which was beneficial to reduce system error from high-throughput sequencing and lower false-positive results. Nevertheless, a prospective study with a larger sample size and more validation techniques is ongoing.

Conclusions
The present study investigated potential circRNAmediated ceRNA interplays using sample-matched whole transcriptome profiles between PMOP patients and healthy individuals. This PMOP-specific dysregulated ceRNA network might provide a comprehensive understanding of ceRNA-mediated gene regulation in the pathogenesis of PMOP and lay a firm foundation for exploring promising diagnostic biomarkers and novel treatment targets for PMOP.