Osteoporosis (OP) is an age-dependent metabolic bone disease characterized by a decrease in bone mineral density (BMD) and destruction of bone microstructure, which can lead to elevated bone fragility and fractures (Pisani et al., 2016). There are approximately 18.3% with OP worldwide (Salari et al., 2021). The aging of the population increases the disease prevalence, so treating OP becomes a major financial issue. The disease is primarily associated with postmenopausal women. However, approximately 25% of OP cases occur in men over 50 years of age. Its incidence is thought to be underestimated because men are under-screened for OP compared to women (Salari et al., 2021), this results in men having more complications of the disease such as vertebral or hip fractures, and they have a higher mortality rate from osteoporotic fractures compared to women (Rinonapoli et al., 2021; Wright et al., 2014). In Russia, OP is diagnosed in about 34% of women and 27% of men when a densitometric examination is performed on a random sampling (Zaigrova et al., 2017; Sozen et al., 2017; Lesnyak et al., 2018). Osteoporosis is characterized by increased bone resorption, prevailing over mineralization and anabolic processes, which leads to impaired bone strength and elevated risk of fractures, even from minor falls or impacts (Foger-Samwald et al., 2020). Genomic and multi-omics projects showed that primary osteoporosis is strongly associated with hereditary factors (Makitie et al., 2019), which determine up to 65% of the variability in bone mineral density levels (Howard et al., 1998).

At the same time, quite a few DNA loci with a high risk effect were identified (Rivadeneira et al., 2009). A meta-analysis of a genome-wide association search (GWAS) identified about 56 loci correlating with low BMD and 14 loci associated with fracture risk (Estrada et al., 2012). Morris et al. (2019) identified more than 1,000 conditionally independent signals at genome-wide significance (p < 6.6 × 10−9) mapping to 515 loci. However, replication of GEFOS/GENOMOS results on an independent sampling from the Volga-Ural region of Russia (VUR) did not confirm most of the identified DNA markers (Estrada et al., 2012). VUR is an ethnically heterogeneous region on the border of Europe and Asia, separated by the Ural Mountains. Turkic, Finno-Ugric, and Slavic peoples live here, having in their gene pool European and Mongoloid components in different proportions (Khusnutdinova et al., 1999).

Aberrations of epigenetic mechanisms are significantly correlated with elevated fracture risk and decreased BMD (Khusnutdinova et al., 1999; Xu et al., 2021). This is supported by studies of DNA methylation patterns (Visconti et al., 2021) and the microRNA pool in blood or bone tissue in patients with OP (Huai et al., 2020). Non-coding RNAs (ncRNAs) act synergistically to regulate cell proliferation, osteoclastogenesis, osteogenesis, autophagy, and other processes. MicroRNA and long ncRNA control multiple signaling pathways and gene expression. Therefore, understanding how binding site polymorphism affects the effect of RNA interference expands the possibilities for developing targeted therapies tailored to the affinity of microRNAs for a specific mRNA target when delivering RNA vectors (Lin et al., 2021). As a therapeutic precedent, gene therapy for musculoskeletal diseases has been well demonstrated in a murine ovariectomy model over the past 2 decades (Baltzer et al., 2001). The development of RNA therapies for parandontal remodeling through the regulation of microRNA activity affecting orthodontic tooth movement seems very promising (Chen and Zhang, 2023). Thus, microRNAs are promising markers for the development of diagnostic strategies and targeted therapy for primary OP (Wang and Grainger, 2012). However, there is a knowledge gap regarding the role of microRNA binding sites in forming the risk of primary OP.

Because each mRNA differs in its ability to bind to microRNAs, some do so more efficiently than others, and this is associated with sequence variations in the 3′-UTR untranslated regions of mature mRNAs (Plotnikova et al., 2019). In addition, different allelic variants have an affinity for different types of microRNAs, altering the qualitative composition and effect(s) of DNA-RNA interactions (Cloonan, 2015). Thus, polymorphism of the 3′-UTR sequence of mRNA changes the binding affinity to microRNAs and thereby may alter the regulation of target mRNA translation or induce mRNA degradation, leading to differential expression of target genes (Zhang and Wang, 2017). Previously, using the method of polygenic risk assessment, the authors developed models predicting the risk of fracture development in women from the Volga-Ural region of Russia, which, along with other loci, included polymorphic variants of microRNA binding sites (Yalaev et al., 2022). The first results in this direction were obtained by Lei et al.(2011) in a full genomic study of microRNA binding sites, where statistically significant associations of 7 polymorphic variants–rs6854081, rs1712, rs10518716, rs17054320, rs10793442, rs10098470, and rs2745426 with femoral neck BMD were revealed. The most significant associations with OP were described for rs1712 of the FBX05 gene and rs6854081 of the FGF2 gene (Lei et al., 2011). Another study found associations between rs735890 and lumbar low BMD in the elderly, probably due to changes in microRNA target sites (Amjadi-Moheb et al., 2018). Ahn et al. (2020) investigated microRNA binding site polymorphic variants in mRNA of genes involved in vitamin B metabolism, finding that combinations of risk alleles of the rs9426 (CD320), rs10418 (TCN2), rs1051296 (SLC19A1), and rs16862199 (SLC19A2) loci were associated with postmenopausal OP and compression fractures in women of Korean ancestry. Thus, studies show a significant role of polymorphic variants of microRNA binding sites in elevating the risk of osteoporotic fractures and decreasing BMD. However, studies in this direction are sporadic and need to expand the number of studied genes. In addition, replication of the previously obtained results of some authors on different cohorts of patients with primary OP from the other populations and ethnically differentiated groups is no less relevant.

The aim of this study was to search for associations of polymorphic variants of microRNA binding sites in the mRNA of genes involved in the metabolism of connective tissue in general and bone tissue in particular: rs1061947 (COL1A1), rs1031820 (COL11A1), rs9659030 (COL11A1), rs11540149 (VDR), rs6854081 (FGF2), rs1042673 (SOX9), rs10793442 (ZNF239), rs10098470 (TPD52), rs1054204 (SPARC), rs1712 (FBXO5), and rs5854 (MMP1), as well as associations of polymorphic variants rs2910164 and rs11614913 of the miR-146a and miR-196a genes with primary OP and development of clinical and genetic models. It was done to predict the risk of fractures and low BMD in general and standard localizations in men and women from the Volga-Ural region of Russia.

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.”

Table 1. Characteristics of the patient cohort by phenotypic subgroups.

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).

Table 2. Characteristics of polymorphic variants of microRNA binding sites in mRNAs.

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 I2 (the proportion of variability due to the heterogeneity of samplings) (Higgins and Thompson, 2002). When the I2 value was less than 30%, heterogeneity was assessed as mild; when I2 was between 30% and 50%, it was assessed as moderate; and when I2>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).

Associations were searched using the nonparametric χ2 criterion in the total sampling, taking into account ethnicity and sex, as well as the localization of the phenotypic features of OP: fractures and low BMD. The “typical fractures” group included fractures of the femoral neck, spine, and radial bone, and the “general fractures” group included individuals with fractures in any parts of the skeleton without regard to their localization. In addition, analysis was performed taking into account fractures or a low BMD level in individual localizations. The BMD level was categorized according to the T-criterion (low BMD level at −2.5 SD and below, up to −0.9 SD–normal mineral density level). The identified statistically significant associations are presented in Tables 5, 6, full data in supplementary data.

The Hardy-Weinberg equilibrium test revealed some deviation at the rs1712 locus (p = 0.014) in the male sampling, while equilibrium was maintained in the control group without fractures (p = 0.159). Therefore, all polymorphic loci were eligible for the association study. The minor allele frequency ranged from 8.4% at the rs1031820 locus to 46.8% at the rs1054204 locus (Table 3).

Table 3. Characterization of the studied polymorphic loci.

The authors conducted a sequential analysis of associations with OP in general–fractures in typical localizations and/or low BMD, as well as with general fractures and various localizations without regard to BMD, low BMD in general and various localizations without regard to the presence of fractures in the total sample. Then, it was conducted in cohorts of men and women separately, taking into account their ethnicity. Certain regularities were found: some associations were characterized by an increase in the level of statistical significance in the analysis of associations taking into account the fracture localization and BMD level, as well as ethnicity and sex (Tables 4, 5).

Table 4. Identified associations by pairwise comparison in the total sampling of men and women.

Table 5. Identified associations by pairwise comparison in men and women separately.

In the total sampling, the highest number of associations was found for loci rs10098470 (TPD52) and rs11540149 (VDR), which are associated with both fractures and low BMD alone and in combination with a high level of statistical significance (Table 4).

The association of T locus rs1164913 (miR-196a), A locus rs11540149 (VDR), and A locus rs10098470 (TPD52) alleles with OP was found in a pooled sampling of men and women. However, when the cohorts of men and women were analyzed separately, the identified associations persisted only in men, except for rs1164913, and there was a new association of the G allele of locus rs1054204 (SPARC) with OP in general in men. When the sample was further stratified into subgroups based on the presence of fractures and low BMD levels separately, the rs1164913 locus was associated with general fractures and the rs1054204 locus was associated with low femoral neck BMD in men (Table 5). This indicates that primary OP is a genetically heterogeneous disease with different contributions of risky marker variants according to the sex.

This is also true for ethnicity, when cohorts are more homogeneous due to belonging to the same ethnic group, the association becomes more significant. For example, the A allele of the rs11540149 locus (VDR) showed an association with fractures in the pooled sampling of men and women overall (OR = 5.2), in Russians (OR = 7.24), and with typical fractures–OR = 6.85 (overall) and OR = 8.20 (in Russians).

When dividing the sampling by sex, we observed the persistence of the identified associations mainly in men despite the fact that their effect was lower than in women (Table 4). For example, the OR for the association of the C allele of the rs1712 locus (FBX05) with fractures overall in men reached 2.32 compared with 1.4 in the pooled sampling of men and women. However, no statistically significant associations were found with general fractures in women. The C allele of the rs1712 locus (FBX05) was associated with general fractures, as well as typical osteoporotic fractures and radial fractures, with OR increasing from 1.34 to 1.74, respectively.

Figure 1 presents a Venn diagram showing the contribution of the studied polymorphic variants of microRNA binding sites in target genes to fracture risk and low BMD individually and common markers for the two endophenotypes of primary OP by ethnicity and sex.

Figure 1. Genes involved in the formation of osteoporotic fractures and low BMD in men and women from the Volga-Ural region of Russia, taking into account ethnic characteristics. For easier perception, each phenotype localisation and ethnic specificity characteristics are highlighted in different colour.

Polymorphic variants of the VDR and TPD52 genes are common risk markers of fractures and low BMD for men and women, FBX05 and COL11A1 – for men. At the same time, the SOX9, MMP1, ZNF239, COL1A1, FGF2, miR-146a, and miR-196a gene loci are characterized by associations for separate groups depending on sex and ethnicity, as well as the localization of the pathological process (Figure 1).

In the total sampling of men and women, the association of the polymorphic variant rs10098470 of the TPD52 gene with lumbar fractures in Russians (χ2 = 23.600, OR = 7.77, p = 1.183e−06) reaches the highest level of statistical significance (Table 3), which is also associated with a low level of BMD in general, as well as femoral neck in the pooled sampling, but with a lower level of significance. When the sampling was divided by sex, the associations with lumbar fractures remained in cohorts of men and women of Russian ethnicity.

In the sampling of women, variant A of the SOX9 gene locus rs1042673, associated with lumbar fractures in women of Tatar ethnicity, has the highest odds ratio (OR = 7.00). At the same time, the G allele is associated with low BMD in women of Russian origin with a significantly lower odds ratio (OR = 1.61). The risk allele A of the locus rs10793442 of the ZNF239 gene, which is probably a specific marker of radial fractures in women of Tatar ethnicity, reaches statistical significance.

In women, allele G of locus rs6854081 of the FGF2 gene is associated with femoral neck fractures, with OR increasing from 4.167 to 5.399 in women of Russian ethnicity. Probably, there are population differences in the distribution of allele frequencies, which confirms the relevance of the search for ethnospecific markers of primary OP.

The loci associated with OP traits in men and women separately contribute to the notion of sex differences in the mechanisms of fracture risk and low BMD. The results require validation in independent samplings.

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 (I2>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 (I2 = 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.

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).

Figure 2. Prognostic models for assessing the risk of osteoporosis in general (A) and osteoporotic fractures (B) in women.

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).

Figure 3. Prognostic models for osteoporosis risk assessment in general (A) and osteoporotic fractures (B) in men.

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).

Figure 4. Predictive models for assessing the risk of osteoporosis in general (A) and of osteoporotic fractures (B) in a pooled sampling of men and women.

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.

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).

Table 6. Predicted microRNAs, their binding free energy, and binding site localization depending on the polymorphic variant rs6854081 in the mRNA transcript of FGF2 gene 002,006.6 at alternative C allele of FGF2.

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 7. Predicted microRNAs, their binding free energy, and binding site localization depending on the polymorphic variant rs10098470 in the mRNA transcript of the TPD52 gene NM_001287144.2 at alternative allele C of the TPD52 gene.

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.

Table 8. Predicted microRNAs, their binding free energy, and binding site localization depending on the rs1712 polymorphic variant in FBXO5 gene mRNA transcript NM_012177.5 at alternative C allele of FBXO5.

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).

Table 9. Predicted microRNAs, their binding free energy, and binding site localization as a function of the rs1712 polymorphic variant in the mRNA transcript of the FBXO5 NM_001142522.3 gene at alternative C allele of FBXO5.

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).

Table 10. Predicted microRNAs, their binding free energy, and binding site localization as a function of the rs1712 polymorphic variant in the mRNA transcript of the FBXO5 NM_001142522.3 gene at the alternative C allele of FBXO5.

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.

Each microRNA has certain chemical affinity, which is affinity for the target mRNA in the regulation of which it is involved. The efficiency of their interaction is determined by different thermodynamic and chemical properties of these macromolecules. The affinity of mRNA and microRNA affects the regulation of gene expression–the higher it is, the stronger is the effect of RNA-induced gene silencing (RNA interference), i.e., suppression of gene expression (up to complete degradation of mRNA with cessation of translated protein expression) (Brien et al., 2018).

Several factors are known to influence the affinity of mRNAs and microRNAs:

1. Complementarity of the canonical interaction between the “seed”-region of the microRNA (2–8 nucleotides at the 5′-end) and the 3′-untranslated region (UTR) in the mRNA binding site of the target gene. The higher is this complementarity, the more effective is the RISC complex of RNA interference (or discrete microRNA without the RISC complex) (Bartel, 2009).

2. Non-canonical interaction of the “seed”-region of microRNA with the coding (CDS) and 5′-UTR regions of mRNA. According to CLASH data, most of the interactions (about 58%) lie in the CDS region, 38% lie in the 3′-UTR, and 4% are in the 5′-UTR in mRNA. However, non-canonical interactions are not as well studied as in the case of microRNA interactions with the 3′-UTR of mRNAs. Besides, their functional significance in mRNA upregulation was not established. Therefore, their role in influencing the regulation of gene expression needs additional functional studies (Helwak et al., 2013).

3. Spatial configuration of the secondary structure of mRNA (the double helix of RNA is formed by two linear molecules connected to each other along the entire length by hydrogen bonds). It is known that complex conformational changes in the spatial orientation of mRNA molecules affect the efficiency of interaction between mRNA and microRNAs that have affinity for binding sites in different regions of mRNA (Yang et al., 2020).

Moreover, SNPs of binding sites can affect the secondary structure of mRNA, potentially limiting the accessibility of the microRNA binding site to the binding centers of the RNA-induced gene silencing complex, thereby, affecting both the composition of microRNAs and their spatial accessibility for binding (Rykova et al., 2022).

The field of research on the influence of polymorphisms in 3′UTR-regions of mRNA on the risk of various diseases is actively developing. However, due to the huge number of microRNA variations in humans (more than 2,500), as well as the huge representation of microRNA binding sites in the genome, and taking into account the presence of multiple mRNA transcripts for each gene, an unprecedented number of possible and complex variants of functional consequences of microRNA-mRNA interactions appears. To reveal the mechanisms and significance of these interactions, there are various bioinformatic tools that aim to establish the role of aberrations of these interactions in the risk of disease development, including osteoporosis. Therefore, an urgent task is to search for the most functionally significant changes in predicted mRNA-microRNA interactions and their role in disease pathogenesis depending on polymorphic variants of microRNA binding sites, for which a statistically significant association with pathology phenotypes were shown in the studied patient cohorts (Chhichholiya et al., 2021).

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.