Novel quantitative pigmentation phenotyping enhances genetic association, epistasis, and prediction of human eye colour

Success of genetic association and the prediction of phenotypic traits from DNA are known to depend on the accuracy of phenotype characterization, amongst other parameters. To overcome limitations in the characterization of human iris pigmentation, we introduce a fully automated approach that specifies the areal proportions proposed to represent differing pigmentation types, such as pheomelanin, eumelanin, and non-pigmented areas within the iris. We demonstrate the utility of this approach using high-resolution digital eye imagery and genotype data from 12 selected SNPs from over 3000 European samples of seven populations that are part of the EUREYE study. In comparison to previous quantification approaches, (1) we achieved an overall improvement in eye colour phenotyping, which provides a better separation of manually defined eye colour categories. (2) Single nucleotide polymorphisms (SNPs) known to be involved in human eye colour variation showed stronger associations with our approach. (3) We found new and confirmed previously noted SNP-SNP interactions. (4) We increased SNP-based prediction accuracy of quantitative eye colour. Our findings exemplify that precise quantification using the perceived biological basis of pigmentation leads to enhanced genetic association and prediction of eye colour. We expect our approach to deliver new pigmentation genes when applied to genome-wide association testing.

Scientific RepoRts | 7:43359 | DOI: 10.1038/srep43359 is reflected by the stroma of the iris, and the eye colour is perceived as grey to blue through Tyndall scattering 3 . As with many traits, the nature of human eye colour variation is continuous, spanning from the lightest shades of grey or blue to the darkest shades of brown or black 4 . Dark eye colour reflects the ancestral state in humans linked to their commonly believed origin in Africa, while light eye colour is assumed to be derived; shaped by positive selection perhaps due to sexual selection during European history 5 .
Several gene-mapping studies on eye colour were previously conducted by using manually defined phenotype categories [6][7][8][9][10][11][12][13][14] , inevitably oversimplifying the continuous nature of human eye colour variation. Although this incomplete use of the underlying basics of eye colour variation reduces the thoroughness of such studies, several eye colour genes were previously identified with this simplified phenotyping approach [5][6][7][8][9][10][11][12][13] . Moreover, it provides an accurate prediction from DNA with reasonably high accuracies for at least the extreme categories of blue and brown demonstrated via the IrisPlex system 15 , a system consisting of only six single nucleotide polymorphisms (SNPs) from six genes. Consequently, eye colour was one of the first externally visible characteristic for which the concept of Forensic DNA Phenotyping (FDP) 16,17 was put into practice 15 , later followed by hair colour 18 and most recently by skin colour 19 . However, as described elsewhere 20 there is a desire to move pigmentation colour prediction from the currently applied categorical level to the continuous level. As prerequisite, this requires an understanding of the genes that determine eye colour in its fully continuous spectrum as well as methodology that allows the capture of continuous eye colour as accurately and completely as possible.
The first quantitative approach to measure eye colour was proposed by Frudakis 21 , who introduced two quantities representing the iris colour properties, i.e. the iris colour score, and the melanin index as derived from average luminosity (L) and colour reflectance values (C) from selected boxes in digital eye-imagery. The melanin index can be directly related to the amount of melanin that is known to decrease to extremely low levels (even complete absence) in blue eyes 22,23 . In the hue (H) and saturation (S) measurements introduced to eye colour quantification by Liu et al. 24 , H defines the colour itself, which can be related to the type of melanin having more red or yellow components. The S value describes the richness of a certain colour (defined in H) that is supposed to correlate with the amount of eu-or pheo-melanin. The V value is usually discarded, as it is supposed to rather represent the brightness due to different lighting conditions 24,25 . Digital quantification of eye images of thousands of Europeans using the H-S colour space and its use in a genome-wide association study allowed the identification of three new eye colour genes, not previously identified when using categorical eye colour 24 . This study clearly demonstrated the increase of power to find new genes when moving pigmentation phenotyping from the classical categorical approach to a quantitative approach. Beleza et al. 26 averaged and normalized B-G values from the RGB-space (red, green, blue, value) and proposed a T-index quantity to describe the amount of melanin per iris. Recently, the CIE-L * a * b* values have been used in place of H, S values 27 to quantify iris colour, where the L value describes the lightness, the a* the red/green, and b* the yellow/blue component, of the colour respectively. However, taking a quantity that is averaged over the iris as previous methods 21,23,24,27 have done may obscure the different mixture proportions of pigments. An alternative has been proposed by Anderson et al. 28 , which includes a clustering of the segmented iris into blue and brown pixels deriving a ratio score (PIE score). Different types of pigmentation (i.e. eumelanin or pheomelanin), however, are not distinguished with this approach.
Here, we introduce an improvement of quantitative pigmentation phenotyping based on an automated segmentation of the iris followed by a measurement of the digital equivalents of eumelanin, pheomelanin, and total absence of any pigment in the iris. This clustering is based on manually predefined (assumed) image segments that depict eumelanin, pheomelanin and nonpigmented areas. We exemplified the advantage of this approach by an empirical analysis of high-resolution eye images from over 3000 individuals sampled from seven European countries (EUREYE study). By using genotypes of 12 SNPs previously involved in human eye colour variation that we generated in the same individuals, we demonstrate the impact of this novel pigmentation phenotyping approach on genetic association, epistasis, and prediction.

Results and Discussion
Detection and segmentation of the iris in digital imagery. Prior to colour assessment, the iris needed to be segmented from the pupil and sclera. Several approaches have been proposed previously for the segmentation of the iris in digital eye imagery [28][29][30][31][32] . The utility of a certain approach depends strongly on the properties of the given image data. In the majority of imagery available to us, (i) the pupil was always centred in the middle (Fig. 1a), (ii) the iris was fully visible (those images where it was not were excluded from the analyses), and iii) eye lashes rarely overlapped with the iris. Because of these features, we followed a previously proposed 30 two-step procedure, which we implemented in Matlab (R2007a). First, we used a Canny filter 33 to distinguish the rim confining the iris and then applied the Hough transformation 34 to detect the iris circle. To reduce the number of multiple solutions, we constrained the results of the Hough transformation on those circles only that were centred in the middle of the image. As a result, we were able to maintain a high ratio of correctly segmented irides (> 90%). Falsely segmented irides, for example when the pupil was extremely dilated, or eyes were closed by chance, were curated manually. It may be that in other types of image data taken under less normalized conditions (i.e. DSR camera systems with macro lenses), other approaches 28 might be more sensitive to apply for iris segmentation.
Quantification of iris colour from digital imagery. The developed method assesses each pixel according to its digital classification of pheomelanin, eumelanin, and non-pigmentation, using a machine learning approach (Supplementary Figure S1). The colour information of each image pixel was available as a red-green-blue (RGB) triplet, which we first transformed into a hue-saturation-value (HSV) triplet 35 . We then used a support vector machine classifier 36 with a quadratic kernel to assign each pixel within the HSV space to one of the manually defined classifications representing the total absence of pigment (non-pigmentation), pheomelanin, and eumelanin. To define the distribution of the respective classes for the training of the support vector machine (prior to the sample phenotyping procedure), we used a set of 10 randomly chosen images of different eye colours to manually Scientific RepoRts | 7:43359 | DOI: 10.1038/srep43359 label areas that most obviously contain the two types of pigments and their absence (Supplementary Figure S1a). The distributions of these three types were well separated in at least one dimension from the selected training images (Supplementary Figure S1b), which implies that the assignment of a pixel from the iris into one of the three types of pigment outcomes can be achieved with good precision. Note that the clustering outcomes were robust against the choice of the kernel and choice of colour model (RGB or HSV, data not shown).
We finally considered the proportion of the clustered pixels relative to the segmented iris as our quantitative eye colour phenotype, reflecting the equivalent amounts of non-pigmented, pheomelanin, and eumelanin areas per iris (see Fig. 1c), the sum of which equals to one.
Comparing the new eye colour quantification method with previous methods. One important motivation behind quantitative measurements is to capture the information about complex phenotypes with a small number of variables. We use a correlation analysis to reveal how well different eye colour quantifications maintain the information about our estimated abundance of different types of melanin. We calculated quantitative iris colour phenotypes using five previously proposed methods: (i) the mean H, S values from the HSV space 24 (ii) the mean luminosity value and colour score 21 , (iii) the components of the L * a * b* space 27 , (iv) the PIE score 28 and (v) the T-index 26 . Table 1 provides information about the relationship between the considered eye colour quantifiers. The amount of non-pigmentation we find mostly positively correlated with b* (r = 0.93, P < 1e-300), followed by S (− 0.74, P < 1e-300). The more melanin is distributed on the iris (fewer non-pigmentation) the lower will be the b* indicating the yellow components contained in the melanin and the higher is the a* value for the red components of the melanin. Therefore we find the amount of pheomelanin was most strongly correlated with a* (r = 0.76, P < 1e-300), followed by b* (r = − 0.68, P < 1e-300). The amount of eumelanin estimated with our approach was most strongly positively correlated with L 21 , (r = 0.71, P < 1e-300), followed by S (r = 0.69, P < 1e-300). The colour score C, which summarizes the components a* and b* 21 , shows a weaker correlation with eumelanin (r = − 0.52, P = 1.1e-286) and much weaker positive correlation with pheomelanin (r = 0.14, P = 3.6e-14), which indicates a reduced utility in predicting the respective phenotypes from colour components, which is discussed in the next paragraph.
Finally, the amount of pheomelanin estimated with our approach showed only weak correlation with the estimated amount of eumelanin (r = − 0.06, P = 0.001), confirming our expectation that these two measurements represent different biological phenotypes.
Comparing quantitative eye phenotypes with categorical eye phenotypes. A good quantitative colour measure is expected to distinguish well among manually defined eye colour categories. To test for this, we plotted the manually derived eye colour categories (blue, intermediate, brown) in the different components of the described eye colour quantities (see Fig. 2). To demonstrate the quality in separating the distributions of two selected categories, we used the Hellinger distance (HD, Table 2). The HD describes the distance between two probability distributions, which is zero if both are perfectly overlapping and one in case they are entirely separated. Note that the PIE-score and T-index are only one-dimensional quantities. The manually graded eye colour categories blue and brown were best differentiated by the non-pigmentation/eumelanin-space generated by our new method (Fig. 2c, HD = 0.923), followed by the non-pigmentation/pheomelanin space from our method (Fig. 2b, HD = 0.911). The categories blue and brown can be nearly completely separated with our approach. The categories blue and intermediate were also best separated in the non-pigmentation/pheomelanin space from our method (Fig. 2b, HD = 0.743). However, brown and intermediate were best separated by the L-a* space (Fig. 2g, HD = 0.545), followed closely by the pheomelanin/eumelanin space from our approach (0.495). A large overlap of the intermediate category with both, blue and brown, was evident, likely caused by well-known difficulties in manual assignment of eye colours to the intermediate category.
Considering all pairs of categories, the colour score-luminance space performed least accurately in comparison to all other quantitative methods used. The summary of the a* and b* components most likely obscures important information about the different types of melanin. The one-dimensional PIE score and T-index

Impact of improved eye colour phenotyping on strength of genetic eye colour association.
Next, we investigated the impact of our eye colour phenotyping approach, as well as other previously used methods, on the strength of genetic eye colour association. For this, we used 12 SNPs from 11 genes highlighted in previous genetic eye colour studies mostly using categorical phenotypes 15,18,24,[37][38][39][40] . We analysed these 12 SNPs in the same 3,087 individuals from whom we used the eye images for eye colour quantification. We employed partial correlation controlling for age, sex, and the sampling population (see Methods), and report the amount of genetically explained eye colour variance as measured by R 2 ( Table 3). Note that stronger association or higher effect size 41 (as measured by R 2 in our study) is beside the allele frequency -an important factor controlling for the power of discovering an unknown causal variant 42 . Due to the SNPs HERC2 rs12913832 and HERC2 rs1129038 being in very high linkage disequilibrium (LD) (R 2 = 98%, P < 1e-300), we only describe the results for HERC2 rs12913832 (for complete results, see Table 3). Among all quantitative measurements, the phenotypic eye colour variance explained by any of these 11 SNPs was highest for the amount of non-pigmentation as measured with our new approach (HERC2 rs12913832 R 2 = 49.7%, P < 1e-300) followed by S (HERC2 rs12913832 R 2 = 40.7%, P < e-300), and eumelanin measured with our approach (HERC2 rs12913832 R 2 = 32.1%, P = 1.3e-261), and H (HERC2 rs12913832 R 2 = 26.1%, P < 7e-205). The explained variance was considerably less of a* (HERC2 rs12913832 R 2 = 15.4%, P = 1.e-300) and pheomelanin as measured with our approach (HERC2 rs12913832 R 2 = 14.6%, P = 1.8e-107). HERC2 rs12913832 (and similarly the LD SNP HERC2 rs1129038, Table 3) individually explained considerably more variation in continuous eye colour than any other SNP tested, which is in agreement with our previous study based on HS colour space 24 ). This also is in line with previous categorical eye colour studies, where HERC2 rs1129038 contributed most 12,15,38 .
In contrast to Liu et al. 22 , we could not replicate a statistically significant association of NPLOC4 rs9894429, DSCR9 rs2835630, and LYST rs3768056 (border-line significant association with S, P = 0.06) with any of the quantitative eye colour phenotypes tested including H and S used by Liu et al. 24 (Table 3). The lack of association noted here could be explained by their small effect size together with the smaller sample size used here (N = 3,087) relative to the previous study (N = 5,951) 22 . Alternative factors could be image and individual age differences between EUREYE samples used here and the on average elder Rotterdam Study samples used previously (note that the effects are measured as single contributions to the phenotype, for combined analysis see Supplementary Table S2). Age has been previously shown to be a significant predictor for eye colour 24 . Epistatic effects detected with the detailed quantitative eye colour phenotypes. Next, we studied the impact of the detailed eye colour phenotyping as achieved with our new method, as well as previously used methods, on epistatic effects between pairs of SNPs that deviate from additivity while correcting for sex and age (see Methods and Table 4). For this, we ignored the interaction between HERC2 rs1129038 and HERC2 rs12913832 due to their strong LD (R 2 = 0.98, P < 1e-300, Supplementary Table S3). We found two novel interactions between SLC24A4 rs12896399 and SLC45A2 rs16891982 strongest in pheomelanin (P = 8.7e-04) as well as LYST rs3768056 and DSCR9 rs2835630 solely evident in eumelanin (P = 2.6e-2). Moreover, we confirmed several interaction pairs reported in previous studies using quantitative or categorical eye colour phenotypes, namely: HERC2 rs12913832 and SLC24A4 rs12896399 previously observed for S 24 as well as blue vs. non-blue 43 , observed here in PIEscore (P = 2.3e-16), b* (P = 3.3e-4), S (P = 4.2e-3), non-pigmentation (P = 5.4e-13), and pheomelanin (P = 8e-9); HERC2 rs12913832 and SCL45A2 rs16891982 previously described in blue vs. non-blue 44 , observed here in S (P = 1.3e-12), pheomelanin (P = 1.2e-16), non-pigmentation (P = 1.8e-20), S (P = 1.3e-12), Colour score (P = 3.5e-7), PIEscore (P = 3.6e-21), and b* (P = 3.8e-18); HERC2 rs12913832 and TYPR1 rs1325127 previously described in hazel vs. non-hazel 43 found only in pheomelanin (P = 7.8e-4); HERC2 rs12913832 and IRF4 rs12203592 previously observed in H and S 24 , found here in H (P = 1.9e-3), a* (P = 2.6e-4), and pheomelanin (P = 3.8e-3). We exemplified the interaction between HERC2 rs12913832 and SLC45A2 rs16891982 in Supplementary Figure S2. From the marginal distribution of HERC2 rs12913832 (SNP2 in Supplementary Figure S2) the dominant effect for allele T can be observed as the presence of a T allele strongly decreases the amount of non-pigmentation. The interaction partner SLC45A2 rs16891982 shows a recessive behaviour for the G allele that increases the amount of non-pigmentation only for the homozygote genotype GG. From the combined distributions a masking effect of HERC2 rs12913832 can be observed. The effect of the homozygote GG in SLC45A2 rs16891982 becomes much more evident if HERC2 rs12913832 is homozygote for CC (R 2 = 0.037, P = 3.7e-12). The effect of SLC45A2 rs16891982 on the amount of non-pigmentation is almost neutralized with HERC2 rs12913832 heterozygote CT or homozygote TT (R 2 = 0.003 P = 1.5e-02). Thus it appears that the allelic states of SLC45A2 rs16891982 will lighten/darken only those irises that are already stated to be "blue" by HERC2 rs12913832 = CC. A similar masking effect can be observed for the interaction pair HERC2 rs12913832 and SLC24A4 rs12896399 (Supplementary Figure S3). We similarly investigated the interaction between HERC2 rs12913832 and TYRP1 rs1325137 (Supplementary Figure S4) and observe a slight reinforcing effect of TYRP1 rs1325137 on the amount of quantified pheomelanin. Having a C allele in HERC2 rs12913832 decreases the amount of pheomelanin by the state of TYRP1 rs1325137 (R 2 = 0.06, P = 7.6e-3). In contrast, the state of TYRP1 rs1325137 tends to increase the amount of pheomelanin if HERC2 rs12913832 has at least one T allele (R 2 = 0.02, P = 6.9e-2). Condensing all SNP-SNP interactions for certain quantitative eye colour estimates, the following SNP pairs were found to be most relevant: IRF4 rs12203592 × HERC2 rs12913832, SLC24A4 rs12896399 × SLC45A2 rs16891982, HERC2 rs1291383) × SLC24A4 rs12896399, HERC2 rs12913832 × TYRP1 rs1325127, HERC2 rs12913832 × SLC45A2 rs16891982. Consequently, we included them in the prediction analyses considered as a separate model.

SNP-based prediction of quantitative eye colour phenotypes.
We conducted a formal prediction analysis of continuous eye colour expressed by the various quantitative measures with (model 1, Supplementary Table S4) or without considering the SNP-SNP interaction terms previously described (model 2, Supplementary Table S5). For this, we divided our samples into a model building set (N = 2,087) and a model-validation set (N = 1,000). As prediction accuracy measure we used the percentage of phenotypic variance predicted by the respective model (mean coefficient of determination R 2 ). Model 1 (without considering interactions) explained 52% of the phenotypic variance of non-pigmented area, 16% of pheomelanin, and 33% of eumelanin (Table 5), while H and S were predicted as 24% and 42%, respectively. PIE-score can be predicted as 45%, T-index as 31%, a* and b*, as 42% and 30% respectively. Model 2 (with considering interactions) provided  Table 2. Hellinger distances between pairs of colour categories from two dimensional colour subspaces (see Fig. 2). The higher the value close to one, the better is separability of the respective clusters in the quantitative colour space. Bold values denote the components that separate pairs of categories best.
increased prediction of 55% (i.e., by 3%) and 20% (i.e., by 4%) for non-pigmentation and pheomelanin, respectively, with eumelanin remaining at 33%, while for S the prediction increases to 44% (i.e., by 1%) and remained the same for H at 24%. Hence, by including the interaction terms in the prediction model, we observed an overall increase in the predictability of our newly derived quantitative eye colour measures. The strongest increase was seen for pheomelanin (i.e., by 4%) that notably is involved in intermediate eye colours. The following interactions were mostly relevant for the prediction of pheomelanin: HERC2 rs12913832 and SLC45A2 rs16891982 (beta = 0.73, P = 6.12e-14), HERC2 rs12913832 and SLC24A4 rs12896399 (beta = 0.32, P = 5.72e-9), as well as IRF4 rs12203592 and HERC2 rs12913832 (beta = 0.39, P = 9.56e-7), (see Supplementary Table S5). Notably, when constraining the quantitative eye colour predictions on the 6 SNPs included in the IrisPlex system for categorical eye colour prediction 15 , we observed only a minor loss of information relative to the full 12-SNP model with interaction, and varying degrees of effect relative to the 12-SNP model without interaction (Supplementary Table S6). Based on the 6 IrisPlex SNPs with considering interaction (see Supplementary Table S7 for betas), we estimated 53% for the amount of non-pigmentation (2% less than the full 12-SNP model), 18% for the amount of pheomelanin (2% less), and 33% for the amount of eumelanin (same as the 12-SNP model). The observed slight increase of predictability of pheomelanin with the 12-SNPs relative to the 6 IrisPlex SNPs underlines the added value of the additional 6 SNPs in predicting non-blue and non-brown eye colours, which is notably lower than the gain of power by our novel quantification method.
To demonstrate the potential of the 12 SNPs for predicting quantitative eye colour, in Fig. 3 we plotted the observed eye colour phenotypes as measured with our new approach, and their DNA-predictions from 12-SNP genotypes. We already find a very good prediction of the proportions of melanin-types, which is more informative with respect to iris colour than a plain category. Hence, intermediate iris colours can be visually predicted from DNA by the quantitative amounts of pheomelanin, and eumelanin, which would provide a better overall representation of eye colour rather than one single intermediate category. Knowledge about the distribution of the types of melanin on the iris will allow more realistic prediction of iris colour and structure, which can be hypothetically achieved by applying this approach to sub-areas of the iris as applied in Edwards et al. 27 .

Conclusions
In summary, we have developed an automated computational approach that separates the iris from digital eye images and quantifies iris pigmentation by estimating digital equivalents of eumelanin, pheomelanin, and non-pigmentation. By applying this new approach to high-resolution eye imagery of thousands of Europeans, we demonstrated that it mostly outperformed previous eye colour quantification methods. When only low-resolution imagery is available that does not allow for a distinct clustering of pixels into melanin types, average H and S values are recommended for usage.
The power to detect a causal variant in a genome wide association study is known to depend mainly on its effect and frequency in the population 42 . By using these detailed quantitative eye colour phenotypes, we noted that SNPs    previously involved in human eye colour variation showed stronger associations, revealed new and confirmed previously noted SNP-SNP interactions, and increased DNA-based prediction compared to quantitative eye colour phenotypes established by other methods. Overall, our findings imply that using our approach for detailed quantitative pigmentation phenotyping in future genome-wide association studies will likely deliver new pigmentation genes and new pigmentation predictive DNA variants, which is relevant for medical, evolutionary, and forensic genetics.

Materials andMethods
Ethics statement. The study was approved by national ethical committees and met the criteria of the Helsinki declaration. The ethics committees of each of the institutional review boards of the following collaborating eye study centres gave approval: Department of Epidemiology and Biostatistics, National Institute for Health Development,   Subjects, images, and genotyping. Eye images and DNA samples were collected as part of the EUREYE study. EUREYE is a population-based study of age related macular degeneration (AMD) in seven centers located across Europe. Participants were recruited from random sampling of the population aged over 65 years in Bergen (Norway), Tallinn (Estonia), Belfast (UK), Paris Creteil (France), Verona (Italy), Thessaloniki (Greece), and Alicante (Spain). Participants were interviewed by fieldworkers, underwent an eye examination including digital capture of the iris, and provided a blood sample for DNA analysis. A detailed description of the EUREYE study including image collection can be found elsewhere 38,45,46 . In brief, iris photography involved illumination of the anterior segment of each eye to show the colour of the iris with a flash intensity of 25-36 mW using a Topcon TRC 50 EX camera (http://www.topconmedical.com/categories/imaging-retinalcameras.htm) under normalized conditions. We manually excluded samples where the image quality was limited due to partially closed eyelids or largely dilated pupils that prevented a full view of the iris. DNA samples were genotyped using the SNaPShot technology (Life Technologies) via a single multiplex assay targeting 12-eye colour SNPs. We used 6 SNPs that were previously identified to predict categorical eye colour in genome wide studies 12,24,37 ) and which are included in the IrisPlex and HIrisPlex DNA test systems for categorical eye colour prediction 15,18,[37][38][39] : HERC2 rs12913832, OCA2 rs1800407, SLC24A4 rs12896399, SLC45A2 (MATP) rs16891982, TYR rs1393350, and IRF4 rs12203592. The 6 additional SNPs used were identified to be eye colour predictive in previous studies 24,40 : rs2070959 (UGT1A6), rs9894429 (NPLOC4), rs1129038 (HERC2), rs3768056 (LYST), rs2835630 (DSCR9), and rs1325127 (TYRP1) (see Supplementary Table S1 for primer sequences). Our analyses in the present manuscript were based on 3087 samples with suitable iris pictures and complete genotype data available.
Statistical analysis. Deviations from Hardy Weinberg equilibrium (HWE) were calculated using methods described elsewhere 47 . Among other causes, HWE can be violated because of genotyping error or population substructure. The latter is indicated by an increased informativeness of ancestry 48 . The statistical properties of the 12 SNPs from the 3087 samples are summarized in Supplementary Table S1. Three SNPs showed significant deviations from the Hardy-Weinberg equilibrium, namely HERC2 rs1129038 (P HWE = 8.25e-38), HERC2 rs12913832 (P HWE = 2.38e-37), and SLC45A2 rs16891982 (P HWE = 2e-16). In the absence of evidence for genotyping errors, we expect this being caused by the samples coming were from seven European populations as these three SNPs displayed elevated informativeness of ancestry (In) values (see Supplementary Table S1 for HWE and In values). Linkage disequilibrium (LD) between pairs of SNPs was estimated by means of R 2 . Regression analysis, Hellinger distance analysis, and interaction analysis were performed in Matlab (R2007a). For association and prediction analysis, linear models were fitted for both, individual SNPs and a combined model including all SNPs in Matlab (R2007a) while controlling for age, sex and population id. Presence of interactions between pairs of SNPs was tested by comparing nested linear models with and without the interaction term by means of an F-test. DNA-based eye colour prediction was achieved using general linear models Y~b 0 + b 1 X 1 + b 2 X 2 + b 3 X 3 + b 4 age + b 5 sex + b 6 pop, where X 1 and X 2 denote the genotypes of SNP 1 and SNP 2 respectively and X 3 denotes the interaction term on a multiplicative scale. Genotypes were coded additively as counts of minor alleles. P-value for the interaction term (b 3 ) were used and corrected for multiple testing using Bonferroni correction for 12 interaction pairs (396 tests applied). Linear models were fitted to predict the quantitative amounts of melanin from genotypic variation. To study the importance of genetic interaction, we added the most important interaction identified to the predictor using only main effects of SNPs. Quality of prediction was assessed by means of adjusted coefficient of determination (R 2 ) of the models. We used each colour quantification values separately as response variable. To account for over-fitting, we provide the mean R 2 from cross validation as estimated from 100 randomized train/holdout experiments where 2/3 of the samples were used as a training set (N = 2087) and 1/3 as a validation set (N = 1000).
The computational method to perform statistical analyses and extract quantitative eye colour from digital images is available on request.