Advertisement

Genome-Wide Association Shows that Pigmentation Genes Play a Role in Skin Aging

Open ArchivePublished:May 11, 2017DOI:https://doi.org/10.1016/j.jid.2017.04.026
      Loss of fine skin patterning is a sign of both aging and photoaging. Studies investigating the genetic contribution to skin patterning offer an opportunity to better understand a trait that influences both physical appearance and risk of keratinocyte skin cancer. We undertook a meta-analysis of genome-wide association studies of a measure of skin pattern (microtopography score) damage in 1,671 twin pairs and 1,745 singletons (N = 5,087) drawn from three independent cohorts. We identified that rs185146 near SLC45A2 is associated with a skin aging trait at genome-wide significance (P = 4.1 × 10–9); to our knowledge this is previously unreported. We also confirm previously identified loci, rs12203592 near IRF4 (P = 8.8 × 10–13) and rs4268748 near MC1R (P = 1.2 × 10–15). At all three loci we highlight putative functionally relevant SNPs. There are a number of red hair/low pigmentation alleles of MC1R; we found that together these MC1R alleles explained 4.1% of variance in skin pattern damage. We also show that skin aging and reported experience of sunburns was proportional to the degree of penetrance for red hair of alleles of MC1R. Our work has uncovered genetic contributions to skin aging and confirmed previous findings, showing that pigmentation is a critical determinant of skin aging.

      Abbreviations:

      BG6 (Beagley-Gibson six-point rating system), eQTL (expression quantitative trait locus), GWAS (genome-wide association study), LD (linkage disequilibrium), QC (quality control), SL (solar lentigines), SNP (single-nucleotide polymorphism), TEST (Twins Eyes Study in Tasmania)

      Introduction

      Skin aging, characterized by epidermal thinning and reduced DNA repair, results from both external and internal influences. Major extrinsic factors include exposure to UVR, causing photoaging, and smoking. Specifically, photoaging is typified by elastosis and damage to the complex extracellular support matrix of the dermis (for review, see
      • Bilac C.
      • Sahin M.T.
      • Ozturkcan S.
      Chronic actinic damage of facial skin.
      ). This photodamage alters mechanical properties in the skin, leading to deterioration and loss of regularity of the fine reticular patterning of the skin surface (
      • Lavker R.M.
      • Kwong F.
      • Kligman A.M.
      Changes in skin surface patterns with age.
      ), as well as wrinkles and sagging. Because UVR also induces DNA damage, the degree of photoaging correlates with skin carcinogenesis (
      • Gandini S.
      • Sera F.
      • Cattaruzza M.S.
      • Pasquini P.
      • Zanetti R.
      • Masini C.
      • et al.
      Meta-analysis of risk factors for cutaneous melanoma: III. Family history, actinic damage and phenotypic factors.
      ,
      • Green A.C.
      Premature ageing of the skin in a Queensland population.
      ,
      • Holman C.D.
      • Armstrong B.K.
      • Evans P.R.
      • Lumsden G.J.
      • Dallimore K.J.
      • Meehan C.J.
      • et al.
      Relationship of solar keratosis and history of skin cancer to objective measures of actinic skin damage.
      ). Thus, heritable factors that modify skin aging, particularly photoaging, have implications for both personal appearance and risk of skin cancer.
      To date, a number of genome-wide association studies (GWASs) have analyzed various measures of photoaging or overall skin aging. These include perceived facial age (
      • Liu F.
      • Hamer Merel A.
      • Deelen J.
      • Lall Japal S.
      • Jacobs L.
      • van Heemst D.
      • et al.
      The MC1R gene and youthful looks.
      ), facial photoaging (
      • Le Clerc S.
      • Taing L.
      • Ezzedine K.
      • Latreille J.
      • Delaneau O.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women points out a putative role of the STXBP5L gene in facial photoaging.
      ), and the combined impact of intrinsic and extrinsic factors on facial aging (
      • Chang A.L.
      • Atzmon G.
      • Bergman A.
      • Brugmann S.
      • Atwood S.X.
      • Chang H.Y.
      • et al.
      Identification of genes promoting skin youthfulness by genome-wide association study.
      ). Additional GWASs have been conducted on traits related to photoaging such as actinic keratosis and facial pigmented spots (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ,
      • Jacobs L.C.
      • Liu F.
      • Pardo L.M.
      • Hofman A.
      • Uitterlinden A.G.
      • Kayser M.
      • et al.
      IRF4, MC1R and TYR genes are risk factors for actinic keratosis independent of skin color.
      ,
      • Laville V.
      • Le Clerc S.
      • Ezzedine K.
      • Jdid R.
      • Taing L.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women suggests the involvement of HLA genes in the severity of facial solar lentigines.
      ). Together, these have identified genetic variants on chromosomes 3, 6, 16, 20, and X.
      Here, we report on a meta-analysis of five GWASs drawn from three independent cohorts (Methods) analyzed using the six-point Beagley and Gibson (BG6) microtopography scoring system of skin patterning regularity and complexity (
      • Beagley J.
      • Gibson I.M.
      Changes in skin condition in relation to degree of exposure to ultraviolet light.
      ,
      • Holman C.D.
      • Armstrong B.K.
      • Evans P.R.
      • Lumsden G.J.
      • Dallimore K.J.
      • Meehan C.J.
      • et al.
      Relationship of solar keratosis and history of skin cancer to objective measures of actinic skin damage.
      ). BG6 score strongly correlates with amounts of dermal solar elastosis (
      • Battistutta D.
      • Pandeya N.
      • Strutton G.M.
      • Fourtanier A.
      • Tison S.
      • Green A.C.
      Skin surface topography grading is a valid measure of skin photoaging.
      ,
      • Fritschi L.
      • Battistutta D.
      • Strutton G.M.
      • Green A.
      A non-invasive measure of photoageing.
      ,
      • Seddon J.M.
      • Egan K.M.
      • Zhang Y.
      • Gelles E.J.
      • Glynn R.J.
      • Tucker C.A.
      • et al.
      Evaluation of skin microtopography as a measure of ultraviolet exposure.
      ), more so at younger ages (
      • Hughes M.C.
      • Strutton G.M.
      • Fourtanier A.
      • Green A.C.
      Validation of skin surface microtopography as a measure of skin photoaging in a subtropical population aged 40 and over.
      ). It is correlated with self-reported sun exposure history and skin phototype, and sunscreen can protect against the age-related increase in BG6 score (
      • Hughes M.C.
      • Williams G.M.
      • Baker P.
      • Green A.C.
      Sunscreen and prevention of skin aging: a randomized trial.
      ). At age 12 years, BG6 is highly heritable (86%), with heritability decreasing to 62% in adults (
      • Shekar S.N.
      • Luciano M.
      • Duffy D.L.
      • Martin N.G.
      Genetic and environmental influences on skin pattern deterioration.
      ), suggesting that BG6 is a useful measure to explore the genetics of photoaging. Here, we aim to confirm previous GWAS findings, identify additional loci, and characterize the interaction of MC1R red hair alleles with BG6 and various sun exposure behaviors.

      Results

      Skin patterning GWAS meta-analysis

      We performed a meta-analysis of five BG6 GWASs drawn from three independent cohorts, representing 1,671 twin pairs and 1,745 singletons (N = 5,087) (see Methods and Table 1). To account for differences in age, genotyping array, and phenotype, a total of five BG6 GWASs (Table 1) were conducted using imputed genotype dosage scores in Merlin-offline (
      • Abecasis G.R.
      • Cherny S.S.
      • Cookson W.O.
      • Cardon L.R.
      Merlin—rapid analysis of dense genetic maps using sparse gene flow trees.
      ), a method that appropriately models the family structure in related samples, including twins (see Methods and Figures 1 and 2). GWAS quantile-quantile plots indicated that there was no evidence for population stratification (see Supplementary Figure S1 online). Single-nucleotide polymorphisms (SNPs) at or near SLC45A2 reached genome-wide significance with a measurement of skin aging, as did SNPs near IRF4 and MC1R (Table 2; Supplementary Figures S2, S3, and S4 online). Full results for all SNPs with P less than 1 × 10–5 are reported in Supplementary Table S1 online. Conditioning on the top SNPs did not show additional associations (see Supplementary Methods, Supplementary Table S2, and Supplementary Figure S5 online).
      Table 1Population Description
      DatasetAdolescentAdultTESTRaine
      Phase 1Phase 2
      MZ, pairs38419460190
      DZ, pairs63429246420
      Singleton twins, n1105517080
      Siblings, n409167060
      Unrelated, n0000820
      Total, n2,5551,194382136820
      Total families, n1,06150419771820
      Mean age in years (SD)15.3 (1.9)14.4 (2.0)44.6 (9.8)9.6 (3.0)22.2 (0.6)
      Female:male1,304:1,251649:545228:15469:67412:408
      Abbreviations: DZ, dizygotic twins; SD, standard deviation; MZ, monozygotic twins; TEST, Twin Eyes Study in Tasmania.
      Figure 1
      Figure 1Quantile-quantile plot of –log10 of fixed-effects P-values derived from the meta-analysis of BG6 GWASs. The 95% confidence interval for the expected –log10(P) is indicated by black lines. BG6, Beagley-Gibson six-point rating system; GIF, genomic inflation factor; GWAS, genome-wide association study.
      Figure 2
      Figure 2Manhattan plot of –log10 fixed-effects P-values derived from the meta-analysis of BG6 GWASs. The negative log10 of meta-analysis P-values are plotted by chromosome. The orange line indicates genome-wide significance (P = 5 × 10–8) while the light blue line is P = 5 × 10–5. The top SNP and gene are indicated for loci reaching genome-wide significance. BG6, Beagley-Gibson six-point rating system; GWAS, genome-wide association study; SNP, single nucleotide polymorphism.
      Table 2Loci with SNPs with P-values less than 5 × 10–8 after the fixed effects meta-analysis
      We have reported the top SNP, and if known, additional functional SNPs at the same loci. The full results, including random effects estimates and individual genome-wide association study data, for all SNPs with P-values less than 1 × 10–5 can be found in Supplementary Table S1. Note that Table 1 results are given for the minor allele for ease of interpretation; as indicated in Supplementary Table S1, effect sizes may be reported for the other allele. SNP positions are Hg19. A1 Frequency is the frequency of allele 1 in the combined dataset, weighted by sample number. min rˆ2 is the minimum imputation quality score across contributing genome-wide association studies. Beta is the change in BG6 score for each copy of the A1 allele. A negative effect size indicates that allele 1 is associated with reduced BG6 scores and reduced skin aging. I2 is the measure of heterogeneity of effect sizes across the included studies, and a low I2 score indicates that the signals are homogenous.
      Hg19 PositionSNPGeneA1 Effect AlleleA2A1 Frequencymin rˆ2Fixed PFixed βI2
      5:33952106rs185146SLC45A2CT0.060.524.1 × 10–9–0.250
      6:396321rs12203592IRF4AG0.210.908.8 × 10–130.130
      rs12203592 did not pass quality control in the Raine study; all other SNPs in Table 1 are reported for all contributing studies.
      16:90026512rs4268748MC1RCT0.290.951.2 × 10–150.131.8
      16:89986117rs1805007MC1RTC0.090.831.2 × 10–100.180
      16:89986144rs1805008MC1RTC0.080.981.1 × 10–50.1224.6
      Abbreviations: A, adenine; C, cytosine; G, guanine; min, minimum; SNP, single-nucleotide polymorphism; T, thymine.
      1 We have reported the top SNP, and if known, additional functional SNPs at the same loci. The full results, including random effects estimates and individual genome-wide association study data, for all SNPs with P-values less than 1 × 10–5 can be found in Supplementary Table S1. Note that Table 1 results are given for the minor allele for ease of interpretation; as indicated in Supplementary Table S1, effect sizes may be reported for the other allele. SNP positions are Hg19. A1 Frequency is the frequency of allele 1 in the combined dataset, weighted by sample number. min rˆ2 is the minimum imputation quality score across contributing genome-wide association studies. Beta is the change in BG6 score for each copy of the A1 allele. A negative effect size indicates that allele 1 is associated with reduced BG6 scores and reduced skin aging. I2 is the measure of heterogeneity of effect sizes across the included studies, and a low I2 score indicates that the signals are homogenous.
      2 rs12203592 did not pass quality control in the Raine study; all other SNPs in Table 1 are reported for all contributing studies.
      In addition to BG6, a number of other skin aging-related traits have been analyzed by GWAS (
      • Chang A.L.
      • Atzmon G.
      • Bergman A.
      • Brugmann S.
      • Atwood S.X.
      • Chang H.Y.
      • et al.
      Identification of genes promoting skin youthfulness by genome-wide association study.
      ,
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ,
      • Jacobs L.C.
      • Liu F.
      • Pardo L.M.
      • Hofman A.
      • Uitterlinden A.G.
      • Kayser M.
      • et al.
      IRF4, MC1R and TYR genes are risk factors for actinic keratosis independent of skin color.
      ,
      • Laville V.
      • Le Clerc S.
      • Ezzedine K.
      • Jdid R.
      • Taing L.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women suggests the involvement of HLA genes in the severity of facial solar lentigines.
      ,
      • Le Clerc S.
      • Taing L.
      • Ezzedine K.
      • Latreille J.
      • Delaneau O.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women points out a putative role of the STXBP5L gene in facial photoaging.
      ,
      • Liu F.
      • Hamer Merel A.
      • Deelen J.
      • Lall Japal S.
      • Jacobs L.
      • van Heemst D.
      • et al.
      The MC1R gene and youthful looks.
      ). Supplementary Table S3 online summarizes the BG6 association for those SNPs previously reported as reaching genome-wide significance with other skin aging-related traits. Of note, the top SNP from the chromosome 3 region reported by
      • Le Clerc S.
      • Taing L.
      • Ezzedine K.
      • Latreille J.
      • Delaneau O.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women points out a putative role of the STXBP5L gene in facial photoaging.
      as associated with facial photoaging; rs322458 has a P-value of 0.07 for BG6. More promising is the nearby SNP rs470647, noted as a possible functional SNP tagged by rs3322458 (BG6 P = 3.8 × 10–3).

      Gene-based approach

      We used the gene-based association test VEGAS2 to test for multiple independent signals within a gene that are individually not genome-wide significant (
      • Mishra A.
      • Macgregor S.
      VEGAS2: software for more flexible gene-based testing.
      ) (see Supplementary Methods). All results with a P-value less than 5 × 10–5 are reported in Supplementary Table S4 online. Although many genes at the 16q24.3 locus have significant gene-based P-values (P < 2.084 × 10–6), because of extensive linkage disequilibrium (LD) across this region, this is driven by red hair alleles in MC1R.

      Functional annotation

      Because the top SNP at each locus may not be the functional one, meta-analysis results were filtered to P-values less than 1 × 10–7 and homogenous effects sizes across contributing studies (I2 < 31%), leaving 102 SNPs mapping to the three genome-wide significant regions on chromosomes 5, 6, and 16 (see Supplementary Methods and Supplementary Table S5 online). Our annotation work did not identify any additional regulatory or functional SNPs at 5p13.2/SLC45A2 beyond that our peak SNPs are in LD (r2 = 0.73) with a known missense rs16891982 (F374L) variant. The L374 allele, rare in European-derived populations, is associated with darker skin (
      • Graf J.
      • Hodgson R.
      • van Daal A.
      Single nucleotide polymorphisms in the MATP gene are associated with normal human pigmentation variation.
      ) and here with reduced BG6 score, indicating less photoaging (Table 2). The association at 6p25.3, near IRF4, is strongest at rs12203592, which lies in a melanocyte enhancer and is a weak IRF4 expression quantitative trait locus (eQTL) in whole blood and lymphocytes (Genotype-Tissue Expression [GTEx] project [
      GTEx Consortium
      The Genotype-Tissue Expression (GTEx) project.
      ] P = 6.0 × 10–7and P = 5.0 × 10–7, respectively). rs12203592 has been shown to modulate IRF4 regulation (
      • Praetorius C.
      • Grill C.
      • Stacey S.N.
      • Metcalf A.M.
      • Gorkin D.U.
      • Robinson K.C.
      • et al.
      A polymorphism in IRF4 affects human pigmentation through a tyrosinase-dependent MITF/TFAP2A pathway.
      ). In the MuTHER dataset (
      • Nica A.C.
      • Parts L.
      • Glass D.
      • Nisbet J.
      • Barrett A.
      • Sekowska M.
      • et al.
      The architecture of gene regulatory variation across multiple human tissues: the MuTHER study.
      ,
      • Grundberg E.
      • Small K.S.
      • Hedman A.K.
      • Nica A.C.
      • Buil A.
      • Keildson S.
      • et al.
      Mapping cis- and trans-regulatory effects across multiple tissues in twins.
      ), 17 BG6-associated SNPs at 16q24.3 are eQTLs for seven genes (see Supplementary Methods and Supplementary Table S6 online), and an overlapping list of 16 SNPs are eQTLs for 10 genes in the Westra data (
      • Westra H.J.
      • Peters M.J.
      • Esko T.
      • Yaghootkar H.
      • Schurmann C.
      • Kettunen J.
      • et al.
      Systematic identification of trans eQTLs as putative drivers of known disease associations.
      ) (see Supplementary Table S7 online). As reported by
      • Grundberg E.
      • Small K.S.
      • Hedman A.K.
      • Nica A.C.
      • Buil A.
      • Keildson S.
      • et al.
      Mapping cis- and trans-regulatory effects across multiple tissues in twins.
      these include the MC1R red hair SNP rs1805007 as an eQTL for DBNDD1. Together, these show direct impacts on pigmentation via nonsynonymous changes to SLC45A2 and MC1R, as well as expression modulation of IRF4. Annotation also highlights that MC1R red hair alleles may also be associated with altered expression in additional genes, which may further influence skin patterning and photoaging.

      Impact of MC1R alleles on sun exposure effects

      Our GWAS meta-analysis provides evidence that red hair alleles of MC1R are associated with skin patterning. 1,246 of the adolescent twins used in the adolescent phase 1 GWAS (Table 1) have been previously genotyped for nine low-frequency–coding MC1R alleles (
      • Duffy D.L.
      • Box N.F.
      • Chen W.
      • Palmer J.S.
      • Montgomery G.W.
      • James M.R.
      • et al.
      Interactive effects of MC1R and OCA2 on melanoma risk phenotypes.
      ). To improve our statistical power, we grouped these MC1R alleles by their penetrance for red hair and tested them for association with a range of sun exposure behaviors; a summary of these measures and their correlation can be found in the Supplementary Methods and Supplementary Tables S8 and S9 online. Male and female participants were analyzed separately given the influence of sex on skin patterning (see Supplementary Methods).
      There were significant differences in BG6 for MC1R compound genotypes (having any combination of two red hair alleles at MC1R) in both sexes (P < 0.001, Figure 3), with a significant trend for increasing BG6 score and greater red hair penetrance (P < 0.001). Together MC1R red hair alleles explain at least 4.12% of BG6 variation. rs12203592 in IRF4 and rs16891982 in SLC45A2 explain a further 2.95% and 0.26% of BG6 variation respectively. Sunburn risk increased for compound genotypes in line with the degree of penetrance for red hair (Figure 4). Additionally, we explored the impact of the MC1R genotype and thus vulnerability to skin damage on sun exposure behavior and use of sun protection (see Supplementary Methods). We concluded that there were no differences in sun-exposed hours (UV index justified, see Supplementary Figure S6 online) or sun protection behavior between the MC1R genotype groups (see Supplementary Figure S7 online). The genome-wide significant SNP rs12203592 in IRF4 did not exhibit any association with sun exposure or protection traits (data not shown).
      Figure 3
      Figure 3Adolescent BG6 scores by MC1R compound genotype. Error bars are mean ± standard error of the mean. The y-axis is BG6 score and x-axis is MC1R compound genotype. R indicates strongly penetrant red hair alleles (D84E, R151C, R160W, D294H), r weakly penetrant (V60L, V92M, R163Q), and + is the consensus haplotype. For each sex, the first F test result represents the test for significant difference in mean BG6 score across the compound genotype categories. A second F statistic and P-value is provided for the test for trend across all seven compound genotypes. BG6, Beagley-Gibson six-point rating system; r, weakly penetrant red hair alleles; R, strongly penetrant red hair alleles.
      Figure 4
      Figure 4Self-report sunburn scores by compound MC1R genotype. Error bars are mean ± standard error of the mean. The y-axis is Sunburn Score () and x-axis is MC1R compound genotype. R indicates strongly penetrant red hair alleles (D84E, R151C, R160W, D294H), r weakly penetrant (V60L, V92M, R163Q), and + is the consensus haplotype. For each sex, the first F test result represents the test for significant difference in mean BG6 score across the compound genotype categories. A second F statistic and P-value is provided for the test for trend across all seven compound genotypes. r, weakly penetrant red hair alleles; R, strongly penetrant red hair alleles.

      Discussion

      We show that SLC45A2 SNPs in LD with the nonsynonymous rs16891982 (F374L) (
      • Graf J.
      • Hodgson R.
      • van Daal A.
      Single nucleotide polymorphisms in the MATP gene are associated with normal human pigmentation variation.
      ) are genome-wide significantly associated with a measure of skin aging, BG6; to our knowledge, this is previously unreported. We also replicate previous associations at IRF4 and MC1R (see Supplementary Table S2) (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ,
      • Jacobs L.C.
      • Liu F.
      • Pardo L.M.
      • Hofman A.
      • Uitterlinden A.G.
      • Kayser M.
      • et al.
      IRF4, MC1R and TYR genes are risk factors for actinic keratosis independent of skin color.
      ,
      • Liu F.
      • Hamer Merel A.
      • Deelen J.
      • Lall Japal S.
      • Jacobs L.
      • van Heemst D.
      • et al.
      The MC1R gene and youthful looks.
      ). The associations include well-known alleles affecting pigmentation (
      • Duffy D.L.
      • Box N.F.
      • Chen W.
      • Palmer J.S.
      • Montgomery G.W.
      • James M.R.
      • et al.
      Interactive effects of MC1R and OCA2 on melanoma risk phenotypes.
      ). For each variant, the allele correlated with greater skin pigmentation was associated with reduced photoaging. A functional SNP rs470647 previously reported to be associated with facial photoaging (
      • Le Clerc S.
      • Taing L.
      • Ezzedine K.
      • Latreille J.
      • Delaneau O.
      • Labib T.
      • et al.
      A genome-wide association study in Caucasian women points out a putative role of the STXBP5L gene in facial photoaging.
      ) was associated (P = 3.8 × 10–3) with BG6. The detection of overlapping SNPs using a range of different measures of skin aging is promising (see Supplementary Table S3) and confirms that there are germline effects on multiple skin aging components. The combined effects of the seven MC1R red hair alleles polymorphic in our population explained at least 4.12% of BG6 variation, representing a gene of large effect compared with most quantitative trait loci.
      Replication across independent datasets is critical for confidence in GWAS findings. Data contributing to this meta-analysis were drawn from three different Australian cohorts and analyzed as five discrete GWASs to account for differences in genotyping array. BG6-associated SNPs at SLC45A2 and MC1R had P-values of less than 0.05 in three or more GWAS sets drawn from two independent cohorts; however, SNPs at IRF4 did not pass the imputation threshold in the Raine study (
      • Straker L.M.
      • Hall G.L.
      • Mountain J.
      • Howie E.K.
      • White E.
      • McArdle N.
      • et al.
      Rationale, design and methods for the 22 year follow-up of the Western Australian Pregnancy Cohort (Raine) Study.
      ) (see Supplementary Table S1 online). Nevertheless, the fact that we replicated skin aging association at MC1R and IRF4 suggests that our observations are robust. Although conditional analysis did not show additional genome-wide significant associations (see Supplementary Figure S5), a second red hair allele, rs1805008 (R160W) was also associated with BG6 at P = 1.1 × 10–5 (see Supplementary Table S1), suggesting that, like many other complex traits, additional variation remains to be discovered with larger sample sizes.
      The three loci associated with BG6 are well characterized, giving direct insights into how they influence fine skin patterning. Animal models suggest that SLC45A2 low-activity alleles impair the processing and transportation of melanogenic proteins to melanosomes (
      • Costin G.E.
      • Valencia J.C.
      • Vieira W.D.
      • Lamoreux M.L.
      • Hearing V.J.
      Tyrosinase processing and intracellular trafficking is disrupted in mouse primary melanocytes carrying the underwhite (uw) mutation. A model for oculocutaneous albinism (OCA) type 4.
      ), possibly via altered ability to maintain melanosome pH (
      • Dooley C.M.
      • Schwarz H.
      • Mueller K.P.
      • Mongera A.
      • Konantz M.
      • Neuhauss S.C.
      • et al.
      Slc45a2 and V-ATPase are regulators of melanosomal pH homeostasis in zebrafish, providing a mechanism for human pigment evolution and disease.
      ). As well as regulating pigmentation together with SLC45A2, MC1R signaling enhances melanocyte UVR resistance by rapidly up-regulating antioxidant defenses (
      • Song X.
      • Mosby N.
      • Yang J.
      • Xu A.
      • Abdel-Malek Z.
      • Kadekaro A.L.
      alpha-MSH activates immediate defense responses to UV-induced oxidative stress in human melanocytes.
      ) independent of melanin synthesis (
      • Maresca V.
      • Flori E.
      • Bellei B.
      • Aspite N.
      • Kovacs D.
      • Picardo M.
      MC1R stimulation by alpha-MSH induces catalase and promotes its re-distribution to the cell periphery and dendrites.
      ). Melanocytes can act as a sink for reactive oxygen species generated by keratinocytes in vitro (
      • Pelle E.
      • Mammone T.
      • Maes D.
      • Frenkel K.
      Keratinocytes act as a source of reactive oxygen species by transferring hydrogen peroxide to melanocytes.
      ), and MC1R red hair alleles have a reduced capacity to mediate reactive oxygen species protection (
      • Kadekaro A.L.
      • Leachman S.
      • Kavanagh R.J.
      • Swope V.
      • Cassidy P.
      • Supp D.
      • et al.
      Melanocortin 1 receptor genotype: an important determinant of the damage response of melanocytes to ultraviolet radiation.
      ). MC1R-null “red hair” mice are at higher risk of melanoma and oxidative DNA damage in the absence of UVR, suggesting that the increased synthesis of pheomelanin is itself pro-oxidant (
      • Mitra D.
      • Luo X.
      • Morgan A.
      • Wang J.
      • Hoang M.P.
      • Lo J.
      • et al.
      An ultraviolet-radiation-independent pathway to melanoma carcinogenesis in the red hair/fair skin background.
      ). In parallel to animal models, there is some evidence that MC1R red hair alleles are associated with higher risk for keratinocyte skin cancers in a high pigmentation background (
      • Bastiaens M.T.
      • ter Huurne J.A.
      • Kielich C.
      • Gruis N.A.
      • Westendorp R.G.
      • Vermeer B.J.
      • et al.
      Melanocortin-1 receptor gene variants determine the risk of nonmelanoma skin cancer independently of fair skin and red hair.
      ,
      • Box N.F.
      • Duffy D.L.
      • Irving R.E.
      • Russell A.
      • Chen W.
      • Griffyths L.R.
      • et al.
      Melanocortin-1 receptor genotype is a risk factor for basal and squamous cell carcinoma.
      ). Skin cancer risk is likely influenced by the fact that individuals carrying MC1R variant genotypes experience more sunburns (Figure 4).
      • Liu F.
      • Hamer Merel A.
      • Deelen J.
      • Lall Japal S.
      • Jacobs L.
      • van Heemst D.
      • et al.
      The MC1R gene and youthful looks.
      reported that any influence of MC1R on perceived aging is independent of wrinkling; here we show that MC1R SNPs are associated with higher BG6 scores, which correlates with greater photoaging, dermal solar elastosis, and wrinkling. This suggests that these genetic variants are not acting independently of wrinkling.
      We have used VEGAS2 to identify multiple genes near MC1R as potentially associated with BG6 (see Supplementary Table S4). MC1R red hair alleles, and SNPs in LD, are eQTLs for genes other than MC1R (
      • Grundberg E.
      • Small K.S.
      • Hedman A.K.
      • Nica A.C.
      • Buil A.
      • Keildson S.
      • et al.
      Mapping cis- and trans-regulatory effects across multiple tissues in twins.
      ,
      • Westra H.J.
      • Peters M.J.
      • Esko T.
      • Yaghootkar H.
      • Schurmann C.
      • Kettunen J.
      • et al.
      Systematic identification of trans eQTLs as putative drivers of known disease associations.
      ). However, it is also possible that the observation of additional genes via VEGAS2, and MC1R red hair alleles acting as eQTLs, are simply a product of the strong selection on the MC1R locus in Europeans and resulting long-range LD.
      Three SNPs at 6p25.3 are genome-wide significantly associated with BG6, with the lowest P-value at rs12203592 (P = 1.0 × 10–13) in intron four of IRF4 (Table 2). Variants in IRF4 including the T allele of rs12203592 (associated with higher BG6 or more severe skin aging) have been associated with lower pigmentation (but not red hair), poorer tanning ability, greater freckling, and in adults lower mole nevus counts and a higher risk of melanoma and keratinocyte skin cancers (
      • Duffy D.L.
      • Iles M.M.
      • Glass D.
      • Zhu G.
      • Barrett J.H.
      • Hoiom V.
      • et al.
      IRF4 variants have age-specific effects on nevus count and predispose to melanoma.
      ,
      • 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.
      ,
      • Han J.
      • Qureshi A.A.
      • Nan H.
      • Zhang J.
      • Song Y.
      • Guo Q.
      • et al.
      A germline variant in the interferon regulatory factor 4 gene as a novel skin cancer risk locus.
      ,
      • 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.
      ,
      • Zhang M.
      • Song F.
      • Liang L.
      • Nan H.
      • Zhang J.
      • Liu H.
      • et al.
      Genome-wide association studies identify several new loci associated with pigmentation traits and skin cancer risk in European Americans.
      ). The same high photoaging T allele has been shown to reduce activity of an IRF4 enhancer in melanocytes (
      • Praetorius C.
      • Grill C.
      • Stacey S.N.
      • Metcalf A.M.
      • Gorkin D.U.
      • Robinson K.C.
      • et al.
      A polymorphism in IRF4 affects human pigmentation through a tyrosinase-dependent MITF/TFAP2A pathway.
      ). The direction of effect for the T allele of rs12203592 on nevus count and melanoma risk may be reversed in adolescents, with the T allele associated with lower rates of melanoma and higher mole counts (
      • Duffy D.L.
      • Iles M.M.
      • Glass D.
      • Zhu G.
      • Barrett J.H.
      • Hoiom V.
      • et al.
      IRF4 variants have age-specific effects on nevus count and predispose to melanoma.
      ,
      • Kvaskoff M.
      • Whiteman D.C.
      • Zhao Z.Z.
      • Montgomery G.W.
      • Martin N.G.
      • Hayward N.K.
      • et al.
      Polymorphisms in nevus-associated genes MTAP, PLA2G6, and IRF4 and the risk of invasive cutaneous melanoma.
      ). We do not see evidence of a reversal in the BG6 association, the direction of effect of rs12203592 is consistent in adolescents and adults.
      Of particular interest is the overlap in genes associated with two different phenotypes, skin patterning and pigmented spots/lesions (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ). Both studies report the same functional SNP in IRF4, rs12203592. The top BG6 SNP at MC1R is rs4268748, whereas for pigmented lesions it is rs35063026, and rs35063026 is also associated with BG6 (P = 3.5 × 10–10; see Supplementary Table S1). Although the LD between MC1R SNPs rs4268748 and rs35063026 is modest (r2 = 0.22; 1000 Genomes phase 1 version 3 [
      • Abecasis G.R.
      • Altshuler D.
      • Auton A.
      • Brooks L.D.
      • Durbin R.M.
      • et al.
      1000 Genomes Project Consortium
      A map of human genome variation from population-scale sequencing.
      ]), rs35063026 is in LD with the red hair allele rs1805007 (r2 = 0.83), and conditioning on red hair alleles at MC1R abolished the association of rs35063026 with pigmented lesions (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ). The minor allele of rs4268748 is modestly in LD with the minor red hair allele of both rs1805007 (r2 = 0.24) and rs1805008 (r2 = 0.17), likely explaining why it is shows the strongest association with BG6. Conditioning on rs4268748 essentially abolishes the BG6 association of rs1805007 (see Supplementary Table S2; P = 0.03) and rs1805008 (P = 0.4). In summary, in this study the low skin pigmentation minor alleles of functional SNPs in IRF4 and MC1R are associated with higher BG6 scores (i.e., more severe photoaging). The same alleles are associated with a greater percentage of the face covered with pigmented lesions (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ).
      In
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      , the pigmented spots/lesions measured included solar lentigines (SLs) and seborrheic keratoses. SLs are pigmented lesions that contain increased numbers of melanocytes and increased pigmentation of basal keratinocytes, and they are a consequence of chronic UVR damage to the epidermis (
      • Cario-Andre M.
      • Lepreux S.
      • Pain C.
      • Nizard C.
      • Noblesse E.
      • Taieb A.
      Perilesional vs. lesional skin changes in senile lentigo.
      ,
      • Praetorius C.
      • Sturm R.A.
      • Steingrimsson E.
      Sun-induced freckling: ephelides and solar lentigines.
      ). Seborrheic keratoses are common benign skin tumors appearing with age, usually pigmented, and the evidence that they are associated with sun exposure and UVR damage is less conclusive (
      • Hafner C.
      • Vogt T.
      Seborrheic keratosis.
      ). Furthermore, seborrheic keratoses are seen in all skin types, whereas SLs are more common in fairer skin. Skin patterning, as measured by BG6, and epidermal dyspigmentation including SLs, are influenced by both chronological aging and photodamage. Although the sites of the tests in the two studies differ, both are assessing an outcome of sun damage (and natural aging) on a highly sun exposed site, and thus we find overlapping genes that mediate UVR protection (MC1R, IRF4). It is worth noting that, in the case of pigmented lesions, this association is at least partly independent of skin color (
      • Jacobs L.C.
      • Hamer M.A.
      • Gunn D.A.
      • Deelen J.
      • Lall J.S.
      • van Heemst D.
      • et al.
      A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
      ), and the relationship between skin color and SL is complex (
      • Praetorius C.
      • Sturm R.A.
      • Steingrimsson E.
      Sun-induced freckling: ephelides and solar lentigines.
      ), so nonpigmentary roles of these genes may also be involved (e.g.
      • Mitra D.
      • Luo X.
      • Morgan A.
      • Wang J.
      • Hoang M.P.
      • Lo J.
      • et al.
      An ultraviolet-radiation-independent pathway to melanoma carcinogenesis in the red hair/fair skin background.
      ).
      In conclusion, we have identified a (to our knowledge previously unreported) genome-wide significant association with SLC45A2 and confirmed that functional SNPs in IRF4 and MC1R are associated with BG6, a measure of wrinkling and photoaging. At each locus it is the lower-pigmentation alleles that are associated with more severe photoaging. Low-pigmentation alleles of SLC45A2 and MC1R have been extensively associated with increased risk for melanoma and keratinocyte skin cancers (
      • Bastiaens M.T.
      • ter Huurne J.A.
      • Kielich C.
      • Gruis N.A.
      • Westendorp R.G.
      • Vermeer B.J.
      • et al.
      Melanocortin-1 receptor gene variants determine the risk of nonmelanoma skin cancer independently of fair skin and red hair.
      ,
      • Chatzinasiou F.
      • Lill C.M.
      • Kypreou K.
      • Stefanaki I.
      • Nicolaou V.
      • Spyrou G.
      • et al.
      Comprehensive field synopsis and systematic meta-analyses of genetic association studies in cutaneous melanoma.
      ,
      • Han J.
      • Kraft P.
      • Colditz G.A.
      • Wong J.
      • Hunter D.J.
      Melanocortin 1 receptor variants and skin cancer risk.
      ). MC1R red hair alleles, even in the heterozygous state, are associated with a higher burden of somatic mutations in melanoma (
      • Robles-Espinoza C.D.
      • Roberts N.D.
      • Chen S.
      • Leacy F.P.
      • Alexandrov L.B.
      • Pornputtapong N.
      • et al.
      Germline MC1R status influences somatic mutation burden in melanoma.
      ). Photoaging is itself associated with keratinocyte skin cancers. Thus, although photoaging and general changes in appearance as we age are of interest from a quality-of-life and cosmetic point of view, the genes discovered through BG6, which can be measured easily and quantitatively in large-scale populations, may improve our understanding of cancer pathogenesis.

      Materials and Methods

      Samples

      Three independent cohorts were used for this analysis. For the first, adolescent twins and their siblings were recruited in Brisbane, Australia, as a part of the Brisbane Twin Nevus Study (
      • McGregor B.
      • Pfitzner J.
      • Zhu G.
      • Grace M.
      • Eldridge A.
      • Pearson J.
      • et al.
      Genetic and environmental contributions to size, color, shape, and other characteristics of melanocytic naevi in a sample of adolescent twins.
      ,
      • Zhu G.
      • Duffy D.L.
      • Eldridge A.
      • Grace M.
      • Mayne C.
      • O'Gorman L.
      • et al.
      A major quantitative-trait locus for mole density is linked to the familial melanoma gene CDKN2A: a maximum-likelihood combined linkage and association analysis in twins and their sibs.
      ,
      • Zhu G.
      • Montgomery G.W.
      • James M.R.
      • Trent J.M.
      • Hayward N.K.
      • Martin N.G.
      • et al.
      A genome-wide scan for naevus count: linkage to CDKN2A and to other chromosome regions.
      ), totaling 3,759 adolescent twins and their siblings from 1,565 families; these were genotyped in two phases of 2,555 and 1,194 (Table 1). An unrelated set of adult twins comprising a total of 382 people from 197 families was drawn from an alcohol use study (
      • Heath A.C.
      • Bucholz K.K.
      • Madden P.A.
      • Dinwiddie S.H.
      • Slutske W.S.
      • Bierut L.J.
      • et al.
      Genetic and environmental contributions to alcohol dependence risk in a national twin sample: consistency of findings in women and men.
      ,
      • Whitfield J.B.
      • Cullen L.M.
      • Jazwinska E.C.
      • Powell L.W.
      • Heath A.C.
      • Zhu G.
      • et al.
      Effects of HFE C282Y and H63D polymorphisms and polygenic background on iron stores in a large community sample of twins.
      ) collected as a part of the same cohort. The second cohort, the Twins Eyes Study in Tasmania (TEST) (
      • Mackey D.A.
      • Mackinnon J.R.
      • Brown S.A.
      • Kearns L.S.
      • Ruddle J.B.
      • Sanfilippo P.G.
      • et al.
      Twins eye study in Tasmania (TEST): rationale and methodology to recruit and examine twins.
      ), contributed phenotype and genotype data from 136 twins. For the third cohort, 820 Raine participants were drawn from the 22-year follow-up (the Raine Eye Health Study, conducted between March 2012 and July 2014) of the ongoing prospective Western Australian Pregnancy Cohort Study funded by the Raine Medical Research Foundation (
      • Straker L.M.
      • Hall G.L.
      • Mountain J.
      • Howie E.K.
      • White E.
      • McArdle N.
      • et al.
      Rationale, design and methods for the 22 year follow-up of the Western Australian Pregnancy Cohort (Raine) Study.
      ).

      Phenotype collection

      For all participants in all cohorts, imprints of the back of the left hand were made using Affinis light-body silicone elastomer (Coltène AG, Altstätten, Switzerland) (
      • Barnes I.E.
      Techniques for the replication of skin surfaces. A new method.
      ,
      • Sarkany I.
      A method for studying the microtopography of the skin.
      ,
      • Sarkany I.
      • Caron G.A.
      Microtopography of the human skin. Studies with metal-shadowed replicas from plastic impressions.
      ). Imprints for adolescent twins were made close to their 12th, 14th, and 16th birthdays; BG6 scores at age 12 were used preferentially, with missing data supplemented with data from age 14 years over age 16 years, because trait heritability decreases with age. Coincident with the skin measurements, a sun exposure questionnaire was administered for 2,204 of the phase 1 adolescent twins (see Supplementary Methods) (
      • Shekar S.N.
      • Luciano M.
      • Duffy D.L.
      • Martin N.G.
      Genetic and environmental influences on skin pattern deterioration.
      ,
      • Zhu G.
      • Duffy D.L.
      • Eldridge A.
      • Grace M.
      • Mayne C.
      • O'Gorman L.
      • et al.
      A major quantitative-trait locus for mole density is linked to the familial melanoma gene CDKN2A: a maximum-likelihood combined linkage and association analysis in twins and their sibs.
      ). Imprints for nontwin siblings were collected during their twin sibling’s 12th birthday visit. Imprints for the adult sample were collected at a single visit (
      • Shekar S.N.
      • Duffy D.L.
      • Montgomery G.W.
      • Martin N.G.
      A genome scan for epidermal skin pattern in adolescent twins reveals suggestive linkage on 12p13.31.
      , Shekaret al., 2005), as were the Raine and TEST sets.
      A single rater (T. Luong) used a low-power dissection microscope to rate all skin imprints according to the six categories of the Beagley-Gibson rating system, with higher scores reflecting greater disruption of the normal pattern of the lines on the epidermis (
      • Holman C.D.
      • Armstrong B.K.
      • Evans P.R.
      • Lumsden G.J.
      • Dallimore K.J.
      • Meehan C.J.
      • et al.
      Relationship of solar keratosis and history of skin cancer to objective measures of actinic skin damage.
      ). The test-retest correlation for 50 BG6 measurements rescored a second time after 9 months was 0.87 (95% confidence interval = 0.66–0.97) (
      • Shekar S.N.
      • Duffy D.L.
      • Montgomery G.W.
      • Martin N.G.
      A genome scan for epidermal skin pattern in adolescent twins reveals suggestive linkage on 12p13.31.
      ).

      Ethical statement

      Approval to undertake this study was obtained from the Human Research Ethics Committee of the QIMR Berghofer Medical Research Institute. The Raine Study and TEST protocols were approved by the Human Research Ethics Committee of the University of Western Australia and the University of Tasmania, respectively. All studies conformed to the Declaration of Helsinki protocols, with participants giving their written informed consent.

      Genotyping and imputation

      Adolescent samples were genotyped in two phases. Adolescent phase 1 and TEST samples were genotyped using Human610-Quadv1_B arrays (Illumina, San Diego, CA) at deCODE Genetics (Reykjavik, Iceland). For imputation they were merged with a larger set of individuals genotyped on Human610-Quadv1_B, and Human660W-Quad_v1_C arrays (Illumina), and filtered to a common set of 482,586 markers after quality control (QC). Adolescent phase 2 samples were genotyped on HumanOmniExpress-12v1-1_A (Illumina), HumanOmni25M-8v1-1_B (Illumina), and HumanCoreExome 12v1-0_C arrays (Illumina) at the Diamantina Institute, University of Queensland (Queensland, Australia) and filtered to a common set of 233,546 SNPs that passed QC in all sets.
      Adult samples were genotyped on 317K, HumanCNV370-Quadv3_C (Illumina), Human610-Quadv1_B, and Human660W-Quad_v1_C arrays by the Center for Inherited Disease Research (Johns Hopkins University, Baltimore, MD), the University of Helsinki (Helsinki, Finland), and deCODE. For imputation they were merged with a large set of individuals genotyped on the same chips and filtered to a common set of 276,755 SNPs after QC.
      Adolescent, adult, and TEST SNP data were excluded if they failed standard SNP filters in GenomeStudio (Illumina) as recommended by Illumina, and further cleaned out if BeadStudio GenCall (Illumina) was less than 0.7, minor allele frequency was less than 1%, Hardy-Weinberg equilibrium P-value was less than 10–6, or call rate was less than 95%. X chromosome SNPs were filtered if their heterozygosity was greater than 1% in males. Samples with missingness of 5% or greater were dropped, and principal component analysis in Eigenstrat (
      • Price A.L.
      • Patterson N.J.
      • Plenge R.M.
      • Weinblatt M.E.
      • Shadick N.A.
      • Reich D.
      Principal components analysis corrects for stratification in genome-wide association studies.
      ) was used to exclude samples with a principal component 1/principal component 2 value of more than ± 6 standard deviations from 1000 Genomes Projects phase 1 European populations (
      • Abecasis G.R.
      • Altshuler D.
      • Auton A.
      • Brooks L.D.
      • Durbin R.M.
      • et al.
      1000 Genomes Project Consortium
      A map of human genome variation from population-scale sequencing.
      ). Identity by descent estimates in PLINK version 1.07 (
      • Purcell S.
      • Neale B.
      • Todd-Brown K.
      • Thomas L.
      • Ferreira M.A.
      • Bender D.
      • et al.
      PLINK: a tool set for whole-genome association and population-based linkage analysis.
      ) were used to confirm pedigree structure. Because the analyses for these four sets were done in Merlin-offline, which explicitly models familial relationships, samples were not filtered for relatedness.
      Adolescent, adult, and TEST samples were then imputed to the November 23, 2010 release of the 1000 Genomes Phase 1 version 3 (
      • Abecasis G.R.
      • Altshuler D.
      • Auton A.
      • Brooks L.D.
      • Durbin R.M.
      • et al.
      1000 Genomes Project Consortium
      A map of human genome variation from population-scale sequencing.
      ) using minimac (July 17, 2013 version) (
      • Fuchsberger C.
      • Abecasis G.R.
      • Hinds D.A.
      minimac2: faster genotype imputation.
      ,
      • Howie B.
      • Fuchsberger C.
      • Stephens M.
      • Marchini J.
      • Abecasis G.R.
      Fast and accurate genotype imputation in genome-wide association studies through pre-phasing.
      ) with markers with ≤1 copy of the minor alleles removed. Imputation was performed in two phases based on chip overlap: (i) adolescent phase 1 and TEST; (ii) adult and adolescent phase 2. After imputation the four sets were analyzed separately. For genotyped SNPs original genotype calls were used; for imputed SNPs imputation dosage scores were used.
      Raine participants were genotyped on the 660 Quad Array (Illumina) at the Centre for Applied Genomics (Toronto, Ontario, Canada). As part of QC, individuals were excluded if they had identity by descent of π greater than 0.1875 with another participant or greater than 3% missingness. SNPs were filtered out at Hardy-Weinberg equilibrium P-value greater than 5.7 × 10–7, missingness greater than 5%, or a minor allele frequency less than 0.01. Potential population stratification was controlled for via a principal component analysis using Eigenstrat (
      • Price A.L.
      • Patterson N.J.
      • Plenge R.M.
      • Weinblatt M.E.
      • Shadick N.A.
      • Reich D.
      Principal components analysis corrects for stratification in genome-wide association studies.
      ). Imputation of 22 autosomes was performed using MaCH, version 2.3.0 software (
      • Li Y.
      • Willer C.J.
      • Ding J.
      • Scheet P.
      • Abecasis G.R.
      MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes.
      ) and the November 23, 2010, version of the 1000 Genome Project phase 1 European reference panel (
      • Abecasis G.R.
      • Altshuler D.
      • Auton A.
      • Brooks L.D.
      • Durbin R.M.
      • et al.
      1000 Genomes Project Consortium
      A map of human genome variation from population-scale sequencing.
      ).

      Genome-wide association scans

      For the adolescent phases 1 and 2, adult, and TEST samples, GWASs were performed separately using 6-point BG6 scores as a continuous variable and imputed dosage genotypes using the Merlin-offline software (http://www.sph.umich.edu/csg/abecasis/Merlin/) under an additive model that explicitly models relatedness and produces unbiased estimates with standard errors corrected for the sample size inflation due to familial relationship (
      • Chen W.M.
      • Abecasis G.R.
      Family-based association tests for genomewide association scans.
      ). Sex, age, age2, sex × age, and sex × age2 were included as covariates.
      For Raine, a linear regression analysis was performed in R (

      R Core Team. R: a language and environment for statistical computing. Vienna, Austria; 2014.

      ) using the ProbABEL package (
      • Aulchenko Y.S.
      • Struchalin M.V.
      • van Duijn C.M.
      ProbABEL package for genome-wide association analysis of imputed data.
      ). The model was adjusted for age, sex, age2, age × sex, age2 × sex, and the five principal components that accounted for the population stratification.

      Postimputation QC and meta-analysis

      Imputation quality was assessed by the ratio of imputed genotype dosage variance to the variance of the allele frequency under Hardy-Weinberg equilibrium (MaCH rˆ2) (
      • Li Y.
      • Willer C.J.
      • Ding J.
      • Scheet P.
      • Abecasis G.R.
      MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes.
      ). For adolescent phases 1 and 2, TEST, and adult sets, SNPs with a, rˆ2 less than 0.5 or minor allele frequency less than 0.05 in any set were excluded; for Raine, SNPs with imputation quality less than 0.5 were excluded; retained SNPs were meta-analyzed using PLINK version 1.9 (https://www.cog-genomics.org/plink2) (
      • Chang C.C.
      • Chow C.C.
      • Tellier L.C.
      • Vattikuti S.
      • Purcell S.M.
      • Lee J.J.
      Second-generation PLINK: rising to the challenge of larger and richer datasets.
      ), with variant βcoefficients weighted by their inverse variance. Between-study heterogeneity was assessed using the I2 value (
      • DerSimonian R.
      • Laird N.
      Meta-analysis in clinical trials.
      ).

      Conflict of Interest

      The authors state no conflicts of interest.

      Supplementary Material

      References

        • Abecasis G.R.
        • Altshuler D.
        • Auton A.
        • Brooks L.D.
        • Durbin R.M.
        • et al.
        • 1000 Genomes Project Consortium
        A map of human genome variation from population-scale sequencing.
        Nature. 2010; 467: 1061-1073
        • Abecasis G.R.
        • Cherny S.S.
        • Cookson W.O.
        • Cardon L.R.
        Merlin—rapid analysis of dense genetic maps using sparse gene flow trees.
        Nat Genet. 2002; 30: 97-101
        • Aulchenko Y.S.
        • Struchalin M.V.
        • van Duijn C.M.
        ProbABEL package for genome-wide association analysis of imputed data.
        BMC Bioinformatics. 2010; 11: 134
        • Barnes I.E.
        Techniques for the replication of skin surfaces. A new method.
        Br J Dermatol. 1973; 89: 277-283
        • Bastiaens M.T.
        • ter Huurne J.A.
        • Kielich C.
        • Gruis N.A.
        • Westendorp R.G.
        • Vermeer B.J.
        • et al.
        Melanocortin-1 receptor gene variants determine the risk of nonmelanoma skin cancer independently of fair skin and red hair.
        Am J Hum Genet. 2001; 68: 884-894
        • Battistutta D.
        • Pandeya N.
        • Strutton G.M.
        • Fourtanier A.
        • Tison S.
        • Green A.C.
        Skin surface topography grading is a valid measure of skin photoaging.
        Photodermatol Photoimmunol Photomed. 2006; 22: 39-45
        • Beagley J.
        • Gibson I.M.
        Changes in skin condition in relation to degree of exposure to ultraviolet light.
        Western Australia Institute of Technology School of Biology, Perth1980
        • Bilac C.
        • Sahin M.T.
        • Ozturkcan S.
        Chronic actinic damage of facial skin.
        Clin Dermatol. 2014; 32: 752-762
        • Box N.F.
        • Duffy D.L.
        • Irving R.E.
        • Russell A.
        • Chen W.
        • Griffyths L.R.
        • et al.
        Melanocortin-1 receptor genotype is a risk factor for basal and squamous cell carcinoma.
        J Invest Dermatol. 2001; 116: 224-229
        • Cario-Andre M.
        • Lepreux S.
        • Pain C.
        • Nizard C.
        • Noblesse E.
        • Taieb A.
        Perilesional vs. lesional skin changes in senile lentigo.
        J Cutan Pathol. 2004; 31: 441-447
        • Chang A.L.
        • Atzmon G.
        • Bergman A.
        • Brugmann S.
        • Atwood S.X.
        • Chang H.Y.
        • et al.
        Identification of genes promoting skin youthfulness by genome-wide association study.
        J Invest Dermatol. 2014; 134: 651-657
        • Chang C.C.
        • Chow C.C.
        • Tellier L.C.
        • Vattikuti S.
        • Purcell S.M.
        • Lee J.J.
        Second-generation PLINK: rising to the challenge of larger and richer datasets.
        GigaScience. 2015; 4: 7
        • Chatzinasiou F.
        • Lill C.M.
        • Kypreou K.
        • Stefanaki I.
        • Nicolaou V.
        • Spyrou G.
        • et al.
        Comprehensive field synopsis and systematic meta-analyses of genetic association studies in cutaneous melanoma.
        J Natl Cancer Inst. 2011; 103: 1227-1235
        • Chen W.M.
        • Abecasis G.R.
        Family-based association tests for genomewide association scans.
        Am J Hum Genet. 2007; 81: 913-926
        • Costin G.E.
        • Valencia J.C.
        • Vieira W.D.
        • Lamoreux M.L.
        • Hearing V.J.
        Tyrosinase processing and intracellular trafficking is disrupted in mouse primary melanocytes carrying the underwhite (uw) mutation. A model for oculocutaneous albinism (OCA) type 4.
        J Cell Sci. 2003; 116: 3203-3212
        • DerSimonian R.
        • Laird N.
        Meta-analysis in clinical trials.
        Control Clin Trials. 1986; 7: 177-188
        • Dooley C.M.
        • Schwarz H.
        • Mueller K.P.
        • Mongera A.
        • Konantz M.
        • Neuhauss S.C.
        • et al.
        Slc45a2 and V-ATPase are regulators of melanosomal pH homeostasis in zebrafish, providing a mechanism for human pigment evolution and disease.
        Pigment Cell Melanoma Res. 2013; 26: 205-217
        • Duffy D.L.
        • Box N.F.
        • Chen W.
        • Palmer J.S.
        • Montgomery G.W.
        • James M.R.
        • et al.
        Interactive effects of MC1R and OCA2 on melanoma risk phenotypes.
        Hum Mol Genet. 2004; 13: 447-461
        • Duffy D.L.
        • Iles M.M.
        • Glass D.
        • Zhu G.
        • Barrett J.H.
        • Hoiom V.
        • et al.
        IRF4 variants have age-specific effects on nevus count and predispose to melanoma.
        Am J Hum Genet. 2010; 87: 6-16
        • Fritschi L.
        • Battistutta D.
        • Strutton G.M.
        • Green A.
        A non-invasive measure of photoageing.
        Int J Epidemiol. 1995; 24: 150-154
        • Fuchsberger C.
        • Abecasis G.R.
        • Hinds D.A.
        minimac2: faster genotype imputation.
        Bioinformatics. 2015; 31: 782-784
        • Gandini S.
        • Sera F.
        • Cattaruzza M.S.
        • Pasquini P.
        • Zanetti R.
        • Masini C.
        • et al.
        Meta-analysis of risk factors for cutaneous melanoma: III. Family history, actinic damage and phenotypic factors.
        Eur J Cancer. 2005; 41: 2040-2059
        • Graf J.
        • Hodgson R.
        • van Daal A.
        Single nucleotide polymorphisms in the MATP gene are associated with normal human pigmentation variation.
        Hum Mutat. 2005; 25: 278-284
        • Green A.C.
        Premature ageing of the skin in a Queensland population.
        Med J Aust. 1991; 155 (7–8): 473-474
        • Grundberg E.
        • Small K.S.
        • Hedman A.K.
        • Nica A.C.
        • Buil A.
        • Keildson S.
        • et al.
        Mapping cis- and trans-regulatory effects across multiple tissues in twins.
        Nat Genet. 2012; 44: 1084-1089
        • GTEx Consortium
        The Genotype-Tissue Expression (GTEx) project.
        Nat Genet. 2013; 45: 580-585
        • Hafner C.
        • Vogt T.
        Seborrheic keratosis.
        J Dtsch Dermatol Ges. 2008; 6: 664-677
        • Han J.
        • Kraft P.
        • Colditz G.A.
        • Wong J.
        • Hunter D.J.
        Melanocortin 1 receptor variants and skin cancer risk.
        Int J Cancer. 2006; 119: 1976-1984
        • 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.
        PLoS Genet. 2008; 4: e1000074
        • Han J.
        • Qureshi A.A.
        • Nan H.
        • Zhang J.
        • Song Y.
        • Guo Q.
        • et al.
        A germline variant in the interferon regulatory factor 4 gene as a novel skin cancer risk locus.
        Cancer Res. 2011; 71: 1533-1539
        • Heath A.C.
        • Bucholz K.K.
        • Madden P.A.
        • Dinwiddie S.H.
        • Slutske W.S.
        • Bierut L.J.
        • et al.
        Genetic and environmental contributions to alcohol dependence risk in a national twin sample: consistency of findings in women and men.
        Psychol Med. 1997; 27: 1381-1396
        • Holman C.D.
        • Armstrong B.K.
        • Evans P.R.
        • Lumsden G.J.
        • Dallimore K.J.
        • Meehan C.J.
        • et al.
        Relationship of solar keratosis and history of skin cancer to objective measures of actinic skin damage.
        Br J Dermatol. 1984; 110: 129-138
        • Howie B.
        • Fuchsberger C.
        • Stephens M.
        • Marchini J.
        • Abecasis G.R.
        Fast and accurate genotype imputation in genome-wide association studies through pre-phasing.
        Nat Genet. 2012; 44: 955-959
        • Hughes M.C.
        • Strutton G.M.
        • Fourtanier A.
        • Green A.C.
        Validation of skin surface microtopography as a measure of skin photoaging in a subtropical population aged 40 and over.
        Photodermatol Photoimmunol Photomed. 2012; 28: 153-158
        • Hughes M.C.
        • Williams G.M.
        • Baker P.
        • Green A.C.
        Sunscreen and prevention of skin aging: a randomized trial.
        Ann Intern Med. 2013; 158: 781-790
        • Jacobs L.C.
        • Hamer M.A.
        • Gunn D.A.
        • Deelen J.
        • Lall J.S.
        • van Heemst D.
        • et al.
        A genome-wide association study identifies the skin color genes IRF4, MC1R, ASIP, and BNC2 influencing facial pigmented spots.
        J Invest Dermatol. 2015; 135: 1735-1742
        • Jacobs L.C.
        • Liu F.
        • Pardo L.M.
        • Hofman A.
        • Uitterlinden A.G.
        • Kayser M.
        • et al.
        IRF4, MC1R and TYR genes are risk factors for actinic keratosis independent of skin color.
        Hum Mol Genet. 2015; 24: 3296-3303
        • Kadekaro A.L.
        • Leachman S.
        • Kavanagh R.J.
        • Swope V.
        • Cassidy P.
        • Supp D.
        • et al.
        Melanocortin 1 receptor genotype: an important determinant of the damage response of melanocytes to ultraviolet radiation.
        FASEB J. 2010; 24: 3850-3860
        • Kvaskoff M.
        • Whiteman D.C.
        • Zhao Z.Z.
        • Montgomery G.W.
        • Martin N.G.
        • Hayward N.K.
        • et al.
        Polymorphisms in nevus-associated genes MTAP, PLA2G6, and IRF4 and the risk of invasive cutaneous melanoma.
        Twin Res Hum Genet. 2011; 14: 422-432
        • Laville V.
        • Le Clerc S.
        • Ezzedine K.
        • Jdid R.
        • Taing L.
        • Labib T.
        • et al.
        A genome-wide association study in Caucasian women suggests the involvement of HLA genes in the severity of facial solar lentigines.
        Pigment Cell Melanoma Res. 2016; 29: 550-558
        • Lavker R.M.
        • Kwong F.
        • Kligman A.M.
        Changes in skin surface patterns with age.
        J Gerontol. 1980; 35: 348-354
        • Le Clerc S.
        • Taing L.
        • Ezzedine K.
        • Latreille J.
        • Delaneau O.
        • Labib T.
        • et al.
        A genome-wide association study in Caucasian women points out a putative role of the STXBP5L gene in facial photoaging.
        J Invest Dermatol. 2013; 133: 929-935
        • Li Y.
        • Willer C.J.
        • Ding J.
        • Scheet P.
        • Abecasis G.R.
        MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes.
        Genet Epidemiol. 2010; 34: 816-834
        • Liu F.
        • Hamer Merel A.
        • Deelen J.
        • Lall Japal S.
        • Jacobs L.
        • van Heemst D.
        • et al.
        The MC1R gene and youthful looks.
        Curr Biol. 2016; 26: 1-8
        • Mackey D.A.
        • Mackinnon J.R.
        • Brown S.A.
        • Kearns L.S.
        • Ruddle J.B.
        • Sanfilippo P.G.
        • et al.
        Twins eye study in Tasmania (TEST): rationale and methodology to recruit and examine twins.
        Twin Res Hum Genet. 2009; 12: 441-454
        • Maresca V.
        • Flori E.
        • Bellei B.
        • Aspite N.
        • Kovacs D.
        • Picardo M.
        MC1R stimulation by alpha-MSH induces catalase and promotes its re-distribution to the cell periphery and dendrites.
        Pigment Cell Melanoma Res. 2010; 23: 263-275
        • McGregor B.
        • Pfitzner J.
        • Zhu G.
        • Grace M.
        • Eldridge A.
        • Pearson J.
        • et al.
        Genetic and environmental contributions to size, color, shape, and other characteristics of melanocytic naevi in a sample of adolescent twins.
        Genet Epidemiol. 1999; 16: 40-53
        • Mishra A.
        • Macgregor S.
        VEGAS2: software for more flexible gene-based testing.
        Twin Res Hum Genet. 2015; 18: 86-91
        • Mitra D.
        • Luo X.
        • Morgan A.
        • Wang J.
        • Hoang M.P.
        • Lo J.
        • et al.
        An ultraviolet-radiation-independent pathway to melanoma carcinogenesis in the red hair/fair skin background.
        Nature. 2012; 491: 449-453
        • Nica A.C.
        • Parts L.
        • Glass D.
        • Nisbet J.
        • Barrett A.
        • Sekowska M.
        • et al.
        The architecture of gene regulatory variation across multiple human tissues: the MuTHER study.
        PLoS Genet. 2011; 7: e1002003
        • Pelle E.
        • Mammone T.
        • Maes D.
        • Frenkel K.
        Keratinocytes act as a source of reactive oxygen species by transferring hydrogen peroxide to melanocytes.
        J Invest Dermatol. 2005; 124: 793-797
        • Praetorius C.
        • Grill C.
        • Stacey S.N.
        • Metcalf A.M.
        • Gorkin D.U.
        • Robinson K.C.
        • et al.
        A polymorphism in IRF4 affects human pigmentation through a tyrosinase-dependent MITF/TFAP2A pathway.
        Cell. 2013; 155: 1022-1033
        • Praetorius C.
        • Sturm R.A.
        • Steingrimsson E.
        Sun-induced freckling: ephelides and solar lentigines.
        Pigment Cell Melanoma Res. 2014; 27: 339-350
        • Price A.L.
        • Patterson N.J.
        • Plenge R.M.
        • Weinblatt M.E.
        • Shadick N.A.
        • Reich D.
        Principal components analysis corrects for stratification in genome-wide association studies.
        Nat Genet. 2006; 38: 904-909
        • Purcell S.
        • Neale B.
        • Todd-Brown K.
        • Thomas L.
        • Ferreira M.A.
        • Bender D.
        • et al.
        PLINK: a tool set for whole-genome association and population-based linkage analysis.
        Am J Hum Genet. 2007; 81: 559-575
      1. R Core Team. R: a language and environment for statistical computing. Vienna, Austria; 2014.

        • Robles-Espinoza C.D.
        • Roberts N.D.
        • Chen S.
        • Leacy F.P.
        • Alexandrov L.B.
        • Pornputtapong N.
        • et al.
        Germline MC1R status influences somatic mutation burden in melanoma.
        Nat Comm. 2016; 7: 12064
        • Sarkany I.
        A method for studying the microtopography of the skin.
        Br J Dermatol. 1962; 74: 254-259
        • Sarkany I.
        • Caron G.A.
        Microtopography of the human skin. Studies with metal-shadowed replicas from plastic impressions.
        J Anat. 1965; 99: 359-364
        • Seddon J.M.
        • Egan K.M.
        • Zhang Y.
        • Gelles E.J.
        • Glynn R.J.
        • Tucker C.A.
        • et al.
        Evaluation of skin microtopography as a measure of ultraviolet exposure.
        Invest Ophthalmol Vis Sci. 1992; 33: 1903-1908
        • Shekar S.N.
        • Duffy D.L.
        • Montgomery G.W.
        • Martin N.G.
        A genome scan for epidermal skin pattern in adolescent twins reveals suggestive linkage on 12p13.31.
        J Invest Dermatol. 2006; 126: 277-282
        • Shekar S.N.
        • Luciano M.
        • Duffy D.L.
        • Martin N.G.
        Genetic and environmental influences on skin pattern deterioration.
        J Invest Dermatol. 2005; 125: 1119-1129
        • Song X.
        • Mosby N.
        • Yang J.
        • Xu A.
        • Abdel-Malek Z.
        • Kadekaro A.L.
        alpha-MSH activates immediate defense responses to UV-induced oxidative stress in human melanocytes.
        Pigment Cell Melanoma Res. 2009; 22: 809-818
        • Straker L.M.
        • Hall G.L.
        • Mountain J.
        • Howie E.K.
        • White E.
        • McArdle N.
        • et al.
        Rationale, design and methods for the 22 year follow-up of the Western Australian Pregnancy Cohort (Raine) Study.
        BMC Public Health. 2015; 15: 663
        • 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.
        Nat Genet. 2007; 39: 1443-1452
        • Westra H.J.
        • Peters M.J.
        • Esko T.
        • Yaghootkar H.
        • Schurmann C.
        • Kettunen J.
        • et al.
        Systematic identification of trans eQTLs as putative drivers of known disease associations.
        Nat Genet. 2013; 45: 1238-1243
        • Whitfield J.B.
        • Cullen L.M.
        • Jazwinska E.C.
        • Powell L.W.
        • Heath A.C.
        • Zhu G.
        • et al.
        Effects of HFE C282Y and H63D polymorphisms and polygenic background on iron stores in a large community sample of twins.
        Am J Hum Genet. 2000; 66: 1246-1258
        • Zhang M.
        • Song F.
        • Liang L.
        • Nan H.
        • Zhang J.
        • Liu H.
        • et al.
        Genome-wide association studies identify several new loci associated with pigmentation traits and skin cancer risk in European Americans.
        Hum Mol Genet. 2013; 22: 2948-2959
        • Zhu G.
        • Duffy D.L.
        • Eldridge A.
        • Grace M.
        • Mayne C.
        • O'Gorman L.
        • et al.
        A major quantitative-trait locus for mole density is linked to the familial melanoma gene CDKN2A: a maximum-likelihood combined linkage and association analysis in twins and their sibs.
        Am J Hum Genet. 1999; 65: 483-492
        • Zhu G.
        • Montgomery G.W.
        • James M.R.
        • Trent J.M.
        • Hayward N.K.
        • Martin N.G.
        • et al.
        A genome-wide scan for naevus count: linkage to CDKN2A and to other chromosome regions.
        Eur J Hum Genet. 2007; 15: 94-102