Circular RNA hsa_circ_0076690 acts as a prognostic biomarker in osteoporosis and regulates osteogenic differentiation of hBMSCs via sponging miR-152

CircRNA-seq analysis revealed that hsa_circ_0076690 and hsa_circ_0111433 were down-regulated in OP patients

The high quality clean reads were obtained from five pairs of OP and control samples by Illumina Hiseq 3000 and bioinformatics preprocessing. In total, 271 differentially expressed circRNAs were identified between OP and control samples, with |log₂FC| ≥ 1 and P-value < 0.05 (Supplementary Table 2). The expression profiles of abnormally expressed circRNAs were demonstrated by volcano plot in Figure 1A. Among these circRNAs, we noticed that hsa_circ_0076690 and hsa_circ_0111433 were significantly down-regulated in OP patients compared to control (fold change: -20.77 and -16.29, respectively). To explore the role of hsa_circ_0076690 and hsa_circ_0111433 during osteoporosis pathogenesis, we validated the expression of two circRNAs in 57 pairs of OP and control samples by qRT-PCR assay. As shown in Figure 1B, hsa_circ_0076690 expression was significantly decreased in OP group (P-value < 0.01) compared to control. The expression level of hsa_circ_0111433 in OP showed a similar trend with hsa_circ_0076690 (P-value < 0.01, Figure 1C).

Figure 1. (A) Differentially expressed circRNAs were demonstrated by volcano plot (red: up-regulation green: down-regulation). (B, C) The down-regulation of hsa_circ_0076690 and hsa_circ_0111433 was validated by qRT-PCR in 57 pairs of OP and control samples.

Hsa_circ_0076690 acts as a potential diagnostic biomarker for OP

To examine the clinical significance of hsa_circ_0076690 and hsa_circ_0111433 for OP, the Pearson correlation analysis was conducted to measure the correlation between circRNAs expression and clinical parameters of OP patients. The statistical results showed that both the expression of hsa_circ_0076690 and hsa_circ_0111433 was significantly correlated with bone mineral density (BMD) and T-score (P-value < 0.05, r < -0.7, respectively), while had no significant correlation with age (r = -0.032 and 0.198, respectively) and Body Mass Index (BMI, r = 0.064 and 0.130, respectively) (Table 1). Subsequently, ROC analysis was conducted to investigate the accuracy of the two circRNAs serving as biomarkers for OP. As shown in Figure 2A, the AUC value was 0.8299 (P-value < 0.001) with 79% sensitivity and 85% specificity for hsa_circ_0076690, while hsa_circ_0111433 did not show an acceptable diagnostic value in clinical samples (data not shown).

Figure 2. (A) ROC curve analysis of hsa_circ_0076690 in OP patients. (B) Predicted binding sites of hsa_circ_0076690 and miR-152. (C, D) Dual-luciferase reporter and RIP assay showed the direct binding of hsa_circ_0076690 and miR-152.

The potential regulatory axis of hsa_circ_0076690/miR-152/RUNX2

To investigate the molecular mechanism of hsa_circ_0076690 underlying osteoporosis pathogenesis, bioinformatics analyses was conducted. As predicted by miRanda and TargetScan software, hsa-miR-152-5p and hsa-miR-3678-3p were identified to have potential binding sites of hsa_circ_0076690 sequence (Figure 2B). Dual-luciferase reporter assay revealed that miR-152 could significantly inhibit the luciferase intensity in the hsa_circ_0076690 wild type compared to control and had no effect on hsa_circ_0076690 mutant (Figure 2C). Meanwhile, the luciferase intensity was not significantly reduced after co-transfection of miR-3678 and hsa_circ_0076690 wild type reporter (data not shown). Moreover, the RIP experiment indicated that hsa_circ_0076690 pulled down with anti-AGO2 was significantly enriched in miR-152 mimics group compared to control (P-value < 0.01, Figure 2D), indicating that hsa_circ_0076690 may function as a sponge for miR-152.

We further predicted the downstream targets of miR-152 by bioinformatics algorithm and found that RUNX2 was a potential target (Figure 3A). A subsequent dual-luciferase reporter assay revealed that the luciferase activity was repressed by miR-152 mimics in RUNX2 wild type, whereas the mutated luciferase reporter showed no changes (Figure 3B). In addition, qRT-PCR demonstrated that miR-152 expression was significantly increased in OP samples compared to control (Figure 3C), while RUNX2 expression was decreased (Figure 3D).

Figure 3. (A) Predicted binding sites of miR-152 and RUNX2. (B) Dual-luciferase reporter assay showed RUNX2 was a target of miR-152. (C, D) The up-regulation of miR-152 and down-regulation of RUNX2 were validated by qRT-PCR in 57 pairs of OP and control samples.

Hsa_circ_0076690 regulated osteogenic differentiation of hBMSCs through targeting miR-152

During the osteogenic differentiation of hBMSCs at 0, 7, 14 and 21 day, osteogenic-associated biomarkers were measured (Supplementary Figure 1A). We noticed that the expression of hsa_circ_0076690 and RUNX2 was accumulated (Figure 4A), on the contrary, miR-152 expression was gradually decreased in a time-dependent manner during osteogenic differentiation of hBMSCs (Figure 4A).

Figure 4. (A) The expression of hsa_circ_0076690, miR-152 and RUNX2 were measured by qRT-PCR during osteogenic differentiation of hBMSCs. (B) MiR-152 and RUNX2 expression levels were measured by qRT-PCR after transfected with hsa_circ_0076690 in hBMSCs. (C) RUNX2 expression was inhibited by miR-152 mimics and remained constant under miR-152/hsa_circ_0076690 co-overexpression. (D) The expression level of osteogenic-associated biomarkers and ALP activity were detected after transfected with hsa_circ_0076690, miR-152 mimics and miR-152/hsa_circ_0076690 co-overexpression in hBMSCs.

To further investigate the function of hsa_circ_0076690/miR-152/RUNX2 axis during osteogenic differentiation, qRT-PCR assays were performed in transfected hBMSCs (Supplementary Figure 1B). As a result, over-expression of hsa_circ_0076690 lead to the decreased expression of miR-152 and the increased expression of RUNX2 in hBMSCs, while si-hsa_circ_0076690 showed the opposite result (Figure 4B). In addition, the expression of RUNX2 was significantly inhibited by miR-152 mimics and remained constant under miR-152/hsa_circ_0076690 co-overexpression condition (Figure 4C). Moreover, the expression level of osteogenic-associated biomarkers and ALP activity were detected to measure the differentiation level of transfected hBMSCs. The results showed that osteogenesis differentiation was promoted by hsa_circ_0076690 over-expression and remain unchanged in miR-152/hsa_circ_0076690 co-overexpression condition (Figure 4D, Supplementary Figure 1C).

Study subjects and ethics statement

In total, 57 pairs of osteoporosis patients and non-osteoporosis healthy subjects were enrolled in this study. None of the subjects had metabolic diseases that affect bone metabolism or received medical treatment before. Clinical parameters were recorded from all the participants, including age, body mass index (BMI), BMD, T-score, β-crosslaps (β-CROSSL), N-terminal osteocalcin (OSTEOC) and Total Procollagen Type1 Intact N-terminal Propeptide (TPINP), etc (Supplementary Table 1). This study was approved by Ethics Committee of Qilu Hospital of Shandong University and was in accordance with Helsinki Declaration. The written informed consents were signed by all participants.

CircRNA-seq analysis

Five pairs of serum or plasma samples from osteoporosis patients and healthy controls were subjected to circRNA sequencing using Illumina Hiseq 3000 (Illumina, CA, USA). Firstly, total RNAs were extracted from samples using (TRIzol CA, USA) following the instruction. After measuring the quantity and quality of RNAs, Rnase R was used to remove linear transcripts. Subsequently, cDNA library was built for circRNA sequencing following the instruction of TruSeq RNA Sample Prep Kit (Illumina, CA, USA). FastQC and cutadapt software were used to obtain high quality clean reads for further analysis. GRCh38 genome and circBase were used for circRNAs annotation. EdgeR package was used to filter the differentially expressed circRNAs between osteoporosis and control group (log₂FC ≥ 1 and P-value < 0.01).

Bioinformatics and statistical analysis

Differentially expressed circRNAs between osteoporosis and control group was demonstrated by Volcano plot, which was drawn by ggplot2 package. The diagnostic values of candidate circRNAs were measured by ROC curve analysis and the areas under the curve (AUC) were calculated by R. Target prediction was performed by TargetScan and miRanda software. All experiments data are representative of independent triplicate experiments and the statistical analysis was performed using R. Statistically significant was set to P-value < 0.05.

Osteogenic differentiation of hBMSCs

The hBMSC cells were purchased from Ribio company (Beijing, China) cultured in culture medium supplemented with 12% fetal bovine serum (24-well plates, 37 °C). The hBMSCs induction process was conducted by supplying the cells with osteogenic differentiation medium and changed every three days. HBMSC cells were cultured for twenty-one days and the osteogenic associated biomarkers including ALP, Osteocalcin, β-catenin, BMP2 and RUNX2 were detected to measure the osteogenic differentiation process.

Cell transfection and ALP measurement

HBMSC cells were cultured in six-well plates and transfected with Lipofectamine 3000 for 6h according to the manufacturer’s instruction. The negative controls (NC), miR-152 mimics and si-hsa_circ_0076690 were purchased from Genechem company (Shanghai, China). The ALP activity was measured after transfected for 14 days using a commercial kit according to the instruction. Absorption was measured at 405 nm.

Quantitative real time-polymerase chain reaction (qRT-PCR)

Total RNAs were extracted from osteoporosis samples, control samples and transfected hBMSCs using TRIzol reagent (Invitrogen, CA, USA) according to the protocol. The reverse-transcription reaction was conducted using Reverse Transcriptase kit (Invitrogen, CA, USA). The Real-time PCR was performed to detect the expression of candidate circRNAs in tissue samples and specific molecule in transfected hBMSC cells. A SYBR Green PCR kit (Takara, Japan) was used following the instrument: denaturation for 180 s at 93 °C and followed by 30 cycles at a Real-Time PCR Applied Biosystems. The relative expression level was determined by 2^-ΔΔCt method, GAPDH and U6 were used as references.

Western blot and Alizarin red staining (ARS)

The hBMSC cells from different groups were lysed with protein-extraction reagent (Roche, Mannheim, Germany) following the instrument. The BCA Protein Assay Kit (Beyotime, China) was used to measure protein concentration. Then proteins were separated by 12% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to PVDF membrane to incubate with primary antibody. Finally, the membrane was washed (three times) and further incubated with secondary antibody. The bands were visualized by ECL (Millipore) following the instrument. ARS staining was performed at 14d and 21d during osteogenic differentiation by treating cells with 0.1% ARS-Tris-Hcl (1h, 37°C), then the mixture was washed by distilled water and ARS level was examined by microscope.

Dual-luciferase reporter assay and RNA immune-precipitation

The potential binding sequences between RNA molecules were synthesized and cloned into pGL3-Basic luciferase vector. Luciferase plasmid and miR-152 NC/mimics were transiently transfected into cells with Lipofectamine 3000 for 24h following the protocol. The luciferase activities were then measured using an applied system (Promega, USA). The RNA immunoprecipitation (RIP) experiment was conducted in transfected cells using an RNA Binding Protein Kit (Millipore, USA) and human anti-Ago2 antibody following the protocols. IgG was used as negative control and binding output was measured by qRT-PCR.