Risk of generalized vitiligo is associated with the common 55R-94A-247H variant haplotype of GZMB (encoding granzyme B).

In generalized vitiligo (GV), patches of depigmented skin and overlying hair (Picardo and Taïeb, 2010) result from autoimmune destruction of melanocytes by skin-homing auto-reactive cytotoxic T-lymphocytes (CTLs) (Ogg ). In two genome-wide association studies (GWAS) to identify GV susceptibility loci in European-derived white (EUR) patients (GWAS1, GWAS2; Jin , 2012) we detected association with GZMB, encoding granzyme B, a serine protease marker of activated CTLs (Sattar et al., 2003). In concert with perforin, GZMB mediates direct and caspase-mediated apoptosis of target cells (Granville, 2010; Ewen ), as well as proteolytic cleavage of autoantigens, creating or exposing autoimmune epitopes that may initiate or propagate the autoimmune process (Darrah and Rosen, 2010).

Meta-analysis of GWAS1, GWAS2, and replication study data for SNPs genotyped in the GZMB region (Jin ) showed greatest association with rs8192917-C (P = 5.60 × 10−9, OR 1.22), a common non-synonymous SNP (R55Q) that is in very strong linkage disequilibrium with two other common non-synonymous SNPs, rs11539752 (P94A; r2 = 0.99) and rs2236338 (Y247H; r2 = 0.93). Together, these three variants define two predominant haplotypic multi-variant GZMB polypeptide isoforms 55Q-94P-247Y (termed “QPY”) versus 55R-94A-247H (“RAH”) (McIlroy ; Zaitsu ). To carry out higher-resolution analysis of association of GV with SNPs in the GZMB region, we first imputed genotypes for all SNPs in the 1000 Genomes Project phase 1 integrated variant call set (June 21, 2012 release) across the 160 kb region (chr14:25045130-25204958) showing even nominal association with GV in the combined GWAS1 and GWAS2 dataset (Jin et al., 2010, 2012), from 101.5 kb upstream of GZMB through 55 kb downstream. After quality control procedures, there were 464 genotyped or imputed SNPs in 7202 total subjects (Supplementary Table S1 online). In this analysis, the most significant SNP was rs10909625 (P = 1.03 × 10−8), a synonymous variant at codon K80 (Table 1). To determine which SNPs in the GZMB region might be causal for GV susceptibility, we then carried out logistic regression, individually testing each of 463 SNPs conditional on rs10909625. This analysis indicated there is only one primary association signal in the GZMB region, represented by nine SNPs (Table 1) whose effects cannot be distinguished (Supplementary Table 2) due to very strong linkage disequilibrium (Supplementary Figure 1). These include rs8192917 (R55Q) and rs11539752 (P94Q); none of the other seven primary associated SNPs (rs6573910, rs6573911, rs45628336, rs113822535, rs45442494, rs1126639, rs10909625) are predicted to have functional consequences. In contrast, the associations of rs2236338 (Y247H), as well as that of another potentially functional SNP rs2273844 (stop codon within the 5′ untranslated region), were found to be secondary.

Table 1

Association of GV with SNPs in the GZMB region of chromosome 14q12

SNP	nt	A1	A2	A1freq	GWAS1 P	GWAS1OR	GWAS2 P	GWAS2OR	CMH P	CMHOR	CMHL95	CMHU95	Call type	Function
rs2236338	25100282	G	A	0.24	4.18 × 10−5	1.25	2.96 × 10−2	1.20	3.32 × 10−6	1.23	1.13	1.35	genotyped	Y247H
rs6573910	25100882	T	C	0.25	3.15 × 10−6	1.28	1.82 × 10−2	1.22	1.84 × 10−7	1.27	1.16	1.38	imputed
rs6573911	25100933	T	C	0.25	3.15 × 10−6	1.28	1.82 × 10−2	1.22	1.84 × 10−7	1.27	1.16	1.38	imputed
rs45628336	25101440	T	C	0.24	1.48 × 10−6	1.30	1.88 × 10−2	1.22	9.69 × 10−8	1.27	1.17	1.39	imputed
rs113822535	25101465	T	G	0.24	9.11 × 10−7	1.30	2.08 × 10−2	1.22	6.60 × 10−8	1.28	1.17	1.40	imputed
rs45442494	25101475	T	G	0.24	6.20 × 10−7	1.31	1.84 × 10−2	1.22	4.01 × 10−8	1.28	1.17	1.40	imputed
rs1126639	25101548	A	G	0.24	2.38 × 10−6	1.29	2.19 × 10−2	1.21	1.74 × 10−7	1.27	1.16	1.38	imputed
rs11539752	25101589	C	G	0.24	2.98 × 10−6	1.29	2.32 × 10−2	1.21	2.27 × 10−7	1.26	1.16	1.38	imputed	P94A
rs10909625	25101629	C	T	0.24	9.15 × 10−7	1.31	1.22 × 10−2	1.24	3.92 × 10−8	1.28	1.17	1.40	imputed
rs8192917	25102160	C	T	0.25	2.83 × 10−6	1.29	2.39 × 10−2	1.21	2.25 × 10−7	1.26	1.16	1.38	genotyped	R55Q
rs2273844	25103414	A	G	0.24	4.81 × 10−6	1.28	3.00 × 10−2	1.20	4.78 × 10−7	1.26	1.15	1.37	genotyped	Stop gained

Data are shown for nine genotyped or imputed SNPs (Supplementary Table S1) whose effects could not be distinguished by logistic regression analysis, due to linkage disequilibrium (Supplementary Table S2), as well as for two additional potential functional SNPs, rs2236338 (Y247H) and rs2273844 (stop codon within 5′ untranslated region). nt, nucleotide position on chromosome 14 (alleles denoted on forward strand)

A1, effect allele (here, high-risk)

A2, reference allele

A1 freq, A1 allele frequency in summary 7202 GV

cases and controls

GWAS 1, data from 1388 GV cases and 2586 controls ()

GWAS2, data from 418 GV cases and 2810 controls ()

OR, odds ratio

CMH, Cochran-Mantel-Haenszel meta-analysis

L95, lower limit of 95% confidence interval; U95, upper limit of 95% confidence interval. Note that subject quality control was carried out using the combined GWAS1+GWAS2 dataset, and subjects in GWAS1 who were related to subjects in GWAS2 (pi-hat > 0.05

4 GV cases, 43 controls) were removed. Cochran-Mantel-Haenszel analysis here does not include replication studies of GWAS1 or GWAS2, which could not be imputed

therefore, P-values shown here for some SNPs are less significant here than reported in Jin et al., 2012 because of lower sample size.

Analysis of the haplotypes defined by the three non-synonymous SNPs rs8192917-rs11539752-rs2236338 similarly indicated that principal association with GV is with the haplotype comprised of rs8192917-C—rs11539752-C. As shown in Table 2, this constitutes the higher-risk haplotype (OR 1.27; P = 4.42 × 10−8), compared to the haplotype that also includes the high-risk G allele of rs2236338 (OR 1.26; P = 3.26 × 10−7). Nevertheless, it is difficult to completely exclude any effect of rs2236338, as the haplotypes containing this SNP (C-C-G and C-C-A) are also associated with elevated risk for GV.

Table 2

Analysis of haplotypes of GZMB non-synonymous SNPs rs8192917 (R55Q), rs11539752 (P94A), and rs2236338 (Y247H)

SNPHaplotype	Amino AcidHaplotype	F_A	F_U	P	OR
3-SNP Haplotypes (rs8192917-rs11539752-rs2236338)
TGA	55Q-94P-227Y	0.716	0.760	1.08E-07	ref
CCG	55R-94A-247H	0.270	0.228	3.26E-07	1.26
CCA	55-R-94A-247Y	0.009	0.005	0.02587	1.69
TGG	55Q-94P-247H	0.006	0.006	0.5314	0.9
2-SNP Haplotypes (rs8192917-rs11539752)
TG	55Q-94P	0.7211	0.7663	4.42E-08	ref
CC	55R-94A	0.2789	0.2337	4.42E-08	1.27

P values compare the frequency of the designated haplotype with all other haplotypes in cases versus controls. F_A, frequency in cases

F_U, frequency in controls

OR, odds ratio compared to the reference (ref) haplotype.

Genotypes of a number of uncommon (MAF < 0.01) potentially functional SNPs in the GZMB region could not be imputed with high confidence (MaCH r2 < 0.3). To assess whether association of GV with GZMB might result from rare deleterious variants “hitchhiking” on the most associated haplotype of common SNPs, we therefore carried out next-generation DNA re-sequencing of a total 8.3 kb across the GZMB locus, spanning 3.2 kb upstream of the gene, the five exons and intervening sequences (with a single 783 nt gap within intervening sequence 1), to 3.1 kb downstream, in 114 unrelated EUR GV patients from GWAS1 (Jin ). As shown in Supplementary Table S3 online, in these 228 alleles we observed a total of 41 variants. In addition to the three aforementioned common non-synonymous SNPs, we observed one additional rare non-synonymous SNP, rs185979723 (M225L) in a single heterozygote, predicted to be functionally benign. There were no nonsense, frameshift, or splice junction mutations, and no indels or other apparent rearrangements. We also observed five previously unreported variants intergenic and intronic variants, none of which were of predicted functional significance. Our findings thus suggest a major role for the common rs8192917-C—rs11539752-C (rs2236338-G) haplotype in genetic association of GZMB with GV, with no evidence for a major contributation attributable to rare causal variants.

The mechanism by which the GZMB rs8192917-C—rs11539752-C—(rs2236338-G) haplotype increases GV risk is not yet known, as functional correlates of the corresponding 55R-94A-(247H) and 55Q-94P-(247Y) GZMB protein variant isoforms remain uncertain. Individually, the R55Q, A94P, and Y247H substitutions are each predicted to be benign (not shown); however, in silico functional predictions of individual amino acid variants are of unknown validity in the context of such multi-variant haplotypes. We speculate that the 55R-94A-247H and 55R-94A-247Y GZMB polypeptides may both increase risk of GV, the former perhaps to a greater degree than the latter. The GZMB QPY and RAH alleles exhibit similar expression, stability, and proteolytic activity of the corresponding alternative GZMB polypeptide isoforms (McIlroy ). However, there have been conflicting claims regarding their possible differential effects on CTL effector functions. Homozygotes for the rs8192917-T (55Q) allele have been reported to express a greater proportion of GZMB-positive cells after mitogen stimulation (Girnita et al., 2009), while two studies reported reduced apoptotic induction by RAH versus QPY GZMB (McIlroy ; Gaafar et al., 2009), and another found no difference (Sun ). Regardless, we cannot exclude the alternative possibility that association of all GZMB coding variants with GV is secondary to LD with an as-yet unidentified causal non-coding SNP.

GZMB has not been genetically associated with other autoimmune diseases, specifically including juvenile idiopathic arthritis (Donn ) and Behçet’s disease (Kucuksezer ). This suggests the possibility that GZMB may be relatively specific for melanocyte-directed autoimmune susceptibility. In addition to the effector function of GZMB in target cell killing by CTLs, perhaps particularly relevant to GV may be the role played by GZMB cleavage in activating autoantigens in autoimmune disease (Darrah and Rosen, 2010). GZMB cleaves target proteins at aspartic acid residues (Odake ), including both caspases associated with induction of apoptosis and most target cell autoantigens, potentially enhancing presentation of cryptic epitopes and leading to activation of self reactive T-cells (Sercarz ). GZMB cleavage sites of GV melanocyte autoantigens have not yet been defined experimentally. However, in silico analysis of the principal human GV autoantigens using GRASVM 1.0 (Wee et al., 2011) predicts one high-probability GZMB cleavage site in each of TYR and TYRP1, three in TYRP2, two in SOX10, one in PMEL, and none in MC1R. These findings thus suggest that in vivo GZMB cleavage of melanocyte proteins that constitute GV autoantigens may contribute to initiation and propagation of autoimmunity directed against melanocytes.

Supplementary Figure 1. Linkage disequilibrium (r2) of 11 GZMB SNPs from Table 1

Supplementary Table S1. Association of GZMB region SNPs with GV

Supplementary Table S2. Association of GV with most significant and potential functional SNPs conditional on rs10909625

Supplementary Table S3. GZMB Sequence Variants Observed in 114 unrelated GV patient

Supplementary Table S4. PCR Primers for DNA Sequencing of the GZMB Region