SNPs in bone-related miRNAs are associated with the osteoporotic phenotype

Biogenesis and function of microRNAs can be influenced by genetic variants in the pri-miRNA sequences leading to phenotypic variability. This study aims to identify single nucleotide polymorphisms (SNPs) affecting the expression levels of bone-related mature microRNAs and thus, triggering an osteoporotic phenotype. An association analysis of SNPs located in pri-miRNA sequences with bone mineral density (BMD) was performed in the OSTEOMED2 cohort (n = 2183). Functional studies were performed for assessing the role of BMD-associated miRNAs in bone cells. Two SNPs, rs6430498 in the miR-3679 and rs12512664 in the miR-4274, were significantly associated with femoral neck BMD. Further, we measured these BMD-associated microRNAs in trabecular bone from osteoporotic hip fractures comparing to non-osteoporotic bone by qPCR. Both microRNAs were found overexpressed in fractured bone. Increased matrix mineralization was observed after miR-3679-3p inhibition in human osteoblastic cells. Finally, genotypes of rs6430498 and rs12512664 were correlated with expression levels of miR-3679 and miR-4274, respectively, in osteoblasts. In both cases, the allele that generated higher microRNA expression levels was associated with lower BMD values. In conclusion, two osteoblast-expressed microRNAs, miR-3679 and miR-4274, were associated with BMD; their overexpression could contribute to the osteoporotic phenotype. These findings open new areas for the study of bone disorders.

miRNAs have been extensively studied in bone research, particularly their relationship to osteoporosis [1][2][3] . These studies clearly showed altered miRNAs profiling in serum from patients with osteoporosis, as well as in bone tissue after osteoporotic (OP) fracture. However, these miRNA expression signatures observed in patients with osteoporosis do not provide evidence of causality because the altered pattern could be a consequence of the disease or even unrelated to the pathogenesis.
Another approach in miRNAs studies is the association analysis between one SNP within a candidate miRNA (miR-SNP) or in a miRNA target site and one disease related-outcome. In this case, the associated variant is likely involved in the pathophysiology or confers susceptibility to develop the disease. Moreover, many evidences suggest that the genetics of complex traits are attributable to genetic variations that modulate gene expression, rather than the variations resulting in protein structure changes 4 . However, functional assays are needed in order to elucidate the role of the associated variants in the pathophysiology of the disease since the SNP could be in linkage disequilibrium with the true functional variant.
The aim of this study was to identify SNPs within candidate miRNAs in order to perform an association study between those SNPs and bone mineral density (BMD), the main outcome used to define osteoporosis. First, we searched for miR-SNPs in the primary miRNA transcript (pri-miRNA), which has a hairpin structure with a terminal loop and two single-stranded flanking regions. The pri-miRNAs are recognized and cleaved by the Drosha and DCGR8 complex, resulting in a shorter structure called pre-miRNA 5 . In this step, the sequences at the unpaired flanking arms and within the hairpin double-stranded stem structure are crucial to correct binding and cleavage by the Drosha-DCGR8 complex. Thus, the existence of genetic variants within the pri-miRNA sequences could lead to an alteration of the hairpin structure, affecting molecular processing and the underlying miRNA maturation 6 . Changes in miRNA maturation would trigger changes in miRNA abundancy, and consequently a deregulation of the expression levels of target genes. Supporting this idea, large-scale in silico analyses of SNPs in human miRNA genes have demonstrated lower SNP densities in the miRNA sequence than their flanking regions or the human genome 7 . Hence, our study was based on the detection and subsequent genetic association analysis of putative functional miR-SNPs. Furthermore, associated miR-SNPs were explored in bone cells in order to validate the association with the OP phenotype.

Results
An overall overview of methodology and results is schematized in Fig. 1.
Association analysis with BMD. The first approach used in our study was to identify functional variants within microRNAs involved in bone metabolism. The minor allele frequency (MAF) of many of the miR-SNPs found in databases has not been assessed in European population; therefore, a validation step was necessary. In this regard, SNPs in microRNA sequences were firstly validated as polymorphic in patients from BARCOS subcohort 8 . Of the 53 miR-SNPs tested, 14 were validated in this cohort. However, only 5 SNPs had a MAF > 0.01 and were finally genotyped in the OSTEOMED2 cohort for their association analysis with LS (lumbar spine) BMD and FN (femoral neck) BMD (Table 1) Quantification of miRNA expression levels in total FN bone samples. The anthropometric features of the OP and Control groups were shown in Table 2. Using the Mann-Whitney U test, no statistical differences in any features were observed between both groups.  MiR-3679 and miR-4274, which harbored BMD-associated SNPs in their pre-miRNA sequence (Fig. 2), were quantified by qPCR in order to compare the expression levels between OP and non-OP bone samples.
MiR-3679-5p was not detected in total bone samples or in cultured human osteoblast cells, but miR-3679-3p was present in both sample types. Therefore, only mature miR-3679-3p was assayed in the qPCR. Both miRNAs miR-3679-3p and miR-4274 were significantly overexpressed in the OP samples (Table 3).
Other bone miRNAs, miR-631 and miR-574-5p 1,9 , were also tested in order to preclude a non-specific phenomenon related to a general upregulation of gene expression. No differences were found between osteoporotic and non-osteoporotic bone samples (p = 0.626 and p = 0.183 respectively).
Association analysis of miRNA expression levels with genotypes of BMD-associated SNPs. Human primary osteoblasts (n = 38) were cultured for DNA and RNA extraction and sorted by   Table 3. miRNA expression levels, comparison between osteoporotic and non-osteoporotic bone samples. Abbreviations: RQ = Relative quantification; IQR = Interquartile Range.
genotype for both rs6430498 and rs12512664. Association between expression levels of mature miRNA miR-3679-3p and miR-4274 with its corresponding own genotype was analyzed. A significant correlation was observed between miRNA levels and the genetic variant (Fig. 3). The A allele for both SNPs was associated with higher expression of each corresponding miRNA (miR-3674; log-additive model; p value = 0.015; recessive model; p value = 0.03, and miR-4274; recessive model; p value = 0.013). Additionally, in order to corroborate that the differences among expression levels are due to genotypes "per se" and not for other cellular circumstances, other bone-related miRNAs 1,9 were checked in these cells. Expression levels of miR-320a, and miR-22-3p were measured and no differences were found irrespective of genotypes. Moreover, as a sensitivity test, we analyzed the miR-3679-3p and miR-4274 expression levels separated according to the exchanged genotype amongst both miRNAs (rs12512664 and rs6430498 respectively). Once again, no differences were found in both miRNA analyzed (miR-3679-3p; rs12512664: p-value = 0.713 and miR-4274; rs6430498: p-value = 0.645).

Bioinformatic analysis.
A subtle effect of the BMD-associated SNPs on secondary pre-miRNA structure was detected. The genetic variants rs6430498 and rs12512664 did not provoke evident changes in the loop formation on the RNAstructure web server (Fig. 2).
In-vitro assessment of osteoblast activity and matrix mineralization. hOBs were treated with miR-3679-3p and miR-4274 mimics or inhibitors and their respective controls for investigating their role in osteoblast function (Supplemental Fig. 1).
For miR-4274 transfections, no changes in any parameter evaluated were found.

Discussion
Two osteoblast-related microRNAs, miR-3679 and miR-4274, were found to be associated with femoral neck BMD and overexpressed in OP fractured bone. These miRNAs harbor genetic variants in their pre-miRNA sequences that were correlated with the miRNA expression levels in human osteoblasts. Finally, inhibition of miR-3679-3p in human osteoblastic cells increased matrix mineralization. The lower density of genetic variants in miRNA genes compared to other non-coding genomic regions suggests that SNPs can have a remarkable biological role in pri-miRNAs by differential regulation of their target genes 10 . Hence, we chose genetic variants within the pri-miRNA sequences for their putative functionality. Moreover, the BMD-associated miR-SNPs in our study were located in the pre-miRNA molecules, which have been described to exhibit strong selective constraints 11,12 . SNPs within the pri-miRNA sequence can affect the miRNA maturation in multiple aspects, modifying the hairpin structure or changing the binding affinity of biogenesis enzymes to the miRNA hairpin. Moreover, variants can affect alternative cleavage sites of biogenesis enzymes, producing novel isomiRs or changing their frequency. These changes in maturation can result in altered expression of the functional miRNA, resulting in deregulation of target genes. Hence, several studies have suggested that functional SNPs in pri-miRNAs affect the processing and expression levels of mature miRNAs 13 .
Two mir-SNPs, rs6430498 and rs12512664, within the pre-miRNA sequence of miR-3679 and miR-4234, respectively, were associated with FN BMD in the OSTEOMED2 cohort. Bioinformatic analysis of the hairpin structure of these BMD-associated miRNAs showed no striking effect on its secondary structure caused by allelic changes but we cannot rule out that the expression was altered by affecting the affinity binding of the processing enzymes or their accessory proteins. In this regard, other studies have demonstrated that miR-SNPs in pre-miRNA sequences lead to altered levels of mature miRNA in vivo even though changes in its secondary structure are not predicted by RNAfold program 14 .
In our study the AA genotype for both rs6430498 and rs12512664 showed higher expression levels of miR-3679-3p and miR-4274 in hOBs, respectively. Moreover, these miRNAs were overexpressed in OP bone samples. Consistently with miR-3679-3p overexpression in OP fractured bone, transfection results suggest that miR-3679-3p negatively regulates matrix mineralization by osteoblasts which could lead to an OP phenotype. Further supporting these results, the A allele for both SNPs was associated with lower FN BMD levels.
In overall, predicted targets and pathways involving miR-3679-3p belong to cell cycle, differentiation and tissue development, suggesting a dysfunction in bone remodeling in the osteoporotic bone. In agreement with these data, an impaired mineralization was observed when miR-3679-3p was overexpressed in osteoblasts. In contrast, analysis of targeted genes and pathways for miR-4274 suggests that this miRNA could be involved, broadly, in vesicular and membrane trafficking, cell migration and adhesion, and cell metabolism. That would explain why changes in the assayed parameters were not detected after miR-4274 transfection in osteoblasts.
Altogether, these findings suggest that these variants can play a role in bone metabolism. Since information about these two miRNAs is very scarce, this study could be considered as a preliminary approach presenting miR-3679-3p and miR-4274 within the complex network of bone regulation.
Our study is strengthened by the extremely careful in control of potentially confounding characteristics among bone samples in terms of age, sex, BMI and metabolic diseases highly prevalent in elderly population 15,16 . Moreover, we excluded from the study samples from patients treated with hormones and other anti-osteoporotic drugs that could alter the miRNA expression in bone cells 17 . These strict inclusion criteria support that miRNAs detected in our study can be more reliably considered involved in the osteoporotic phenotype.
One limitation is the use of osteoarthritic samples as control group. Due to obvious ethical reasons the collection of bone from healthy individuals is not allowed. However, in an attempt to minimize this potential limitation, we obtained the bone samples from a location distant from the interface between bone and cartilage and, therefore, as far away as possible from the osteoarthritic lesion. Moreover, although both miRNAs have found overexpressed in OP bone tissue, it cannot be attributed to the SNP regulation in these samples since this could not be analyzed due to lack of bone sample's availability. Other limitation of the study is the need to replicate the association results in other cohorts because these miR-SNPs lacked linkage information with other tag-SNPs related to bone phenotypes. However, this is the first time that two genetic variants in human osteoblast-related miRNAs have been associated with FN BMD in an extensive cohort of postmenopausal women, and this association was functionally demonstrated in bone samples.
In conclusion, two putative functional SNPs in the pre-miRNA sequence of miR-3679 and miR-4274 were associated with femoral neck BMD. In both cases, the allele associated with lower BMD was correlated with higher expression levels of mature miRNA in human osteoblastic cells. Moreover, both miRNAs were overexpressed in OP fractured bone. Our results open new exploratory avenues for future studies in the bone field.

Methods
Ethics Statement. The study protocols for obtaining DNA from blood samples, bone tissue samples and primary osteoblasts were approved by the CEIC-Parc de Salut MAR, the coordinating center (Registry number 2010/3882/I). Methods were carried out in accordance with the approved guidelines by the CEIC-Parc de Salut MAR. The approved protocols were explained to potential study participants, who provided written informed consent before being included in the study. All studies were carried out in accordance with the principles of the Declaration of Helsinki as revised in 2008.

Study subjects.
Genetic association analysis was performed in an observational, clinical cohort study (OSTEOMED2) of patients recruited from 14 medical centers in several regions of Spain. All patients were consecutive, unselected, postmenopausal women attended in an outpatient clinic. Patients were prospectively recruited regardless of BMD value (Table 4). Exclusion criteria were any history of metabolic or endocrine disease, chronic renal failure, chronic liver disease, malignancy (except superficial skin cancer), Paget's disease of bone, malabsorption syndrome, and any anti-OP or bone-affecting treatment. In addition, women with early menopause (before the age of 40) were excluded from analysis. BMD (g/cm 2 ) was measured at the lumbar spine (LS) L2-L4 and at the non-dominant femoral neck (FN) using the dual-energy X-ray densitometers available in each participating center (Supplemental Table 1).
Whole bone samples for qPCR miRNA quantification were obtained from the transcervical region of the femoral neck of postmenopausal women undergoing hip replacement due to either OP fracture (n = 10) or osteoarthritis in the absence of osteoporosis (n = 10) (T-score measurements [mean ± SD]: 0.616 ± 0.523) ( Table 2). Bone was taken from a location distant from the interface between bone and cartilage and, therefore, as far away as possible from the osteoarthritic lesion. Again, none of the participants had a history of metabolic or endocrine disease, chronic renal failure, chronic liver disease, malignancy, Paget's disease of bone, malabsorption syndrome, hormone replacement therapy, anti-resorptive or anabolic agents, oral corticosteroids, anti-epileptic drugs, or treatment with lithium, heparin, or warfarin.
Only those SNPs with a minor allele frequency (MAF) in Utah residents with ancestry from northern and western Europe (CEU) > 0.01 were included in the study; SNPs lacking published MAF were validated in participants from the BARCOS cohort 8 (included in OSTEOMED2) by means of PCR-RFLP, genotyping or Sanger sequencing.
Polymorphism genotyping. DNA was obtained from peripheral blood collected in EDTA tubes and genotyping was performed in LGC Genomics platform using KASPar v4.0 genotyping systems. To ensure genotyping quality, a random sample (5% of the total number of samples) was also genotyped in a separate control plate. There was 100% concordance between these results. Human osteoblasts culture. Human primary osteoblasts (hOB) were obtained from trabecular bone of patients who underwent knee or hip replacement due to osteoarthritis. Bony tissue was cut up into small pieces,   Table 4. Baseline characteristics of the OSTEOMED2 cohort. Abbreviations: BMI = body mass index; BMD = bone mineral density; LS = lumbar spine; FN = femoral neck.
For RNA extraction of total bone, a piece of tissue was cut out and added to QIAzol Lysis Reagent (Qiagen), then homogenized for 5 min using the TissueLyser system. Chloroform was added to each sample, followed by centrifugation for 15 min (12000 g). The upper water phase was collected and the extraction continued according to manufacturer's instructions. The quality of the total RNA was verified by an Agilent 2100 Bioanalyzer profile. The concentration of the purified RNA was analyzed on a spectrophotometer (Nanodrop, Thermo Fisher Scientific Inc). miRNA quantification by qPCR. Using the miScript II RT kit (Qiagen), 1 µg of total RNA was reverse-transcribed in 20 μl reactions. cDNA from total bone was diluted x8 and cDNA from hOBs was diluted x16; in both cases, 2 µl were assayed in 10 µl PCR reactions in 384-well plates using MiScript Syber Green PCR kit according to the protocol. The sequence of the mature miRNAs selected, according to the mirBase web site, was used as a forward primer and the Universal primer as a reverse. U6 snRNA was used as the reference gene for normalization. Amplification was performed in a QuantStudio 12 K Flex Real-Time PCR (Applied Biosystems), and the Expression Suite software was used both for determination of relative quantification (RQ) (by 2 −ΔΔCt method) and for melting curve analysis.
Bioinformatics analyses of BMD-associated miRNAs. Gene target prediction was assessed using the following six programs: PicTar (http://pictar.mdc-berlin.de), TargetScan Human (http://www.targetscan. org), miRDB (http://mirdb.org), MiRanda (http://www.microrna.org), DIANA-TarBase (http://diana.imis. athena-innovation.gr), and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw). The DIANA-mirPath web-based computational tool 26 was used to identify molecular pathways potentially altered by the intersection of miRNAs differentially expressed in fractured bone. miRNASNP database 10 , RNAstructure (http://rna.urmc.rochester.edu/RNAstructureWeb) and RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) were used to predict the effect of the variants on the miRNA secondary structure. These predictions are based on the assumption that variants can destabilize the hairpin and therefore reduce the mature miRNA levels. Alizarin Red quantification. To induce osteoblastic mineralization, hOBs (n = 4) were cultured during 28 days with hOB medium supplemented with 5 mM β-glicerophosphate (Sigma-Aldrich, St Louis, MO, USA). Cells were transfected with both mimics and inhibitors of miR-3679-3p and miR-4274, and their corresponding controls at day 1, and day 14 after seeding. During the cell culture time period, the medium was changed once every 3 days. At day 28, cells were stained with alizarin red to quantify mineralized nodules. At this time point, media was removed from the cell monolayer and gently washed 3 times with PBS. The cells were then fixed in 10% buffered formalin for 10 minutes at room temperature. Fixative was removed and cultures were washed in PBS. The cell layer was stained with 2% Alizarin-S (Sigma-Aldrich, St Louis, MO, USA) at ~pH 4.2 for 20 minutes. Cell preparations were washed with PBS to eliminate nonspecific staining. To quantify calcium deposition, the dye was leached by the addition of 10% cetylpyridinium chloride until all of dye had been drawn. Optical density was then quantified by spectrophotometry at 550 nm, using 10% cetylpyridinium chloride as a blank reference. Statistical methods. Hardy-Weinberg equilibrium (HWE) was calculated using the chi-square test.

Cell transfection.
Multiple linear regression models were fitted to assess the association between genotyped SNPs and BMD. Log-additive, dominant and recessive models were tested for each SNP association analysis. Potential confounders considered for adjustment were densitometer devices, body mass index (BMI) and age. Association analyses were performed using R software version 2.13.2 with the SNPassoc, foreign, gdata and multtest packages.
Mann-Whitney U test was performed for OP and non-OP group comparisons in the SPSS v.12.0 for Windows. This test was also used for treatment comparisons between miRNA mimics or inhibitors and their respective controls.
Linear regression was used to analyze the association between miRNA expression levels and sample genotypes. Log-additive, dominant and recessive models were tested for each miRNA analysis. Correlation analyses were performed using R software version 2.13.2 with the SNPassoc and foreign packages.
All analyses were two-tailed, and p-values < 0.05 were considered significant.