Interaction of Wnt pathway related variants with type 2 diabetes in a Chinese Han population.

Aims. Epistasis from gene set based on the function-related genes may confer to the susceptibility of type 2 diabetes (T2D). The Wnt pathway has been reported to play an important role in the pathogenesis of T2D. Here we applied tag SNPs to explore the association between epistasis among genes from Wnt and T2D in the Han Chinese population. Methods. Variants of fourteen genes selected from Wnt pathways were performed to analyze epistasis. Gene-gene interactions in case-control samples were identified by generalized multifactor dimensionality reduction (GMDR) method. We performed a case-controlled association analysis on a total of 1,026 individual with T2D and 1,157 controls via tag SNPs in Wnt pathway. Results. In single-locus analysis, SNPs in four genes were significantly associated with T2D adjusted for multiple testing (rs7903146(C) in TCF7L2, p = 3.21∗10(-3), OR = 1.39, 95% CI [1.31-1.47], rs12904944(G) in SMAD3, p = 2.51∗10(-3), OR = 1.39, 95% CI [1.31-1.47], rs2273368(C) in WNT2B, p = 4.46∗10(-3), OR = 1.23, 95% CI [1.11-1.32], rs6902123(C) in PPARD, p = 1.14∗10(-2), OR = 1.40, 95% CI [1.32-1.48]). The haplotype TGC constructed by TCF7L2 (rs7903146), DKK1 (rs2241529) and BTRC (rs4436485) showed a significant association with T2D (OR = 0.750, 95% CI [0.579-0.972], P = 0.03). For epistasis analysis, the optimized combination was the two locus model of WNT2B rs2273368 and TCF7L2rs7903146, which had the maximum cross-validation consistency. This was 9 out of 10 for the sign test at 0.0107 level. The best combination increased the risk of T2D by 1.47 times (95% CI [1.13-1.91], p = 0.0039). Conclusions. Epistasis between TCF7L2 and WNT2B is associated with the susceptibility of T2D in a Han Chinese population. Our results were compatible with the idea of the complex nature of T2D that would have been missed using conventional tools.

Introduction

Genome-wide-association studies (GWAs) detect half a million to a million markers at once. According to comparative analysis of different GWAs on the same disease, genuine susceptible loci could rank at moderate levels of significance rather than high (Zeggini et al., 2008; Baum et al., 2008). And multiple testing corrections may neglect many modest association signals (Wang et al., 2011). To compensate for the shortcomings of GWAs, rather than focusing on individual SNPs, epistasis from gene sets based on the function-related genes, which jointly take into account multiple variants, have been developed (Zhou et al., 2012). Emerging reports suggest that epistasis is responsible for a proportion of complex human disease (McClellan & King, 2010).

Pathway-analysis of Type 2 Diabetes (T2D) based on the Welcome Trust Case Control Consortium (WTCCC) shows that Wnt pathway is associated with T2D (Zeggini et al., 2007). Wnt signaling pathway is initially recognized in colon cancer research (Peifer & Polakis, 2000). The heterodimerisation of free β-cat with one of the four members of the TCF family, TCF7L2 (TCF3) being themajor partner of β-cat in the intestinal epithelia, form the bipartite transcription factor β-catenin/T cell factor (β-cat/TCF), which is the major factor that effects the Wnt pathway(s) (Doble & Woodgett, 2003). In the absence of WNT, the free β-cat is phosphorylated by one complex, then destroyed. Following the Wnt ligands binding, the destructive complex is disrupted, and then β-cat enters the nucleus and formulates the β-cat/TCF downstream target genes. Wnt pathway is involved in the secretion of the incretin hormone glucagon-like peptide-1 (GLP-1) by the intestinal endocrine L-cells (Yi, Brubaker & Jin, 2005). Thus, alteration in this pathway could reduce GLP-1 levels which, in turn, could affect insulin production. GLP-1, accompanied by insulin, plays a vital role in blood glucose homeostasis. It has been postulated that variants in this pathway may confer a predisposition of T2D by altering GLP-1 levels (Smith, 2007). The role of epistasis among genes from Wnt pathway to genetic susceptibility to T2D remains unclear.

Based on the pathway-analysis of T2D from WTCCC, our study aims at the exploration of the contribution of epistasis from variants among Wnt pathway to the susceptibility to T2D. Selected SNPs of 14 genes from Wnt pathway were explored in unrelated T2D patients and non-diabetic control subjects, then epistasis among these genes were examined. Here we present the results in Chinese Han individuals.

Method

Ethics statement

The study was conducted with the approval from the Ethics Committee of Beijing Tongren Hospital, Capital Medical University, approval number 2013.12.24. Written informed consent was obtained from each participant.

Study populations and genotype

The Beijing Diabetic Complications Survey is an ongoing database that was established for studying the genetic basis of T2D. The study subjects included 1,026 T2D cases and 1,157 controls from northern China, who were acquired from this database (Yang et al., 2010; Zhou et al., 2010b). Single nucleotide polymorphisms (SNPs) from the Wnt pathway were assessed. Cases were chosen based on the presence of type 2 diabetes. Control groups were individuals who had fasting plasma glucose <5.6 mmol/l, 2 h post-OGTT plasma glucose <7.8 mmol/l. Blood pressure, serum creatinine, and ages that were matched in all subjects of both groups. Genotyping was completed using a MassARRAY platform (Sequenom, San Diego, California, USA).

Gene and SNP selection

Individual level GWA data of T2D from the WTCCC were completed to perform pathway-analysis using the Gene Set Enrichment Algorithm (GSEA) (Zeggini et al., 2007; Wang, Li & Bucan, 2007). A total of 390,025 autosomal SNPs were genotyped across 1,924 case subjects and 2,938 population-based control subjects. All human KEGG pathways were used in the pathway-analysis (Elbers et al., 2009). A total of 139 genes of Wnt pathway are included from KEGG. The most strongly associated SNP from a gene ±200 kb of flanking sequence either side is used to represent each gene in the genome.

Final results showed that Wnt pathway was associated with T2D. Top SNPs of fourteen genes reaching p < 0.01 were included. If selected SNPs were scare or homomorphism based on HapMap CHB data, tags SNPs wereincluded from the selected genes according to the following principles: (1) threshold of a minor allele frequency (MAF) >0.05. (2) linkage disequilibrium (LD) between SNPs with the threshold of r2 > 0.4 and D > 0.7 according to Haploview (v.4.2) based on HapMap CHB data sets (v.3, release 2).Datum was excluded if the allele call rate was less than 95%, or p value of the Hardy-Weinberg equilibrium (HWE) was less than 0.001. (3) functional relevance and importance.

The included SNPs were rs1796390 in PRICKLE1, rs2273368 in WNT2B, rs33830 in FBXW11, rs7903146 in TCF7L2, rs6902123 in PPARD, rs2242445 in CSNK2A2, rs3860318 in CTNNBIP1, rs4436485 in BTRC, rs12904944 in SMAD3, rs1049612 in CCND2, rs2656272 in SOX17, rs3810765 in SFRP2, rs3741571 in FZD10 and rs2241529 in DKK1.

Statistic methods

SPSS version 19.0 for Windows (SPSS Inc., Chicago, Illinois, USA) was used to perform statistical analysis. Means ± SD were used to show the continuous data. The between-group data for categorical datum were analyzed by the X2 test. Allele frequencies were calculated from the genotypes of the subjects. Chi-squared-test or Fisher’s exact test were performed to compare the differences in allele and genotype frequencies between cases and controls.

To evaluate the interaction among a large number of loci in our research, we used the MDR (multifactor dimensionality reduction) method to identify high order gene–gene interactions in our samples. The MDR method was described in our paper (Zhou et al., 2012). All possible SNP interactions were tested using 10-fold cross-validation in an exhaustive search. The MDR method was described in detail elsewhere (Dong et al., 2003). Briefly, in an analysis of n-way interactions, the n-dimensional space formed by all possible combinations of values (classes) of a given set of n variables (in this case, SNPs) was reduced to a single dimension by reclassifying each class as either high risk or low risk according to the relative proportion of cases to controls in that class. This statistical method included a 10-fold cross-validation and permutation-testing procedure to minimize false positive results by multiple examinations of the data. With 10-fold cross-validation, the data were divided into 10 equal parts, and the model was developed on 9/10 of the data (training set) and then tested on 1/10 of the remaining data (testing set). This was repeated for each possible 9/10 and 1/10 of the data, and the resulting 10 prediction errors were averaged. The combination with the lowest prediction error was reported. Finally, hypothesis testing for this selected model could then be performed by evaluating the consistency of the model across cross-validation data sets. Under the null hypothesis that no association was derived from 1,000 permutations, the average cross-validation consistency from the observed data was compared to the distribution of average consistencies. The null hypothesis was rejected when the P-value derived from the permutations was ≤0.05.

Generalized multifactor dimensionality reduction (GMDR) was based on the score of a generalized linear model, of which the original MDR method was a special case (Lou et al., 2007). Then logistic regression models were employed to confirm the best combination of loci identified in the GMDR. Armitage’s test for trend calculated power analysis of case-control samples (alpha has been set at 0.05) (Slager & Schaid, 2001).

Results

Population characteristics

The characteristics of the population were provided. Table 1 illustrated the detailed clinical features of the validation cohorts. Age, blood pressure, and BMI were comparable between the two groups. Significance was noted for total cholesterol, low-density lipoprotein, high-density lipoprotein and triglyceride (p < 0.01). No deviation from the Hardy–Weinberg Equilibrium was observed for all studied polymorphisms in controls (p > 0.05).

Table 1

Clinical characteristics of the cohorts.

Age, body mass index (BMI), systolic blood pressure (SBP), diastolic bloodpressure (DBP), TC, HDL-c, LDL-c, Cr and fast plasma glucose (FPG) values are given as mean (SD); TG values are given as median (range).

Characteristic	Type 2 diabetes	Controls
Number	1,026	1,157
Men/Women (n/n)	523/503	569/588
Age (years)	59.4 ± 9.3	59.6 ± 9.4
BMI (kg/m2)	26.5 ± 3.73	25.8 ± 3.45
Waist to hip ratio	0.91 ± 0.07	0.86 ± 0.08
Fasting blood glucose (mmol/l)	7.82 (6.69–9.74)	4.83 (4.56–5.46)
TG (mmol/l)	1.75 (1.26–2.55)	1.19 (0.76–1.66)
TC (mmol/l)	4.84 (4.27–5.52)	4.34 (3.83–4.87)
HDL-c (mmol/l)	1.13 (0.96–1.31)	1.19 (1.02–1.36)
LDL-c (mmol/l)	3.09 (2.52–3.63)	2.79 (2.31–3.26)
SBP (mmHg)	137 (126–151)	129 (119–144)
DBP (mmHg)	81 (79–90)	81 (74–89)

Single locus analysis

The results of the single-locus analyses for the case-control population indicated that three SNPs were significantly associated with T2D (Table 2, p < 0.01). After adjusting of significance threshold for multiple testing (p < 0.01/14 = 0.00357), two SNPs were still significantly associated with T2D (rs7903146C in TCF7L2 p = 3.21∗10−3, OR = 1.39, 95% CI [1.31–1.47], rs12904944G in SMAD3, p = 2.51∗10−3, OR = 1.39, 95% CI [1.31–1.47]).

Table 2

The study results for selected SNPs in Wnt-signalling pathway.

Gene	SNP	Chr.	Position	Location	Risk/non–risk allele	Control risk frequency	T2D risk frequency	P
PRICKLE1	rs1796390	12	42878604	intron	T/C	0.294	0.311	0.08342
WNT2B	rs2273368	1	112521149	intron	C/T	0.453	0.492	0.00446
FBXW11	rs33830	5	171721166	intron	A/T	0.30	0.323	0.05348
TCF7L2	rs7903146	10	112998590	intron	C/T	0.047	0.090	0.003209
PPARD	rs6902123	6	35362644	intron	C/T	0.06	0.893	0.01141
CSNK2A2	rs2242445	16	58165299	intron	G/T	0.144	0.168	0.0582
CTNNBIP1	rs3860318	1	9864887	intron	A/T	0.188	0.199	0.1372
BTRC	rs4436485	10	101485604	intron	G/C	0.118	0.135	0.09461
SMAD3	rs12904944	15	67069436	intron	G/A	0.42	0.473	0.002511
CCND2	rs1049612	12	4303596	3′ region	A/G	0.10	0.129	0.04518
SOX17	rs2656272	8	54349340	intron	A/G	0.38	0.398	0.08346
SFRP2	rs3810765	4	153788328	intron	G/A	0.45	0.476	0.05414
FZD10	rs3741571	12	130152232	transcript variant	T/C	0.05	0.062	0.1568
DKK1	rs2241529	10	52314997	exon	A/G	0.332	0.358	0.04561

Interaction among genes in Wnt pathway

The contribution of combinations of the variants to the etiology of T2D was assessed by GMDR, which showed that WNT2B rs2273368 and TCF7L2 rs7903146 was the best combination. The optimum had the highest level of testing balance accuracy (0.5620), and the maximum cross-validation consistency, this was 9 out of 10 for the sign test at 0.0107 level (see Table 3). The two-locus genotype combinations were classified into high- or low-risk groups. In the chi-squared test, the OR of high-risk combination of the two-locus model increased the risk of T2D by 1.47 times (95% CI [1.13–1.91], p = 0.0039). To further elucidate the contribution of interactions on the predisposition of T2D, the logistic analysis was performed on the basis of two variants identified by GMDR. No significant multiplicative deviation of gene–gene interaction was shown in the logistic regression model.

Table 3

GMDR results of multi-locus interaction with T2D.

No. of loci with T2D	Best model*	Training balance accuracy	Testing balance accuracy	Cross-validation consistency	Sign test (p)
1	TCF7L2	0.5464	0.5430	6	7 (0.1719)
2	WNT2B, TCF7L2	0.5619	0.5620	10	9 (0.0107)
3	WNT2B, TCF7L2, PPARD	0.5350	0.5429	8	5 (0.6230)
4	PRICKLE1, WNT2B, TCF7L2, PPARD	0.5276	0.5429	8	6 (0.3770)
5	PRICKLE1, WNT2B, TCF7L2, PPARD, CSNK2A2	0.5295	0.5444	7	6 (0.3770)
6	PRICKLE1,WNT2B, TCF7L2, PPARD, CSNK2A2, DKK1	0.5307	0.5460	7	5 (0.6230)

Notes.

Adjusted by age and sex.

Since GMDR was based on MDR, analysis from MDR was also performed. The accuracy of the two-locus model without co-variable adjustment was similar to that with adjustment (see Table 4).

Table 4

Comparison of best models from MDR and GMDR for prediction of T2D.

	Training balance accuracy	Testing balance accuracy	Cross-validation consistency	Sign test (p)
MDR	0.5619	0.5620	10	0.0107
GMDR	0.5619	0.5619	10	0.0107

Discussion

The potential for gene–gene interaction has been proposed to be one of the possible reasons for the so-called “missing heritability” of T2D (Yang et al., 2010). Contribution of epistasis among variants from gene sets based on the function-related genes still remains unclear. Epistasis from gene sets that were evaluated may provide a cumulative contribution and new insights into the biology of T2D. Our results evidenced the interaction among Wnt pathway related genes and the correlation of susceptibility to T2D.

Single locus analysis indicated significant association with T2D in our research (rs7903146 in TCF7L2, rs12904944 in SMAD3). The first GWAS on T2D showed strong signal for TCF7L2 (Sladek et al., 2007). TCF7L2 harbors common genetic variants and a strong effect on T2D (Zhou et al., 2014). Our study substantiated the association between variant of TCF7L2 and T2D which was validated in another study with Chinese Han patients (Wang et al., 2013). The SNP rs7903146 (TCF7L2) has been associated with early-onset T2D in some populations (Iwata et al., 2012). Since the absence of age at diagnosis, we could not analyze association between early-onset T2D and susceptible genes.There is a novel association between SMAD3—rs12904944 and T2D in our study. As an important metabolic regulator, SMAD3 is linked to the pathogenesis of T2D (Tan et al., 2011). MAF of rs12904944 in our research is similar to that of CHB population in HapMap, while higher than that of Europeans population in HapMap (MAF: 0.42 vs. 0.40 vs. 0.33). The underlying mechanisms responsible for the association between these variants and T2D remain to be found.

In recent years, more and more studies have focused on gene–gene interactions (Moore & Williams, 2009; Jiang, Barmada & Visweswaran, 2010), which are partly due to the appearance of new statistical theory. GMDR is proved to be a useful statistical tool to detect gene–gene interactions while avoiding “the dimension curse” (Ritchie et al., 2001). In our study, GMDR suggested that the combination between WNT2B rs2273368 and TCF7L2 rs7903146 was the best model. However, the interaction was not validated by the logistic regression analysis. The possible reason for these inconsistent results is that GMDR did not detect the interaction defined by “deviation from the multiplicative” in the logistic regression model. The significant results from GMDR only show that the combination of different loci may increase the risk of T2D and the interaction among the loci may refer to either multiplicative, deviation from the multiplicative or departure from additivity (Zhou et al., 2012). WNT2b, which activates the canonical Wnt pathway, is highly unregulated in T2D (Lee et al., 2008). The association between WNT2B and T2D was novel in our research, this significance disappeared after multiple adjusting. The mechanism of epistasis between TCF7L2 and WNT2B need to be further explored.

Variants conferring modest disease risk may not reveal themselves in multiple underpowered GWA studies, but can be readily identified by epistasis from gene sets based on function-related genes in a single study. Our results validated this strategy. Analysis of epistasis from one gene set can provide complementary information to conventional single-marker analysis (Wang, Li & Hakonarson, 2010). Therefore, genes that interact with others in one gene set can serve as candidates for further replication and functional studies (Zhou et al., 2010a).

In conclusion, our study supports epistasis association with the predisposition of T2D between TCF7L2 and WNT2B. We believe these findings will guide further investigations in gene detection, and hopefully lead to a better understanding of T2D pathogenesis.

Click here for additional data file.