Longitudinal Changes of Circulating miRNAs During Bisphosphonate and Teriparatide Treatment in an Animal Model of Postmenopausal Osteoporosis.

MicroRNAs regulate bone homeostasis, and circulating microRNAs have been proposed as novel bone biomarkers. The effect of anti-osteoporotic treatment on circulating microRNAs has not been described in detail. Therefore, we performed a comprehensive analysis of microRNA serum levels in ovariectomized (OVX) and sham-operated (SHAM) rats over 12 weeks of antiresorptive or osteoanabolic treatment. Forty-two Sprague Dawley rats underwent SHAM surgery (n = 10) or ovariectomy (n = 32). After 8 weeks, OVX rats were randomized to antiresorptive treatment with zoledronate (n = 11), osteoanabolic treatment with teriparatide (n = 11), or vehicle treatment (n = 10). Serum samples were collected at weeks 8, 12, 16, and 20 after surgery. A total of 91 microRNAs were analyzed by RT-qPCR in serum samples collected at week 20. Based on the results, 29 microRNAs were selected for longitudinal analysis at all four study time points. Changes in bone mineral density and microstructure were followed up by in vivo micro-CT and ex vivo nano-CT. Ovariectomy resulted in the loss of trabecular bone, which was reversed by osteoanabolic and antiresorptive treatment. Differential expression analysis identified 11 circulating miRNAs that were significantly regulated after treatment. For example, miR-107 and miR-31-5p increased in vehicle-treated OVX animals, whereas they decreased during teriparatide treatment. Additional miRNAs were identified that showed significant correlations to bone microstructure or bone miRNA expression, including miR-203a-3p, which exhibited a significant negative correlation to vertebral and tibial trabecular bone volume fraction (%). Longitudinal analysis confirmed eight microRNAs with significant changes in serum over time that were prevented by teriparatide and zoledronate treatment (miR-34a-5p, miR-31-5p, miR-30d-3p, miR-378a-5p) or teriparatide treatment only (miR-375-3p, miR-183-5p, miR-203a-3p, miR-203b-3p). Gene target network analysis identified WNT and Notch signaling as the main signaling pathways controlled by these miRNAs. Thus, ovariectomy results in time-dependent deregulation of circulating miRNAs compared with SHAM animals. Anti-osteoporotic treatments can rescue this effect, showing that bone-related miRNAs might act as novel biomarkers for treatment monitoring.

Introduction

Osteoporosis is a worldwide health concern affecting more than 200 million people.( ) Osteoporosis‐related fractures are known to predict further fracture risk.( ) Moreover, increased morbidity, mortality, impaired quality of life, and vast health care cost were reported repeatedly.( ) Anti‐osteoporotic treatments are used to decrease fracture risk in osteoporosis patients by controlling bone turnover and improving bone mineral density, microstructure, geometry, and material properties. In recent years, options for clinicians regarding anti‐osteoporotic compounds have increased. Although bisphosphonates such as zoledronic acid (ZOL) are still recommended as first‐line antiresorptive therapy,( ) teriparatide (TPD), a recombinant human parathyroid hormone‐fragment (rhPTH1‐34), has received increased attention and is currently especially considered for osteoporosis patients at high fracture risk.( ) It has to be noticed that the mechanisms of action are considerably different between antiresorptive and osteoanabolic agents. Bisphosphonates such as ZOL inhibit osteoclasts and thereby bone resorption mainly via the mevalonate pathway.( ) On the other hand, osteoanabolic TPD stimulates osteoblasts and, thus, bone formation. Numerous mechanisms including protein kinases, MAP‐kinase, phospholipases, as well as the WNT‐signaling pathway are involved in this process.( ) However, the exact mechanism of TPD is not fully understood yet. The follow‐up and treatment response of anti‐osteoporotic agents are usually measured by dual‐energy X‐ray absorptiometry (DXA) analyses and established bone turnover markers. However, bone mineral density (BMD) changes only account for 4% to 74% of the fracture risk reduction that can be achieved under anti‐osteoporotic treatment,( ) and repeated, short‐term DXA measurements are not recommended. Bone turnover biomarkers such as C‐terminal cross‐linking telopeptide of type I collagen (CTX) or procollagen type 1 N‐propeptide (PINP) were reported to predict BMD changes( ) and estimate fracture risk( ) at cohort level. These established bone biomarkers can monitor adherence and partial efficacy of therapy at an individual patient level but do not give exact information on the pathophysiology of disease.

MicroRNAs (miRNAs) are a new class of biomarkers. miRNAs are small, non‐coding RNAs, which act as posttranscriptional regulators of gene expression. Several studies have focused on miRNAs in bone metabolism. Significant differences in serum miRNA levels were found in patients with osteoporosis when compared with healthy controls,( ) as well as in recent( ) and manifest osteoporotic fractures.( ) Patients with vertebral fractures show elevated miRNA serum levels when compared with patients without fractures and healthy controls.( ) Thus, there is evidence for the association of circulating miRNAs to various bone diseases, but only a few studies have examined the impact of anti‐osteoporotic treatment on miRNA levels. For example, the influence of alendronate, an oral bisphosphonate, on miRNA expression in bone cells was investigated. A significant upregulation of miR‐182‐5p attributable to ovariectomy (OVX) was found in femoral bone tissue and alendronate treatment resulted in a reduction of miR‐182‐5p tissue levels, which was found to be a regulator of osteoblast differentiation and apoptosis.( ) This raises the question whether antiresorptive agents such as ALN could directly or indirectly influence osteoblast function through miRNA regulation in humans. The role of alendronate on the differential expression of miRNAs, their target genes, as well as osteoclastogenesis was further determined by Jie and colleagues.( ) Their data suggest that bisphosphonates regulate the miR‐101‐3p/Rap1b signaling pathway and thereby osteoclast differentiation.( ) On the other hand, in vitro data suggest TPD induces osteogenic differentiation via regulation of the miR‐375/Runx2 axis in human marrow mesenchymal cells (hMSCs).( ) Cross‐sectional comparisons of miRNAs in serum of postmenopausal osteoporosis patients reported changes after 3 months of treatment with recombinant parathyroid hormone (TPD)( ) and after sequential treatment with TPD or ZOL followed by denosumab.( ) However, the changes of bone‐related miRNAs during anti‐osteoporotic treatment in a longitudinal manner is presently not sufficiently understood.

Recently, we have initiated the TAMIBAT Study (Time‐dependent Analysis of microRNAs and Bone Microstructure under Consideration of Anti‐Osteoporotic Treatment) to investigate miRNA transcription in femoral head bone tissue after ZOL and TPD treatment, respectively.( ) We observed both common and distinct effects of antiresorptive and osteoanabolic treatments on bone tissue miRNA levels: first, TPD treatment had a stronger impact on miRNA transcription in bone (46 differentially regulated miRNAs) compared with ZOL treatment (10 differentially regulated miRNAs). Interestingly, the two anti‐osteoporotic treatments partially reversed distinct effects of OVX‐related changes in miRNA transcription. For example, treatment with ZOL but not TPD resulted in downregulation of OVX‐induced miR‐205, whereas both treatments led to a reduction of miR‐203b that was also upregulated due to OVX.

Conversely, miR‐203a, which is a potent inhibitor of osteogenic differentiation,( ) was observed to be upregulated during OVX‐induced bone loss and downregulated only by TPD but not ZOL treatment. Interestingly, the levels of miR‐203a in bone and serum were significantly correlated, which provides evidence for potential application of circulating miRNA biomarkers for monitoring anti‐osteoporotic treatment.

Our initial analysis in the TAMIBAT study had focused on the regulation of miRNA transcription in bone tissue in response to antiresorptive and osteoanabolic treatment. In the next step, we aimed to investigate the utility of miRNAs as bone biomarkers by systematically analyzing the impact of anti‐osteoporotic therapies on the circulating (serum) levels of miRNAs. Hence, the primary objective of this study was to evaluate the time‐dependent changes of circulating miRNAs in serum of ovariectomized rats during ZOL and TPD treatment, respectively. Secondary objectives included the longitudinal changes in bone microstructure, assessed in vivo and ex vivo, as well as the association of miRNAs to parameters of trabecular bone microstructure.

Materials and Methods

The TAMIBAT study was conducted at the Ludwig Boltzmann Institute for Traumatology in cooperation with the Ludwig Boltzmann Institute of Osteology, the Medical University of Vienna, and TAmiRNA GmbH (all Vienna, Austria).

Study design

We conducted an animal‐study in OVX rats, an established model for postmenopausal osteoporosis, and tested the effects of OVX as well as osteoanabolic and antiresorptive treatment on serum miRNA levels and their association to bone microstructure. The study duration was 20 weeks in total with a treatment period of 12 weeks for all groups. Six‐month‐old Sprague Dawley rats (n = 42) were included in the study. Ovariectomy was performed in 32 rats in order to reach postmenopausal state (OVX). Ten animals underwent sham‐operation (SHAM). After 8 weeks, animals in the OVX group were randomized either to vehicle (VEH; n = 10), TPD (n = 11), or ZOL (n = 11) therapy. A single dose of zoledronic acid 100 μg/kg body weight s.c. was injected into animals in the ZOL group (Novartis Pharma AG, Basel, Switzerland). In the TPD group, the agent was administered s.c. periodically (weekly dose 210 μg/kg body weight [42 μg/kg body weight per day over 5 days], Eli Lilly and Company, Indianapolis, IN, USA). Animals in the VEH and SHAM group received vehicle injections (0.9% saline with 20 mM NaH2PO4 in 0.9% NaCl, 3 mg/mL mannitol, dosed at a volume of 1 mL/kg body weight) s.c. at the same interval as TPD, and thus acted as control groups.

For the present study, blood samples were collected at weeks 8, 12, 16, and 20 and analyzed. Blood sampling was performed from the rat's lateral tail vein under short inhalation anesthesia (isoflurane 1% to 3%, oxygen). To avoid an impact of food intake on bone metabolism, blood was taken in the early afternoon after a 6‐hour fasting interval. Clotting time for the serum samples was 30 to 45 minutes at room temperature. Afterward, samples were immediately centrifuged at 2500g for 15 minutes at 20°C and subsequently frozen at –80°C.

Bone microstructure analysis in vivo was performed at weeks 8 and 12 using micro‐CT. At week 20 (end of the study), ex vivo scans were conducted using nano‐CT. The study design is shown in Fig. 1.

Fig 1

Study design of the TAMIBAT study. (A) Forty‐two Sprague Dawley rats were ovariectomized (OVX; n = 32) or underwent sham‐surgery (SHAM; n = 10) at the age of 6 months. After 8 weeks, rats in the OVX group were randomized either to placebo (vehicle solution [VEH]; n = 10), teriparatide (TPD; n = 11), or zoledronic acid (ZOL; n = 11) therapy. Animals in the SHAM group received vehicle solution (SHAM; n = 10). Serum for miRNA measurements was taken at weeks 8, 12, 16, and 20. In vivo micro‐CT analyses were performed at week 8 and week 12. At the end of study (week 20), ex vivo nano‐CT was performed in all animals. (B) The first report( ) on the TAMIBAT study was focused on cross‐sectional characterization of bone tissue microRNA expression, whereas the present report systematically assessed the levels of bone‐related circulating microRNAs in serum. First, a panel of 91 microRNAs was analyzed cross‐sectionally in serum at the study endpoint (week 20). Based on differential expression analysis and correlation analysis using NanoCT and NGS data, a focused set of 29 microRNAs was selected for longitudinal serum analysis throughout the whole study.

Animals

In total, 42 six‐month‐old female Sprague Dawley rats were included in the study. Thirty‐two of them were ovariectomized and 10 sham‐operated. All animals were sourced by Charles River Laboratories (Gottingen, Germany). Two rats were kept per cage (makrolon type 4 cages). Ready‐made feed (Rat Maintenance, Ssniff GmbH, Soest, Germany) and water were fed ad libitum. The number of animals per group was determined according to previously published studies,( , ) with a difference in the group average of 15% (90% power). Bedding was obtained from Abedd GmbH (Kalnciems, Latvia; Midi Chips). Gnawing blocks and tunnels were provided as cage enrichment.

Before all study‐related procedures, the study design and the experimental protocols were approved by the Municipal Government of Vienna in accordance with Austrian law and the Guide for the Care and Use of Laboratory Animals as defined by the National Institute of Health (revised 2011).

Bone microstructure assessment

Assessment of trabecular bone microstructure as well as BMD in vivo was performed at weeks 8 and 12 at the vertebral body L4 (SCANCO vivaCT 75, SCANCO Medical AG, Brüttisellen, Switzerland). In vivo scans were operated at 114 μA, 70 kVp, 500 projections/180° with an integration time of 300 ms and reconstructed to an isotropic resolution of 45 μm. For in vivo analyses, the trabecular bone volume fraction (BV/TV, %) as well as BMD (mgHA/cm3) were chosen to demonstrate the effects of OVX, sham‐operation, as well as TPD, ZOL, and placebo treatment. Measurements were performed under short inhalation anesthesia, using the same anesthesia protocol as for blood draw.

Ex vivo scans were assessed at week 20 (end of the study) using a high‐resolution μCT imaging system (nano‐CT, μCT 50, SCANCO Medical AG). Scans were performed at 90 kVp, 200 μA, 1000 projections/180° with an integration time of 500 ms and an isotropic nominal resolution of 10 μm.

The tibia and the lumbar spine (L4) were used for the evaluation of trabecular bone microstructure. Trabecular microstructure parameters included the BV/TV (%), the mean trabecular number (Tb.N, mm−1), the mean trabecular thickness (Tb.Th, mean, μm), and mean trabecular spacing (Tb.Sp, mean, μm).

Image segmentation and analyses were done by Definiens Developer XD 2.1.1 (Definiens AG, Munich, Germany), Fiji, and BoneJ plugin, respectively. The procedure is described in detail elsewhere.( )

MicroRNA analysis serum

RNA extraction

Total RNA was extracted from 200 μL serum using the miRNeasy Mini Kit (Qiagen, Hilden, Germany). Samples were thawed on ice and centrifuged at 12,000g for 5 minutes to remove any cellular debris. For each sample, 200 μL of serum were mixed with 1000 μL Qiazol and 1 μL of a mix of three synthetic spike‐in controls (Qiagen). After a 10‐minute incubation at room temperature, 200 μL chloroform were added to the lysates followed by cooled centrifugation at 12,000g for 15 minutes at 4°C. Precisely 650 μL of the upper aqueous phase were mixed with 7 μL glycogen (50 mg/mL) to enhance precipitation. Samples were transferred to a miRNeasy mini column where RNA was precipitated with 750 μL ethanol followed by automated washing with RPE and RWT buffer in a QiaCube liquid handling robot. Finally, total RNA was eluted in 30 μL nuclease‐free water and stored at −80°C to await further analysis.

MicroRNA reverse‐transcription quantitative PCR analysis in serum RNA (RT‐qPCR)

Starting from total RNA samples, cDNA was synthesized using the miRCURY LNA RT kit (Qiagen). Reaction conditions were set in accordance with the manufacturer's specifications. In total, 2 μL of total RNA were used per 10 μL reverse transcription (RT) reaction. To monitor RT efficiency and presence of impurities with inhibitory activity, a synthetic RNA spike‐in (cel‐miR‐39‐3p) was added to the RT reaction. PCR amplification was performed in a 384‐well plate format in a Roche LC480 II instrument (Roche, Mannheim, Germany) using miRCURY LNA SYBR Green PCR kit (Qiagen) with the following settings: 95°C for 10 minutes, 45 cycles of 95°C for 10 seconds and 60°C for 60 seconds, followed by melting curve analysis. To calculate the cycle of quantification values (Cq‐values), the second derivative method was used. Spike‐in control values were used for monitoring data quality. Serum samples with low RNA recovery assessed by spike‐in levels were not included in the analysis (5 of 42 serum samples were excluded based on poor spike‐in signals) (Supplemental Fig. S1). Endogenous miRNAs with an average raw Cq‐value of <35 were included in the analysis (80 of 91 microRNAs met this criterium) and were normalized to the RNA spike‐in controls by subtracting the individual miRNA Cq‐value from the RNA Spike‐In Cq‐value, thus obtaining delta‐Cq (dCq) values that were used for the analysis. Hemolysis was assessed in all samples using the ratio of miR‐23a‐3p versus miR‐451a and applying a cut‐off of a delta Cq > 7 indicates a high risk of hemolysis.( )

Statistical analysis

Statistical analysis was performed using GraphPad Prism version 8.0.0 (GraphPad Software, San Diego, CA, USA). Treatment and OVX‐related changes of BMD and bone microstructure analyzed via micro‐CT between week 8 and week 12 were analyzed by two‐way ANOVA in conjunction with Bonferroni post hoc test. Statistical differences in bone microstructure, which were analyzed via nano‐CT at week 20, were assessed using one‐way ANOVA followed by Tukey's multiple comparison post hoc test.

For longitudinal analysis of microRNA levels in serum, two‐way ANOVA followed by Tukey's multiple comparison post hoc test was performed. Correlation analysis of serum microRNA levels with different bone microstructure parameters (n = 37) and bone microRNA levels (n = 15) was conducted using linear regression analysis and calculation of Pearson correlation coefficients. For correlation with serum qPCR data, bone nano‐CT data and bone NGS data were log2 transformed.

Data in all figures were visualized using box plots with individual data points indicated, for longitudinal analysis of serum microRNA levels mean values between time points were connected.

Multiple comparison adjusted p values reaching p < .05 were indicated in all figures.

microRNA target network construction

A total of 12 microRNAs identified as significantly regulated (Table 5) were used for constructing a target network using the online tool miRnet 2.0 (accessible via mirnet.ca).( ) Genes listed in miRTarBase v8.0 were selected for network construction and the degree filter was set to 4; hence, only target nodes with at least four connections remained in the network. The reactome database( ) was used for pathway enrichment, using all genes identified in the network and hypergeometrical testing.

Table 5

Results of Two‐Way ANOVA on Longitudinal Data Set

microRNA	Two‐way ANOVA (mixed effects analysis)
	Time point	Treatment	Time point × Treatment
miR‐103a‐30	p = .013
miR‐126‐5p	p = .001
miR‐183‐5p		p = .031	p = .007
miR‐203a‐3p			p = .010
miR‐30d‐3p	p = .001
miR‐34a‐5p	p = .016
miR‐375	p = .015	p = .004
miR‐378a‐5p	p = .005
miR‐31‐5p	p = .003
miR‐203b‐3p	p = .019
miR‐455‐3p	p < .0001		p = .007
miR‐582‐5p	p = .036

Two‐way ANOVA–derived p values reaching a threshold of p < .05 for the factor timepoint, treatment, and the interaction effect between timepoint and treatment are listed.

Results

Study design for TAMIBAT

We have previously described the impact of OVX treatment and therapy on tissue levels of miRNAs in femoral head bone tissue.( ) Here, we intended to systematically characterize the circulating levels of bone‐related miRNAs in serum at multiple time points. First, a cross‐sectional analysis of 91 miRNAs was performed in serum samples collected at the endpoint (week 20) of the study in order to identify miRNAs with highest differences between treated and untreated OVX rats. Based on these data, we narrowed the number of miRNAs to 29, which were analyzed in all 128 serum samples collected at weeks 8, 12, 16, and 20 (Fig. 1).

BMD and BV/TV change during the early phase of anti‐osteoporotic treatment (week 8 to week 12)

To confirm appropriate response in the animals to the treatments, in vivo micro‐CT analyses were performed at week 8 and week 12. The analysis showed a significant increase in BMD from week 8 to week 12 in TPD (p < .0001), while BMD for SHAM (p > .99), VEH (p > .99), and ZOL (p = .48) did not change over time (Fig. 2).

Fig 2

Changes in bone mineral density (BMD; mgHA/cm3) and ratio of bone volume to tissue volume (BV/TV; %) analyzed by micro‐CT between week 8 and week 12 at the tibia. Significant increases in (A) BMD and (B) BV/TV were observed in the teriparatide‐treated group (TPD). In the OVX‐placebo group (VEH), BV/TV but not BMD significantly declined over time. No differences between weeks 8 and 12 were found in sham‐operated (SHAM) and zoledronic acid (ZOL)‐treated rats. Multiple comparisons Bonferroni corrected p values reaching a threshold of p < .05 for the contrast week 12 versus week 8 in each group are indicated.

The BV/TV (%) increased significantly between week 8 and week 12 in TPD‐treated animals (p < .0001). A significant reduction in BV/TV was found in the VEH group (p = .037). BV/TV remained stable in SHAM and ZOL animals (p > .99 and p = .74, respectively) (Fig. 2).

Bone microstructure analysis at study endpoint (week 20)

Trabecular bone microstructure was assessed by ex vivo nano‐CT at the study endpoint (week 20) at two sites, the lumbar spine (L4) and tibia (Fig. 3). At the lumbar spine, BV/TV was lowest in the VEH group. ZOL resulted in a partial rescue of BV/TV loss compared with the VEH group, but values in ZOL remained significantly lower than in SHAM. TPD‐treated animals showed the highest BV/TV values with significant differences to SHAM, ZOL, and VEH. The high BV/TV values in TPD‐treated animals corresponded to high Tb.Th but not Tb.N. ZOL therapy did not influence Tb.N or Tb.Th. However, a lower Tb.Sp was observed in ZOL in comparison to VEH‐treated rats.

Fig 3

Trabecular bone microstructure analyzed by nano‐CT at the study end at 20 weeks, corresponding to 12 weeks of treatment. Trabecular bone microstructure assessed ex vivo by nano‐CT (A) at the tibia and (B) the 4th vertebral body. Figures indicate the ratio of bone volume to tissue volume (BV/TV; %), trabecular number (Tb.N; 1/mm), trabecular thickness (Tb.Th; mm), trabecular separation (Tb.Sp; mm) for the TAMIBAT rat cohort (SHAM, n = 10; VEH, n = 10, TPD, n = 11; ZOL, n = 11). Tukey's multiple comparisons adjusted derived p values reaching a threshold of p < .05 are indicated.

At the tibia, TPD, but not ZOL, restored BV/TV when compared with SHAM (Fig. 3). As observed at the lumbar spine, (i) Tb.N was decreased in all OVX groups, independently of type of treatment, (ii) TPD resulted in highest Tb.Th values, and (iii) VEH in highest Tb.Sp.

Bone loss and anti‐osteoporotic therapy induces cross‐sectional changes in miRNA serum levels at study endpoint (week 20)

A panel of 91 miRNAs was designed based on our previously reported small RNA sequencing data in femoral head bone tissue,( ) with the rational to focus on bone/serum correlative miRNAs. We analyzed this panel of miRNAs in serum samples from all rats included in the TAMIBAT study (n = 42) at the study endpoint (week 20). After filtering of low‐abundant miRNAs and hemolytic samples, 80 miRNAs and 37 samples were included in the analysis. Threshold for p value (p < .05) and log2 fold change > |1| were set to identify deregulated miRNAs (Fig. 4, Table 1).

Fig 4

Cross‐sectional serum microRNA analysis at the study end (week 20). (A–C) Volcano plots showing log2 transformed fold change (x axis) and p values (y axis) for 80 microRNAs included in the analysis and three group contrasts (SHAM, n = 9; VEH, n = 9; TPD, n = 9; ZOL, n = 10).

Table 1

Cross‐Sectional Regulation of Circulating miRNAs in Ovariectomized Animals and After 12 Weeks of Treatment With Teriparatide (TPD) and Zoledronate (ZOL) at the Study End (Week 20)

microRNA	VEH (n = 9) versus SHAM (n = 9)		TPD (n = 9) versus VEH (n = 9)		ZOL (n = 10) versus VEH (n = 9)
	Log2FC	p Value	Log2FC	p Value	Log2FC	p Value
miR‐31‐59	1.07	.031	–0.98	.046	–0.30	.556
miR‐107	1.02	.098	–1.06	.242	–0.86	.086
miR‐375	–0.41	.119	0.16	.532	0.32	.281
miR‐203a‐3p	0.64	.123	–0.59	.148	0.01	.986
miR‐34a‐5p	0.95	.143	–0.91	.270	–0.35	.610
miR‐205	0.73	.228	–0.27	.566	–0.33	.531
miR‐21‐5p	0.32	.235	–0.19	.555	0.06	.832
miR‐103a‐3p	0.74	.24	–0.42	.52	–0.03	.96
miR‐363‐3p	0.61	.25	–0.35	.57	–0.43	.38
miR‐181d‐5p	0.89	.26	–0.55	.55	–0.34	.63
miR‐378a‐5p	0.67	.285	–1.26	.062	0.08	.889
miR‐708‐5p	0.70	.321	0.18	.790	0.18	.788
miR‐30d‐3p	0.42	.340	–0.50	.319	0.08	.842
Let‐7b‐5p	0.40	.38	–0.14	.69	0.15	.73
miR‐126‐5p	0.17	.386	–0.01	.979	0.28	.226
miR‐324‐5p	0.46	.410	–0.12	.849	–0.07	.889
miR‐466c‐5p	0.41	.453	–0.69	.454	–0.33	.520
miR‐582‐5p	–0.28	.543	0.54	.203	1.18	.010
miR‐214‐3p	–0.26	.568	–0.03	.929	0.88	.049
miR‐188‐5p	0.30	.572	–1.22	.192	0.16	.777
miR‐152‐3p	–0.25	.589	–0.26	.528	0.65	.102
miR‐139‐5p	0.09	.746	0.17	.595	0.38	.209
miR‐183‐5p	0.13	.760	–0.36	.330	0.29	.380
miR‐455‐3p	–0.16	.805	0.07	.903	1.13	.064
miR‐335‐5p	–0.23	.812	–0.10	.912	1.82	.056
miR‐143‐3p	–0.13	.834	0.34	.557	1.29	.024
miR‐29b‐3p	0.04	.93	–0.11	.81	0.62	.2
miR‐203b‐3p	0.03	.967	–1.73	.059	0.54	.374
miR‐204‐5p	0.00	.995	–0.14	.814	1.15	.017

p < .1 are indicated in bold.

In vehicle‐treated OVX animals miR‐31‐5p (log2FC = 1.07, p < .05) and miR‐107 (log2FC = 1.02, p < .1) were found to be upregulated in comparison to SHAM controls. Treatment with teriparatide was able to reverse both effects and resulted in a downregulation of miR‐31‐5p (log2FC = −0.98, p < .05) and miR‐107 (log2FC = −1.06, p = .24) in TPD versus VEH.

TPD treatment furthermore resulted in the downregulation of miR‐188‐5p (log2FC = −1.22, p = .19), miR‐378a‐5p (log2FC = −1.26, p < .1), and miR‐203b‐3p (log2FC = −1.73, p < .1) in comparison to VEH.

OVX animals undergoing treatment with zoledronate showed higher serum levels for miR‐582‐5p (log2FC = 1.18, p < .05), miR‐143‐3p (log2FC = 1.29, p < .05), miR‐204‐5p (log2FC = 1.15, p < .05), miR‐322‐5p (log2FC = 0.60, p < .05), miR‐214‐3p (log2FC = 0.88, p < .05), miR‐335‐5p (log2FC = 1.82, p < .05), and miR‐455‐3p (log2FC = 1.13, p < .1) in comparison to the vehicle control.

miRNA serum levels are correlated to bone microstructure parameters

miRNA serum levels were correlated to trabecular bone microstructure parameters, assessed by nano‐CT at week 20. At the vertebral body (Table 2), significant correlations between miR‐31‐5p and BV/TV (r = −0.35, p = .036) as well as miR‐378a‐5p and Tb.Th (r = −0.38, p = −.020) were found. Tb.N was significantly correlated to miR‐708‐5p (r = −0.38, p = .030) and miR‐139‐5p (r = −0.36, p = .030). The most significant findings were observed for miR‐203 family members: miR‐203a‐3p serum levels correlated negatively to BV/TV (r = −0.34, p = .040), Tb.N (r = −0.37, p = .027), Tb.Th (r = −0.33, p = .049), and positively (by trend) to Tb.Sp (r = 0.30, p = .076). miR‐203b‐3p levels were negatively correlated to BV/TV (r = −0.36, p = .028) and Tb.Th (r = −0.50, p = .002).

Table 2

Correlation Analysis (Pearson) of Serum microRNA Levels at the Study End (Week 20) With Lumbar Spine (L4) Bone Microarchitecture

L4	BV/TV		Tb.Th		Tb.N		Tb.Sp
microRNA	PCC	p Value	PCC	p Value	PCC	p Value	PCC	p Value
miR‐139‐5p	0.03	.883	–0.03	.877	–0.36	.030	–0.02	.920
miR‐203a‐3p	–0.34	.040	–0.33	.049	–0.37	.027	0.30	.076
miR‐203b‐3p	–0.36	.028	–0.50	.002	–0.16	.362	0.22	.183
miR‐31‐5p	–0.35	.036	–0.0.30	.067	–0.21	.207	0.30	.070
miR‐378a‐5p	–0.30	.068	–0.38	.020	–0.04	.829	0.20	.243
miR‐708‐5p	–0.15	.431	–0.15	.408	–0.38	.030	0.12	.524

BV/TV = trabecular bone volume fraction; Tb.Th = trabecular thickness; Tb.N = trabecular number; Tb.Sp = trabecular spacing; PCC = Pearson correlation coefficient.

N = 37; p < 0.05 are indicated in bold.

Similar correlations between circulating miRNAs and bone microstructure parameters were observed at the tibia (Table 3). miR‐31‐5p and BV/TV were significantly correlated (r = −0.40, p = .015) and miR‐378‐5p was correlated to Tb.Th (r = −0.39, p = .017). Significant correlations were also found between miR‐203b‐3p and Tb.Th (r = −0.40, p = .014) as well as for miR‐455‐3p and Tb.Th (r = −0.36, p = .031).

Table 3

Correlation Analysis (Pearson) of Serum microRNA Levels at the Study End (Week 20) With Tibia Bone Microarchitecture

Tibia	BV/TV		Tb.Th		Tb.N		Tb.Sp
microRNA	PCC	p Value	PCC	p Value	PCC	p Value	PCC	p Value
miR‐203b‐3p	–0.19	.263	–0.40	.014	0.03	.849	0.03	.858
miR‐31‐5p	–0.40	.015	–0.29	.082	–0.24	.146	0.30	.070
miR‐378a‐5p	–0.29	.083	–0.39	.017	–0.08	.626	0.14	.417
miR‐455‐3p	–0.04	.801	–0.36	.031	0.27	.115	–0.20	.254

BV/TV = trabecular bone volume fraction; Tb.Th = trabecular thickness; Tb.N = trabecular number; Tb.Sp = trabecular spacing; PCC = Pearson correlation coefficient.

N = 37; p < .05 are indicated in bold.

Associations between miRNA bone tissue expression and serum levels

miRNA serum levels were correlated with miRNA expression in femoral head bone tissue, which had previously been assessed by small RNA sequencing.( ) Serum and bone levels of miR‐139‐5p were significantly correlated (r = 0.62, p = .013), whereas miR‐582‐5p (r = 0.49, p = .077), miR‐203a‐3p (r = 0.45, p = .096), and miR‐375‐3p (r = 0.46, p = .084) showed trending correlations between levels in the circulation and bone tissue (Table 4).

Table 4

Correlation Analysis (Pearson) of Serum microRNA Levels at the Study End (Week 20) With Bone microRNA Expression

microRNA	Bone microRNA expression
	PCC	p Value
miR‐139‐5p	0.62	.013
miR‐203a‐3p	0.45	.096
miR‐375	0.46	.084
miR‐582‐5p	0.49	.077

PCC = Pearson correlation coefficient.

N = 15; p < .05 are indicated in bold.

Longitudinal changes in serum microRNA levels over 12 weeks of antiresorptive or osteoanabolic treatment

On the basis of differential expression analysis as well as correlation analysis of serum miRNA levels with different bone parameters, a panel consisting of 29 of the 80 analyzed miRNAs was selected for longitudinal analysis in 128 serum samples collected at week 8, week 12, and week 16 (Table 1). After quality control, 118 serum samples and all 29 miRNAs were included in the analysis and two‐way ANOVA was used to test for significant effects of time point, treatment, or time point x treatment (Table 5).

Indeed, serum levels of 12 miRNAs were significantly influenced by the time point (10 miRNAs) or the treatment (1 miRNA) or showed a significant interaction between the two factors time point and treatment (3 miRNAs) (Table 5). The 12 miRNAs identified were subjected to statistical comparison of all four individual time points within each treatment group. Eight miRNAs showed significant differences between individual time points that were dependent on the respective treatment group (Figs. 5 and 6).

Fig 5

MicroRNAs regulated by ovariectomy‐induced osteoporosis and rescued by osteoanabolic and antiresorptive treatment. Four miRNAs (miR‐30d‐3p, miR‐31‐5p, miR‐34a‐5p, and miR‐378a‐5p) significantly increased in vehicle‐treated OVX rats throughout the study from week 8 to week 20. Treatment of OVX rats with teriparatide and zoledronate, respectively, was able to prevent this upregulation.

Fig 6

MicroRNAs regulated by ovariectomy‐induced osteoporosis and rescued by osteoanabolic but not antiresorptive treatment. (A) Two miRNAs (miR‐375‐3p and miR‐183‐5p) behaved similarly in SHAM and TPD and VEH and ZOL throughout the study; hence, either an upregulation was observable in VEH and ZOL that was prevented in TPD (miR‐375‐3p) or a downregulation was observed in SHAM and TPD (miR‐183‐5p). (B) Two miRNAs (miR‐203a‐3p and miR‐203b‐3p) showed opposing trends in VEH (upregulated) and TPD (downregulated). For miR‐203a‐3p, a similar trending downregulation was observed in ZOL from week 8 to week 16 that was followed by an upregulation from week 16 to week 20.

The levels of miR‐30d‐3p (log2FC = 1.62, p < .05), miR‐34a‐5p (log2FC = 2.45, p < .05), miR‐31‐5p (log2FC = 1.21, p < .05), and miR‐378a‐5p (log2FC = 2.25, p < .01) significantly increased in VEH‐treated OVX rats from week 8 to week 20 (Fig. 5). Within the same time frame, serum levels for these four miRNAs remained stable under TPD and ZOL therapy; hence, osteoanabolic as well as antiresorptive therapy were able to prevent this effect.

Longitudinal changes in serum levels between week 12 and week 20 were identified for miR‐375‐3p and miR‐183‐5p only in the VEH and ZOL groups, but not SHAM or TPD‐treated animals (Fig. 6). miR‐375‐3p levels increased between week 12 and week 20 in VEH (log2FC = 0.70, p < 0.05) and ZOL (log2FC = 0.66, p < .05), but remained unchanged in TPD and SHAM. Conversely miR‐183‐5p serum levels decreased between week 12 and week 20 in TPD (log2FC = −1.29, p < .01) and SHAM (log2FC = −1.17, p < .05), whereas levels in VEH and ZOL followed a different trajectory during the whole study and had unchanged levels comparing week 8 and week 20.

miR‐203a‐3p and miR‐203b‐3p are known inhibitors of osteoblast differentiation and their transcription is upregulated in bone tissue of OVX animals.( ) Serum levels of both miRNAs are correlated; however, miR‐203a‐3p is higher abundant compared with miR‐203b‐3p. By taking advantage of the multiple time points in the longitudinal analysis, we observed a trend toward upregulation in serum of VEH animals for both miRNAs (log2FC = 0.45, p = n.s., log2FC = 1.76, p < .1, respectively) and downregulation in response to TPD (log2FC = −0.82, p = n.s., log2FC = −0.36, p = n.s., respectively) (Fig. 6). Furthermore, we observed that a single injection with zoledronate caused a comparable downregulation of both miRNAs until week 16; however, between weeks 16 and 20, a significant rebounding increase was observed (log2FC = 1.45, p < .05).

Gene target and network analysis reveal WNT and NOTCH signaling as potential targets for identified miRNA candidates

Twelve miRNAs selected by two‐way ANOVA (Table 5) were mapped to their experimentally verified target genes using mirnet 2.0. The resulting list of gene~miRNA interactions was filtered for genes with at least four miRNA interactions, narrowing down this list to seven genes (Fig. 7). Gene set enrichment analysis based on the reactome classification identified several potentially relevant pathways, including NOTCH and WNT signaling (Fig. 7).

Fig 7

Target analysis of differentially expressed microRNAs (Table 5) using miRnet. (A) Twelve microRNAs identified as significantly influenced by time point and/or treatment in two‐way ANOVA analysis were used for constructing a target network using the online tool miRnet (accessible via mirnet.ca). The degree filter for constructing the network was set to 4; hence, only target nodes with at least four connections remained in the network. (B) Pathway enrichment using the Reactome database identified 18 pathways significantly enriched with an adjusted p value <.05. NOTCH and Wnt signaling were found as the top enriched pathways via the targets MYC and AGO3 (highlighted in yellow).

Discussion

In the present study, we have used an established animal model for postmenopausal osteoporosis to demonstrate that bone loss as well as two common forms of anti‐osteoporotic therapies impact serum miRNA levels throughout disease progression and treatment response.

Based on an unbiased small RNA sequencing analysis of femoral head samples and subsequent cross‐sectional analysis of matched serum samples, we selected a panel of 29 miRNAs for analysis at four different time points. We found that serum levels of 12 of 29 miRNAs changed significantly over time (10 miRNAs) or depending on the type of treatment (1 miRNA) or showed a significant interaction between the two factors time point and treatment (3 miRNAs). These 12 miRNAs were used for statistical post hoc analysis of all individual time points, which revealed eight miRNAs with distinct changes over time as a consequence of OVX‐induced osteoporosis. Importantly, these miRNAs showed differential responses to either osteoanabolic treatment with TPD or antiresorptive treatment using ZOL.

Progression of OVX‐induced osteoporosis led to an increase of miR‐34a‐5p, miR‐378a‐5p, miR‐30d‐3p, and miR‐31‐5p, whereas levels for all four miRNAs remained stable throughout the study under TPD and ZOL treatment.

Interestingly, the increase of miR‐34a‐5p in the circulation of vehicle‐treated OVX rats was also observed in femoral head tissue and reversed by TPD treatment.( ) miR‐34a was previously reported to interact with several targets within the WNT signaling pathway( ) and shown to be capable of antagonizing WNT signaling by targeting WNT1 mRNA.( , , ) This connection is further underscored by a recent study that characterized circulating miRNA profiles in individuals with decreased WNT signaling due to a WNT1 mutation. Serum levels of both miR‐34a‐5p and miR‐31‐5p were found to be decreased and suggested that decreased WNT signaling in these patients might result in a compensatory reduction of miRNAs that otherwise would repress mRNAs of the WNT pathway and thus further decrease WNT signaling.( )

For miR‐31‐5p, we observed the most significant cross‐sectional effect in serum at week 20 (Fig. 4) that was further confirmed in the longitudinal analysis by a significant increase in the circulation of OVX rats between week 8 and week 20. During the same time frame, miR‐31‐5p levels remained constant in TPD and ZOL (Fig. 5). Several former studies have reported an important role of miR‐31‐5p in different aspects of bone biology: overexpression of miR‐31 suppressed osteogenesis of hMSCs, while its downregulation improved it.( ) Senescent endothelial cells upregulate miR‐31 transcription and can transfer miR‐31 via extracellular vesicles to hMSCs, where osteoblastogenesis is impaired through downregulation of frizzled 3 (FZD3), inhibiting WNT signaling. This effect could be rescued by antagonizing miR‐31‐5p.( ) In line, inhibition of miR‐31‐5p using antagomiRs resulted in an increase of osterix (Osx) in hMSCs and pre‐osteoblastic cells.( ) We have recently observed significant correlations of miR‐31‐5p and the corresponding mRNAs Dlx5 (distal‐less homeobox 5), Runx2 (Runt‐related transcription factor 2), and Osx (our unpublished data). All of them are key regulators of bone formation,( ) highlighting the importance of miR‐31‐5p in bone formation and the pathophysiology of osteoporosis. Thus, the measurement of miR‐31‐5p in serum might indeed have an impact for clinical diagnosis, prognosis, and follow‐up of treatment. Because inhibition of miR‐31‐5p also dampens osteoclastogenesis,( ) it might represent a therapeutic target in osteoporotic treatment targeting both the anabolic and catabolic arms of bone metabolism.

miR‐378a‐5p and miR‐30d‐3p serum levels increased over time as a result of ovariectomy‐induced osteoporosis. This effect was rescued by treatment with TPD and ZOL. We also observed a moderate increase of circulating miR‐30a‐5p in postmenopausal women with recently sustained osteoporotic fractures.( ) Additionally, miR‐30d‐3p was found to be mechanosensitive and upregulated in an in vitro model of unloading conditions and with inhibitory effects on osteogenic differentiation via targeting RUNX2,( ) as well as to be increased in postmenopausal osteoporotic women.( ) Reduced plasma levels were found after exercise in postmenopausal women,( ) which is well known to positively impact on bone metabolism.( ) Taken together, it seems that higher levels of miR‐30d‐3p are associated with bone loss and lower levels with increased bone formation.

Similarly, we report that miR‐203a‐3p, miR‐203b‐3p, and miR‐375‐3p levels increase in serum with progression of osteoporosis but are rescued by TPD. ZOL treatment showed similar effects until week 16, but then miR‐203b‐3p levels increased at week 20. Because ZOL was only administered once (at week 8), this could lead to a reduced efficacy of zoledronic acid after 8 to 12 weeks of treatment compared with the daily administered teriparatide treatment. Clinical evidence suggests that therapeutic effects of ZOL may last for several years after administration;( ) hence, circulating miRNAs could be of clinical value for monitoring treatment efficacy over time and identifying intervals in between treatments. The steady increase of miR‐203 serum levels in OVX animals is in line with our previous finding that the transcription of both miR‐203a and miR‐203b in bone tissue is highly upregulated in osteoporotic animals compared with control. The fact that TPD treatment reduces miR‐203a to the baseline level observed in non‐osteoporotic control animals provides evidence for the efficacy of TPD therapy. In human patients, miR‐203a is upregulated in serum of diabetic and non‐diabetic postmenopausal women with osteoporosis and a history of low‐traumatic fractures( ) and fully in line with our data.

Functionally, miR‐203a‐3p is known to be involved at different steps in bone formation. For example, the differentiation of osteoblasts is slowed by miR‐203a through direct silencing of Dlx5 and RUNX2,( , ) as well as of Smad genes in the context of diabetes‐associated osteogenesis genes.( )

Yet another miRNA responding to anti‐osteoporotic treatment is miR‐375, and again several links to bone biology are documented. Postmenopausal women with incident vertebral fractures were shown to have higher levels of miR‐375 in serum in comparison to non‐fractured controls, ) and miR‐375 was shown to be capable of inhibiting osteogenesis via targeting LRP5 and β‐catenin.( )

The level of miR‐183‐5p does not change with induction of osteoporosis, whereas it decreases with time in non‐ovariectomized rats. Interestingly, TPD treatment results in a similar decrease, whereas ZOL does not. In vitro studies revealed a role for miR‐183‐5p in osteoclastogenesis and to be upregulated by stimulation with RANKL, whereas inhibition of this miRNA resulted in reduced osteoclastogenesis.( ) In addition, this miRNA was also found to be increased in extracellular vesicles derived from bone marrow interstitial fluid from aged mice.( )

To assess whether these 12 miRNAs might cooperate to impact gene expression, we looked for common experimentally confirmed target mRNAs. Surprisingly, 8 of 12 miRNAs were found to target either AGO3 or MYC (also known as c‐MYC). Pathway enrichment using the reactome database revealed significant associations of both genes with NOTCH and Wnt signaling, both with key roles in bone remodeling.

Interestingly, it was previously shown that myeloid‐specific MYC depletion increases bone mass and protects from OVX‐induced osteoporosis by impairing osteoclastogenesis.( ) On the other hand, a recent study showed that miR‐203a and miR‐203b are transcriptionally responsive to c‐MYC.( ) Further, miR‐34a‐5p, miR‐375, and miR‐30d are known to be capable of antagonizing c‐MYC (Fig. 7).( , , ) Thus, the expression of MYC possibly represents an important connecting node for the observed changes in serum and prompts further investigation.

So far, only a few studies have been published dealing with anti‐osteoporotic treatment and miRNAs. Zarecki and colleagues observed that selected bone‐related miRNAs were upregulated in osteoporosis patients with vertebral fractures compared with non‐fractured patients and healthy controls. They also found a non‐significant treatment effect with lower levels in those treated with oral bisphosphonates. miR‐375 was lower in oral bisphosphonate pretreated subjects compared with untreated subjects.( ) In the present study, miR‐375‐3p was rescued by TPD but not ZOL.

Moreover, TPD therapy led to a different expression in postmenopausal osteoporosis patients. miR‐133‐3p and miR‐33‐3p reached statistical significance in a previously published study but were not examined in our study.( )

Bone microstructure is one of the main components of bone strength. The association between deteriorations of bone microarchitecture and fracture risk has been reported repeatedly.( , ) Teriparatide, an osteoanabolic agent, increases trabecular and cortical bone microstructure indices.( , ) Using high‐resolution micro‐CTs, we found a significant increase in BMD and BV/TV during TPD treatment. At the end of the study, Ct.Th (data not shown) as well as Tb.Th were highest in TPD‐treated animals. miR‐31‐5p, miR‐203b‐3p, and miR‐378a‐5p, which were the most significantly downregulated miRNAs under TPD treatment, showed a significant negative correlation to trabecular bone microstructure. These data indicate that circulating miRNA levels might reflect the improvement of bone microstructure during osteoanabolic treatment. Interestingly, also miR‐203a‐3p was highly correlated to indices of bone microstructure. This miRNA was found previously to be the most significantly upregulated miRNA in bone tissue in ovariectomized untreated rats.( ) TPD treatment reversed this effect in bone and in serum. The correlation between levels in bone tissue and serum levels suggest miRNAs to be valuably predictive of bone‐specific processes in blood.

As expected, weaker effects on bone microstructure were found in the zoledronate treatment group. Still, it partially rescued bone loss. In the cross‐sectionally upregulated miRNAs in ZOL compared with untreated control at week 20, miR‐582‐5p, miR‐204‐5p, and miR‐214‐3p (Fig. 4) were not correlated to bone microstructure, suggesting no beneficial role for these miRNAs in the prediction of bone microstructure.

Presently, the follow‐up of osteoporosis and the decision for continuation as well as discontinuation of therapy is based on clinical risk factors, incident fractures as well as changes in BMD and/or BTM.( ) However, BMD changes during antiresorptive therapy only partly explain the risk reduction for vertebral and non‐vertebral fractures.( ) Moreover, the increase of areal BMD during TPD treatment does not sufficiently reflect the magnitude of gain in bone organic matrix, mineral content, or bone microstructure assessed by high‐resolution techniques.( ) Additionally, the least significant change of BMD changes has to be considered.( ) Therefore, new and noninvasive tools for treatment monitoring are required. Changes in circulating miRNAs during anti‐osteoporotic treatment, as observed in this study, could reflect treatment response and non‐response on specific therapy in patients with osteoporosis in future.

In summary, the effects observed in the present study of OVX‐induced osteoporosis and anti‐osteoporotic treatment on bone loss and treatment response that were followed by quantitative changes of bone‐related miRNAs in serum suggest potential clinical value for treatment monitoring, and clinical investigations are needed to confirm our preclinical results.

Disclosures

RK, JF, GL, PH, XF, and JZ state that they have no conflicts of interest regarding this project. MH and JG are co‐founders of TAmiRNA. MH and MW are employed by TAmiRNA GmbH. JG, PP, and HR are scientific advisors to TAmiRNA. MH and JG hold patents related to the use of circulating microRNAs as biomarkers for bone diseases.

AUTHOR CONTRIBUTIONS

Moritz Weigl: Conceptualization; data curation; formal analysis; investigation; methodology; visualization; writing‐original draft; writing‐review & editing. Roland Kocijan: Conceptualization; data curation; formal analysis; methodology; writing‐original draft; writing‐review & editing. James Ferguson: Investigation; methodology; writing‐original draft. Gabriele Leinfellner: Investigation; methodology; writing‐original draft. Patrick Heimel: Data curation; methodology; writing‐original draft. Xaver Feichtinger: Writing‐original draft. Peter Pietschmann: Conceptualization; writing‐original draft. Johannes Grillari: Conceptualization; formal analysis; writing‐original draft. Jochen Zwerina: Conceptualization; writing‐original draft. Heinz Redl: Conceptualization; writing‐original draft.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1002/jbmr.4276.

Supplemental Fig. S1. Quality control of circulating microRNA analysis in rodent serum samples. Raw Cq‐values obtained for RNA spike‐in control (UniSp4) and cDNA spike‐in control (cel‐miR‐39) are plotted for all 118 serum samples included in the data analysis.

Click here for additional data file.