Common variants in FOXP1 are associated with generalized vitiligo.

In a recent genome-wide association study of generalized vitiligo, we identified ten confirmed susceptibility loci. By testing additional loci that showed suggestive association in the genome-wide study, using two replication cohorts of European descent, we observed replicated association of generalized vitiligo with variants at 3p13 encompassing FOXP1 (rs17008723, combined P=1.04x10(-8)) and with variants at 6q27 encompassing CCR6 (rs6902119, combined P=3.94x10(-7)).

Generalized vitiligo (GV) is a common, complex autoimmune disease in which patchy depigmentation of skin and hair results from loss of melanocytes from involved areas1, and which is epidemiologically associated with several other autoimmune diseases2. A number of potentially contributory genes have been suggested for GV on the basis of candidate gene association and genetic linkage studies, though few of these have received consistent support3.

We recently carried out a genome-wide association (GWA) study of GV patients and families of European descent, identifying 10 different loci that contribute to GV risk4. In addition, seven other loci showed suggestive association with GV in the initial GWA analysis (Supplementary Table 1), defined as nominal P values < 10-4 for multiple SNPs clustered across a contiguous genomic region. At 3p13, we observed suggestive association of 9 SNPs spanning nt 71,505,650-71,571,667; most significant was rs17008713 (P = 3.70 × 10-6, OR = 1.32), located within FOXP1 (Figure 1a). At 3q13.13, we observed suggestive association of 15 SNPs spanning nt 107,078,487-108,580,778; most significant was rs2603127 (P = 2.69 × 10-7, OR = 1.34), located within MYH15. At 6q27 we observed suggestive association of 17 SNPs spanning nt 166,054,922-167814784; most significant was rs6902119 (P = 5.72 × 10-5, OR = 1.21), upstream of CCR6 (Figure 1b). At 7p21.3 we observed suggestive association of three SNPs spanning nt 8,176,301-8,185,089; most significant was rs2192346 (P = 6.59 × 10-5, OR 1.26), located within ICA1. At 9q22.33 we observed suggestive association of 10 SNPs spanning nt100,951,838-101,049,252; most significant was rs7868451 (P = 8.37 × 10-5, OR = 1.22), located within TBC1D2. At 12q13.2 we observed suggestive association of 6 SNPs spanning nt 56,368,078-56,491,880; most significant was rs1701704 (P = 1.66 × 10-7, OR = 1.30), located upstream of IKZF4. At 12q24.12 we observed suggestive association of 18 SNPs spanning nt 110,557,312-113,039,943; most significant was rs3184504 (P = 6.91 × 10-6, OR = 1.24), located within SH2B3.

Figure 1

Newly replicated associations in GV. Upper panel shows genomic control-corrected PLINK association results from the GWA scan for genotyped (black) and imputed (blue) SNPs on the y axis versus chromosomal nucleotide position (GRCh37) on the x axis surrounding (a) FOXP1 and (b) CCR6. Red squares indicate Cochran-Mantel-Haenszel combined P-values for the most strongly associated SNP in each locus. LD patterns for SNPs across the regions are shown below. Arrows indicate gene positions and transcriptional orientation. Lower panel shows pairwise r2 values for LD; darker boxes indicate stronger disequilibrium for SNPs in the upper panel.

To test replication of association of these seven candidate signals, we genotyped the most significant SNP in each locus (Supplementary Table 2) in two independent replication cohorts of European descent: Replication 1 consisted of 647 unrelated GV cases and 1056 non-GV controls (principally spouses of GV cases) and Replication 2 consisted of 183 simplex GV trios and 332 GV multiplex families. SNP rs17008713 could not be genotyped for technical reasons; accordingly, we imputed genotypes for nearby SNP rs17008723 (imputed genotype r2 = 0.995), which is in almost complete LD with rs17008713 in the GWA dataset (r2 = 0.99), and in the replication study we therefore genotyped SNP rs17008723. SH2B3 region SNP rs3184504 and IKZF4 region SNP rs1701704 deviated significantly from Hardy-Weinberg equilibrium in the Replication 1 controls; therefore, these SNPs were excluded from further analysis. Case-control association statistics were calculated using the Cochran-Armitage trend test implemented in PLINK5 and family-based association statistics were calculated using FBAT6. We carried out a combined meta-analysis of the two replication studies using a Cochran-Mantel-Haenszel test with cases and controls from Replication 1, and cases and pseudocontrols derived from Replication 2, and an overall combined meta-analysis of the two replication studies plus the GWA study. We considered as conservative joint criteria for confirmed association: 1) identification of the same high-risk allele in the GWA study and both of the replication studies; 2) nominally significant (P < 0.05) association in both replication cohorts or significant association in one and marginal association in the other replication cohort; 3) significant (Bonferroni uncorrected P < 1.0 × 10-2; 0.05 / 5) combined replication stage 1 + 2 P value; and 4) a genome-wide significant (P < 5 × 10-8) (ref. 7) overall combined P value.

As shown in Table 1, we replicated association for FOXP1 region SNP rs17008723 (risk allele G; combined replication stage 1 + 2 P = 1.36 × 10-3; overall combined P = 1.04 × 10-8, OR = 1.33; Figure 1a) and for CCR6 upstream SNP rs6902119 (risk allele C; combined replication P = 3.79 × 10-3; overall combined P = 3.94 × 10-7, OR = 1.23; Figure 1b), although the latter did not achieve the genome-wide significance threshold in the combined analysis. MYH15 SNP rs2603127 showed a consistent high-risk allele across the three study cohorts, and near genome-wide significant association in the combined analysis (P = 5.36 × 10-8), but showed marginal or no association in the replication cohorts and marginal combined replication P values and so was considered unconfirmed. SNPs at two loci, ICA1 and TBC1D2, failed to even show consistent high-risk alleles across the three study cohorts and thus are considered most likely excluded.

Table 1

Replication analysis of novel candidate GV susceptibility loci

		GWA Study				Replication 1		Replication 2		Meta-analysis Replication 1 + Replication 2		Meta-analysis GWA + Replication 1 + Replication 2
SNP	Risk allele	AFRA, cases	AFRA, controls	PLINK P	OR	P	OR	P	OR	P	OR	P	OR
Replicated
rs17008723	G	0.214	0.171	3.65 × 10-6	1.32	1.47 × 10-3	1.35	0.017	1.26	1.36 × 10-3	1.36	1.04 × 10-8	1.33
rs6902119	C	0.495	0.446	5.72 × 10-5	1.21	2.92 × 10-3	1.28	0.058	1.09	3.79 × 10-3	1.25	3.94 × 10-7	1.23
Uncertain
rs2603127	A	0.247	0.213	2.69 × 10-7	1.34	0.404	1.09	0.078	1.00	0.085	1.19	5.36 × 10-8	1.30

SNPs rs17008723, rs2603127, and rs6902119 were located at nt 71,573,135 in the FOXP1 region of 3p13; nt 167,505,791 in the CCR6 region of 6q27; and nt 108,243,551 in the MYH15 region of 3q13.13, respectively. SNP nucleotide positions are from GRCh37 and genes in close proximity to the designated SNP are denoted. AFRA, allele frequency of the risk allele. EIGENSTRAT GWA P-values for SNPs rs17008723, rs6902119, and rs2603127 were 6.25 × 10-6, 7.23 × 10-5, and 4.50 × 10-7, respectively. PLINK and EIGENSTRAT GWA test statistics were calculated and adjusted for the genomic inflation factor 1.048 as described in the Supplementary Methods. The Bonferroni adjusted significance threshold for the combined replication stage 1 + 2 meta-analysis was P < 1.00 × 10-2, and the significance threshold for the overall combined GWA + replication stage 1 + 2 meta-analysis was P < 5 × 10-8. FOXP1 SNP rs17008713 from the GWA study could not be genotyped in the replication study for technical reasons; therefore, genotypes for nearby SNP rs17008723 (r2 = 0.99) were imputed, and rs17008723 was genotyped in the replication study. Data for SNPs rs2192346 and rs7868451, for which association was excluded, are shown in Supplementary Table 3.

FOXP3 and CCR6 both encode proteins that play important roles in immune regulation. FOXP1 encodes a widely expressed transcription factor that is essential for the development of B-cells8, quiescent naïve T-cells9, and monocytes10, and is paralogous to FOXP3, which encodes a transcriptional regulator of regulatory T-cell development and function and is the defective gene in the IPEX multiple autoimmune disease syndrome11. Foxp1 conditional knock-out mice have premature cell-autonomous hyper-activation of early thymocytes to CD4+ and CD8+ mature T-cells with effector functions9.

CCR6 encodes a receptor for macrophage inflammatory protein-3a (CCL20) expressed on unactivated memory B- and T-cells, T-helper 17 cells, and some dendritic cells, and plays a key role in B-cell differentiation and tissue specific migration of dendritic and T cells during inflammatory and immunological responses12. Another CCR6 SNP, rs2301436, has been associated with inflammatory bowel disease13, with the same high-risk allele observed in the GV GWA dataset (P = 2.27 × 10-4). SNPs rs6902119 and rs2301436 are in moderate LD (r2 = 0.64), and logistic regression analysis indicated that association of GV with rs2301436 might be secondary to LD with rs6902119 (P = 0.0712). Recently, we described a small GWA study of GV in an isolated Romanian founder population14, identifying association in that group with another SNP in 6q27, rs13208776, located 1.44 Mb distal to rs6902119, within SMOC2. We find no association between GV and SNPs in the vicinity of rs13208776 in the present study, and no apparent long-range LD between rs6902119 and rs13208776 (r2 = 0); nevertheless, we cannot exclude the possibility that in the Romanian founder population a variant near rs13208776 might influence expression of CCR6 at a distance. Interestingly, these vitiligo-associated 6q27 SNPs are in close proximity to IDDM8 (http://www.t1dbase.org), a locus with linkage and association with type I diabetes15-18 and rheumatoid arthritis19, autoimmune diseases that are epidemiologically associated with GV2.

Our findings provide additional evidence that variation in genes encoding proteins with roles in the immune system contribute to susceptibility towards GV. Moreover, many of these genes also contribute to other autoimmune diseases, particularly those with which GV is epidemiologically associated. While each of the GV susceptibility loci thus far identified accounts for only a small increase in relative risk, the biological pathways they highlight provide insights into the pathogenesis of GV and other autoimmune diseases, and may afford relatively tractable targets for therapeutic intervention.