The case-control study included 701 postmenopausal women (mean age = 61.95 ± 7.94) and 501 men (mean age 62 ± 10.8) who underwent medical examination between 2004 and 2011 at the City Clinical Hospitals No. 5, No. 21, and No. 22 in Ufa and Regional Clinical Hospital No. 1. in Yekaterinburg. Ethnic composition of women: Russians–516, Tatars–185. Ethnic composition of men: Russians–470, Tatars–31. The sampling of patients consisted of people with primary OP; the control group included people without fractures and with a normal level of BMD. In the group of women, the number of patients with fractures in general is 280, with fractures in typical localization–160, with low BMD–324. In the group of men, the number of patients with fractures in general is 145, with fractures in typical localization–83, with low BMD–304. The subgroups and their detailed characteristics are summarized in Table 1. The authors completely excluded the presence of any family groups and relatives on the basis of questionnaires and family history data. Exclusion criteria included a history of alcohol and drug abuse, long-term use of glucocorticoids, hormone replacement therapy, current treatment of acute diseases, and the presence of chronic diseases that affect bone metabolism. The level of BMD was measured by Dual-energy X-ray absorptiometry (DEXA) using a Hologic QDR 4500/A DXA system (United States) in standard localizations (femoral neck and lumbar spine) [DEXA system in Ufa and Yekaterinburg demonstrated a high level of reproducibility (coefficient of variation CV<2.6%) and a low level of error (ε<1.5%)]. The general sampling was divided according to the T-criterion–a T value between +2.5 and −0.9 standard deviations (SD) characterized normal BMD, values between −1.0 and −2.5 SD–osteopenia, and values less than −2.5–OP (according to the World Health Organization’s recommendations). The presence of osteoporotic fractures in standard localizations (axial part of the femur, lumbar spine) in general and separately, as well as in combination with any other skeletal fractures, was taken into account. Each participant signed an informed consent form for participation in the study in accordance with the standards of the World Medical Association Declaration of Helsinki “Ethical Principles for Scientific Medical Research Involving Human Subjects.”

DNA was isolated using the phenol-chloroform extraction method from peripheral blood leukocytes according to the protocol of Mathew (1985). The quality of isolated DNA was checked using a NanoDrop 1,000 spectrophotometer (Thermo Scientific, United States). DNA concentration was measured using a Qubit 4 fluorimeter (Thermo Scientific, United States). Genotyping of the studied samples was performed using KASP™ (Kompetitive Allele Specific PCR) technology in “real-time” by the endpoint (Semagn et al., 2013) on the QuantStudio 12K Flex Real-Time platform. Target microRNA loci were selected using the microRNA database (http://compbio.uthsc.edu/miRSNP, access: 1.02.2024) (Table 2).

The PLINK software (v. 1.09) was used to search for associations of alleles and genotypes with case-control fractures and a low BMD level using Pearson’s concordance criterion. The degree of association was assessed in odds ratio (OR) values using the formula: OR = (a*d)/(b*c), where a is the frequency of the trait in the sampling of patients, b is the frequency of the trait in the control sampling, c is the sum of the frequencies of the other traits in the sampling of patients, and d is the sum of the frequencies of the other traits in the control sampling. Tests were performed for the two-sided significance level; differences at p < 0.05 were considered statistically significant.

The PLINK program (v. 1.09) was used for meta-analysis of the results. To calculate the mean OR and significance level, fixed effects models (Mantel-Haenszel method) and random effects models (Dersimonian-Laird method) were considered. To assess the statistical heterogeneity of different samples, Cochran’s Q test was used; differences at p < 0.1 were considered statistically significant. The level of heterogeneity was determined using the statistical criterion I² (the proportion of variability due to the heterogeneity of samplings) (Higgins and Thompson, 2002). When the I² value was less than 30%, heterogeneity was assessed as mild; when I² was between 30% and 50%, it was assessed as moderate; and when I²>50%, it was assessed as highly heterogeneous.

Logistic regression analysis was performed using the MedCalc software (v. 22.016). The disease phenotype (presence/absence of fracture, low BMD/normal BMD, combined conditions) served as the dependent variable, and risk alleles of polymorphic variants and clinical features (BMI, BMD level) served as independent variables (predictors). In the course of logistic regression analysis, several dozens of regression equations were obtained, which the equations with the highest values of statistical significance indices were selected from. To determine the quality of the obtained prediction model, ROC (Receiver Operating Characteristic) analysis was used, the sensitivity and specificity of the model were determined, as well as the AUC (Area Under Curve)–an indicator of the area under the ROC curve (Wray et al., 2010).

To evaluate the functional significance of polymorphic variants associated with the risky phenotype, the authors in silico predicted microRNAs that potentially interact with mRNAs of genes, wherein the studied high-risk loci are localized. For this purpose, we used miRBase, NCBI, and IntaRNA 2.0 bioinformatics tool for prediction of mRNA-microRNA interactions to analyze changes in the qualitative composition of microRNAs depending on allelic variants of SNP that reached a high level of significance within the association analysis. Three files were prepared for this purpose:

1. Original reference sequences of gene mRNA transcripts in FASTA format from NCBI repositories.

Further, IntaRNA computational algorithms and these files were used to predict all microRNAs that had affinity for mRNA transcripts separately for the reference sequence of known transcripts and the sequences of transcripts with an alternative risk allele. After that, we excluded all those microRNAs which minimum free energy of hybridization was >-8 from all microRNAs that were predicted for these transcripts. After filtering, we matched and compared the lists of these microRNAs for the original mRNA transcript sequence and the new transcript sequence with the alternative allele.

Due to ethnic heterogeneity and the presence of sexual dimorphism in OP, the authors assessed the genetic heterogeneity of the studied sampling using meta-analysis tools based on the PLINK software. Meta-analysis provides the optimal opportunity to find effects that are specific to the whole cohort. Heterogeneity of information (inherited genetic background, covariation of genotypes and phenotypes) affects not only statistical power but also increases the potential for false positives (Brien et al., 2018). Fixed-effects models assume that there is a common true genetic effect in all association studies, and any variation is explained by random error; on the other hand, random-effects models assume that there are different effect sizes in the association tests, and any differences, i.e., heterogeneity, are due to real population differences. This allows cohort information to be harmonized as much as possible. A pooled sampling of men and women of Russian and Tatar ethnicity was assessed, and then, men and women were assessed separately.

For the rs10098470 (TPD52) and rs11540149 (VDR) loci, the authors found high heterogeneity in the pooled sampling of men and women (I²>50%) and applied randomized p values (p(R)) for the random-effects model, which did not show the significance of these loci for low-BMD and fracture risk in the pooled sampling, general fractures, without regard to the BMD level. For the rs10098470 (TPD52) locus, heterogeneity is apparently due to the fact that the A allele is a marker of fractures and low BMD in combination and separately in Russians in general and for men and women separately, and the associations with fractures and low BMD of specific skeletal sites did not always coincide. The same patterns are characteristic of the A allele of the rs11540149 locus (VDR), which is associated with general fractures and standard localizations mainly in men, except for lumbar fractures and a low level of lumbar BMD in women, all associations were stronger in representatives of Russian ethnicity.

In the low-BMD sample, meta-analysis revealed a high level of heterogeneity for loci rs1031820, rs10793442, and rs1042673 (I² = 73.04, 54.24, and 63.08, respectively). The G allele rs1031820 of the COL11A1 gene was associated with lumbar fractures in the pooled sampling of men and men of Russian ethnicity, and the A allele was associated with low femoral neck BMD in the pooled group of men and women of Tatar origin due to differences in the frequency of the minor allele A, which was 3.5 times higher in men than in the sampling of women (0.301 and 0.084, respectively). The high level of heterogeneity at the locus rs10793442 of the ZNF239 gene can probably be explained by the association with BMD only in Russians. For the SOX9 gene rs1042673 locus, heterogeneity is observed in a sampling of women; this locus is associated with a low level of BMD in women of Russian origin. However, the risky allele in this case is allele A, whereas allele G is associated with lumbar fractures in women of Tatar ethnicity. It is necessary to expand the sample to confirm the obtained associations, as well as to replicate them in an independent sampling.

Thus, heterogeneity in ethnicity, sex, and fracture-associated and low BMD of individual localizations for loci rs10098470 (TPD52) and rs11540149 (VDR) in the studied cohorts of women and men from the Volga-Ural region was revealed, where the majority of the sampling consisted of Russians, belonging to the Slavic group, and Tatars, a Turkic group of the Altai language family, having diversity of ethnogenesis, diet and religion, but living on a common territory for a long time.

The conducted studies yielded molecular genetic risk markers for the development of OP in general, as well as its various phenotypes. However, the analysis of individual risk factors does not take into account their interaction, which is necessary for predicting the risk of developing the studied pathology. It is also of interest to evaluate the joint influence of genetic and clinical factors, such as body mass index (BMI) and BMD levels of the lumbar spine and femoral neck. In order to study the mutual influence of clinical and molecular genetic factors on the risk of OP and to create diagnostic algorithms, the authors performed statistical processing of the obtained results by logistic regression and ROC analysis methods with the calculation of the area under the curve (AUC), the values of which range from 0.6 to 0.7 indicating the average diagnostic value of the model; 0.7–0.8–good; 0.8–0.9–very good.

Predictive models were developed to assess the risk of OP in general, as well as osteoporotic fractures in women and men individually, and in a pooled sampling. The model for assessing the risk of OP in women included a low BMD in standard localizations, as well as loci rs6854081 (FGF2) and rs1054204 (SPARC); the model for predicting osteoporotic fractures included a low BMD and locus rs1712 of the FBXO5 gene (Figure 2). The models were characterized by the highest level of predictive significance (AUC = 0.909 and AUC = 0.830, respectively). It should be noted that the exclusion of genetic loci from the analysis reduced the predictive value of the models (AUC = 0.844 and AUC = 0.805, respectively).

The model for OP risk assessment in men included a decreased level of BMD in standard localizations and the locus rs1712 (FBXO5) and rs11540149 (VDR), while the model for fractures included a low level of BMD and the loci rs1712 (FBXO5), rs1042673 (SOX9), and rs10793442 of the ZNF239 gene (Figure 3). The models were characterized by a sufficient level of predictive significance (AUC = 0.764 and AUC = 0.716, respectively), there was a tendency of decreasing predictive value when genetic predictors were excluded (AUC = 0.740 and AUC = 0.649, respectively).

The OP risk assessment model in the pooled sampling of men and women also included decreased BMD levels in standard localizations, as well as loci rs1712 (FBXO5), rs6854081 (FGF2), rs1054204 (SPARC), rs10793442 (ZNF239), and the model for fractures–decreased BMD levels and loci rs1712 (FBXO5), rs10793442 (ZNF239), and rs11614913 of the miR-196a gene (Figure 4). The models were characterized by a sufficient level of predictive significance (AUC = 0.838 and AUC = 0.771, respectively), the trend of decreasing predictive value when genetic predictors were excluded was maintained (AUC = 0.808 and AUC = 0.720, respectively).

Thus, the conducted logistic regression and ROC analyses made it possible to identify several models with high prognostic significance for predicting the risk of primary OP and osteoporotic fractures in men and women in general and separately, which confirmed the role of microRNA binding site loci in the FBXO5, FGF2, SPARC, ZNF239, and miR-196a genes in the formation of the disease phenotype obtained by association analysis as well. The loci rs10098470 (TPD52) and rs11540149 (VDR), for which high heterogeneity was detected, were not included in the models with high statistical significance.

Based on the association study data of microRNA binding site polymorphic variants rs6854081, rs10098470 and rs1712, which showed the highest association signal with osteoporosis phenotypes, we performed a comparative analysis of predicted microRNAs with potential affinity for mRNA of FGF2, TPD52 and FBXO5 genes. It was done to test how binding site variation at the studied loci affected affinity for microRNAs, whether the qualitative composition of predicted microRNAs changed depending on the presence of the alternative allele, and whether the localization of microRNA binding to mRNA changed in the 3′UTR region of mRNA.

We found that in silico, at risk allele C, a novel microRNA hsa-miR-6838-3p was predicted in transcript 002,006.6, which had no affinity for the reference mRNA sequence. Through analysis of the entire mRNA transcript sequence, we found that a new binding site in the 3′UTR, which previously had no affinity for the reference sequence, was predicted for the microRNAs hsa-miR-34a-5p, hsa-miR-380-3p, hsa-miR-654-3p, hsa-miR-3649, and hsa-miR-6770-3p (Table 6).

It was revealed that at risk allele C in the 3′UTR of TPD52 gene in transcript NM_001287144.2, three microRNAs lost affinity to the reference sequence: hsa-miR-548d-3p, hsa-miR-873-3p, hsa-miR-1537-3p. The microRNA hsa-miR-6753-3p is also predicted to have affinity for the 3′UTR mRNA with high binding energy (−13.49) (Table 7). Analysis of predicted transcript interactions with NM_005079.4 microRNAs showed that these microRNAs retained affinity for mRNA, but the binding site changed to a position upstream of the 3′UTR region. Analysis of transcript NM_001025253.3 in the 3′UTR predicts the appearance of hsa-miR-6753-3p microRNA affinity.

Table 8 shows the results of comparison of all predicted microRNAs, which were filtered based on the appearance or loss of affinity for mRNA or on the change of the binding site in the 3′UTR of the FBXO5 gene. Comparative analysis of the predicted microRNAs for transcript NM_012177.5 revealed several changes. At risk allele A, affinity for the new microRNAs hsa-miR-105-5p, hsa-miR-379-5p, hsa-miR-552-3p, and hsa-miR-617 appeared, and affinity for hsa-miR-187-5p, hsa-miR-196b-3p, and hsa-miR-1273c was completely lost (Table 8). In addition, a change in the microRNA binding site with localization in the 3′UTR with high binding energy was detected for microRNAs hsa-miR-139-5p, hsa-miR-3125, hsa-miR-122b-5p and hsa-miR-5579-3p.

Comparative analysis of predicted microRNAs for transcript NM_001142522.3 revealed that the risk allele resulted in a new affinity for the microRNAs hsa-miR-105-5p, hsa-miR-379-5p, hsa-miR-552-3p, hsa-miR-617, and hsa-miR-3085-5p, and loss of affinity for hsa-miR-187-5p, hsa-miR-382-3p, hsa-miR-196b-3p, hsa-miR-549a-3p, hsa-miR-2277-3p, and hsa-miR-1273c (Table 9).

In addition, changes in the localization of binding sites were noted. Thus, for microRNAs hsa-miR-139-5p, hsa-miR-6878-5p, hsa-miR-3125, hsa-miR-3127-5p, and hsa-miR-3945, the binding site moved to the 3′UTR-region of the gene with high binding energy. At the same time, the microRNA binding sites hsa-miR-146b-5p, hsa-miR-181d-3p, hsa-miR-498-3p, hsa-miR-371b-5p, and hsa-miR-5194 binding site moved upstream beyond the 3′UTR (Table 10).

The point is that an alternative allele of the binding site can lead to a change in the composition of microRNAs with potential affinity for mRNA. In turn, a change in the composition of microRNAs can lead to a functional change in protein production activity and play a biological role in the dysregulation of the activity of genes involved in bone tissue metabolism and homeostasis.

Fibroblast growth factor type 2 (FGF2) is a highly pleiotropic member of a large family of growth factors with a broad spectrum of activity, including proliferation and self-renewal of human pluripotent stem cells. The FGF2 gene does not have alternative splicing (Nickle et al., 2023). Instead, isoforms of the FGF2 protein are expressed from a single mRNA: a high molecular weight 34 kDa HMW isoform localized to the cell nucleus and a low molecular weight 18 kDa LMW isoform predominantly localized to the membrane as a ligand for fibroblast growth factor receptor (FGFR). HMW has an inhibitory effect on mineralization and LMW promotes bone formation (Coffin et al., 2018).

Since the role of aberrant HMW expression in low bone mineralization was proven, we focused on the transcript encoding the HMW isoform of FGF2. Using IntaRNA algorithms (Mann et al., 2017), we evaluated predictions of microRNA and NM_002006.6 (HMW mRNA) transcript interactions and found that the risk allele of the rs6854081 polymorphic variant in the 3′UTR in silico results in affinity for FGF2 gene mRNA with a novel microRNA–hsa-miR-6838-3p.

Hsa-miR-6838-3p is overexpressed in patients with osteosarcoma (Wang et al., 2022). In addition, decreased levels of miR-6838-5p increase the level of FGF2 gene expression (Zhang et al., 2023). This was found in a study by Zhang et al. (2023), who found that Circ_0001667 RNA was associated with breast cancer and angiogenesis processes involving miR-6838-5p. Hsa-miR-34a-5p with affinity for in the 3′UTR region of the mRNA is associated with the induction of osteogenic differentiation of bone marrow mesenchymal stem cells. This microRNA modulates bone metabolism by targeting HDAC1 and stimulating ER-α transcription (Sun et al., 2023). The interaction of miR-6838-3p and miR-34a-5p with FGF2 gene mRNA may be associated with decreased FGF2 expression; however, additional studies are needed to test this hypothesis.

In the study by Lei et al. (2011), rs6854081 (FGF2) is a functionally investigated locus. They not only found an association of this variant with low BMD, but also identified aberrant expression of the FGF2 gene in women with osteoporosis, with the rs6854081 locus being associated with altered osteoclast differentiation (Lei et al., 2011). We reproduced a statistically significant association of this locus with OP: in our sample, the G allele of rs6854081 is a predictor in a prognostic model of the risk of low BMD formation and is also associated with femoral neck fractures in men and women. Studies show that the FGF2 gene transcript with the major T allele of the rs6854081 locus has a lower expression level of monocytes and is associated with a high level of BMD, which only confirms the potentially risky role of the C allele.

D52 is one of the members of the tumor protein family (Abe et al., 2021). Several studies identified the TPD52 gene sequence by its increased expression in human breast cancer tissue and in lung cancer cell lines. Several studies identified overexpression of the TPD52 gene in osteosarcoma cell lines (Abe et al., 2021). About 80% of the TPD52 transcript was initially identified as a 3′-untranslated region, so a large number of microRNAs involved in the regulation of TPD52 gene expression was identified (Han et al., 2015). This gene has 29 transcripts (splice variants), 208 orthologs, and 3 paralogs. Eight of them encode different protein isoforms. The polymorphic variant rs10098470 is localized in the TPD52 gene encoding tumor protein D52. The function of this protein in the pathogenesis of osteoporosis was not studied. Nevertheless, the variant rs10098470 was still identified in the GWAS study with a high level of significance (Lei et al., 2011). Studies showed that TPD52 family proteins were not involved in osteoclast differentiation. Still, they indicated that the protein played an important role in terminal chondrocyte differentiation during endochondral ossification, maintaining the pre-hypertrophic state of chondrocytes (Ito et al., 2017). Hence, given the replicated association in two cohorts [ours and Lei et al. (2011)], and considering the important role of TPD52 protein in ossification, the rs10098470 variant deserves further functional studies. Using IntaRNA algorithms, we evaluated predictions of microRNA and NM_002006.6 transcript (HMW mRNA) interactions and found that the risk allele of the rs6854081 polymorphic variant in the 3′UTR in silico results in affinity for FGF2 gene mRNA with a novel microRNA–hsa-miR-6838-3p.

FBXO5 is an F-box 5 protein that constitutes one of the four subunits of the ubiquitin protein-ligase complex called SCF (SKP1-cullin-F-box), functioning in a phosphorylation-dependent ubiquitination process. FBXO5 gene expression is upregulated in many human cancers and is associated with chromosomal instability. A marked increase in mineralization was observed in groups overexpressing both transcripts compared to controls. Knockdown of FBXO5 resulted in decreased migration and osteogenic differentiation functions. In contrast, overexpression of FBXO5 enhanced the potential for migration and osteogenic differentiation. These data suggest that FBXO5 plays a positive role in the regulation of migration and osteogenic differentiation (Liu et al., 2018). In the work of Lei et al. (2011), the most significant level of association with the level of femoral neck BMD in women was revealed for the rs1712 variant of the FBX05 gene (p = 2.5e⁻³). We replicated this result. However, the difference is that in our work, a significant association was found only in men, who had the C allele of the rs1712 locus (FBX05) associated with fractures of all localizations except the femoral neck, as well as with a low level of BMD in general. In the PolymiRTS database, microRNA miR-549a has affinity for this site, which level is elevated in Huntington’s chorea and decreased in fibroblasts in patients with nonsyndromal cleft palate.

At present, the role of the FBX05 gene in the pathogenesis of osteoporosis remains largely unknown, but the reproducibility and high significance in OP risk indicate the high prognostic potential of the rs1712 locus. We assessed whether the affinity for microRNAs of transcripts of this gene changes depending on the rs1712 in silico variant and found that the risk allele was associated with a significant change in microRNAs potentially interacting with the gene’s mRNA. Among the microRNAs for which a new affinity appeared, miR-617 is of interest, as this microRNA is involved in the regulation of osteogenic cell differentiation (Ge et al., 2020). Among the microRNAs that lost affinity, the most interesting is miR-382, which stimulates osteogenic differentiation in MG63 cells (Su et al., 2023). Further functional studies are needed.

Thus, within the framework of the study, for the first time, new regularities were revealed in the little-studied area of epigenetic factors of primary OP–the role of polymorphic variants of microRNA binding sites in the risk of fracture development and the formation of low BMD. The study of the influence of polymorphic variants of microRNA binding sites on the risk of OP is important because the discovered patterns suggest a high role of polymiRTS in increasing the risk of osteoporotic fractures, and therefore can be of potentially great value for the development of methods of early diagnosis of OP to prevent fractures.