Investigating DNA methylation as a potential mediator between pigmentation genes, pigmentary traits and skin cancer.

Pigmentation characteristics are well-known risk factors for skin cancer. Polymorphisms in pigmentation genes have been associated with these traits and with the risk of malignancy. However, the functional relationship between genetic variation and disease is still unclear. This study aims to assess whether pigmentation SNPs are associated with pigmentary traits and skin cancer via DNA methylation (DNAm). Using a meta-GWAS of whole-blood DNAm from 36 European cohorts (N = 27,750; the Genetics of DNA Methylation Consortium, GoDMC), we found that 19 out of 27 SNPs in 10 pigmentation genes were associated with 391 DNAm sites across 30 genomic regions. We examined the effect of 25 selected DNAm sites on pigmentation traits, sun exposure phenotypes and skin cancer and on gene expression in whole blood. We uncovered an association of DNAm site cg07402062 with red hair in the Avon Longitudinal Study of Parents and Children (ALSPAC). We also found that the expression of ASIP and CDK10 was associated with hair colour, melanoma and basal cell carcinoma. Our results indicate that DNAm and expression of pigmentation genes may play a role as potential mediators of the relationship between genetic variants, pigmentation phenotypes and skin cancer and thus deserve further scrutiny.

Several genetic variants in pigmentation‐related genes have been associated with skin cancer, particularly in populations of European ancestry. DNA methylation was investigated as a potential link between pigmentation SNPs, pigmentation phenotypes and skin cancer. DNA methylation sites near MC1R and ASIP were strongly associated with pigmentation SNPs, gene expression, hair colour, melanoma and basal cell carcinoma in an analysis that included data from the Genetics of DNA Methylation consortium, the Avon Longitudinal Study of Parents and Children and the UK Biobank. Further evaluation of the role played by DNA methylation in these relationships is warranted.

INTRODUCTION

The incidence of skin cancer, comprising malignant melanoma, basal cell carcinoma (BCC) and squamous cell carcinoma (SCC), has increased rapidly in the past decades (Cahoon et al., 2015; Hu et al., 2009; Rees et al., 2014). Melanoma is the most aggressive of skin cancers although it has a low incidence, while non‐melanoma skin cancer (NMSC) shows high incidence yet considerably lower mortality rates compared to melanoma (Apalla et al., 2017). To date, recognized risk factors for melanoma and NMSC include fair skin, light‐coloured eyes, red hair, freckles and melanocytic naevi as well as genetic variants, some of which underlie these skin pigmentation and sun sensitivity phenotypes (Gordon, 2013). Even though there are missense and nonsense polymorphisms in pigmentation genes strongly associated with skin cancer, particularly in the gene MC1R (Nasti & Timares, 2015), it is not fully understood how non‐coding variation in these genes relates to malignancy.

DNA methylation (DNAm) is an epigenetic modification with a prospective role in cancer aetiology. The association of global whole‐blood DNA hypomethylation with cancer is well known and has also been described for melanoma (Cappetta et al., 2015; Shen et al., 2017). In addition, recent studies have shown the presence of hypomethylation in skin biopsy samples that were exposed to sunlight or artificial ultraviolet radiation (Grönniger et al., 2010; Holzscheck et al., 2020; Vandiver et al., 2015). Inter‐individual DNAm variation at specific sites, measured in peripheral blood, has been uncovered as a predictor of a number of complex trait risk factors as well as all‐cause mortality (McCartney et al., 2018). DNAm has the potential to be used as a molecular biomarker to diagnose disease or to assess prognosis in those affected by disease and is a promising candidate leading towards the realization of personalized medicine (Nikolouzakis et al., 2020).

This study investigated whether DNAm plays a role in the relationship between genetic variants, pigmentation‐related skin cancer risk factors and skin cancer. First, we examined the effect of genetic variation on DNAm levels in peripheral blood across 10 regions robustly associated with pigmentation traits and skin cancer, in 36 cohorts of European descent. We then analysed the association of pigmentation SNP‐associated DNAm sites with sun exposure and pigmentation phenotypes in participants of the Avon Longitudinal Study of Parents and Children (ALSPAC). Finally, using summary data‐based Mendelian randomization (SMR), we explored whether SNP‐related DNAm was likely to causally underlie the expression of genes associated with pigmentation phenotypes and skin cancer. The purpose of this work was to carry out an integrative analysis to evaluate the joint contribution of genetic and epigenetic risk factors to skin cancer susceptibility. As this is the first instance of such an analysis involving genes associated with pigmentation traits and skin cancer, the clinical relevance of our findings will depend on subsequent replication and further investigation of the genes and pathways that we described here.

MATERIALS AND METHODS

The different stages of this study are depicted in Figure 1.

FIGURE 1

Stages of analysis implemented and datasets used in this study. From the 391 unique DNA methylation (DNAm) sites associated with pigmentation SNPs, 25 were selected for further analysis as depicted. GoDMC, Genetics of DNA methylation consortium; ALSPAC, Avon Longitudinal Study of Parents and Children; CIT, causal inference test; SMR, summary data‐based Mendelian randomization; eQTL browser, blood eQTL browser; CAGE, Cap Analysis of Gene Expression; GTEx, Genotype‐Tissue Expression consortium; UKB, UK Biobank

Identification of DNAm sites associated with pigmentation SNPs

The genetics of DNA methylation consortium (GoDMC) was created to study the genetic basis of DNAm variation and bring together resources and researchers with expertise in the epigenetics field (www.godmc.org.uk). One of its aims was to carry out a meta‐GWAS of DNAm, as measured on Illumina 450k or EPIC Beadchips. Results from a meta‐GWAS involving 36 European cohorts (N = 27,750) were used here. We provide a brief description of the analysis in the Methods S1.

Using GoDMC data, we searched for DNAm sites that were strongly associated (p < 1 × 10‐5) with well‐known pigmentation‐related SNPs previously identified via genome‐wide association studies (GWAS) or candidate gene studies (reviewed in reference Pavan & Sturm, 2019), located within the genes (or their surrounding regions) ASIP (rs1015362, rs4911414, rs619865), BNC2 (rs2153271), IRF4 (rs12203592, rs12210050), HERC2 (rs12913832), MC1R (rs1110400, rs11547464, rs11648785, rs1805005, rs1805006, rs1805007, rs1805008, rs1805009, rs2228479, rs258322, rs4785763, rs885479), OCA2 (rs1800401, rs1800407), SLC24A4 (rs12896399), SLC24A5 (rs1426654), SLC45A2 (rs16891982, rs28777) and TYR (rs1042602, rs1393350).

We uncovered 874 strong SNP‐DNAm associations across 30 genomic regions, which included 391 unique DNAm sites. No data were available for SNPs in SLC24A5 and SLC45A2, neither for SNPs rs1110400, rs11547464, rs1805006 and rs1805009 in MC1R and rs1393350 in TYR. Some of these polymorphisms have low minor allele frequencies (MAF) and GoDMC only included SNPs with MAF > 1%. The other SNPs may not have been part of the genotyping platforms, successfully imputed, or strongly associated with any DNAm site (p > 1 × 10‐5).

Variance in DNAm explained by SNPs was estimated as 2*β2*MAF*(1 − MAF), where β is the effect size and MAF is the minor allele frequency.

Linkage disequilibrium (LD) r 2 values were obtained with LDlink (https://ldlink.nci.nih.gov/).

DNA methylation and pigmentation/sun exposure phenotypes in the Avon Longitudinal Study of Parents and Children (ALSPAC)

Cohort description and information regarding the collection of DNAm data in ALSPAC can be found in the Methods S1.

Regression analysis

The pigmentation/sun exposure phenotypes evaluated included the following: skin reflectance, freckles, moles, sunburning, tanning ability, hair colour and eye colour. We examined DNAm at the nearest time point with respect to when the phenotypes were measured. We assessed the association of cord blood DNAm and pigmentation traits measured before age 7 [i.e. skin reflectance (49 months), freckles (49 and/or 61 months), red hair (15 and 54 months), eye colour (54 months), tanning ability (69 months)]; the association of DNAm in childhood (~7 years old) with sunburning from birth to age 12, and total number of moles at 15 years old; and the association of DNAm in adolescence (15–17 years old) with total number of moles at 15 years old, red hair at 18 years old and tanning ability at ~25 years old. Description of most of these phenotypes has been provided in an earlier work (Bonilla et al., 2014). Data on hair colour and tanning ability in young adults were collected in a recent ALSPAC questionnaire (Life@25+), using the same scales employed in past questionnaires (Bonilla et al., 2014).

We tested the association of SNP‐associated DNAm sites with pigmentation and sun exposure phenotypes using t tests, one‐way ANOVA tests, linear and logistic regression models. Regressions were adjusted for age, sex and the first 10 genetic principal components to account for population stratification. Pairwise correlation between DNAm sites was examined in ALSPAC children at age 7. All analyses were carried out with the statistical package Stata v15.

Mediation analysis

The causal inference test (CIT) approach (Millstein et al., 2009), implemented in the R package ‘cit’, was employed to investigate the causal direction between pigmentation SNPs, DNAm and red hair colour. Regressions were adjusted for sex, age and the top 10 genetic principal components.

Summary data‐based Mendelian randomization (SMR)

In order to investigate the potential for DNAm at selected sites to be causally related to gene expression and phenotypes such as pigmentation characteristics and skin cancer, we carried out a Mendelian randomization analysis with colocalization, implemented in the summary data‐based Mendelian randomization (SMR) method (Zhu et al., 2016). We used the platform Complex Trait Genetics Virtual Lab (CTG‐VL, https://genoma.io/; Cuellar‐Partida et al., 2019) to run the statistical package SMR (Lloyd‐Jones et al., 2017; Qi et al., 2018; Zhu et al., 2016). This method uses summary data from GWAS, mQTL or eQTL studies to distinguish causal or pleiotropic associations between DNAm and gene expression, or between either of these and a phenotype, from a situation where the traits are caused by different genetic variants that are in strong LD (see figure 1b of Zhu et al., 2016). If the former is true, that is a single genetic variant underlies DNAm and the other tested phenotype, the traits are said to be colocalized. The different SMR analyses run are described below. Unless otherwise reported, the SMR analyses p‐value thresholds were as follows: p SMR < 5 × 10‐7/5x10‐6 and p HEIDI ≥ 0.05.

SMR of whole‐blood DNA methylation and gene expression

Peripheral blood DNAm summary data were extracted from McRae et al. (McRae et al., 2018) while summary data on gene expression in blood were obtained from three different sources: the Genotype‐Tissue Expression consortium (GTEx, https://gtexportal.org/home/), Cap Analysis of Gene Expression (CAGE; Lloyd‐Jones et al., 2017), and the blood eQTL browser (https://genenetwork.nl/bloodeqtlbrowser/; Westra et al., 2013). All expression datasets used the hg19 genome assembly.

SMR of gene expression in skin, pigmentation/sun exposure traits and skin cancer

We ran SMR to assess the association of gene expression in sun‐exposed and sun‐unexposed skin, obtained from GTEx, with pigmentary and skin cancer phenotypes.

Data were obtained by CTG‐VL from the UK Biobank and trait definitions can be found in the biobank website (https://www.ukbiobank.ac.uk/). Pigmentation characteristics analysed were skin colour (UKB ID#1717, N = 356,530), ease of skin tanning (#1727, N = 353,697), childhood sunburn occasions (#1737, N = 269,734), and black (N = 15,809) and blonde (N = 41,178) hair colour (#1747, N = 360,270). We also considered diagnosed malignant melanoma (#ICD10:C43, N = 1,672, total N = 361,194), self‐reported malignant melanoma (N = 2,898), self‐reported basal cell carcinoma (N = 3,441), self‐reported squamous cell carcinoma (N = 449; all #20001, total N = 361,141) and having melanocytic naevi (#ICD10:D22, N = 3,501 cases, total N = 361,194).

SMR of whole‐blood DNA methylation, pigmentation/sun exposure traits and skin cancer

We also investigated the potential colocalization of genetic variants underlying DNAm with pigmentation traits and skin cancer. DNAm utilized in this analysis was measured in blood (McRae et al., 2018).

Heritability of the DNA methylation sites associated with pigmentation SNPs

We checked the heritability of DNAm sites using the resource provided by the Complex Disease Epigenetics Group (www.epigenomicslab.com/online‐data‐resources/; Hannon, Knox, et al., 2018) that employs twin data to report the variance in whole‐blood DNAm explained by an additive genetic component, a shared environmental component and a unique environmental component.

RESULTS

Identification of DNA methylation sites associated with pigmentation SNPs

We summarized available functional information on the analysed SNPs in Table 1.

TABLE 1

Pigmentation SNPs investigated in GoDMC in association with DNA methylation (DNAm) changes

SNP	Chromosome	Position	Gene a	Effect allele b	Other allele	EAF c	DNAm sites associated	Other information d
rs12203592	6p25.3	396321	IRF4	T	C	0.188	3	eQTL (sun‐exposed skin)
rs12210050	6p25.3	475489	IRF4	T	C	0.173	7
rs2153271	9p22.2	16864521	BNC2	T	C	0.592	4	eQTL (whole blood)
rs1042602	11q14.3	89011046	TYR	A	C	0.262	5	NP_000363.1:p.Ser192Tyr
rs12896399	14q32.12	92773663	SLC24A4	T	G	0.468	6
rs1800407	15q13.1	28230318	OCA2	T	C	0.081	1	NP_000266.2:p.Arg419Gln/ NP_001287913.1:p.Arg395Gln
rs1800401	15q13.1	28260053	OCA2	G	A	0.946	2	NP_000266.2:p.Arg305Trp
rs12913832	15q13.1	28365618	HERC2	G	A	0.744	3	eQTL (whole blood)
rs258322	16q24.3	89755903	CDK10	A	G	0.118	98	eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs1805005	16q24.3	89985844	MC1R	T	G	0.121	26	rhc (NP_002377.4:p.Val60Leu)
rs2228479	16q24.3	89985940	MC1R	A	G	0.089	163	rhc (NP_002377.4:p.Val92Met)/eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs1805007	16q24.3	89986117	MC1R	T	C	0.086	97	RHC (NP_002377.4:p.Arg151Cys)/ eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs1805008	16q24.3	89986144	MC1R	T	C	0.083	97	RHC (NP_002377.4:p.Arg160Trp)/ eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs885479	16q24.3	89986154	MC1R	A	G	0.105	110	rhc (NP_002377.4:p.Arg163Gln)/ eQTL (sun‐exposed skin)
rs4785763	16q24.3	90066936	AFG3L1P	A	C	0.332	104	eQTL/ eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs11648785	16q24.3	90084561	DBNDD1	C	T	0.689	67
rs4911414	20q11.22	32729444	ASIP	T	G	0.342	24
rs1015362	20q11.22	32738612	ASIP	C	T	0.727	27	eQTL (sun‐exposed and sun‐unexposed skin and whole blood)
rs619865	20q11.22	33867697	EIF6	A	G	0.109	29
rs16891982	5p13.2	33951693	SLC45A2	G	C	0.938	n/a	NP_001012527.2:p.Leu374Phe
rs28777	5p13.2	33958959	SLC45A2	A	C	0.956	n/a
rs1393350	11q14.3	89011046	TYR	A	G	0.244	n/a
rs1426654	15q21.1	48426484	SLC24A5	A	G	0.997	n/a	NP_995322.1:p.Thr111Ala
rs1805006	16q24.3	89985918	MC1R	A	C	0.010	n/a	RHC (NP_002377.4:p.Asp84Glu)
rs11547464	16q24.3	89986091	MC1R	A	G	0.009	n/a	RHC (NP_002377.4:p.Arg142His)
rs1110400	16q24.3	89986130	MC1R	C	T	0.008	n/a	RHC (NP_002377.4:p.Ile155Thr)
rs1805009	16q24.3	89986546	MC1R	C	G	0.008	n/a	RHC (NP_002377.4:p.Asp294His)

In this study, for simplicity, SNPs on chromosome 16q24.3 are considered part of the MC1R genetic region, and SNPs on chromosome 20q11.22 are considered part of the ASIP genetic region.

The effect allele is the allele associated with a lighter skin, hair and eye colour, and skin cancer susceptibility.

Effect allele frequency from GoDMC for SNPs with associated DNAm sites. Effect allele frequency from 1,000 Genomes for SNPs not available in GoDMC.

RHC = high penetrance red hair colour variant, rhc = low penetrance red hair colour variant. eQTL information obtained from the GTEx Portal v8 release.

Pigmentation SNPs that were strongly associated with DNAm sites in GoDMC are shown in Table S1.

We selected DNAm sites with the most reliable effects for further analysis as follows: sites associated with at least six of the MC1R region SNPs, three of the ASIP region SNPs or showing the strongest association with pigmentation SNPs in the other genes tested. Additionally, chosen DNAm sites had to display consistent associations with the allele that increased fair pigmentation for a minimum of five out of eight SNPs in MC1R and two out of three SNPs in ASIP (i.e. the allele increasing fair pigmentation had to always increase or always decrease DNAm at the same site for a majority of the SNPs tested in that gene). Finally, DNAm sites had to be associated with SNPs in ALSPAC participants consistently across the time points where DNAm was measured (i.e. birth, childhood, adolescence, pregnancy, middle age). Based on these criteria, we followed up a total of 25 DNAm sites: 17 in MC1R, 3 in ASIP and one each in IRF4, BNC2, TYR, SLC24A4 and HERC2. DNAm sites that were selected for further analysis are shown in Table 2.

TABLE 2

Pigmentation SNP‐associated DNA methylation (DNAm) sites selected for follow‐up

DNAm site	Gene region	Gene	chr	chr position	# Associated SNPs	Associated in ALSPAC at all time points?
cg26114043	MC1R a	INTU	4	128544375	7	Y
cg09806625	IRF4	EXOC2	6	611523	1	Y
cg03291755	BNC2	BNC2	9	16868891	1	Y
cg05041596	TYR	NAALAD2	11	89867385	1	Y
cg04136915	SLC24A4	intergenic	14	92721383	1	Y
cg14091419	HERC2	HERC2	15	28356810	1	N b
cg05504729	MC1R	ANKRD11	16	89399490	6	Y
cg08845973	MC1R	SPG7	16	89592071	6	Y
cg09738481	MC1R	CPNE7	16	89653835	6	Y
cg01097406	MC1R	intergenic	16	89675127	7	Y
cg04240660	MC1R	CHMP1A	16	89714849	6	Y
cg03605463	MC1R	LOC284241	16	89740564	6	Y
cg04287289	MC1R	FANCA	16	89883240	6	Y
cg26513180	MC1R	FANCA	16	89883248	7	Y
cg07402062	MC1R	SPIRE2	16	89894098	7	Y
cg21285383	MC1R	SPIRE2	16	89894308	7	Y
cg03388025	MC1R	SPIRE2	16	89894329	7	Y
cg01768446	MC1R	MC1R	16	89982419	8	Y
cg01023759	MC1R	MC1R	16	89983766	8	Y
cg01883217	MC1R	DEF8	16	90015832	8	Y
cg05185784	MC1R	DEF8	16	90016020	7	Y
cg15936718	MC1R	GAS8	16	90088801	7	Y
cg01901788	ASIP	MAP1LC3A	20	33145848	3	Y
cg08999081	ASIP	PIGU	20	33150536	3	Y
cg21132536	ASIP	ACSS2	20	33465180	3	Y

cg26114043 is located on chromosome 4 but it is strongly associated with SNPs on chromosome 16 (MC1R region).

cg14091419 was associated with HERC2 SNP rs12913832 in pregnancy and middle‐age only.

MC1R

In this gene region, we obtained results for red hair colour loss‐of‐function variants (RHC) R151C (rs1805007) and R160W (rs1805008), and missense variants (rhc) V60L (rs1805005), V92M (rs2228479) and R163Q (rs885479) and for non‐coding polymorphisms found in GWAS to be associated with pigmentation and melanoma risk (rs11648785, rs258322 and rs4785763).

MC1R non‐coding polymorphisms and functional variants exhibited the largest number of associations with DNAm sites of all the pigmentation genes examined, ranging from 26 (rs1805005) to 163 (rs2228479; Table 1). Three hundred and twelve unique sites were strongly associated with at least one of the MC1R SNPs tested. Of these, 131 were associated with just one genetic variant. There was overlap between DNAm sites associated with different MC1R SNPs (Table S2). Twenty‐five DNAm sites were associated with a minimum of 6 of the SNPs with results, while 5 were associated with all 8 of them. However, not all of these DNAm sites showed a consistent direction of association with the alleles that increase the odds of having a fair skin phenotype and being affected by melanoma. For instance, full consistency was achieved for MC1R DNAm site cg01097406 (i.e. all alleles that are associated with fair skin decrease DNAm at this position) and cg07130392, both associated with 7 of the 8 SNPs tested, and for cg08845973 and cg09738481, associated with 6 SNPs.

Interestingly, DNAm sites associated with both RHC variants, which are located in the first exon of the gene SPIRE2 (i.e. cg03388025, cg04287289, cg07402062, cg21285383, cg26513180), exhibited effects of similar magnitude and direction (Figure S1). However, for these same DNAm sites, the association with rhc variants rs2228479 and rs885479 occurred in the opposite direction to that of the RHCs (Table S1).

Although most associations uncovered were between SNPs and DNAm sites in cis, MC1R showed several strong associations in trans (Figure 2). Besides associations with more distant DNAm sites on chromosome 16 (>1 Mb), we identified associations with sites on chromosomes 1, 2, 3, 4, 5, 8, 10, 11, 12 and 13 (Table S1, Figure 2). One of the trans sites most frequently associated with MC1R SNPs (7 of them, largely in low LD) was cg26114043 located in the gene inturned planar cell polarity protein (INTU) on chromosome 4 position 128,544,375. This gene is related to the processes of keratinocyte differentiation and hair follicle morphogenesis (https://www.ncbi.nlm.nih.gov/gene/27152).

FIGURE 2

Circos plot showing the relationship between SNPs in pigmentation genes and DNA methylation (DNAm) sites in cis and trans. The external track (light blue) depicts the effect of each pigmentation SNP on DNAm (negative in red, positive in blue). The internal track (pink) depicts the R 2 (variance explained) values. Information extracted from Table S1

The maximum variance in DNAm explained by a pigmentation SNP was 24.4% shown by the pair rs2228479‐cg09569215. Amongst the RHC, the highest R2 was found for rs1805007, which explained 17.8% of variation at cg08854185 (Table S1).

Other pigmentation genes

The number of SNP–DNAm site associations for the other pigmentation genes examined ranged from 1 (rs1800407 in OCA2) to 29 (rs619865 in ASIP; Table 1). There were 127 unique DNAm sites associated with at least one SNP in these 7 genes, while 102 DNAm sites were associated with only one polymorphism.

Five different DNAm sites in the ASIP region were associated with three SNPs each, whereas one site was associated with 3 SNPs in the OCA2/HERC2 region (Table S2). There was no DNAm site overlap between SNPs rs12210050 and rs12203592 in IRF4. The only fully consistent association between SNPs and DNAm sites was with cg01901788, located near ASIP (i.e. for all 3 SNPs the alleles associated with a fair skin phenotype increased DNAm at this site).

As expected, most associations were in cis, although there were three associations in trans in ASIP (distant chromosome 20) and one in BNC2 (chromosome 1; Table S1, Figure 2). Maximum variance in DNAm explained by SNPs was: 8.9% (ASIP rs4911414‐cg08999081), 8.6% (BNC2 rs2153271‐cg03291755), 5.3% (OCA2 rs1800407‐cg20906524), 3.0% (TYR rs1042602‐cg25941151), 2.7% (IRF4 rs12210050‐cg26187313), 1.7% (HERC2 rs12913832‐cg20906524) and 1.3% (SLC24A4 rs12896399‐cg16993582; Table S1).

DNA methylation sites and pigmentation/sun exposure phenotypes in ALSPAC

Pairwise correlation coefficients (Pearson's) and p‐values for DNAm sites measured in ALSPAC children at age 7 are shown in Table S3a–c. There was moderate to strong correlation between DNAm sites in the region that extends from cg04287289 to cg03388025 in MC1R (Table S3b), and some moderate correlation was detected between DNAm sites in unlinked genetic regions (Table S3a).

When we investigated the association of the selected DNAm sites with pigmentation and sun exposure traits in ALSPAC children, we found the strongest associations with red hair colour. The most robust association with this phenotype was that of cg07402062, where redhead children showed lower levels of DNAm at this site than non‐redhead children (Table S4). This effect was seen with cord blood DNAm and red hair at 15 and 54 months (after correction for multiple testing), as well as with DNAm in adolescence and red hair at 18 years old (before correction for multiple testing). Associations with the other traits tested were either weaker or absent.

The causal inference test (CIT) showed inconclusive results for SNPs rs1805007, rs1805008, rs258322, and rs4785763 (Table S5). There was evidence for a mediating effect of DNAm at cg07402062 between these pigmentation SNPs and red hair colour, as well as for red hair as a mediator of the relationship between the pigmentation SNPs and DNAm at cg07402062. On the other hand, for variants rs2228479 and rs885479, we found evidence of a causal effect of DNAm at cg07402062 on red hair colour, but not the other way around. Finally, rs11648785 appeared to influence both traits independently.

Summary data‐based Mendelian randomization (SMR) of whole‐blood DNA methylation and gene expression

Out of the 25 DNAm sites under analysis, four (cg01023759, cg09738481, cg14091419 and cg26114043) were not present in any of the expression datasets used. Since the data sources only included variation in cis, cg26114043, located on chromosome 4, was left out. The 21 remaining sites showed associations with gene expression to varying degrees which were, for the most part, consistent between datasets. Of note were the associations of correlated DNAm sites cg03388025, cg04287289, cg07402062, cg21285383 and cg26513180 with the expression of SPIRE2, cg08845973 with SPG7, cg01901788 with NCOA6 and cg21132536 with GGT7. Results of this analysis, including total numbers of associations per DNAm site and dataset, are shown in Table S6.

Summary data‐based Mendelian randomization (SMR) of gene expression in skin, pigmentation/sun exposure traits and skin cancer

We found no evidence of a pleiotropic or causal association between gene expression at 16q24.3, 20q11.22, or any of the other pigmentation regions included in this study, and ease of skin tanning, skin colour, childhood sunburn occasions, melanocytic naevi and SCC, in sun‐exposed and sun‐unexposed skin.

There was evidence of an association with the same genetic variant between black hair colour and the expression of genes ASIP and FAM83C (20q11.22) in sun‐exposed skin and of CDK10 (16q24.3), ASIP, NCOA6 (20q11.22) and EXOC2 (6p25.3) in sun‐unexposed skin, whereas blonde hair colour was associated with the expression of ASIP and FAM83C in sun‐exposed skin and with NCOA6 and EXOC2 in sun‐unexposed skin. In addition, ASIP gene expression was associated with self‐reported and diagnosed malignant melanoma, and BCC in both skin tissues, and CDK10 was associated with self‐reported and diagnosed malignant melanoma, and BCC in unexposed skin (Table S7).

When the HEIDI p‐value threshold was Bonferroni‐corrected for multiple testing (instead of being set as ≥ 0.05), additional genes showed potential pleiotropic/causal associations with the traits of interest, including CPNE7, DBNDD1, DEF8, MC1R, SNAI3‐AS1, SPATA2L, URAHP in 16q24.3, and EIF6, GGT7, UQCC1 in 20q22.11.

Additional summary data‐based Mendelian randomization (SMR) of whole‐blood DNA methylation and gene expression

Because the DNAm sites we selected were associated with the expression of just 5 of the genes identified in the previous step (i.e. DBNDD1, DEF8, GGT7, NCOA6, SPATA2L), we searched for additional DNAm sites that showed an association with the expression of the genes reported above. For CDK10, there was no common colocalized DNAm site across all the datasets, rather associations were evident with several DNAm sites in the 16q24.3 region. For ASIP and FAM83C, there were no data available or no evidence of association with DNAm sites, and EXOC2 showed two overlapping associated DNAm sites between the CAGE and GTEx datasets. Results for these and other genes are shown in Table S8.

Summary data‐based Mendelian randomization (SMR) of whole‐blood DNA methylation, pigmentation/sun exposure traits and skin cancer

We uncovered 696 DNAm sites associated and colocalized with pigmentary and skin cancer traits (Table S9), but of these, only a reduced number (n = 18) were associated with gene expression colocalizing with these phenotypes (see Table S8). From this set of 18 DNAm sites, 9 were identified in our original analysis in GoDMC but only one (cg01901788) was on the list we followed up (the others were associated with three SNPs or less, Table S2). SNPs associated with these “mediating” 9 DNAm sites were as follows: rs11648785, rs1805008 (R160W), rs2228479 (V92M) and rs4785763 in the MC1R region; and rs1015362, rs4911414 and rs619865 in the ASIP region (Table 3).

TABLE 3

DNA methylation sites identified in GoDMC that colocalized with the expression of genes which, in turn, colocalized with pigmentation/sun exposure/skin cancer phenotypes. These DNAm sites also colocalized with pigmentation/sun exposure/skin cancer phenotypes in an independent test

DNAm sites	Chromosome	Associated pigmentation SNPs	Colocalized genes	Colocalized phenotypes a
cg01901788	20q11.22	rs1015362, rs4911414, rs619865	NCOA6	skin colour, ease of skin tanning, blonde hair, black hair
cg04270048	20q11.22	rs4911414, rs619865	NCOA6	skin colour, ease of skin tanning, childhood sunburn, blonde hair, black hair
cg16810054	20q11.22	rs619865	NCOA6	blonde hair, black hair
cg04378830	16q24.3	rs2228479	DBNDD1, MC1R	blonde hair, black hair, self‐reported BCC
cg04752812	16q24.3	rs2228479	DBNDD1, DEF8	blonde hair, skin colour, ease of skin tanning, self‐reported BCC
cg09569215	16q24.3	rs11648785, rs2228479, rs4785763	MC1R	blonde hair, black hair
cg10062109	16q24.3	rs1805008, rs4785763	CDK10	blonde hair, black hair, self‐reported melanoma, skin colour, ease of skin tanning, childhood sunburn, diagnosed melanoma, self‐reported BCC
cg11215013	16q24.3	rs2228479	SPATA2L	skin colour, black hair, self‐reported melanoma
cg26598918	16q24.3	rs1805008, rs4785763	CDK10	blonde hair, black hair, skin colour, ease of skin tanning, childhood sunburn, diagnosed melanoma, self‐reported BCC

In bold: DNAm site followed up in this study.

Phenotypes in red colocalized with the DNAm site shown in the first column, phenotypes in blue colocalized with the expression of the gene/s listed in the fourth column, and phenotypes in green colocalized with both, the DNAm site and gene expression. See Tables S6–S9 for the colocalization results.

Heritability

Nine (38%) of 24 DNAm sites followed up showed a heritability value of 80% or more (additive genetic contribution), almost all of these are located in MC1R. DNAm site cg01023759 was unavailable in this dataset. Median heritability was 63% (IQR 31%, 85%; Table S10). The five correlated DNAm sites associated with SPIRE2 expression and also with SNPs rs1805007, rs1805008, rs2228479 and rs885479 showed a heritability higher than 63%.

DISCUSSION

In this study, we reported the existence of strong associations between pigmentation/sun exposure/skin cancer‐associated SNPs and DNAm variation at specific sites, mostly in cis, with a few cases in trans. This appears to be especially marked in the MC1R gene region, which could be related to the high CpG content present in this region (Martinez‐Cadenas et al., 2013). The SNP‐DNAm site associations were observed for functional variants (RHC and rhc) as well as for non‐coding polymorphisms within MC1R.

For a set of selected DNAm sites (those associated with several pigmentation SNPs in the same direction), we assessed their effect on pigmentation/sun exposure phenotypes measured in ALSPAC participants. Although sample sizes were modest, we detected a strong inverse association of red hair with cg07402062, which is located in the gene SPIRE2 (MC1R region). Remarkably, a recent GWAS of hair colour has identified the strongest genetic signal of association with red hair (versus brown/black hair) at SNP rs34357723 (T allele associated with red hair; LD with rs1805007, r 2 = .27), a signal that remained relevant even after adjustment for MC1R coding variants (Morgan et al., 2018). Rs34357723 is an intronic variant of SPIRE2, and GoDMC data indicate that it is also a strong mQTL for cg07402062, explaining ~16% of DNAm at this site (T allele beta = −0.82, SE = 0.01). The CIT analysis suggested that a mediation effect of cg07402062 DNAm between rs34357723 and red hair colour is compatible with ALSPAC data (see Table S5), although replication of these findings is needed.

Galanter and colleagues reported that DNAm at cg07402062 and cg21285383 differed between ethnic groups and that these differences could be accounted for by differences in genetic ancestry (Galanter et al., 2017), which would also be consistent with hair colour variation between groups.

The high heritability detected for these two DNAm sites, along with cg03388025, cg04287289 and cg26513180, and their association with RHC and rhc variants suggests a close relationship between SPIRE2 DNAm status and the MC1R gene structure.

SPIRE2 was recently shown to participate in myosin Va‐dependent transport of melanosomes in melanocytes via the generation of actin tracks (Alzahofi et al., 2018). However, despite SPIRE2 expression colocalizing with DNAm at cg07402062 and correlated sites, we detected no evidence of colocalization with hair colour GWAS signals, and thus, we cannot yet conclude that this gene causally underlies the phenotype.

Using summary data and two‐sample Mendelian randomization implemented in SMR, we determined that the selected DNAm sites affected expression of a number of genes in cis, but there was limited overlap between this set of genes and those that showed an association of expression with pigmentation traits and skin cancer. Given that the DNAm expression analysis was carried out using whole‐blood data, and the expression‐phenotype analysis used skin data, the differences observed may be due to this discrepancy.

Previous reports on DNAm in pigmentation/skin cancer‐associated regions include a study by Roos and colleagues (Roos et al., 2017) who showed that DNAm in healthy human skin was associated in cis with SNPs underlying melanoma risk discovered by GWAS, including two polymorphisms in MC1R (rs258322, rs4785763, which we tested) and one in the ASIP region (rs910873, not considered by us but in LD with rs619865). SNPs in SLC45A2 (rs35390) and TYR (rs1393350, rs1847134) were also investigated but yielded no evidence of association with DNAm, a finding similar to ours. Roos et al. also found that none of the DNAm sites associated with these SNPs was associated with naevus count in their epigenome‐wide association study (EWAS). Likewise, in our SMR analysis presence/absence of melanocytic naevi failed to colocalize with any DNAm site.

Another study by Budden and Bowden described the regulation of MC1R expression by a CpG island overlapping a potential enhancer in the MC1R gene (Budden & Bowden, 2018). These authors found that this CpG island is unmethylated in melanocytes but highly methylated in other skin cell types. In addition, most of the melanoma cell lines and tumours examined in their study were methylated in this region. Further analysis indicated that differential methylation of the CpG island in tumours mainly affected the expression of MC1R, and less so of TUBB3. Budden and Bowden also reported that there was weak evidence of increasing CpG island methylation levels with increasing numbers of red hair colour alleles (R) in melanoma tumours from The Cancer Genome Atlas (TCGA) dataset (Budden & Bowden, 2018).

While seven of the 12 MC1R SNPs that we evaluated were located within the CpG island identified by Budden and Bowden, they were strongly associated with DNAm sites located nearby, although outside the island/enhancer region, making direct comparisons of results difficult.

Lastly, Martinez‐Cadenas and colleagues (Martinez‐Cadenas et al., 2013) postulated that the considerable nucleotide diversity evident at MC1R was probably the result of a high mutation rate that occurred as a consequence of the elevated CpG concentration present in the region. Chromosome 16q is also one of the richest regions in gene and GC content, which has been found to be correlated (Zeng et al., 2018). In this study, polymorphisms in MC1R were associated with DNAm variation at substantially more sites than is the case in other pigmentation genes, with lower GC content, tested (for instance, MC1R region SNPs are associated with 95 DNAm sites each on average, whereas ASIP region SNPs were associated with an average of 27 DNAm sites each, see Table 1).

From an evolutionary point of view, it is interesting to consider the potential for epigenetic modifications in pigmentation genes as a mechanism of adaptation to variable levels of ultraviolet radiation exposure, providing phenotypic flexibility in a changing environment, similar to what has been suggested for adaptation to high altitude (Julian, 2017). Since there is evidence that pigmentation genes, such as MC1R and ASIP, have been shaped by natural selection (Martinez‐Cadenas et al., 2013; Norton et al., 2007; Quillen et al., 2019), polymorphic variants in them may act, in part, by affecting how they are epigenetically regulated to optimize a response to an environment with reduced sunlight.

Limitations

Our analysis was conducted on DNA derived from whole‐blood rather than skin tissue, and any pigmentation SNP‐DNAm site associations may not be as strong or present in skin. We used heritability estimates provided by Hannon, Gorrie‐Stone, et al. (2018) to assess the likelihood of the DNAm sites examined to vary in a similar way in different tissues. These authors reported a low contribution on average of additive genetic effects on DNAm across the genome (~16%) and a higher heritability of the sites that were correlated across tissues (median 71%). Since most of the DNAm sites (16/24) we analysed had heritability values > 50%, it is possible that the associations that we detected with pigmentation SNPs are also there in skin. Moreover, for two MC1R polymorphisms, Roos et al. (2017) already described an association in healthy skin with some of the same DNAm sites identified in our study. Although not examining mQTLs, Yuan et al. (2019) also reported an association of DNAm sites in 16q24.3 within the genes SPG7, DPEP1, FANCA and SPIRE2, which included cg04287289 and cg26513180, with ethnicity in placental tissue. However, the small overlap between the genes whose whole‐blood expression was affected by the DNAm sites we tested and those whose expression in skin was associated with pigmentation and skin cancer phenotypes could be suggesting that there are in fact differences between tissues.

Unfortunately, mQTL data on populations of non‐European ancestry are limited. In the only study involving one of these populations that we identified, Pierce and colleagues reported an association of 5 of our pigmentation DNAm sites with SNPs in cis in the population of Bangladesh (Pierce et al., 2018). While this is an indication that similar genetic–epigenetic relationships are likely to occur in other groups, knowledge of the role of these chromosomal regions on pigmentation in non‐White populations is still scarce (Pavan & Sturm, 2019).

Each of the methods we used has their own particular limitations, so it is advisable to be cautious in the interpretation of our findings.

The CIT method requires that the analysis is carried out in one sample where the SNP, exposure and outcome have all been measured, and therefore, statistical power is constrained by the size of this sample. It may also be affected by collider bias and phenotype measurement error. In order to protect against collider bias, it is essential to adjust for all confounders of the exposure‐outcome relationship. While we adjusted the regression models for age, sex and genetic principal components, unknown confounders may be still playing a role (Hemani et al., 2017).

The SMR method was conceived to try to tell apart pleiotropy from LD in the event of an association between gene expression and a complex trait, and thus to interpret results under a model of pleiotropy rather than causality, even though it could be consistent with the latter as well (Zhu et al., 2016). To tease out a causal effect, additional analyses are necessary, such as two‐step Mendelian randomization (Relton & Davey Smith, 2012) using independent (in cis or trans) mQTLs. However, independent instrumental variables for the DNAm sites in this study have been hard to find since there is usually considerable LD with the original SNPs. As a consequence, SMR results are mostly intended to be used to highlight genes for follow‐up studies and not to establish unequivocal causal relationships (Zhu et al., 2016).

In summary, our results suggest that DNAm may lie on the causal pathway leading from pigmentation‐associated SNPs to pigmentation/sun exposure traits and eventually skin cancer, possibly via regulation of gene expression, although more evidence to support these findings is needed. This study puts forward a set of genes that could be prioritized in upcoming investigations (e.g. CDK10 or SPIRE2) to delve deeper into the role played by the variability in DNAm as a potentially causal risk factor for skin cancer. We have reported for the first time the existence of DNAm sites strongly influenced by genetic variants across pigmentation genes and demonstrated the importance of integrating different sources of variation and of uncovering the interaction of markers within and between chromosomal regions (for instance, on chromosome 16 and chromosome 4) to obtain a better picture of pathways leading to disease. However, there is still some way to go before human genetic and epigenetic information can be introduced fully into the dermatology clinic for the provision of personalized medicine.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

Figure S1

Click here for additional data file.

Tables S1–S10

Click here for additional data file.

Methods S1

Click here for additional data file.