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

Abstract

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.

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 microRNA (miRNA) sponges, serving as scaffolds in the assembly of protein complexes, regulating the expression of parental genes, and modulating alternate splicing, as well as RNA-protein 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 closed-loop 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 circRNA 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 circRNAs (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.

Materials and methods

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 real-time quantitative polymerase chain reaction (qRT-PCR). The demographic and clinical characteristics of the six participants are summarized in Table 1.

Table 1 Demographic characteristics of the participants

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.

Whole transcriptome sequencing

RNA high throughput sequencing was performed Cloud-Seq Biotech Ltd. Co. (Shanghai, China). Briefly, total RNA was used after removing ribosomal RNAs (rRNAs) with NEBNext rRNA Depletion Kit (New England Biolabs, Inc., Massachusetts, USA) and constructing RNA libraries with NEBNext® Ultra™ II Directional RNA Library Prep Kit (New England Biolabs, Inc., Massachusetts, USA). Libraries were controlled for quality and quantified using the BioAnalyzer 2100 system (Agilent Technologies, Inc., USA). Libraries (10 pM) were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters and finally sequenced for 150 cycles on the Illumina HiSeq sequencer (HiSeq 4000; Illumina, Inc.).

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₂(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 pathway-enrichment analysis was also performed according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Target miRNA prediction

The circMic software (http://www.bioinf.com.cn/) was referred to predict specific microRNAs (miRNAs) of corresponding DECs. Interactions between miRNAs and circRNAs were evaluated with miRanda (http://www.microRNA.org/). The interaction between dysregulated mRNA and miRNA were predicted with TargetScan release 7.0 (http://www.targetscan.org/) and presented via miRNet 2.0 platform (https://www.mirnet.ca/).

circRNA-miRNA-mRNA network prediction

The overlapping miRNAs between the predicted circRNAs 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.

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.

Result

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 circRNAs 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).

Fig. 1

Overview of circRNAs in the PMOP and the control groups. (a) Expressions were normal in both the PMOP group (samples 1-3) and the control group (samples 4-6). (b) Length distribution of all circRNAs. (c) circRNAs were unevenly located in all chromosomes. (d) circRNAs were categorized into five types

Table 2 Overview of sequencing and quality control of samples

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 sense-overlapping 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.

Fig. 2

Comparison of differentially expressed circRNAs (DECs) between the PMOP and the control groups. (a) Distribution of DECs was exhibited via volcano plot map and compared between the PMOP group (the former three columns) and control group (the latter three columns). (b) Each row represented one circRNA and each column represented each sample. (c) DECs with different functions were exhibited

Table 3 TOP 15 upregulated and downregulated circRNAs in PMOP patients

qRT-PCR validation of DECs

To verify the reliability of sequencing data, four upregulated circRNAs (circ_0000471, circ_0008139, circ_0001824, circ_0008345) and four downregulated circRNAs (circ_0112054, circ_0000443, circ_0077548, circ_0001395) were randomly selected. The expression levels of these dysregulated circRNAs were following the sequencing data (Supplementary Figure 1).

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 circRNAs 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). Furthermore, the top 20 KEGG pathways in the dysregulated circRNAs are shown in Fig. 3c and d (Supplementary Tables 4 and 5). The results revealed that significantly upregulated circRNAs were mainly associated with chemokine signaling pathway (hsa04062), bacterial invasion of epithelial cells (hsa05100), and RAS signaling pathway (hsa04014), while significantly downregulated circRNAs were mainly associated with insulin resistance (hsa04931), tight junction (hsa04530), and transcriptional misregulation (hsa05202).

Fig. 3

GO annotation and KEGG pathway analysis for DECs. GO analysis results of (a) upregulated and (b) downregulated circRNAs in terms of “biological process,” “cellular component,” and “molecular function.” TOP 10 pathway analysis results of (c) upregulated and (d) downregulated circRNAs based on the KEGG database

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).

Fig. 4

Differentially expressed mRNAs (DEMs) and prediction of miRNA-mRNA interactions. DEMs between the PMOP and the control groups were assessed and exhibited via (a) heat map, (b) volcano plot map, and (c) scatter plot analysis. Predicted interactions between (d) miRNAs (blue square) and upregulated mRNAs (red circle) and (e) downregulated mRNAs

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).

Fig. 5

The ceRNA network was constructed based on circRNA-miRNA-mRNA interactions. Square, triangle, and circle nodes represent circRNA, miRNA, and mRNA, respectively. Red color indicates upregulation and blue color indicates downregulation

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, circRNAs 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 sub-types, 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 factor-related 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 circRNA 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 circRNAs 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 circRNA-mediated 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.