Global miRNA expression and correlation with mRNA levels in primary human bone cells.

MicroRNAs (miRNAs) are important post-transcriptional regulators that have recently introduced an additional level of intricacy to our understanding of gene regulation. The aim of this study was to investigate miRNA-mRNA interactions that may be relevant for bone metabolism by assessing correlations and interindividual variability in miRNA levels as well as global correlations between miRNA and mRNA levels in a large cohort of primary human osteoblasts (HOBs) obtained during orthopedic surgery in otherwise healthy individuals. We identified differential expression (DE) of 24 miRNAs, and found 9 miRNAs exhibiting DE between males and females. We identified hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b and their target genes as important modulators of bone metabolism. Further, we used an integrated analysis of global miRNA-mRNA correlations, mRNA-expression profiling, DE, bioinformatics analysis, and functional studies to identify novel target genes for miRNAs with the potential to regulate osteoblast differentiation and extracellular matrix production. Functional studies by overexpression and knockdown of miRNAs showed that, the differentially expressed miRNAs hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b target genes highly relevant to bone metabolism, e.g., collagen, type I, α1 (COL1A1), osteonectin (SPARC), Runt-related transcription factor 2 (RUNX2), osteocalcin (BGLAP), and frizzled-related protein (FRZB). These miRNAs orchestrate the activities of key regulators of osteoblast differentiation and extracellular matrix proteins by their convergent action on target genes and pathways to control the skeletal gene expression.

INTRODUCTION

MicroRNAs (miRNAs) are small (∼22 nt long), single-stranded noncoding RNAs which regulate gene expression primarily by binding to partially complementary sites in the 3′ UTR of mRNAs and inhibit gene expression either by facilitating mRNA degradation or by inhibition of translation (Bartel 2004; Doench and Sharp 2004; Engels and Hutvagner 2006; Rajewsky 2006; Hu and Coller 2012; Wilczynska and Bushell 2015). Recent evidence from Bartel et al., however, suggests that the predominant mechanism by which miRNAs regulate protein levels of target genes is mRNA stability rather than on translation (Guo et al. 2010). Recent studies indicate that over one-third of all human genes are regulated by microRNAs (Davis and Hata 2009). Studies suggest that a single gene can be targeted by a cluster of miRNAs, but also that a single miRNA can target many protein coding genes. As the field of miRNA studies has evolved, the actions of miRNAs have been shown to be surprisingly complex, and intricate networks of factors stringently control miRNA processing and biological activities.

Osteoporosis is a common disease, characterized by diminished bone mineral density (BMD) and bone strength and consequently an increased risk for fractures. The genetic contribution to osteoporosis is well known and much effort has been made to identify genetic variants associated with the susceptibility to osteoporotic fractures; both for the discovery of novel drug targets but also for the identification of risk factors.

Expression of different subsets of bone-regulating miRNAs (with the proposed name osteomiRs) during different stages of osteoblast maturation and osteoclastogenesis have been reported, regulating the timing and progression of cellular differentiation programs (Lian et al. 2012).The miRNA families; hsa-miR-29, hsa-miR-let-7, and hsa-miR-26 influence bone formation by targeting a number of collagens and other extracellular matrix proteins (Zhao et al. 2014). It has been shown that hsa-miR-125b down-regulates cell proliferation by inhibiting osteoblast differentiation (Miyazono et al. 2005). Hsa-miR-29b is a positive regulator of osteoblast differentiation, targeting bone matrix RNAs, including COL1A1 (Kapinas et al. 2010) and COL3A1, down-regulating inhibitors such as histone deacetylase (HDAC) 4, TGFβ3, AcvR2A, and CTNNBIP1 (catenin β-interacting protein I) in rat and mouse cells (Li et al. 2009). Studies on human mesenchymal stem cells (MSC) revealed that hsa-miR-31 negatively regulates Osterix expression, which in turn targets BGLAP and COL1A1 genes during early osteoblastic differentiation, thereby making hsa-miR-31 a potentially important regulator of bone mineralization (Baglìo et al. 2013). Recent studies identified nine up-regulated miRNAs in osteoporotic patients, suggesting that alterations in the levels of circulating miRNAs in serum is associated with either increased osteoclastogenesis or inhibited osteoblast differentiation, and that circulating miRNA levels could be used as novel biomarkers for diagnostic purposes (Seeliger et al. 2014).

RESULTS

Global mRNA expression

The global mRNA-expression profiles from primary human osteoblasts (HOBs) from 95 individuals were obtained using the Illumina HumRef-8v2 arrays. Microarray data have been deposited in the Gene Expression Omnibus (GEO) (www.ncbi.nlm.nih.gov/geo, accession no. GSE15678).

microRNA profiling in human bone samples by microarray

Normalization and data filtering of the 757 probes resulted in 251 high quality probe signals passing predetermined quality criteria, and these were used in downstream analyses (Supplemental Table 1a). In the initial analysis of data, we determined relative levels of miRNA expression and differentially expressed (DE) miRNAs, defined as having M-values with SD > 0.5, resulting in 24 DE miRNAs overall (Fig. 1A; Supplemental Table 1a). When assessing gender differences, 9 of 251 miRNA had significantly different expression levels between males and females in two-tailed t-test with a P-value cut-off at <0.05 (Fig. 1B; Supplemental Table 1b). The gender difference is further exemplified in Supplemental Figure 2 using principal component analysis (PCA) showing separation of female and male clusters. The annotation miRPlus in tables and figures refers to Exiqon in-licensed human sequences not yet annotated in miRBase (latest release version 21: June 2014).

FIGURE 1.

Heatmaps representing miRNA expression. (A) Heatmap of miRNA expression by LNA microarray in human osteoblasts from eight males and 11 females. Out of the 251 expressed miRNAs, a heatmap was generated for the 24 miRNAs differentially expressed by log2 (Hy3/Hy5) ratios with a variation across all samples having SD > 0.50. For comparison, expression levels for hsa-miR-125b are shown below the 24 DE miRNAs. (B) Heatmap depicting variation in miRNA expression between females and males. Out of the expressed miRNA, nine exhibited differential expression levels between males and females based on log2 (Hy3/Hy5) ratios significantly different between males and females using a two-tailed t-test and P-values <0.05. Clustering was made and heatmaps were processed and displayed using MultiExperiment Viewer.

microRNA profiling in human bone samples by real-time RT-PCR

From the 251 expressed miRNAs identified from the LNA microarray experiment a subset of miRNAs with readily available assays for verification was chosen. The miRNAs were selected based on the LNA array data with either DE between genders with a P < 0.05 (hsa-miR-29b, hsa-miR-99a, and hsa-miR-140-3p), and/or overall DE with an SD > 0.5 (hsa-miR-30c2, hsa-miR-503, hsa-miR-31, hsa-miR-335, hsa-miR-22, and hsa-miR-198), or based on high Hy3/Hy5 average intensity on the LNA array (hsa-miR-125b, data not shown), and for use as controls (hsa-miR-191 and hsa-let-7a, see below). The expression profiles of these 12 miRNAs were successfully quantitated in HOBs from 95 individuals using TaqMan MicroRNA Assays. There was no correlation between miRNA levels and age of the bone cell donor (P = 0.23–0.97, bivariate analysis by age). The miRNAs exhibited some degree of correlation overall, less pronounced for hsa-miR-335, hsa-miR-198, and hsa-miR-29b. The correlation pattern differed slightly between males and females (Supplemental Table 2). The validation by TaqMan analyses showed that 9/10 miRNAs (hsa-miR-125b, hsa-miR-29b, hsa-miR-30c2, hsa-miR-31, hsa-miR-99a, hsa-miR-140-3p, hsa-miR-198, hsa-miR-335, and hsa-miR-503) exhibited DE with overall SD 0.9–17.2. Four out of 10 tested miRNAs showed the same trend for differences in expression between males and females as was previously observed in the microarray analysis; hsa-miR-125b (P = 0.089), hsa-miR-29b (P = 0.113), hsa-miR-198 (P = 0.064), with hsa-miR-31 (P = 0.049) reaching statistical significance (Fig. 2B).

FIGURE 2.

miRNA-expression profiles showing segregation between females and males. (A) Venn diagram depicting the number of correlations miRNA versus mRNA with P < 0.05, and the number of shared correlations between males, females, and the overall data set. (B) TaqMan results for the subset selected for verification of miRNA LNA array data. The selection was based on overall differential expression and/or differential expression between males and females. Bars represent means with SD. P-values for T-tests for independent samples of males versus females are indicated above the respective bars.

Correlations between miRNA- and mRNA levels

Analysis revealed large numbers of significant correlations at the P < 0.05 level with preponderance for significant correlations in females (Fig. 2A). Using an FDR cut-off of 0.05, the number of remaining correlations was reduced (Table 1), but the pattern of more correlations being observed for females remained. The number of positive correlations was generally higher than negative correlations, with a mean ratio for positive versus negative correlations of 1.5 (range 0.9–2.7). hsa-miR-29b showed the highest number of significant correlations followed by an intermediate group consisting of hsa-miR-335, hsa-miR-140-3p, hsa-miR-99a, hsa-miR-31, hsa-miR-22, and hsa-miR-125.

TABLE 1.

miRNA:mRNA correlations

Quantile–quantile (Q–Q) plots performed on the 10 tested miRNAs in the overall data set showed substantial egression of the observed findings from the null hypothesis of no correlation, indicative of significant correlations (Supplemental Fig. 3). In addition, an overrepresentation analysis showed that mRNAs having a predicted hsa-miR-29b target site in their 3′ UTR were fourfold enriched among those mRNA most significantly correlated with miR-29b (FDR < 0.01) (Supplemental Fig. 4). This further strengthens the case that the significant correlations we observed for hsa-miR-29b are indeed true primary effects due to miRNA–mRNA interactions.

A number of genes involved in bone metabolism exhibited a statistically significantly correlation between miRNA and mRNA expression. COL1A1 (P = 0.0001) and SPARC (P = 0.001) showed a negative correlation to hsa-miR-29b (Fig. 3B), RUNX2 (P = 0.0148), FRZB (P = 0.003), and BGLAP (P = 0.0325) were negatively correlated with hsa-miR-30c2 (Fig. 4B) and FRZB (P = 0.0117) and BGLAP (P = 0.0450) were negatively correlated with hsa-miR-125b (Fig. 5B).

FIGURE 3.

hsa-miR-29b negatively regulates Col1A1 and SPARC in primary human bone cells. (A) Predicted binding sites for miR-29b in Col1A1 and SPARC 3′ UTR by the TargetScan, miRanda, and PicTar softwares. (B) Pearson correlation scatter plots showing a negative correlation between hsa-miR-29b and COL1A1 (P = 0.001) and with SPARC (P = 0.001). (C) RT-PCR analysis of COL1A1 and SPARC mRNA levels normalized to GAPDH. (D) The amount of COL1A1 protein was assessed by Western blot analysis with β-actin as a loading control (left panel) and SPARC protein secreted into conditioned media was quantified by Western blot. Mean ± SD * is significantly different from negative control (NC) (right panel). (E) ALP staining was visualized at 10× magnification. (F) Deposition of calcium was visualized (left panel) and quantified (right panel), by Alizarin red stain in HOBs after they were treated in osteogenic medium but for 12 d.

FIGURE 4.

hsa-miR-30c2 negatively regulates RUNX2, FRZB, and BGLAP in primary human bone cells. (A) Predicted binding sites for hsa-miR-30c2 in RUNX2, FRZB, and BGLAP 3′ UTR by the TargetScan, miRanda, and PicTar softwares. (B) Pearson correlation scatter plots showing a negative correlation between hsa-miR-30c2 RUNX2 (P = 0.0148), FRZB (P = 0.003), and BGLAP (P = 0.0325). (C) RT-PCR analysis of RUNX2, FRZB, and BGLAP mRNA levels normalized to GAPDH. (D) The amount of RUNX2 and FRZB was assessed by Western blot analysis with β-actin as a loading control (left panel) and BGLAP protein secreted into conditioned media was quantified by Western blot. Mean ± SD * is significantly different from negative control (NC) (right panel). (E) ALP staining was visualized at 10× magnification. (F) Deposition of calcium was visualized (left panel) and quantified (right panel), by Alizarin red stain in HOBs after they were treated in osteogenic medium but for 12 d.

FIGURE 5.

hsa-miR-125b negatively regulates BGLAP and FRZB in primary human bone cells. (A) Predicted binding sites for hsa-miR-125b in BGLAP and FRZB 3′ UTR by the TargetScan, miRanda, and PicTar softwares. (B) Pearson correlation scatter plots showing a negative correlation between hsa-miR-125b and BGLAP (P = 0.0450) and FRZB (P = 0.0117). (C) RT-PCR analysis of BGLAP and FRZB mRNA levels normalized to GAPDH. (D) The amount of FRZB protein was assessed by Western blot analysis with β-actin as a loading control (left panel) and BGLAP protein secreted into conditioned media was quantified by Western blot. Mean ± SD * is significantly different from negative control (NC) (right panel). (E) ALP staining was visualized at 10× magnification. (F) Deposition of calcium was visualized (left panel) and quantified (right panel), by Alizarin red stain in HOBs after they were treated in osteogenic medium but for 12 d.

Bioinformatic pathway analysis

The miRNA versus mRNA correlation data set was interrogated by the Ingenuity Pathway Analysis software (IPA, Ingenuity Systems). Analyses were performed for the overall data set, and for male and females assessed separately. IPA analysis on the miRNA versus mRNA correlation data set yielded several pathways with known or putative relations to bone metabolism reaching a P-value <0.05. Although the majority of the significant pathways were unique for each miRNA, one pathway reaching statistical significance was the Wnt/ β-catenin signaling pathway, for hsa-miR-29b, hsa-miR-125b, and hsa-miR-30c2 (Supplemental Fig. 5 for hsa-miR-29b).

hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b regulate osteogenesis

Functional studies of hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b were performed by overexpression and knockdown studies in HOBs. The putative binding sites in the 3′ UTRs of hsa-miR- 29b in COL1A1 and SPARC mRNA (Fig. 3A), of hsa-miR-30c2 RUNX2, FRZB, and BGLAP (Fig. 4A), of hsa-miR-125b FRZB and BGLAP (Fig. 5A), were predicted by the softwares TargetScan, miRanda, and PicTar. The overexpression studies show that augmented amount of hsa-miR-29b repressed Col1A1 and SPARC, whereas hsa-mir-29b inhibitor increased Col1A1 and SPARC at the mRNA and protein level. Augmented amounts of hsa-miR-30c2 repressed RUNX2, FRZB, and BGLAP, whereas hsa-mir-30c2 inhibitor elevated the amount of RUNX2, FRZB, and BGLAP. Augmented amounts of hsa-miR-125b repressed FRZB and BGLAP, whereas hsa-miR-125b inhibitor elevated the amount of FRZB and BGLAP as assessed by RT-PCR (Figs. 3–5C) and Western blotting (Figs. 3–5D). In vitro matrix mineralization levels was visualized and quantified by ALP and Alizarin red staining, showing increased mineralization with hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b inhibitors and a decreased mineralization by hsa-miR-29b, hsa-miR-30c2, and hsa-miR-mimics (Figs. 3–5E,F).

DISCUSSION

The human genome is estimated to encode >1000 miRNAs which may regulate up to ∼60% of human genes (Lian et al. 2012). In recent years, a number of studies have identified unique miRNA signatures for several diseases (Philippidou et al. 2010). Advances in miRNA basic research, discovery of new miRNAs and their target networks have been paralleled by developments and improvements in miRNA detection techniques. In our study, we have used sensitive LNA miRNA microarrays and specific qPCR protocols to assess the miRNAome in primary human bone cells.

To explore whether microRNAs have a role in human bone dynamics it is necessary to determine which miRNAs are active within bone. So far, the function of a number of miRNAs has been characterized in bone; e.g., hsa-miR-125b, hsa-miR-29b, and hsa-miR-30a have been associated with osteoblast differentiation in human multipotent mesenchymal stromal cells (MSCs) (Gao et al. 2011). In the present study, the miRNA-expression profile of primary human bone cells was analyzed, revealing expression of ∼250 miRNAs in HOBs. We then further validated differential expression of a subset of 10 miRNAs in 95 individuals by qPCR, and correlated the miRNA levels with global mRNA levels in these same individuals and found that hsa-miR-29b had the largest number of significant correlations to mRNA levels.

Validation of miRNA targets is an ongoing process, and the TarBase v6.0 database (Vergoulis et al. 2012) currently contains ∼65,000 experimentally validated miRNA–gene interactions. A number of these validated targets are present in our data set. The predicted targets of the miRNAs in this study were also verified using the softwares TargetScan, miRanda, and PicTar. hsa-mir-29b, hsa-mir-30c2, and hsa-mir-125b were differentially expressed in our data set and for these miRNAs, our miRNA–mRNA correlation studies show a number of statistically significant correlations for genes involved in bone metabolism, giving these three miRNAs high priority for assessing their effects on osteoblast differentiation and mineralization.

hsa-miR-29b has been shown to target the most abundant collagenous matrix gene COL1A1 (Li et al. 2009), and hsa-miR-29a and -29c to target the most abundant noncollagenous matrix gene SPARC (Kapinas et al. 2009). Our data now show a novel negative correlation between hsa-miR-29b and SPARC mRNA levels, a finding which we then further confirmed by functional studies at the protein level. Synthesis of collagenous and noncollagenous proteins is crucial for the formation of extracellular matrix by the osteoblast. hsa-miR-29b suppresses the expression of collagen proteins allowing the collagen fibril matrix to mature for mineral deposition (Fig. 6). Previous studies have shown that collagen proteins are not suppressed in immature osteoblasts but attenuated only in mature osteoblasts (Li et al. 2009). As osteoblasts differentiate, the expression and synthesis of extracellular matrix genes are reduced (Kalajzic et al. 2005). In our functional studies on HOBs and matrix proteins, we confirm that COL1A1 shows a significant negative correlation and also show a novel negative correlation between and SPARC and hsa-miR-29b both on the mRNA and protein level.

FIGURE 6.

Regulation of osteoblast differentiation by miRNAs. MiRNAs influence different stages of osteoblast differentiation from mesenchymal stem cell to osteocyte. MiRNAs target important transcription factors and genes and affect regulation during each step of differentiation. MiRNAs may have a positive or negative effect depending on the target and stage of osteoblast differentiation. Runt-related transcription factor 2 (RUNX2), osteocalcin (BGLAP), and frizzled-related protein (FRZB).

RUNX2 is a key regulator of osteogenesis and it is subjected to multiple levels of regulation of miRNAs. During early stages of osteoblast differentiation, RUNX2 is up-regulated, activating major matrix genes, while RUNX2 is down-regulated in mature osteoblasts (Fig. 6). It is well established that hsa-miR-30c mediates inhibition of osteoblastic differentiation by targeting RUNX2 (Zhang et al. 2011). In our functional studies on HOBs, we confirm a significant negative effect of hsa-miR-30c2 on RUNX2. In addition, our results also show a novel effect of hsa-miR-30c2 and hsa-miR-125b with negative correlations to expression levels of FRZB and the late-stage osteoblast marker BGLAP. Our pathway analysis showed many statistically significant pathways that were unique for each miRNA, while one pathway that reached statistical significance shared by many miRNA was the Wnt/β-catenin signaling pathway. MiRNAs modulate Wnt signaling by targeting and suppressing negative regulators of the Wnt pathway, e.g., FRZB/sFRPs. Our functional results show that hsa-miR-30c2 and hsa-miR-125b are negatively correlated to FRZB/sFRPs indicating an activating effect on the Wnt/β-catenin pathway. Furthermore, our results support previous findings that hsa-miR-29b is a key down-regulator of Wnt antagonists like Dkk1, Kremen2, sFRP2, and WIF1 (Zhou et al. 2011). The data also show a positive correlation of hsa-miR-29b with Wnt1, one of the activators of canonical Wnt signaling pathway and β-catenin which activates target gene transcription (Li et al. 2011).

Compared with the paradigm that miRNAs are repressive in nature, our miRNA versus mRNA correlation data showed a higher number of significant positive correlations than negative correlations, which is intriguing. Though unexpected, the detection of positive correlations between miRNA and miRNA has been reported in many studies (Vasudevan et al. 2007; Nunez-Iglesias et al. 2010). Gene regulatory networks often include feedback motifs, and miRNAs are interwoven into such complex regulatory networks which coordinate transcriptional regulation in signaling networks (Inui et al. 2010). The overrepresentation of positive correlations in our data could be due to downstream effects, perhaps due to inhibited repression or perturbations of feedback loops (Martinez et al. 2008; Sun et al. 2012).

The aim of this study was to identify miRNAs involved in bone metabolism, using primary human osteoblasts (HOBs). By using human primary cells rather than rat or mouse, or rodent cell lines, we hope to circumvent issues regarding species or strain-specific effects. Our study is unique in demonstrating that expression of miRNAs exhibits interindividual variation in primary bone cells from a large cohort and also in showing differential expression of miRNAs between males and females, which has not been shown before in studies of bone cells. Although the initial findings regarding gender differences found in 20 subjects by LNA array showed the same trend when analyzing by TaqMan in HOBs from 95 subjects, only miR-31 remained significant. The difference in miRNA levels between cells from males and females was observed in culture conditions that were identical for all cells, which indicates either an underlying gender difference, or that cells from males and females have been primed epigenetically or by other means to express miRNAs differently in response to changes in hormonal milieu. These findings are intriguing, and although our data were obtained by two different quantification techniques, findings need to be replicated in larger experimental settings for putative gender differences to be further elucidated.

In conclusion, this study shows that several miRNAs exhibiting interindividual variation act in important metabolic pathways, with hsa-miR-29b, hsa-miR-30c2, and hsa-miR-125b identified as affecting osteoblast differentiation and mineralization by targeting genes highly relevant to bone metabolism. These miRNAs play a regulatory role affecting osteoblast differentiation and extracellular matrix production. In this study, we used an integrated analysis of global mRNA—miRNA correlations, bioinformatics analysis, and functional studies to identify novel target genes for miRNAs with the potential to regulate osteoblast differentiation.

MATERIALS AND METHODS

Bone cell culture and transfection

Primary human osteoblast (HOB) cells were isolated from human trabecular bone excised from the unaffected distal region of bone specimens obtained in conjunction with prosthesis surgery. Specimens were collected from 104 donors from a cohort of patients with osteoarthritis undergoing total hip or knee replacement, with no other bone-related pathologies reported (Grundberg et al. 2009, 2011). The bone chips obtained from each donor were thoroughly minced and washed with PBS. The minced bone chips were cultured in medium containing α-MEM (Sigma-Aldrich) supplemented with 2 mmol/L L-glutamine, 100 units/µL penicillin, 100 mg/µL streptomycin, and 10% fetal bovine serum (Sigma-Aldrich) at 37°C with 5% CO2 until confluence was reached. The culture medium was changed twice weekly. The study was approved by the local ethics committee (Ethical approval # Ups 03-561).

HOBs from seven individuals were seeded at a density of 35,000 cells from the second passage in a 24-well cell culture plate to achieve 60%–80% confluence at the day of transfection. To over express and inhibit the function of hsa-miR-29b, hsa-miR-125b, and hsa-miR-30c2 the cells were transfected with mirVana hsa-miR-29b mimic (Ambion, catalog no. 4464066), mirVana hsa-miR-29b inhibitor (Ambion, catalog no. 4464084), mirVana hsa-miR-125b mimic (Ambion, catalog no. 4464066), mirVana hsa-miR-125b inhibitor (Ambion, catalog no. 4464084), mirVana hsa-miR-30c mimic (Ambion, catalog no. 4464066), mirVana hsa-miR-30c inhibitor (Ambion, catalog no. 4464084), and with mirVana miRNA inhibitor negative control #1 (Ambion, catalog no. 4464077) and mirVana miRNA mimic negative control #1 (Ambion, catalog no. 4464061) at 40 nM concentrations with Magnet Assisted Transfection (MATra-si) reagent (IBA GmbH) according to the manufacturers’ protocols. Each transfection was performed in triplicate. Post-transfection, cells were incubated for 48 h.

RNA isolation

Total RNA was isolated from the HOB cells by adding 500 µL of lysis buffer (Qiagen). The cell lysates were homogenized with QIAshredder (Qiagen, Germany). RNA was extracted from the cell lysates with the RNeasy Mini Kit (Qiagen). High RNA quality was confirmed for all samples with an Agilent 2100 BioAnalyzer (Agilent Technologies, Palo Alto, CA), where samples with RNA Integrity Number (RIN) above seven were taken forward for analysis (Schroeder et al. 2006). The concentrations were determined with NanoDrop ND-1000 (NanoDrop Technologies). For RNA in the LNA microarray pipeline, RIN values were 6.2–9.5, and OD 260/280 between 1.95 and 2.03. Of the original 104 samples, a total of 95 samples from 50 males and 45 females (Supplemental Fig. 1) had RIN values above threshold and were carried forward for further analyses.

microRNA profiling in human bone samples by LNA microarray

Total RNA from 11 females and 8 males were labeled with Hy3 and common reference with Hy5 fluorescent label, using the miRCURY LNA Array power labeling kit (Exiqon). The miRCURY LNA array microarray slides were scanned using the Agilent G2565BA Microarray Scanner System (Agilent Technologies, Inc.) and the image analysis was carried out using the ImaGene 8.0 software (BioDiscovery, Inc.). The quantified signals were background corrected using Normexp with offset value 10, (Ritchie et al. 2007) and normalized using the global LOESS (LOcally WEighted Scatterplot Smoothing) nonlinear regression algorithm (details provided in Supplemental Methods), whereafter data were subjected to unsupervised hierarchical clustering (Eisen et al. 1998) using MeV 4.9 release (MultiExperiment Viewer).

Real-time RT-PCR

Total RNA was reverse transcribed enabling miRNA specific cDNA synthesis for the 10 different human miRNAs and 2 miRNA controls. RT-PCR performed on a 7500 Fast Real-Time PCR System and the comparative quantitation 2−ΔΔCT method (Livak and Schmittgen 2001) was used to compare differences in cycle number thresholds for samples normalized for endogenous controls (details provided in Supplemental Methods).

mRNA-expression profiling of the 95 samples, each with three biological replicates, was performed using the Illumina HumRef-8v2 BeadChips (Illumina, Inc.) where 200 ng of total RNA was processed according to the protocol supplied by Illumina as previously reported (Grundberg et al. 2009) The biological replicates were always separated and hybridized on different BeadChips in order to avoid that results are confounded by technical effects. Any sample or array that showed an abnormal distribution of signal intensities as visualized in box plots was removed or rehybridized. The raw data were imported to Bioconductor (Gentleman et al. 2004) using the R 2.5.0 package for variance-stabilizing transformation and robust spline normalization to obtain normalized mean-expression values.

Western blot analysis

HOBs were washed twice with ice-cold PBS and total cell lysates was prepared using RIPA lysis buffer (50 mM Tris–HCl, pH 8.0, with 150 mM NaCl, 1.0% Igepal CA-630 (NP-40), 0.5% sodium deoxycholate, 0.1% SDS, supplemented with 1.0% protease inhibitor cocktail (Sigma-Aldrich). Lysates were incubated on ice for 30 min and centrifuged at 10,000-RPM for 20 min to collect supernatant. For estimation of osteocalcin and osteonectin protein levels, HOBs were grown and in conditioned media (serum free) and after transfection, the medium was harvested to quantify the secreted protein, osteocalcin and osteonectin. Conditioned media was precipitated with 4 volume of cold acetone and resuspended in RIPA buffer. Protein concentration was quantified by using Coomassie Plus—The Better Bradford Assay Reagent (Thermo Scientific). Twenty micrograms of soluble protein was subjected to SDS-PAGE and the separated proteins were transferred to polyvinylidene difluoride membrane (Millipore). The protein bands were probed with primary antibodies against Col1A1 (1:1000 dilution, a kind gift of Dr. Larry Fisher [NIDCR, National Institutes of Health]), BGLAP (1:1000 dilution, GeneTex), SPARC (ab89218, 1:1000 dilution; Abcam), FRZB (1:1000 dilution, Sigma), RUNX2 (1:250 dilution, Sigma), and ACTB (1:1000 dilution; Cell Signaling Technology). Anti-Rabbit-HRP conjugated secondary antibodies (1:3000 dilution; R&D Systems) were used to detect the primary antibodies, followed by the target protein visualization with EMD Millipore Immobilon Western Chemiluminescent HRP Substrate (ECL). Images were acquired using LI-COR Odyssey Fc Dual-Mode Imaging system (LI-COR Biosciences) and Image Studio Software.

Data analysis

RT-PCR-expression normalization

Assessment of the stability of the miRNA gene expression was performed using softwares geNorm, version 3.5 (Vandesompele et al. 2002) and Normfinder (Andersen et al. 2004). The program geNorm is a Visual Basic application tool for Microsoft Excel. The program selects from a panel of candidate reference genes the two most stable genes or a combination of multiple stable genes for normalization (Fu et al. 2009). The NormFinder (freely available on the Internet http://www.mdl.dk) is a Microsoft Excel add-in and calculates the stability values of the individual candidate reference genes for normalization (Vandesompele et al. 2002). A low stability value indicating a low combined intra- and inter-group variation proves high expression stability. Using this approach, the most stable single gene is calculated.

Downstream data analysis

Downstream data analysis was carried out with R statistical computing framework version 2.11.1 (http://www.R-project.org) (The R Foundation for Statistical Computing 2010) using tools from Bioconductor (Gentleman et al. 2004). The association between paired miRNA and mRNA profiles was tested by Pearson correlation coefficients. False Discovery Rate (FDR) and P-values <0.05 were also computed. The R software package adjusts P-values generated in multiple hypotheses testing of gene expression data. The software applies multiple testing procedures that control the False Discovery Rate (FDR) criterion introduced by Benjamini and Hochberg (1995). For visualization of miRNA and mRNA correlations, the Statistica v10 software was used (StatSoft, Inc.).

Overrepresentation analysis of predicted miRNA targets

Predictions of genes carrying conserved miRNA target sites in their UTRs were downloaded from the Targetscan website (Lewis et al. 2005). We then further analyzed hsa-miR-29b since it had the largest number of significant correlations to mRNA levels. (details provided in Supplemental Methods).

Bioinformatic analysis of correlations—ingenuity pathway analysis

The ingenuity pathway analysis (IPA; Ingenuity Systems, www.ingenuity.com), a web-delivered application, was used to explore molecular functions and relevant networks and gene targets (Mori et al. 2009), where the significant correlations were associated and/or overrepresented (details provided in Supplemental Methods).

Statistical analysis

Correlations between miRNA fold-change values and mRNA-expression values were determined using the R statistical software 3.11.1 package in Bioconductor (Release 3.0). The Student's t-test was used to determine statistical differences between pairs of groups. Two-way analysis of variance (ANOVA) was used to evaluate the statistical significance for comparisons within groups. P < 0.05 (two sided) was considered as statistically significant. P-values obtained were adjusted for false discovery rate due to multiple testing correction (Benjamini and Hochberg 1995). For generation of correlation graphs and statistical analysis, Statistica v12 software (StatSoft, Inc.) and GraphPad Prism software package (version 6.0) were used.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.