Prevalence of polymorphisms in OPG, RANKL and RANK as potential markers for Charcot arthropathy development

Charcot arthropathy is one of the most serious complications of diabetic foot syndrome that leads to amputation of the affected limb. Since there is no cure for Charcot arthropathy, early diagnosis and implementation preventive care are the best available treatment. However, diagnosis is hindered by obscure clinical picture of the disease and lack of molecular markers for its early detection. Results of recent research suggest that OPG-RANKL-RANK axis regulating bone metabolism can be associated with Charcot arthropathy and that SNPs in OPG gene are associated with the disease. Here we report the results of comprehensive analysis of ten SNPs in OPG, RANKL and RANK genes in 260 subjects divided into diabetes, neuropathy and Charcot arthropathy groups. Besides genotype analysis we performed linkage disequilibrium and hierarchical clustering to obtain information about correlation between SNPs. Our results show that OPG 245T/G (rs3134069) and OPG 1217C/T (rs3102734) polymorphisms co-occur in patients with Charcot arthropathy (r2 = 0.99). Moreover, hierarchical clustering revealed a characteristic profile of all SNPs in Charcot arthropathy and neuropathy, which is distinct from control group. Our results suggest that analysis of multiple SNPs can be used as potential marker of Charcot arthropathy and provide insight into possible molecular mechanisms of its development.

RANK. Binding of RANKL to OPG results in the inhibition of bone resorption and stimulates bone mass building. The dynamic equilibrium between RANKL and OPG concentrations is crucial for normal bone metabolism 12 . Alternatively, the imbalance of the RANKL/OPG ratio could lead to an uncontrolled loss of bone mass 13 .
Diabetes-induced distortion of the RANKL/OPG ratio can directly affect bone formation and regeneration, i.e. processes that are likely impaired in CN 7 . Moreover, altered levels of RANKL and OPG have been linked to the development of vascular calcification, which has been proposed as one of the causes of this disease 9,14 . A recent study indicated that the observed bone destruction in CN could be partially caused by an increase in the RANKL gene expression and subsequent local inflammation in the affected limb 15 . It has been also proposed that an increased bone resorption and calcification of blood vessels could lead to the development of Charcot arthropathy or intensify the development of this disease 9,15,16 . Therefore, increased RANKL/OPG ratio can be both, a symptom and a cause in the development of the Charcot arthropathy.
Recent studies revealed CN-specific occurrence of certain polymorphisms in OPG, RANKL and RANK genes [17][18][19] . Therefore, the aim of this study was to determine the usefulness of selected SNPs as potential genetic markers of CN, allowing to assess the risk of CN development in patients with diabetes. We have analysed ten polymorphisms in the OPG, RANKL and RANK genes in patients with CN, diabetic neuropathy (N) and diabetes without symptoms of neuropathy or CN (D) (blood samples were collected and analysed from 77 (CN), 77 (N) and 106 (D) patients). The study design was focused on finding difference in the specific SNP occurrence between patients with the diabetic neuropathy and Charcot arthropathy, in order to find genetic markers of the disease. In addition, hierarchical clustering analysis of all analysed polymorphisms in patients provided general overview of genetic variants associated with Charcot arthropathy. In parallel to genetic analysis, serum level of RANKL and OPG was measured to gain comprehensive view on Charcot arthropathy in relation to neuropathy and diabetes. Obtained results show that combination of the genetic analysis with assessment of the RANKL and OPG serum levels could be used as a prediction marker for the development of CN in diabetic patients.

Results
Analysis of selected polymorphisms in OPG, RANKL and RANK genes. In this study, we analysed ten polymorphisms in three genes: OPG (five variants), RANKL (three variants) and RANK (two variants), in samples collected from 260 patients. Results of this analysis are summarized in Table 1. For three out of ten analysed loci, namely OPG 950T/C, RANK 421C/T and RANK 575C/T we found no significant changes in the genotype and allele frequencies across the studied groups. Out of remaining SNPs, we observed the strongest change for the C allele frequency at OPG 1181C/G in patients with CN. In patients with CN, the frequency of the CC genotype was higher than in patients in the N and D groups (CN-36%, N-15%, D-23%, Fig. 1A). Statistical significance of observed differences was tested by the chi-square test (p < 0.001 and p < 0.05 for CN vs. D, and CN vs. N, respectively). Similarly, we observed a bias in allele and genotype frequency towards N and CN groups for the OPG 6890A/C allele (Table 1 and Fig. 1A). At this locus the frequency of the C allele was increased in N and CN groups, what corresponds with higher occurrence of the CC genotype. In case of OPG 245T/G and OPG 1217C/T alleles the analysis revealed a higher frequency of the G and C alleles, respectively (Table 1). Interestingly, genotypes distribution at the OPG 245T/G and 1217C/T SNPs was identical among the phenotypic groups (Fig. 1A). The TT homozygote was the most dominant in both loci for all groups, whereas the GG and CC homozygotes (at OPG 245T/G and 1217C/T, respectively) were detected only in the diabetes group. At both loci, we observed that the frequency of heterozygotes is the highest in the N group. The analysis of RANKL 290C/T, 643C/T and 693G/C polymorphisms also revealed differences in the allele and genotype frequencies among the studied groups (Table 1 and Fig. 1B). The T allele in RANKL 290C/T and 643C/T, as well as the C allele in 693G/C, were more frequent in the N and CN groups. It correlates with the increased frequency of the corresponding homozygotes in those groups (TT for 290C/T and 643C/T, and CC for 693G/C). Statistical analysis of the genotype frequencies indicated that N and CN groups differ significantly from the D group (Fig. 1B).
Linkage disequilibrium analysis. The analysis of SNPs revealed that some alleles and genotypes have similar patterns of distribution in studied groups. In order to determine if these patterns are specific for a given condition, we calculated linkage disequilibrium (LD) for the analysed SNPs in the OPG and RANKL genes (Figs 2 and 3, respectively). For OPG gene polymorphisms we observed the highest LD between OPG 245T/G and OPG 1217C/T for the CN group (r 2 = 0.99) and lower values for D and N groups (r 2 = 0.58 and r 2 = 0.34 respectively). A weaker association was observed between OPG 1181G/C and 950T/C in CN group (r 2 = 0.51), although it was still higher than in D and N groups (r 2 = 0.24 and r 2 = 0.04 respectively).
The same analysis was performed for the RANKL gene polymorphisms and has revealed the highest LD to be between RANKL 693G/C and 290C/T polymorphisms in the CN group (r 2 = 0.89, Fig. 3) but in other groups the their association was also high (r 2 = 0.60 and r 2 = 0.72 in D and N group respectively). In case of RANKL 693G/C and 643C/T we observed the lowest association in CN group (r 2 = 0.52) while it was the highest in D group (r 2 = 0.69). The LD analysis has been also performed between OPG and RANKL polymorphisms but did not show disequilibrium at any of the studied groups.
Hierarchical clustering of SNPs in the OPG, RANKL and RANK genes. Human diseases rarely depend on single genetic variants and are likely a cumulative effect of multiple changes in different genes. Therefore analysis of the occurrence of multiple SNPs in patients is more suitable for finding genetic basis of a disease. In order to obtain a general overview of genotype distribution in the Charcot arthropathy, we have performed hierarchical clustering (HC) analysis for all analysed polymorphisms in the studied group of patients (Fig. 4A). The genotype data has been transformed to generate a matrix where each row represents a patient and each column represents a SNP. Both, rows and columns, have been sorted on the basis of the genotypes occurrence similarity, resulting in clusters depicted as dendrograms. The major advantage of HC is that it shows similarities between patients (rows) and SNPs (columns). Consistently with the LD analysis (Fig. 2B), OPG 245T/G and OPG 1217C/T have nearly identical pattern of distribution. However, RANK 421C/T and OPG 6890A/C also cluster together with OPG 245T/G and 1217C/T, which was not detected by pairwise comparison using only LD. All three RANKL SNPs analysed here, have been grouped together and patterns of genotypes at these positions fit well with the three distinguished clusters. Analysis of the patient study group designation revealed that the three clusters are comprised in different proportions of D, N and CN patients (Fig. 4B). The first cluster is composed mostly of the D patients (40.85%), whereas in the second and third cluster, the D patients comprise 26.83% and 20.56%, respectively. The same proportion of N and CN patients (36.59%) characterizes the second cluster. The highest number of patients with CN has been observed in the third cluster (50.47%) (Fig. 4B).  Cytokine concentration in the blood serum. It has been reported that in the CN patients an increased ratio of RANKL to OPG in the blood serum is observed in respect to healthy subjects. Therefore, we have determined levels of these cytokines in the blood serum of the three examined groups of patients, using commercial ELISA assays. We observed increased levels of RANKL in CN and N groups comparing to D group (1.01 ± 1.45 pmol/l, 2.66 ± 1.74 pmol/l and 0.5 ± 0.43 pmol/l respectively, Fig. 5A). Surprisingly, the highest concentration of the protein was in N group, what has not been observed before. Observed differences in OPG concentration were statistically significant as determined by T-test (p < 0.01 for CN vs. N and CN vs. D, p < 0.001 for N vs. D).
The OPG levels were higher in patients in the CN and N groups, when comparing to the D group (7.36 ± 4.1 pmol/l in CN, 6.29 ± 1.68 pmol/l in N and 4.77 ± 2.38 pmol/l in D, Fig. 5A). These differences were statistically significant as indicated by T-test (p < 0.01 for N vs. D and p < 0.001 for CN vs. D), while there was no difference between CN and D groups. The proportion of RANKL to OPG was the highest in the N group (0.42) and was approximately 4-fold higher than in D group (0.10), whereas in the CN group (0.14) the ratio of RANKL to OPG was approximately 1.5-fold higher than in the D control (Fig. 5B).

Discussion
The pathogenesis of Charcot arthropathy is still unknown but research on the OPG/RANKL/RANK cytokines axis has highlighted this system's role in the bone-associated diseases. In this work, we have performed a comprehensive genetic analysis of ten polymorphic loci in the 260 patients, divided into three groups (106 CN, 77 N and 77 D). To our knowledge, this is the largest group of patients studied in the context of Charcot arthropathy and neuropathy in diabetes. In parallel to classical genotyping, we present here hierarchical clustering as a useful analysis in the studies of SNP association with human disease. Our results suggest association of the RANKL gene polymorphisms with N (neuropathy) and CN (Charcot arthropathy). Moreover, we suggest that increased ratio of RANKL/OPG in the blood serum is specific for neuropathy and can be a factor leading to the development of Charcot arthropathy.
To date, studies addressing a possible link between genotype and Charcot arthropathy have been limited to the analysis of the OPG gene polymorphisms 17,18 . In the study presented here, we have analysed additional SNPs that have been associated with altered bone metabolism in postmenopausal women 20,21 . In total, out of ten selected SNPs three did not show any association with either of the studied conditions. Our results indicate that OPG 950T/C, and RANK 421C/T and 575C/T, which have been previously analysed in the context of Paget's disease and osteoporosis, are not associated with CN 21,22 . For the remaining seven SNPs, four in the OPG gene (245T/G, 1181G/C, 1217C/T and 6890A/C) and three in the RANKL gene (290C/T, 643C/T and 693G/C), we have found statistically significant differences in the genotype distribution in the studied groups.
Consistently with previous reports, the strongest association with CN was observed for the OPG 1181G/C polymorphism 17,18 . Moreover, although the hierarchical clustering did not show any characteristic OPG 1181G/C genotype pattern for a particular patient cluster, the majority of patients with the CC genotype belongs to the CN and N groups (Fig. 4A). Interestingly, OPG 1181G/C is located in the coding sequence of this gene and point mutations result in the codon change from lysine into asparagine ( Fig. 2A). This amino acid substitution could potentially influence the protein activity, explaining its association with CN. Intriguing is also the co-occurrence of OPG 1181G/C and 950T/C in CN. The later polymorphism is located upstream from transcription start site and its role in the disease and in association with OPG 1181G/C is unknown. It is also challenging to propose a mechanism explaining observed association between OPG 245T/G and 1217C/T with CN and N (Fig. 2B). Both SNPs are located in the noncoding region of this gene, with OPG 245T/G located in the putative promoter region of the OPG gene and 1217C/T located in the first intron ( Fig. 2A). This location, especially that of OPG 245T/G, would seem to suggest altered gene expression as the underlying cause of the observed phenotypes, however, since little is known about molecular mechanisms underlying OPG gene expression, it is unknown if SNP variants at OPG 245 and 1217 can influence this process. Yet, both loci show strong linkage disequilibrium and very similar pattern of genotype distribution in the studied patients (Figs 2B and 4A).
In contrast to previous studies 17, 18 , we aimed to find differences between CN, N and the D group, which was used as the control group. This distinction allowed for finding that OPG 6890A/C is associated with both, CN and N. However, as in the case of OPG 1217C/T, this SNP is located in the intron sequence and it is unclear how this position contributes to the development of CN and N. Furthermore, OPG 6890A/C was clustered together with OPG 1217C/T, OPG 245T/G and RANK 421C/T polymorphisms, which did not show distribution patterns (Fig. 4A, upper dendrogram).
Our results also show that RANKL 290C/T, 643C/T and 693G/C polymorphisms are associated with both, CN and N (Fig. 1B). All three polymorphisms show similar distribution across all patients and form a separate cluster in the HC analysis (Fig. 4A, upper dendrogram). The general distribution pattern of these polymorphisms fits well with cluster 1 and clusters 2 and 3. Cluster 1, where diabetic patients are the most abundant group (40.85%, Fig. 4B) is mostly composed of genotypes CC, GG, CC for RANKL 643C/T, 693G/C and 290C/T polymorphisms respectively. On the contrary, in clusters 2 and 3 where Charcot patients comprise 36.59% and 50.47% respectively, heterozygotes and homozygotes TT, CC and TT for aforementioned SNPs are the most frequent. This suggests that all studied RANKL polymorphisms can be associated with either CN alone or CN and N. Noteworthy, all three SNPs are located in the putative promoter region of RANKL what raises the possibility that sequence variation at these loci affects the gene expression. As revealed by LD analysis, there is an association between RANKL 693G/C and 290C/T co-occurrence in CN (Fig. 3B). On the contrary, the other pair of RANKL SNPs, 643C/T and 693G/C shows decreased association in CN in respect to D and N groups. However, it is only a hypothesis and further studies are required to determine if SNPs combination affects the expression of the RANKL gene. It has been proposed that bone loss observed in CN is a result of inflammation-induced RANKL expression 22 . It would be interesting to test whether SNPs in the promoter region enhance stimulation of RANKL expression by the tumour necrosis factor α (TNF-α).
Although it is unknown if SNPs can influence RNAKL expression, it has been reported that CN patients show increased levels of RANKL protein in the blood serum 23 . In agreement with this report, we found increased level of RANKL in patients with N and CN. However, in our analysis RANKL concentration is mostly increased in the blood serum of N group, resulting in nearly 4-fold higher RANKL/OPG ratio than in the diabetic patients ( Fig. 5A and B). It is likely that the elevated level of RANKL during the neuropathy stage is the major factor responsible for the bone loss observed in the Charcot arthropathy. Interestingly, we observed that the neuropathy group in our study showed lower RANKL/OPG ratio, due to high level of OPG which can counterbalance bone loss triggered by RANKL. Currently it is unknown if this observation is of clinical significance and can be a marker of neuropathy or reflects a stage during Charcot arthorpathy development.
Based on the results reported here, it is possible that at least two molecular mechanisms underlie the development of CN. The first one involves the OPG protein, which is required for quenching serum-soluble RANKL. Changes in the protein sequence can affect its activity and result in increased bone resorption due to RANKL-induced osteoclastogenesis. A second mechanism can be directly related to the RANKL gene expression, which can be influenced by the DNA sequence in the proximity of the promoter. SNPs in this position could directly alter the RANKL gene expression, or could enhance response to the TNF-α-mediated RANKL expression activation. Both mechanisms are likely and are supported by observed, increased levels of RANKL in the blood serum. However, regulation of human genes expression rarely depends solely on promoter activity and involves additional levels of regulation, including among others DNA methylation and post-transcriptional regulation by miRNA, which may underlie the phenotypes observed in diabetes and its complications 24 . Recently, it has been shown that overall DNA methylation level is decreased in fibroblasts derived from diabetic foot patients 25 . As Charcot arthropathy is a consequence of diabetic foot syndrome, it is likely that its development is also related with changed methylation pattern. Although results of methylation status in fibroblasts indicate that it affects genes associated with wound healing 25 , it would be interesting to analyse methylation status in bone cells and correlate methylation patterns with gene expression analysis to obtain a comprehensive view on the epigenetic mechanism of Charcot arthropathy development. Integrating transcriptome analysis in such study should include miRNA analysis, since non-coding RNAs are often factors in human disease 26 . Currently there is no data on changes in transcriptome in Charcot arthropathy or on activity of non-coding RNAs in this particular disease. However, there are reports indicating that individual miRNAs can affect levels of RANKL and RANK mRNAs. It has been shown that upregulation of miR-18a results in decreased level of RANKL mRNA in hepatocellular carcinoma 27 . Recently, miR-503 has been shown to regulate RANK levels and down regulation of miR-503 results in bone loss due to increased osteoclastogenesis 28 . Considering the above, it is likely that miRNAs play a role in the development of Charcot arthropathy and testing this hypothesis should be one of the next steps in dissecting mechanisms of this condition.
In summary, our work provides an insight into genetic changes associated with Charcot arthropathy. Analysis of loci presented in this work could be potentially used to predict risk of Charcot arthropathy development in patients with diabetes, especially in multiplex analysis. We also suggest possible molecular mechanisms underlying Charcot arthropathy development. However, verification of these mechanisms will require additional studies, integrating SNPs, methylation pattern and transcriptome analyses which could lead to better understanding and preventing Charcot arthropathy.

Methods
Study participant's characteristics. This study was performed with the approval of the Medical University in Gdansk Research Ethics Committee and all the participants gave written informed consent. It was carried out in accordance with the Declaration of Helsinki ethical principles on human studies and all the experiments was conducted with the approved guidelines and regulations. The study group consisted of patients of the Diabetology Clinique at the Medical University of Gdansk. Patients were recruited to the study on the basis of diagnosed Type I or Type II diabetes, in accordance with applicable criteria of the Polish Diabetes Association.
Patients were classified into study groups based on the occurrence of clinical features, characteristic of Charcot arthropathy, neuropathy or diabetes. Qualification for Charcot arthropathy group required the presence of the following symptoms: unilateral redness, edema of the foot, increased foot temperature compared to the other limb and radiological changes. Other diseases that have similar symptoms, like gout, deep vein thrombosis or inflammation of the bone, were excluded. Patients were classified into the neuropathy group based on the sensorimotor neuropathy examination. Patients diagnosed with other potential causes of damage to the peripheral nervous system, such as Vitamin B12 deficiency, use of neurotoxic drugs, hereditary neuropathy, abuse of alcohol or heavy metal poisoning, have been excluded. A group of people with diabetes consisted of patients without neurological changes or other disqualifying conditions. SNPs Genotyping. Three ml of blood were collected from the cubital vein into the EDTA containing tubes (BD Vacutainer K 2 EDTA) from each patient. After collecting the blood, the serum was separated from whole