17 variants interaction of Wnt/β-catenin pathway associated with development of osteonecrosis of femoral head in Chinese Han population

The genes of Wnt/β-catenin pathway may have potential roles in fat accumulation of Non-traumatic osteonecrosis of the femoral head (ONFH), but the effects of their variants in the pathway on ONFH development have been remained unclear. To explore the potential roles of the variants in the development of ONFH, we completed the investigation of the paired interactions as well as their related biological functions of 17 variants of GSK3β, LRP5, and FRP4 genes etc. in the pathway. The genotyping of the 17 variants were finished by MASS ARRAY PLATFORM in a 560 ONFH case–control system. The association of variants interactions with ONFH risk and clinical traits was evaluated by logistic regression analysis etc. and bioinformatics technology. The results showed that the genotype, allele frequency, and genetic models of Gsk3β rs334558 (G/A), SFRP4 rs1052981 (A/G), and LRP5 rs312778 (T/C) were significantly associated with the increased and decreased ONFH risk and clinical traits, respectively (P < 0.001–0.0002). Particularly, the paired interactions of six variants as well as eight variants also showed statistically increased and decreased ONFH risk, bilateral hip lesions risk and stage IV risk of ONFH, respectively (P < 0.044–0.004). Our results not only at the first time simultaneously showed exact serum lipid disorder and abnormal platelet function of ONFH in the same study system with the 17 variants polymorphisms of Wnt/β-catenin pathway but also shed light on the variants closely intervening the lipid disorder and abnormal coagulation of ONFH.

Non-traumatic osteonecrosis of the femoral head (ONFH) is a complex bone disorder related to genetic and environmental factors 1 .Although the pathogenesis of ONFH remains unclear, decreased osteogenesis together with increased adipogenesis is a common feature found in the bone marrow of osteonecrotic lesion since the transdifferentiation of BMSCs is regulated by the Wnt/β-catenin pathway [2][3][4] .It has been showed that critical genes of Wnt/β-catenin pathway play important roles in fat accumulation of the lesion bone marrow cavity of ONFH 5,6 , but it has rarely been reported that the effects of paired interactions of multiple variants in Wnt/β-catenin pathway genes on ONFH occurrence and development.The relevant reports suggest that chronic corticosteroid use and excessive alcohol consumption are the two main environmental risk factors of ONFH, and multiple genes are associated with ONFH risk 7,8 .In particular, some evidences revealed that the imbalance between osteogenesis and adipogenesis regulated by Wnt/β-catenin pathway played essential roles in the development of ONFH 9,10 .Moreover, much adipose tissue are often observed in the focus of ONFH so that lipid metabolism disorder has long been recognized as the core pathogenesis of ONFH 11,12 .In the view of the relationship between the multiple variants interaction of key genes in the Wnt/β-catenin pathway and the occurrence and development of ONFH has remained obscure, we selected the optimized 17 variants of GSK3β, LRP5, SFRP4, EMDR1, and

Statistical analysis
The Hardy-Weinberg Equilibrium (HWE) test and the genotype and allele distribution in the ONFH and control groups were performed for each variant by the Shesis software platform (http:// analy sis.bio-x.cn/ myAna lysis.php) 14 .
Multiple genetic models were used to assess the association of the 17 variants with ONFH risk by logistic regression analysis with adjusting for age and sex.The 2 × 2 gene interaction of multiplicative model among the 17 variants with ONFH risk, clinical stages, and hip lesions of ONFH were analyzed by logistical regression analyses, respectively.The generalized multifactor dimensionality reduction method (GMDR, version 0.7) were used to test high dimensional gene-gene interactions with ONFH risk.
The student's t-test was used to analyze the differences of age, BMI, serum lipid levels, and platelet indices between ONFH and control groups.ONE-WAY ANOVA was used to finish the different analysis of the serum lipid levels, platelet indices, and the age at onset (years old) among the genotype of each variant.Measurement data are presented as mean ± standard deviation.The χ 2 analysis was used to assess differences of gender, etiological classification, hip lesions and clinical stages among the genotypes of each variant.The above statistical analysis was performed using SPSS software (version 26, SPSS, Chicago, IL, USA).Statistical significance is defined as P < 0.05.

Functional annotation and eQTLs analysis of the 17 variants in Wnt/β-catenin pathway
The HAPLOREG v4.1 database (https:// pubs.broad insti tute.org/ mamma ls/ haplo reg/ haplo reg.php) was employed for exploring functional annotations of the 17 variants.The genetic variants can be visualized along with protein binding information by searching this bioinformatical online tool.To examine the potential functions of the significant variants, we extracted gene expression data from the GTEX database (https:// gtexp ortal.org/ home/) to investigate the tissue-specific eQTL.The research data of GTEX have explored the differences in gene expression of different genotypes in various tissues and have indicated them with statistically significant P-values (P < 0.05).

Association of the genotypes and allele frequencies of 17 variants in the Wnt/β-catenin pathway with ONFH risk
Genotype and allele frequency analysis in the 17 variants are represented in Table 1, Supplementary Table 2.

Association of the genotypes of 17 variants in the Wnt/β-catenin pathway with the clinical traits of ONFH
The results were shown in Table 2, Supplementary Table 3.The age at onset of the AA carriers of Gsk3β rs3755557 (T/A) was significantly younger than that of TT carriers in ONFH patients (P = 0.001), as well as the AA carriers of steroid-induced ONFH was statistically higher than that of idiopathic ONFH (P = 0.026).Besides, the bilateral hip lesions (9.4%) of the CC carriers of SFRP4 rs2598116 (A/C) was significantly higher than that (1.6%) of unilateral hip lesions, P = 0.026.The female (26.2%) of the GG carriers of SFRP4rs1802073 (T/G) presented statistically higher than that of the male, P = 0.009.

Lipid levels and platelet parameters between ONFH and Control groups and their association with the 17 variants polymorphism of Wnt/β-catenin pathway
Compared to control group, the decrease of HDL-c and the increace of LDL-c and LDL-c /HDL-c in ONFH patients showed statistical significance, P < 0.001, P = 0.048, and P < 0.001, respectively, see Fig. 2. The TC levels in the ONFH patients was lower than that of healthy controls because of a abnormal result for older age in healthy controls.Meanwhile, the PCT, MPV, and PDW levels of ONFH patients were also significantly increased (all P < 0.001), shown in Fig. 3.
Association of 17 variants genotypes in Wnt/β-catenin pathway with lipid levels and platelet parameters are shown in Table 3.The LDL-c /HDL-c of the TT carriers of Gsk3βs2037547(C/T) was significantly higher than that of the heterozygous CT carriers, P = 0.03, and the TC and HDL-c level in the AA carriers of Gsk3βrs334558(G/A) were significantly decreased, compared to those of the major homozygous GG carriers, P = 0.002, P = 0.005, respectively.And the TC level in the AA carriers and the LDLc/HDLc of the GG carrier of Gsk3βrs3732361 (A/G) were statistically higher than that of the heterozygous AG carriers, P = 0.033, P = 0.032, respectively.The TC level and LDL-c/HDL-c in the AA carrier of Gsk3β rs6438552(G/A) were decreased or increased, compared to that of the GG carrier and GA carriers, P = 0.040, 0.050, respectively.The HDLc level of the TT carrier in the SFRP4 rs1802073(T/G)showed significantly higher than that of heterozygous TG, P = 0.010.
Moreover, the PLT of the TT carriers of LOC105375236 rs1721400(C/T) was significantly decreased, compared to that of the CC carriers, P = 0.032, while the PDW(%)in the AA carries of Gsk3βrs334558(G/A) was obviously lower than that of the GG carrier, P < 0.001.The PDW(%)in the AA carriers of Gsk3βrs3732361 (A/G) and GG carriers of Gsk3βrs6438552(G/A) were also significantly lower than that of the heterozygous carriers, P = 0.006, P = 0.029, respectively.And the MPV(fl) in the TT carriers of LRP5rs312778(T/C) and in the TT carries of SFRP4 rs1802074(C/T) all revealed statistically higher than that of heterozygous TC and CT carriers, P = 0.018, P = 0.040, respectively.

Potential regulatory role and eQTL gene expression analysis in the 17 variants of Wnt/ β-catenin pathway
To explore the functional effects and their potential regulatory roles of the 17 variants, we performed a functional prediction of statistically significant variants by using the HaploReg database.The results exhibited that rs334558, rs3732361, rs3755557 and rs6438552 of Gsk3β and rs312778 of LRP5 revealed their potential biological functions, as shown in Table 4.The five variants all are located in 5 prime UTR of Gsk3β and LRP5 gene, respectively.As motif variants, they play a switch role in controlling transcription factor binding to the gene and affecting its transcription and mRNA expression.Besides, the variants, also as promoter histone marks and DNAse function, exert multiple effects.The eQTLs analysis results also show that the 5 variants were closely associated with the modulation of gene expression because they are identified as eQTLs in different tissues expressions with significant p-values.The significant effects of different genotypes on the gene expressions in the above five eQTs variants are shown in Fig. 4. The expression of minor homozygous GG genetype of eQTLrs 334558 of Gsk3β is the highest among three genotypes in muscle skeletal, while the major homozygous AA genotype in eQTLrs 2732361 of Gsk3β has the highest expressions in subcutaneous adipose tissue.Moreover, the expressions of minor homozygous AA and heterozygous TA genetypes of eQTLrs 3755557 of Gsk3β are the higher than that of major homozygous TT genotype in muscle skeletal.And the minor homozygous GG genetype of eQTLrs 6438552 of Gsk3β and the minor homozygous CC genetype of eQTLrs312778 of LRP5 also reveal the highest expressions among three genotypes in whole blood and adrenal grand, respectively.

Discussion
The multiple functions of Wnt signaling pathway include mitogen stimulation, cell fate determination, and cell differentiation 15 .Of them, Wnt/β-catenin pathway affects bone formation and remodeling processes 16,17 .
The human Wnt family is composed of 19 highly conserved genes that regulate gene expression, cell behavior, and cell polarity.GSK3β located in 3q13.33 is a protein coding gene, and it is the initiation factor of the Wnt/β-catenin pathway and as a downstream protein directly interacting with AKT in the PI3K/Akt signaling pathway 18 .Considering the critical roles of GSK3β in bone metabolism, we focused on the effects of its five variants on the development of ONFH.Our results showed that the AA genotype and the A allele frequency of Gsk3βrs334558 (G/A) were significantly associated with the increased ONFH risk, as well as the genetic models of co-dominants, dominants, and recessives of the variant all were related to the increased risk.Moreover, the paired interaction of Gsk3β rs334558 with LRP5 rs2306862, with LRP5 rs3736228, and with SFRP4 rs1802074 also showed a significantly increased ONFH risk, respectively, whereas the paired interactions between Gsk3β rs3755557 and SFRP4 rs2598116 revealed a decreased ONFH risk.ONFH is a complex disease involved in genes and environmental risk factors.Even if the identified variants from genome-wide association study(GWAS) often show only modest effects on the disease risk, leading to the "missing heritability".The gene-gene interactions (epistasis) analysis can solve a part of this "missingness" to thereby elucidating their effect on complex diseases 19 .The paired interactions among the variants of GSK3β gene itself or with other gene variants strongly revealed the contribution of epistatic effects to ONFH risk, that is, recovered the missing heritability of ONFH.GWASs has  www.nature.com/scientificreports/been reported that the "missingness" is attributed to genetic heterogeneity, epistasis (gene-gene interaction), and gene environment interaction 20,21 .Moreover, it has been investigated that epistasis could drive the evolution of recombination frequencies among genes on the same chromosome, thereby altering gene order 22 .The gene-gene interactions include synthetic, suppressive, and epistatic types.When Gsk3β rs334558, together with the other multiple variants from residing in Wnt signaling pathway, may produce a new trait with synthetic and epistatic types, increased or decreased ONFH risk 23 .
Our results enough revealed that Gsk3β rs334558, as an essential variant, involved in the development of ONFH.What is the possible mechanism of this variant leading to ONFH? First, rs334558 is a 5 prime UTR variant, as motif variant, it is a switch that controls transcription factor binding to Gsk3β gene and affecting its transcription.Second, rs334558, as a eQTL variant, may play potential regulatory role in Gsk3β gene expression.Besides, the variant, also as a promoter histone mark and DNAse function, exerts multiple effects.In the HaploReg database, we performed a functional prediction of the significant variants.And the results showed rs334558, rs3732361, rs3755557 and rs6438552 of Gsk3β all exhibited potential biological functions in Table 3. Association of the genotypes of 17 variants in Wnt/β-catenin pathway with lipid levels and platelet indices of ONFH.TG, triglyceride(0.56-1.71mmol/L) # ; TC, total cholesterol (2.90-5.17mmol/L) # ; HDL-c, high-density lipoprotein cholesterol (0.90-1.68 mmol/L) # ; LDL-c, low-density lipoprotein cholesterol (2.70-3.40mmol/L) # .PLT, platelets (125-350/L × 10 9 ) # ; PCT, Plateletcrit (0.108-0.282%) # ; MPV, Mean platelet volume (7.0-11.0↑fl)# ; PDW, platelet distribution width (9.0-17.0%)# .# Reference value; P: one-way ANOVA (or Tukey test).gene modulation.The variants were identified in different tissues expression with significant p-values.Using the GTEx database,we also find that the four variants regulate the functional importance of Gsk3β gene expression as eQTLs.The variants potential functions of GSK3β discovered by the databases analysis may explain their molecular mechanism during the development of ONFH.SFRP4 acts as soluble modulator of Wnt signaling and exerts critical roles in regulating bone formation and resorption as the main receptor for Wnt [24][25][26] .In our evaluation results of the six variants of SFRP4, the rs1802074 showed the increased ONFH risk by the epistasis with Gsk3β rs334558, and the increased bilateral hip lesions risk in the CC carriers of SFRP4 rs2598116 (A/C).The recessives model of SFRP4 rs1052981 was reported to relate to increased ONFH risk and knee OA in women 27 .Previous results had also revealed a significant association between SFRP4 rs1802073 and BMD in Danish and Belgian men 28 , and the variant also connected with elevated serum lipid levels after acitretin treatment 29 .Our results also showed that the female percentage of minor homozygous GG of the variant was statistically higher than that of the male.The close association of the SFRP4 variants with ONFH risk and its clinical traits indicated their potential biological roles in the development of ONFH.

Gene
LRP5 also acts as a co-receptor with Fz protein family members for transducing signals by Wnt proteins.It plays a key role in skeletal homeostasis, and the many diseases related bone density are caused by mutations in this gene.This study involved the four variants of LRP5, and both the co-dominant and dominant model of LRP5 rs312778 (T/C) showed a protective role in decreased ONFH risk, while the LRP5 rs2306862 and rs3736228 connected with the increased ONFH risk by the paired interaction with Gsk3β rs334558.With using the HaploReg database, the LRP5 rs312778 was also found to played the roles of promoter histone marks, DNAse, and proteins bound.By GTEx database, the variant was further identified as an eQTL of regulating gene expression, exhibited its potential biological roles in the development of ONFH.As Wnt co-receptor, LRP5 has been identified as a crucial protein for mechanical signaling in bone, and the modifiers of LRP5 activity play important roles in mediating signaling efficiency 30 .LRP5 gene has been focused on GWAS and has confirmed a close relationship of its genotypes with osteoporosis 31,32 .Therefore, the associations between the variants of LRP5 gene and BMD or osteoporotic fractures have been extensively explored.Rs4988300 and rs634008 of LRP5 gene also associated with bone phenotypes in the elderly with OP or fractures 33 .To our knowledge, the association of the four variants of LRP5 with the development of ONFH was for the first time reported in this study.
LOC105375236 is a lncRNA located in 7p14.1, and the rs1721400(C/T), previously thought to be a variant of GSK3β gene, and now belonging to LOC105375236, is a genetic downstream transcript variant.In a polygenic inheritance GWAS study, the intergenic variant rs1524058 locates between LOC105375236 and STARD3NLwas identified as one of pleiotropic variants related to all four traits 34 .In our present study, the one of best model composed of 7 variants containing rs1721400 (C/T) in GMDR analysis connected with ONFH risk tendency.The protein encoded by EPDR1 gene is a type II transmembrane protein similar to two families of cell adhesion molecules 35 .EPDR1was reported to be closely related to the developmental processes of various tumors 36-38.Our results revealed that the paired interactions between EPDR1 rs16879765(C/T) and Gsk3β rs334558 significantly increased stage IV risk of ONFH.And we have also introduced the GMDR approach to our study and the results revealed four best models related to ONFH risk were set up.Even though failed to reach significance (all P-values = 0.0547), GMDR analysis is a useful tool to predict ONFH risk in expanded samples.www.nature.com/scientificreports/Our results also revealed that the age at onset of AA carriers of Gsk3β rs3755557 (T/A) was significantly younger than that of TT carriers, as well as the proportion of steroid-induced ONFH of the AA carriers was statistically higher than that of idiopathic ONFH.Moreover, the paid interaction between Gsk3β rs3755557 and LRP5 rs2306862, LRP5 rs3736228, and LRP5 rs556442 showed significantly decreased bilateral hip lesions risk, while the interactions between Gsk3β rs334558 and EPDR1 rs16879765 revealed the increased stage IV risk of ONFH, which indicate that the interaction of the variants strongly involved in the development of ONFH.
Fat accumulation is often observed in the lesion bone marrow cavity of ONFH in the animal experiments as well as the clinical studies 39,40 .Here, we clearly showed not only the increase of serum LDLc, LDLc/HDLc and the decrease of HDLc but also the rise of PCT, MPV, and PDW parameters of ONFH patients, which at the first time provided forceful evidence that there were simultaneous lipid metabolic disorders and coagulation disturbance in the same research system of focusing on the roles of 17 variants in Wnt/β-catenin pathway in ONFH development.Fat embolism and coagulation disorders are considered the two possible causes of interruption of blood supply in ONFH 41 .It is well-known that platelet indices are closely associated with platelet function.MPV indicates early platelet activation, and the increased MPV shows that the release of thromboxane A2 from large platelets is more than from smaller ones 42 .PCT (PLT × MPV) reflects the volume of blood occupied by platelets.As an index of the size uniformity of PLTs, PDW is the marker of platelet activity, and the increased PDW level are associated with abnormal thrombosis 43 .Especially, the variants in Wnt /β-catenin pathway were significantly involved in the lipid metabolic disorders and abnormal coagulation of ONFH by different genotypes and genetic model, which strongly indicates that the possible roles of the variants in the hyperlipidemia as well as the hypercoagulation of ONFH.The pathological features of ONFH are affected by impaired MSCs differentiation, and the Wnt/β-catenin pathway plays a vital role in the lineage specification during MSCs differentiation. .Some studies also show steroid-and alcohol-induced ONFH from the inhibition of MSCs pathway.Our previous works found that a molecular axis driven by circRNA CDR1as-miR-7-5p-WNT5B and a molecular axis driven by lncRNA HOTAIR-miRNA-217-DKK1 played essential roles in the transdifferentiation between osteogenesi and adipogenesis of BMSCs in ONFH 44,45 .Both the molecular axes vigorously promoted the adipogenesis and inhibited the osteogenesis of the BMSCs.Especially significant, their ending protein molecules, WNT5B and DKK1, respectively, also are as key proteins of Wnt/β-catenin pathway, which presented the evidence of the abnormal molecular damage network driven by non-coding RNAs and eventually collected the Wnt/β-catenin pathway, elucidated the molecular mechanism of ONFH fat accumulation and the molecular targets for early prevention and treatment of ONFH.The evidences were also highly consistent with the results of multiple variants in Wnt/β-catenin pathway involving in the development of ONFH.Therefore, the variants in Wnt/β-catenin pathway facilitate specific biological events, the lipid and coagulation disorder, the differentiation of BMSCs, and the abnormal molecular damage network leading to the development of ONFH.
Our study have several limitations.First, the samples size of 560 case-control was smaller, which may have limited our statistical power to detect slight differences between groups, and further expanded samples investigation will improve our statistical results.Second, there was a difference in age and gender between ONFH and control group due to the older age and more female of control group, which may also affect the statistical results because of age and gender as covariates of analysis.Future research will overcome the bias through precise selection of enrolling participants.Third, although the serum lipid levels and platelet parameters between ONFH and control group revealed statistical significance, the actual difference values between two groups were smaller.A accurate grouping or stratification of enrolled ONFH patients in future research may improve such studies.Fourth, the molecular biological roles of these variants in the development of ONFH remain to be deeply identified by multiple complex verification experiment.
Fifth, although our study revealed the effects of Wnt pathway variants on ONFH development involving in lipid metabolism disorders and platelet parameters changes, the exact roles and its mechanisms of Wnt /β-catenin pathway in lipid metabolism and platelet regulation of ONFH need to be in-depth explored in future study..

Conclusions
Our results at the first time revealed that the genotypes, allele frequency, and genetic models of Gsk3β rs334558, SFRP4 rs1052981, and LRP5 rs312778 in Wnt/β-catenin pathway were associated with increased and decreased ONFH risk, respectively, and the paired interactions of six variants and eight variants in the pathway were statistically connected with the increased or decreased ONFH risk, respectively.The interactions of six paired variants were also related to the decreased risk of bilateral hip lesions of ONFH, while the paired interactions of two variants closely connected with the increased stage IV risk of ONFH.In particular, we showed that the crucial variants in Wnt/β-catenin pathway were involved in the lipid disorder and abnormal coagulation of ONFH and played significant roles during the occurrence and development of ONFH.

Figure 1 .
Figure 1.Association of the paired gene-gene interactions of 17 variants in Wnt/β-catenin pathway with the risk, hip lesions, and clinical stages of ONFH.The chromosomes are arranged end by end in a clockwise direction.The exterior of the circle is chromosome number and scale.The inner circle represents the gene where the variant is located.The innermost part is ID of the variant.The lines connecting two variants represents paired gene-gene interactions.The red lines represent interactions with OR > 1 while the blue lines represent interactions with OR < 1.In addition, the color shade and lines thickness indicate whether the interaction is statistical significance: the thinnest and lightest lines represent interactions with no statistical significance (P > 0.1), the middle ones represent interactions with 0.05 < P < 0.1, and the thickest and the darkest lines represent interactions with P < 0.05.(a) The paired gene-gene interactions of 17 variants in Wnt/β-catenin pathway with the ONFH risk.(b) The paired gene-gene interactions of 17 variants in Wnt/β-catenin pathway with the unilateral and bilateral hip lesions of ONFH.(c) The paired gene-gene interactions of 17 variants in Wnt/β-catenin pathway with the clinical stages of ONFH.Gsk3β, glycogen synthase kinase 3 beta; LRP5, LDL receptor related protein 5; SFRP4, secreted frizzled related protein 4; EPDR1, ependymin related 1; LOC105375236, uncharacterized LOC105375236, LncRNA.

Figure 4 .
Figure 4.The effects of genotypes of eQTLs variants in Wnt/β-catenin pathway on gene expression in different tissues.