Distinct Subsets of Noncoding RNAs Are Strongly Associated With BMD and Fracture, Studied in Weight‐Bearing and Non–Weight‐Bearing Human Bone

We investigated mechanisms resulting in low bone mineral density (BMD) and susceptibility to fracture by comparing noncoding RNAs (ncRNAs) in biopsies of non–weight‐bearing (NWB) iliac (n = 84) and weight bearing (WB) femoral (n = 18) postmenopausal bone across BMDs varying from normal (T‐score > −1.0) to osteoporotic (T‐score ≤ −2.5). Global bone ncRNA concentrations were determined by PCR and microchip analyses. Association with BMD or fracture, adjusted by age and body mass index, were calculated using linear and logistic regression and least absolute shrinkage and selection operator (Lasso) analysis. At 10% false discovery rate (FDR), 75 iliac bone ncRNAs and 94 femoral bone ncRNAs were associated with total hip BMD. Eight of the ncRNAs were common for the two sites, but five of them (miR‐484, miR‐328‐3p, miR‐27a‐5p, miR‐28‐3p, and miR‐409‐3p) correlated positively to BMD in femoral bone, but negatively in iliac bone. Of predicted pathways recognized in bone metabolism, ECM‐receptor interaction and proteoglycans in cancer emerged at both sites, whereas fatty acid metabolism and focal adhesion were only identified in iliac bone. Lasso analysis and cross‐validations identified sets of nine bone ncRNAs correlating strongly with adjusted total hip BMD in both femoral and iliac bone. Twenty‐eight iliac ncRNAs were associated with risk of fracture (FDR < 0.1). The small nucleolar RNAs, RNU44 and RNU48, have a function in stabilization of ribosomal RNAs (rRNAs), and their association with fracture and BMD suggest that aberrant processing of rRNAs may be involved in development of osteoporosis. Cis‐eQTL (expressed quantitative trait loci) analysis of the iliac bone biopsies identified two loci associated with microRNAs (miRNAs), one previously identified in a heel‐BMD genomewide association study (GWAS). In this comprehensive investigation of the skeletal genetic background in postmenopausal women, we identified functional bone ncRNAs associated to fracture and BMD, representing distinct subsets in WB and NWB skeletal sites. © 2020 The Authors. Journal of Bone and Mineral Research published by American Society for Bone and Mineral Research.

Introduction P ostmenopausal osteoporosis is the most common disease in women above 50 years of age, leading to fragile bone and fractures even after mild skeletal trauma. (1) About 40% of all US white women will suffer at least one clinically apparent fragility fracture during their lifetime. (2) The consequences of fractures, commonly occurring in the spinal vertebrae, wrists, and hip, pose a serious health risk to the patients and have a profound social and economic impact. (3) High and low bone mineral densities (BMDs) are typically associated with strong and fragile bones, respectively. Reduced BMD leads to osteoporosis, defined by the World Health Organization (WHO) as a BMD of 2.5 standard deviations (SDs) or more below the young adult mean (T-score ≤ −2.5). BMD is highly heritable. Mothers with low BMD have children with significantly lower BMD more frequently than the general female population. (4) Also, studies involving twins or kin relationships have demonstrated that individual variations in BMD are up to 80% genetically determined. (5)(6)(7) A recent genomewide association study (GWAS), (8) the most extensive so far, identified 518 loci accounting for approximately 20% of heel BMD variation. Some of the discrepancies between that GWAS and family studies may be explained by different epigenetics between osteoporotic and healthy subjects, as indicated by recent DNA methylation studies and difficulties detecting disease associated genetic variants due to minor contribution. (9)(10)(11) MicroRNAs (miRNAs) are considered part of the possible epigenetic mechanisms contributing to osteoporosis risk, and their regulatory effect are heritable. (12,13) miRNAs are 19-nucleotide (nt)-long to 24-nt-long noncoding molecules that are distinct from, but related to, small interfering RNAs (siRNAs), and have been identified in a variety of organisms. (14) They are transcribed from distinct genes and regulate the translation and/or degradation of specific messenger RNAs (mRNAs) or other transcripts, thereby modulating many critical cellular signaling pathways and functions. (14) The miRBase database (release 22; Manchester University, Manchester, UK; http:// www.mirbase.org/) contains well over 2500 mature human miRNAs. (15) Despite the critical importance of miRNAs in the regulation of various cell functions, their role in osteoporosis has not yet been investigated in depth in well-characterized human bone.
In addition to miRNAs, thousands of other noncoding RNAs (ncRNAs) are transcribed from the human genome, and may regulate gene expression at transcription, RNA processing, and translation levels. (16) These untranslated RNAs include small nucleolar RNAs (snoRNAs, SNORDs, SNORAs), most of which are pieces of introns (~70 nt) that have separate functions after being excised through exonucleolytic processing. (17) The snoR-NAs are essential for pre-ribosomal RNA (pre-rNA) processing, catalyze nucleotide modifications, and may also serve to chaperone the correct RNA fold for rRNA processing and ribosome assembly. (18) We aimed to establish a putative association of ncRNAs with BMD or fracture by performing global gene expression profiling of bone biopsies from postmenopausal white women with bone mass varying from healthy to osteoporotic with or without prior fragility fractures. mRNA levels vary considerably between weight-bearing (WB) spine and non-weight-bearing (NWB) iliac bone in the same individuals. (19,20) To make the present study representative for different skeletal sites, we included NWB iliac and WB femoral bone where the most serious fractures occur.

Ethics statement
The study was approved by the Norwegian Regional Ethical Committee (REK no 2010/2539, Norway), all volunteers gave their written informed consent, and sampling and procedures were according to the Act of Biobanking in Norway.

Participants
The cohort of iliac bone donors was enlisted after completing a questionnaire that included questions on lifestyle factors and was deemed representative of the Oslo-based Norwegian ethnic female population aged 50 to 86 years. Nonrelated postmenopausal ethnic Norwegian women (50 to 86 years, n = 300) were recruited from the outpatient clinic of Lovisenberg Diaconal Hospital (Oslo, Norway). Of these, 177 were excluded from the study due to medication or diseases other than primary osteoporosis that are known to alter bone metabolism. Altogether, 100 women were enrolled in the study, and underwent transiliac bone biopsies of the anterior superior iliac crest. Eighty-four biopsies underwent gene expression analysis. The biochemical and physiological parameters of these women have been described. (11) Moreover, site-specific BMD was evaluated with Lunar Prodigy DEXA (GE Lunar, Madison, WI, USA) following the manufacturer's instructions. The precision of the instrument for measuring the lumbar spine (L 2 -L 4 ) and hip BMD was 1.7% and 1.1%, respectively.
Femoral bone intertrochanteric biopsies were obtained from postmenopausal women (n = 18) with a wide BMD range, ie, from healthy to osteoporotic, who were undergoing hip replacement surgery ( Table 1). None of the participants had medication or diseases other than primary osteoporosis known to affect bone metabolism. Anthropometric data, including BMD, for each of the 18 donors are shown in Table S1.
Purification and quantification of ncRNAs Total RNA was isolated from bone biopsies as described (11) and quality checked using Agilent BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA).
For analysis of the iliac ncRNAs, total RNA (45 ng) was reversetranscribed using Megaplex TM reverse transcription (RT) primersets A and B (Applied Biosystems, Foster City, CA, USA). TaqMan miRNA LDA Arrays A and B (Applied Biosystems) were used for profiling 758 different ncRNAs. The data were analyzed using SDS 2.3 software (Applied Biosystems, Foster City, CA, USA) and the comparative threshold cycle (ΔΔCt) method. We used a restricted mean of expressed miRNAs, and normalized between two plates with the same sample and performed a normalization across samples.
For mature miRNA and ncRNA microchip expression profiling of femoral bone intertrochanteric RNA, 300 ng total RNA was used for biotin labeling by a Genisphere FlashTag HSR kit (Genisphere, Hatfield, PA, USA) according to the manufacturer's protocol. Labeled miRNAs were hybridized to the GeneChip miRNA 4.0 Array (Affymetrix, Santa Clara, CA, USA). Affymetrix CEL files (containing probe intensities) were imported into Partek Genomics Suite software (Partek, St. Louis, MO, USA). Robust microarray analysis (RMA) was applied for normalization.

Correlation studies
The Benjamini and Hochberg procedure (21) was used for controlling the false discovery rate (FDR). The correlation coefficient was computed between variation in the expression levels of each mature miRNA/ncRNA expressed in 95% of the biopsies and the BMD for all subjects (n = 84). For each RNA, the nullhypothesis of zero correlation was tested against the two-tailed alternative hypothesis. Correlation coefficients were computed for BMD adjusted for age (Z-score) and BMI, first by regressing the BMD Z-score on BMI values (with intercept), and then defining the residuals as the BMI-adjusted Z-scores.

Lasso analysis and cross-validation
In addition to identifying the individual transcripts showing significant correlation with BMD, we also selected a set of transcripts that together accounted for a substantial percentage of the BMD variation. The set was selected using the Lasso method for variable selection and shrinkage in regression models. (22) Cross-validation was performed to estimate the penalty parameter of the Lasso method. The Lasso analysis, including the crossvalidation step, was performed for BMD Z-scores adjusted for BMI in total hip, femoral neck, and spine (L 2 -L 4 ) using the Rpackage lars (22) (http://cran.r-project.org/web/packages/lars/ index.html) (23) and the R-package glmnet (https://cran.r-project. org/web/packages/glmnet/index.html). The R 2 value for each estimated model was also computed. To compare the performance of the selected sets of transcripts to a random selection of transcripts, the R 2 value of an ordinary least squares regression with the selected set was compared to the R 2 values of an ordinary least squares regression for each of 10,000 randomly selected sets of the same. The p value was calculated as the proportion of the random sets with a higher R 2 value than the Lassoselected set.

Calculation of fracture associated iliac bone ncRNAs
To test whether the expressions of the transcripts were related to fracture, a logistic regression model was fitted for each of the 303 transcripts from iliac bone, with fracture (yes/no) as the response variable and transcript, age, and BMI as explanatory variables. The Wald test p values were corrected for multiple testing by using the Benjamini and Hochberg procedure. (21) In this analysis, 47 iliac bone donors without fracture and 36 with previous fracture were included. Out of the 36, 23 had at least one vertebral fracture and 27 had at least one nonvertebral fracture (wrist or femoral neck). Associations with fracture were performed only for iliac bone donors because only two of the femoral bone donors did not have fracture.
Cis-eQTL (expressed quantitative trait loci) analysis of iliac bone samples Genomewide genotyping of the iliac bone donors was performed using the Affymetrix Axiom Biobank Genotyping Array (Affymetrix, Santa Clara, CA, USA) (~700,000 SNPs assessed), (26) followed by imputation to the haplotype reference panel (HRC 1.0). SNPs with minor allele frequency (MAF) > 0.05 and imputation quality (R2) > 0.3 were considered for further analysis. Combining the available genotyped patient data and iliac bone ncRNA expression data of patients yielded data of 77 patients in total. We then used a two-step approach to assess the association of SNPs with the expression of each ncRNA. First, using linear models we adjusted the ncRNA expression for either age or age and BMI, and then defined the model residuals as either age-adjusted RNA expressions or age-and BMI-adjusted RNA expression. These model residuals were scaled using zero mean unit variance standardization. Second, we performed cis-eQTL analyses, using the Rvtests software package with single variant tests, (27) for each adjusted ncRNA expression. SNPs within 500 kilobases (kb) of the ncRNA transcription start or termination sites were surveyed. Values of p of association were adjusted for multiple testing (FDR < 0.05). Verification of expression of relevant miRNAs in osteoblasts/osteocytes About 95% of bone cells constitute osteoblasts/osteocytes, but the bone biopsies also contained bone marrow, including nonbone cell types. Thus, it was necessary to examine whether the BMD associated miRNAs were expressed in osteoblasts/osteocytes. For this purpose we took advantage of data from De-Ugarte and colleagues, (28) who performed microarray profiling of primary osteoblasts obtained from bone biopsies taken during knee replacement due to osteoarthritis (n = 4, NCBI GEO repository, accession number GSE74211).

Cohorts and experimental strategy
We used bone biopsies from two cohorts of postmenopausal women (age range, 50 to 86 years), and statistical analyses were carried out as outlined in Fig. 1.

Demographic, clinical, and laboratory parameters of the iliac bone donors
Demographic, clinical, and laboratory parameters of the iliac bone donors, have been described. (11) Briefly, the mean AE SD age and BMI of the 84 subjects were 64.6 AE 9.6 years and 24.2 AE 3.6 kg/m 2 , respectively. They had a wide range of BMD of the spine (L 2 -L 4 ; T-score range: −6.2 to 2.9), femoral neck (−4.1 to 1.3), and total hip (−4.1 to 2.0), spanning osteoporotic to normal BMD. All participants had normal clinical, biochemical, and nutritional status. The mean levels of parathyroid hormone (PTH), vitamin K, Ca 2+ , phosphate, and 25(OH) vitamin D, as well as the bone remodeling parameters osteocalcin, pyridinoline cross-linked carboxyterminal telopeptide of type I collagen (1CTP) and bone-specific alkaline phosphatase were all within normal reference ranges. (11) Quantitative histomorphometry in an additional seven women (T-score total hip range: −3.3 to 1.1) showed no change in osteocyte number between normal and osteoporotic trabecular bone from iliac crest. (11) Demographic, clinical and laboratory parameters of the femoral bone donors Data on the femoral bone donors are presented in Table 1. This cohort consisted of postmenopausal women (n = 18) with a wide range in BMD, from healthy to osteoporotic, undergoing hip replacement surgery.
Identification of BMD-associated ncRNAs in iliac or femoral bone biopsies Table 2 summarizes significant correlations found between gene expression signal strength of ncRNAs originating from either iliac or femoral bone and BMD variation in the same donors. For iliac bone, 303 ncRNAs were detected in >95% of the samples and were used for all remaining statistical analyses. For femoral bone, analyzed by microarrays, ncRNA transcripts with log 2 signal <5 across all samples were omitted.
In iliac bone, when accounting for multiple testing, applying FDRs of varying stringency (1%, 5%, 10%), 36, 58, and 75 transcripts, respectively, correlated with BMI-adjusted total hip Zscore (age adjusted T-score) ( Table 2). Table 2 also shows the corresponding numbers for the femoral neck and L 2 -L 4 BMIadjusted Z-scores. In femoral bone, no ncRNAs were identified for total hip BMI-adjusted Z-score at 1% or 5% FDR, but 94 were identified at 10% FDR. Table 2 shows that in femoral bone as well as in iliac bone, more transcripts were correlated to BMD at the total hip as compared to BMD at the femoral neck and the lumbar regions after adjustment for FDR. The correlated transcripts identified for the different skeletal sites at ≤10% FDR are presented in Table S2 (iliac bone) and Table S3 (femoral bone). The relative amounts of the 36 iliac ncRNAs at 1% FDR are presented in Fig. S1.

Lasso analysis
Lasso analysis is a regression-based analytical method commonly used when variables outnumber samples. (29,30) We used Lasso analysis as a separate independent statistical method on both cohorts of bone donors to identify sets among all detected transcripts that would explain a high proportion of BMD. The Lasso analysis selected nine ncRNAs in both iliac and femoral bone. The degree of explained total hip BMD variation after adjustment for age and BMI was calculated from an ordinary least squares regression model with the selected ncRNAs, and was 53% (R 2 = 0.53, p = 1.0E-4) and 88% (R 2 = 0.88, p = 3.0E-5) for iliac bone and femoral bone, respectively ( Table 3). Absence of ncRNAs selected for femoral ncRNAs associated with total femoral neck and spine is probably due to weaker associations and few samples, because the number of selected variables depends on sample size.
Because it is possible that the Lasso-selected sets of ncRNAs could explain a large proportion of the variation in BMD just by chance, we compared the results (R 2 = 0.53 and R 2 = 0.88) to the R 2 of the ordinary least squares regressions for 10,000 randomly selected sets of ncRNAs with set sizes of nine. On average, the random sets resulted in R 2 = 0.26 and R 2 = 0.74 for iliac bone and femoral bone, respectively, which is clearly lower than what was found as variance explained of total hip BMD for the initial The table shows number of BMD correlated ncRNAs derived from iliac bone (A) and femoral bone (B) after correction for variation in age and BMI at indicated sites with unadjusted p value and with different levels of FDR. BMI = body mass index.  R 2 values were calculated for ordinary least squares regression models as described in Subjects and Methods. In parentheses, the number of ncRNAs selected by the Lasso analysis and used for calculation of R 2 values. As control, average R 2 and p values were calculated from 10,000 random sets of the same size as those selected by Lasso for the respective anatomical sites. Lasso did not select femoral ncRNAs associated with total femoral neck nor spine. This is probably due to weaker associations and because of few samples as the number of selected variables is dependent on sample size.
Lasso-selected sets. None of the random sets for iliac bone (total hip, total femoral neck, L 2 -L 4 ) had R 2 as large as the sets identified by Lasso analysis.
Seven of the nine Lasso-identified iliac bone ncRNAs were among the 36 that were significantly correlated (<1% FDR) to adjusted total hip BMD (Tables S2 and S4A). Furthermore, eight of the nine Lasso-identified femoral ncRNAs were among the set of 94 ncRNAs correlating (<10% FDR) to total hip BMD adjusted Z-score (Tables S3 and S4B). Table S4A lists all of the selected Lasso-identified iliac ncRNAs for age-and BMI-adjusted BMD for total hip (n = 9), total femoral neck (n = 7), and L 2 -L 4 (n = 33). Table S4B shows the corresponding Lasso-identified femoral bone ncRNAs. Notably, Lasso identified no common ncRNAs between iliac and femoral bone.
Femoral bone BMD-associated miRNAs were also verified by PCR (Table S5), showing correlation r values in the same direction as found in Lasso analysis, although with p values >.05, possibly due to the limited number of samples. Although the nine femoral ncRNAs identified by Lasso and microarray analysis could explain 88% (R 2 = 0.88) of the variation in adjusted total hip BMD, PCR and multiple linear regression analyses of seven of these femoral bone ncRNAs could explain 44% (R 2 = 0.44) of the variation in the adjusted total hip BMD, thus supporting microarray and Lasso findings. Relative levels of the PCR tested miRNAs are presented in Fig. S3.

Identification of fracture associated ncRNAs
Using a logistic regression model (Subjects and Methods) we identified 29 ncRNAs to be associated with any previous fracture in the cohort of iliac bone donors at adjusted p < .05. The topmost fracture associated ncRNAs (adjusted p < .1) are listed in Table S6. Out of the 36 women with fracture 23 had at least one vertebral compression fracture.
Of the 29 fracture-associated iliac miRNAs, only miR-16-1-3p and miR-23a-3p were not among the ncRNAs correlated to BMD in Table S2. To test whether they were associated with hip bone geometry they were correlated to hip structure parameters in n = 80 of the same bone donors. Hip bone geometry in these donors has been described. (31) Interestingly, miR-16-1-3p showed nominally significant correlation with neck shaft angle (r = −0.22, p = .050) and section modulus of the femoral neck (r = 0.22, p = .050).
The PCA plot (Fig. 3A) is based on all detected iliac ncRNAs, and indicates that iliac bone donors with fracture has a common ncRNA expression profile that at least partly separates them from donors without fracture. The woman representing the leftmost blue outlier ball had nonvertebral as well as vertebral fractures, thus suggesting an underlying disease, although all biochemical markers were normal. The PCA plot in Fig. 3B is also based on all iliac ncRNAs and shows grouping of donors with T-score ≤ −2.5 (blue balls) as compared to donors with normal BMD (Tscore > −1, red balls).

Functional annotation clustering of predicted miRNA targets in iliac and femoral bone
Gene enrichment analysis showed that the majority of the identified pathways are specific to either iliac or femoral bone ( Table 4). As expected, extracellular matrix (ECM)-receptor interaction, being central in bone metabolism, was the second most significantly enriched pathway at both skeletal sites. Moreover, fatty acid metabolism was more prominent in iliac bone than in femoral bone.
Cis-eQTL analysis of iliac bone samples DNA from the iliac bone biopsies has previously undergone gene profiling using the Affymetrix UK Biobank Axiom Array. (26) To evaluate whether SNPs within a 1-megabyte (MB) window of the ncRNA-encoded regions influenced the ncRNAs levels, we performed two cis-eQTL analyses, adjusting for age and for age and BMI (Subjects and Methods). At q < 0.05, rs12187909 was associated with miR-449b-5p (p = 4.14E-5) and rs74918612 associated with miR-331-3 (p = 7.59E-4). The former came out as suggestively correlated when adjusting ncRNA expressions only for age, and the latter when adjusting ncRNA expressions for age and BMI (Table 5). Tables S7 to S11 include detailed results from the cis-eQTL analysis, including all ncRNA-associated SNPs adjusted for age or age and BMI with q values. These tables also present p values for the same SNPs from the GWAS by Morris and colleagues (8) on estimated heel bone mineral density (eBMD). Figures S4 and S5 show the squared coefficients of correlation (R 2 ) between different SNPs, as curated from the cis-eQTL analysis results with age as predictor, where their respective q values were both <0.1 for miR-449b-5p and miR-331-3p.

Discussion
We present the hitherto most comprehensive study characterizing functional ncRNA expression patterns in NWB iliac and WB femoral bone in well-characterized cohorts of osteoporotic and healthy postmenopausal white women. Comparisons between the BMD-correlated ncRNAs expressed at the two sites revealed major differences with minor overlap and distinct, overrepresented signaling pathways. Lasso analysis of the bone ncRNAs enabled the identification of small sets of ncRNAs explaining a significant proportion of BMD variation measured at different skeletal site. We also identified miRNAs associated with fracture or SNPs.
The common pathways affected by BMD-associated miRNAs in both iliac and femoral bone include that for ECM-receptor interaction and proteoglycans in cancer, both of which are obviously important in all types of bone. The three other common pathways, Amoebiasis, Glioma, and Lysine degradation, have no obvious function in bone and may be the result of miRNA targets being common to several functional pathways. For example, the glioma pathway includes miRNA targets such as mammalian target of rapamycin (MTOR), platelet-derived growth factor receptor alpha (PDGFRA), and epidermal growth factor (EGF), which have established functions in bone metabolism. (32)(33)(34) The ECM-receptor interaction pathway involves interactions between various cell-associated integrins  and matrix collagens which are central actors in bone metabolism. This pathway overlaps with the focal adhesion pathway, ranked tenth in iliac bone but ranked as not significant in femoral bone. Notably, the focal adhesion pathway is linked to the Wnt pathway, which is central to bone metabolism. (35) Other iliac bone pathways related to bone matrix include that for mucin O-glycan biosynthesis and glycosaminoglycan biosynthesis, which are of obvious relevance to BMD. The pathways specific to iliac bone also include that for fatty acid biosynthesis and fatty acid metabolism. It is well known that bone marrow stem cells (BMSCs) tend to differentiate toward adipocytes rather than osteoblasts in osteoporotic bone, (36) and an inverse association between BMD and marrow adipose tissue has been demonstrated in different populations. (37) A larger study with higher power might be able to detect these pathways in WB femoral bone.
It is unclear how several of the identified pathways targeted specifically by femoral bone miRNAs are associated with BMD, with the exception of the Hippo signaling pathway. The Hippo signaling pathway is important for stem cell function, bone development, and bone remodeling via interactions with, eg, the transforming growth factor beta (TGF-β)/SMAD pathway, Wnt signaling, and the transcription factor RUNX2 (RUNX family transcription factor 2). (38) Supporting the relevance of these findings in bone tissue, 22 of 46 ncRNAs have been shown to be important for bone cell differentiation and function in cell and animal models (Table S5). Figure S2 shows that some ncRNAs were positively correlated to BMD when expressed in WB bone, but inversely when expressed in NWB bone. Transcript levels have previously been shown to vary widely between WB and NWB bone. For example, the transcription factor ZIC1 was shown to be expressed 200-fold higher in male bone biopsies from WB lumbar spine as compared to NWB iliac crest in the same donors. (20) In the same study, more than 4000 genes were identified as differentially expressed between the two sites. Thus, the observed inverse associations may reflect another aspect of differential bone metabolism in WB and NWB bone.
We searched the literature to identify whether the total hip BMD associated ncRNAs from iliac bone and femoral bone have been experimentally shown to impact bone metabolism. The relevant studies are summarized in Table S12. Briefly, 22 of the 46 ncRNAs have been associated with bone metabolism and/or osteogenic cell differentiation; eg, miR-92a-1-5p influences osteogenic differentiation in vitro by regulating β-catenin (39) and reduces body and skull length in miR-92a −/− mice (40) and inhibition of miR-92a-1-5p enhances fracture healing. (41) The miR-15a-5p is part of the deleted in lymphocytic leukemia 2 (DLEU2) gene, which is highly correlated to BMD in postmenopausal women (11) and mice (own results unpublished). In addition, miR-149-5p may be of clinical importance, because polymorphisms within its encoded primary RNA sequence have been associated with vertebral fractures in postmenopausal Korean women. (42) Most of the miRNAs in Table S7 affect functions attributed to bone cell differentiation and mineralization involving both osteoblast and osteoclast lineages confirming their skeletal functionality.
The bone biopsies contain also non-bone cells, and it may be argued that some of the BMD-associated ncRNAs are not originating from bone cells. However, all but five of the 36 topmost femoral bone miRNAs significant at 10% FDR (Table S3) when correlated to BMI=adjusted total hip Z-score were also produced in cultured differentiated normal human osteoblasts (publicly available data described by De-Ugarte and colleagues (28) ; not shown). Notably, the five miRNAs not found were not covered by the GPL20999 miRCURY LNA miRNA Array used to profile the osteoblast miRNAs.
Notably, separate sets of small nucleolar RNAs (SNORAS, SNORDS, or snRNAs) from iliac and femoral bone were correlated with BMD. In iliac bone, SNORD44 and SNORD48, which direct 2 0 O-ribose methylation of 18S rRNA and 28S rRNA, respectively, (43)(44)(45) were negatively correlated with BMD, whereas 28SrRNA was positively correlated with BMD (Table S2). Interestingly, 28SrRNA showed an inverse correlation with SNORD44 and SNORD48 (r = −0.57 and r = −0.81; p = 1.2E−8 and p = 1.0E−15, respectively). rRNA methylation is important for rRNA processing and function, and failure in rRNA methylation has been linked to several diseases. (46) The connection between BMD and rRNA levels is supported by the identification of functional RUNX2 binding sites within DNA repeats encoding rRNA (47) and the fact that RUNX2 suppresses rRNA gene transcription during osteoblast lineage progression. (48,49) In iliac bone, 28S rRNA was significantly correlated to nearly all BMD-correlated miRNAs at 1% FDR (not shown). This may be due to the expansion segments of 28S rRNA, which have been suggested to constitute an important component of miRNA balance by binding and decreasing ("sponging") the availability of GC-rich miRNAs and thereby aid the conservation of GC-rich mRNAs. (50) Iliac RNU11 (RNA, U11 small nuclear) was positively correlated to BMD and has functions such as 5 0 splice site recognition at constitutive splicing, activation of U2-dependent alternative splicing, and regulating the U12-dependent spliceosome, (51) indicating that altered splicing may be an important mechanism affected in osteoporosis.
In femoral bone at 10% FDR, another set of total hip BMDcorrelated nucleolar RNAs was found, consisting of seven SNORDS catalyzing 2 0 -O-ribose methylation of rRNA and five SNORAS guiding rRNA pseudourididylation (52) (Table S3). The majority of snoRNAs are encoded in the introns of host genes, (53) and their expression levels are thus correlated with their host gene transcripts and regulated similarly. Furthermore, crosslinking studies have indicated that SNORDS/SNORAS can also modify mRNAs, (54) thereby possibly affecting the level of bone cell transcripts related to skeletal metabolism.
Of the 29 fracture-associated ncRNAs (Table S6), it is striking that all but one (miR-181a-5p) increase risk of fracture at increased expression levels. Furthermore, miR-16-1-3p (adj. p = .0356) and miR-23a-3p (adj. p = .0745) were not among any of the ncRNAs correlated to BMD in Table S2, and may therefore be particularly associated with structure, as also indicated by the association of miR-16-1-3p with hip structure.
Because of the limited availability of human bone biopsies, only a few similar, smaller ncRNA profiling studies have been performed, and none has been performed on NWB iliac bone. In a smaller study De-Ugarte and colleagues (28) determined that miR-320a and miR-483-5p are increased in osteoporotic bone. miR-320a was replicated in the present study, whereas miR-483-5p was detected but did not reach significance. A more comprehensive study by Seeliger and colleagues (91) indicated that miR-21, miR-23a, miR-24, miR-25, miR-100, and miR-125b are upregulated in osteoporosis using unadjusted statistics. They were also detected in the femoral bone biopsies in the present study, but did not reach significance when correlated to BMD adjusted for variation in age and BMI. Garmilla-Ezquerra and colleagues (92) compared expression of 760 miRNAs in explanted osteoblasts from the femoral heads from patients with osteoporotic hip fractures and used patients with osteoarthritis as surrogate controls. They found that miR-187 was increased in the controls, whereas miR-518f was increased in the osteoporosis group. Neither were identified in our femoral bone biopsies using multivariate analyses.
It is tempting to speculate that mechanically sensitive miRNAs would be associated with BMD and/or fracture, especially in WB femoral bone, but very few of those identified in cell culture studies were associated with BMD or fracture in the current study. A recent study (93) identified 10 miRNAs that changed expression when mouse MLOY4 osteocytes were subjected to mechanical strain. However, only two of the identified miRNAs were identified as BMD associated in our study. Of those, miR-29b-3p correlated positively to BMI-and age-adjusted total hip BMD when expressed in both femoral and iliac bone (Fig. S2) but showed reduced expression in MLOY4 cells subjected to mechanical stress. In the same study (93) miR-574-3p showed increased expression at mechanical stress, but this miRNA was inversely correlated to BMD when expressed in iliac bone (Table S2). Thus, it is hard to conclude on the association between mechanically sensitive miR-NAs and their levels in WB and NWB bone.
Cis-eQTL analyses of iliac bone identified two loci/miRNAs of interest. The locus near the miR-331-3p coding sequence was also identified in the hitherto most comprehensive GWAS on heel eBMD. (8) However, none of the SNPs identified in the present study were associated with eBMD in the GWAS by Morris and colleagues, (8) nor were they in LD with any of the eBMDassociated variants.
Morris and colleagues (8) did not attempt to explain the effect of genomic variation on miRNA expressions but our findings may offer an explanation. Given the biological background of miRNA mechanism of action, an argument can be given for such a nondirect association of miRNA with eBMD variants genomewide based on the present results. Polymorphisms influencing miRNA variation are not influencing BMD variation directly due to epigenetic regulation. It would be desirable to perform a replicate study with a larger cohort, but acquisition of additional biological samples from well-characterized participants will be difficult to obtain.
Most of the BMD-correlated ncRNAs in iliac bone, and in femoral bone, represent novel bone biology. The use of nontraditional statistical methods, such as Lasso, allowed us to select ncRNA sets maximizing the BMD variance explained. The power of this technique was further confirmed by a permutation analysis. This and our previous study (11) confine the large number of transcripts expressed in bone biopsies (summing thousands) to a nadir of nine ncRNAs identified through the Lasso analysis and a limited set of mRNAs, which statistically are strongly associated with BMD. These molecular components in bone, as presented in Tables 3 and S4, admittedly at the risk of oversimplification, present an initial schematic outline of an emerging picture of ncRNAs having important roles for normal bone remodeling in NWB iliac bone and WB femoral bone.
Our study has some limitations. The number of participants is limited and the study includes only postmenopausal women of Northwestern European background. Especially for femoral data, power is a limiting factor in terms of the detected number of BMD associated ncRNAs. With a larger cohort, power would be stronger and more BMD associated ncRNA would probably be detected. Nevertheless, our rigorous statistical analyses ensures validity of our findings. The bone biopsies contain marrow, hence the results may be somewhat influenced by transcripts in this tissue. RNA sequencing of the samples indicated that marrow cells of hematopoietic and immunological lineage constituted a minor part (about 20%) of the total cell population contained (not shown). Furthermore, it is expected that marrow constituents across iliac and femoral bone biopsies would be similar. Thus, the detected differences in ncRNA expression between these two sites are most likely attributable to differences in bone cell expression levels. Twenty-two of 46 ncRNAs have been ascribed functional bone cell significance in different experimental systems. Also, the most significantly BMDcorrelated miRNAs have been found to be expressed in primary cultured human osteoblasts. We used different platforms for profiling the iliac and femoral bone ncRNAs, but all ncRNAs profiled in iliac bone were covered by microarrays used to profile femoral bone ncRNAs, thus supporting the validity of our findings. Still, we acknowledge that detection methodologies (microarrays with detection of biotin labeled ncRNAs versus TaqMan assays involving detection of PCR products after cleavage of quenching hybridizing probes) and quantification analyses (eg, different types of sample normalization, software, and thresholds) could have an influence on the results. However, several papers show that in general there is a very good concurrence in studies where miRNA levels was measured by both GeneChip miRNA 4.0 Array and TaqMan PCR analysis. (94)(95)(96) The results invite to a more focused bone transcriptome and functional analyses and thus represent a promising foundation for pursuing studies of miRNAs and other ncRNAs preparing the ground for future gene therapy of a common and devastating disease.

Disclosures
All authors have no conflicts of interest to declare. Jahre Foundation for the Promotion of Science; the Lovisenberg Diaconal Hospital research fund; and the Research Council of Norway. This work is part of the European Union project OSTEO-GENE (no. FP6-502491). CM-G and FR are supported by the Netherlands Organization for Health Research and Development (ZonMw VIDI 016.136.367). The present work is dedicated all the Norwegian women who voluntarily participated and made this study possible. We are indebted to the Lovisenberg Diaconal Hospital, Department of Clinical Chemistry for blood sampling and to the hospital's medical outpatient clinic, for access to the surgical facility and patient care assistance. We also acknowledge Diakonhjemmet Hospital, Oslo, for assistance in the collection of femoral bone biopsies. Bioengineer Anne Runningen at The Osteoporosis Laboratory, Lovisenberg Diaconal Hospital, is especially gratefully acknowledged. All primary data are freely available from Dryad (https://datadryad.org/).
Authors' roles: KMG and SR designed the study. KMG, VTG, HV, EL, and SR were responsible for patient inclusion and collection of bone biopsies. OKO, LR, VTG, and SR were involved in RNA purification and quantification of bone ncRNAs. CM-C, ES, VP, and FR contributed to the eQTL analysis. C-CG, VP, MH, and SR performed the biostatistical analyses. C-CG, MY, TPU, SR, VP, and KMG contributed to the writing of the manuscript. All authors reviewed the manuscript, added appropriate revisions, agreed to submission for publication, and approved the final version. The corresponding authors attest that all listed authors meet authorship criteria and that no others meeting the criteria have been omitted.