Genomic dissection of 43 serum urate-associated loci provides multiple insights into molecular mechanisms of urate control

Abstract

High serum urate is a prerequisite for gout and associated with metabolic disease. Genome-wide association studies (GWAS) have reported dozens of loci associated with serum urate control; however, there has been little progress in understanding the molecular basis of the associated loci. Here, we employed trans-ancestral meta-analysis using data from European and East Asian populations to identify 10 new loci for serum urate levels. Genome-wide colocalization with cis-expression quantitative trait loci (eQTL) identified a further five new candidate loci. By cis- and trans-eQTL colocalization analysis, we identified 34 and 20 genes, respectively, where the causal eQTL variant has a high likelihood that it is shared with the serum urate-associated locus. One new locus identified was SLC22A9 that encodes organic anion transporter 7 (OAT7). We demonstrate that OAT7 is a very weak urate-butyrate exchanger. Newly implicated genes identified in the eQTL analysis include those encoding proteins that make up the dystrophin complex, a scaffold for signaling proteins and transporters at the cell membrane; MLXIP that, with the previously identified MLXIPL, is a transcription factor that may regulate serum urate via the pentose–phosphate pathway and MRPS7 and IDH2 that encode proteins necessary for mitochondrial function. Functional fine mapping identified six loci (RREB1, INHBC, HLF, UBE2Q2, SFMBT1 and HNF4G) with colocalized eQTL containing putative causal SNPs. This systematic analysis of serum urate GWAS loci identified candidate causal genes at 24 loci and a network of previously unidentified genes likely involved in control of serum urate levels, further illuminating the molecular mechanisms of urate control.

Introduction

Elevated serum urate (hyperuricemia) is causal of gout, an inflammatory arthritis increasing in prevalence worldwide (1,2). Monosodium urate crystals, which form in individuals with hyperuricemia, can activate the NLRP3-inflammasome of resident macrophages to mediate an IL-1β-stimulated gout flare (3). Long established genome-wide association studies (GWAS) (4,5) have reported 28 loci associated with serum urate levels in European and East Asian sample sets with a more recent East Asian study reporting an additional eight loci (6). The loci of strongest effect are dominated by renal and gut transporters of urate, with two loci (SLC2A9 and ABCG2) together explaining up to 5% of variance in serum urate levels in Europeans (4). Most of these 36 loci also associate with gout in multiple ancestral groups (4,7–9). There has, however, been little progress on understanding the molecular basis of the association for the various loci. Probable causal genes have been identified at only about one-fifth of the 36 loci (10–12), with strong evidence for causality for variants identified at ABCG2 (rs2231142; Q141K) and PDZK1 (rs1967017) (11,13–17).

There are a number of resources and analytical techniques that can be applied to the summary statistics of GWAS. Differences in underlying linkage disequilibrium (LD) structure between ancestral groups can be leveraged to amplify signal of association at shared causal variants (18,19). The epigenomics roadmap and ENCODE projects have generated a large resource of cell- and organ-specific regulatory regions (20,21). This information can be used to discover the cell-type-specific regulatory regions that are known to be overrepresented in the heritability of a typical complex trait. Variants in regulatory regions identified by the epigenomics roadmap, and ENCODE can be further analyzed with functional annotation fine-mapping tools to identify candidate causal variants. Once credible sets of causal variants have been identified, expression quantitative trait loci (eQTL) sample sets [e.g. Genotype-Tissue Expression (GTEx) project (22)] can be used to translate from causal variants to affected genes, thus informing the design of functional experiments for insights into molecular pathogenic pathways. Because sample sizes for eQTL studies are relatively modest (<1000), colocalization analyses of GWAS and eQTL data have remained primarily focused on cis-eQTL. However, recent methods that integrate high-resolution genomic interaction data with eQTL data can reduce the number of trans-eQTL investigated substantially (23,24), although a limitation of this filtering approach is that it excludes trans-eQTL not mediated by genomic interactions (25). Despite this limitation, integrating genomics interaction-filtered trans-eQTL signals with GWAS allows expansion of our view of how GWAS associations underpin gene expression (24).

In this study, we integrated these analytical approaches with the summary statistics of two serum urate GWAS from European and East Asian individuals (4,5). By meta-analysis and colocalization analysis of serum urate and eQTL signals, we identified 15 new serum urate loci, identified 54 candidate causal genes connected to 29 loci, revealed the cell types that are enriched in serum urate heritability and used this functional information to identify credible sets of causal variants using trans-ancestral fine mapping.

Figure 1

Analysis flowchart.

Results

Trans-ancestral meta-analysis identifies 10 new loci associated with serum urate levels

The analysis approach for this study is summarized in Figure 1. Z-scores were imputed into the European (4) and East Asian (5) summary statistics using reference haplotypes from the Phase 3 1000 Genomes release and combined by meta-analysis (Fig. 2 and Supplementary Material, Fig. S1). Study-specific results revealed three new loci at chromosome 11 in the East Asian sample set (Chr11, 63.2–67.2 Mb, SLC22A9, PLA2G16, AIP) in addition to those reported as genome-wide significant in the original GWAS (5) (Table 1; Fig. 2 and Supplementary Material, Fig. S2). All loci reported in the original GWAS reports (4,5) were also detected in the trans-ancestral meta-analysis. However, the separate signals at SLC22A11 and SLC22A12 (Chr11, 64.4 Mb) reported by Köttgen et al. (4) are reported as one signal in the trans-ancestral meta-analysis and an additional signal was detected at Chr11 65.4 Mb (RELA). The trans-ancestral meta-analysis identified seven new loci (Chr4, 81.2 Mb, FGF5; Chr5, 40.0 Mb, LINC00603; Chr6, 32.7 Mb, HLA-DQB1; Chr9, 33.2 Mb, B4GALT1; Chr10, 60.3 Mb, BICC1; Chr11, 63.9 Mb, FLRT1; Chr11, 119.2 Mb, USP2) (Fig. 2 and Supplementary Material, Fig. S3). Of the 10 new loci identified (seven from the trans-ancestral meta-analysis and three in the East Asian-specific analysis), five mapped within an extended Chr11 locus (63.2–67.2 Mb) that encompassed the previously identified SLC22A11, SLC22A12 and OVOL1/RELA loci (4,5). In the East Asian GWAS, the peak marker falls outside the RELA locus (Supplementary Material, Fig. S4). On closer inspection of the association signal from the region within and surrounding the RELA locus, it is clear that the causal variants in the East Asian population are not the same as in the European population (Supplementary Material, Fig. S4). On Chr6, given the association of the HLA-DQB1 locus with T-cell-mediated autoimmunity (5), we also investigated if the lead HLA-DQB1 SNP (rs2858330) was associated with other phenotypes using Phenoscanner (26). This analysis identified 23 associations at P < 5E-8, including associations with autoimmune disorders, such as ulcerative colitis, rheumatoid arthritis and asthma. This indicates that the HLA-DQB1 signal in the serum urate GWAS is not distinct from the association of this region with autoimmunity. The 35 loci found in Europeans explain 6.9% of variance in age and sex-adjusted serum urate levels. In summary, a total of 38 loci associated with serum urate concentration at a genome-wide level of significance were identified by this analysis.

Conditional analysis identifies eight additional variants associated with serum urate

Using the European summary statistics, a conditional and joint analysis was performed with the objective of identifying independent genetic effects. Conditional and joint analysis identified an additional four genome-wide significant associations at SLC2A9 and one at each of ABCG2, TMEM171, SLC16A9 and SLC22A11/A12 (Table 1). The conditional analysis was limited to five independent associations at each locus; therefore, it remains possible that there are additional unidentified associations at these loci. In a joint model, these five loci explained an additional 0.54% of the variance of age- and sex-adjusted serum urate levels (Table 1).

Population-specific associations with serum urate levels

LocusZoom plots from each population were visually compared to the trans-ancestral meta-analysis to identify population-specific and shared patterns of association. For 16 loci (ABCG2, B4GALT1, BCAS3, FGF5, BICC1, HFN4G, IGF1R, INHBB, NFAT5, PDZK1, QRICH2, SLC16A9, SLC17A1, TMEM171, TRIM46 and UBE2Q2), the pattern of association was consistent between the East Asian and European GWAS, suggesting strong similarity between the underlying haplotypic structure and casual variant(s) (Supplementary Material, Fig. S5). The MAF locus contains two association signals in the East Asian population, one that is shared with the European population and one that is specific to the East Asian population (Supplementary Material, Fig. S6) (12). The lead SNP for the East Asian population is rs889472; this variant is common in both European (C-allele = 0.38) and East Asian (C-allele = 0.60) individuals from the 1000 Genomes Project; yet there is no serum urate association signal in the European population. Two other East Asian-specific signals were identified on chromosome 11 near the SLC22A9 and PLA2G16 genes in addition to the previously mentioned East Asian-specific signal at the RELA locus. These loci, in combination with the conditionally independent and trans-ancestral associations, mean that there are seven independent associations on chromosome 11 between 63.1 and 67.3 Mb.

Table 1

SNPs associated with serum urate concentrations by meta-analysis individuals of European and East Asian ancestry

SNP	Chr:bp	Closest gene	A1a	A2	Freq. A1 EUR	Freq. A1 EAS	βKöttgenb, se, P	βOkada, se, P	βMeta, se, P	% varc
Previously reported loci (Köttgen et al.)
rs1471633	1:145723739	PDZK1	A	C	0.49	0.84	0.061, 0.005, 1.5E-29	0.047, 0.023, 0.045	0.060, 0.005, 8.9E-29	0.11
rs11264341	1:155151493	TRIM46	C	T	0.59	0.31	−0.048, 0.006, 3.2E-16	−0.056, 0.019, 2.7E-03	−0.049, 0.006, 4.1E-18	0.07
rs1260326	2:27730940	GCKR	T	C	0.41	0.48	0.077, 0.006, 1.3E-44	0.052, 0.013, 1.0E-04	0.073, 0.005, 1.1E-46	0.17
rs17050272	2:121306440	INHBB	G	A	0.55	0.54	−0.037, 0.006, 1.2E-09	−0.035, 0.014, 0.014	−0.037, 0.006, 5.7E-11	0.04
rs11894371	2:148575872	ACVR2A	A	C	0.67	0.54	−0.034, 0.006, 4.5E-09	−0.004, 0.013, 0.76	−0.029, 0.005, 4.0E-08	0.03
rs6770152	3:53100214	SFMBT1	G	T	0.41	0.60	0.048, 0.006, 9.1E-18	0.022, 0.014, 0.14	0.045, 0.005, 2.0E-17	0.07
rs11722228	4:9915741	SLC2A9	C	T	0.65	0.73	0.206, 0.006, 6.9E-294	0.183, 0.014, 1.2E-36	0.309, 0.000, 2.2E-251	3.10
rs4148155	4:89054667	ABCG2	A	G	0.91	0.71	−0.221, 0.009, 3.6E-129	−0.168, 0.015, 2.3E-30	−0.226, 0.009, 2.3E-135	0.51
rs17632159	5:72431482	TMEM171	G	C	0.71	0.74	0.038, 0.006, 3.4E-10	0.039, 0.015, 9.4E-03	0.031, 0.006, 5.0E-07	0.04
rs675209	6:7102084	RREB1	T	C	0.27	0.91	0.063, 0.006, 5.3E-24	0.043, 0.026, 0.10	0.062, 0.006, 4.0E-23	0.10
rs1165213	6:25799676	SLC17A1	G	A	0.44	0.18	−0.094, 0.005, 3.8E-66	−0.070, 0.018, 8.9E-05	−0.080, 0.005, 4.5E-67	0.26
rs729761	6:43804571	VEGFA	T	G	0.27	0.13	−0.046, 0.006, 2.0E-13	−0.022, 0.020, 0.27	−0.044, 0.006, 7.6E-13	0.05
rs1178977	7:72857049	BAZ1B	A	G	0.81	0.90	0.050, 0.007, 4.4E-13	0.049, 0.022, 2.5E-03	0.050, 0.007, 5.1E-14	0.05
rs17786744	8:23777006	STC1	A	G	0.57	0.65	−0.031, 0.005, 1.6E-08	−0.014, 0.015, 0.37	−0.029, 0.005, 3.1E-08	0.03
rs2941484	8:76478768	HNF4G	C	T	0.57	0.70	0.049, 0.006, 5.4E-19	0.049, 0.013, 3.1E-04	0.049, 0.005, 7.9E-22	0.07
rs10994856	10:52645248	A1CF	G	A	0.80	0.95	0.053, 0.007, 2.3E-13	0.033, 0.032, 0.30	0.052, 0.007, 1.1E-12	0.06
rs1171617	10:61467182	SLC16A9	G	T	0.24	0.001	−0.073, 0.007, 2.4E-25	-	-	0.12
rs10897518	11:64360705	SLC22A12	C	T	0.29	0.78	0.070, 0.006, 1.2E-33	0.249, 0.018, 2.9E-45	0.074, 0.006, 4.0E-36	0.13
rs12289836	11:65436888	RELA	A	G	0.65	0.78	−0.044, 0.006, 3.5E-15	−0.068, 0.020, 6.1E-04	−0.046, 0.005, 8.9E-18	0.05
rs3741414	12:57844049	INHBC	C	T	0.81	0.93	−0.071, 0.007, 8.3E-24	−0.030, 0.025, 0.24	−0.068, 0.007, 4.1E-22	0.10
rs653178	12:112007756	ATXN2	C	T	0.47	0.003	0.036, 0.005, 2.0E-11	-	-	0.04
rs1976748	15:76160951	UBE2Q2	A	G	0.50	0.37	−0.037, 0.006, 1.3E-11	−0.028, 0.014, 0.046	−0.036, 0.005, 2.7E-12	0.05
rs6598541	15:99271135	IGF1R	A	G	0.35	0.43	0.044, 0.006, 2.7E-14	0.049, 0.014, 3.1E-04	0.044, 0.005, 3.7E-17	0.04
rs33063	16:69640217	NFAT5	A	G	0.14	0.08	0.042, 0.008, 6.1E-08	0.096, 0.026, 1.7E-04	0.046, 0.007, 9.9E-11	0.03
rs11150189	16:79734227	MAF	A	G	0.65	0.72	0.032, 0.006, 2.4E-08	0.054, 0.014, 2.1E-04	0.035, 0.005, 4.0E-11	0.03
rs7224610	17:53364788	HLF	C	A	0.43	0.14	0.038, 0.006, 4.3E-12	0.030, 0.017, 0.083	0.037, 0.005, 1.9E-12	0.04
rs9895661	17:59456589	BCAS3	C	T	0.19	0.53	−0.045, 0.008, 1.7E-09	−0.053, 0.015, 6.1E-04	−0.047, 0.007, 5.3E-12	0.04
rs164009	17:74283669	QRICH2	A	G	0.62	0.35	0.029, 0.006, 2.1E-07	0.027, 0.014, 0.065	0.028, 0.005, 4.0E-08	0.02
New loci
rs11099098	4:81169912	FGF5	G	T	0.72	0.61	0.033, 0.006, 1.4E-07	0.039, 0.014, 5.5E-03	0.034, 0.006, 2.9E-09	0.03
rs7706096	5:39994900	LINC00603	G	A	0.41	0.49	0.028, 0.005, 3.8E-07	0.029, 0.013, 0.031	0.028, 0.005, 3.5E-08	0.02
rs2858330	6:32658715	HLA-DQB1	T	C	0.49	0.25	0.026, 0.005, 1.9E-06	0.043, 0.015, 5.4E-03	0.027, 0.005, 4.2E-08	0.02
rs10813960	9:33180362	B4GALT1	C	T	0.73	0.46	0.033, 0.006, 1.9E-07	0.040, 0.013, 2.9E-03	0.035, 0.006, 2.3E-09	0.03
rs1649053	10:60321487	BICC1	T	C	0.61	0.77	−0.027, 0.006, 9.8E-07	−0.051, 0.016, 1.6E-03	−0.029, 0.005, 8.9E-09	0.02
rs11231463	11:63184455	SLC22A9	A	G	0.99	0.94	-	0.312, 0.029, 4.3E-27	-	0.80
rs7928514	11:63360114	PLA2G16	G	A	0.90	0.90	-	0.115, 0.021, 5.5E-10	-	0.17
rs641811	11:63869596	FLRT1	G	A	0.82	0.78	0.030, 0.007, 7.7E-06	0.072, 0.014, 5.9E-07	0.037, 0.006, 1.0E-09	0.02
rs11227805	11:67246757	AIP	C	T	0.79	0.88	-	0.141, 0.023, 5.5E-10	-	0.25
rs2195525	11:119235404	USP2	C	T	0.46	0.17	0.031, 0.006, 5.3E-08	0.035, 0.017, 0.040	0.031, 0.005, 6.3E-09	0.03
Independent signals							βjoint, se, P
rs10939614	4:9926613	SLC2A9	T	C	0.34	0.10	0.179, 0.008, 9.2E-118	-	-	0.90
rs4447861	4:9953940	SLC2A9	C	T	0.59	0.91	0.143, 0.009, 8.2E-57	-	-	0.50
rs12499240	4:10103890	SLC2A9	T	C	0.27	0.38	0.062, 0.007, 3.12E-19	-	-	0.07
rs4698031	4:10315921	SLC2A9	A	G	0.82	0.87	0.121, 0.009, 4.18E-18	-	-	0.38
rs2622629	4:89094064	ABCG2	T	C	0.69	0.44	−0.056, 0.006, 3.5E-22	-	-	0.08
rs575416	5:72437534	TMEM171	A	G	0.90	0.99	−0.045, 0.009, 7.7E-07	-	-	0.02
rs1171606	10:61434519	SLC16A9	G	A	0.72	0.40	−0.036, 0.006, 2.2E-08	-	-	0.03
rs2269730	11:64423831	SLC22A12	G	A	0.80	0.55	−0.050, 0.007, 7.1E-14	-	-	0.05

^aA1 is the effect allele.

^bNot identical to Köttgen et al. (4) or Okada et al. (5) P-value. This P-value is from z-score adjusted by LD-score intercept in both Okada and Köttgen. β in mg/dl.

^cVariance explained is estimated where possible using the European data set. The variance explained for the conditionally independent signals is estimated from the joint model. For all other loci, the marginal effects are used.

Cis-eQTL colocalization analysis identifies 34 candidate causal genes at 24 serum urate loci

To connect the serum urate associations with the genes they influence, we utilized publicly available expression data provided by the GTEx consortium (Methods) and performed colocalization with COLOC (27) (Supplementary Material, Fig. S7, Table 4 and Supplementary Material, Table S1). This method attempts to identify whether the causal variant is the same in both the eQTL and GWAS signal indicating a putative causal mechanism, whereby the variant alters gene expression (transcript levels) and expression influences the trait—in this case serum urate levels. This approach provides further support for the loci identified by the trans-ancestral meta-analysis by linking the serum urate signals into the biological process of gene regulation, although we note that the data in GTEx come predominantly from people of European (83%) and African (13%) ancestry meaning that we will not detect colocalization with East Asian-specific variants. For 24 of the serum urate GWAS loci strong evidence for colocalization (posterior probability of colocalization (PPC) > 0.8) was seen with 34 cis-eQTL (Table 4; Supplementary Material, Fig. S7). The 24 loci included five loci identified by inclusion of sub-genome-wide significant GWAS loci in the analysis (DHRS9, RAI14, MLXIP, IDH2 and MRPS7) (Table 2). For 15 of the previously identified Köttgen et al. (4) GWAS loci, there are colocalized cis-eQTL (PDZK1, TRIM46, INHBB, SFMBT1, BAZ1B, SLC16A9, INHBC, UBE2Q2, IGF1R, MAF, QRICH2, HLF, ABCG2, SLC17A1 and SLC22A12). Of the 10 new loci discovered as genome-wide significant in the trans-ancestral meta-analysis, colocalized eQTL were identified at four loci (HLA-DQB1, B4GALT1, RELA and FGF). The kidney-specific analysis (kidney cortex from GTEx v8 and nephrotic tubulointerstitial) revealed four additional colocalization events; FGF5, ITIH4 (at SFMBT1), TMEM151A (at RELA) and HLF (Table 4 and Supplementary Material, Fig. S8).

Supplementary Material, Table S1 presents posterior probabilities for colocalization for all GTEx tissues with a heatmap in Figure 3. Of the individual loci, it was notable that PDZK1 and MLXIP had colocalization scores largely restricted to the gut, with HLF, SLC16A9, FGF5 and TMEM151A restricted to the kidney and SFBMT1 to both the gut and kidney. The gut and kidney were noticeable for having colocalizations, but not the liver, emphasizing that genetic control of excretion is central to regulation of serum urate levels.

For the five loci (SLC2A9, ABCG2, TMEM171, SLC16A9 and SLC22A11/A12) where we identified additional independent associations, colocalization analysis was performed. This analysis revealed seven colocalization events, of which six were for genes not otherwise detected as cis-eQTLs (there was an additional eQTL for SLC16A9). For SLC22A11/A12, the signal arising from rs10897518/rs2078267 colocalized with cis-eQTLs for SLC22A11 and MEN1, and the signal arising from rs2269730 colocalized with cis-eQTLs for NRXN2 and SLC22A12 (Supplementary Material, Fig. S9). Notably the SLC22A11 and SLC22A12 eQTLs were detected in the kidney. For the ABCG2 locus, the lead GWAS variant is the established rs2231142 (p.Gln141Lys) variant; however, the independent association marked by rs2622629 colocalizes with a cis-eQTL for ABCG2.

Figure 2

Manhattan plots showing −log₁₀(P) for all SNPs of the European, East Asian and trans-ancestral GWAS ordered by chromosomal position. (A) Manhattan plot of the European GWAS. (B) Manhattan plot of the East Asian GWAS. (C) Manhattan plot of the trans-ancestral GWAS. SNPs within previously identified serum urate loci are colored light green. SNPs located within novel serum urate loci are colored orange. For the 10 new genome-wide significant loci identified by trans-ancestral meta-analysis, the closest gene to the lead SNP is indicated.

Figure 3

Heatmap of cis-eQTL colocalizations. Drawn using the pheatmap package using Supplementary Material, Table S1 data, with row and column clustering. Note the significant amount of ‘not available’ data for ITIH4, GALNTL5 and FGF5. STIMATE = TMEM110.

CoDeS3D analysis integration with GTEx and colocalization for identification of trans-eQTL

To identify candidate causal genes that are activated in trans by eQTLs at GWAS loci (i.e. trans-eQTLs), we pre-screened for SNP-gene physical connectivity using the ‘Contextualize Developmental SNPs using 3D Information’ (CoDeS3D) algorithm and then tested for colocalization with serum urate GWAS signals (Table 4). This identified 20 trans-eQTL signals that co-localized (PPC > 0.8) with 15 GWAS loci (Supplementary Material, Fig. S7). Of the 20 genes with colocalized trans-eQTL we identified, only two had evidence within the gene (P < 5 × 10⁻⁰⁴) for a signal of association with serum urate by GWAS (Supplementary Material, Fig. S10)—UTRN in the Köttgen et al. dataset (lead variant rs4896735, P = 2 × 10⁻⁰⁴) and DMD (rs1718043; P = 9 × 10⁻⁰⁵) in the Kanai et al. (6) dataset. The DMD and UTRN genes encode components of the dystrophin complex. A trans-eQTL for MAPK6 (also known as ERK3) colocalized with the serum urate GWAS signal at the IDH2 locus. Interestingly, the MAPK6 locus is associated with serum urate levels in response to allopurinol in gout by GWAS (rs62015197, P = 8 × 10⁻⁰⁷) (28).

Ten serum urate loci (SLC16A9, BAZ1B, QRICH2, UBE2Q2, INHBB, INHBC, HLF, DHRS9, MLXIP and IDH2) exhibited both cis- and trans-eQTL of which the latter three had been identified by the genome-wide colocalization analysis. At SLC16A9, the signal is different between the cis- and trans-eQTL (Supplementary Material, Fig. S7), with all of the GWAS signal present in the cis-eQTL whereas only the signal associated with the lead GWAS SNP was evident in the trans-eQTL. Differential cis- and trans-eQTL signals are reminiscent of the situation at the serum urate-associated cis- and trans-eQTL signals at the MAFTRR locus (12).

Table 2

Loci identified by genome-wide cis-eQTL colocalization

SNP	Chr:bp	Closest gene	A1a	A2	Freq. A1 EUR	Freq. A1 EAS	βKöttgena, se, P	βOkada, se, P	βMeta, se, P	% varb
Cis-eQTL co-localized loci
rs3815574	2:169963330	DHRS9	A	C	0.52	0.37	0.025, 0.005, 4.0E-06	0.013, 0.016, 0.41	0.023, 0.005, 5.3E-06	0.02
rs461660	5:34657025	RAI14	A	C	0.55	0.38	0.026, 0.006, 5.4E-06	0.007, 0.014, 0.63	0.023, 0.005, 1.4E-05	0.02
rs7953704	12:121517820	MLXIP	A	G	0.47	0.51	−0.028, 0.006, 3.5E-07	−0.003, 0.013, 0.84	−0.024, 0.005, 2.1E-06	0.02
rs8024386	15:90670526	IDH2	A	C	0.74	0.79	−0.029, 0.006, 3.9E-06	−0.018 0.020, 0.37	−0.028, 0.006, 4.3E-06	0.02
rs4788878	17:73260298	MRPS7	A	G	0.18	0.04	−0.034, 0.007, 2.4E-06	−0.000, 0.034, 0.99	−0.032, 0.007, 1.5E-05	0.02

^aA1 is the effect allele.

^bVariance explained is estimated where possible using the European data set. The variance explained for the conditionally independent signals is estimated from the joint model. For all other loci, the marginal effects are used.

Replication in Kanai et al.

While this work was being finalized, a serum urate GWAS comprising 109 029 Japanese individuals [of whom 18 519 overlapped with the Okada et al. (5) study] was published (6) allowing an opportunity to validate our findings. Seven of the 10 new loci we identified replicated (P < 0.003) in the Kanai et al. (6) study, and there was P < 0.004 evidence to support the candidacy of two (RAI14 and MLXIP) of the five loci identified by the genome-wide cis-eQTL colocalization analysis (Table 3). The nine loci included two (FGF5 and BICC1) of the total 27 genome-wide significant signals reported by Kanai et al.—of the remaining 25 loci identified by Kanai et al. (6) 17 had previously been reported by others (4,29,30) and eight were new (the gene containing the lead SNP or the flanking genes at each locus: RNF115 (rs12123298), USP23 (rs7570707), UNCX (rs4724828), TP53INP1 (rs7835379), EMX2/RAB11F1P2 (rs1886603), SBF2 (rs2220970), MPPED2/DCDC5 (rs963837), GNAS (rs6026578)).

Table 3

Validation in serum urate and gout datasets

SNP	Chr: bp	Closest gene	A1a	A2	βMeta, se, P	βKanai, se, P	ORGoutUKBB, se, Pb	ORGoutJapan, se, P	ORGoutChina, se, P
New loci
rs11099098	4: 81169912	FGF5	G	T	0.034, 0.006, 2.9E-09	0.025, 0.004, 2.8E-08	1.03, 0.019, 0.14	1.11, 0.065, 0.12	-
rs7706096	5: 39994900	LINC00603	G	A	0.028, 0.005, 3.5E-08	0.005, 0.004, 0.23	1.00, 0.017, 0.90	1.00, 0.061, 0.97	0.89, 0.098, 0.24
rs2858330	6: 32658715	HLA-DQB1	T	C	0.027, 0.005, 4.2E-08	-c	1.01, 0.017, 0.46	-d
rs10813960	9: 33180362	B4GALT1	C	T	0.035, 0.006, 2.3E-09	0.017, 0.004, 5.8E-05	1.04, 0.019, 0.070	0.90, 0.066, 0.10	1.04, 0.096, 0.68
rs1649053	10: 60321487	BICC1	T	C	−0.029, 0.005, 8.9E-09	−0.020, 0.005, 5.3E-05	0.95, 0.017, 1.6E-03	0.99, 0.095, 0.92	0.97, 0.112, 0.83
rs11231463	11: 63184455	SLC22A9	A	G	-	0.207, 0.009, 6.6E-123	1.05, 0.098, 0.61	1.39, 0.133, 0.013	1.41, 0.282, 0.22
rs7928514	11: 63360114	PLA2G16	G	A	-	0.086, 0.006, 3.8E-54	0.98, 0.027, 0.47	1.31, 0.082, 9.7E-04	-
rs641811	11: 63869596	FLRT1	G	A	0.037, 0.006, 1.0E-09	0.079, 0.004, 5.0E-70	1.07, 0.021, 1.2E-03	0.99, 0.066, 0.82	-
rs11227805	11: 67246757	AIP	C	T	-	0.085, 0.007, 2.6E-34	0.97, 0.021, 0.21	1.01, 0.101, 0.89	-
rs2195525	11: 119235404	USP2	C	T	0.031, 0.005, 6.3E-09	0.013, 0.005, 0.019	1.06, 0.017, 1.7E-03	1.12, 0.078, 0.13	-
Independent signals
rs10939614	4: 9926613	SLC2A9	T	C	-	-	1.24, 0.017, 1.7E-34	1.30, 0.103, 0.013	1.37, 0.135, 0.021
rs4447861	4: 9953940	SLC2A9	C	T	-	-	1.15, 0.017, 3.5E-15	0.78, 0.116, 0.031	0.89, 0.149, 0.42
rs12499240	4: 10103890	SLC2A9	T	C	-	-	0.95, 0.019, 7.5E-03	1.09, 0.063, 0.19	0.87, 0.100, 0.18
rs4698031	4: 10315921	SLC2A9	A	G	-	-	1.58, 0.024, 3.6E-78	1.09, 0.085, 0.33	1.45, 0.156, 0.017
rs2622629	4: 89094064	ABCG2	T	C	-	-	0.90, 0.018, 3.9E-09	0.80, 0.063, 3.3E-04	0.96, 0.093, 0.64
rs575416	5: 72437534	TMEM171	A	G	-	-	0.90, 0.027, 1.4E-04	1.27, 0.18, 0.19	-
rs1171606	10: 61434519	SLC16A9	G	A	-	-	0.95, 0.019, 3.2E-03	0.91, 0.061, 0.13	1.00, 0.096, 0.99
rs2269730	11: 64423831	SLC22A12	G	A	-	-	0.83, 0.021, 9.8E-19	0.84, 0.062, 3.5E-03	0.79, 0.096, 0.016
Cis-eQTL co-localized loci
rs3815574	2: 169963330	DHRS9	A	C	0.023, 0.005, 5.3E-06	0.000, 0.005, 0.98	1.06, 0.017, 8.1E-04	0.98, 0.067, 0.76	0.87, 0.095, 0.16
rs461660	5:34657025	RAI14	A	C	0.023, 0.005, 1.4E-05	0.013, 0.005, 3.6E-03	1.07, 0.017, 1.0E-04	0.91, 0.065, 0.17	0.97, 0.099, 0.78
rs7953704	12: 121517820	MLXIP	A	G	−0.024, 0.005, 2.1E-06	−0.017, 0.004, 3.3E-05	0.90, 0.017, 3.1E-09	0.90, 0.061, 0.091	0.90, 0.094, 0.27
rs8024386e	15:90670526	IDH2	A	C	−0.028, 0.006, 4.3E-06	−0.007, 0.006, 0.25	0.89, 0.020, 8.9E-09	0.89, 0.081, 0.16	1.10, 0.120, 0.44
rs4788878	17: 73260298	MRPS7	A	G	−0.032, 0.007, 1.5E-05	0.006, 0.010, 0.57	0.95, 0.023 0.027	1.10, 0.169, 0.57	0.75, 0.241, 0.22

^aA1 is the effect allele

^bAge- and sex-adjusted. At the independent signals data were extracted directly from GWAS summary statistics with no modeling to determine whether or not the signals were independent with gout as outcome.

^cHLA-DQ data were not reported by Kanai et al. (6)

^dThe variant was not imputed in the Japanese data set.

^eSurrogate rs7175469 was used in the Chinese gout data set—r² = 0.97, D′ = 1.0.

Testing for association with gout

To further validate the urate signals, we tested the independent signals at eight existing loci, 10 new loci with genome-wide significance in the trans-ancestral meta-analysis and five loci discovered by genome-wide colocalization (albeit at a sub-genome-wide level of significance) with eQTL (Table 3) for association with gout in European (UK Biobank) (31), Chinese (32) and Japanese (30) sample sets. The BICC1, FLRT1 and USP2 loci replicated (P ≤ 1.6 × 10⁻⁰³) in the European dataset in a directionally consistent fashion (i.e. the urate-increasing allele associated with an increased risk of gout). The SLC22A9 and PLA2G16 loci replicated (P ≤ 0.013) in the Japanese dataset also in a directionally consistent fashion. All eight additional variants identified in the European serum urate data set by conditional analysis (Table 1) were replicated (P ≤ 3.3 × 10⁻⁰³), in the European gout data set (Table 3). For the five loci identified by colocalization with eQTL, their candidacy as serum urate-controlling loci was supported by evidence for association with gout (P ≤ 0.027) in the European UK Biobank gout data set, with IDH2 and MLXIPL at a genome-wide level of significance (P < 5.0 × 10⁻⁰⁸). All had an OR for gout consistent with the direction of effect on serum urate levels. None of these five loci were associated with gout in the Chinese or Japanese sample sets.

Functional partitioning of the heritability of serum urate levels

To understand the functional categories that contribute most to the heritability of serum urate level, we used LD score regression to functionally partition the SNP heritability of the European serum urate GWAS (Fig. 4 and Supplementary Material, Fig. S11; Supplementary Material, Tables S2 and S3). Functional partitioning of serum urate SNP heritability according to cell type revealed significant enrichments in the kidney (P = 3.2 × 10⁻⁰⁸), the gastrointestinal tract (P = 5.3 × 10⁻⁵) and the liver (P = 3.4 × 10⁻⁰³). A refined analysis of 218 functional annotations, which contribute to the larger cell type groups, revealed 11 significant annotations: four histone marks in the kidney H3K27ac (P = 1.2 × 10⁻⁰⁷), H3K9Ac (P = 1.5 × 10⁻⁰⁶), H3K4me3 (P = 9.6 × 10⁻⁰⁶) and H3K4me1 (P = 2.5 × 10⁻⁰⁵) and two histone marks in the gastrointestinal tract—H3K27ac (P = 2.5 × 10⁻⁰⁶) and H3K4me1 (P = 1.6 × 10⁻⁰⁵). These histone marks are characteristic of transcriptional activation and consistent with active expression of nearby genes.

Table 4

Serum urate associated loci with colocalized GTEx eQTL

Locus	Lead GWAS variant	Colocalized eQTL gene	PPCa	Tissue(s)	Direction (allele, βSU mg/dl, β Expression, PExpression)
Cis-eQTL (genome-wide significant by trans-ancestral meta-analysis)
PDZK1	rs1471633	PDZK1	0.98	Colon - transverse, small intestine	A, 0.061, 0.58, 1.3E-11 (colon - transverse)
TRIM46	rs11264341	MUC1	0.93	Adipose subcutaneous, artery aorta, esophagus mucosa, esophagus muscularis, testis, whole blood, kidney tubulointerstitial	T, −0.048, 0.34, 7.7E-16 (esophagus—mucosa)
		GBAP1	0.98	Skin	T, −0.048, −0.35, 2.7E-08
		FAM189B	0.92	Heart atrial appendage	T, −0.048, −0.23, 1.1E-05
INHBB	rs17050272	INHBB	0.80	Lung	G, −0.037, 0.23, 2.2E-05
FGF5	rs11099098	FGF5	0.98	Kidney cortex, kidney tubulointerstitial	T, −0.038, 0.62, 1.1E-07 (kidney cortex)
SFMBT1	rs6770152	TMEM110	0.86	Adipose subcutaneous, skin	T, −0.048, 0.24, 1.1E-05 (adipose subcutaneous)
		SFMBT1	0.97	Colon transverse, kidney tubulointerstitial	T, −0.048,0.45, 1.9E-10 (colon transverse)
		ITIH4	0.94	Kidney cortex	T, −0.048, 0.36, 1.4E-05
HLA-DQB1	rs2858330	HLA-DQA2	0.82	Prostate	T, 0.026, −0.79, 2.8E-12
BAZ1B	rs1178977	MLXIPL	0.88	Adipose visceral, transformed fibroblasts	T, 0.050, −0.48, 1.2E-08 (transformed fibroblasts)
B4GALT1	rs10813960	B4GALT1	0.81	EBV transformed lymphocytes, esophagus mucosa	T, −0.033, −0.47, 1.9E-04 (EBV transformed lymphocytes)
SLC16A9	rs1171617	SLC16A9	0.86	Artery aorta, thyroid, kidney tubulointerstitial	T, 0.073, −0.28, 2.2E-06 (thyroid)
RELA	rs12289836	OVOL1-AS1	0.88	Thyroid, caudate basal ganglia, cortex	A, −0.043, 0.55, 3.7E-06 (basal ganglia)b
		TMEM151A	0.80	Kidney cortex	A, −0.043, −0.39, 0.04
INHBC	rs3741414	R3HDM2	0.84	Transformed fibroblasts	T, −0.071, 0.21, 2.2E-05
UBE2Q2	rs1976748	UBE2Q2	0.91	Dorsolateral prefrontal cortex	A, −0.037, −0.10, 8.5E-30
IGF1R	rs6598541	IGF1R	0.83	Heart left ventricle	A, 0.044, −0.33, 1.7E-07
MAF	rs11150189	MAFTRR	0.84	Colon sigmoid, pancreas	A, 0.032, 0.52, 1.2E-05 (colon sigmoid)
HLF	rs7224610	HLF	0.85	Kidney tubulointerstitial	C, 0.038, −0.75, 2.5E-15
QRICH2	rs164009	UBALD2	0.90	Esophagus muscularis, caudate basal ganglia, whole blood	A, 0.029, 0.23, 7.7E-07 (esophagus muscularis)
		PRPSAP1	0.85	Anterior cingulate cortex	A, 0.029, 0.62, 1.7E-06
Cis-eQTL (conditional association signals)
ABCG2	rs2622629	ABCG2	0.85	Esophagus mucosa	T, −0.056, −0.22, 3.3E-03
SLC17A1	rs1359231	SLC17A4	0.82	Pancreas	G, −0.044, −0.49, 1.0E-4
SLC16A9	rs1171617	SLC16A9	0.89	Artery aorta, thyroid, kidney tubulointerstitial	T, −0.061, −0.28, 2.2E-06 (thyroid)
SLC22A12	rs2078267	SLC22A11	0.82	Kidney tubulointerstitial	C, 0.074, 0.19, 8.6E-04
		MEN1	0.84	Whole blood	C, 0.074, 0.088, 4.9E-05
SLC22A12	rs2269730	SLC22A12	0.94	Kidney tubulointerstitial	G, −0.050, −0.26, 2.4E-04
		NRXN2	0.85	Transformed fibroblasts	G, −0.050, −0.29, 5.8E-11
Cis-eQTL (sub genome-wide significant by trans-ancestral meta-analysis)
DHRS9	rs3815574	DHRS9	0.83	Whole blood	A, 0.024, −0.30, 1.7E-26
RAI14	rs461660	RAI14	0.83	Thyroid	A, 0.026, 0.24, 3.5E-07
MLXIP	rs7953704	MLXIP	0.83	Small intestine	A, −0.028, 0.48, 9.6E-07
IDH2	rs8024386	IDH2	0.87	Atrial appendage	A, −0.029, 0.35, 1.1E-06
MRPS7	rs4788878	GGA3 MRPS7	0.83 0.87	Thyroid dorsolateral prefrontal cortex, colon transverse, transformed fibroblasts, pancreas	A, −0.034, 0.23, 1.9E-08 A, −0.034, 0.07, 2.5E-12 (dorsolateral prefrontal cortex)
Trans-eQTL
NFAT5	rs33063	AIF1L	0.92	Brain substantia nigra	A, 0.042, 0.39, 1.9E-05
		CACNA2D3	0.99	Basal ganglia	A, 0.042, 0.23, 4.2E-06
		STIM1	0.99	Basal ganglia	A, 0.042, 0.26, 5.6E-07
SLC16A9	rs1171617	ANKS1B	0.98	Testis	T, 0.073, 0.17, 1.8E-05
		DSCAM	0.98	Brain hypothalamus	T, 0.073, 0.51, 1.9E-05
BAZ1B	rs1178977	RNF24	0.98	Brain cortex	A, 0.050, −0.43, 1.3E-05
QRICH2	rs164009	PPP3R1	0.98	Heart left ventricle	A, 0.029, −0.18, 7.0E-06
INHBB	rs17050272	CHAC2	0.97	Basal ganglia	A, 0.037, −0.36, 1.3E-05
		ZNF804A	0.99	Brain anterior cingulate cortex	A, 0.037, 0.23, 1.1E-05
UBE2Q2	rs1976748	COL11A1	0.95	Colon transverse	A, −0.037, 0.23, 6.0E-06
HNF4G	rs2941484	CSMD2	0.99	Brain cerebellum	T, 0.049, −0.40, 5.2E-06
INHBC	rs3741414	SPIN1	0.98	Brain frontal cortex	T, −0.071, −0.22, 1.7E-05
DHRS9	rs3815574	JHDM1D	0.94	Brain hippocampus	A, 0.024, 0.32, 1.3E-05
IDH2	rs8024386	MAPK6 ZBTB20	0.86 0.89	Brain amygdala testis	A, −0.029, −0.32, 1.2E-04 A, −0.029, 0.10 1.8E-05
RREB1	rs675209	UTRN	0.99	Brain putamen basal ganglia	T, 0.063, 0.35, 2.1E-05
HLF	rs7224610	DMD	0.94	Brain nucleus accumbens basal ganglia	A, −0.038, −0.36, 1.9E-05
VEGFA	rs729761	CLPS	1.00	Brain cerebellar hemisphere	T, −0.046, −0.59, 3.8E-06
MLXIP	rs7953704	NDUFA12	0.95	Brain putamen basal ganglia	A, −0.028, −0.32, 2.0E-05
BCAS3	rs9895661	TMEM117	0.99	Prostate	T, 0.045, 0.48, 3.9E-06

^aPosterior probability of colocalization.

^bData from proxy variant rs642803.

Figure 4

Tissue-focused functional heritability enrichments. Tissue-focused functional heritability enrichments for serum urate levels. The color of each bar indicates whether heritability was depleted or enriched within a particular cell-type group (only CNS was depleted). The –log10 P-value for the enrichment in each cell-type group is on the Y-axis. These enrichments were generated using LD-Score functional partitioning of the European GWAS summary statistics.

We next sought to assess the proportion of heritability in serum urate levels attributable to regulatory variation across the GTEx tissues. We first assessed the heritability surrounding genes that are specifically expressed in certain tissue using LD-score regression (Methods). This revealed a significant enrichment of the regions near genes that are specifically expressed in the kidney (P = 9.1 × 10⁻⁰⁴) (Supplementary Material, Fig. S12), further highlighting the central role of the kidney in serum urate control. We observed a modest 1.5-fold enrichment overall for GTEx cis-eQTLs in the serum urate GWAS (SNP proportion = 0.36, proportion of h² = 0.54; enrichment P = 5.5 × 10⁻⁰⁴). This analysis was repeated using annotations derived from each tissue, with no tissue being dominantly enriched in the serum urate GWAS (Supplementary Material, Fig. S13). For this analysis, we also included the kidney cortex eQTL data from GTeX v8, but did not observe a significant enrichment, possibly because this specific dataset comprised only 73 individuals. Overall, this suggests that cis-eQTLs that contribute to serum urate heritability are distributed throughout many tissues.

Trans-ancestral functional fine mapping identifies putative causal variants

We sought to leverage both the functional enrichments and linkage disequilibrium differences between the populations to identify candidate causal variants at each locus associated with serum urate levels. To this end, we performed trans-ancestral fine mapping with PAINTOR (Probabilistic Annotation INtegraTOR) using the kidney, gastrointestinal tract and liver cell type group annotations as functional priors. Restricting to variants with posterior probabilities >0.8, when analyzing only the European GWAS the 90% causal credible sets had on average 129 SNPs. With the addition of the East Asian GWAS data, the set size reduced to an average of 56 SNPs and functional annotations reduced the average credible set size to 41. Of the 36 loci used in this analysis (the 28 reported by Köttgen et al. (4) and the 10 new genome-wide significant loci reported here, excluding RELA and HLA-DQB1), 14 loci had seven or fewer causal variants in their 80% causal credible set (Table 5 and Supplementary Material, Table S4). The combination of both the functional annotations and East Asian GWAS data significantly improves our ability to identify the causal variants for loci associated with serum urate levels.

SLC2A9 is a complex locus with a very strong effect on serum urate levels and multiple independent genetic effects (33,34). A subset of the lead urate SNPs at SLC2A9 with PAINTOR posterior probabilities of 1.0 overlap putative regulatory elements (Table 5 and Supplementary Material, Table S5). One of these urate-associated variants at SLC2A9, rs11723382, is also among the maximally associated cis-eQTL variants for RP11-448G15.1 (transformed lymphocytes) (RP11-448G15.1 is a lncRNA located within the second intron of SLC2A9) and disrupts two predicted motifs Hmx and Nkx2 (Supplementary Material, Fig. S1 and Table 5) (35). This eQTL was not identified in our COLOC analysis; however, visual inspection of the RP11-448G15.1 eQTL and SLC2A9 GWAS signal indicates that the signals coincide and suggests RP11-448G15.1 expression is likely important for serum urate control.

At SLC22A12/NRXN2, four of the seven putative causal variants (Table 5) are in LD (R² > 0.6) with the maximal trans-eQTL variant for RNF169 identified by CoDeS3D. Visual inspection of the RNF169 trans-eQTL and the serum urate signal at the SLC22A12/NRXN2 locus indicates that these signals overlap (Supplementary Material, Fig. S14). rs2277311, an intronic variant located within NRXN2, is the most likely candidate of these variants to have regulatory function. rs2277311 has promoter, enhancer and DNase signatures and the urate-decreasing A-allele disrupts a predicted HiC1 motif (Table 5 and Supplementary Material, Table S5) (35).

Six loci (RREB1, INHBC, HLF, UBE2Q2, SFMBT1 and HNF4G) with PAINTOR candidate causal SNPs (PP > 0.8) also have colocalized eQTL (Tables 4 and 5). These loci represent good candidates for follow-up analyses of regulatory function (11,12). The lead urate variant at the HLF locus, rs7224610 (PAINTOR posterior probability = 0.92) is intronic, has enhancer signatures, is bound by multiple transcription factors including POL2 (Supplementary Material, Table S4) and is amongst the maximally associated trans-eQTL variants for DMD (encodes dystrophin). rs675209 at RREB1 is the maximal trans-eQTL variant for UTRN (encodes utrophin), overlaps enhancer signatures in six tissues and alters 8 transcription factor binding motifs (Table 5 and Supplementary Material, Table S5). The variants at HNF4G, SFMBT1, UBE2Q2 and INHBC do not overlap putative regulatory elements (Supplementary Material, Table S5). Although rs13264750 (HNF4G), rs2115779 (SFBMT1) and rs9870898 (upstream of SFMBT1) are predicted to change 8 binding motifs including HNF1 (Table 5).

Table 5

Putative credible causal SNP set in the combined European – East Asian analysis identified with PAINTOR

Locusa	Chr:pos	SNP ID	Posterior probb	Proportion of total PP	European Z-scorec	East Asian Z-score	Motifs changed (Haploreg)
SLC2A9	4:9915741	rs11722228	1.000	0.200	−36.73	−12.67	None
SLC2A9	4:9946095	rs4697701	1.000	0.200	59.52	9.52	E2a, Mxi1
SLC2A9	4:9954660	rs11723382	1.000	0.200	38.02	7.24	Hmx, Nkx2
SLC2A9	4:9981997	rs13145758	1.000	0.200	−57.11	−2.66	Egr1, GCNF, HNF4
SLC2A9	4:9982330	rs13125646	1.000	0.200	−49.28	−2.66	None
ABCG2	4:88917735	rs17013705	1.000	0.143	4.44	5.10	Irf, TATA, TCF12
ABCG2	4:88944511	rs2725227	1.000	0.143	−7.22	−2.52	Gfi1, TCF11, MafG
ABCG2	4:88960528	rs2725217	1.000	0.143	−17.56	−9.73	CDP7, Dbx1, HNF1, Mef2, Pouf1, TATA
ABCG2	4:88973427	rs2725210	1.000	0.143	−15.27	−4.69	Fox, FoxA, FoxC1, FoxJ2, FoxF2, FoxK1, PLZF
ABCG2	4:88999222	rs2728126	1.000	0.143	−15.65	−8.72	EWRS1, FLI1
ABCG2	4:89052323	rs2231142	1.000	0.143	−24.26	−11.43	GR, Irf
ABCG2	4:89098731	rs9631715	1.000	0.143	15.34	7.24	GATA, SREBP
SLC17A1	6:25785295	rs6909187	1.000	0.250	15.65	−0.77	FoxC1, HDAC2, HMG-IY, Pou2F2
SLC17A1	6:25786993	rs3799344	1.000	0.250	15.36	2.10	Eomes, Pax6, TBX5
UBE2Q2	15:76194286	rs335685	1.000	0.348	5.25	0.13	CEBPB, DMRT1, FoxA, Nanog, Nkx6, Pou2F2, Pou3F4, STAT, TATA
SFMBT1	3:53092375	rs9870898	0.996	0.331	5.15	0.63	AhR, GR, HES1, HNF1, Pax4
SFMBT1	3:53026384	rs2564938	0.996	0.331	−6.86	−1.46	None
SFMBT1	3:53026714	rs2115779	0.996	0.331	5.38	1.63	Hsf, Ptflb
SLC22A12	11:64333296	rs1783811	0.994	0.143	7.91	5.00	Mef2, TAL1, ZID
UBE2Q2	15:76160951	rs1976748	0.991	0.348	−6.78	−2.00	Arid3a, Sox, TCF4
AIP	11:67246757	rs11227805	0.991	0.953	NA	6.22	None
SLC22A12	11:64358241	rs11602903	0.990	0.143	11.16	13.45	BAF155
SLC22A12	11:64387932	rs2277311	0.990	0.143	10.79	14.29	Hic1
SLC22A12	11:64338228	rs11231822	0.990	0.143	12.88	−2.62	CTCF, ERalphaA, Lmo2, Nanog
SLC22A12	11:64419217	rs502571	0.990	0.143	−8.41	−10.70	Mrg1, Hoxa9, TAL1
SLC22A12	11:64474752	rs2957564	0.990	0.143	−8.45	−10.94	GR, LUN1
SLC22A12	11:64622502	rs2007521	0.990	0.143	−1.50	−6.90	AP1, CTCF, Ets
SLC17A1	6:25798932	rs1165215	0.986	0.250	−16.73	−3.92	None
SLC17A1	6:26125342	rs129128	0.986	0.250	5.24	−0.09	None
SLC22A9	11:63170736	rs7925182	0.981	0.330	−1.26	−0.94	Barx1
USP2	11:119235404	rs2195525	0.978	0.881	5.46	2.06	AP1, Foxa, STAT
RREB1	6:7102084	rs675209	0.964	0.581	10.13	1.63	CCNT2, Ets, MZF1, NRSF, STAT, VDR, Zfp281, Zfp740
HNF4G	8:76401359	rs13264750	0.964	0.252	0.44	−2.78	Pou3f2
HLF	17:53364788	rs7224610	0.920	0.754	6.95	1.74	None
GCKR	2:27730940	rs1260326	0.868	0.398	14.00	3.89	NRSF
INHBC	12:57807114	rs540730	0.851	0.192	−9.44	0.75	None
SLC22A9	11:63859120	rs11231454	0.810	0.273	−0.31	10.01	GATA
SLC22A9	11:63171309	rs12281229	0.810	0.273	−1.01	−3.09	Cdx, Dbx1, Fox, FoxA, fOXc1, FoxD3, FoxF1, FoxI1 FoxJ1, Foxj2, FoxK1, FoxL1, Foxo, FoxP1, HDAC2, HNF1, Hlx1, HoxD8, Mef2, NF-Y, Ncx, Pbx-1, Pbx3, TATA
PLAG2G16	11:63360114	rs7928514	0.751	0.736	NA	5.51	LBP-1
BCAS3	17:59456589	rs9895661	0.622	0.555	−6.04	−3.44	AhR_1, Eomes, NRSF, Pax-8 TBX5
INHBB	2:121310269	rs6706968	0.610	0.555	5.97	2.59	Myb, SP1
PDZK1	1:145723739	rs1471633	0.711	0.380	11.31	2.01	Cdx2, Hoxa9, Hoxc10, Hoxc9, Irf, PRDM1, Pax-4, VDR
ATXN2	12:112007756	rs653178	0.621	0.307	6.72	NA	Esr2
UBE2Q2	15:76232422	rs10444856	0.696	0.243	−6.11	−2.02	None
HNF4G	8:76401202	rs13264686	0.734	0.192	−1.91	−2.72	None
HNF4G	8:76401202	rs16939082	0.733	0.192	−1.91	−2.62	CEBPB,NF-AT1,STAT
HNF4G	8:76401202	rs2941456	0.733	0.192	−4.74	1.48	Cdx2, Foxp1, Hoxa9, Hoxb13, Hoxb9, Mrg1::Hoxa9, Sox_5, TATA
MAF	16:79800762	rs7190102	0.660	0.177	−1.48	−0.40	ELF1, HDAC2, PU.1, Pou2f2, SPIB, p300
FLRT1	11:63859120	rs883798	0.755	0.153	0.94	6.14	GATA,HDAC2,SRF
INHBC	12:57807114	rs3923885	0.625	0.141	3.25	1.84	BDP1,Smad4
INHBC	12:57807114	rs12313306	0.610	0.137	9.76	0.57	ERalpha-a, Ik-1, SP1, ZBRK1
INHBC	12:57807114	rs7309659	0.607	0.137	4.76	1.76	GATA
INHBC	12:57807114	rs12315434	0.607	0.137	9.54	0.45	GATA, HDAC2,TCF12

^aThe top panel comprises the 38 variants with PP > 0.80; the bottom panel comprises 16 other variants in the 80 variants with the highest locus-specific proportion of PP (out of 31 215 analyzed by PAINTOR) and PP > 0.60.

^bPAINTORv3 allows for more than one causal variant meaning that the total probability will be greater than one.

^cZ-scores are reported because effect sizes are not available for imputed variants.

SLC22A9

SLC22A9 encodes organic anion transporter 7 (OAT7). OAT7, expressed only in the liver, is a relatively poorly characterized member of the OAT family (36) that includes urate secretory transporters OAT1–3 and the urate reuptake transporter OAT4 (encoded by SLC22A11) (37). RT-PCR screening of human cell lines indicated expression in HepG2 cells (Fig. 5). OAT7 exhibited modest uricosuric-sensitive urate uptake when expressed in Xenopus oocytes (Fig. 5). Pre-injection of oocytes with butyrate, but not other anions (data not shown), led to a modest trans-activation of urate transport, consistent with urate-butyrate exchange.

Figure 5

Expression analysis and functional expression of SLC22A9 (OAT7). (A) RT-PCR of SLC22A9/OAT7 expression in the human PTC-05 proximal tubular cell line, HEK-293 T cells and HepG2 hepatic cells. All three cell lines are positive for GAPDH but SLC22A9/OAT7 is unique to HepG2. (B) OAT7 is a weak urate transporter. Xenopus oocytes were microinjected with water (control cells) or cRNA for OAT1 or OAT7. OAT7-expressing cells have a very modest urate transport activity that is increased by prior microinjection with butyrate, to ‘trans-activate’ urate-butyrate exchange. This transport activity is inhibited by the uricosurics tranilast and benzbromarone, each at a concentration of 100 μM; DMSO, the diluent for tranilast and benzbromarone, has no effect on urate transport. *Refers to P < 0.001 compared to OAT7-expressing cells without butyrate pre-injection and water control cells. Data shown are from a single representative experiment.

Discussion

Identification of new loci associated with serum urate

The new loci identified here as associated with serum urate levels can be ranked according to the strength of genetic evidence according to two criteria; a genome-wide significant association with serum urate levels, replication in gout and/or replication in a Japanese serum urate GWAS (6). In addition to the strength of the genetic evidence, eight of these loci were co-localized with at least one eQTL signal, which identifies a putative causal gene and provides further evidence of a genuine association with serum urate levels. Of the 15 novel loci, eight (FGF5, BICC1, PLA2G16, B4GALT1, SLC22A9, AIP, FLRT1 and USP2) were genome-wide significant and replicated in gout (Table 1) or the Kanai et al. urate dataset (6). The HLA-DQB1 locus was genome-wide significant and a putative causal gene HLA-DQA2 was identified. The DHRS9, MLXIP, MRPS7, RAI14 and IDH2 loci were only of suggestive association in the trans-ancestral meta-analysis but the colocalization analysis provided strong evidence that they participate in a causal pathway. Of these five loci, we were able to replicate the association at MLXIP and RAI14 in both Kanai et al.’s (6) serum urate analysis and gout in the UK Biobank dataset, and for IDH2, MRPS7 and DHRS9, there was supporting evidence from association with gout in the UK Biobank dataset. Overall, the evidence that these five loci, identified solely by colocalization of GWAS signal with an eQTL signal, have a true association with serum urate is strong and provide empirical support for our genome-wide co-localization approach using sub-genome wide significant GWAS signals. Overall, we identified 14 novel loci that we are confident are unlikely to represent false positive associations (FGF5, B4GALT1, PLA2G16, SLC22A9, FLRT1, USP2, BICC1, DHRS9, RAI14, IDH2, MLXIP, AIP, MRPS7 and HLA-DQB1). The remaining locus LINC00603 was identified only as genome-wide significant in the trans-ancestral meta-analysis.

A total of seven loci (three new, one independent signal, three previously reported) are concentrated in a 4 Mb segment of Chr 11 (63.2–67.2 Mb). In the previous Okada et al. and Kanai et al. East Asian and Japanese GWAS (5,6) these loci were reported as a single locus. Köttgen et al. (4) reported three loci in this region (SLC22A11, SLC22A12 and OVOL1). The Chr11 region is clearly of importance for serum urate control, and there are more genome-wide associated loci in East Asian populations than in Europeans. At the RELA locus, the causal variants in the East Asian population are not the same as in the European population (although we note that the Okada et al. (5) East Asian RELA signal is based entirely on imputed SNPs). Notably, the effect sizes of SLC22A9 and SLC22A12 (change in urate of 0.31 and 0.25 mg/dl per allele, respectively) are larger in East Asian populations than SLC2A9 and ABCG2 (0.18 and 0.17 mg/dl, respectively). In comparison the effect sizes in Europeans for SLC2A9, ABCG2 and SLC22A12 are 0.21, 0.22 and 0.07 mg/dl, respectively, [note that the lead SLC22A9 SNP rs11231463 is uncommon in Europeans (1.1%)]. In Europeans, SLC22A12 is the seventh strongest signal in serum urate after SLC2A9, ABCG2, GCKR, SLC17A1, SLC16A9 and INHBC. At the SLC22A12 locus, only the conditional analysis detected cis-eQTLs to each of SLC22A11 (encoding OAT4) and SLC22A12 (encoding URAT1), in additional to cis-eQTLs to other genes in the region (MEN1 and NRXN2). There is extensive signal in the SLC22A12 region suggesting that (as at other loci such as SLC2A9, ABCG2 and SLC17A1) the presence of multiple signals obscured the initial cis-eQTL analysis.

Recently, a separate trans-ancestral meta-analysis of the Köttgen et al. (4) and a new serum urate GWAS of 121 745 Japanese individuals (that encompassed all the individuals in the Kanai et al. (6) study) was published (38). This study discovered 59 loci, of which 22 are newly reported beyond those reported in the Köttgen et al. and Kanai et al. studies (4,6). Of the 22, three overlapped with the 15 newly identified loci in the study reported here (HLA-DQB1, B4GALT1, USP2). Both the Kanai et al. (6) and Nakatochi et al. (38) studies reported this segment of Chr11 as a single locus. Here, we dissected the Chr11 63.2–64.4 Mb segment and identified four loci in this region, including SLC22A9 (encoding OAT7), the locus with the largest effect size on serum urate in the Japanese population (Table 1). Comparing to the very recent Tin et al. (39) serum urate trans-ancestral GWAS comprising 457 690 individuals, we identified one novel locus (MRSP7 from the genome-wide cis-eQTL colocalization analysis).

Assigning causality to reported GWAS loci

We identified 54 genes with strong evidence for colocalization with a serum urate association signal (34 from the cis-eQTL analysis and 20 from the trans-eQTL analysis). Candidate causal genes at eight loci deserve brief mention (in addition to those discussed in more detail later). First, MUC1 encodes mucin-1 (CD227), a membrane protein with excessive O-glycosylation in the extracellular domain that protects from pathogens. Mutations in MUC1 cause autosomal dominant tubulointerstitial kidney disease (40), suggesting that regulation of this gene could influence serum urate levels via an effect on the structure and function of the kidney tubule. Second, IGF1R encodes the insulin-like growth factor-1 receptor, with the eQTL implicating IGF-1 signaling and resultant anabolic processes in urate control. Third, SLC16A9 encodes mono-carboxylate transporter 9 and the urate GWAS signal is also associated with DL-carnitine and propionyl-L-carnitine levels, which are both strongly associated with serum urate levels (41). Kidneys reabsorb carnitine from the urinary filtrate by a sodium-dependent transport mechanism (42), possibly influencing urate levels indirectly as a result of the secondary sodium dependency of urate transport (37). Fourth, B4GALT1 encodes β-1,4-galactosyltransferase 1, a Golgi apparatus membrane-bound glycoprotein. This implicates sugar modification of proteins (e.g. urate transporters) in serum urate control, either by regulating their level of expression and/or activity. Fifth, PRPSAP1 has a cis-eQTL at the QRICH2 locus—PRPSAP1 (encoding phosphoribosyl pyrophosphate synthetase-associated protein 1) is a strong candidate gene. As a negative regulator of phosphoribosyl pyrophosphate synthetase that catalyzes the formation of phosphoribosyl pyrophosphate from ATP and ribose-5-phosphate in the purine salvage pathway decreased expression of PRPSAP1 would be predicted to contribute to increased urate levels. However, our data are not consistent with this hypothesis—rs164009_A associated with increased PRPSAP1 expression and increased urate levels. Sixth, a very strong colocalized trans-eQTL for CHAC2 was identified at INHBB. CHAC2 is a y-glutamyl cyclotransferase involved in glutathione homeostasis (43,44). Proximal tubule cells contain high levels of glutathione, which is transported in and out of the kidney via OAT1/3, MRP2/4 and OAT10 (45,46). Specifically, glutathione serves as a counter ion for urate reabsorption via OAT10, releasing glutathione into the lumen (46). Thus, it could be predicted that changes in CHAC2 expression would disrupt glutathione homeostasis altering urate secretion/reabsorption in the kidney. Seventh, there was evidence for a regulatory effect at ABCG2, which, in certain genotype combinations, may amplify the effects of the established causal variant at that locus (p.Gln141Lys (rs2231142)) (13). We did not find evidence of this at GCKR where p.Leu446Pro (rs1260326) is the maximally associated variant (Supplementary Material, Fig. S5). Finally, the extensive linkage disequilibrium at the SLC17A1 locus, for which there is also evidence for a role of missense variants (47), complicates deconvolution of the signal. However, the conditional analysis did detect a cis-eQTL for SLC17A4.

The dystrophin complex

Of the 20 spatially supported trans-eQTL that colocalize with European GWAS serum urate signals, two genes, DMD and UTRN (trans-eQTL at HLF and RREB, respectively), also have serum urate association signals in cis. There was a sub-genome-wide signal of association at the UTRN locus in the European serum urate GWAS data (rs4896735; P = 2.0 × 10⁻⁰⁴) (4) and a similar signal has been reported in an Indian serum urate GWAS study (rs12206002; P < 10⁻⁴; not in LD with rs4896735) (48). DMD associated with serum urate levels in the Japanese serum urate GWAS sample set (rs1718043; P = 8.8 × 10⁻⁰⁵) (6). UTRN and DMD are components of the dystrophin complex and the urate-raising alleles at these trans-eQTL increase expression of UTRN and DMD (Table 2). The canonical function of the dystrophin complex is well defined from its role in Duchenne muscular dystrophy and is crucial for stabilization of the plasma membrane in muscle cells (49). However, syntrophins within the dystrophin complex also act as scaffolding for transporters (e.g. ABCA1 (50)) and ion channels via PDZ domains, reminiscent of the PDZK1 interaction with urate transporters (37). Isoforms of the proteins within the dystrophin complex have segment-specific distribution in the mouse nephron (51); thus, it is possible that expression changes in the components of this complex in the kidney could alter the function of renal transporters that influence serum urate levels.

OAT7

An East Asian-specific genome-wide significant signal near the gene encoding OAT7, SLC22A9, was confirmed. Ideally, we would have performed a colocalization analysis to assess whether this genetic association may be influencing the expression of SLC22A9. However, because SLC22A9 is specifically expressed in the liver and brain and no East Asian eQTL are currently available for those tissues, this could not be performed. In lieu of providing genetic evidence that this association influences the expression of SLC22A9, we sought to evaluate whether OAT7 transported urate. Our data suggest that OAT7 is a very weak urate transporter in the presence of the various anions tested as exchangers (glutarate, α-ketoglutarate, butyrate, β-hydroxybutyrate). It is possible that OAT7 may function as a more efficient urate transporter in the presence of the appropriate (as yet unidentified) exchanging anion. OAT7 is a hepatic transport protein that exchanges, for the short-chain fatty acid butyrate, sulphyl conjugates, xenobiotics and steroid hormones and is not inhibited by established inhibitors and substrates of other organic anion transporters such as probenecid, paraaminohippurate, non-steroidal anti-inflammatory drugs and diuretics (36). We found that urate transport mediated by OAT7 is inhibited by the uricosuric drugs benzbromarone and tranilast, which inhibit multiple other urate transporters (52). Three uncommon missense variants that influence the ability of OAT7 to transport pravastatin by either causing the protein to be retained intracellularly or reducing protein levels at the plasma membrane have been reported (53), all at a frequency <1% in East Asian. HNF4α plays a key role in the transactivation of the SLC22A9 promoter (53), an interesting observation given that HNF4α is also required for expression of the gene encoding the urate transportosome-stabilizing molecule PDZK1 in the liver (11), and is implicated in control of serum urate levels via the MAFTRR locus (12). Very recently, initially implicated by GWAS, a functional variant of HNF4α (p.Thr139Ile) has been experimentally demonstrated to play a role in serum urate levels via transactivation of the promoter of ABCG2 (39).

Colocalization analysis assigns causation to variants at MLXIPL and MLXIP

We identified the paralogs MLXIPL and MLXIP as the putative causal genes at the BAZ1B and MLXIP loci, respectively. These genes encode the ChREBP and MondoA proteins, which are glucose-sensitive transcription factors involved in energy metabolism—including glycolytic targets and glycolysis (54–56). These proteins form heterodimers with the Mlx protein, and both of these proteins are activated by high levels of intracellular glucose-6-phosphate—a product of the first step of the glycolysis and pentose phosphate pathways. Increased activity of the pentose phosphate pathway leads to the production of ribose-5-phosphate, thus stimulating de novo purine nucleotide synthesis. The resulting nucleotides are ultimately catabolized into urate if they are not otherwise utilized. In Drosophila at least, the ChREBP/Mondo-Mlx complex is responsible for the majority of transcriptional changes that result from glucose consumption, including the pentose phosphate pathway (54). The colocalization results reveal that the serum urate-increasing variants at both loci decrease expression of MLXIPL and MLXIP. Taken together, this suggests a possible mechanism whereby the decreased basal expression of ChREBP and MondoA results in increased activity of the pentose phosphate pathway and therefore higher levels of serum urate.

MRPS7 and IDH2 and mitochondrial function

MRPS7 is putatively involved in serum urate control via mitochondrial processes. Of relevance, reduced relative mitochondrial DNA copy number is associated with gout (57). The association signal at the MRPS7 locus colocalized with gene expression of MRPS7 and GGA3. MRPS7 encodes the mitochondrial ribosomal protein S7, which is required for the assembly of the small ribosomal subunit of the mitochondria. A whole exome study revealed that a non-synonymous mutation in MRPS7 (p.Met184Val), which destabilizes the protein and reduces expression, results in impaired mitochondrial protein synthesis and impaired mitochondrial function (58). The patients in this study presented with congenital sensorineural and significant hepatic and renal impairment, consistent with a role for reduced MRPS7 activity in renal function. Our findings show that the urate-increasing G-allele decreases the expression of MRPS7 (Table 4), consistent with the hypothesis generated by the p.Met184Val phenotype. Also implicating mitochondrial function is IDH2, which encodes isocitrate dehydrogenase that catalyzes the decarboxylation of isocitrate to 2-oxyglutarate in the citric acid cycle. The urate-increasing allele associates with reduced expression of IDH2. Somatic mutations in IDH2 are implicated in a range of diseases including cancers such as glioma and acute myeloid leukemia (where an inhibitor is in phase III clinical trial (59)) and the tumor syndromes Ollier disease and Maffucci syndrome (60). Understanding the molecular mechanism of urate control by the MRPS7 and IDH2 loci locus could lead to insights into the mitochondrial processes that influence serum urate levels.

Trans-ancestral functional fine mapping identifies putative causal variants

To connect GWAS loci where we identified candidate causal genes to an underlying causal variant, we performed trans-ancestral fine mapping with PAINTOR using the kidney, gastrointestinal tract and liver cell type group annotations as functional priors. We identified six loci (RREB1, INHBC, HLF, UBE2Q2, SFMBT1 and HNF4G) that had colocalized eQTL and contained SNPs with high posterior probabilities of causality (>0.8). Two additional loci SLC2A9 and SLC22A12 also contained SNPs with high posterior probabilities of causality (>0.8) that were cis- and trans-eQTL for RP11-448G15.1 and RNF169, respectively. Seven loci (AIP, USP2, RREB1, HLF, PLA2G16, BACS3 and INHBB) had single variants with posterior probabilities >0.6 and that explained over half of the total posterior probability at their respective loci. Many of these SNPs overlapped annotated regulatory regions of the genome (Table 5). These candidate causal variants and genes provide a starting point for understanding how these variants alter serum urate levels. The power of this approach is illustrated in our prior work on the PDZK1 locus (11). In Ref. (11), we experimentally confirmed that PDZK1 was the causal gene, with rs1967017 (one of the two candidate causal variants identified with posterior probabilities >0.25 (Supplementary Material, Table S1) being a highly likely causal variant via altering a binding site for hepatocyte nuclear factor 4α. We have also applied a similar approach to the MAF locus (12). MAF is a complex locus with population-specific signals, and for one of these signals, we experimentally demonstrated that the effect on urate arises from one of two SNPs within a kidney specific enhancer that is co-expressed with MAF and HNF4A in the developing proximal tubule. This study also identified colocalized eQTL for two long intergenic non-coding RNAs MAFTRR and LINC01229 that regulate MAF expression in cis, and other genes implicated in urate metabolism in trans (12). These studies highlight the power of initially combining colocalization analyses and fine mapping using prior information to determine the molecular mechanisms that underlie GWAS signals.

In conclusion, we have identified 15 new GWAS signals associated with serum urate levels. By cis-eQTL colocalization we identified 34 candidate causal genes and by trans-eQTL analysis we implicated a further 20 genes in the molecular control of serum urate levels. Highlighted insights into molecular mechanisms come from identification of the protein encoded by SLC22A9 (OAT7) to be a urate transporter, the implication of mitochondrial function via MRPS7, the identification of MLXIP (alongside the already identified MLXIPL) and intriguing data genetically implicating the dystrophin complex in control of serum urate levels.

Materials and Methods

Data preparation and quality control

Summary statistics from the Global Urate Genetics Consortium (GUGC) meta-analysis of GWAS data consisting of 110 238 individuals of European ancestry (4) (http://metabolomics.helmholtz-muenchen.de/gugc/), and a meta-analysis consisting of 21 417 individuals of East Asian ancestry (5) were utilized. For both datasets, the following quality control procedure was followed. Firstly, we removed any SNPs that were not present in the Phase 3 release of the 1000 Genomes for the representative populations (EUR and EAS), or where the alleles were not identical between this summary data and the 1000 Genomes (e.g. the alleles were G/T in the GUGC meta-analysis and T/A in the 1000 Genomes dataset) (61). The effective sample size for each SNP was calculated using the Genome-wide Complex Trait Analysis (GCTA, v1.25.2) toolkit (62) and SNPs with effective sample sizes >2 standard deviations from the mean were excluded. Finally, SNPs with a minor allele frequency (MAF) of less than 0.01 were excluded.

Trans-ancestral meta-analysis

ImpG (v1.0) was used to impute Z-scores into the European and East Asian summary statistics. For the reference haplotypes, the Phase 3 release of the 1000 Genomes project was used (61), and only bi-allelic SNP markers having a minor allele frequency greater than 0.01 in the relevant population were included. All imputed markers with a predicted R² of less than 0.8 were removed. Meta-analysis was performed by summing the Z-scores and weighting by sample size. We used the sample sizes that were reported in the original GWAS summary statistics. For the imputed SNPs, the sample size was estimated as the median of the sample size of SNPs that were not removed in the filtering step above. To provide an adjustment for inflated test statistics, the LD-score intercept in the original summary statistics files was calculated using LD-score regression (63). This intercept adjusts the test statistics for confounding, such as cryptic relatedness, but in contrast to genomic control will not remove inflation caused by a true polygenic signal.

Independent regions were identified using the following protocol. Firstly, SNPs that were genome-wide significant (P < 5 × 10⁻⁰⁸) were padded 50 kb either side of the SNP position, and all overlapping regions were clumped together. Secondly, the maximal R² > 0.6 for the most significant SNP in each of these regions was calculated for each population. Finally, the maximal regions from the P-value clumping and LD approach were created, and any overlapping regions were merged. SNPs that were not present in both datasets were also analyzed, and for those SNPs, the LD was only calculated in the relevant population. Based on their proximity to stronger signals four loci (Chr4/rs114188639/CLNK1, Chr8/rs2927238/HNF4G, Chr11/rs641811/FLRT1, Chr11/rs117595559/VPS51) were visually examined by LocusZoom and subjected to conditional analysis—of these only rs641811/FLRT1 was concluded to be independent of the nearby signal. For all lead SNPs, the meta-analysis effect estimate was calculated using the inverse variance method, and when there was no effect estimate, it was estimated from the Z-score using the following Eq. (1).

$$\begin{equation} {S}_x=\!\!\!\!{\displaystyle \begin{array}{c}{\hat{\beta}}_x={Z}_x{S}_x\\ 1/\sqrt{2p\left(1-p\right)\left(\hat{n}+{z}_x^2\right)}\\{}\hat{n}= median(n)\end{array}} \end{equation}$$

(1)

where Z_x is the Z-score for the xth SNP, and S_x is an estimate of the standard error that requires only the Z-score, allele frequency (p) and sample size (n).

Conditional analysis

A conditional and joint analysis of the European summary statistics for all genome-wide significant regions identified by meta-analysis was done. This was not done on the East Asian summary statistics owing to the lack of availability of both an LD matrix and a reference haplotype set of sufficient size. The genotypic data from the UK Biobank was used as the reference for the LD, and to improve computational efficiency, only a random 15% (22872) of samples from the first release were included. To prevent any LD differences between 1000 Genomes and UK Biobank from influencing our conditional analysis, we restricted these analyses only to SNPs that were provided in the Köttgen et al. summary statistics. Because the GCTA-COJO module (64) was not designed to utilize dosage matrices, we performed this analysis using our own software, Correlation-based conditional analysis (COCO: https://github.com/theboocock/coco), based on the methods presented in GCTA-COJO (64), was designed to perform conditional and joint analysis from summary statistics with some minor alterations to use LD correlation matrices as input. To discover conditional associations, the coco pipeline implemented a forward stepwise selection using a residual-based regression. First, SNPs were ranked on marginal test statistics, then the top SNP was selected and the result of extracting the residuals from this model and performing a regression with every other SNP was estimated. These test statistics were then ranked. If the new top SNP passed the P-value threshold it was added to a joint model with the other selected SNP, which was used as the new model for residual extraction. This process was then repeated until no SNPs passed the significance threshold. In practice, we restricted the maximum number of selected SNPs at a locus to five [there is evidence for multiple signals at SLC2A9 (10)], and we did not consider any pairs of SNPs having an R² > 0.9. To ensure that the method was working correctly, simple phenotypes were simulated and it was verified that COCO yields almost identical results to the lm function in the R programming language.

A mathematical explanation of the method is given as follows. We assume we have mean centered genotypes in a matrix X. To perform GWAS, we generate a marginal statistic for each variant individually [Eq. (2)].

$$\begin{equation} {\hat{\beta}}_{gwas}={\left(\mathit{\operatorname{diag}}\left(X\hbox{'}X\right)\right)}^{-1}x\hbox{'}y \end{equation}$$

(2)

Using substitution into the ordinary least squares equations, we can convert these marginal effects into joint effects and also calculate the standard error [Eq. (3)].

$$\begin{equation} {\displaystyle \begin{array}{ccc}{\hat{\beta}}_{jo\mathit{\operatorname{int}}} & =& {\left(X\hbox{'}X\right)}^{-1}\mathit{\operatorname{diag}}\left(X\hbox{'}X\right){\hat{\beta}}_{gwas}\\[6pt] {}\mathit{\operatorname{var}}\left(\hat{\beta}\right) & =& {\sigma}_j^2{\left(X\hbox{'}X\right)}^{-1}\\[8pt] {}{\hat{\sigma}}_j^2 & =& \frac{y\hbox{'}y-{\hat{\beta}}_{jo\mathit{\operatorname{int}}}\mathit{\operatorname{diag}}\left(X\hbox{'}X\right){\hat{\beta}}_{gwas}}{n-N-1}\end{array}} \end{equation}$$

(3)

where N is equal to the number of SNPs in the joint model. n is equal to the effective sample size [Eq. (4)], described in more detail in Ref. (64).

$$\begin{equation} n= Median\left(\frac{Var(Y)- Var\left({X}_{,j}\right)\cdot{\beta}_j^2}{Var\left({X}_{,j}\right)\cdot Var\left({\beta}_j\right)}\right) \end{equation}$$

(4)

The calculation for the variance is given later [Eq. (9)]. The genotypic variance for the jth SNP Var(X_,j) is estimated from the UK Biobank dosages. β_j is the effect estimate, and Var(B_j) is the square of the standard error from the GWAS summary statistics.

Finally, we can approximate a regression of the residuals from a joint model [Eq. (5)].

$$\begin{equation} {\displaystyle \begin{array}{l}{\hat{\beta}}_{resid1}={\hat{\beta}}_{gwas1}-{\left(X{\hbox{'}}_1{X}_1\right)}^{-1}\left(X{\hbox{'}}_1{X}_2\right){\left(X{\hbox{'}}_2{X}_2\right)}^{-1}\mathit{\operatorname{diag}}\left(X{\hbox{'}}_2{X}_2\right){\hat{\beta}}_{gwas2}\\[8pt] {}\mathit{\operatorname{var}}\left({\hat{\beta}}_{resid1}\right)=\frac{\left(y\hbox{'}y-\hat{\beta}{\hbox{'}}_{jo\mathit{\operatorname{int}}2} daig\left(X{\hbox{'}}_2X\right){\hat{\beta}}_{gwas2}-{\hat{\beta}}_{resid1}\mathit{\operatorname{diag}}\left(X{\hbox{'}}_1X\right){\hat{\beta}}_{resid1}\right)}{n-N-1}{\left(X{\hbox{'}}_1X\right)}^{-1}\end{array}} \end{equation}$$

(5)

where X₁ is the genotype matrix of SNPs to be regressed on by the residuals, X₂ is the genotype matrix of the joint model and N is equal to the number of SNPs in the residual model. In practice, the data matrix X is unavailable as summary statistics were used, but it is possible to approximate this matrix using the LD structure from a reference panel [Eq. (6)].

$$\begin{eqnarray} \mathit{\operatorname{diag}}{\left(R\hbox{'}R\right)}_{11} &=& \left(n-1\right)\mathit{\operatorname{var}}\left({R}_1\right)\nonumber\\[6pt] X\hbox{'}X &\approx & \mathit{\operatorname{diag}}\left(R\hbox{'}R\right) cor\left(R\hbox{'}R\right) daig\left(R\hbox{'}R\right) \end{eqnarray}$$

(6)

where R is the reference genotype matrix, sigma is the LD matrix for the locus and the diagonal of R’R is modified to be equal to the sample size of the SNP minus one in the GWAS multiplied by the genotypic variance of the SNP observed in the reference panel. Because the data were generated from dosages and not hard-called genotypes, using the observed genotypic variance in the reference panel would have accounted for some of the uncertainty introduced by imputation.

To calculate the effective number of hypothesis tests, Eigen value decomposition was performed on the SNP correlation matrix for each region, using data from the European individuals from the 1000 Genomes Project. The number of hypotheses tested per region was calculated as the number of Eigen values that were required to explain 0.995 of the total sum of the Eigen values [Eq. (7)]. The total number of hypotheses tested in the conditional analysis was taken as the sum of the per region hypothesis counts (65). This revealed that in the focused conditional analysis, we were performing approximately 5443 hypothesis tests. The multiple-testing threshold for our conditional analysis was therefore determined to be 9.2 × 10⁻⁰⁶ (0.05/5443).

$$\begin{equation} \frac{\sum_{m=1}^{M_{gao}}{\lambda}_m}{\sum_{m=1}^M{\lambda}_m}\ge c \end{equation}$$

(7)

The following equation was used to calculate variance explained by each SNP in the meta-analysis and joint analysis [Eq. (8)].

$$\begin{equation} {q}^2=\frac{\mathrm{Var}(X){\beta}^2}{\mathrm{Var}(Y)} \end{equation}$$

(8)

where Var(X) was the variance for each SNP, calculated as 2p(1—P) with P as the allele frequency. β was the effect estimate, and Var(Y) was calculated as the pooled variance estimate provided by the GCTA software when performing the conditional analysis. This pooled variance was calculated using the equation below [Eq. (9)].

$$\begin{equation} {s}_p^2=\frac{\sum_{i=1}^k\left({n}_i-1\right){s}_i^2}{\sum_{i=1}^k\left({n}_i-1\right)} \end{equation}$$

(9)

where s_i² was the phenotypic variance estimated by the GCTA software for each chromosome, and n_i was the number of SNPs on each chromosome. This revealed that the empirical variance of sex adjusted serum urate was 1.624. Unadjusted variance in serum urate was also calculated using the equation above, where |${s}_i^2$| was replaced with the variance in serum urate for each study, and n_i was replaced with the number of participants in each study. This analysis revealed that the empirical variance of unadjusted serum urate was 1.964.

Heritability and functional enrichments

LD Score regression was used to partition SNP heritability of serum urate (63). An estimate was generated using LD Score for the amount of heritability explained by all SNPs additively in the Köttgen et al. (4) meta-analysis. We also performed functional annotation-partitioned LD score regression to determine which cell type groups and cell types contribute significantly to the heritability of serum urate (66). The comprehensive set of functional annotations that were released with partitioned LD Score regression (https://data.broadinstitute.org/alkesgroup/LDSCORE/) were used. This compares the results to a baseline model that contains annotations such as evolutionary conservation and pooled cell type annotations such as DNase1 hypersensitivity. We calculated the enrichment of specifically expressed genes using LD-score regression (67). To assess the enrichment of eQTLs, we performed LD score regression using annotations derived from the GTEx eQTL data with the AllCisQTL annotation used to calculate the overall enrichment of eQTLs across all tissues (68). We separately generated this annotation for kidney cortex data from GTEx v8. We calculated P-values and Bonferroni-corrected thresholds by dividing by the total number of tests within each of the cell type group and cell-type-specific analyses, noting that this is a conservative adjustment because the annotations are correlated. Benjamini–Hochberg false discovery rate (FDR) adjusted P-values were also calculated (69). All results were visualized using ggplot2 (70).

Functional trans-ancestral fine mapping with PAINTOR

PAINTOR (v3.0) (71) was initially used to fine map the 38 loci associated at a genome-wide level of significance in this study with serum urate within the separate European and East Asian GWAS. This initial analysis revealed that both the RELA and HLA loci were inappropriate loci for trans-ancestral fine mapping. For the RELA locus, the association signal between the European and East Asian GWAS clearly involves different causal variants. For HLA-DQB1, the large number of SNPs in the region resulted in computational errors in the PAINTOR software. Both of the loci were excluded from all additional PAINTOR analyses. Cell-type groups that were significant in the LD score regression analysis were used with PAINTOR. To assess how much the East Asian GWAS and these functional annotations improved serum urate fine mapping, the average size of the 90% causal credible sets in three analyses was calculated.

1. European GWAS

2. European and East Asian GWAS

3. European and East Asian GWAS and functional annotations

Given the difficulty in comparing results from loci such as SLC2A9, ABCG2, SLC17A1 and SLC22A12 that have strong signal over a relatively large region and extensive linkage disequilibrium (resulting in a larger total posterior probability) to loci with weaker signal over a smaller region and less linkage disequilibrium (a lower total posterior probability), we also calculated the proportion of total posterior probability for each SNP at their respective loci.

Cis-eQTL analysis

We used COLOC (72) to colocalize the urate-associated loci with publicly available eQTL data from the Genotype Tissue Expression Project (GTEx v6p) (73) and Dorsal Lateral Prefrontal Cortex (DLFPC) eQTL summary statistics from the CommonMind Consortium (27). We also performed colocalization using eQTL summary statistics for the kidney cortex from GTEx v8 (74) and micro-dissected tubulointerstitial and glomerular transcriptomes from patients with nephrotic syndrome (75). COLOC is a Bayesian method that compares four different statistical models at a locus. These models are: no causal variant in the GWAS or the eQTL region; a causal variant in either the GWAS or the eQTL region, but not both; different causal variants in the GWAS and the eQTL region; or a shared causal variant in the GWAS and the eQTL region. All the cis-eQTL regions from a GTEx tissue were merged with the genome-wide European serum urate GWAS data. Genes that were annotated as novel transcripts were removed. For learning the priors, each cis-eQTL region was treated as independent and the likelihood was maximized using the Nelder–Mead algorithm. Genes that had a PPC greater than 0.8 were considered to have a shared causal variant with serum urate. We did not restrict our analysis only to the genome-wide significant loci, which made it possible to identify novel serum urate loci. If multiple tissues supported colocalization at probability > 0.8 the posterior probability was averaged. Colocalization analysis was also performed for all conditional associations that we identified. In this case, we generated association traces for each independent association by first regressing out any other independent associations. To generate these traces using only summary statistics, we used COCO with the serum urate GWAS and LD information from the UK Biobank.

Trans-eQTL analysis

The Contextualize Developmental SNPs using 3D Information (CoDeS3D) algorithm (GitHub, https://github.com/alcamerone/codes3d) (76) was used to identify long-distance regulatory relationships for serum urate-associated SNPs. This analysis leverages known spatial associations from Hi-C databases (77) and gene expression associations (eQTL data from the GTEx catalogue (72)) to assess regulatory connections. Although it would be ideal to use to use chromatin interaction and eQTL data from the same source, we used heterogeneous sources of data (Supplementary Material, Table S6) because a single data source that comprises both types is currently not available. However, the presence of eQTL between SNP and gene pairs from non-paired sources demonstrates the robustness of the process. Briefly, SNPs were mapped onto Hi-C restriction fragments, the genes that physically interact with these restriction fragments identified and collated (SNP-gene spatial pairs). SNP-gene pairs were screened through GTEx to identify eQTL. The FDR was calculated using a stepwise Benjamini–Hochberg correction procedure and incorporated the number of tests and eQTL value list. An FDR value of < 0.05 was accepted as statistically significant (76). COLOC was then used to co-localize trans-eQTL with serum urate GWAS signals. A limitation of trans-eQTL approaches is that there is a high false-positive rate (78)—to mitigate this, replication would be ideal. However, this is not possible owing to the lack of availability of a clinical resource of sufficient size with Hi-C, RNA-Seq and genotype data.

Gout case-control sample sets for replicating serum urate associations

The Japanese gout data set, generated as previously described (30), consisted of 945 male gout patients and 1213 male controls, where gout was clinically ascertained. The Chinese data set, generated as previously described (32), consisted of 1255 male gout cases and 1848 male control where gout was clinically ascertained according to the American College of Rheumatology diagnostic criteria. The European gout dataset was generated from 7342 gout patients and 352 534 controls of European ancestry from the UK Biobank (31), where gout was ascertained by self-report of physician-diagnosed gout or use of urate-lowering therapy (79). Gout association in UK Biobank was tested using logistic regression, adjusted by age, sex and the first 10 principal components (out of 40).

SLC22A9—cell lines, RNA extraction and RT-PCR

Human kidney proximal tubule epithelial cell line (PTC-05) was obtained from Ulrich Hopfer (Case Western Reserve University, Cleveland, Ohio) and grown (37°C in a 5% CO₂) on type IV collagen-coated Petri dish in a 1:1 mixture of DMEM and HAM’S F12 media containing 5 mm glucose, 10% fetal bovine serum (FBS), 2 mm glutamine, 1 mm pyruvate, 5 μg/ml transferrin, 5 μg/ml insulin, 10 ng/ml human epidermal growth factor, 4 μg/ml dexamethasone, 15 mm HEPES (pH 7. 4), 0.06% NaHCO₃, 10 ng/ml interferon-gamma, 50 μM ascorbic acid, 20 nM sodium selenite (Na₂SeO₃), 1 nM triiodothyronine (T3) and penicillin (50 units/ml)/streptomycin (50 μg/ml). Human embryonic kidney HEK293 cells (ATCC) and human hepatocellular carcinoma HEPG2 cells (ATCC) were grown (37°C in a 5% CO₂) and maintained in Dulbecco’s Modified Eagle’s Medium (DMEM) and Eagle’s Minimum Essential Media (EMEM), respectively, supplemented with 4.5 g/l glucose, 2 mm glutamine, 1 mm sodium pyruvate, 10% FBS and penicillin (50 units/ml)/streptomycin (50 μg/ml).

Total RNA from human cell lines (PTC-05, HepG2 and HEK-293 T) was extracted using spin columns with the RNeasy Mini Kit (QIAGEN, GmbH, Germany) following the manufacturer’s instructions. Approximately, 2 μg of DNase-treated total RNA, isolated from cells, was primed with poly-dT and random hexamers and then reverse-transcribed using AMV reverse transcriptase (New England Biolabs, Ipswich, MA). An equal amount of cDNA was used for PCR amplification of OAT7 and GAPDH cDNAs using the following primers, followed by electrophoresis.

hOAT7-2S [sense] 5′-CAACCTCAATGGCCTTTCAGGACCTCCTGG-3′.

hOAT7-3A [antisense] 5′-GCCTGGAATCTGTGTGTTGCCCACTCGG-3′.

hGAPDH-1S [sense] 5′-CGGAGTCAACGGATTTGGTCGTATTG-3′.

hGAPDH-1A [antisense] 5′-GACTGTGGTCATGAGTCCTTCCACGA-3′.

SLC22A9—urate transport analysis of OAT7

Studies using Xenopus laevis oocytes were performed in accordance with the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the U.S. National Institutes of Health and were approved by the Institution’s Animal Care and Use Committee. Mature female Xenopus laevis frogs (NASCO, Fort Atkinson, MI) were subjected to partial ovariectomy under tricaine (SIGMA, St Louis, MO) anesthesia (0.17% for 15–20 min) as described previously (52). A small incision was made in the abdomen and a lobe of ovary was removed. Subsequently, the oocytes were pre-washed for 20 min in Ca²⁺-free ND96 medium (96 mm NaCl, 2 mm KCl, 1 mm MgCl₂, 5 mm HEPES, pH 7.4) to remove blood and damaged tissue. Oocytes were then defolliculated by treatment with 3.5 mg/ml of collagenase enzyme (Roche, Indianapolis, IN) in Ca²⁺-free ND96 medium for about 120 min with gentle agitation at room temperature (25°C). Subsequent to this treatment, oocytes were washed three times with ND96 medium and incubated (16–18°C) in isotonic Ca²⁺-containing ND96 medium (96 mm NaCl, 2.0 mm KCl, 1.8 mm CaCl₂, 1.0 mm MgCl₂ and 5 mm Hepes, pH 7.4) supplemented with 2.5 mm pyruvate and gentamycin (10 μg/ml).

For expression of OAT7 and OAT1 in Xenopus laevis oocytes, their respective full-length cDNAs were cloned into the pGEMHE vector, wherein the cDNA insert is flanked by the Xenopus laevis β-globin 5′-UTR and 3′-UTR (80). These constructs were linearized and cRNAs were synthesized in vitro using T7 RNA polymerase (mMESSAGE mMACHINE; Ambion, Austin, TX) following the supplier’s protocol. Isopropanol-precipitated, in vitro transcribed capped cRNAs were washed twice with 70% ethanol, the cRNA pellet was dried and then dissolved in sterile nuclease-free water. The yield and integrity of the capped cRNA samples was assessed by spectroscopy (at 260 nm) and 1% agarose–formaldehyde gel electrophoresis, respectively. All cRNA samples were stored frozen in aliquots at −80°C until used.

About 18 hours after isolation, oocytes were microinjected with 50 nl of sterile water, 50 mm Tris pH 7.4, or 50 nl of a cRNA solution in 50 mm Tris buffer (pH 7.4) containing 25 ng of the indicated cRNA using fine-tipped micropipettes by a microinjector (World Precision Instrument Inc. Sarasota, FL). The microinjected oocytes were then incubated in isotonic ND96 medium (pH 7.4) containing 1.8 mm CaCl₂, 2.5 mm pyruvate, gentamycin (10 μg/ml) at 16–18°C for approximately 48 h to allow expression of protein from microinjected cRNA.

For [¹⁴C]-urate (specific activity: 50 mCi/mmol) uptake experiments in Xenopus laevis oocytes, oocytes expressing proteins as indicated (OAT7 and OAT1) were washed four times with ND96 medium (96 mm NaCl, 2.0 mm KCl, 1.8 mm CaCl₂, 1.0 mm MgCl₂ and 5 mm Hepes, pH 7.4) without pyruvate and gentamycin. OAT7 functions as a butyrate exchanger (36), therefore OAT7-expressing oocytes were microinjected with 50 nl of 100 mm butyrate to optimize urate transport by ‘trans-activation’ (52). After approximately 60 min of starvation, oocytes were preincubated in the ND96 uptake medium for 30 min before incubation (25°C, in a horizontal shaker-incubator) in the uptake medium containing [¹⁴C]-urate (40 μM). After 60 min of incubation in the uptake medium, oocytes (20 per group) were washed three times with ice-cold uptake medium to remove external adhering radioisotope. OAT7-expressing oocytes were then exposed to DMSO (diluent for uricosurics) or the uricosuric drugs tranilast and benzbromarone, as indicated. The radioisotope content of each individual oocyte was measured by scintillation counter following solubilization in 0.3 ml of 10% (v/v) SDS and addition of 2.5 ml of scintillation fluid (Ecoscint). All uptake experiments included at least 20 oocytes in each experimental group; statistical significance was defined as two-tailed P < 0.05, and results were reported as means ± S.E. Statistical analyses including linear regressions and significance were determined by Student’s t-test using SigmaPlot software.

Data Availability

The summary statistics for the trans-ancestral meta-analysis are available for download (DOI: 10.5281/zenodo.3366490) as are the summary statistics for the genome-wide colocalization analysis (https://zenodo.org/record/3590410#.Xf_nKNbYryw). All data used in this research to generate or further analyze these summary statistics were sourced from third parties. Publicly available data include the European serum urate GWAS data (http://metabolomics.helmholtzmuenchen.de/gugc/), the GTEx eQTL data (https://gtexportal.org), the 1000 Genomes data (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/) and the functional annotations for partitioned LD Score regression (https://data.broadinstitute.org/alkesgroup/LDSCORE/). The East Asian serum urate GWAS data are available on request for researchers who meet the criteria from Yuki Okada (Osaka University, [email protected]u.ac.jp). The UK Biobank data are available upon request for researchers who meet the criteria from the UK Biobank Access Management System (https://bbams.ndph.ox.ac.uk/ams/). All other relevant data are within the manuscript and its supporting information files.

Abbreviations

GWASgenome-wide association studies

eQTL expression quantitative trait loci

OAT7 organic anion transporter 7

GTEx project Genotype-Tissue Expression

CoDeS3D Contextualize Developmental SNPs using 3D Information’

PAINTOR Probabilistic Annotation INtegraTOR

SNP single nucleotide polymorphism

Acknowledgements

This research has been conducted using the UK Biobank Resource under Application Number 12611. Akiyoshi Nakayama is thanked for statistical analytic support. The Health Research Council of New Zealand is acknowledged for funding support (grant 14/527).

Conflict of Interest statement. The authors have declared that no competing interests exist.

Author notes

James Boocock, Megan Leask and Tony R. Merriman contributed equally to the work.