Introduction

Gout has been a great concern as the most prevalent arthritis associated with impaired quality of life.1, 2 Gout is developed with a preliminary increase of serum uric acid level and a subsequent deposition of monosodium urate crystals in and around joints. The serum uric acid level is, therefore, considered the most critical risk factor for susceptibility to gout. As inefficient excretion of uric acid is more influential in elevating serum uric acid level and predisposing people to gout than its overproduction, renal mechanisms have received an overwhelming attention in terms of explaining gout susceptibility.3 Especially, recent genetic studies have primarily focused on genes regulating renal urate transport.4 For example, genetic associations of serum uric acid level and/or gout susceptibility were repeatedly observed with nucleotide sequence variants of the gene encoding GLUT9 (SLC2A9), a urate transporter in proximal renal tubules.5, 6, 7 Genetic associations were also identified with other transporters such as ABCG2,5, 8 NPT1 (SLC17A1),9 SLC17A35 and URAT1 (SLC22CA12).10, 11 Phenotypic variability for gout susceptibility has been quite limitedly explained by a small number of genetic factors identified so far,12 especially considering its large heritability (0.63).13 A small number of genetic associations have been discovered even in recent genome-wide association studies (GWAS).14 This was because many false-negative genetic associations were produced by a multiple test with hundreds of thousands of nucleotide markers.

The current study aimed to conduct a GWAS for gout susceptibility in a large-scale cohort, to identify functional associations of cis-regulatory regions with gout susceptibility, to compare the functional associations with gout susceptibility with those with serum uric acid level, and to show a genetic architecture for gout susceptibility with interaction effects among genes within each functional association.

Materials and methods

Subjects and data

The Korea Association REsource (KARE) Analysis Consortium was established to conduct large-scale GWAS for understanding of human genetic basis. The consortium has cohort data of 10 038 unrelated Korean individuals collected by the Korean Genome Epidemiology Study. The data collection was initiated in 2001, and thereafter, follow-up examinations for each participant have been conducted every 2 years. Genotypic data of the KARE were obtained using the Affymetrix Genome-Wide Human SNP Array 5.0 (Affymetrix, Santa Clara, CA, USA). A total of 8842 individuals were obtained after screening by genotype calling and quality control. For details, see Cho et al.15 However, four individuals without phenotype or covariate information were excluded in the current analysis. An underlying set of unphased genotypes for each individual in the cohort was imputed with the Japanese and Chinese HapMap phase 2 haplotype panel using IMPUTE software program (version 2).16

A total of 520 out of 8834 subjects were self-reported patients with gout, and they were assumed to be diagnosed with gout by rheumatologists according to the American College of Rheumatology diagnostic criteria.17 The other subjects in the cohort were all used as controls. The patients with gout were excluded in association analysis for serum uric acid level. Phenotypic characteristics of the patients were compared with those of controls in Supplementary Table 1.

Single locus association

Genotypic associations of each single-nucleotide polymorphism with susceptibility to gout and with serum uric acid level were tested with an additive model. The analytical model included age, gender and body mass index as covariates. Threshold of false-positive error in the significance test was 0.05, and the Bonferonni correction was introduced to control for occurrence of false-positive genetic associations. The association analyses were conducted using PLINK (version 1.06, http://pngu.mgh.harvard.edu/purcell/plink) and SPSS (version 12.0, SPSS, Chicago, IL, USA) software programs.

Genome-wide functional prediction

Genomic regions enrichment analysis was conducted to analyze functional associations of cis-regulatory regions identified by localized measurements of DNA-binding events using GREAT (http://great.stanford.edu/).18 The input genomic regions were determined on the basis of single-nucleotide polymorphisms with P<0.01 resulted from the GWAS. Biological processes provided by the Gene Ontology (http://www.geneontology.org/) were utilized as input annotation terms. The analysis incorporated distal binding sites up to 1 Mb, as well as proximal binding sites up to upstream 5 Kb and downstream 1 Kb. Significance of enrichment was determined by both a binomial test over genomic regions and a hypergeometric test over genes. To control incorrect rejections of null hypothesis, the significance thresholds were corrected by false discovery rate.

Entropy decomposition

Phenotypic variability of gout as a complex disease might be largely explained by epistasis. The epistatic effects among genes within each biological process were estimated in the current study. We first selected the optimal genes within each biological process for susceptibility to gout using multifactor dimensionality reduction (MDR) analysis and then estimated pairwise epistatic estimates using entropy decomposition.

The MDR collapses high-order dimensions of genotypic data produced with multiple loci into a single dimension, which enables to infer epistasis in a relatively small sample size by grouping genotypes.19 Thus, the MDR could identify multiple loci and their genotype combinations that were associated with gout susceptibility. A 10-fold cross-validation strategy was applied to the current MDR. Case and control data were randomly and equally partitioned into 10 pieces. Nine pieces were used as a training data set to build a genetic model for predicting susceptibility to gout, and the rest one piece was used as a testing data set to test the model established by using the training set. This analysis was conducted with each of the 10 possible partitions. The best n-locus model was selected on the basis of minimum classification error with 10-fold cross validation. Each of combined genotypes produced with every possible n-locus was determined as to whether it was a high or low risk, based on the case–control ratio. All the procedures explained above were replicated 100 times through shuffling data sets. The best models with 1-, 2-, 3- and 4-loci were determined with the maximum average estimates of cross-validation consistency and testing accuracy for the replicates. All the variants in the best models with 1-, 2-, 3- and 4-loci were selected for estimating their additive and interactive genetic effects. The MDR analysis was conducted with all the loci of each function identified by the genomic regions enrichment analysis using the MDR software package (http://www.multifactordimensionalityreduction.org).

Pairwise interactions among the selected loci associated with a susceptibility to gout were estimated as the portions of entropy removed for interaction information of pair variants after adjusting marginal entropy removed by each individual variant.20 Therefore, interaction analysis reduced uncertainty about either of two attributes with the knowledge of the other attribute. The interaction effect was partitioned into synergistic and redundant effects by the sign of its value. The entropy decomposition analysis was conducted utilizing the freely available Orange software (http://www.ailab.si/orange).

Results

The GWAS revealed no significant nucleotide sequence variants associated with gout susceptibility or with serum uric acid level by Bonferroni multiple test (P>1.42 × 10−7, Figure 1). Previously identified associations of gout susceptibility or serum uric acid level were replicated in the current study with nucleotide sequence variants in solute carrier family genes such as GLUT9, NPT1, PDZK1 and SLC17A3 (P<0.05, Table 1).

Figure 1
figure 1

Manhattan plot for genome-wide associations with gout susceptibility (a) and serum uric acid level (b). The y axis represents negative log-transformed P-values for their associations with nucleotide polymorphisms at their relative genomic positions by chromosome represented along the x axis. No nucleotide polymorphisms were identified for association with gout susceptibility or with serum uric acid level by Bonferroni multiple test (P<1.42 × 10−7). A full color version of this figure is available at the Journal of Human Genetics journal online.

Table 1 Associations of SNPs within previously identified genes with gout susceptibility and serum uric acid level

Genome-region set analysis revealed a variety of biological processes associated (false-discovery-rate-corrected P<0.05) with each of gout susceptibility (Table 2) and serum uric acid level (Table 3). The identified biological processes associated with gout susceptibility differed from those with serum uric acid level except for regulation of chemotaxis.

Table 2 Biological process associated with gout susceptibility using genomic cis-regulatory region set of proximal and distal binding sites
Table 3 Biological process associated with serum uric acid level using genomic cis-regulatory region set of proximal and distal binding sites

Individual and interactive genetic effects of nucleotide sequence variants were estimated for each of the identified biological processes using entropy decomposition (Figure 2). The entropy decomposition analysis included all the loci of the best 1-, 2-, 3- and 4-locus models selected by MDR (Supplementary Table 3). Some interactions were observed with entropy estimates larger than those of individual effects. For example, rs7035561 in AGTPBP1 gene and rs9960723 in NETO1 gene showed a strong synergistic epistasis, which explained phenotypic variability approximately four times larger than their individual effects. Furthermore, although marginal effect of rs9308056 in NPY5R was negligible, its interaction with the rs7035561 was considerably large. In general, strong epistatic associations were synergistic in regulations of neurological system process, cellular metabolic process and neurotransmitter transport, whereas redundant epistatic associations were observed in regulations of leukocyte migration and chemotaxis.

Figure 2
figure 2

Marginal and interactive genetic effects for gout susceptibility using entropy decomposition. The estimate of marginal genetic effect is presented in each box, and the estimate of interactive genetic effect is presented along with the line connecting the corresponding variants. The effect size of epistasis is proportional to the width of the corresponding line. The black line indicates synergistic (positive) epistasis, and the gray line indicates redundant (negative) epistasis. The effects were estimated within each biological process: leukocyte migration (a), chemotaxis (b), neurotransmitter transport (c), neurological system process (d) and cellular metabolic process (e).

Discussion

We conducted a GWAS for gout susceptibility and serum uric acid level with a large-cohort data set in the current study, but any nucleotide sequence variants were not identified by employing Bonferroni multiple test. In spite of their high heritability (0.63),13 complexity of the traits might lead to hard detection for genetic factors with small effects. A large number of false-negative associations could be produced by such a conservative testing. We applied a genomic regions enrichment analysis using variants with a marginal association to overcome the problem. Furthermore, we could exclude potential biases by incorporating two mutually beneficial tests: one was a binomial test explicitly accounting for the sizes of the regulatory domains of genes, and the other was a hypergeometric test counting each gene only once. As a result, a variety of annotation terms were identified with statistical enrichment for their functional associations with gout susceptibility. The annotation terms turned out to be quite different from those associated with serum uric acid level. This kind of a large divergence was unexpected, considering the serum uric acid level as the underlying cause of gout. Nevertheless, the annotation terms were identified for association with gout susceptibility, and specifically with serum uric acid level. For example, the annotation terms included ‘transforming growth factor (TGF)β receptor signaling pathway’ for serum uric acid level. This was explained well with the role of TGFβ1 as a promoter in synthesizing extracellular matrix. The change of extracellular matrix declines renal function and increases serum uric acid level.21 On the other hand, gout susceptibility was not associated with the term of ‘TGFβ receptor signaling pathway’, but significantly associated with the annotation term of ‘leukocyte migration’. Leukocyte adhesion in monosodium urate crystal-induced inflammation can be suppressed by TGFβ1, and the suppression would lead to gout.22 The discrepancy between the annotation terms associated with gout susceptibility and with serum uric acid level implied heterogenous underlying genetic factors for the gout and its underlying cause. The different genetic associations might explain that many people with hyperuricemia would be highly unlikely to develop gout. That is, gout susceptibility might be largely explained by a variety of genetic factors for inducing chronic inflammation after urate deposition.

Phenotypic variability of gout as a complex disease might be largely explained by interaction effects between genes.23 Including the interaction effects is essential in understanding genetic architecture for gout susceptibility, but may present a serious drawback of containing a large number of all the possible interaction effects in the analytical model. To overcome the problem, we employed a parsimonious interaction model with selected genes within each of the biological process terms in the current study. As a result, large interactive effects of the genetic factors for gout susceptibility were observed using an entropy-based approach. Especially, some interaction effects explained larger phenotypic variability than their marginal effects. For example, neurological system process-related genes showed their strong synergistic epistasis for gout susceptibility.

One thing to be noted is a potential heterogeneity of gout genetic architecture by gender. The current study was not balanced by gender, and thus, a larger portion of genetic factors was attributed to genetic mechanism for susceptibility to gout in males. However, this unbalance in the current study was not too extreme with a larger female ratio of gout (21%; Supplementary Table 1) than that in usual (5%).24 The larger proportion might be explained partially by having a larger female ratio (10%)25 of gout in Korea, and more importantly by including only healthy subjects on the starting point of the cohort, that is, male subjects with gout at that point might have been excluded.

In conclusion, the current study provided the first evidence that genetic factors for gout susceptibility considerably differed from those for serum uric acid level. This implied that gout susceptibility was greatly attributed to genetic factors irrelevant to pathogenesis of uric acid, and thus understanding genetic architecture of gout susceptibility would also require research endeavors on identification of such genetic factors. Interaction effects between genes turned out to be not negligible for explaining phenotypic variability for gout susceptibility. To validate our findings, the current study calls for rigorous replications in other ethnic populations. Furthermore, comparing subjects with gout with those with asymptomatic hyperuricemia would confirm the heterogeneous genetic factors for gout susceptibility and serum uric acid level, which eventually help understand decomposed genetic factors for such complex disease.