Genome-wide association study identifies 12 new susceptibility loci for primary biliary cirrhosis.

In addition to the HLA locus, six genetic risk factors for primary biliary cirrhosis (PBC) have been identified in recent genome-wide association studies (GWAS). To identify additional loci, we carried out a GWAS using 1,840 cases from the UK PBC Consortium and 5,163 UK population controls as part of the Wellcome Trust Case Control Consortium 3 (WTCCC3). We followed up 28 loci in an additional UK cohort of 620 PBC cases and 2,514 population controls. We identified 12 new susceptibility loci (at a genome-wide significance level of P < 5 × 10⁻⁸) and replicated all previously associated loci. We identified three further new loci in a meta-analysis of data from our study and previously published GWAS results. New candidate genes include STAT4, DENND1B, CD80, IL7R, CXCR5, TNFRSF1A, CLEC16A and NFKB1. This study has considerably expanded our knowledge of the genetic architecture of PBC.

Primary biliary cirrhosis is a chronic, autoimmune liver disease characterized by non-suppurative destructive cholangitis and highly specific auto-antibodies to pyruvate dehydrogenase complex (PDC)1. It is an important cause of chronic liver disease and a well-established indication for liver transplantation. In the UK, the prevalence of PBC is approximately 35 per 100,000 adults and 94 per 100,000 women aged 40 years or older2. PBC has a sibling relative risk of ~10, suggesting a substantial genetic contribution to PBC3.

Previous studies have established that PBC is associated with HLA-DR8, with odds ratio ranging from 2.4 to 3.3 depending on the population examined4. To date, six non-HLA loci have been associated with PBC at a genome-wide level of significance (P<5×10−8): IL12A (3q25)5, IL12RB2 (1p31)5, IRF5/TNPO3 (7q32)6,7, ORMDL3/IKZF3 (17q12)6,7, MMEL1 (1p36)6, and SPIB (19q13)7. To identify additional PBC risk loci we conducted a GWAS in a large cohort of UK PBC cases and population controls as part of the Wellcome Trust Case Control Consortium 3.

PBC cases were drawn from the UK PBC Consortium, which consists of 142 NHS Trusts in the UK, including all liver transplant centers. All cases were of self-declared British or Irish ancestry. PBC cases were genotyped using the Illumina 660W-Quad array. UK population controls were genotyped using the Illumina Human1M-Duo by the Wellcome Trust Case Control Consortium 2 and consisted of individuals from the 1958 British Birth Cohort and National Blood Service8. Following stringent quality control (see Online Methods), 507,467 SNPs were available across 1,840 cases and 5,163 ‘historical’ population controls (see Supplementary Tables 1 and 2). The quantile-quantile plot of the case-control chi-square test statistics demonstrates a substantial excess of significant associations in the tail of the distribution, even after removal of known loci, which cannot be ascribed to overall inflation of the distribution (genomic control λ=1.09; see Supplementary Figure 1).

We identified 34 loci where one or more SNPs showed at least suggestive evidence for association (P<1×10−5), including six of the seven previously associated loci (Table 1). We found weaker evidence for replication at the seventh previously associated locus 1p36 (containing MMEL1; P=4×10−3). 28 loci achieving at least suggestive significance and not previously associated with PBC at genome-wide significance were followed up by genotyping 46 SNPs in an independent panel of 620 cases from the UK PBC Consortium (Supplementary Table 3). A comparison set of 2,514 ‘historical’ UK population controls previously genotyped by TwinsUK using the Illumina HumanHap610 array was obtained (see Online Methods). Twelve of these 28 loci were significant at P<0.05 in the replication study, and P<5×10−8 in the combined analysis of the discovery and replication cohorts, and thus represent novel PBC associations (Table 2). Two of these loci (2q32 and 1q31) showed suggestive evidence of association in a previous PBC GWAS study5. Follow-up genotyping in large independent panels of cases and controls from a range of ethnicities is needed to further characterize these loci.

Table 1

Association results and in silico analyses for the 7 previously confirmed primary biliary cirrhosis risk loci.

CHR	SNP	RISKALLELE	LEFT/RIGHT REGION(MB)	CANDIDATE GENEa(Number of genes in region)	AIb	nsSNPc	GWAS Cohort
							Control RAFd	Case RAFd	P-value	OR(95% CI)
1p36	rs10752747	T	2.39 - 2.78	MMEL1 (6)	○		0.339	0.367	2.65×10−3	1.13(1.04 - 1.22)
1p31	rs17129789	C	67.53 – 67.71	IL12RB2 (2)	○		0.177	0.247	9.48×10−20	1.52(1.39 - 1.67)
3q25	rs485499	T	160.96 - 161.3	IL12A (2)	●		0.574	0.651	2.29×10−16	1.38(1.28 - 1.50)
6p21	rs7774434	C	26.21 - 33.74	Many (MHC)	○	●	0.379	0.494	3.86×10−34	1.60(1.48 - 1.73)
7q32	rs12531711	G	128.33 - 128.57	IRF5 (2)	●		0.107	0.159	8.90×10−17	1.58(1.41 - 1.76)
17q12	rs7208487	T	34.61 - 35.49	ORMDL3 (23)	●		0.840	0.874	7.89×10−7	1.32(1.18 - 1.48)
19q13	rs3745516	A	55.52 - 55.73	SPIB (11)		●	0.226	0.287	1.63×10−13	1.38(1.32 – 1.44)

The putative candidate gene represents the strongest candidate within the region based on available evidence, but does not preclude the existence of other plausible candidate genes within the region. The number of genes is based upon the RefSeq gene track.

Open circles indicate the locus has been previously associated with another autoimmune disease; a filled circle indicates the same candidate gene has also been suggested.

Filled circles indicate a non-synonymous SNP in LD with our top SNP was identified in the candidate gene.

RAF indicates the risk allele frequency.

Table 2

Association results and in silico analyses for 12 newly confirmed primary biliary cirrhosis risk loci.

CHR	SNP	RISKALLELE	LEFT/RIGHTREGION (MB)	CANDIDATEGENEa(Number ofgenes in region)	AIb	GRAILc	nsSNPd	GWAS cohort				Replication cohort				Combined samples
								ControlRAFe	CaseRAFe	P-value	OR(95% CI)	ControlRAFe	CaseRAFe	P-value	OR(95% CI)	P-value	OR(95% CI)
1q31	rs12134279	T	195.58 - 196.21	DENND1B (4)	○	●		0.202	0.250	1.07×10−9	1.32(1.21 - 1.44)	0.201	0.263	2.57×10−6	1.42(1.23 - 1.65)	2.06×10−14	1.34(1.25 - 1.45)
2q32	rs10931468	A	190.77 - 191.61	STAT4 (7)	●	●		0.119	0.164	2.55×10−12	1.46(1.31 - 1.62)	0.120	0.183	2.64×10−9	1.64(1.39 - 1.94)	2.35×10−19	1.50(1.37 - 1.64)
3q13	rs2293370	G	120.58 - 120.79	CD80 (6)	●	●		0.804	0.853	7.70×10−11	1.41(1.27 - 1.56)	0.809	0.835	0.036	1.19(1.01 - 1.41)	2.53×10−11	1.35(1.23 - 1.47)
4q24	rs7665090	C	103.61 - 104.24	NFKB1 (7)				0.524	0.572	5.33×10−7	1.21(1.13 - 1.31)	0.513	0.593	5.50×10−7	1.38(1.22 - 1.57)	4.06×10−12	1.26(1.18 - 1.34)
5p13	rs860413	A	35.74 - 36.08	IL7R (5)	●	●	●	0.719	0.773	3.09×10−10	1.33(1.22 - 1.45)	0.729	0.769	4.50×10−3	1.24(1.07 - 1.43)	1.02×10−11	1.30(1.21 - 1.40)
7p14	rs6974491	A	37.32 - 37.41	(0)				0.170	0.205	3.39×10−6	1.25(1.14 - 1.38)	0.177	0.215	2.40×10−3	1.27(1.09 - 1.49)	4.44×10−8	1.25(1.16 - 1.36)
11q23	rs6421571	C	117.82 - 118.30	CXCR5 (10)	○			0.809	0.855	3.53×10−10	1.40(1.26 - 1.55)	0.810	0.847	2.10×10−3	1.30(1.10 - 1.55)	2.69×10−12	1.37(1.25 - 1.50)
12p13	rs1800693	C	6.29 - 6.33	TNFRSF1A (3)	●	●		0.401	0.452	5.51×10−8	1.23(1.14 - 1.33)	0.403	0.445	8.70×10−3	1.18(1.04 - 1.34)	1.80×10−9	1.22(1.14 - 1.30)
14q24	rs911263	T	67.34 - 67.98	RAD51L1 (2)				0.712	0.764	1.68×10−9	1.31(1.20 - 1.43)	0.717	0.760	2.30×10−3	1.25(1.08 - 1.45)	1.76×10−11	1.29(1.20 - 1.39)
16p13	rs12924729	G	10.92 - 11.22	CLEC16A (3)	●			0.679	0.737	7.68×10−11	1.32(1.21 - 1.44)	0.680	0.718	8.80×10−3	1.20(1.05 - 1.38)	2.95×10−12	1.29(1.20 - 1.38)
16q24	rs11117432	G	84.55 - 84.58	(0)	○			0.760	0.808	1.20×10−6	1.26(1.15 - 1.39)	0.774	0.838	9.52×10−7	1.52(1.28 - 1.79)	4.66×10−11	1.31(1.21 - 1.43)
22q13	rs968451	T	37.87 - 38.19	MAP3K7IP1(3)				0.194	0.233	4.31×10−7	1.27(1.16 - 1.39)	0.193	0.237	6.45×10−4	1.30(1.12 - 1.51)	1.08×10−9	1.27(1.18 - 1.38)

PBC loci that meet genome-wide significance P<5×10−8 in the combined analysis and P<0.05 in the replication cohort. GWAS and replication cohort data for the replicated SNPs were merged using PLINK.

The putative candidate gene represents the strongest candidate within the region based on available evidence, but does not preclude the existence of other plausible candidate genes within the region. The number of genes is based upon the RefSeq gene track.

Open circles indicate the locus as been previously associated with another autoimmune disease; a filled circle indicates the same candidate gene has also been suggested.

Filled circles indicate the gene was identified by GRAIL as the most plausible functional candidate in the region (Ptext < 0.01). GRAIL results are not provided for previously confirmed loci as these were used as seeds in the analysis.

Filled circles indicate a non-synonymous SNP in LD (r2 >0.8) with our top SNP was identified in the candidate gene.

RAF indicates the risk allele frequency.

To identify additional risk loci, we combined summary statistics from our discovery cohort with those from the two datasets included in the previously published meta-analysis of PBC GWAS7. Three further novel loci reached genome-wide significance (see Table 3). This included one locus (14q32) that just failed to achieve genome-wide significance in our combined analysis of discovery and replication cohorts (P = 1.69×10−7), but did so with the addition of data from the study by Liu et al.7 (P = 2.61×10−13). As SNPs at the two other loci (3p24, 11q13) were not genotyped in our replication cohort, and the loci were identified based on summary statistics alone, genotyping using an independent technology in additional cohorts is needed to fully validate these associations. A combined GWAS meta-analysis is still warranted because we were only able to meta-analyze the top 100 SNPs from the Liu et al. study. Genome-wide imputation using HapMap3 reference panels did not identify any further genome-wide significant loci (see Online Methods, Supplementary Figure 2, and Supplementary Table 4), although for some loci imputed SNPs provided stronger evidence of association than the genotyped SNPs. No statistically significant gene-gene interactions were detected between associated loci, or after fitting an HLA-risk model (see Online Methods and Supplementary Table 5).

Table 3

Genomic regions reaching genome-wide significance after meta-analysis with Liu et al. (2010) data.

CHR	SNP	RISK ALLELE	LEFT/RIGHTREGION (MB)	CANDIDATEGENEa(Number of genes inregion)	Discovery sample				Liu et al. 2010		Meta-analysis
					ControlRAFb	CaseRAFb	P-value	OR (95% CI)	P-value	ORc	P-value	OR (95% CI)
3p24	rs1372072	A	16.82 - 17.13	PLCL2 (1)	0.365	0.400	1.38×10−4	1.16(1.08 – 1.25)	1.52×10−5	1.27	2.28×10−8	1.20(1.12 – 1.27)
11q13	rs538147	G	63.60 - 64.04	RPS6KA4 (20)	0.606	0.647	1.01×10−5	1.19(1.10 – 1.29)	7.72×10−6	1.28	2.06×10−10	1.23(1.15 – 1.31)
14q32	rs8017161	A	102.54 - 102.68	TNFAIP2 (3)	0.396	0.439	4.71×10−6	1.20(1.11 – 1.29)	4.86×10−7	1.31	2.61×10−13	1.22(1.16 – 1.27)

The putative candidate gene represents the strongest candidate within the region based on available evidence, but does not preclude the existence of other plausible candidate genes within the region. The number of genes is based upon the RefSeq gene track.

RAF indicates the risk allele frequency.

Liu et al. (2010) do not provide confidence intervals for the odds ratios (ORs) estimated from their meta-analysis for these SNPs. None of these genomic regions have previously been associated at genome-wide significance with another autoimmune disease. GRAIL failed to identify any strong candidate genes within these regions and no nsSNPs were identified in high LD (r2>0.8) with the most associated SNP at each locus.

We found evidence for a second independent association at the 3q25 locus containing IL12A and SCHIP1, as did Liu et al7 (see Supplementary Table 6). Three SNPs, located between IL12A and SCHIP1, remained genome-wide significant following a conditional logistic regression adjusting for the most significant SNP in the region (rs485499). Conducting the same analysis using the imputed data identified a further two SNPs reaching genome-wide significance. These five SNPs are all in linkage disequilibrium (LD; r2>0.2) with each other, but none are in LD with rs485499. They are located downstream of SCHIP1, but upstream of IL12A (see Supplementary Figure 3). Fine mapping of this locus is needed to determine whether these association signals implicate independent variants affecting the same gene, or two different genes.

We identified plausible candidate genes within associated loci via manual curation, supported by evidence from: a) previous GWAS findings for other autoimmune diseases; b) GRAIL9, a literature-mining tool that identifies non-random, evidence-based links between genes; c) identification of non-synonymous SNPs in 1000 genomes data that are in LD (r2>0.8) with the most associated genotyped SNP in each locus; d) identification of eQTL within associated loci that are in LD (r2>0.8) with the most associated SNP at that locus, using data from Dixon et al.10 (see Online Methods and Supplementary Tables 7-10). Even in aggregate these analyses do not identify the gene(s) containing causal variants but they allow us to identify potential candidate genes for future follow-up studies. Supplementary Figure 4 shows all genes within each of the associated loci.

The results from the GWAS of PBC conducted to date provide additional support for the involvement of three pathways previously implicated in the pathogenesis of PBC: NF-κB signaling, T-cell differentiation, and Toll-like receptor (TLR) and Tumor Necrosis Factor (TNF) signaling.

We identified several loci containing genes involved in activation of NF-κB, a transcription factor which regulates expression of many genes involved in the immune response and is highly activated in other autoimmune disorders such as rheumatoid arthritis, multiple sclerosis, and asthma11. Its importance in PBC is suggested by evidence that NF-κB modulates the balance of survival and apoptosis in activated hepatic stellate cells12, and NF-κB p50 −/− mice show aggressive hepatic inflammation and fibrosis13. The locus we identified at 4q24 contains the NFKB1 gene itself, and we identified genes in pathways leading to NF-κB activation at four other loci: 22q13 (TAB1), 12p13 (TNFRSF1A), 3q13 (CD80), and 11q13 (RPS6KA4).

Loci identified to date suggest a role for T-lymphocyte differentiation in the development of PBC. TH1 immune responses have been implicated in many autoimmune diseases14 and may be involved in development of autoreactive T-cells, consistent with the putative role of PDC-specific autoreactive TH1 cells in the pathogenesis of human PBC and animal disease models15. IL-12 signaling promotes TH1-type immune responses by driving differentiation of activated, naïve T-cells to TH1 cells16 and three loci containing genes involved in IL-12 signaling have been identified for PBC: 3q25 (IL12A) and 1p31 (IL12RB2) by Hirschfield et al.5, and 2q32 (STAT4) in this study. These results provide further support for the TH1 hypothesis regarding PBC development.

Activation of TLR signaling, and its downstream effectors such as TNFα, is well described in PBC17. The 7q32 locus, identified by Liu et al.7, contains IRF5 which is activated in response to TLR-signaling and leads to selective expression of TNFα. We identified a locus at 11q13 containing RPS6KA4, which suppresses TLR-dependent cytokine production18. TNFα is an activating factor for a number of intracellular pathways that determine the fate of hepatocytes, and thus plays a key role in liver homeostasis19. We identified three loci containing genes in TNFα signaling pathways: 12p13 (TNFRSF1A), 1q31 (DENND1B), and 14q32 (TNFAIP2). TNFRSF1A is one of two receptors for TNFα, and TNFRSF1A −/− mice show attenuated liver fibrosis when compared to wild-type mice after administration of a potent hepatotoxin20. DENND1B interacts directly with TNFRSF1A21 and has previously been associated with asthma22. TNFα signaling also directly induces TNFAIP2 expression23.

In summary, this is the first report in a new series of GWAS undertaken by the WTCCC3. Twelve novel PBC risk loci have been identified in this study of >7,000 European samples, making this the largest GWAS of PBC to date. In addition, a further three loci achieve genome-wide significance following meta-analysis with published data. For many of the associated loci we have identified plausible candidate genes that support the involvement of the innate and adaptive immune systems in PBC etiology, particularly signaling via the NF-κB, TLR, and TNF pathways, although these findings require confirmation through fine-mapping, gene-expression and functional studies.