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

Prediction of human pigmentation traits, one of the most differentiable externally visible characteristics among individuals, from biological samples represents a useful tool in the field of forensic DNA phenotyping. In spite of freckling being a relatively common pigmentation characteristic in Europeans, little is known about the genetic basis of this largely genetically determined phenotype in southern European populations. In this work, we explored the predictive capacity of eight freckle and sunlight sensitivity-related genes in 458 individuals (266 non-freckled controls and 192 freckled cases) from Spain. Four loci were associated with freckling (MC1R, IRF4, ASIP and BNC2), and female sex was also found to be a predictive factor for having a freckling phenotype in our population. After identifying the most informative genetic variants responsible for human ephelides occurrence in our sample set, we developed a DNA-based freckle prediction model using a multivariate regression approach. Once developed, the capabilities of the prediction model were tested by a repeated 10-fold crossvalidation approach. The proportion of correctly predicted individuals using the DNA-based freckle prediction model was 74.13%. The implementation of sex into the DNA-based freckle prediction model slightly improved the overall prediction accuracy by 2.19% (76.32%). Further evaluation of the newly-generated prediction model was performed by assessing the model’s performance in a new cohort of 212 Spanish individuals, reaching a classification success rate of 74.61%. Validation of this prediction model may be carried out in larger populations, including samples from different European populations. Further research to validate and improve this newlygenerated freckle prediction model will be needed before its forensic application. Together with DNA tests already validated for eye and hair colour prediction, this freckle prediction model may lead to a substantially more detailed physical description of unknown individuals from DNA found at the crime scene.

Only variants with protein sequence alterations are shown. Previous in vitro or in silico functional analysis were used to classify MC1R variants into three categories. In the absence of previous published data, impact on receptor function was predicted by the PolyPhen programme. Genetic variants with PolyPhen scores lower than 0.50 were considered as p alleles (pseudoalleles with similar function compared to wild-type), while variants with scores from 0.50 to 0.95 were classified as r alleles (possibly damaging).

Association analysis
Association analyses were performed using the R software (http://www.R-project.org). All analyses performed were two-sided, and a significance level of 0.05 was considered for rejection of the null-hypothesis. The conservative Bonferroni correction was used to adjust the significance level for multiple testing (P-value < 4.54 × 10 −3 = 0.05/11). Unknown and missing values were excluded at each specific analysis.
For each polymorphism studied, Fisher's exact test was used to check for deviations from Hardy-Weinberg equilibrium (HWE) among controls. Minor allele frequencies (MAFs) for freckled and non-freckled individuals were estimated from our population.

Gan-Or et al. (2016)
Associations between the genotyped polymorphisms and the presence of ephelides were assessed according to the co-dominant model of inheritance via binary logistic regression. Genotype-related odds ratios (ORs), their corresponding 95% confidence intervals (CIs) and associated P-values were estimated. The proportion of the total variance in freckling explained by each genetic variant was estimated using Nagelkerke pseudo-R 2 statistic (R 2 ). This statistic parameter was used to rank the genetic variants included in the study based on their contribution to the freckling phenotype.

Prediction model
A multivariate logistic regression approach was applied to build the prediction models. The DNA-based freckle prediction model was developed using backward elimination of genetic variants based on the Akaike information criterion (AIC), being the optimal model the one with the smaller AIC valuebest balance between goodness-of-fit and parsimony. For each iteration, the lowest predictor in the variable set is excluded. Then, the model is rebuilt, used to predict again freckle occurrence, and the quality of the newly generated model is re-tested. This model-building process is repeated until all remaining polymorphisms included in the prediction model have a statistically significant contribution, and the estimated information loss is minimised.
To assess the influence of sex in the freckling phenotype, the DNA-based prediction model was additionally adjusted by including sex, as well as the interaction between sex and each genetic variant as covariates in the multivariate logistic regression. As above, the optimal sex-adjusted prediction model was developed based on the AIC variable selection criteria.
Finally, statistical interactions between genetic variants were examined using the MDR software (http://www.epistasis.org/). The multifactor dimensionality reduction method is a powerful strategy for detecting and interpreting statistical locus-locus epistasis [34]. Dendrogram interaction graphs provided by MDR illustrate the presence, strength, and nature of epistatic effects. Pairwise genetic interactions were tested by including interactions into the logistic regression model.

Model evaluation 2.7.1 Internal evaluation phasecross-validation
In order to evaluate the predictive ability of the freckle prediction models, we applied a repeated 10-fold cross-validation approach. Briefly, individuals were randomly divided into 10 equal data subsets. For each round of cross-validation, nine data subsets were used for performing the analysis (training set) while the remaining subset was used for validating the analysis (testing set). This process was run 10 times, using a different subset as testing set each time. This entire procedure was repeated 200 times, using different random partitions of the original sample to protect against chance divisions of the dataset. Data partition was carried out using the 'caret package' for R software, and the repeated cross-validation approach was performed with custom programmes written in R.
Receiver operating characteristic (ROC) curve analysis was adopted to evaluate the cross-validated performance of the prediction models. The Youden Index method was used to set the optimal cut-off point -the point on the curve at which (sensitivity + specificity−1) is maximised. For each prediction model, basic prediction accuracy parameters -including the area under the ROC curve (AUC), sensitivity, specificity and accuracy -were calculated from the confusion table reporting the number of true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN). The final values reported were averaged from the 200 × 10 individual cross-validation curves to produce a single estimation. All calculations and plots were performed using the 'ROCR package' for R software.

External evaluation phaseindependent population
In order to carry out the external evaluation, only the genetic variants that were significantly associated with freckles in the first phase of the study -IRF4, ASIP, BNC2 and MC1R -were analysed in the newly collected independent samples.
This independent dataset was used to assess the performance of the newly-generated freckle prediction model regarding the correct classification of individuals as freckled or non-freckled. Subsequently, a comparison between predicted and observed data was performed, and the basic accuracy parameters (AUC, sensitivity, specificity and accuracy) were calculated from the confusion table by using the 'pROC package' for R software.

Results
Our original sample set comprised 266 non-freckled and 192 freckled individuals of Spanish ancestry ( Table 2). The study cohort included 260 females (56.52%) and 198 males (43.04%). Interestingly, the prevalence of the freckling phenotype was remarkably higher in females than in males (68.75% vs 31.25%, P-value = 9.00 × 10 −9 ). As expected, freckled individuals were more likely to have skin phototypes I-II and childhood sunburns than non-freckled individuals, although no significant differences between them were found regarding sun exposure habits ( Table 2). For each polymorphism, minor allele frequencies in freckled and non-freckled individuals from our population are detailed in Table 3. Evidence of departure from Hardy-Weinberg equilibrium among controls was observed for two SNPs: rs12896399 in SLC24A4 gene and rs16891982 in SLC45A2 -as seen in previous studies involving these human pigmentation genes [35][36][37]. Out of 458 individuals studied, 277 (60.48%) carried at least one non-synonymous MC1R variant, including 135 (54.43%) in a total of 248 non-freckled individuals, and 142 (75.93%) in a total of 187 freckled individuals. All genetic variants found in the MC1R gene are listed and classified in Table 1. To our knowledge, two mutations (G32R and L135R) have never been previously described, and they add up to the growing list of rare MC1R natural variants. The frequencies of R and r alleles in both freckled and non-freckled individuals are detailed in Table 3. Seventeen samples were discarded due to unsuccessful MC1R sequencing. Table 3 Genotypic association with freckles in the Spanish population.
alt-text: Table 3 Gene The proportion of total variance in freckling phenotype explained by the genetic variants was estimated using R 2 .
c The genetic variants were ranked according to their importance in freckling phenotype.
Univariate association analyses were performed to assess the independent effects of each genetic variant on freckling occurrence in childhood (Table 3). R variants of the MC1R gene were the alleles most strongly associated with freckling (P-value = 7.05 × 10 −8 ), explaining 6.91% of the variance. The rs12203592 polymorphism in the IRF4 gene presented the second strongest association with the presence of ephelides (P-value = 1.81 × 10 −5 ), explaining 5.92% of the variance. Weaker significant associations were also observed for r variants (P-value = 0.039) and rs4911442 in ASIP (P-value = 0.044), explaining around 1% of the variance (R 2 = 1.09% and R 2 = 1.01%, respectively). No associations were found for the remaining six polymorphisms, although rs2153271 in BNC2 showed a marginal significant association with freckling (P-value = 0.067).
A prediction model was then constructed based on multinomial logistic regression. To design the optimal model, we performed a step-wise analysis by iteratively excluding the lowest predictor from the multinomial logistic regression model. The genetic variants were ranked according to their impact on freckling. The DNA-based prediction model included four loci (MC1R, IRF4, ASIP and BNC2), which together explained about 30% of the freckling phenotype variance in our population (Table 4). A repeated 10-fold cross-validation approach was used to perform an internal validation of the DNA-based prediction model by ROC analysis, which tests the power of the prediction model [38]. The prediction accuracy of the tested DNA-based prediction model was 74.13% in the population, with an AUC of 0.771, a specificity of 82.00%, and a sensitivity of 63.51% (Fig. 1A). Table 4 Multivariate logistic regression testing freckle association with genetic variants and sex.
alt-text: Table 4 DNA-based prediction model Sex-adjusted prediction model Optimal model selection procedure was based on the Akaike information criterion (AIC).
a P-value for the multivariate logistic regression association analysis. Bold indicates significant results at Bonferroni threshold of 4.50E-03.
The inclusion of sex on the freckle prediction model appeared to increase the prediction accuracy by 2.19%, to a total of 76.32% (Fig. 1B). Slight increases in AUC (0.781), specificity (83.80%) and sensitivity (66.21%) values were observed when sex was taken into account. The results of the multivariate association analysis confirmed sex as an important variable in the estimation of the freckling phenotype (P-value = 2.70 × 10 −3 ). No significant alt-text: Fig. 1 interaction between sex and other genetic predictors was found.
An external validation of the newly-generated prediction model, assessing the model generalisability [39], was also performed by collecting an independent cohort of 212 individuals -though only 193 individuals successfully genotyped for all loci were included in the external validation. The prediction accuracy of the freckle prediction model was 74.61% in the independent cohort, with an AUC of 0.809, a specificity of 89.12%, and a sensitivity of 59.79%. A total of 49 individuals were incorrectly predicted (misclassification rate of 25.39%) ( Table 5). We also evaluated the prediction capacity of the sex-adjusted prediction model to determine the freckling phenotype in adulthood. Taking into account that freckles often disappear after adolescence or during young adulthood, the success rate of the prediction model was relatively good in our population (69.63%), with a high true negative rate (specificity of 80.70%) but a low true positive rate (sensitivity of 58.91%) (Fig. 2). Out of all six freckle predictors included in the model, only the r variants in MC1R did not show a significant contribution to freckling in adulthood (P-value > 0.05).

Fig. 2
Area under the receiver operating characteristic curves for the sex-adjusted freckle prediction model in adulthood. Individuals from both the first (model building and internal validation) and the second phase (external validation) of the study were pooled to perform this analysis. Using Youden Index, optimal cut-off was selected for maximising the sensitivity and specificity of the prediction model to classify individuals. AUC, area under the curve.
alt-text: Fig. 2 MDR analysis for locus-locus interaction detection indicated that the best freckle prediction model was composed by four genetic variables and sex (BNC2 rs2153271 was not included). The locus-locus interaction analysis showed a redundant interaction between R variants in the MC1R gene and rs12203592 in the IRF4 gene, denoted by the blue lines connecting these two genetic determinants in the dendrogram (Fig. S2 in Supplementary material in online version at DOI: 10.1016/j.fsigen.2017.11.013). However, the inclusion of this genetic interaction in our freckle prediction model did not significantly modify the prediction capacity using a repeated 10-fold cross-validation approach (AUC = 0.768, specificity = 83.20%, sensitivity = 61.27%, and accuracy = 75.45%).

Discussion
In the current study, eight freckle and sunlight sensitivity-related genes were genotyped in 458 individuals from Spain, with the intention of analysing their putative implication in the appearance of ephelides, and its potential for future forensic applications. However, only the 438 individuals successfully genotyped for all loci were included in the development and cross-validation of the prediction model. Furthermore, an external evaluation of the freckle prediction model was performed on an extra 212 sample independent dataset by genotyping the genetic variants associated with freckling in the original population.
In the last few years, the genetic basis of ephelides has been adequately studied in populations of North European origin [11,[16][17][18]. However, little is known about the genetic determinants of this pigmentation characteristic in southern European populations. Our results showed an association between ephelides and genetic variants in the IRF4, MC1R and ASIP genes in a Mediterranean population, confirming previous studies performed in North European populations [11,17,18]. R variants of the MC1R gene ha ve been previously acknowledged as the most relevant locus associated with both freckling and sun sensitivity. The abnormal function of this receptor leads to a higher phaeomelanin to eumelanin ratio, commonly resulting in the well-known freckle-generating RHC phenotype. IRF4 has also been associated with freckling in several studies [8,18]. This interferon regulatory factor cooperates with MITF to activate the expression of tyrosinase in melanocytes, a function which seems to be impaired in carriers of the rs12203592*T derived allele [40]. The other main freckle-associated loc us in this study, ASIP, antagonises the activation of MC1R, leading to a down-regulation of eumelanogenesis and an up-regulation of phae omelanogenesis. Variants in this gene, also previously linked to the RHC phenotype, may therefore have effects similar to variants of the MC1R gene [16,17].
The impact of BNC2 rs2153271 on freckling was detected only after applying a multivariate association approach. The identification of the association between freckling and rs2153271, located in an intron of the BNC2 gene, was recently discovered in a GWAS study performed in a population of North European ancestry [16], and was additionally correlated with acquired facial pigmented spots during aging [19].
Although rs1042602 and rs1393350, in the TYR gene, have been previously associated with freckling in two European populations of Icelandic and Dutch origin [18], a lack of association between these two SNPs and ephelides occurrence was observed in our Spanish population. Accordingly, no association with freckling was previously observed for a TYR melanoma-associated genetic variant (rs1126809) in a sample set comprising three Mediterranean populations from France, Spain and Italy [36]. However, that melanoma case-control study identified a moderate association between rs16891982 in SLC45A2 and the absence of ephelides [36]. The inheritance of a set of genetic variants in SLC45A2 was also associated with freckling in a population from Brazil [41]. However, we are not able to verify this association in the current study, perhaps due to the limited number of individuals included in the sample set.
This association study allows us to develop a freckle prediction model based on five genetic predictors: R variants in MC1R, IRF4 rs12203592, r variants in MC1R, ASIP rs4911442 and BNC2 rs2153271 (listed in order, from greater to lower genetic contribution). Predicting externally visible human traits from genotypes represents a potential valuable tool in forensic investigations [1]. However, DNA-based phenotyping needs a very demanding approach since externally visible characteristics are typically influenced by several genetic as well as environmental factors [42]. To date, DNA-based human eye colour prediction is the most accurate, advanced and applied test in forensic applications [3][4][5]43]. The IrisPlex system, a robust prediction model based on six eye colour SNP predictors, has been validated in numerous studies for forensic eye colour prediction [4,5,43]. Recently, a multiplex genotyping assay, called HIrisPlex, has been developed for simultaneous hair and eye colour prediction [6], being suitable and sufficiently sensitive for forensic use approval [44]. Currently, forensic genome-wide association studies are focused on increasing the genetic determinants of different phenotypes such as body height, hair shape, baldness, and facial variation, in order to reach more detailed predictions of an unknown person's appearance from the analysis of his/her DNA [16,[45][46][47][48].
As a whole, the DNA-based freckle prediction model presents with high specificity (82.00%), but low sensitivity (63.51%). The low sensitivity levels achieved could be due to the fact that one third of the individuals carrying R variants did not display freckles, suggesting that other genetic determinants are also important in the appearance of freckles. Hopefully, future studies with increased sample sizes may add to the genetic factors influencing freckle occurrence, so that the model's prediction potential may be significantly improved.
Interestingly, we found a much higher prevalence of ephelides in females than in males. Previous studies have also stated discrepancies in human pigmentation and sunlight sensitivity traits between sexes [19,35,[49][50][51][52].
Notably, higher freckle prevalence in females was previously observed in a GWAS study performed in an northern European population [18]. Additionally, Jacobs and cols. (2015) showed that females presented a much higher occurrence of facial sun spots than males, being the total variance explained by sex higher than any of the genetic variants studied [19].
s to to i n The inclusion of sex in the IrisPlex model has been proposed to improve prediction accuracy mainly for intermediate eye colours [35,50,53,54]. Adding sex as a covariate in our DNA-based freckle prediction model slightly increased the prediction performance for ephelides prevalence. In particular, it presented with both higher specificity (83.80%) and sensitivity (66.21%) if compared to the prediction model based on the five freckle genetic predictors alone. Overall, this sex-adjusted freckle prediction model could explain a high proportion of the freckling phenotypic variance (R 2 = 32.92%), which is quite large compared to other human complex traits [55].
The main purpose of a prediction model is to provide correct classification for new individuals, being therefore external validation a crucial phase of the model development process. Assessing both reproducibility (internal validation) and generalisability (external validation) of the developed prediction model are necessary in order to avoid overfitted models [56]. However, application of a prediction model should only be considered after proving adequate accuracies in an independent dataset. For this reason, the discrimination ability of the newly-generated prediction model was assessed by comparing the observed and predicted outcomes in a new dataset. As usually noted in most external validations, a slight decrease of the overall prediction accuracy was observed (from 76.32% to 74.61%).
In terms of prediction power, our freckle prediction model provided a similar accuracy than the one obtained for green/hazel eye colour prediction using IrisPlex (AUC = 0.809 and AUC = 0.76, respectively), although significantly lower than the accuracy level achieved for brown (AUC = 0.93) and blue eye colour prediction (AUC = 0.91) [5]. HIrisPlex also presented equivalent prediction power for hair colour prediction, with an average accuracy of 73% [44], compared to 74.61% of our freckle prediction model. Notice that our results have been reached after a validation approach with a relatively small independent sample, while both IrisPlex and HIrisPlex prediction models have been validated in independent large population samples.
After our preliminary validation, the freckle prediction model developed in this study has shown to display limited success, but it also appears promising for future applications in forensics. However, further research is needed to increase the correct prediction rate (sensitivity) of the model, since the percentage of freckled individuals predicted as non-freckled is considerably high (false negative rate of 40.21%).
For practical considerations, knowing that ephelides tend to disappear during aging due to undetermined reasons, we tested the accuracy of our prediction model to determinate the presence of ephelides in adulthood. The Since the sample set used to build and evaluate this freckle prediction model was relatively small, the ability to accurately predict the presence of ephelides may be further improved by using larger studies. It is likely that phenotype prediction accuracy may be substantially improved by including undiscovered genetic variants at novel loci or by taking into account possible gene interactions affecting the freckling phenotype. Also, DNA samples from different populations of European origin should be included in further studies in order to validate this DNA predictive test for future forensic applications.
It is also important to note that predicting complex phenotypic traits from DNA studies remains a difficult task even if all genetic loci involved -and the interactions among them -are taken into account, since different environmental factors may always have a considerable effect (for example, the effect of UV exposure on freckle occurrence).

Conclusions
The main genetic variants involved in freckle appearance in the Spanish population are located in the MC1R and IRF4 genes, with minor contributions from ASIP, BNC2 and perhaps other as yet unknown freckle genes.
However, the influence on freckling of different MC1R variants (R or r) is substantially different. As a result, the DNA-based model for freckle prediction developed in this work considers five genetic determinants: R variants of the MC1R gene, IRF4 rs12203592, r variants of the MC1R gene, ASIP rs4911442 and BNC2 rs2153271 -in order of greatest to lowest contribution -reaching a cross-validated prediction accuracy of up to 74.13%, a more than respectable percentage. When sex is added to the model, the cross-validated prediction accuracy reached is even higher, growing up to 76.32%.
Furthermore, the newly-generated freckle model was tested in an independent cohort reaching an acceptable specificity level of 89.12%. However, the sensitivity of the model is certainly improvable, attaining a percentage of around 60%.
For a more detailed pigmentation phenotype prediction, future research may focus on designing multiplex genetic analysis for simultaneously predicting different externally visible traits. Perhaps the inclusion of two more SNPs in the HIrisPlex system, BNC2 rs2153271 and ASIP rs4911442, could increase the forensic potential of this DNA-based prediction test to include freckle occurrence (as well as the current eye and hair colour prediction).

E-Extra
For practical considerations, knowing that ephelides tend to disappear during aging due to undetermined reasons, we tested the accuracy of our prediction model to determin e the presence of ephelides in adulthood. The model presented an AUC of 0.756, a slightly lower value compared to the AUC value obtained for freckle prediction in childhood (0.781 for internal cross-validation and 0.809 for external validation).
No significant contribution to freckling prediction in adulthood was noted for r variants, suggesting that other unknown predictors could have an independent additional impact on the freckling phenotype in adulthood. Interestingly, after adjustment by IRF4 rs12203592 (the most associated genetic predictor for freckling in adulthood), the minor allele of SLC45A2 rs16891982 was weakly associated with the absence of freckles in adulthood (OR = 0.64 (0.43-0.94), P-value = 0.023). As revealed by the MDR analysis, the effect of this SLC45A2 variant may not be strongly influenced by a locus-locus interaction, since there is no interaction between SLC45A2 and the rest of the freckle loci (Fig. S2B).

E-component
Under the supervision of a professional, each participant completed a standardised questionnaire to collect information on sex, age, pigmentation traits, history of childhood sunburns, Fitzpatrick's skin type classification, and sun exposure habits. Details of ephelides occurrence both during the infancy or adolescence periods and in adulthood were collected in the questionnaire (an illustration is shown in Fig. S1).
MDR analysis for locus-locus interaction detection indicated that the best freckle prediction model was composed by four genetic variables and sex (BNC2 rs2153271 was not included). The locus-locus interaction analysis showed a redundant interaction between R variants in the MC1R gene and rs12203592 in the IRF4 gene, denoted by the blue lines connecting these two genetic determinants in the dendrogram (Fig. S2). However, the inclusion of this genetic interaction in our freckle prediction model did not significantly modify the prediction capacity using a repeated 10-fold cross-validation approach (AUC = 0.768, specificity = 83.20%, sensitivity = 61.27%, and accuracy = 75.45%).
at Multimedia Component 1