Association between miRNA Target Sites and Incidence of Primary Osteoarthritis in Women from Volga-Ural Region of Russia: A Case-Control Study

Over the past decades, numerous studies on the genetic markers of osteoarthritis (OA) have been conducted. MiRNA targets sites are a promising new area of research. In this study, we analyzed the polymorphic variants in 3′ UTR regions of COL1A1, COL11A1, ADAMTS5, MMP1, MMP13, SOX9, GDF5, FGF2, FGFR1, and FGFRL1 genes to examine the association between miRNA target site alteration and the incidence of OA in women from the Volga-Ural region of Russia using competitive allele-specific PCR. The T allele of the rs9659030 was associated with generalized OA (OR = 2.0), whereas the C allele of the rs229069 was associated with total OA (OR = 1.43). The T allele of the rs13317 was associated with the total OA (OR = 1.67). After Benjamini-Hochberg correction, only rs13317 remained statistically significant. According to ethnic heterogeneity, associations between the T allele (rs1061237) with OA in women of Russian descent (OR = 1.77), the G allele (rs6854081) in women of Tatar descent (OR = 4.78), the C allele (rs229069) and the T allele (rs73611720) in women of mixed descent and other ethnic groups (OR = 2.25 and OR = 3.02, respectively) were identified. All associations remained statistically significant after Benjamini-Hochberg correction. Together, this study identified miRNA target sites as a genetic marker for the development of OA in various ethnic groups.


Introduction
Osteoarthritis (OA) is one of the most common joint diseases, often causing persistent pain and immobilization [1,2]. The genetic contribution to OA varies from 40% to 65%, depending on the presence of the pathology in the immediate family, the location of the affected joint, sex, and ethnicity of the patients [3][4][5][6]. A large number of studies using both the gene-candidate approach and the genome-wide association study (GWAS) have been conducted, but the existing findings still fail to elucidate the molecular pathogenesis of the disease [7]. The study of miRNA target sites is a promising new direction in molecular diagnostics. MiRNAs are single-stranded non-coding RNAs that are about 19-25 bps long and inhibit protein expression by directly binding to the 3 UTR region of target mRNAs. Recent studies have shown that miRNAs are also involved in the development and progression of OA [8]. A large-scale analysis of the expression of 365 miRNAs using microarrays revealed nine miRNAs (miR-16, miR-22, miR-23b, miR-30b, miR-103, miR-223, miR-377, miR-483, and miR-509) that were associated with upregulation of the genes potentially involved in cartilage homeostasis in patients with OA compared to normal cartilage, and seven miRNAs (miR-25, miR-26a, miR-29a, miR-140, miR-210, miR-337, and miR-373) that were associated with the downregulation of the PPARA, BMP7, and MMP13 genes [9]. Recent studies of healthy cartilages and cartilages with knee OA have shown that miR-125b suppresses aggrecanase gene expression (ADAMTS4) and IL-1β activation in cartilage with OA [10]. Park et al. demonstrated that miR-558 reduces the cytochrome-C oxidase-2 (COX2) gene expression through the activation of IL-1β and stimulates the catabolism of chondrocytes in OA [11]. Even small variations in the 3 UTR of the target sequence can alter the miRNA binding, leading to changes in these key regulatory interactions [12], so that even a single nucleotide difference can have profound effects [13]. Aberrations of miRNA expression patterns have been shown to be involved in oncogenesis, systemic diseases of connective tissue, osteoporosis, and a number of other conditions [14][15][16]. As more and more miRNAs are revealed to play a role in chondrogenesis, cartilage homeostasis, and OA pathogenesis, it is imperative to identify the key miRNAs and their targets that could potentially affect the disease progression. The aim of our work was to study the microRNA binding sites in patients with OA of different localizations in various ethnic groups of the Volga-Ural region of Russia.

Materials and Methods
Patient samples. A total of 417 women (mean age 51.67 ± 11.5) from inpatient care units of Ufa (Republic of Bashkortostan, Russia) hospitals participated in the study between January 2013 and August 2017. The patients were examined to diagnose for osteoarthritis according to the criteria of the American College of Rheumatologists (1995) with X-ray confirmation. Overall, 356 women with OA were recruited and divided into three groups as follows: (1) subgroup 1 included 84 women affected by generalized OA, (2) subgroup 2 included 197 women affected by knee OA, and (3) subgroup 3 included 75 women affected by hip OA. Patients with oncological pathology, systemic diseases of the connective tissue (rheumatoid arthritis, systemic lupus erythematosus, etc.), signs of an active inflammatory process of both infectious and non-infectious etiology, or traumatic joint injuries in the anamnesis were excluded from the study. A total of 161 normal healthy individuals were recruited to form the age-matched control group. The protocol was approved by the Ethics Review Committee of the Bashkir State Medical University, and signed informed consent was obtained from all of the participants. The ethnic composition of the sample was as follows: 144 Russian (Slavic group of the Indo-European language family), 159 Tatar (Turkic branch of the Altaic language family), and 114 mixed and representatives of small ethnic groups (Table 1). Patient and public involvement participants were not involved in designing the research question, conducting the study, or in the interpretation or writing of the results. There are no plans to involve participants or relevant patient communities in dissemination of results. Results are not disseminated to study participants.
DNA extraction and genotyping. Genomic DNA was extracted from peripheral blood leucocytes by phenol-chloroform extraction (Mathew, 1984). Genotyping was performed using a competitive allele-specific PCR (KASP)-patented technology by LGC-Genomics.
Selection of the miRNA target loci was carried out using the database of the National Center for Biotechnological Information (https://www.ncbi.nlm.nih.gov/, accessed on 1 November 2020), Ensemble GenomeBrowser (www.ensembl.org, accessed on 1 November 2020) and the base of polymorphisms of microRNA target sites (http://compbio.uthsc. edu/miRSNP, accessed on 1 November 2020). The candidate gene selection was based on the involvement in connective tissue metabolism or disease development ( Table 2).
hsa-miR-3919 hsa-miR-3925-5p hsa-miR-7153-5p Studies of their polymorphic variants in the coding part do not fully explain the development of pathology [17][18][19][20][21]. The genotyping results were obtained using QuantStudio 12K FlexReal-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). Genotyp-ing was performed by the laboratory staff who were "blind" to the status of the samples. Further, 10% of the samples were tested twice to validate the genotyping results with 100% reproducibility. Two authors independently reviewed the genotyping results, data entry, and statistical analyses.
Statistical analysis. A standard χ 2 -analysis was used to examine the differences in allelic frequencies and genotype distributions between the patients with OA and controls using the Statistica software version 6.0 (http://statsoft.ru/, accessed on 1 November 2020). Considering the type I error caused by multiple testing, p values were adjusted by calculating the FDR value using the Benjamini-Hochberg method (https://tools.carbocation.com/ FDR, accessed on 1 November 2020). Hardy-Weinberg equilibrium and the heterozygosity of the studied loci were estimated using the Haploview 2.0. The OR and 95% CI were calculated using the reported risk allele as a reference. Statistical significance was considered at p < 0.05.

Genetic Association Analysis
A total of 12 polymorphic loci in the 3 UTR region of 17 genes, involved in the formation and degradation of connective tissue, were genotyped. Hardy-Weinberg equilibrium was maintained for all the loci except for rs9978597. When considering the control sample independently, the indicator corresponded to equilibrium values ( Table 2). The distributions of alleles and genotypes for the selected loci are available from the Supplementary File (Tables S1-S14).

Ethnic Analysis
Due to the profound ethnic heterogeneity of residents of the Volga-Ural region of Russia [22], we analyzed the distribution of alleles and genotypes in Russian, Tatar, and mixed individuals combined with small ethnic groups, which included Bashkir, Chuvash, Kabardian, Mari, Armenian, Azerbaijani, Belarusian, and Ossetian. The following associations between a given allele and incidence of total OA in women of a certain descent were identified: T allele (rs1061237) was enriched in women of Russian descent (OR = 1.77; 95% CI 1.07 to 2.94), the G allele (rs6854081) was enriched in Tatar women (OR = 4.78; 95% CI 1.89 to 12.02), the C allele (rs229069), and T allele (rs73611720) prevailed in the groups of mixed ancestry and small ethnic groups combined (OR = 2.25; 95% CI 1.30 to 3.89 and OR = 3.02; 95% CI 1.38 to 6.60, respectively) ( Table 4). All associations remained statistically significant after Benjamini-Hochberg correction.  Subsequently, we performed meta-analysis to verify the collected data and identified general and ethnospecific trends (Table 5).
In rs6854081 and rs4836732 loci, high genetic heterogeneity of the studied sample was found (I 2 = 79.81 and 77.43, respectively). These loci did not produce statistically significant values in regard to the incidence of the total OA in women but were identified as ethnically specific markers of OA progression in women of Tatar descent. High heterogeneity was also inherent in rs7639618, rs1061237, and rs9350591 loci, which were distinguished as the markers of the OA progression in women of Russian descent. However, these associations did not remain statistically significant after adjusting for multiple comparisons. The association of the rs229069 is noteworthy. It was associated with total OA and, according to the high ethnic specificity and genetic heterogeneity, with total OA only in mixed and other ethnic groups at the same time. Most likely, the rs229069 is an ethnospecific marker in Bashkir women, since the number of women of Bashkir descent prevailed in this cohort. The obtained associations must be confirmed in future studies with a larger sample size and a cohort restricted to Bashkir women exclusively.

Discussion
It is generally recognized that OA is a complex, multi-factorial disease and genetic factors have been suggested to play an important role in the development and progression of cartilage degeneration. More than 50 candidate genes related to OA have been identified by genetic association studies. Unfortunately, today, science is faced with a problem of missing heritability. One way to solve this problem is to study epigenetic regulatory mechanisms, including miRNAs and their target sites. However, the most common weakness is difficulty in replicating the previous association signals. Furthermore, studies of microRNA binding sites to date are not very numerous.

rs9659030 (COL11A1)
When the A variant is present in the seed site, miR-7-1-3p and miR-7-2-3p fail to bind to the respective genes, whereas the presence of the G variant in the seed site creates new binding sites for miR-495-3p and miR-5688, respectively. The role of these miRNAs in the pathogenesis have not been examined yet, but these observations are, nevertheless, compelling. A recent EVA study by Raine et al. analyzed the allelic expression imbalance (AEI) of the COL11A1 gene using three polymorphic loci: rs9659030, rs1676486, and rs2615977. In this study, a decrease in the expression of the COL11A1 gene was detected in the presence of a rare T allele at the polymorphic locus rs1676486 (p.Ser1419Thr), which has been shown to be in linkage disequilibrium with rs9659030. AEI at rs1676486 locus is a risk factor for lumbar disc herniation (LDH), but not for OA. Therefore, this suggests that the two diseases share a common functional phenotype, namely, AEI of COL11A1, but this appears to be a disease risk only in LDH [23]. Another noteworthy study of allelic imbalance (AI) was done by den Hollander et al. Their findings were based on the AI effect observed in individuals heterozygous for both rs2615977 and rs9659030 loci. Two independent COL11A1 SNPs (D = 0.2, r 2 = 0.01) with profoundly significant meta-ϕ AI were revealed: the former was the abovementioned rs9659030 locus with meta-ϕ AI of 0.65 (P = 1.2 × 10−25) and the rs2229783 locus with meta-ϕ AI of 0.53 (P = 3.7 × 10−8) was the latter. Notably, in contrast to the results obtained by Raine et al., the extent and consistency of AI for rs9659030 were considerably higher than those for rs2229783 and rs1676486 loci [24]. There were no statistically significant associations at rs9659030 locus in regard to the development of glaucoma and myopia in the Chinese population [25,26]. However, further studies are required to confirm these findings.

rs229069 (ADAMTS5)
Unfortunately, limited literature is available on this polymorphic variant. The G allele disrupts the interaction of miR-3144-3p and miR-875-5p with the target mRNA, whereas the C allele creates a novel binding site for let-7d-3p and let-7e-3p miRNAs, which suggests a change in the regulation of expression of the ADAMTS5 gene. Our findings regarding the association of this locus with the development of OA are a significant contribution to the understanding of the molecular pathogenesis of the disease.

rs13317 (FGFR1)
The A allele at 3 UTR of rs13317 harbors binding sites for both miR-3128 and miR-4470, while miRNAs associated with the G allele at this locus are not currently defined. However, this polymorphic variant was found to be associated with the pathologies of the nervous system and bone tissue development; an association with the anomalies in the development of the facial skeleton (overall head size and midfacial development, the ratio of the distance between the eyes) was observed in patients from Japan and Korea [27]. The CC genotype occurred 50% as often in patients with ossification of the posterior longitudinal ligament, as compared to the control group in Koreans [28]. In the study of 203 patients with an injury of the rotator cuff of the shoulder, an association between the T allele and the presence of pathology was revealed (P = 0.08/0.02; OR 2.67; CC overrepresented in control patients) [29]. In a study of non-union (NU) fractures, 167 patients with long bone fractures, 101 with uneventful healing, and 66 presenting aseptic non-unions were investigated. A significant association of FGFR1 (rs13317) C allele with NU was observed. Carrying the C allele increased twice the risk of NU, in contrast to TT genotype, which acts as a protective factor [30].

rs6854081 (FGF2)
Fibroblast growth factor is one of the key links in bone metabolism. It is expressed in most mesenchymal and bone-related cells, including chondrocytes, osteoblasts, adipocytes and osteoloclasts. FGF2 is expressed in the embryos of the extremities at the developmental stage and controls the growth and formation of limb patterns. In bone cells, by regulating bone formation induced by parathyroid hormone (PTH) and bone morphogenetic protein 2 (BMP2) promotes the differentiation of bone marrow stromal cells into osteoblast [31][32][33]. It was found that rs6854081 and rs1048201 in the 3 -UTR of FGF2 are associated with bone mineral density (BMD). The minor G allele of rs6854081 was associated with low BMD in a cohort of white women (n = 2725). The study demonstrated that the FGF2 transcript containing the T allele of rs6854081 had a lower expression level in monocytes and B cells, and the allele was associated with higher BMD. T allele could have a higher binding affinity for miR-146a/b, therefore, inhibiting FGF2 translation. Suppression of FGF2 will reduce osteoclastogenesis, which in turn will lead to an increase in BMD [34]. In another study, the minor T allele of rs1048201 was significantly associated with spinal BMD in the Han people (n = 2339), it may act by affecting binding of hsa-miR-196a-3p [35]. Study of the rs6854081, rs1048201 and rs7683093 of FGF2 in one more Chinese ethnic cohort-Zhuang, revealed significant relationship of these loci with the bone density of of the femoral neck [36]. However, a study only by the Asian could not reproduce the relationship between target loci and BMD, suggesting that the effects of rs6854081, rs1048201, rs7683093 on bone regulation may be ethnic specific [37]. In this regard, results of our study on associations of the rs6554081 in particular ethnic groups are especially important.

Limitations of the Study
A limitation is that the sample size was relatively small and the study was conducted only among female patients; polymorphic target sites were not evaluated for the level of microRNA expression.

Conclusions
The current study demonstrated an association between the T allele of rs13317 in FGFR1 and incidence of the total OA and knee OA in women from Volga-Ural region of Russia. The G allele of rs6854081 in FGF2 was a risk factor for OA development for women of Tatar descent, the C allele of rs1061237 in COL1A1 was a risk factor for OA development for Russian women, and the T allele of rs229069 in ADAMTS5 and the T allele of rs73611720 in GDF5 were risk factors for OA development for women of mixed descent.
Author Contributions: A.T. conducted the clinical examination of the patients and statistical data processing, D.S. carried out the molecular genetic studies, V.P. and H.G. participated in the study design and coordination, R.K. participated in the study design and assisted in the drafting of the manuscript. All authors have read and agreed to the published version of the manuscript.

Informed Consent Statement:
The study was approved by the Ethics Review Committee of the Bashkir State Medical University, and signed informed consent was obtained from all of the participants.

Data Availability Statement:
The additional data are available in Supplementary.

Conflicts of Interest:
The authors declare that they have no competing interest.