Association of Variant Interactions in RANK,RANKL,OPG,TRAF6,and NFATC1 Genes with the Development of Osteonecrosis of the Femoral Head

Abstract

Multiple gene polymorphisms have been demonstrated to correlate with the susceptibility to osteonecrosis of the femoral head (ONFH). However, as a complex disease induced by multiple genes, the development of ONFH has rarely been reported to involve in gene interaction. In this study, we first explored the association of 10 variants interactions in receptor activator of nuclear factor-kappa B (RANK), RANK ligand (RANKL), osteoprotegerin (OPG), tumor necrosis factor receptor-associated factor 6 (TRAF6), and nuclear factor of activated T cells cytoplasmic 1 (NFATC1) genes with the development and clinical phenotypes of ONFH in a 377 ONFH case–control study with using Mass ARRAY^® platform. Our results showed that not only a total of 6 interactional variants in the paired 10 variants interactions were significantly associated with the development of ONFH (OPG rs2073617 and NFATC1 rs754093, p < 0.019; OPG rs2073618 and NFATC1 rs754093, p < 0.008; OPG rs2073617 and RANKL rs1054016, p < 0.039, respectively) but also a total of 4 paired interactional variants were found to involve significantly in the increased risk of bilateral hip lesions in ONFH (OPG rs2073617 and TRAF6 rs5030411, p = 0.044; RANK rs884205 and TRAF6 rs5030411, p = 0.045, respectively). Moreover, the results from generalized multifactor dimensionality reduction also showed that the five best models were identified and associated significantly with ONFH risk, p = 0.001, 0.01, 0.01, 0.01, and 0.01, respectively. Our results first suggest that the variants in RANK/RANKL/OPG pathway genes affected the development of ONFH in gene interaction manner through the interaction of the paired variants and multiple variants.

Introduction

Osteonecrosis of the femoral head (ONFH) is thought as a complex disease caused by genetic and environmental factors (Song et al., 2017; Hao et al., 2018). The prevalence of ONFH has been increasing during the past decades throughout the world, and ∼20,000 people in the United States, and 2200 people in Japan are newly diagnosed as ONFH annually, respectively (Kubo et al., 2016; Adesina et al., 2017). The estimated ONFH cases were 8.12 million among Chinese people 15 years of age and over (Zhao et al., 2015). ONFH has a significant economic impact due to largely affecting persons in the prime of life (peak age 35 years). It has been shown that a series of genetic polymorphisms were involved in the development of ONFH and that excessive alcohol and steroid are recognized as major environmental risk factors (Zhao et al., 2015). These genes, forming linkage disequilibrium (LD) or affecting the expression of other susceptible genes, may play a synergistic (Song et al., 2017) role in the development of ONFH. However, the effect of the gene interactions among these susceptible genes on the development of ONFH remains unclear. Osteoprotegerin (OPG)/receptor activator of nuclear factor-kappa B ligand (RANKL)/receptor activator of nuclear factor-kappa B (RANK) signaling pathway plays a vital role in osteoclastogenesis to maintain bone balance (Wu et al., 2017). RANKL, which is expressed primarily by osteoblasts and bone marrow stromal cells and occurs membrane bound on the surface of these cells, exerts an important effect on the activation of osteoclasts (Kong et al., 1999). RANK, the receptor for RANKL, is significant for osteolysis (Li et al., 2000). The expression of OPG and RANKL has been shown to relate to the regulation of the osteoclast differentiation controlled by osteoblasts (Goto et al., 2011).

Increasing evidences have also shown the association of susceptibility genes with the development of ONFH (Baek et al., 2017; Song et al., 2017). But the genetic contribution of the identified variants is limited and a large amount of genetic variants still remain unexplained, which is known as “missing heritability.” Several innovative approaches have been explored to find the missing heritability. Gene–gene interactions, also called epistasis, are one approach that has the potential to explore the missing genetic contribution (Meng et al., 2017). An epistasis analysis is capable of finding genes with nonapparent or weak effects, when the gene–gene interactions are significant and have the potential to discover novel genes (Zuk et al., 2012). Our previous study has revealed that a few genetic variants in RANKL/RANK/OPG signaling pathway are associated with the development of ONFH (Song et al., 2017). However, it has never been reported on the effects of gene–gene interactions in this pathway on the development of ONFH. To explore the potential missing heritability in the pathway genes, we further analyzed the roles of gene–gene interactions in 10 variants of RANKL/RANK/OPG pathway genes in the development of ONFH in a 377 ONFH case–control study.

Materials and Methods

Participants

This is a retrospective case–control study with Level III of evidence. The case–control study included 200 unrelated patients with ONFH and 177 unrelated health control subjects. Clinical information of the study participants is shown in Table 1. Two hundred ONFH patients were consecutively enrolled at the Department of Orthopedics, the Second Hospital of Jilin University, (Changchun, China) from March, 2014 to June, 2015 in the study. Diagnostic and exclusion criteria, etiological classification, clinical stage, and unilateral or bilateral hip lesions for ONFH were followed as per reference (Song et al., 2017). Moreover, 177 unrelated health control subjects who were age and sex matched for the ONFH group were consecutively enrolled at the Health Examination Center of the Second Hospital of Jilin University, (Changchun, China) from October 2014 to December 2014.

Table 1.

Clinical Information of the Study Participants

Characteristics	Value	p
Age, median (range), year
Control	49 (26–83)	>0.05
ONFH	53 (21–84)
Sex (male/female)
Control	112/65	>0.05
ONFH	132/68
CHO (mM)
Control	3.978 ± 1.077	>0.05
ONFH	4.212 ± 1.648
TG (mM)
Control	1.161 ± 0.595	<0.01
ONFH	1.379 ± 0.616
HDL-c (mM)
Control	1.074 ± 0.321	<0.0001
ONFH	0.894 ± 0.497
LDL-c (mM)
Control	2.554 ± 0.823	<0.005
ONFH	2.835 ± 1.082
LDL-c/HDL-c
Control	2.494 ± 0.798	<0.0001
ONFH	3.395 ± 1.299
Clinical stage of ONFH (n, %)
Stage II	10 (5.6)
Stage III	54 (30.2)
Stage IV	115 (64.2)
Etiological classification of ONFH (n, %)
Alcohol-induced	71 (39.7)
Idiopathic	64 (34.0)
Steroid-induced	47 (26.3)
Unilateral or bilateral hip lesions of ONFH (n, %)
Unilateral	76 (42.5)
Bilateral	103 (57.5)

HDL-c, high-density lipoprotein cholesterol; LDL-c, low-density lipoprotein cholesterol; ONFH, osteonecrosis of the femoral head.

Health control subjects were also defined following references (Song et al., 2017). All of the 377 participants were Han Chinese from northeast China and gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the 2008 Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the Second Hospital of Jilin University, Changchun, China (2014-34). The serum TC, TG, high-density lipoprotein cholesterol (HDL-c), and low-density lipoprotein cholesterol (LDL-c) level of all participants were detected by automatic biochemistry analyzer (AU680; Beckman Coulter, Inc., Brea, CA).

Genomic DNA extraction, variant selection, and genotyping

Approximately 2 mL of venous blood was collected from all of the participants after a minimum of 10 h fasting. Genomic DNA was extracted following reference [12]. The database: http://gvs.gs.washington.edu/GVS138/was used to select variants of the genes according to reference (Song et al., 2017). The 10 variants from the five genes were selected based on LD analysis by HapMap (www.hapmap.org/index.html.en). The list of 10 variants is shown in Table 2.

Table 2.

Basic Information and Primers for PCR and Sequencing of 10 Variants in OPG, RANK, RANKL, TRAF 6, and NFATC 1 Genes

Gene	Chr	Variant ID	Allele (+/−) ^a	Positon	Function	PCR primer	Sequencing primer
OPG	8q24	rs2073617	C/T (+)	NM_002546.3:c.-223C>T	5′UTR	5′ACGTTGGATGTTTCCGCCCCAGCCCTGAAA3′	AAAGCCCTGAAAGCGTTAA
		rs2073617	C/T (+)	NM_002546.3:c.-223C>T	5′UTR	5′ACGTTGGATGTTTTTAGGAAGCTTGGGCGG3′	AAAGCCCTGAAAGCGTTAA
		rs2073618	C/G (+)	NM_002546.3:c.9C>G	Missense	5′ACGTTGGATGGACTTACCACGAGCGCGCA3′	GGGGAGAGCGCGCAGCACAGCAA
		rs2073618	C/G (+)	NM_002546.3:c.9C>G	Missense	5′ACGTTGGATGTCCAAGCCCCTGAGGTTTC3′	GGGGAGAGCGCGCAGCACAGCAA
RANK	18q22.1	rs884205	G/T (−)	NM_003839.3:c.2590A>C	3′UTR	5′ACGTTGGATGTCTTTCACTGGGTACCCAGG3′	TGGTACCCAGGCTGTATTTTAA
RANK	18q22.1	rs884205	G/T (−)	NM_003839.3:c.2590A>C	3′UTR	5′ACGTTGGATGATTCCCTGGACCACAGACTA3′	TGGTACCCAGGCTGTATTTTAA
RANKL	13q14	rs7984870	G/C (+)	NM_003701.3:c.-1958G>	Promoter	5′ACGTTGGATGAACAACTAAACCCACAATG3′	ATTTTACTATAGTGTCTTCCATT
		rs7984870	G/C (+)	NM_003701.3:c.-1958G>	Promoter	5′ACGTTGGATGTGAGGATCAAGATGTCAAGG3′	ATTTTACTATAGTGTCTTCCATT
		rs9525641	T/C (+)	NM_003701.3:c.-416T>C	Promoter	5′ACGTTGGATGCTGAGGTGCAATCCTTATCC3′	GGTGTCCTCTGCGTCTTC
		rs9525641	T/C (+)	NM_003701.3:c.-416T>C	Promoter	5′ACGTTGGATGTTCTGGAGGTCCAAGAGATG3′	GGTGTCCTCTGCGTCTTC
		rs1054016	G/T (+)	NM_003701.3:c.948G>T	Promoter	5′ACGTTGGATGCTTTTTTGCTGATAAACCAC3′	CACTAATAACTTTTCTATTTTGAAAAAA
		rs1054016	G/T (+)	NM_003701.3:c.948G>T	Promoter	5′ACGTTGGATGTGTTTTTCCACAAGTGCCGC3′	CACTAATAACTTTTCTATTTTGAAAAAA
NFATC1	18q23	rs754093	T/G (+)	NM_001278669.1:c.2251T>G	3′UTR	5′ACGTTGGATGAAGACCATACTACAGCCAGC3′	TCAGCGTGCCCAGCAAGGCCTACAC
		rs754093	T/G (+)	NM_001278669.1:c.2251T>G	3′UTR	5′ACGTTGGATGTGGCATCAGGGTGCTTCTCT3′	TCAGCGTGCCCAGCAAGGCCTACAC
		rs9518	T/C (+)	NM_001278669.1:c.1229T>C	Missense	5′ACGTTGGATGTGTACAGAGTGCCCAGCAAG3′	CCCCGCCGGCTTCCCGCCC
		rs9518	T/C (+)	NM_001278669.1:c.1229T>C	Missense	5′ACGTTGGATGCCTGCGTTGTGTGTTTTCAG3′	CCCCGCCGGCTTCCCGCCC
TRAF6	11p12	rs5030411	C/T (−)	NM_004620.3:c.-2178C>T	Promoter	5′ACGTTGGATGAACGCTTTGACCTAGGAATG3′	GGGAACTAGGAATGATATGCATCCTTCT
		rs5030411	C/T (−)	NM_004620.3:c.-2178C>T	Promoter	5′ACGTTGGATGTTCCCTTTGCAATCAGGCAG3′	GGGAACTAGGAATGATATGCATCCTTCT
		rs5030416	A/C (−)	NM_004620.3:c.-913A>C	Promoter	5′ACGTTGGATGTCTGCATATTTGGTTCATCG3′	TGGTTCATCGTCCACC
		rs5030416	A/C (−)	NM_004620.3:c.-913A>C	Promoter	5′ACGTTGGATGGCATTCTGGCTTCAGGATAC3′	TGGTTCATCGTCCACC

SNP alleles are reported in forward (+) or reverse (−) orientation with respect to the genome.

Chr, chromosome; NFATC1, nuclear factor of activated T cell cytoplasmic 1; OPG, osteoprotegerin; RANK, receptor activator of nuclear factor-kappa B; RANKL, receptor activator of nuclear factor-kappa B ligand; TRAF6, tumor necrosis factor receptor-associated factor 6; UTR, untranslated regions.

Primers for polymerase chain reaction and sequencing were designed by Sequenom Assay Design 3.1 software (Sequenom, San Diego, CA) following the manufacturer's protocols, shown in Table 2. Mass ARRAY^® platform (Sequenom Analyzer 4, Inc., San Diego, CA) was used to genotype the 10 variants of RANK, RANKL, OPG, tumor necrosis factor receptor-associated factor 6 (TRAF6), and nuclear factor of activated T cells cytoplasmic 1 (NFATC1) genes in 200 ONFH patients and 177 healthy controls. The genotyping success rates for the 10 variants were >95%, respectively.

Statistical analysis

Shesis software (http://analysis.bio-x.cn/SHEsisMain.htm) was used to analyze the Hardy–Weinberg equilibrium, genotype, and allele frequency between the ONFH and control groups. The association of paired gene–gene interaction of the 10 variants in RANKL/RANK/OPG pathway genes with the risk, clinical stages, and hip lesions of ONFH were estimated using binary logistic regression model by SPSS 21.0 software, respectively. Case and control individuals were defined as a “1” and “0,” clinical stage III and stage IV individuals of ONFH were defined as a “0” and “1,” unilateral and bilateral hip lesion individuals of ONFH were defined as a “0” and “1,” respectively. Every variant genotype was defined from 0 to 2 indicating the number of minor alleles in an individual. Each of the two variants that was multiplicative pairwise was an independent variable; case–control status, stage IV–stage III status, and bilateral–unilateral hip lesion status were as dependent variables for building the logistic regression model. p < 0.05 was considered as statistical significance. The regression models were adjusted, including age and sex.

Generalized multifactor dimensionality reduction (GMDR, GMDR Software Beta 0.9, www.ssg.uab.edu/gmdr/) method was used to examine all possible genetic models for multiplicative gene–gene interaction effects. Statistical significance was assessed using 1000-fold permutation testing. The value of Training or Testing Accuracy is between 0.50 and 1.00, indicating model prediction degree. Sign test (p) displays the number of testing accuracies greater than 0.5 and, in parentheses, the p-value computed using the nonparametric sign test. The crossvalidation (CV) consistency score was a measure of the number of times, in a particular crossvalidated run, a given attribute combination was selected as the best model.

Results

Genotypes and allele frequencies of the 10 variants and serum liquid levels between ONFH and control groups

Genotypes and allele frequencies of the 10 variants through Mass ARRAY platform in 377 samples are shown in Table 3. χ² test results revealed that the genotypes and allele frequencies of the 10 variants failed to show statistical significance between ONFH and control groups in spite of a marginal increase of NFATC1 rs9518 CC frequency in the ONFH group (p = 0.057). The TG and LDL-c levels of serum in the ONFH group were significantly increased, compared with the control group (p = 0.01, p = 0.005, respectively), but the serum HDL-c level of the ONFH group was dramatically lower than the control group (p = 0.0001), as shown in Table 1.

Table 3.

Genotypes and Allele Frequencies of 10 Variants in OPG, RANK, RANKL, TRAF 6, and NFATC 1 Genes Between ONFH Patients and Controls

Gene	Variant ID	Group	Genotype (frequency)			MAF	HWE ^a	p^b	Allele (frequency)		p^b	OR (95% CI)
OPG	rs2073617 (C/T)		CC	CT	TT				C	T
		Control	33 (0.186)	85 (0.480)	59 (0.333)	0.427	0.807	0.990	151 (0.427)	203 (0.573)	0.902	0.982
		Control	33 (0.186)	85 (0.480)	59 (0.333)	0.427	0.807		151 (0.427)	203 (0.573)		(0.735–1.312)
		Case	36 (0.181)	96 (0.482)	67 (0.337)	0.422	0.875		168 (0.422)	230 (0.578)		(0.735–1.312)
	rs2073618 (G/C)		CC	GC	GG				C	G
		Control	14 (0.079)	70 (0.395)	93 (0.525)	0.277	0.870	0.726	98 (0.277)	256 (0.723)	0.833	0.966
		Case	12 (0.060)	84 (0.420)	104 (0.520)	0.270	0.355	0.726	108 (0.270)	292 (0.730)	0.833	(0.701–1.332)
RANK	rs884205 (G/T)		GG	GT	TT				G	T
		Control	106 (0.599)	64 (0.362)	7 (0.040)	0.220	0.486	0.280	276 (0.780)	78 (0.220)	0.128	0.767
		Control	106 (0.599)	64 (0.362)	7 (0.040)	0.220	0.486		276 (0.780)	78 (0.220)		(0.545–1.080)
		Case	94 (0.516)	78 (0.429)	10 (0.055)	0.277	0.229		266 (0.731)	98 (0.269)		(0.545–1.080)
RANKL	rs1054016 (G/T)		GG	GT	TT				G	T
		Control	49 (0.277)	91 (0.514)	37 (0.209)	0.466	0.661	0.792	189 (0.534)	165 (0.466)	0.507	1.105
		Case	55 (0.307)	90 (0.503)	34 (0.190)	0.430	0.792	0.792	200 (0.559)	158 (0.441)	0.507	(0.823–1.485)
	rs7984870 (G/C)		CC	CG	GG				C	G
		Control	36 (0.206)	88 (0.503)	51 (0.291)	0.457	0.862	0.653	160 (0.457)	190 (0.543)	0.403	1.141
		Case	34 (0.227)	79 (0.527)	37 (0.247)	0.490	0.510	0.653	147 (0.490)	153 (0.510)	0.403	(0.838–1.554)
	rs9525641 (C/T)		CC	CT	TT				C	T
		Control	36 (0.205)	90 (0.511)	50 (0.284)	0.460	0.698	0.515	162 (0.460)	190 (0.540)	0.849	1.028
		Case	36 (0.182)	113 (0.571)	49 (0.247)	0.467	0.039	0.515	185 (0.467)	211 (0.533)	0.849	(0.771–1.371)
TRAF6	rs5030411 (C/T)		CC	CT	TT				C	T
		Control	12 (0.068)	72 (0.407)	93 (0.525)	0.271	0.699	0.721	96 (0.271)	258 (0.729)	0.588	1.092
		Control	12 (0.068)	72 (0.407)	93 (0.525)	0.271	0.699		96 (0.271)	258 (0.729)		(0.794–1.503)
		Case	18 (0.090)	79 (0.397)	102 (0.513)	0.289	0.632		115 (0.289)	283 (0.711)		(0.794–1.503)
	rs5030416 (A/C)		AA	AC	CC				A	C
		Control	126 (0.712)	48 (0.271)	3 (0.017)	0.153	0.515	0.748	300 (0.847)	54 (0.153)	0.847	1.040
		Case	146 (0.730)	49 (0.245)	5 (0.025)	0.148	0.715	0.748	341 (0.853)	59 (0.147)	0.847	(0.697–1.552)
NFATC1	rs754093 (G/T)		GG	GT	TT				G	T
		Control	26 (0.147)	82 (0.463)	69 (0.390)	0.379	0.838	0.691	134 (0.379)	220 (0.621)	0.453	0.892
		Case	27 (0.138)	84 (0.429)	85 (0.434)	0.350	0.396	0.691	138 (0.352)	254 (0.648)	0.453	(0.662–1.202)
	rs9518 (C/T)		CC	CT	TT				C	T
		Control	1 (0.006)	55 (0.229)	123 (0.695)	0.159	0.061	0.057	55 (0.155)	299 (0.845)	0.789	1.055
		Case	8 (0.040)	49 (0.245)	143 (0.715)	0.163	0.158	0.057	65 (0.163)	335 (0.845)	0.789	(0.713–1.560)

p-values of deviation from Hardy–Weinberg equilibrium between the ONFH group and control group.

χ² test.

CI, confidence interval; OR, odds ratio.

Association of paired gene–gene interactions of 10 variants in the OPG, RANK, RANKL, TRAF6, and NFATC1 genes with ONFH risk

Paired gene–gene interactions of the 10 variants of OPG, RANK, RANKL, TRAF6, and NFATC1 genes related to ONFH risk were completed by using Logistic regression method. A total of six interactional variants were found to associate significantly with ONFH risk, including the paired interaction between OPG rs2073617 and NFATC1 rs754093, which was significantly associated with increased ONFH risk, and odds ratio (OR) (95% confidence interval [CI]) was 2.367 (1.154–4.855), p = 0.019, whereas the paired interactions between OPG rs2073618 and NFATC1 rs754093 as well as between OPG rs2073617 and RANKL rs1054016 were significantly associated with decreased ONFH risk, and OR (95% CI) was 0.268 (0.099–0.711), p = 0.008; 0.282 (0.085–0.939), p = 0.039, respectively. In addition, the paired interactions between RANK rs884205 and NFATC1 rs754093, between RANK rs884205 and NFATC1 rs9518, and between RANKL rs1054016 and NFATC1 rs754093 were marginal changes, their OR (95% CI) were 0.659 (0.427–1.017), p = 0.060; 1.735 (0.945–3.184), p = 0.076; and 0.805 (0.632–1.026), p = 0.080, respectively, as shown in Table 4 and Figure 1.

FIG. 1.

Association of paired variants interactions in the 10 variants of RANK/RANKL/OPG pathway genes with ONFH risk. 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 genes of variants located. The innermost circle is the name of variants. The lines connecting two variants represent paired gene–gene interactions. The thinner gray lines represent interactions with no statistical significance while thicker red lines represent interactions with statistical significance. ONFH, osteonecrosis of the femoral head; OPG, osteoprotegerin; RANKL, receptor activator of nuclear factor-kappa B ligand; RANK, receptor activator of nuclear factor-kappa B; TRAF6, tumor necrosis factor receptor-associated factor6; NFATC1, nuclear factor of activated T cell cytoplasmic 1. Color images are available online.

Table 4.

Association of Paired Gene–Gene Interactions of 10 Variants in OPG, RANK, RANKL, TRAF 6, and NFATC 1 Genes with ONFH Risk

Gene	Variants	OPG		RANK	RANKL			TRAF6		NFATC1
Gene	Variants	rs2073617	rs2073618	rs884205	rs7984870	rs9525641	rs1054016	rs5030411	rs5030416	rs754093	rs9518
OPG	rs2073617	—	0.946 (0.622–1.440)	0.878 (0.411–1.874)	0.97 (0.771–1.229)	0.97 (0.784–1.204)	0.282 (0.085–0.939)	1.617 (0.567–4.614)	0.780 (0.192–3.173)	2.367 (1.154–4.855)	1.135 (0.496–2.599)
	rs2073617	—	0.796	0.736	0.820	0.791	0.039 ^a	0.369	0.728	0.019 ^a	0.764
	rs2073618	0.946 (0.622–1.440)	—	1.451 (0.832–2.357)	1.132 (0.854–1.501)	1.124 (0.862–1.466)	0.446 (0.140–1.423)	1.240 (0.695–2.213)	0.728 (0.344–1.584)	0.268 (0.099–0.711)	1.163 (0.668–2.023)
		0.946 (0.622–1.440)		0.133	0.390	0.389	0.172	0.172	0.435	0.008 ^b	0.593
		0.796		0.133	0.390	0.389	0.172	0.172	0.435	0.008 ^b	0.593
RANK	rs884205	0.878 (0.411–1.874)	1.451 (0.832–2.357)	—	1.270 (0.935–1.725)	1.220 (0.916–1.625)	1.054 (0.369–3.011)	1.482 (0.830–2.647)	0.631 (0.280–1.418)	0.659 (0.427–1.017)	1.735 (0.945–3.184)
RANK	rs884205	0.736	0.133	—	0.126	0.175	0.922	0.183	0.265	0.060	0.076
RANKL	rs7984870	0.97 (0.771–1.229)	1.132 (0.854–1.501)	1.270 (0.935–1.725)	—	1.558 (0.753–3.225),	0.605 (0.285–1.285),	1.284 (0.835–1.978),	0.762 (0.437–1.327)	0.875 (0.652–1.176)	1.124 (0.759–1.664)
		0.820	0.390	0.126		1.558 (0.753–3.225),	0.605 (0.285–1.285),	1.284 (0.835–1.978),	0.762 (0.437–1.327)	0.875 (0.652–1.176)	1.124 (0.759–1.664)
		0.820	0.390	0.126		0.232	0.191	0.257	0.357	0.377	0.561
	rs9525641	0.97 (0.784–1.204)	1.124 (0.862–1.466)	1.220 (0.916–1.625)	1.558 (0.753–3.225)	—	0.978 (0.786–1.216),	1.093 (0.723–1.652)	0.861 (0.507–1.463)	0.839 (0.635–1.108)	1.101 (0.770–1.573)
	rs9525641	0.791	0.389	0.175	0.232	—	0.839	0.674	0.581	0.251	0.598
	rs1054016	0.282 (0.085–0.939)	0.446 (0.140–1.423)	1.054 (0.369–3.011)	0.605 (0.285–1.285),	0.978 (0.786–1.216),	—	1.085 (0.739–1.593)	0.825 (0.494–1.381)	0.805 (0.632–1.026)	1.102 (0.776–1.564)
		0.282 (0.085–0.939)	0.446 (0.140–1.423)	1.054 (0.369–3.011)	0.605 (0.285–1.285),	0.978 (0.786–1.216),		0.676	0.465	0.080	0.587
		0.039 ^a	0.172	0.922	0.191	0.839		0.676	0.465	0.080	0.587
TRAF6	rs5030411	1.617 (0.567–4.614)	1.240 (0.695–2.213)	1.482 (0.830–2.647)	1.284 (0.835–1.978)	1.093 (0.723–1.652)	1.085 (0.739–1.593)	—	1.040 (0.719–1.504)	1.183 (0.781–1.791)	0.829 (0.357–1.928)
	rs5030411	0.369	0.172	0.183	0.257	0.674	0.676	—	0.835	0.428	0.664
	rs5030416	0.780 (0.192–3.173)	0.728 (0.344–1.584)	0.631 (0.280–1.418)	0.762 (0.437–1.327)	0.861 (0.507–1.463)	0.825 (0.494–1.381)	1.040 (0.719–1.504)	—	0.739 (0.369–1.479)	1.592 (0.538–4.711)
	rs5030416	0.728	0.435	0.265	0.357	0.581	0.465	0.835	—	0.392	0.401
NFATC1	rs754093	2.367 (1.154–4.855)	0.268 (0.099–0.711)	0.659 (0.427–1.017)	0.875 (0.652–1.176)	0.839 (0.635–1.108)	0.805 (0.632–1.026)	1.183 (0.781–1.791)	0.739 (0.369–1.479)	—	1.009 (0.668–1.525)
	rs754093	0.019 ^a	0.008 ^b	0.060	0.377	0.251	0.080	0.428	0.392	—	0.965
	rs9518	1.135 (0.496–2.599)	1.163 (0.668–2.023)	1.735 (0.945–3.184)	1.124 (0.759–1.664)	1.101 (0.770–1.573)	1.102 (0.776–1.564)	0.829 (0.357–1.928)	1.592 (0.538–4.711)	1.009 (0.668–1.525)	—
	rs9518	0.764	0.593	0.076	0.561	0.598	0.587	0.664	0.401	0.965	—

Bold values mean statistical significance or marginal changes.

p < 0.05, ^b p < 0.01 represent the interaction between two variants, respectively.

Associations of multidimensional interactions among 10 variants in OPG, RANK, RANKL, TRAF6, and NFATC1 genes with ONFH risk

The number of dimensions was set from 1 to 10. The five best models were identified associated with ONFH risk, including a two dimension containing RANKL rs798487-TRAF6 rs5030416, with testing accuracy of 0.6245, CV consistency of 8/10, and p = 0.001; a three dimension containing RANKrs884205-RANKLrs7984870-NFATC1 rs754193, with testing accuracy of 0.6049, CV consistency of 7/10, and p = 0.01; a five dimension containing OPG rs2073617-RANK rs884205-RANKL rs7984870-TRAF6 rs5030411-NFATC1 rs7540930, with testing accuracy of 0.5795, CV consistency of 9/10, and p = 0.01; a seven dimension containing OPG rs2073617-OPGrs2073618-RANKrs884205-RANKLrs7984870-TRAF6rs5030411-NFATC1rs754093-NFATC1rs9518, with testing accuracy of 0.8151, CV consistency of 10/10, and p = 0.01; and an eight dimension containing OPGrs2073617-OPGrs2073618-RANKrs884205-RANKLrs7984870-TRAF6rs5030411-TRAFrs5030416-NFAC1rs754093-NFAC1-rs9518, with testing accuracy of 0.6420, CV consistency of 10/10, and p = 0.01. Thus, we concluded that the five multiple dimension models were significantly associated with ONFH risk, as shown in Table 5.

Table 5.

Association of Multidimensional Gene–Gene Interactions of 10 Variants in OPG, RANK, RANKL, TRAF 6, and NFATC 1 Genes with ONFH Risk

Model	Training balanced accuracy	Testing balanced accuracy	Sign test (p)	CV consistency
rs7984870	0.6194	0.6194	10 (0.001)	10/10
rs7984870-rs1054016	0.6513	0.6245	10 (0.001) ^b	8/10
rs884205-rs7984870-rs754193	0.6822	0.6049	9 (0.01)^a	7/10
rs2073617-rs884205-rs7984870-rs754193	0.7351	0.5692	8 (0.0547)	7/10
rs2073617-rs884205-rs7984870-rs5030411-rs754193	0.7945	0.5795	9 (0.01) ^a	9/10
rs2073617-rs2073618-rs884205-rs7984870-rs5030411-rs754193	0.8544	0.6196	8 (0.0547)	10/10
rs2073617-rs2073618-rs884205-rs7984870-rs5030411-rs754193-rs9518	0.8950	0.6151	9 (0.01) ^a	10/10
rs2073617-rs2073618-rs884205-rs7984870-rs5030411-rs5030416-rs754193-rs9518	0.9137	0.6420	9 (0.01) ^a	10/10
rs2073617-rs2073618-rs884205-rs7984870-rs1054016-rs5030411-rs5030416-rs754193-rs9518	0.9182	0.6066	8 (0.0547)	10/10
rs2073617-rs2073618-rs884205-rs7984870-rs9525641-rs1054016-rs5030411-rs5030416-rs754193-rs9518	0.9182	0.6055	7 (0.1769)	10/10

Bold values mean statistical significance.

RANK rs884205, RANKL rs1054016, RANKL rs9525641, RANKL rs7984870, OPG rs2073617, OPG rs2073618, TRAF6 rs5030411, TRAF6 rs5030416, NFATC1 rs754093, NFATC1 rs9518.

p = 0.01, ^b p = 0.001 represent the interactions among variants, respectively.

CV, crossvalidation.

Association of paired gene–gene interactions of 10 variants in the OPG, RANK, RANKL, TRAF6, and NFATC1 genes with the clinical stages and unilateral or bilateral hip lesions of ONFH

The association of paired gene–gene interactions of the 10 variants with stage III and stage IV of ONFH as well as with unilateral and bilateral hip lesions of ONFH were also completed by using logistic regression method (SPSS21.0 software), respectively. The paired interactions between rsRANKL rs9525641 and rs7984870, between RANKL rs7984870 and rs1054016, and between OPG rs2073618 and TRAF6 rs5030411 revealed marginal increased stage IV risk of ONFH, and their OR (95% CI) was 1.375 (0.954–1.982), p = 0.088; 1.390 (0.983–1.967), p = 0.063; 1.375 (0.954–1.982), p = 0.088, respectively, whereas the interaction between OPG rs2073618 and TRAF6 rs5030411 revealed marginal decreased stage IV risk of ONFH, 0.6 (0.403–1.059), p = 0.084, as shown in Table 6. A total of four interactional variants were found to associate significantly with the increased risk of bilateral hip lesions of ONFH, including the paired interaction between OPG rs2073617 and TRAF6 rs5030411, 1.544 (1.011–2.348), p = 0.044; and between RANKrs884205 and TRAF6 rs5030411, 1.772 (1.012–3.104), p = 0.045, respectively. Also, the paired interactions between OPG rs2073617 and TRAF6 rs5030416, between RANK rs884205 and TRAF6 rs5030416, between RANKLrs9525641 and TRAF6 rs5030416, between RANKL rs1054016 and TRAF6 rs5030416, between RANKrs884205 and NFATC1 rs754193, and between RANK rs884205 and NFATC1 rs9518 also revealed marginal increased risk of bilateral hip lesions in ONFH, their OR (95% CI) was 1.771 (0.978–3.209), p = 0.059; 2.592 (0.956–7.026), p = 0.061; 1.608 (0.954–2.710), p = 0.074; 1.828 (0.980–3.142), p = 0.058; 1.612 (0.973–2.670), p = 0.064; and 1.951 (0.950–4.009), p = 0.069, respectively, as shown in Table 7.

Table 6.

The Association of the Paired Gene Interactions with Clinical Stages of ONFH

		OPG		RANK	RANKL			TRAF6		NFATC1
		rs2073617	rs2073618	rs884205	rs7984870	rs9525641	rs1054016	rs5030411	rs5030416	rs754193	rs9518
OPG	rs2073617	0.726 (0.456–1.156)	0.873 (0.638–1.195)	0.898 (0.571–1413)	1.107 (0.722–1.698)	0.891 (0.639–1.243)	0.889 (0.609–1.276)	0.830 (0.554–1.243)	0.747 (0.439–1.271)	0.961 (0.686–1.341)	0.792 (0.505–1.243)
	rs2073617	0.177	0.396	0643	0.642	0.497	0.523	0.365	0.287	0.817	0.311
	rs2073618	0.873 (0.638–1.195)	0.716 (0.444–1.303)	0.845 (0.470–1.520)	1.033 (0.631–1.692)	0.910 (0.607–1.363)	0.860 (0.571–1.296)	0.654 (0.403–1.059)	0.580 (0.304–1.109)	1.122 (0.664–1.896)	0.739 (0.402–1.360)
	rs2073618	0.396	0.319	0.574	0.897	0.646	0.472	0.084	0.099	0.667	0.331
RANK	rs884205	0.898 (0.571–1413)	0.845 (0.470–1.520)	1.037 (0.593–1.814)	1.122 (0.653–1.930)	1.100 (0.699–1.733)	1.144 (0.728–1.797)	0.722 (0.442–1.180)	0.681 (0.357–1.298)	1.035 (0.649–1.649)	0.916 (0.517–1.625)
RANK	rs884205	0643	0.574	0.898	0.677	0.680	0.560	0.193	0.243	0.885	0.764
RANKL	rs7984870	1.107 (0.722–1.698)	1.033 (0.631–1.692)	1.122 (0.653–1.930)	1.851 (1.015–3.378)	1.390 (0.983–1.967)	1.375 (0.954–1.982)	0.926 (0.592–1.449)	0.950 (0.547–1.650)	1.424 (0.892–2.275)	1.070 (0.628–1.823)
	rs7984870	0.642	0.897	0.677	0.045	0.063	0.088	0.737	0.855	0.139	0.804
	rs9525641	0.891 (0.639–1.243)	0.910 (0.607–1.363)	1.100 (0.699–1.733)	1.390 (0.983–1.967)	1.385 (0.827–2.320)	1.188 (0.891–1.582)	0.823 (0.566–1.196)	0.773 (0.459–1.208)	1.269 (0.856–1.881)	0.800 (0.537–1.194)
	rs9525641	0.497	0.646	0.680	0.063	0.215	0.241	0.307	0.259	0.235	0.275
	rs1054016	0.889 (0.609–1.276)	0.860 (0.571–1.296)	1.144 (0.728–1.797)	1.375 (0.954–1.982)	1.188 (0.891–1.582)	1.304 (0.780–2.181)	0.859 (0.580–1.273)	0.816 (0.488–1.362)	1.188 (0.781–1.807)	0.784 (0.502–1.225)
	rs1054016	0.523	0.472	0.560	0.088	0.241	0.312	0.449	0.436	0.421	0.285
TRAF6	rs5030411	0.830 (0.554–1.243)	0.654 (0.403–1.059)	0.722 (0.442–1.180)	0.926 (0.592–1.449)	0.823 (0.566–1.196)	0.859 (0.580–1.273)	0.769 (0.470–1.257)	0.791 (0.541–1.155)	0.928 (0.617–1.395)	0.787 (0.456–1.358)
	rs5030411	0.365	0.084	0.193	0.737	0.307	0.449	0.294	0.225	0.719	0.389
	rs5030416	0.747 (0.439–1.271)	0.580 (0.304–1.109)	0.681 (0.357–1.298)	0.950 (0.547–1.650)	0.773 (0.459–1.208)	0.816 (0.488–1.362)	0.791 (0.541–1.155)	0.792 (0.429–1.460)	0.820 (0.462–1.384)	0.607 (0.310–1.191)
	rs5030416	0.287	0.099	0.243	0.855	0.259	0.436	0.225	0.454	0.425	0.147
NFATC1	rs754193	0.961 (0.686–1.341)	1.122 (0.664–1.896)	1.035 (0.649–1.649)	1.424 (0.892–2.275)	1.269 (0.856–1.881)	1.188 (0.781–1.807)	0.928 (0.617–1.395)	0.820 (0.462–1.384)	1.143 (0.708–1.844)	1.140 (0.698–1.862)
	rs754193	0.817	0.667	0.885	0.139	0.235	0.421	0.719	0.425	0.585	0.602
	rs9518	0.792 (0.505–1.243)	0.739 (0.402–1.360)	0.916 (0.517–1.625)	1.070 (0.628–1.823)	0.800 (0.537–1.194)	0.784 (0.502–1.225)	0.787 (0.456–1.358)	0.607 (0.310–1.191)	1.140 (0.698–1.862)	1.047 (0.589–1.861)
	rs9518	0.311	0.331	0.764	0.804	0.275	0.285	0.389	0.147	0.602	0.877

Bold values mean statistical significance or marginal changes.

The data are represented as OR (95% CI), p values (stage IV vs. stage III) by logistic regression. An interaction term was considered significant if p < 0.05. The regression models were adjusted, including age and sex. Clinical stages were only involved in stage III and stage IV.

Table 7.

The Association of the Paired Gene Interactions with Unilateral and Bilateral Hip Lesions of ONFH

		OPG		RANK	RANKL			TRAF6		NFATC1
		rs2073617	rs2073618	rs884205	rs7984870	rs9525641	rs1054016	rs5030411	rs5030416	rs754193	rs9518
OPG	rs2073617	1.230 (0.805–1.878)	1.054 (0.785–1.415)	1.198 (0.766–1.875)	0.879 (0.605–1.279)	0.932 (0.681–1.275)	0.969 (0.688–1.364)	1.544 (1.011–2.348)	1.771 (0.978–3.209)	1.213 (0.870–1.690)	1.332 (0.831–2.134)
	rs2073617	0.338	0.725	0.429	0.501	0.659	0.856	0.044	0.059	0.265	0.234
	rs2073618	1.054 (0.785–1.415)	1.021 (0.623–1.674)	1.206 (0.676–2.153)	0.897 (0.579–1.389)	0.836 (0.572–1.222)	0.855 (0.580–1.260)	1.143 (0.711–1.839)	1.204 (0.629–2.266)	1.221 (0.753–1.980)	1.497 (0.769–2.913)
	rs2073618	0.725	0.933	0.526	0.626	0.355	0.427	0.581	0.566	0.418	0.235
RANK	rs884205	1.198 (0.766–1.875)	1.206 (0.676–2.153)	1.378 (0.815–2.329)	1.387 (0.847–2.269)	1.187 (0.779–1.809)	1.247 (0.821–1.892)	1.772 (1.012–3.104)	2.592 (0.956–7.026)	1.612 (0.973–2.670)	1.951 (0.950–4.009)
RANK	rs884205	0.429	0.526	0.231	0.193	0.425	0.300	0.045	0.061	0.064	0.069
RANKL	rs7984870	0.879 (0.605–1.279)	0.897 (0.579–1.389)	1.387 (0.847–2.269)	1.020 (0.611–1.703)	1.077 (0.828–1.401)	1.153 (0.864–1.537)	1.160 (0.760–1.770)	1.327 (0.763–2.308)	1.151 (0.802–1.652)	1.495 (0.860–2.599)
	rs7984870	0.501	0.626	0.193	0.939	0.519	0.334	0.492	0.317	0.445	0.154
	rs9525641	0.932 (0.681–1.275)	0.836 (0.572–1.222)	1.187 (0.779–1.809)	1.077 (0.828–1.401)	1.032 (0.647–1.648)	1.121 (0.870–1.444)	1.209 (0.833–1.754)	1.608 (0.954–2.710)	1.151 (0.824–1.609)	1.274 (0.829–1.959)
	rs9525641	0.659	0.355	0.425	0.519	0.894	0.377	0.317	0.074	0.408	0.269
	rs1054016	0.969 (0.688–1.364)	0.855 (0.580–1.260)	1.247 (0.821–1.892)	1.153 (0.864–1.537)	1.121 (0.870–1.444)	1.181 (0.739–1.886)	1.294 (0.864–1.938)	1.828 (0.980–3.142)	1.336 (0.909–1.966)	1.367 (0.825–2.265)
	rs1054016	0.856	0.427	0.300	0.334	0.377	0.487	0.211	0.058	0.141	0.225
TRAF6	rs5030411	1.544 (1.011–2.348)	1.143 (0.711–1.839)	1.772 (1.017–3.104)	1.160 (0.760–1.770)	1.209 (0.833–1.754)	1.294 (0.864–1.938)	1.424 (0.896–2.265)	1.350 (0.895–2.036)	1.184 (0.792–1.769)	1.763 (0.863–3.600)
	rs5030411	0.044	0.581	0.045	0.492	0.317	0.211	0.135	0.153	0.411	0.120
	rs5030416	1.771 (0.978–3.209)	1.204 (0.629–2.266)	2.592 (0.956–7.026)	1.327 (0.763–2.308)	1.608 (0.954–2.710)	1.828 (0.980–3.142)	1.350 (0.895–2.036)	1.757 (0.942–3.278)	1.531 (0.820–2.860)	1.910 (0.744–4.900)
	rs5030416	0.059	0.566	0.061	0.317	0.074	0.058	0.153	0.076	0.181	0.178
NFATC1	rs754193	1.213 (0.870–1.690)	1.221 (0.753–1.980)	1.612 (0.973–2.670)	1.151 (0.802–1.652)	1.151 (0.824–1.609)	1.336 (0.909–1.966)	1.184 (0.792–1.769)	1.531 (0.820–2.860)	1.120 (0.723–1.733)	1.140 (0.698–1.862)
	rs754193	0.265	0.418	0.064	0.445	0.408	0.141	0.411	0.181	0.612	0.602
	rs9518	1.332 (0.831–2.134)	1.497 (0.769–2.913)	1.951 (0.950–4.009)	1.495 (0.860–2.599)	1.274 (0.829–1.959)	1.367 (0.825–2.265)	1.763 (0.863–3.600)	1.910 (0.744–4.900)	1.140 (0.698–1.862)	1.493 (0.859–2.594)
	rs9518	0.234	0.235	0.069	0.154	0.269	0.225	0.120	0.178	0.602	0.155

Bold values mean statistical significance or marginal changes.

The data are represented as OR (95% CI), p values (bilateral vs. unilateral hips) by logistic regression. An interaction term was considered significant if p < 0.05. The regression models were adjusted, including age and sex. Hip lesion includes unilateral and bilateral group.

Discussion

Epidemiological evidences have clearly revealed that the two major environmental risk factors, massive steroid treatment and heavy alcohol drinking, play significant roles in the development of ONFH. Therefore, ONFH is further classified into three subgroups, steroid-associated, alcohol-associated, and idiopathic ONFH. Lipid metabolism disorder has been recognized as the key pathogenesis of ONFH in that the typical pathological feature of ONFH shows the accumulation of a large amount of adipose tissue in the damaged bone marrow cavity. Thus, we also analyzed the serum lipid levels between the ONFH and control groups in this study. As we expected, the results showed that serum TG and LDL-c levels of ONFH group showed significant increase (p = 0.01, p = 0.005), respectively but the HDL-c level showed statistical decrease (p = 0.0001) compared with the control group. Until now, the relationship between RANK/RANKL/OPG pathway and lipid metabolism disorder has not been reported. However, a study revealed that the ratio of OPG/RANKL secretion in primary human preosteoblasts increased ninefold of both mRNA and protein when stimulated with adipocyte-secreted factors (Kühn et al., 2012), indicating the possible link between OPG/RANKL pathway and lipid metabolism disorder.

The increasing evidences have also shown that multiple genes are closely associated with the development of ONFH. Genome-wide association study (GWAS) represents a landmark in the genetic study of human complex diseases, including ONFH. GWAS results from steroid-associated osteonecrosis in children with acute lymphoblastic leukemia reported significant loci around glutamate ionotropic receptor NMDA type subunit 3A (GRIN3A), bone morphogenic protein7 (BMP7), long intergenic nonprotein-coding RNA 251 (LINC00251), and PROX1 antisense RNA 1 (PROX1-AS1, Karol et al., 2016). Another GWAS for ONFH in Japanese population with 1602 cases and 60,000 controls identified a 20q12 locus with genome-wide significance and found LINC01370 as a candidate susceptibility gene for ONFH (Sakamoto et al., 2017).

Considering that multiple genetic loci of moderate effect fails to reach genome-wide significance due to the limited power in most genetic studies, we focused on the association of gene–gene interactions among 10 variants in OPG/RANKL/RANK pathway genes with the risk and clinical phenotypes of ONFH in our present study. The results revealed that the genotypes and allele frequencies of the 10 variants in the RANKL, RANK, OPG, TRAF6, and NFATC1 genes between ONFH and control groups showed no statistical significance. However, these variants identified “negative” were found to correlate significantly with ONFH risk through the analysis of paired variant interactions. That is, the paired interaction between OPG rs2073617 and NFATC1 rs754093 was significantly associated with increased ONFH risk, p < 0.019, but the paired interactions between OPG rs2073618 and NFATC1 rs7540930 as well as between OPG rs2073617 and RANKL rs1054016 showed a significant association with decreased ONFH risk, p < 0.008, p < 0.039, respectively.

Moreover, the paired interactions of the 10 variants were also involved in the clinical phenotypes, particularly unilateral or bilateral hip lesions of ONFH. A total of the paired four interactional variants were found to associate significantly with the increased risk of bilateral hip lesions in ONFH, including the paired interaction between OPG rs2073617 and TRAF6 rs5030411 as well as between RANKrs884205 and TRAF6 rs503041, p = 0.044, p = 0.045, respectively. These variants previously failed to show the statistical association with ONFH risk (Li et al., 2016b) because they were only focused on a single variant analysis, which overlooked possible interactions between intervariants. For complex and multiple factor diseases, a single variant did not have a notable association with the development of a complex disease, such as ONFH. Gene–gene interactions take genetic context into account, and on the basis of gene–gene interaction analysis, the missing heritability, partly attributed to the low power to detect genes, was found to have their genetic contributions to ONFH.

Furthermore, the analysis of multidimensional variant interactions of 10 variants in OPG, RANK, RANKL, TRAF6, and NFATC1 genes showed that the five best models identified associated with ONFH risk, including a two-dimensional model of RANKL rs7984870-TRAF6 rs5030416 (p = 0.001), a three-dimensional model of RANK rs884205-RANKL rs7984870-NFATC1 rs7541938 (p = 0.01), a five-dimensional model of OPG rs2073617-RANKrs884205-RANKLrs7984870-TRAF6rs5030411-NFATC1rs7540930 (p = 0.01), a seven-dimensional model of OPG rs2073617-OPGrs2073618-RANKrs884205-RANKLrs7984870-TRAF6rs5030411-NFATC1rs754093-NFATC1rs9518 (p = 0.01), and an eight-dimensional model of OPGrs2073617-OPGrs2073618-RANKrs884205-RANKLrs7984870-TRAF6rs5030411-TRAF6rs5030416-NFAC1rs754093-NFAC1rs9518 (p = 0.01). Gene–gene interaction acts as an indispensable approach of studying how variants influence phenotype. It has been used to discover variants that were associated with platinum-based chemotherapeutic response and toxicity and the analysis strategy of gene–gene or gene–environment interaction may hold potential value to predict chemotherapeutic response and toxicity (Cui et al., 2017).

The effects of gene–gene interactions in RANK/RANKL/OPG pathway on ONFH risk and its clinical phenotypes could be partly explained by the specific functions of these variants. The rs2073617 and rs2073618 are located in 5 UTR and missense of TNFRSF11B on chromosome 8q23–24, respectively. The rs2073617 is only a few base pairs away from a TGF-response area and only 129 bp upstream of the TATA box (Brandstrom et al., 1999). OPG, protein coded by TNFRSF11B, is a member of the TNF-receptor superfamily. Therefore, it has been hypothesized that factors affecting OPGgene regulation may be major genetic factors influencing bone mass and increasing the risk of fractures, osteoporosis (OP), and osteoarthritis (OA, Sun et al., 2014). Paget's disease of bone, known as juvenile ONFH, was early proved to be associated with rs2073617 by genotyping data from a UK population (Beyens et al., 2007).

Evidence suggests also that rs2073617 polymorphism was associated with an increased OP risk in the Chinese population (Li et al., 2017). The RANK/RANKL/OPG pathway plays a central role in the pathogenesis of bone erosions in RA and a result from a meta-analysis involving three French cohorts demonstrated that rs2073618 was significantly associated with erosions in RA (Gu et al., 2015; Ruyssen-Witrand et al., 2016). In present study, the interactions of rs2073617 and rs2073618 with other variants through the paired and multiple variant interactions manner were shown to affect the risk and clinical phenotype of ONFH, which further indicate the significant effects of gene–gene interactions in RANK/RANKL/OPG pathway on the development of ONFH.

The rs754093 is located in missense of NFATC1 gene on chromosome18q23. NFATC1 is a component of the nuclear factor of activated T cell DNA-binding transcription complex. The gene encoding NFATc1 in man consists of 11 exons and span ∼134 kb DNA. NFATc1 expression was notably downregulated by the neutralization of macrophage migration inhibitory factor (MIF), indicating a likely key downstream target of MIF (Gu et al., 2015). A two-staged genotyping study identified the associations between variants in 11 genes (including NFATC1) and volume bone mineral density (vBMD, Yerges et al., 2009). Our previous results also revealed that the minor homozygote (CC) and the heterozygosis (TC) of NFATC1rs9518 were significantly associated with ONFH risk and the hip lesion development of ONFH (Song et al., 2017).

The rs754093 and rs9518 of NFATC1 gene in this study were involved in a seven- and eight-dimensional interaction model associated with ONFH risk, respectively, suggesting that the two variants, together with other genes of the pathway, exert genetic effect on ONFH risk. The rs1054016 and rs7984870 are located in the promoter of RANKL gene on chromosome 13q14. The receptor activator of RANK pathway is involved in bone metastasis of breast cancer. The rs1054016 was genotyped and analyzed with regard to bone metastasis-free survival (BMFS), disease-free survival, and overall survival for a retrospective cohort of 1251 patients and the results showed that rs1054016 seems to influence BMFS, in that patients with two minor alleles had a more favorable prognosis than patients with at least one common allele.

The effect of rs1054016 adds to the evidence that the RANK pathway plays a role in breast cancer pathogenesis and progression with respect to BMFS (Hein et al., 2014). The meta-analysis results also demonstrated that three variants, rs7984870, rs7325635, and rs1054016, and one haplotype of RANKL were significantly associated with anticitrullinated peptide antibodies presence in RA (Ruyssen-Witrand et al., 2016). Remarkably, the five best models of GMDR identified and associated significantly with ONFH risk in the study did include the rs7984870 of RANKL gene, suggesting the significant effect of the variant on ONFH risk through the interactions among multiple variants.

The rs5030416 and rs5030411 are located in promoter of TRAF6 gene on chromosome 11p12. It has been reported that the two variants with their CC haplotype significantly associates with susceptibility to ischemic stroke (IS) in a total of 816 IS cases and 816 controls in southern Chinese Han population, which suggest that TRAF6 gene polymorphisms may be involved in the pathogenesis of IS (Su et al., 2015). The rs884205 is located in 3′UTR of RANK gene on 18q22.1. Rank is a member of the TNF-receptor superfamily and it interacts with various TRAF family proteins and induces the activation of NF-κB and MAPK8/JNK. RA has been thought as the inflammation destruction of synovial and joint tissues.

A report showed that the expression of RANK on CD14⁺ monocyte in RA patients was significantly increased and correlated with osteoclastogenesis, which indicated that high expressed RANK on monocytes may be novel targets for the regulation of bone resorption in RA and OP (Nanke et al., 2016). A recent GWAS in Chinese Han subjects identified that four BMD-loci, including RANK rs884205, were significantly associated with hip osteoporotic fractures (Guo et al., 2012). More significantly, four of the five best models of GMDR in the study that did also involve in RANK rs884205, revealed its potential key roles in the development of ONFH through multiple variant interactions in the pathway.

Although a large number of susceptibility loci have been identified by GWAS, the disease variance explained by these loci was still limited. At least part of the missing heritability may reside in those genetic loci of moderate effect, which failed to reach genome-wide significance due to the limited power in most genetic studies (Li et al., 2016a). However, when the genes located in the same pathway, whose gene products work cooperatively to complete certain molecular and cellular functions, were jointly assessed, the statistical power of these genes can be significantly improved. Furthermore, the genetic molecular etiology of complex diseases, such as ONFH, is more likely correlated with the results of multiple gene interactions residing in the same pathway. Our present results, which were previously proved to be of no statistical significance by means of single variant genotype, showed to significantly associate with ONFH risk and its clinical phenotype through the analysis of both the paired and multiple gene interactions. By collapsing variant statistics to a gene level and further pathway gene level, the missing genetic contribution in OPG/RANK/RANKL pathway genes can be searched through gene–gene interactions.

In this study, several limitations need to be considered. First, the 377 ONFH case–control attributed to smaller samples that may limit our statistical power to detect the interactions of small differences between ONFH and control groups. Second, the lack of gene expression and their function analysis limited the interaction information between genes–functions. Third, the analysis without the markers of bone turnover related to ONFH decreased also the interaction information between genotype–phenotype. Fourth, multiple other than one approaches may be useful to improve the analysis of gene interactions to provide fully supporting evidence. Therefore, the roles of gene interactions in RANK/RANKL/OPG pathway in the development of ONFH remain to be further explored on the basis of gene expression and function, markers of bone turnover, and united gene interaction analysis approaches in a large cohort.

In conclusion, we explored the roles of the 10 variant interactions in the RANKL, RANK, OPG, TRAF6, and NFATC1 genes in the development of ONFH in 377 ONFH cases–controls. Our results were first showed that the paired gene interactions of six interactional variants as well as the five best models of GMDR analysis in the genes of RANK/RANKL/OPG pathway significantly associated with the risk and clinical phenotype of ONFH. These results suggest that the variants of RANK/RANKL/OPG pathway genes affect the development of ONFH in gene interaction manner through the interaction of paired and multiple variants.

Footnotes

Acknowledgments

This research was in part funded by the National Natural Science Foundation of China (Grant No. 81702195); by the Department of Science and Technology of Jilin Province, China (Grant No. 20180520125JH, 20190802001ZG); by the Education Department of Jilin Province (Grant No. JJKH20180103KJ); by the Project of Bethune Youth Foundation of Jilin University, China (Grant No. 2015409); by the Development and Reform Commission of Jilin Province, China (Grant No. 2014G073); and by the project of application demonstration center of precision medicine for molecular diagnosis in Jilin Province (2016–2018, NDRC).

Disclosure Statement

No competing financial interests exist.

References

Adesina

, Brunson

, Keegan

T.H.M.

, and Wun

(2017). Osteonecrosis of the femoral head in sickle cell disease: prevalence, comorbidities, and surgical outcomes in California. Blood Adv, 1, 1287–1295.

Baek

S.H.

, Kim

K.I.

, Yoon

K.S.

, Kim

T.H.

, and Kim

S.Y.

(2017). Genome-wide association scans for idiopathic osteonecrosis of the femoral head in a Korean population. Mol Med Rep, 15, 750–758.

Beyens

, Daroszewska

, de Freitas

, Fransen

, Vanhoenacker

, Verbruggen

, et al. (2007). Identification of sex-specific associationsbetween polymorphisms of the osteoprotegerin gene, TNFRSF11B, and Paget's disease of bone. J Bone Miner Res, 22, 1062–1071.

Brandstrom

, Stiger

, Michaelsson

, Kiss

, Balla

, Nagy

, et al. (1999). Polymorphism in the promoter region of the humen gene for osteoprotegrin: correlation with bone mineral density. Eur J Endocrinol 14, s134.

Cui

J.J.

, Wang

L.Y.

, Zhu

, Gong

W.J.

, Zhou

H.H.

, Liu

Z.Q.

, et al. (2017). Gene-gene and gene-environment interactions influence platinum-based chemotherapy response and toxicity in non-small cell lung cancer patients. Sci Rep 7, Article. ID5082.

Goto

, Hozumi

, Osaki

, Sakamoto

, Yonekura

, Tomita

, et al. (2011). Primary human bone marrow adipocytes support TNF-α-induced osteoclast differentiation and function through RANKL expression. Cytokine, 56, 662–668.

, Santos

, Ngo

, Fan

, Singh

P.P.

, Fingerle-Rowson

, et al. (2015). Macrophage migration inhibitory factor is essential for osteoclastogenic mechanisms in vitro and in vivo mouse model of arthritis. Cytokine, 72, 135–145.

Guo

, Wang

J.T.

, Liu

, Li

, Yang

T.L.

, Zhang

X.W.

, et al. (2012). Are bone mineral density loci associated with hip osteoporotic fractures?. A validation study on previously reported genome-wide association loci in a Chinese population. Genet Mol Res, 11, 202–210.

Hao

, Guo

, Xu

, Qi

, Wang

, Lu

, et al. (2018). The relationship between apolipoprotein genes polymorphisms and susceptibility to osteonecrosis of the femoral head: a meta-analysis. Lipids Health Dis 17, 192.

10.

Hein

, Bayer

C.M.

, Schraude

M.G.

, Häberle

, Heusinger

, Strick

, et al. (2014). Polymorphisms in the RANK/RANKLgenes and their effect on bone specific prognosis in breast cancer patients. 2014, Article. ID842452.

11.

Karol

S.E.

, Mattano

L.A.

Jr. , Yang

, Maloney

K.W.

, Smith

, Liu

, et al. (2016). Genetic risk factors for the development of osteonecrosis in children under age 10 treated for acute lymphoblastic leukemia. Blood, 127, 558–564.

12.

Kong

Y.Y.

, Yoshida

, Sarosi

, Tan

H.L.

, Timms

, Capparelli

, et al. (1999). OPGL is a key regulator of osteoclastogenesis, lymphocyte development and lumpho-nodeorganogenesis. Nature, 397, 315–323.

13.

Kubo

, Ueshima

, Saito

, Ishida

, Arai

, and Fujiwara

(2016). Clinical and basic research on steroid-induced osteonecrosis of the femoral head in Japan. J Orthop Sci, 21, 407–413.

14.

Kühn

M.C.

, Willenberg

H.S.

, Schott

, Papewalis

, Stumpf

, Flohé

, et al. (2012). Adipocyte-secreted factors increase osteoblast proliferation and the OPG/RANKLratio to influence osteoclast formation. Mol Cell Endocrinol, 349, 180–188.

15.

, Sarosi

, Yan

X.Q.

, Morony

, Capparelli

, Tan

H.L.

, et al. (2000). RANK is the intrinsic hematopoietic cell surface receptor that controls osteoclastogenesis and regulation of bone mass and calcium metabolism. Proc Natl Acad Sci U S A, 97, 1566–1571.

16.

, Wei

, and Hakonarson

(2016a). Application of computational methods in genetic study of inflammatory bowel disease. World J Gastroenterol, 22, 949–960.

17.

, Jiang

, and Du

(2017). Association between osteoprotegerin gene T950C polymorphism and osteoporosis risk in the Chinese population: evidence via meta-analysis. PLoS One 12, e0189825.

18.

, Wang

, Guo

, Wang

, Ouyang

, Cao

, et al. (2016b). OPG and RANKL polymorphisms are associated withalcohol-inducedosteonecrosis of the femoral head in the north area of China population in men. Medicine (Baltimore) 95, e3981.

19.

Meng

, Groth

, Quinn

J.R.

, Bisognano

, and Wu

T.T.

(2017). An exploration of gene-gene interactions and their effects on hypertension. Int J Genomics 2017, 7208318.

20.

Nanke

, Kobashigawa

, Yago

, Kawamoto

, Yamanaka

, and Kotake

(2016). RANK expression and osteoclastogenesis in human monocytes in peripheral blood from rheumatoid arthritis patients. Biomed Res Int 2016, 4874195.

21.

Ruyssen-Witrand

, Degboé

, Cantagrel

, Nigon

, Lukas

, Scaramuzzino

, et al. (2016). Association between RANK, RANKL and OPG polymorphisms with ACPA and erosions in rheumatoid arthritis: results from a meta-analysis involving three French cohorts. RMD Open 2, e000226.

22.

Sakamoto

, Yamamoto

, Sugano

, Takahashi

, Watanabe

, Atsumi

, et al. (2017). Genome-wide association study of idiopathic osteonecrosis of the femoral head. Sci Rep 7, 150357.

23.

Song

G.G.

, and Lee

Y.H.

(2017). Association of eNOS polymorphisms with susceptibility to osteonecrosis of the femur head: a meta-analysis. Z Rheumatol, 76, 267–273.

24.

Song

, Du

Z.W.

, Yang

Q.W.

, Ren

, Wang

Q.Y.

, Wang

, et al. (2017). Association of genes variants in RANKL/RANK/OPG signaling pathway with the development of osteonecrosis of the femoral head in Chinese population. Int J Med Sci, 14, 690–697.

25.

, Chen

, Yan

, Liang

, Xie

, Chen

, et al. (2015). Association between TRAF6 gene polymorphisms and susceptibility of ischemic stroke in Southern Chinese Han population. J Mol Neurosci, 57, 386–392.

26.

Sun

, Chen

, Lin

, Yu

, Zhao

, and Wang

(2014). The influence of osteoprotegerin genetic polymorphisms on bone mineral density and osteoporosis in Chinese postmenopausal women. Int Immunopharmacol, 22, 200–203.

27.

, Li

, Pei

, and Hu

(2017). Effects of Brucine on theOPG/RANKL/RANKsignalingpathwayin MDA-MB-231 and MC3T3-E1 cell coculture system. Evid Based Complement Alternat Med 2017, 1693643.

28.

Yerges

, Klei

, Cauley

J.A.

, Roeder

, Kammerer

C.M.

, Moffett

S.P.

, et al. (2009). High-density association study of 383 candidate genes for volumetric BMD at the femoral neck and lumbarspine among older men. J Bone Miner Res, 12, 2039–2049.

29.

Zhao

D.W.

, Yu

, Hu

, Wang

, Yang

, Wang

B.J.

, et al. (2015). Prevalence of nontraumatic osteonecrosis of the femoral head and its associated risk factors in the Chinese population: results from a nationally representative survey. Chin Med J (Engl), 21, 2843–2850.

30.

Zuk

, Hechter

, Sunyaev

S.R.

, and Lander

E.S.

(2012). The mystery of missing heritability: genetic interactions create phantom heritability. Proc Natl Acad Sci U S A, 109, 1193–1198.