Genome-wide association study in almost 195,000 individuals identifies 50 previously unidentified genetic loci for eye color.

Human eye color is highly heritable, but its genetic architecture is not yet fully understood. We report the results of the largest genome-wide association study for eye color to date, involving up to 192,986 European participants from 10 populations. We identify 124 independent associations arising from 61 discrete genomic regions, including 50 previously unidentified. We find evidence for genes involved in melanin pigmentation, but we also find associations with genes involved in iris morphology and structure. Further analyses in 1636 Asian participants from two populations suggest that iris pigmentation variation in Asians is genetically similar to Europeans, albeit with smaller effect sizes. Our findings collectively explain 53.2% (95% confidence interval, 45.4 to 61.0%) of eye color variation using common single-nucleotide polymorphisms. Overall, our study outcomes demonstrate that the genetic complexity of human eye color considerably exceeds previous knowledge and expectations, highlighting eye color as a genetically highly complex human trait.

INTRODUCTION

Eye color is primarily determined by melanin abundance within the iris pigment epithelium, which is greater in brown than in blue eyes (), and both the density and distribution of stromal melanocyte cells (). Ratios of the two forms of melanin, eumelanin and pheomelanin, within the iris as well as light absorption and scattering by extracellular components (Tyndall scattering) are additional factors that give irises their color (). Absolute melanin quantity and the eumelanin:pheomelanin ratio are higher in brown irises (), while blue or green irises have very little of both pigments and relatively more pheomelanin.

European populations, or those with partial European origin, display the largest diversity of iris color, varying from the lightest blue to darkest brown. The prevalence of blue eyes correlates with geographic latitude across Europe and neighboring areas (), likely as a result of human migration, sexual, and possibly natural selection (–). Similarly, eye color variation with varying degrees of brown irises is seen in Asian populations (), although with a much reduced range compared to brown eye color variation in Europeans.

Iris color is highly heritable (). Previous genome-wide association studies (GWASs) have identified various single-nucleotide polymorphisms (SNPs) in and around 10 genes significantly associated with eye color (–), highlighting the polygenic nature of a trait that, in the past, was assumed to be genetically simple (). The strongest genetic influence on eye color is exerted by the neighboring HERC2 and OCA2 genes (, , , ), where a long distance enhancer effect of an intronic SNP in HERC2 was demonstrated interacting with the OCA2 promoter functioning as a molecular switch between light and dark pigmentation (). Previously available genetic knowledge allows accurate prediction of blue and brown eye color, for instance, with a DNA test system based on six SNPs from six genes, including HERC2 and OCA2 (), that has been used in anthropological () and forensic (, ) applications. However, nonblue and nonbrown eye color can be genetically predicted with considerably lower accuracy (, , ), likely because of unknown predictive SNPs and responsible genes. Moreover, the phenotypic variance of eye color not previously explained by GWASs ranged between 26% for a blue versus a brown scale () and 50% for a scale using three categories () across studies. This illustrates the likely presence of yet unknown genes responsible for the noted missing genetic predicting accuracy (, , ) and missing heritability of human eye color ().

To overcome these limitations and better understand the genetics of human eye color, we carried out the largest eye color GWAS to date. Our study involved 157,485 individuals of European ancestry in the discovery stage and an additional 35,501 ancestral European individuals in the replication stage, as well as 1636 Asians (of Han Chinese and Indian ancestry) additionally used for replication purposes, in total 194,622 individuals of different ancestries.

RESULTS

Linear regression results were obtained for 11,532,091 SNPs from 157,485 individuals of European ancestry in the discovery dataset from the personal genetics company 23andMe Inc. The Devlin’s genomic factor () was λ = 1.13, while a linkage disequilibrium (LD) score regression intercept of 1.095 with an [(intercept − 1)/(mean(χ2) − 1)] ratio of 0.19 was obtained, which is consistent with expectations based on the large sample size and polygenic architecture (). In total, 12,192 SNPs were associated at genome-wide significance (P < 5 × 10−8) with eye color. These SNPs clustered in 52 distinct genomic regions, 50 across the autosomes and 2 on chromosome X (Table 1 and Fig. 1). While confirming the 10 genomic regions previously associated with human eye color (, , ), the remaining 42 genomic regions identified represent novel discoveries.

Table 1

Strongest associated SNPs from the 52 genomic regions independently associated with eye color in the European discovery cohort (N = 157,485).

A full list of 115 independently associated SNPs from conditional analysis at these loci can be found in table ST1. SNPs with novel associations for eye color identified here are highlighted in bold. SNPs that were associated with other pigmentation traits, but not iris color, in previous studies are indicated by an asterisk (*). Chr is the chromosome for the given SNP, genome position (Pos) refer to the HG build 37, RS is the rsid for the given SNP, Freq is the allele frequency of the stated reference allele, Beta is the effect size for the stated reference allele, SE is the respective standard error for the beta, Alt allele is the alternative allele for the given SNP, and N is the sample size tested for the respective SNP.

Region	Chr	Pos	RS	Nearest gene	Ref allele	Alt allele	N	Freq	Beta	SE	P
1	1	9166344	rs6693258	GPR157	T	C	157483	0.216	−0.056	0.0091	9.86 × 10−10
2	1	42110888	rs6696511	HIVEP3	T	C	157483	0.382	0.109	0.0077	1.35 × 10−45
3	1	212421629	rs351385	*DTL	G	A	157483	0.580	−0.083	0.0075	2.31 × 10−28
4	1	236035805	rs2385028	LYST	T	C	157483	0.245	0.061	0.0086	2.43 × 10−12
5	2	46233381	rs13016869	PRKCE	G	C	157483	0.208	−0.058	0.0092	2.23 × 10−10
6	2	206950236	rs112747614	INO80D	G	C	157483	0.696	0.054	0.0081	3.69 × 10−11
7	2	219755011	rs121908120	WNT10A	T	A	157483	0.976	0.140	0.0247	1.54 × 10−8
8	2	223483670	rs12614022	FARSB	G	A	157483	0.636	0.048	0.0078	4.26 × 10−10
9	2	239276278	rs74409360	TRAF3IP1	T	C	157483	0.076	−0.128	0.0140	8.09 × 10−20
10	3	42762488	rs3912104	CCDC13	T	A	157483	0.541	0.045	0.0075	1.47 × 10−9
11	3	69980177	rs116359091	*MITF	G	A	157483	0.974	−0.229	0.0246	1.02 × 10−20
12	4	23939399	rs4521336	PPARGC1A	T	C	157483	0.359	0.048	0.0079	1.28 × 10−9
13	4	59359559	rs141318671	No gene	G	A	157483	0.993	−0.318	0.0568	2.09 × 10−8
14	4	90059434	rs6828137	TIGD2	T	G	157483	0.457	−0.041	0.0075	3.93 × 10−8
15	5	311902	rs62330021	PDCD6, AHRR	G	A	157483	0.948	−0.325	0.0171	1.36 × 10−80
16	5	33951693	rs16891982	*SLC45A2	G	C	157483	0.958	−0.654	0.0191	1.97 × 10−255
17	5	40273518	rs348613	DAB2	G	A	157483	0.001	0.731	0.1227	2.54 × 10−9
18	5	123896988	rs72777200	ZNF608	T	C	157483	0.828	0.074	0.0100	1.22 × 10−13
19	5	148216187	rs11957757	ADRB2	G	A	157483	0.547	0.042	0.0075	2.94 × 10−8
20	6	396321	rs12203592	*IRF4	T	C	157483	0.171	−0.385	0.0100	1.61 × 10−321
21	6	10538183	rs6910861	GCNT2	G	A	157483	0.461	0.078	0.0076	1.20 × 10−24
22	6	158841725	rs341147	TULP4	G	A	157483	0.343	−0.045	0.0078	6.45 × 10−9
23	7	45960645	rs2854746	IGFBP3	G	C	157483	0.587	−0.043	0.0076	1.24 × 10−8
24	7	83653553	rs6944702	SEMA3A	T	C	157483	0.512	0.048	0.0075	1.71 × 10−10
25	8	12690997	rs6997494	LONRF1	T	G	157483	0.480	0.048	0.0074	1.07 × 10−10
26	8	42003663	rs12543326	*AP3M2	G	A	157483	0.875	0.065	0.0113	9.07 × 10−9
27	8	81350433	rs147068120	ZBTB10	T	C	157483	0.038	−0.122	0.0198	8.34 × 10−10
28	9	12677471	rs13297008	*TYRP1	G	A	157483	0.388	0.260	0.0077	4.99 × 10−250
29	9	27366436	rs12552712	MOB3B	T	C	157483	0.604	−0.056	0.0077	2.79 × 10−13
30	9	132001056	rs12335410	IER5L	T	C	157483	0.221	−0.058	0.0090	1.28 × 10−10
31	11	68831364	rs72928978	*TPCN2	G	A	157483	0.877	0.142	0.0121	6.42 × 10−32
32	11	89017961	rs1126809	*TYR	G	A	157483	0.721	0.285	0.0083	1.82 × 10−255
33	12	23979791	rs9971729	*SOX5	C	A	157483	0.566	0.082	0.0076	5.03 × 10−27
34	12	92567833	rs790464	BTG1	T	C	157483	0.189	0.060	0.0096	4.50 × 10−10
35	13	74178399	rs2095645	KLF12	G	A	157483	0.672	0.050	0.0079	3.26 × 10−10
36	13	95189401	rs9301973	*DCT	G	A	157483	0.671	0.066	0.0080	1.48 × 10−16
37	14	69236136	rs138777265	ZFP36L1	C	A	157483	0.047	−0.109	0.0176	6.21 × 10−10
38	14	92780387	rs17184180	*SLC24A4	T	A	157483	0.567	0.271	0.0075	1.34 × 10−284
39	15	28211758	rs4778218	*OCA2	G	A	157483	0.834	−0.520	0.0100	<1 × 10−330
40	15	28356859	rs1129038	*HERC2	T	C	157483	0.734	−2.605	0.0058	<1 × 10−330
41	15	48426484	rs1426654	SLC24A5	G	A	157483	0.102	0.684	0.0644	2.72 × 10−26
42	16	1376386	rs761063	UBE2I	G	C	157483	0.277	−0.050	0.0084	2.07 × 10−9
43	17	1966889	rs4790309	SMG6	T	C	157483	0.456	0.041	0.0075	3.71 × 10−8
44	17	67497367	rs3809761	MAP2K6	G	A	157483	0.811	0.059	0.0101	5.08 × 10−9
45	17	79612397	rs6420484	*TSPAN10	G	A	157483	0.645	−0.161	0.0080	7.53 × 10−90
46	19	7581625	rs73488486	ZNF358	T	G	157483	0.097	0.086	0.0128	2.40 × 10−11
47	20	4948248	rs2748901	SLC23A2	G	A	157483	0.463	−0.047	0.0075	3.37 × 10−10
48	21	38568882	rs2835660	TTC3	G	C	157483	0.545	−0.064	0.0075	1.62 × 10−17
49	21	44783287	rs622330	*SIK1	G	A	157483	0.512	0.111	0.0075	3.87 × 10−49
50	22	46369657	rs35051352	WNT7B	G	C	147875	0.549	−0.062	0.0081	3.33 × 10−14
51	23	48118832	rs78542430	SSX1	T	A	157483	0.450	−0.095	0.0062	5.93 × 10−53
52	23	119439335	rs5957354	TMEM255A	T	A	157483	0.861	−0.090	0.0092	5.43 × 10−23

Fig. 1

Manhattan plot of eye color GWAS results in the European discovery cohort (N = 157,485).

All results with P values <5 × 10−8 are indicated in red.

As expected, many of the strongest associations were observed for SNPs within HERC2 (rs1129038, P < 10−330), the strongest eye color–associated region known previously, and genes involved in five of the seven known types of oculocutaneous albinism (OCA): OCA1 (TYR, rs1126809, P = 1.8 × 10−255), OCA2 (OCA2, rs1800407, P < 10−330), OCA3 (TYRP1, rs13297008; P = 5.0 × 10−250), OCA4 (SLC45A2, rs16891982; P = 2.0 × 10−255), and OCA6 (SLC24A5, rs1426654; P = 2.7 × 10−26). HERC2, TYR, OCA2, TYRP1, and SLC45A2 have previously been associated with eye color (, , , , ), while previous studies found association with hair and skin pigmentation for SLC24A5 () and eye color in a South Asian population (), but not with eye color in Europeans as we showed here. No significant eye color associations were found for either the OCA5 (the 4q24 region) or the OCA7 (C10ORF11) locus, the other two genes involved in OCA, in our large discovery sample; however, C10ORF11 has recently been associated with human eyebrow color (). The remaining four previously reported loci () were also strongly associated with eye color in the present study: LYST (rs2385028, P = 2.4 × 10−12), IRF4 (rs12203592, P = 1.6 × 10−321), SLC24A4 (rs17184180, P = 1.3 × 10−284), TSPAN10 (previously reported as NPLOC4, rs6420484; P = 7.5 × 10−90), and TTC3 (rs2835660, P = 1.6 × 10−17).

Among the 41 novel genetic loci identified in our discovery analysis, we detected significant eye color association for SNPs within TPCN2 (rs72928978, P = 6.42 × 10−32) and MITF (rs75114713, P = 8.11 × 10−24). These two genetic loci and five others (near the DTL, AP3M2, SOX5, DCT, and SIK1 genes; Table 1) were associated with hair and skin pigmentation in previous GWASs (highlighted with asterisk in Table 1) (, –), while their eye color association is reported here. A comparison of our results with those from the GWAS Catalog () revealed that 34 (81%) of the 41 novel eye color–associated loci are unique to eye color and not shared with other pigmentation traits (Table 1), which is 36 (69%) when considering all 52 eye color–associated loci we identified here including the 11 previously known loci. The unique novel eye color loci include SNPs in TRAF3IP1 (rs74409360, P = 8.1 × 10−20) and SEMA3A (rs6944702, P = 1.7 × 10−10). TRAF3IP1 has been previously associated with iris furrows, while SEMA3A was associated with iris crypt variation (). These findings suggest that eye color may, at least in part considered by the eye color phenotyping done here, be mediated through structural effects within the iris. Significant association was also observed for SNPs within HIVEP3 (rs6696511, P = 1.3 × 10−45), which was previously associated with refractive error and myopia (). Novel eye color association within PPARGC1A (rs4521336, P = 1.3 × 10−9) and within MAP2K6 (rs3809761, P = 5.1 × 10−9) likely arises from their participation in pigmentation pathways, the PGC1-α protein, the product of the PPARGC1A gene. PGC1-α activates the MITF promoter () acting as a regulator in the pigment pathway of the tanning response, while overexpression of MAP2K6 increases melanocyte dendricity (), allowing greater transportation of melanosomes. We report here eye color associations for genetic loci located on the X chromosome. DNA variants clustering around the SSX1 (rs78542430, P = 5.93 × 10−53) and TMEM255A (rs5957354, P = 5.43 × 10−23) genes were significantly associated with eye color. Little is published about the TMEM255A and SSX1 genes and their potential functional roles relating to the eye, and further studies are required to understand the role of these genes in iris pigmentation. Notably, a recently published large GWAS on hair color in Europeans was the first to discover X-chromosomal loci to be involved in human pigmentation (), albeit with another gene (COL46A, 12 Mbp upstream from nearest associated locus in our study) than we found to be associated with eye color here. Notably, SNPs in or near the MC1R gene, known to be involved in light skin and red hair color, did not show eye color association in our study.

Next, we conducted a conditional analysis to identify SNPs associated with eye color after adjusting for the effect of the main association signal at a respective locus. This analysis highlighted 115 conditionally associated SNPs distributed across the novel and previously known 52 genomic regions (table ST1; a secondary list of 489 high-quality control (QC) scoring, significantly associated SNPs is provided in table ST2) that improve the phenotypic variability explained by genetic factors compared to the lead SNPs. Of the 115 independently associated SNPs, 9 were missense mutations (table ST3) (), which may suggest their causal role.

Using additional independent data from 35,501 individuals of European ancestry enrolled in nine different studies collected by the International Visible Trait Genetics (VisiGen) Consortium and its study partners, we sought to replicate the findings from the discovery analysis. For this, we tested the strongest associated SNP for each of the 52 regions (lead SNPs) identified in the discovery GWAS. To minimize data fragmentation and therefore maximize statistical power, results of regression analyses from each of the nine VisiGen studies were meta-analyzed together, and the outcome of this meta-analysis was used to replicate the discovery stage findings. Because of the differences in sample size, genotyping platforms, and imputation methods across the VisiGen studies, several SNPs highlighted by the discovery GWAS were not available or did not pass quality control in the replication dataset, particularly variants of low minor allele frequency (MAF). Some VisiGen studies had no availability of X chromosome data. Therefore, replication was attempted for 48 lead SNPs from 50 autosomal regions highlighted in the discovery analysis (table ST4). Overall, for 47 (98%) of the 48 lead SNPs, the direction of effect was the same in the replication analysis as previously seen in the discovery analysis (Fig. 2 and table ST4). For 27 (56%) lead SNPs, we obtained significant replication after applying a strict Bonferroni adjustment for multiple testing (0.05/48 = 1.04 × 10−3), including 18 (47%) of the 38 novel lead SNPs. An additional nine (24%) of the novel lead SNPs were associated at nominal (P < 0.05) significance. There was minimal heterogeneity between cohorts for these SNPs; a histogram of the heterogeneity I-scores is provided in fig. S1.

Fig. 2

Comparison of SNP beta and Z scores between the European discovery (23andMe, N = 157,485) and the European replication (VisiGen, N = 35,501) cohorts.

Color codes represent significance in the replication cohort: red = P < 1.04 × 10−3; green = P < 0.05; blue = P > 0.05.

Results from both the discovery and replication studies were then meta-analyzed including a total of 192,986 Europeans from 10 populations. The Devlin’s genomic factor was relatively unchanged in this analysis (λ = 1.14). With the added statistical power, we found genome-wide significant associations for SNPs in an additional nine separate genomic regions (table ST5). All nine genetic loci were novel, not previously associated with human eye color. Three of these, PDE4D (rs62370541, P = 5.8 × 10−9), JAZF1 (rs849142, 3.0 × 10−9), and SOX6 (rs2351061, P = 4.0 × 10−8), have recently been associated with hair color (). PDE4D inhibitors have been shown to increase melanin pigment in the skin of mouse models (), as PDE4D is a target of the melanin-stimulating hormone/cyclic adenosine monophosphate (AMP)/melanocyte-inducing transcription factor (MITF) pathway ().

Next, we compared results obtained from European subjects with association observations from a meta-analysis of 1636 Asian individuals (959 Han Chinese and 677 Indians from Singapore), for which quantitative measurements of the iris color were available (see the Supplementary Methods for information about quantitative phenotyping). Reliable SNP data were available for 44 of the lead autosomal SNPs from the 52 genetic loci identified by conditional analysis in the European discovery cohort. The remaining eight lead SNPs had a MAF smaller than 1% in this cohort and thus were excluded on account of the much smaller sample size. Thirty-one SNPs (70%) had the same direction of effect as the European analysis, and 5 (11%) of these 44 SNPs were significantly associated with eye color in Asians, after adjustment for multiple testing (Bonferroni-adjusted P value: 0.05/44 = 1.13 × 10−3). This included SNPs from two of the newly identified genes GPR157 (rs6693258, P = 3.6 × 10−5) and SIK1 (rs622330, P = 1.7 × 10−4) (table ST6), as well as from the previously known HERC2 gene (rs1129038, P = 2.2 × 10−4), which is the most strongly eye color–associated SNP in Europeans (). For this marker, we observed considerable MAF differences between Europeans and Asians (T allele: 0.05 in Asians and 0.73 in Europeans). The strongest association in the Singaporean cohorts, however, was found within SLC24A5 (rs1426654, P = 1.9 × 10−7) representing the OCA type 6 (OCA6) locus. Notably, OCA6 was first reported in a Chinese family (). However, the HERC2 and SLC24A5 SNPs were the only two variants to display heterogeneity between the two Singaporean cohorts (table ST6), with association primarily being driven by the cohort of South Asian ancestry. This is likely a due to South Asian and European populations being more closely related than East Asian and European populations. The MAFs for rs1426654 were considerably different between Asians (0.80 overall, 0.99 in East Asians, and 0.52 in South Asians) and Europeans (0.10) for the G allele, explaining why this SNP is often used as a DNA marker for biogeographic ancestry (). Previously, SLC24A5 rs1426654 not only explained a considerable proportion of the variation in skin pigmentation between different continental populations () but also was associated with skin and eye color variation within a South Asian population (, ). Despite such drastic differences in allele frequency for some SNP alleles, the shared eye color effect between populations from two different continents demonstrates the value of our multiethnic study.

DNA variants identified via GWAS are markers of statistical association and not necessarily causative. Because of the differences in LD across continental populations, association may not be universally replicable. Therefore, we next assessed the presence of association not just for the lead SNPs in each region but also by testing other SNPs located within the 52 genomic regions identified in the European discovery analysis. Several SNPs within these regions (table ST7) showed evidence for eye color association in both European and Asian populations, with genome-wide significance in Europeans and suggestive levels of genome-wide association (P < 1 × 10−5) in Asians, despite a much smaller sample size. These SNPs also showed no significant heterogeneity between the two Asian cohorts. It is therefore possible that, against the background of much stronger effects of European-only alleles or due to population-specific differences in MAF, some polymorphisms contributing to European eye color variation are also relevant for eye color variability in non-European populations, such as variation of brown eye color in Asians tested here.

In Europeans, the 112 autosomal SNPs identified through conditional analysis (all autosomal SNPs shown in table S1) explained 99.96% (SE = 6.5%, P = 4.8 × 10−279) of the liability scale for blue eyes (against brown eyes) and 38.5% (SE = 5.7%, P = 2.2 × 10−130) for intermediate eyes in the TwinsUK cohort, which was one of the VisiGen cohorts used for replication. Using the same linear scale as the GWAS analysis, these autosomal SNPs explained 53.2% (SE = 4.0%, P = 1.2 × 10−322) of the total phenotypic variation in eye color in TwinsUK.

Last, we performed in silico analyses to explore the putative function of the genetic loci our study highlighted with significant eye color association using the conditional SNPs. Gene set enrichment analysis identified multiple pathways with significant enrichment (table ST8). As expected, this included several pigmentation process pathways, with “Developmental Pigmentation” the most significant (P = 7 × 10−6), followed by “Frizzled Binding” (P = 2.0 × 10−4) and “Melanin Metabolic Process” (P = 2.6 × 10−4). We also examined the potential effects of the identified eye color–associated SNPs on gene expression using data from the GTEx Consortium (). Despite the lack of iris tissue in the GTEx repository, many SNPs showed significant eQTL (expression quantitative trait loci) effects in multiple cell types and tissues (), as seen for the associated SNPs across 38 (79%) of the 48 tissues in the GTEx dataset (table ST9). Most of the strongest effects were seen in nerve and sun-exposed skin tissue (P = 8.47 × 10−62 and P = 4.84 × 10−49, respectively) for rs2835660, where the C allele is significantly associated with a decrease in TTC3 expression. TTC3 is in proximity to DSCR9, a gene whose polymorphisms were previously associated with eye color (). These results implicate TTC3 as a more likely candidate gene influencing eye color at this genetic locus than DSCR9 and that its effect on eye color is likely mediated through variation in gene expression. The lack of iris tissue information in GTEx likely explains the absence of stronger eQTL effects, such as for SNPs with regulatory effects over gene transcription ().

DISCUSSION

We report the results of the largest GWAS for human eye color to date. In addition to confirming the association of SNPs in 11 previously known eye color genes (, , , , ), the identification of 50 novel eye color–associated genetic loci helps explain previously missing heritability of eye color variability in European populations. Moreover, because of the multiethnic design of our study, we demonstrate that several of the genetic loci discovered in Europeans also have an effect on eye color in Asians.

Eight of the genes in or near the loci newly associated with eye color in our study were previously reported for genetic associations with other pigmentation traits, such as hair and skin color, for instance, TPCN2, MITF, and DCT (, , , ). The commonality of associated DNA variants across the three pigmentation traits helps explain why the different pigmentation traits frequently (but not completely) intercorrelate in European populations. While many significant genetic associations are shared between iris color and other pigmentation traits, there are also notable differences. Although DNA variants within the MC1R gene are strongly associated with light skin and red hair color (), no detectable association with eye color was found in our large GWAS, in line with previous albeit smaller-sized GWASs of more limited statistical power (, , ). Similarly, other DNA variants strongly associated with skin and hair color within genes, such as SILV, ASIP, and POMC (), showed no statistically significant effect on eye color in this study, nor in previous studies. Moreover, we also identified 34 genetic loci that were significantly associated with eye color, but for which there is no report of significant association with hair and/or skin color. This is remarkable as the statistical power of the recent GWASs on hair color (, ) and sun sensitivity () were similar to that of our current eye color GWAS. Significant associations for SNPs in/near genes involved in iris structure, such as TRAF3IP1 and SEMA3A, suggest that they exert their effects with changes in Tyndall scattering, rather than through alterations of melanin metabolism. Overall, this demonstrates that although many genes overlap between eye, hair, and skin color, the different human pigmentation traits are not completely determined by the same genes as we showed.

The major strengths of our study compared with previous eye color GWASs arise from the larger sample size, which translated into increased statistical power and also the ability to lower the threshold of MAF for which sufficient power to detect association is available. Rare SNPs are often a source of considerable phenotypic variation (). For instance, seven (6%) of the independently associated SNPs identified by conditional analysis in the discovery cohort had a MAF between 0.1 and 1%. Despite their low frequency, however, five (71%) of these rare SNPs were in the same region as other, more common conditional SNPs that did replicate. The remaining two loci (DAB2 and an intronic region on chromosome 4) that were not formally replicated should therefore be considered only as strong candidates with respect to their association with eye color, pending independent validation in future studies.

Another strength of this work is the inclusion of European and non-European populations. Non-European populations are underrepresented in the GWAS literature in general, including in pigmentation GWASs, but their study is important for the understanding of the genetic basis of human phenotypes (). Although eye color variation is typically attributed to individuals of (at least partial) European descent, or those originating from areas nearby Europe, more subtle variation in brown eyes is also observed in Asian populations without European admixture (). Our results from the Asian cohorts showed remarkable consistence in the genetic architecture of eye color among individuals of different continental ancestries with Asian replication for the two major European genes OCA2 and HERC2. Moreover, our findings also suggest that while a single regulatory variant in HERC2 is responsible for most blue/brown variation in Europeans (), many additional DNA variants across both OCA2 and HERC2 seem to have independent effects. This hypothesis is further supported by our conditional analysis in the European discovery cohort, identifying independent associations spanning ~14 mbp across both genes rather than a concentrated cluster centered at HERC2 rs1129038. This is remarkable given the large eye color variation from the lightest blue to the darkest brown in Europeans, compared with the more limited variation within brown eye color in Asians.

In conclusion, our work has identified numerous novel genetic loci associated with human eye color in Europeans, of which a subset also shows effects in Asians, despite their largely reduced phenotypic eye color variation compared with Europeans. The genetic loci we identified explain the majority (53.2%) of eye color phenotypic variation (classified using a three-category scale) in Europeans and a large proportion of the previously noted missing heritability of eye color. Our findings clearly demonstrate that eye color is a genetically highly complex human trait, similar to hair () and skin color (), as highlighted recently in large European GWASs. The large number of novel eye color–associated genetic loci identified here provide a valuable resource for future functional studies, aiming to understand the molecular mechanisms that explain their eye color association, and for future genetic prediction studies, aiming to improve DNA-based eye color prediction in anthropological and forensic applications.

METHODS

We performed a GWAS using information from two sources: 157,485 research participants of European ancestry recruited among the customer base of the personal genomics company 23andMe Inc. (Sunnyvale, CA, USA) () and a meta-analysis of 35,501 European individuals from nine populations collected by members of the International Visible Trait Genetics (VisiGen) Consortium () and their study partners. As a secondary analysis and for comparative purposes, we included an additional set of 1636 individuals of Asian ancestry from two populations (959 Han Chinese and 677 Indians from Singapore).

Populations and participants

23andMe

Participants

Participants included for this analysis were 23andMe consumers who consented to participate in research and were determined to have greater than 97% European ancestry using genotype clustering. A segmented identity-by-descent algorithm () was then applied to obtain the maximal subset (157,485) of unrelated individuals to be used as our discovery cohort (full details on ancestry and relatedness calculations are provided in the Supplementary Methods, and principal component plots are provided in fig. S2).

Genotyping

Participants were genotyped on one of 23andMe’s own custom-designed SNP arrays. The V4 platform is a fully customized array containing a subset of SNPs from previous versions, while the V3 is a customized platform based on the OmniExpress+ BeadChip, and the V1 and V2 platforms are customized variants of the Illumina HumanHap550+ BeadChip (full descriptions of platforms and genotyping procedure are provided in the Supplementary Methods). Imputation of additional variants was completed using 1000 Genomes phase 1 as reference haplotypes ().

Phenotyping

Phenotyping for iris color in this cohort was derived from questionnaires. Participants were asked to self-categorize eye color into one of seven distinct groups ranging from blue to dark brown based on color matching (full details are provided in the Supplementary Methods). This categorization allowed greater phenotypic resolution than previous three-color categories while maintaining groups distinctive enough to reduce misclassification. These groups were converted into a numerical scale ranging from 0 (blue) to 6 (dark brown) and used as the outcome variable for linear regression–based GWAS analysis.

VisiGen consortium and study partners

Specific recruitment criteria varied for each cohort (Supplementary Materials), with participants from the United Kingdom, The Netherlands, Italy, and European descendants in the United States and Australia.

Principal components analyses examined the genetic ancestry for each cohort and confirmed each was of nonadmixed European ancestry (specific details for each cohort are provided in the Supplementary Methods).

Eye color categorization varied between cohorts, with most of the European cohorts adopting a 3- or 5-point categorical scale (details for each cohorts’ classification are provided in the Supplementary Methods).

Singaporean cohorts

Full GWAS summary statistics were available for 959 participants of Chinese descent recruited at the Singapore Polyclinic and for 677 participants in the Singapore Indian Eye Study (SINDI).

Phenotypes in these cohorts were ascertained using image analysis of digital photographs, better suited to capture variation within these populations (full details are provided in the Supplementary Methods).

Statistical analyses

Association tests

We performed an independent linear regression for each cohort under the assumption of an additive model for allelic effects, with adjustments made for age, sex, and the first five principal components. Adjustment was made also for the genotyping platform (V1 to V4), and additional measures controlling for relatedness were implemented for the VisiGen family cohorts (full details for each cohort are provided in the Supplementary Methods). The use of a linear scale may be limited by the assumption that there is an equal difference between each group; however, as eye color phenotypically follows linear scale across the color spectrum, this choice of statistical model is appropriate and in accordance with previous studies on eye color (, , ).

Meta-analyses

The QC’ed summary statistics from the VisiGen and Singaporean cohorts were pooled into two separate meta-analyses: the first was conducted using data from the nine VisiGen European cohorts, and the second using data from the two Singaporean cohorts. Both meta-analyses were performed using METAL () applying a weighted Z-score approach, given the differing phenotypic scales used across individual cohorts, with genomic control adjustments () made for the VisiGen meta-analysis (see the Supplementary Methods). A final meta-analysis was then conducted between the 23andMe cohort and the European VisiGen cohorts following the same procedures with a weighted Z score approach and genomic control adjustments.

Definition of genomic region and significance

We adopted the customary definition for a GWAS significance threshold of 5 × 10−08. For this work, an “associated region” was a genomic region demarcated by consecutive significantly associated genomic markers, separated by a nonassociated region greater than 1 million base pairs.

LD score regression

LD Hub () was used to perform LD score regression () on the summary statistics from our GWASs to test for possible inflation. This method is more suited to our study than the Devlin () genomic inflation lambda because of high trait polygenicity and the large cohort sizes ().

Conditional analysis

Conditional analysis was conducted using summary statistics with genome-wide complex trait analysis (GCTA) () following the same procedure described by Yang et al. (). Genotypic data from unrelated participants in TwinsUK were used as a reference for LD structure, a distance of 1000 kb was set as an assumption of complete linkage equilibrium, and a collinearity threshold was set at R2 > 0.9.

Variant effect prediction

The predicted effects of variants were performed using the Ensembl Variant Effect Predictor (VEP) ().

Heritability explained by the independent SNPs

Population liability scale heritability () explained by associated SNPs identified in the 23andMe discovery cohort was calculated using an unrelated sample from the TwinsUK cohort using restricted maximum likelihood analysis () under the assumption that the population prevalence (taken from the 23andMe cohort) is 30.3% for blue eyes and 38.3% for intermediate colors.

Pathway analysis

Genetic pathway analysis for iris color was performed with MAGENTA () using the 23andMe discovery SNP association results as input (full details are provided in the Supplementary Methods). Gene set definitions defined by Gene Ontology () were obtained from the Molecular Signatures Database (version MSigDB v6.1) () for this analysis.

Members of consortia and affiliations

23andMe research team

23andMe Inc., Sunnyvale, California, USA—M. Agee, A. Auton, R. K. Bell, K. Bryc, S. L. Elson, P. Fontanillas, K. E. Huber, A. Kleinman, N. K. Litterman, M. H. McIntyre, J. L. Mountain, E. S. Noblin, C. A. M. Northover, S. J. Pitts, O. V. Sazonova, J. F. Shelton, S. Shringarpure, C. Tian, J. Y. Tung, and V. Vacic.

International Visible Trait Genetics Consortium

Department of Twins Research and Genetic Epidemiology, King’s College London, London, United Kingdom—P. G. Hysi and T. D. Spector.

Department of Genetic Identification, Erasmus MC University Medical Center Rotterdam, Rotterdam, The Netherlands—M. Kayser and F. Liu.

QIMR Berghofer Medical Research Institute, Brisbane, Queensland, Australia—D. L. Duffy and N. G. Martin.

University of Queensland Diamantina Institute, University of Queensland, Brisbane, Queensland, Australia—D. M. Evans.