The impact of correlations between pigmentation phenotypes and underlying genotypes on genetic prediction of pigmentation traits

Xiong Z.
Dankova G.
Howe L.J.
Lee M.K.
Hysi P.G.
de Jong M.A.
et al.

Novel genetic loci affecting facial shape variation in humans.

Google Scholar

], hair shape with up to 95 % [

[3]

Medland S.E.
Zhu G.
Martin N.G.

Estimating the heritability of hair curliness in twins of European ancestry.

]; hair and eye colour with up to 99 % [

[4]

Lin B.D.
Mbarek H.
Willemsen G.
Dolan C.V.
Fedko I.O.
Abdellaoui A.
et al.

Heritability and genome-wide association studies for hair color in a dutch twin family based sample.

] estimated heritability values. Various genome-wide association studies (GWASs), the more recent ones with large sample size, have revealed numerous underlying genes per each of several human appearance traits such as eye colour [

5

Sulem P.
Gudbjartsson D.F.
Stacey S.N.
Helgason A.
Rafnar T.
Magnusson K.P.
et al.

Genetic determinants of hair, eye and skin pigmentation in Europeans.

,

6

Kayser M.
Liu F.
Janssens A.C.
Rivadeneira F.
Lao O.
van Duijn K.
et al.

Three genome-wide association studies and a linkage analysis identify HERC2 as a human iris color gene.

,

7

Liu F.
Wollstein A.
Hysi P.G.
Ankra-Badu G.A.
Spector T.D.
Park D.
et al.

Digital quantification of human eye color highlights genetic association of three new loci.

], hair colour [

Sulem P.
Gudbjartsson D.F.
Stacey S.N.
Helgason A.
Rafnar T.
Magnusson K.P.
et al.

Genetic determinants of hair, eye and skin pigmentation in Europeans.

,

Han J.
Kraft P.
Nan H.
Guo Q.
Chen C.
Qureshi A.
et al.

A genome-wide association study identifies novel alleles associated with hair color and skin pigmentation.

,

Hysi P.G.
Valdes A.M.
Liu F.
Furlotte N.A.
Evans D.M.
Bataille V.
et al.

Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability.

]; eyebrow colour [

Peng F.
Zhu G.
Hysi P.G.
Eller R.J.
Chen Y.
Li Y.
et al.

Genome-wide association studies identify multiple genetic loci influencing eyebrow color variation in Europeans.

], skin colour [

Sulem P.
Gudbjartsson D.F.
Stacey S.N.
Helgason A.
Rafnar T.
Magnusson K.P.
et al.

Genetic determinants of hair, eye and skin pigmentation in Europeans.

,

Han J.
Kraft P.
Nan H.
Guo Q.
Chen C.
Qureshi A.
et al.

A genome-wide association study identifies novel alleles associated with hair color and skin pigmentation.

,

[11]

Visconti A.
Duffy D.L.
Liu F.
Zhu G.
Wu W.
Chen Y.
et al.

Genome-wide association study in 176,678 Europeans reveals genetic loci for tanning response to sun exposure.

], hair structure [

[12]

Medland S.E.
Nyholt D.R.
Painter J.N.
McEvoy B.P.
McRae A.F.
Zhu G.
et al.

Common variants in the trichohyalin gene are associated with straight hair in Europeans.

,

[13]

Liu F.
Chen Y.
Zhu G.
Hysi P.G.
Wu S.
Adhikari K.
et al.

Meta-analysis of genome-wide association studies identifies 8 novel loci involved in shape variation of human head hair.

], hair loss in men [

[14]

Hagenaars S.P.
Hill W.D.
Harris S.E.
Ritchie S.J.
Davies G.
Liewald D.C.
et al.

Genetic prediction of male pattern baldness.

], facial shape [

Xiong Z.
Dankova G.
Howe L.J.
Lee M.K.
Hysi P.G.
de Jong M.A.
et al.

Novel genetic loci affecting facial shape variation in humans.

Google Scholar

], body height [

15

Lango Allen H.
Estrada K.
Lettre G.
Berndt S.I.
Weedon M.N.
Rivadeneira F.
et al.

Hundreds of variants clustered in genomic loci and biological pathways affect human height.

,

16

Wood A.R.
Esko T.
Yang J.
Vedantam S.
Pers T.H.
Gustafsson S.
et al.

Defining the role of common variation in the genomic and biological architecture of adult human height.

,

17

Marouli E.
Graff M.
Medina-Gomez C.
Lo K.S.
Wood A.R.
Kjaer T.R.
et al.

Rare and low-frequency coding variants alter human adult height.

], explaining varying portions of their heritability and highlighting their complex genetic nature. Some appearance phenotypes are highly correlated with each other in human populations, such as certain facial shapes [

Xiong Z.
Dankova G.
Howe L.J.
Lee M.K.
Hysi P.G.
de Jong M.A.
et al.

Novel genetic loci affecting facial shape variation in humans.

Google Scholar

], and especially the different human pigmentation traits [

Peng F.
Zhu G.
Hysi P.G.
Eller R.J.
Chen Y.
Li Y.
et al.

Genome-wide association studies identify multiple genetic loci influencing eyebrow color variation in Europeans.

]. It is common knowledge that, for instance, blue eye colour shows a tendency to co-occur with blond hair and light skin colour in Europeans, red hair colour typically co-occurs with pale skin colour and often with freckles in Europeans, and dark skin colour usually co-occurs with black hair and brown eye colour in Africans, certain Asian populations, New Guineans and Australian Aborigines. Moreover, GWAS outcomes have revealed that the different human pigmentation traits share a large (but not complete) proportion of underlying genetic components that explains their phenotypic correlations [

Sulem P.
Gudbjartsson D.F.
Stacey S.N.
Helgason A.
Rafnar T.
Magnusson K.P.
et al.

Genetic determinants of hair, eye and skin pigmentation in Europeans.

,

Han J.
Kraft P.
Nan H.
Guo Q.
Chen C.
Qureshi A.
et al.

A genome-wide association study identifies novel alleles associated with hair color and skin pigmentation.

,

Peng F.
Zhu G.
Hysi P.G.
Eller R.J.
Chen Y.
Li Y.
et al.

Genome-wide association studies identify multiple genetic loci influencing eyebrow color variation in Europeans.

].

Previously, predictive DNA markers were identified and statistical prediction models were developed for all pigmentation-related traits, namely eye colour [

18

Liu F.
van Duijn K.
Vingerling J.R.
Hofman A.
Uitterlinden A.G.
Janssens A.C.
et al.

Eye color and the prediction of complex phenotypes from genotypes.

,

19

Walsh S.
Liu F.
Ballantyne K.N.
van Oven M.
Lao O.
Kayser M.

IrisPlex: a sensitive DNA tool for accurate prediction of blue and brown eye colour in the absence of ancestry information.

,

20

Ruiz Y.
Phillips C.
Gomez-Tato A.
Alvarez-Dios J.
Casares de Cal M.
Cruz R.
et al.

Further development of forensic eye color predictive tests.

], head hair colour [

Hysi P.G.
Valdes A.M.
Liu F.
Furlotte N.A.
Evans D.M.
Bataille V.
et al.

Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability.

,

21

Branicki W.
Liu F.
van Duijn K.
Draus-Barini J.
Pospiech E.
Walsh S.
et al.

Model-based prediction of human hair color using DNA variants.

,

22

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

,

23

Sochtig J.
Phillips C.
Maronas O.
Gomez-Tato A.
Cruz R.
Alvarez-Dios J.
et al.

Exploration of SNP variants affecting hair colour prediction in Europeans.

], skin colour [

24

Maronas O.
Phillips C.
Sochtig J.
Gomez-Tato A.
Cruz R.
Alvarez-Dios J.
et al.

Development of a forensic skin colour predictive test.

,

25

Walsh S.
Chaitanya L.
Breslin K.
Muralidharan C.
Bronikowska A.
Pospiech E.
et al.

Global skin colour prediction from DNA.

,

26

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

], eyebrow colour [

Peng F.
Zhu G.
Hysi P.G.
Eller R.J.
Chen Y.
Li Y.
et al.

Genome-wide association studies identify multiple genetic loci influencing eyebrow color variation in Europeans.

], and freckles [

[27]

Hernando B.
Ibanez M.V.
Deserio-Cuesta J.A.
Soria-Navarro R.
Vilar-Sastre I.
Martinez-Cadenas C.

Genetic determinants of freckle occurrence in the Spanish population: towards ephelides prediction from human DNA samples.

,

[28]

Kukla-Bartoszek M.
Pospiech E.
Wozniak A.
Boron M.
Karlowska-Pik J.
Teisseyre P.
et al.

DNA-based predictive models for the presence of freckles.

]. The recently established HIrisPlex-S system [

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

] currently represents the most complete DNA-based pigmentation prediction tool, allowing simultaneous prediction of eye, head hair, and skin colour from DNA, including low quality and low quantity forensic DNA, based on 41 carefully selected DNA markers and three separate prediction models. HIrisPlex-S reflects an extension of the previously developed IrisPlex system for eye colour [

Walsh S.
Liu F.
Ballantyne K.N.
van Oven M.
Lao O.
Kayser M.

IrisPlex: a sensitive DNA tool for accurate prediction of blue and brown eye colour in the absence of ancestry information.

] and the previous HIrisPlex system for eye and hair colour prediction [

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

]. Genotyping assays of the HIrisPlex-S system are available based on SNaPshot single base extension technology and capillary electrophoresis (CE) [

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

], in addition to two widely used massively parallel sequencing (MPS) technologies: Ion Torrent and Illumina [

[29]

Breslin K.
Wills B.
Ralf A.
Ventayol Garcia M.
Kukla-Bartoszek M.
Pospiech E.
et al.

HIrisPlex-S system for eye, hair, and skin color prediction from DNA: massively parallel sequencing solutions for two common forensically used platforms.

]. Moreover, the VISAGE Consortium recently incorporated the HIrisPlex-S DNA markers together with continental ancestry informative DNA markers to function as all-in-one tools for both MPS platforms separately (Xavier et al. under review; Palencia-Madrid et al. under review). All available genotyping assays of the IrisPlex, the HIrisPlex, and the HIrisPlex-S systems have been forensically validated [

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

,

29

Breslin K.
Wills B.
Ralf A.
Ventayol Garcia M.
Kukla-Bartoszek M.
Pospiech E.
et al.

HIrisPlex-S system for eye, hair, and skin color prediction from DNA: massively parallel sequencing solutions for two common forensically used platforms.

,

30

Walsh S.
Lindenbergh A.
Zuniga S.B.
Sijen T.
de Knijff P.
Kayser M.
et al.

Developmental validation of the IrisPlex system: determination of blue and brown iris colour for forensic intelligence.

,

31

Walsh S.
Chaitanya L.
Clarisse L.
Wirken L.
Draus-Barini J.
Kovatsi L.
et al.

Developmental validation of the HIrisPlex system: DNA-based eye and hair colour prediction for forensic and anthropological usage.

]. The three statistical prediction models, i.e., the IrisPlex model for eye colour [

Walsh S.
Liu F.
Ballantyne K.N.
van Oven M.
Lao O.
Kayser M.

IrisPlex: a sensitive DNA tool for accurate prediction of blue and brown eye colour in the absence of ancestry information.

], the HIrisPlex model for hair colour [

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

], and the HIrisPlex-S model for skin colour prediction [

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

] are all publically available in their most updated versions [

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

] via the website https://hirisplex.erasmusmc.nl/.

These laboratory and/or statistical tools are already of practical relevance for DNA-based pigmentation trait prediction in several different applications such as forensic investigation [

[32]

Kayser M.

Forensic DNA Phenotyping: predicting human appearance from crime scene material for investigative purposes.

,

[33]

Schneider P.M.
Prainsack B.
Kayser M.

The use of forensic DNA phenotyping in predicting appearance and biogeographic ancestry.

], human history inference [

[34]

King T.E.
Fortes G.G.
Balaresque P.
Thomas M.G.
Balding D.
Delser P.M.
et al.

Identification of the remains of king richard III.

], and anthropological research [

[35]

Brace S.
Diekmann Y.
Booth T.J.
van Dorp L.
Faltyskova Z.
Rohland N.
et al.

Ancient genomes indicate population replacement in early neolithic Britain.

], with more applications being expected in the future. Within the concept of Forensic DNA Phenotyping, predicting pigmentation traits of an unknown crime scene sample donor directly from crime scene DNA can provide useful investigative leads to find unknown perpetrators of the crime, in cases without a DNA profile match with a known suspect [

[32]

Kayser M.

Forensic DNA Phenotyping: predicting human appearance from crime scene material for investigative purposes.

,

[33]

Schneider P.M.
Prainsack B.
Kayser M.

The use of forensic DNA phenotyping in predicting appearance and biogeographic ancestry.

]. In human history investigations, DNA-based pigmentation prediction allows us to reveal the pigment of historical individuals from analysing their remains [

[34]

King T.E.
Fortes G.G.
Balaresque P.
Thomas M.G.
Balding D.
Delser P.M.
et al.

Identification of the remains of king richard III.

]. In anthropological and human evolutionary research, genetic pigmentation prediction captures how humans and human populations may have looked in the past, including the distant past [

[35]

Brace S.
Diekmann Y.
Booth T.J.
van Dorp L.
Faltyskova Z.
Rohland N.
et al.

Ancient genomes indicate population replacement in early neolithic Britain.

], and allows a deeper understanding of the evolutionary history of human pigmentation traits [

[36]

Wilde S.
Timpson A.
Kirsanow K.
Kaiser E.
Kayser M.
Unterlander M.
et al.

Direct evidence for positive selection of skin, hair, and eye pigmentation in Europeans during the last 5,000 y.

].

However, all currently available genetic prediction models for pigmentation traits ignore the well-known correlations between the different pigmentation phenotypes and their underlying genotypes. It could be expected, however, that considering phenotype correlations in the genetic prediction modelling may increase the accuracy of DNA-based prediction for eye, hair and skin colour, which to the best of our knowledge has not been reported as of yet. Here, we empirically tested for the impact of correlations between appearance phenotypes and their underlying genotypes on DNA-based appearance prediction using pigmentation traits as classical example. We applied categorical eye, hair and skin colour phenotype data and genotype data for the 41 HIrisPlex-S DNA markers from 762 individuals for whom complete phenotype and genotype information was available. Based on these data, we empirically estimated pigmentation phenotype correlations and their proportions that were attributable to the 41 HIrisPlex-S DNA markers. Next, we performed genetic prediction modelling of eye, hair and skin colour via three different strategies, namely the established approach of predicting phenotypes solely based on genotypes while not considering phenotype correlations, and two novel approaches considering phenotype correlations, either incorporating truly observed correlated phenotypes or DNA-predicted correlated phenotypes as additional predictors, and compared the prediction accuracies of these different models by empirical observation. Finally, we conducted computer simulations, emulating the three different prediction strategies, to better understand the impact of phenotype and genotype correlations on DNA-based phenotype prediction in an effort to interpret the outcomes of the empirical pigmentation prediction modelling we obtained.

2. Materials and methods

2.1 Phenotype and genotype data

The data used here were all from the previous IrisPlex, HIrisPlex and HIrisPlex-S projects (https://hirisplex.erasmusmc.nl/) as described elsewhere [

Walsh S.
Liu F.
Ballantyne K.N.
van Oven M.
Lao O.
Kayser M.

IrisPlex: a sensitive DNA tool for accurate prediction of blue and brown eye colour in the absence of ancestry information.

,

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

,

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

] and represent a subset of 762 individual datasets of different populations from Europe and the US for which a complete pigmentation phenotype profile (i.e., categorical eye, hair, and skin colour) and a complete genotype profile (i.e., all 41 HIrisPlex-S DNA markers) were jointly available (Table 1). Samples had been collected for the purpose of appearance genetic research under written informed consent, and sample collections where approved by the Ethics Committee of the Jagiellonian University (KBET/17/B/2005), the Commission on Bioethics of the Regional Board of Medical Doctors in Krakow (48 KBL/OIL/2008), and by the Indiana University Ethical Institutional Review Board (#1409306349). As previously described in detail [

Walsh S.
Liu F.
Ballantyne K.N.
van Oven M.
Lao O.
Kayser M.

IrisPlex: a sensitive DNA tool for accurate prediction of blue and brown eye colour in the absence of ancestry information.

,

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

,

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

], eye colour was classified into three categories: blue, intermediate, and brown; hair colour into four categories: red, blond, brown, and black; and skin colour into five categories: very pale, pale, intermediate, dark, and dark to black. The 41 HIrisPlex-S DNA markers were described elsewhere [

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

,

Chaitanya L.
Breslin K.
Zuniga S.
Wirken L.
Pospiech E.
Kukla-Bartoszek M.
et al.

The HIrisPlex-S system for eye, hair and skin colour prediction from DNA: introduction and forensic developmental validation.

].

Table 1Characteristics of the study population representing a subset of the HIrisPlex-S dataset with complete data on eye, hair, and skin colour phenotypes and 41 HIrisPlex-S SNP genotypes.

		N	Proportion of trait categories (%)	Males/Females
		762		328/434
Eye colour	Blue	328	43.04	158/170
	Intermediate	97	12.73	37/60
	Brown	337	44.23	133/204
Hair colour	Red	17	2.23	7/10
	Blond	330	43.31	145/185
	Brown	328	43.04	124/204
	Blond	87	11.42	52/35
Skin colour	Very Pale	28	3.67	13/15
	Pale	355	46.59	142/213
	Intermediate	315	41.34	145/170
	Dark	36	4.72	19/17
	Dark to Black	28	3.67	9/19

Open table in a new tab

2.2 Statistical analyses

Since we considered here, for subsequent genetic prediction modelling, a subset of the original data sets previously applied for the initial IrisPlex eye colour, HIrisPlex hair colour and HIrisPlex-S skin colour models, we first evaluated this reduced data subset with respect to sufficient statistical power and concordant SNP effects. To this end, ordinal phenotypes were considered as continuous variables by assigning ascending integer values, i.e. 1, 2 and so forth (from the lightest to the darkest levels), to the ordered trait categories and Z-transformed (

Z = \frac{y - \bar{y}}{s d (y)}, y

is the trait,

\bar{y}

is the mean value of

y

,

s d (y)

is the standard deviation of

y

). Then we tested each SNP for a phenotypic association using a linear regression model adjusted for sex, thereby taking advantage of the robustness of the linear regression approach in the absence of a normally distributed dependent variable. In addition, multivariate analysis of variance (MANOVA) was conducted to test the association between a given SNP and the three pigmentation phenotypes simultaneously. For association testing, two missense MC1R DNA variants rs1805007_T (R151C) and rs1805008_T (R160W), which are known to be involved in human red hair and related pigmentation phenotypes in a compound heterozygote manner [

[37]

Liu F.
Struchalin M.V.
Duijn K.
Hofman A.
Uitterlinden A.G.
Duijn C.
et al.

Detecting low frequent loss-of-function alleles in genome wide association studies with red hair color as example.

] were collapsed into three possible genotypes wt/wt/, wt/R, and R/R, where R is the risk haplotype consisting of at least one minor allele from any of two MC1R variants and wt is the wild-type haplotype consisting of only major alleles. Then, this MC1R compound marker was treated as a discrete bi-allelic marker in the association analysis, while the remaining MC1R SNPs were used as discrete markers.

The 41 HIrisPlex-S DNA markers are located at 11 different genetic loci. For each locus, the associated SNP with the smallest P-value (i.e. the strongest association effect) in a MANOVA was selected and its allelic effects across different phenotypes were investigated. In this analysis, ordinal phenotypes were considered as continuous variables and Z-transformed as introduced before. The effect alleles were selected to have a colour darkening effect and the allelic effect was estimated as E=

2 f (1 - f) β^{2}

, where

f

is the frequency of the effect allele and

β

is the regression coefficient from linear models. The E ranges from 0 for no effect to 1 for fully explaining the phenotype variance and can be interpreted as the genetically explained proportion of the phenotype variance. Ordinal phenotypes (from the lightest to the darkest levels) were considered as continuous variables and Pearson’s correlation coefficient (r) was calculated for each pair of phenotypes.

To quantify the extent to which the correlation between two phenotypes can be explained by their shared explanatory factors under investigation, we derived a statistic

C = \frac{c o r (y_{1}, {\hat{y}}_{1}) c o r ({\hat{y}}_{1}, {\hat{y}}_{2}) c o r ({\hat{y}}_{2}, y_{2})}{c o r (y_{1}, y_{2})}

where

y_{1}

and

y_{2}

represent any two of the three correlated pigmentation phenotypes, and

{\hat{y}}_{1}

and

{\hat{y}}_{2}

represent the predicted values from respective linear models.

C

typically ranges between 0 and 1, although values outside this interval are possible under extreme scenarios. A value of 0 indicates that the considered explanatory factors cannot explain any of the observed phenotypic correlation while 1 represents the case that the observed correlation can be perfectly explained by the considered factors. For example, with values of 0.5 for

C

and 0.8 for the observed phenotypic correlation, 50 % of this correlation can be explained by the considered set of predictors, whereas the remaining 50 % are due to other shared genetic factors, such as other SNPs or other forms of genetic variation, and non-genetic factors, such as age, sex or environmental factors. We then applied this approach to the 41 HIrisPlex-S SNPs constituting the set of considered predictors to assess the proportion of phenotypic correlation that can be explained by these genetic markers.

Moreover, we empirically assessed the statistical properties of

C

by computer simulations in order to ensure that

C

yields an unbiased prediction of the explained correlation. To this end, we decomposed the set of all contributing factors into three mutually exclusive sets (i.e. a variance decomposition by orthogonal factors), namely genetic factors that are shared between phenotypes and known, accessible and therefore included in the model, s₁, genetic and environmental factors that are shared between phenotypes but are unknown, s₂, and genetic and environmental factors that are unique to a particular phenotype, u. We therefore assume that accessible shared environment has a negligible effect in our model. More specific, we generated two correlated traits,

y_{1}

and

y_{2}

, for 1000 individuals by simulating their shared and unique explanatory components with identical effect sizes.

y_{1} = u_{1} + {α s}_{1} + s_{2}

y_{2} = u_{2} + {α s}_{1} + s_{2}

u_{1}, u_{2} i n d e p e n d e n t a n d i d e n t i c a l l y d i s t r i b u t e d (i . i . d .) \sim N (0, \frac{1 + α^{2}}{9});

s_{1}, s_{2} i . i . d . \sim N (0,1);

α = 0.00, 0.33, 0.50, 0.65, 0.82, 1.00, 1.22, 1.53, 2.00, 3.00

where

u_{1}

and

u_{2}

represent explanatory factors having effects unique to

y_{1}

and

y_{2}

, respectively, while

s_{1}

represents accessible genetic factors and

s_{2}

unknown factors, whereas

α

is used to regulate the shared component variance proportion explained by

s_{1}

. For pigmentation traits, u₁ and u₂ may represent genetic and environmental factors having an effect on one pigmentation trait but not on another and vice versa, s₁ may combine all 41 SNPs and s₂ may represent yet to be discovered genetic factors influencing both traits. The expectation of

C

is

\frac{α^{2}}{1 + α^{2}}

(see detailed derivation in Supplementary Materials). Then the estimated

C

was investigated under a range of expectation of

C, E (C) = 0 %, 10 %, 20 %, 30 %, 40 %, 50 %, 60 %, 70 %, 80 %, 90 %

(

c o r r e s p o n d i n g α = 0.00, 0.33, 0.50, 0.65, 0.82, 1.00, 1.22, 1.53, 2.00, 3.00

). The simulation was conducted with 1000 replicates for each

E (C)

.

To evaluate the impact of incorporating additional correlated phenotypes as predictors in the model on prediction accuracy for the targeted phenotype, we compared three strategies by including different sets of predictors. In strategy I, we accessed the base-line prediction accuracies when solely using the 41 HIrisPlex-S DNA markers as the predictors. In strategy II, in addition to the 41 HIrisPlex-S DNA markers, the fitted values of the correlated phenotypes that were predicted from the 41 DNA markers were used as additional predictors. In strategy III, in addition to the 41 HIrisPlex-S DNA markers, the truly observed correlated phenotypes were used as additional predictors. In particular for strategy II and III, for predicting eye colour, DNA-predicted (II) or truly observed (III) hair colour and skin colour phenotypes were included as additional predictors; for predicting hair colour, DNA-predicted (II) or truly observed (III) eye and skin colour phenotypes were included as additional predictors; and when predicting skin colour, DNA-predicted (II) or truly observed (III) eye and hair colour phenotypes were included as additional predictors in the genetic prediction modelling.

All predictions were made via multinomial logistic regression models, while using standard leave-one-out (LOO) cross-validation (CV). However, in contrast to previous IrisPlex, HIrisPlex, and HIrisPlex-S studies that used 6, 22, and 36 of the 41 HIrisPlex-S SNPs, respectively, for eye, hair and skin colour prediction, respectively, we considered all 41 HIrisPlex-S SNPs for the prediction of all three pigmentation traits in the current study in order to test for the complete effect of correlated phenotypes and genotypes. Therefore, the prediction outcomes obtained here for the standard approach not considering correlated phenotypes (strategy I) are not directly comparable with previously reported eye, hair, and skin colour prediction outcomes based on 6 IrisPlex, 22 HIrisPlex, and 36 HIrisPlex-S, respectively. The prediction accuracies were derived using the Area Under the receiver operating characteristic Curves (AUC) as well as other commonly used prediction statistics including sensitivity, specificity, negative predictive value (NPV), and positive predictive value (PPV). All statistical analyses and result visualization were conducted in R version 3.5.3 [https://www.R-project.org/] using following packages: stats of version 3.5.3, nnet of version 7.3-12, pROC of version 1.15.0.

2.3 Computer simulations

To better understand and interpret the results from the empirical prediction modelling described above, we conducted computer simulations as specified in the following. We generated two correlated traits

y_{1}

and

y_{2}

for 1000 individuals by simulating their shared and unique explanatory components with the identical effect sizes,

y_{1} = u_{11} + u_{12} + s_{1} + {α s}_{2}

y_{2} = u_{21} + u_{22} + s_{1} + {α s}_{2}

u_{11}, u_{12}, u_{21}, u_{22}, s_{1}, s_{2} i . i . d . \sim N (0,1)

where the

u_{11}

,

u_{12}

represents explanatory factors having effects unique for

y_{1}

, and the

u_{21}

and

u_{22}

are explanatory factors unique for

y_{2}

, while only

u_{11}

and

u_{21}

are accessible but

u_{12}

and

u_{22}

are inaccessible by investigators. Similarly,

s_{1}

and

s_{2}

are shared explanatory factors having effects on both

y_{1}

and

y_{2}

while only

s_{1}

is accessible by investigators but

s_{2}

is not.

α

is used to regulate the phenotype variance proportion (

V = \frac{α^{2}}{3 + α^{2}}

) explained by the inaccessible factor of

s_{2}

. The variance explained (R [

Xiong Z.
Dankova G.
Howe L.J.
Lee M.K.
Hysi P.G.
de Jong M.A.
et al.

Novel genetic loci affecting facial shape variation in humans.

Google Scholar

]) was derived from linear models fitted using 3 different sets of predictors as specified below,

M o d e l I : y_{1} \sim u_{11} + s_{1}

M o d e l I I : y_{1} \sim u_{11} + s_{1} + {\hat{y}}_{2}

M o d e l I I I : y_{1} \sim u_{11} + s_{1} + y_{2}

where

{\hat{y}}_{2}

is the fitted

y_{2}

using

u_{21} a n d s_{1}

. Thus, Model I mimics a typical genotype-phenotype prediction analysis without considering phenotype and genotype correlations, analogue to strategy I in our empirical prediction analysis. Model II mimics the scenario when additional correlated phenotypes are not available, but were predicted from the same set of pre-selected SNPs and used as additional predictors, analogue to strategy II in our empirical prediction analysis. Model III mimics the scenario when truly observed correlated phenotypes are available for prediction and used as additional predictors, analogue to strategy III in our empirical prediction analysis. All simulations were conducted under

V = 0 %, 10 %, 30 %, 50 %, 70 % (α = 0.00, 0.58, 1.13 1.73, 2.65)

. In addition, we repeated the simulation process by dichotomizing y using the mean value and estimated the AUC values using logistic models, which mimics scenarios of binary trait analysis. All simulations were conducted for 1000 replicates under each investigated model/scenario.

3. Results

3.1 Data suitability check via genetic association testing

Due to the need for complete eye, hair, and skin data and 41 HIrisPlex-SNP data in this study, we used data from 762 HIrisPlex-S subjects for whom such complete phenotype and genotype data was jointly available (Table 1). This dataset represents a subset of the data previously applied to develop the IrisPlex eye colour, HIrisPlex hair colour and HIrisPlex-S skin colour prediction models. Therefore, we first assessed the suitability of this data subset for genetic prediction modelling by testing if the previously reported genetic associations and SNP effects can be replicated in this particular dataset. To this end, we conducted genetic association testing using MANOVA and linear regression for all 41 HIrisPlex-S SNPs regarding eye, hair and skin colour. While with MANOVA all three pigmentation traits were considered in a combined way, with linear regression analysis each pigmentation trait was tested separately. As may be expected, there was a high correlation between the genetic association outcomes from both approaches (Table 2).

Table 2Association between the 41 HIrisPlex-S DNA markers and eye, hair, and skin colour in the study dataset (N = 762).

								MANOVA	Linear Regression
								Eye/hair/skin colour	Eye Colour	Hair Colour	Skin Colour
Nr.	SNP	EA	OA	Freq	MBp	Locus	Gene	P-value	P-value	P-value	P-value
1	rs16891982	C	G	0.12	33.95	5p13.2	SLC45A2	*1.44E-87*	*5.89E-23*	*5.03E-32*	*5.50E-81*
2	rs28777	C	A	0.1	33.96	5p13.2	SLC45A2	*5.90E-69*	*1.31E-19*	*2.02E-25*	*5.03E-65*
3	rs12203592	T	C	0.1	0.4	6p25.3	IRF4	8.48E-08	7.55E-02	1.23E-01	1.32E-05
4	rs4959270	A	C	0.44	0.46	6p25.3	EXOC2	1.16E-02	2.89E-02	1.09E-01	1.93E-03
5	rs683	G	T	0.39	12.71	9p23	TYRP1	*1.53E-13*	2.04E-04	*2.50E-08*	*2.15E-13*
6	rs10756819	G	A	0.4	16.86	9p22.2	BNC2	1.08E-02	8.80E-01	2.64E-01	1.51E-03
7	rs1042602	T	G	0.33	88.91	11q14.3	TYR	6.27E-06	2.48E-03	7.57E-04	1.12E-06
8	rs1393350	T	C	0.21	89.01	11q14.3	TYR	3.57E-02	2.01E-01	3.31E-02	5.66E-03
9	rs1126809	A	G	0.19	89.02	11q14.3	TYR	9.20E-03	1.39E-01	1.59E-03	1.18E-02
10	rs12821256	G	A	0.08	89.33	12q21.33	KITLG	1.54E-03	3.82E-04	4.99E-03	7.44E-03
11	rs12896399	T	G	0.41	92.77	14q32.12	SLC24A4	6.42E-04	7.26E-03	3.63E-02	9.52E-05
12	rs2402130	G	A	0.19	92.8	14q32.12	SLC24A4	5.90E-03	8.20E-03	7.82E-03	3.52E-03
13	rs17128291	C	T	0.19	92.88	14q32.12	SLC24A4	5.06E-02	3.06E-01	8.70E-02	6.99E-03
14	rs1545397	T	A	0.11	28.19	15q13.1	OCA2	3.81E-04	2.98E-04	3.60E-04	6.73E-03
15	rs1800414	C	T	0.01	28.2	15q13.1	OCA2	6.73E-06	3.99E-04	2.82E-06	3.74E-04
16	rs1800407	A	G	0.07	28.23	15q13.1	OCA2	3.65E-04	9.06E-03	4.05E-01	2.36E-02
17	rs12441727	A	G	0.17	28.27	15q13.1	OCA2	1.85E-03	4.60E-04	2.61E-02	5.75E-03
18	rs1470608	A	C	0.19	28.29	15q13.1	OCA2	*2.28E-23*	*1.09E-11*	*1.48E-08*	*2.92E-22*
19	rs1129038	G	A	0.29	28.36	15q13.1	HERC2	*7.02E-46*	*9.06E-32*	*8.29E-19*	*3.36E-31*
20	rs12913832	T	C	0.29	28.37	15q13.1	HERC2	*1.81E-166*	*6.24E-144*	*1.11E-47*	*1.57E-47*
21	rs2238289	C	T	0.17	28.45	15q13.1	HERC2	*6.75E-24*	*3.38E-14*	*7.92E-12*	*3.17E-20*
22	rs6497292	C	T	0.09	28.5	15q13.1	HERC2	*1.81E-20*	*6.56E-10*	*1.18E-09*	*1.62E-19*
23	rs1667394	C	T	0.23	28.53	15q13.1	HERC2	*5.19E-35*	*9.11E-19*	*4.03E-14*	*3.84E-30*
24	rs1426654	G	A	0.07	48.43	15q21.1	SLC24A5	*1.82E-73*	*2.73E-17*	*2.01E-24*	*1.69E-71*
25	rs3114908	T	C	0.32	89.38	16q24.3	ANKRD11	1.36E-01	5.00E-01	1.07E-01	7.25E-01
26	rs3212355	A	G	0	89.98	16q24.3	MC1R	5.42E-01	2.38E-01	7.44E-01	3.73E-01
27	rs312262906	A	---	0	89.99	16q24.3	MC1R	4.27E-02	1.61E-01	2.25E-01	4.41E-01
28	rs1805005	T	G	0.09	89.99	16q24.3	MC1R	1.44E-01	7.12E-01	5.55E-01	2.20E-02
29	rs1805006	A	C	0	89.99	16q24.3	MC1R	6.62E-01	6.48E-01	7.94E-01	3.88E-01
30	rs2228479	A	G	0.1	89.99	16q24.3	MC1R	2.29E-02	4.92E-01	3.09E-01	2.16E-03
31	rs11547464	A	G	0.01	89.99	16q24.3	MC1R	3.61E-01	1.86E-01	1.33E-01	8.87E-01
32	rs1805007	T	C	0.06	89.99	16q24.3	MC1R	4.77E-06	3.74E-01	2.02E-04	1.89E-05
33	rs201326893	A	C	0	89.99	16q24.3	MC1R	5.47E-03	2.37E-02	2.08E-01	5.62E-01
34	rs1110400	C	T	0.01	89.99	16q24.3	MC1R	1.45E-01	3.32E-02	1.94E-01	1.08E-01
35	rs1805008	T	C	0.06	89.99	16q24.3	MC1R	5.59E-05	4.23E-01	1.35E-03	5.07E-05
	Compound	---	---	0.13	---	16q24.3	MC1R	*1.77E-12*	2.00E-01	1.35E-07	*2.03E-10*
36	rs885479	T	C	0.05	89.99	16q24.3	MC1R	5.10E-05	6.96E-03	1.03E-05	1.33E-03
37	rs1805009	C	G	0.01	89.99	16q24.3	MC1R	7.76E-06	8.28E-01	3.11E-04	5.48E-04
38	rs8051733	C	T	0.31	90.02	16q24.3	DEF8	1.82E-01	2.51E-01	8.15E-02	6.78E-02
39	rs6059655	T	C	0.05	32.67	20q11.22	RALY	3.01E-04	8.15E-01	3.95E-03	3.38E-04
40	rs6119471	C	G	0.02	32.79	20q11.22	ASIP	*8.59E-31*	1.26E-06	4.58E-05	*8.20E-33*
41	rs2378249	C	T	0.16	33.22	20q11.22	ASIP/PIGU	3.52E-01	2.96E-01	5.78E-01	4.48E-01

Statistical significance threshold: nominally significant p < 0.05 in italic, significant after Bonferroni correction of multiple testing p < 0.0012 in bold, and genome-wide significant p<5 × 10⁻⁸ in bold and italic.

Open table in a new tab

With MANOVA and by combining the three pigmentation traits, 32 (78 %) of the 41 HIrisPlex-S DNA markers from all the 11 different genetic loci showed nominally significant association (p < 0.05), 22 (54 %) DNA markers from 9 loci showed significant association after Bonferroni correction for multiple testing (p < 0.0012), and 11 (27 %) SNPs from 5 loci even showed significant association on the genome-wide level (p<5 × 10⁻⁸) representing an over-conservative significance threshold in a candidate SNP approach as applied here (Table 2). On chr16, the MC1R compound marker showed stronger association (p = 1.77e-12) than both markers separately i.e., rs1805007 (p = 4.77e-6) and rs1805008 (p = 5.59e-5). For eye colour and hair colour, the most significant genetic association was seen for HERC2 rs12913832 (eye colour p = 6.24 × 10⁻¹⁴⁴, and hair colour p = 1.11 × 10⁻⁴⁷), while for skin colour the most significant association was observed for SLC45A2 rs16891982 (p = 5.5 × 10⁻⁸¹).

With linear regression analysis and by treating the three pigmentation phenotypes separately, we found 19 (46 %) DNA markers from 9 loci significantly associated with eye, hair and skin colour on the nominal significance level, 12 (29 %) DNA markers from 5 loci significantly associated on the Bonferroni significance level, and 9 (22 %) from 3 loci significantly associated on the over-conservative genome-wide significance level. In addition, 34 (83 %) DNA markers from all 11 loci showed nominal significant association for at least one of the three pigmentation phenotypes, 23 (56 %) from 10 loci based on the Bonferroni level, and 11 (27 %) from 5 loci on the over-conservative genome-wide level. Only 7 (17 %) of the 41 HIrisPlex-S DNA markers showed no association with any of the three pigmentation traits in this sample set. These 7 non-significant DNA markers were from two of the 11 genetic loci, MC1R (6 markers) and ASIP/RALY (one marker). However, both of these genetic loci showed significant associations with other SNPs in the HIrisPlex-S marker set.

Next, for each of the 11 genetic loci covered by the 41 HIrisPlex-S DNA markers, the top-associated SNP from MANOVA was investigated with respect to the contribution of its alleles in explaining the observed variance of the three different pigmentation traits (Fig. 1 and Supplementary Table 1). HERC2 s12913832_T explained an extraordinarily large proportion (48 %) of eye colour variance. Except HERC2 rs12913832, the majority of the SNPs tested showed their largest effect on skin colour (1.4 %–25.6 % explained variance). The compound marker in MC1R explained 4 %–6 % of the phenotypic variance of hair and skin colour but had little effect on eye colour. The relatively small effect (4% explained variance) of the MC1R compound marker on hair colour may be caused by the low frequency of red hair (2.2 %) in our dataset and that the remaining causal variants of MC1R were not considered in this compound marker. All of the allelic effects were on the same direction of darkening, except for IRF4 with a non-significant effect on lightening of hair colour, indicating the homogeneous effect of genetic variants on pigmentation phenotypes.

Fig. 1Effects of the 11 top-associated HIrisPlex-S SNPs from the 11 genetic loci covered by the 41 HIrisPlex-S DNA markers on eye, hair, and skin colour (N = 762). Compound represents a collapsed compound heterozygosity marker based on a haplotype analysis of two pre-selected *MC1R* coding DNA variants rs1805007_T (R151C) and rs1805008_T (R160W). Note that this compound *MC1R* marker has a colour lightening effect on pigmentation phenotypes, but this effect was reversely depicted in the figure for convenience reasons. For underlying data see **Supplementary** Table 1 (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

Given that statistically significant associations and genetic effects we observed with this dataset are in broad agreement with findings from previous studies based on larger sample size [

7

Liu F.
Wollstein A.
Hysi P.G.
Ankra-Badu G.A.
Spector T.D.
Park D.
et al.

Digital quantification of human eye color highlights genetic association of three new loci.

,

8

Han J.
Kraft P.
Nan H.
Guo Q.
Chen C.
Qureshi A.
et al.

A genome-wide association study identifies novel alleles associated with hair color and skin pigmentation.

,

9

Hysi P.G.
Valdes A.M.
Liu F.
Furlotte N.A.
Evans D.M.
Bataille V.
et al.

Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability.

,

[11]

Visconti A.
Duffy D.L.
Liu F.
Zhu G.
Wu W.
Chen Y.
et al.

Genome-wide association study in 176,678 Europeans reveals genetic loci for tanning response to sun exposure.

,

[38]

Jacobs L.C.
Wollstein A.
Lao O.
Hofman A.
Klaver C.C.
Uitterlinden A.G.
et al.

Comprehensive candidate gene study highlights UGT1A and BNC2 as new genes determining continuous skin color variation in Europeans.

,

[39]

Liu F.
Visser M.
Duffy D.L.
Hysi P.G.
Jacobs L.C.
Lao O.
et al.

Genetics of skin color variation in Europeans: genome-wide association studies with functional follow-up.

], we deemed this dataset useful for phenotype-genotype correlation analyses and for genetic prediction modelling.

3.2 Phenotype correlations and genetic contributions

A high correlation between eye, hair and skin colour is generally expected in human populations. In our dataset, all three pigmentation traits showed mid-range but highly statistically significant phenotype correlations with each other, namely Pearson correlation coefficients of 0.47 (p = 5.19 × 10⁻⁴³) for eye versus hair colour, 0.41 (p = 1.64 × 10⁻³²) for hair versus skin colour, and 0.36 (p = 9.42 × 10⁻²⁵) for eye versus skin colour (Table 3). That we did not obtain higher phenotype correlations may be explained by sampling errors due to the medium sample size and due to phenotype categorization from underlying quantitative variables in combination with imperfect phenotype quality based on the phenotyping methods applied.

Table 3Quantification of pigmentation phenotype correlations and their proportions contributed by the 41 HIrisPlex-S DNA markers.

	Eye colour	Hair colour	Skin colour
Eye Colour	*	0.76	0.96
Hair Colour	0.47 (5.19 × 10E-43)	*	0.87
Skin Colour	0.36 (9.42E-25)	0.41 (1.64E-32)	*

Below the diagonal line: phenotype correlations estimated by Pearson correlation coefficient and their statistical significance. Above the diagonal line: contribution of the 41 HIrisPlex-S SNPs on the observed pigmentation phenotype correlations (C statistics).

Open table in a new tab

Moreover, given the shared genetic components between the three pigmentation phenotypes seen for several of the 41 HIrisPlex DNA markers in earlier studies [

Sulem P.
Gudbjartsson D.F.
Stacey S.N.
Helgason A.
Rafnar T.
Magnusson K.P.
et al.

Genetic determinants of hair, eye and skin pigmentation in Europeans.

,

Han J.
Kraft P.
Nan H.
Guo Q.
Chen C.
Qureshi A.
et al.

A genome-wide association study identifies novel alleles associated with hair color and skin pigmentation.

], as well as in the present study (Fig. 1), a large contribution of these genetic markers on the phenotype correlations would be expected. To investigate this empirically, we derived a statistical estimate C (see method section for details) to assess the contribution of the 41 HIrisPlex-S DNA markers on the observed pigmentation phenotype correlations. With this analysis, we estimated the proportion of the shared genetic components explained by the 41 HIrisPlex-S DNA markers over all underlying shared genetic and non-genetic components, i.e., i) from unused shared DNA markers in LD with the used markers from the same genetic loci, ii) from unused shared DNA markers from other genetic loci not used here, and iii) from underlying shared non-genetic components such as environmental effects, which may all contribute to the observed phenotype correlations. The statistical property of C was examined via computer simulations, where the estimated C was in the expected range between 0 and 1. It strikingly corresponded to the proportion of the phenotype variance explained by the shared and accessible genetic components over the variance explained by all shared genetic and non-genetic components (Supplementary Figure S1). Based on our empirical data, the 41 HIrisPlex-S SNPs explained 96 % of the phenotype correlation between eye and skin colour, 87 % of the phenotype correlation between hair and skin colour, and 76 % of the phenotype correlation between eye and hair colour (Table 3). This finding demonstrates that a large proportion of the pairwise pigmentation phenotype correlations is explained by the HIrisPlex-S DNA markers used, but there also is a remaining proportion that remains unexplained by these genetic markers that differs between the pairwise phenotype comparisons.

3.3 Empirical impact of phenotype correlations on genetic phenotype prediction

To investigate the impact of the correlations between the pigmentation phenotypes and the underlying genotypes on DNA-based prediction of eye, hair, and skin colour, we performed genetic prediction modelling for the three pigmentation traits based on three different strategies (for details see methods). As evident from Fig. 2 (Supplementary Table 2), the standard genetic prediction model that does not consider pigmentation phenotype correlations (strategy 1), and the novel model incorporating DNA-predicted correlated pigmentation phenotypes as additional predictors (strategy 2), performed almost identical or very similar across all eye and hair colour categories as well as for most skin colour categories. For skin colour (almost) the same AUC values were achieved with both prediction strategies for pale, intermediate and dark, while for very pale the strategy 2 model achieved a slightly lower (0.04) and for dark to black a slightly higher (0.02) AUC. Similar findings were obtained for other prediction accuracy estimates (Supplementary Table 2). Hence, using correlated pigmentation phenotypes that were predicted from a relatively small set of selected DNA markers as additional predictors had no or almost no impact on DNA-based prediction accuracy of eye, hair, and skin colour.

Fig. 2Accuracies of pigmentation phenotype prediction achieved by considering the 41 HIrisplex-S SNPs (N = 762) based on three strategies Strategy I (depicted in red) included the 41 HIrisplex-S SNPs as sole predictors; strategy II (depicted in yellow) included the 41 HIrisPlex-S SNPs as predictors and as additional predictors the two respective other correlated pigmentation phenotypes *predicted* from the 41 DNA markers; strategy III (depicted in blue) included the 41 HIrisplex-S SNPs as predictors and as additional predictors the two respective other correlated *truly observed* pigmentation phenotypes. For underlying data see **Supplementary** Table 2 (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

This was markedly different for the strategy 3 model incorporating truly observed correlated pigmentation phenotypes in addition to the HIrisPlex-S SNPs in the genetic prediction modelling. Compared to the strategy 1 model that only used HIrisPlex-S SNPs, with the strategy 3 model we observed a slight increase in AUC for blue eye colour (0.016), blond hair colour (0.014) and black hair colour (0.016), while a slight decrease in AUC for very pale skin colour (0.02), and a noticeable increase in AUC for intermediate eye colour (0.037) and brown hair colour (0.030). However, a larger increase was seen for dark to black skin colour (0.056), and the AUC increase was largest for dark skin colour (0.176). For the remaining eye, hair and skin colour categories almost the same AUCs were obtained with both prediction strategies. Similar findings were obtained for other prediction accuracy estimates (Supplementary Table 2). The largest effects for dark and dark to black skin colour cannot necessarily be explained by low samples size, since very pale skin colour and red hair colour had similarly low sample size but showed much smaller effects (Table 1).

3.4 Simulated impact of phenotype correlations on genetic phenotype prediction

Next, we performed dedicated computer simulations to better understand and interpret these empirical findings using three models and considering both continuous and categorical phenotype information (for model details see method). These simulations showed for both continuous and categorical phenotypes, that model I without considering phenotype correlations (comparable with the empirical strategy 1 model) and model II with considering DNA-predicted correlated phenotypes (comparable with the empirical strategy 2 model) performed identical (Fig. 3), which agrees with our empirical findings. In contrast, model III with considering truly observed correlated phenotypes (comparable with the empirical strategy 3 model), achieved higher prediction accuracies compared to model I/II, which is in line with our empirical results. This finding suggests that additional information of unobserved shared component (

s_{2}

) was included in model III using truly observed correlated phenotypes, but not in model II using DNA-predicted correlated phenotypes. Moreover, as the variance explained by

s_{2}

increases, more additional information was included and therefore increased improvement in prediction accuracy is seen with model III that uses truly observed correlated phenotypes, which was not seen with model II that uses DNA-predicted correlated phenotypes (Fig. 3).

Fig. 3Results of a computer simulation investigating accuracy of genetic phenotype prediction without and with correlated phenotypes used as additional predictors. Model I (depicted in red) mimics the scenario of typical phenotype prediction from genotypes without considering correlated phenotypes as additional predictors, which corresponds to the strategy 1 empirical model (see red bars in Fig. 2); model II (depicted in yellow) mimics the scenario of correlated phenotypes being used as additional predictors, which, however, were *predicted* using the same set of pre-selected DNA markers that were used for genetic prediction, which corresponds to the strategy 2 empirical model (see yellow bars in Fig. 2); model III (depicted in blue) mimics the scenario of additional correlated *truly observed* phenotypes being used as additional predictors together with the predictive DNA markers, which corresponds to the strategy 3 empirical model (see blue bars in Fig. 2). (A) Prediction accuracy evaluated by R² for continuous phenotypes. (B) Prediction accuracy evaluated by AUC for binary phenotypes (dichotomized from the continuous ones by mean values) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
View Large Image
Figure Viewer
Download Hi-res image
Download (PPT)

4. Discussion

In the present study, we explored, quantified, and now understand via computer simulations, the extent of which correlated appearance phenotypes and underlying genotypes impact on DNA-based appearance prediction using pigmentation traits as example.

We empirically demonstrated mid-range and highly statistically significant correlations between all three pigmentation phenotypes, and quantified a large contribution on these phenotype correlations provided by the 41 HIrisPlex DNA markers. That the 41 HIrisPlex-S SNPs could not explain 100 % of the observed phenotypic correlations is likely explained by phenotyping classifications and additional genetic factors besides the 41 SNPs used here, e.g., SNPs in LD within the same regions and/or other unknown but shared genetic factors, and perhaps by additional non-genetic i.e., environmental factors not considered here. For instance, that considerably lower genetically explained correlations were achieved for the two phenotype comparisons that involve hair colour, compared to the one not involving hair colour, might be influenced by age-dependent hair colour change, which is not covered by the DNA markers used as previously reported for the 22 HIrisPlex hair colour SNPs [

Walsh S.
Liu F.
Wollstein A.
Kovatsi L.
Ralf A.
Kosiniak-Kamysz A.
et al.

The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA.

], and thus not expected to be covered by the total of 41 HIrisPlex-S markers used here including these 22. Although not tested here, we expect that these estimates may decrease when applied to homogeneous population samples (e.g. Europeans only) due to their reduced phenotypic variance.

In our empirical genetic prediction modelling based on the full set of 41 HIrisPlex-S DNA markers, we found no effect on prediction accuracy of DNA-predicted correlated pigmentation phenotypes used as additional predictors. In contrast a noticeable effect, mostly in increasing prediction accuracy, was seen when using truly observed correlated pigmentation phenotypes as additional predictors. The results from our computer simulations suggest that the prediction accuracy improvement achieved by considering truly observed correlated pigmentation phenotypes is explained by the extra genetic (and non-genetic) information that is provided by these correlated phenotypes. Such extra information is not provided when correlated pigmentation phenotypes were predicted from a limited number of DNA markers such as the 41 HIrisPlex-S markers used here. This explains why the prediction accuracies were almost identical or very similar, when comparing the genetic prediction models without considering phenotype correlations, and those that used DNA-predicted correlated phenotypes as additional predictors. Thus, it is not the prediction error of the DNA-predicted phenotypes that does not allow the prediction accuracy to improve, as one may naively assume. Instead, as our computer simulations suggest, it is the incomplete genetic information covered by the DNA markers used to predict the respective correlated pigmentation phenotype that does not provide a prediction improvement. This missing genetic information (together with unknown non-genetic information) is covered by the truly observed correlated phenotypes used as additional predictor, which in turn provides the prediction accuracy increase.

Our results are relevant for a variety of fields where genetic appearance prediction models are used in practise, such as forensic casework, human history studies, and anthropological investigations, all of which typically lack any phenotype information, including phenotypes that are correlated with the phenotype of interest. In particular, our findings imply that it is not advisable in such practical applications to use DNA-predicted correlated phenotypes as extra factors in the prediction modelling, at least not when it comes to pigmentation prediction from a relatively small set of DNA markers such as the HIrisPlex-S DNA markers used here. Moreover, the conclusion from our simulation analysis provide valuable information for future research and development of better prediction models as it demonstrates that using more independently contributing genetic predictors, in addition to the currently used ones, will increase prediction accuracy overall. However, our findings do not provide any information on how many additional DNA predictors would be needed to achieve the prediction improvement seen when using truly-observed correlated phenotypes, which in turn would depend on the independently contributing effect size of the additional DNA predictors.

Theoretical expectations and empirical demonstration of GWASs based on common DNA variants, such as for human pigmentation traits [

Hysi P.G.
Valdes A.M.
Liu F.
Furlotte N.A.
Evans D.M.
Bataille V.
et al.

Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability.

] or body height [

[16]

Wood A.R.
Esko T.
Yang J.
Vedantam S.
Pers T.H.
Gustafsson S.
et al.

Defining the role of common variation in the genomic and biological architecture of adult human height.

], have demonstrated that newly identified significantly associated DNA variants obtained from GWASs with larger sample size have similarly small effects or even smaller effects than those identified in previous GWASs based on smaller sample size. Therefore, it is expected that many more (common) DNA predictors will be needed to further improve the prediction accuracy of already established appearance prediction models compared to the number of DNA markers included in the established models. The first empirical evidence for this notion for human pigmentation traits comes from a recently published large GWAS on hair colour that additionally reports on the use of the many newly discovered DNA-markers for hair colour prediction [

Hysi P.G.
Valdes A.M.
Liu F.
Furlotte N.A.
Evans D.M.
Bataille V.
et al.

Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability.

]. Based on a large European discovery set of 290.891 individuals, this study found 124 genetic loci significantly associated with hair colour, of which 111 were novel and not identified in previous GWASs. Moreover, the study demonstrated that a new hair colour prediction model based on 258 independently associated SNPs, most of them were newly discovered, improved the prediction accuracy for all hair colour categories, except red, compared to the previous HIrisPlex model. However, given the more than 10-times larger number of DNA predictors in the new model compared to the HIrisPlex model, the prediction accuracy increase this new model provided relative to HIrisPlex, i.e., AUC increase of 0.67 to 0.79 for blond, 0.66 to 0.76 for brown, and 0.82 to 0.91 for black, was not overwhelming. Nevertheless, in the era of targeted massively parallel sequencing (MPS), hundreds of SNPs can now be utilized for practical applications including forensic DNA analysis [

[40]

Ralf A.
van Oven M.
Montiel Gonzalez D.
de Knijff P.
van der Beek K.
Wootton S.
et al.

Forensic Y-SNP analysis beyond SNaPshot: high-resolution Y-chromosomal haplogrouping from low quality and quantity DNA using Ion AmpliSeq and targeted massively parallel sequencing.

], and MPS tool development for forensic DNA phenotyping using hundreds of SNPs is currently underway by the VISAGE Consortium.

Although, in our present study, we used correlated pigmentation traits and pigmentation-predictive DNA markers as example, similar findings may be expected for other correlated appearance traits. Such expectation only holds as long as the phenotype correlations are similarly high and the applied predictive DNA markers explain similarly large proportions of the phenotype correlations as seen here for eye, hair and skin colour and the HIrisPelx-S DNA markers. However, the success of DNA-based pigmentation prediction is driven by the presence of major gene effects (e.g. HERC2, SLC45A2 and MC1R) together with minor gene effects, whereas correlated non-pigmentation appearance traits such as facial shape phenotypes are characterized by the absence of major gene effects as far as currently known [