A large-scale genome-wide association meta-analysis for nevus count provides direct insights into the genetics of melanoma

Domain: biology
Published: 2026-03-10
Editorial outcome: minor_revision

Abstract
A greater understanding of the biology of nevi will provide insights into the etiology of melanoma. Our
large-scale meta-analysis of 14 nevus genome-wide association study (GWAS) include 85,965 individuals
of European ancestry. We identify 29 nevus-associated loci (p < 5x10-8), of which 24 have not been
previously reported in a GWAS conducted for nevus count alone. We further identify 255 candidate genes
for nevus loci, including SIKE1 which is involved in immune response regulation. This is of interest
because immune response regulation influences the formation of nevi and melanoma susceptibility. Gene-
set enrichment analyses prioritise immune response-related pathways and cancers that do not have a
pigmentation component (e.g. breast, prostate, and glioma). This suggests that the biology underlying
nevus count captures risk pathways beyond pigmentation that are relevant to melanoma. In sex-specific
analyses we observe higher total-body nevus count in females than in males, however the genetic
architecture is largely shared (genetic correlation = 0.863, 95% CI = 0.453 – 1.273), indicating the
difference may be influenced by environmental and behavioural factors rather than genetics. A nevus
polygenic risk score explains 5% of the variance in nevus count, indicating its potential to enhance
melanoma risk prediction.

Introduction
Melanocytic nevi are small, pigmented benign skin tumours composed of melanocytes. Melanoma is the
deadliest form of skin cancer and arises from the malignant transformation of melanocytes. Melanoma
can develop from a pre-existing nevus or appear as a new lesion on the skin1,2. A high melanocytic nevus
count is the strongest single risk factor for melanoma3. Therefore, understanding the biology of nevi may
also provide insights into the etiology of melanoma.

Nevus count is a highly heritable trait with heritability estimates from twins ranging from 60% to 70% 4.
Previous genome-wide association studies (GWAS) for nevus count (n > 52,000) have identified five
significant loci5. These loci were also associated with melanoma risk5, implying both traits have a shared
genetic risk component. Although combined analyses of both melanoma and nevus count have identified
additional loci that may be pleiotropically associated with both traits5,6, given the high heritability of nevus
count, many causal genetic variants remain to be discovered.

Nevus counts show sex-specific patterns by body site, with males exhibiting higher counts on the head,
neck and trunk, and females often showing higher counts on the limbs, particularly the lower limbs7.
These patterns mirror sex-specific differences observed in the body site distribution of melanoma7, and
both behavioural and genetic factors may contribute to them. Therefore, evaluating potential phenotypic
and genetic differences in nevus count between the sexes provides biological insight into the determinants
of nevus formation, and helps clarify whether sex-specific pathways contribute to downstream melanoma
susceptibility.

In addition to providing insights into cutaneous melanoma, nevus count genetics can improve our
understanding of rarer forms of melanoma. Uveal melanoma is a rare but deadly form of melanoma
originating from ocular melanocytes. Atypical nevi and iris nevi are risk factors for uveal melanoma8–10,
suggesting that the findings from a cutaneous nevus GWAS can be applied to risk factors for uveal
melanoma.

To identify further genetic variants associated with nevus count, we meta-analysed 14 GWASs of nevus
count, comprising over 85,000 individuals, including two previously unpublished datasets. We also report
on the application of post-GWAS functional analyses to identify those genes most strongly involved in
nevus biology and to identify the pathways underlying these genes. We also conducted sex-specific
GWAS meta-analyses of nevus count vary by sex and explored the degree of causal association between
melanoma and nevus count through Mendelian randomisation analysis. Finally, we generated polygenic
scores for nevus count and tested their ability to predict nevus count in adults and to predict iris nevus
count.

Results
Nevus GWAS meta-analysis
We initially estimated the degree to which genotyped SNPs could capture variation in nevus count within
the QSkin I cohort (n = 15,346). The SNP-based heritability of nevus count estimated using GCTA-
GREML was 0.22 (SE = 0.03).

We then performed a large-scale GWAS meta-analysis of nevus count with a total sample size of 85,965
participants. The meta-analysis included variants on both autosomes and the X chromosome (see
Methods); however, no genome-wide significant signals were observed on the X chromosome (Figure 1).
Therefore, subsequent analyses focused exclusively on autosomal findings. The genomic inflation factor
\lambda and LDSC intercept for the GWAS meta-analysis suggested polygenic inheritance (\lambda = 1.068, intercept
= 1.014, ratio = 0.119). In the nevus meta-analysis, 41 independent genome-wide significant SNPs were
assigned to 29 genomic loci (Figure 1). Of these 29 loci, 24 have not previously been reported in a GWAS
conducted solely for nevus count (Table 1). Nineteen of these 24 loci have been previously identified as
associated with both nevus count and melanoma in joint analyses5,6. Among the remaining five loci after
19/24, the regions near MAP3K1 on chromosome 5 and around TPCN2 on chromosome 11 have also been
associated with melanoma and hair colour6. LDSC regression analysis estimated the proportion of the
variance in nevus count captured by the GWAS SNPs in the current GWAS meta-analysis at 0.06 (SE =
0.01); while this is expected to be an underestimate relative to GCTA-GREML as it uses fewer SNPs and
summary data, this estimate of variance is two times higher than that of the previous nevus GWAS meta-
analysis5.

Pleiotropy of nevus loci
To determine whether the significant nevus loci were pleiotropic with melanoma, a pairwise GWAS was
conducted by integrating our results with a previously published melanoma GWAS 6. This approach
assigns a posterior probability of association (PPA), which sums to 1, to four models; models 1 and 2
suggest the gene region is only associated with trait 1 or 2, respectively; model 3 suggests pleiotropy
and model 4 tests an independent association11 (see Methods). Twenty-eight of 29 nevus loci exhibited strong
evidence of pleiotropic regions for both melanoma and nevus count (model 3 PPA > 0.95) (Supplementary
Data 1). Locus, near TPCN2 on chromosome 11, had the highest PPA for model 4 (independence; PPA
0.76), but as no model reached the threshold (> 0.95), we are unable to reach strong conclusions about the
role of this locus.

Functional mapping
We subsequently used FUMA (Functional Mapping and Annotation of GWASs) to identify potential
functional links between nevus-associated variants and candidate genes at genome-wide significant loci
(Methods). Using three mapping approaches (positional, eQTL and chromatin interaction mapping) 248
genes were prioritised across the nevus count loci (Supplementary Data 2). These prioritised genes
included previously reported nevus (e.g. KITLG, MTAP/CDKN2A and IRF4) or melanoma-related genes
(e.g. MC1R, TMEM38B and DCST2)5,6,12. Of 248 genes, 26 overlapped in all three mapping approaches,
and those genes mapped on 14 of 29 nevus-associated genomic loci (Supplementary Data 3). Interestingly,
of these 26 genes, we identified candidate genes for 2/3 novel genomic regions that have not previously
been reported in a nevus GWAS, nevus/melanoma combined analysis or a nevus/pigmentation combined
analysis, marked by rs7537232 (BCAS2 and SIKE1) and rs12902005 (AKAP13 and KLHL25) as indicated
in Table 1. In addition, PTTG1 was identified as the candidate gene for the remaining nevus-associated
novel locus on chromosome 5 through the positional and eQTL mapping.

Gene-based analyses
MAGMA gene-based analysis evaluates candidate genes for association with a trait by combining the
local association test statistics for SNPs within a gene region. This approach can identify additional
significant genetic loci that may not reach genome-wide significance in single SNP tests. MAGMA
identified 33 genome-wide significant genes (P < 2.64x10-6) at 18 genomic loci; three were in addition to
those found in our GWAS meta-analysis (Figure 2 and Supplementary Data 4). Of these 33 genes, 28
overlapped with those prioritised by functional mapping approaches and included well-known nevus-,
melanoma-, or pigmentation-associated genes providing further evidence of their shared genetics. The
three additional loci identified by MAGMA included a novel locus not previously identified for nevus
count or melanoma on chromosome 22 (C220rf26), a locus previously found in nevus/melanoma
combined analysis on chromosome 10 (FAM208B)5, and a locus previously reported for melanoma on
chromosome 3 (MYNN)12. Interestingly, through MAGMA gene analysis we found NRAS as a significant
candidate gene for nevi at the nevus risk locus marked by rs7537232.

Transcriptome-wide association study (TWAS)
To expand on our gene-mapping work, we carried out TWAS to nominate plausible functional genes
linked to nevus count by testing for significant correlation with gene expression. We analysed gene
expression data from three skin-related tissue-based cis-eQTL datasets (sun-exposed, not-sun-exposed,
and fibroblasts) available in GTEx (see Methods), melanocyte-specific eQTL data6 and an additional well-
powered but potentially less tissue-relevant whole blood dataset from the eQTLGen consortium. After
applying multiple-testing corrections followed by colocalisation analysis, the five TWAS analyses
identified a range of promising candidate genes. Across all TWAS analyses, candidate genes were
identified at 12 genome-wide significant nevus loci and two additional loci on chromosome 10 (LRMDA
and GDI2) (Supplementary Data 5.1 and 5.2). LRDMA was significant in non-sun-exposed skin tissue and
showed increased predicted gene expression associated with higher nevus count (Z-score = 5.061). GDI2
was significant across sun-exposed skin, fibroblasts and whole blood tissues with uniformly negative Z-
scores (Z-score = -5.157, -5.505 and -4.950, respectively), indicating a consistent inverse association
across tissues. Further, GDI2 was replicated in MAGMA gene-based analysis (Supplementary Data 4).
The consistency of signals across tissues, together with its replication across methods, provides strong
evidence for implicating GDI2 in nevus count.

Pathway analyses
In the MAGMA gene-set analysis, three gene sets involved in the regulation of the microtubule
cytoskeleton during radial glia-guided migration in the cerebral cortex, development of astrocytoma in
patients with Neurofibromatosis Type 1 (NF1) and regulation of the catalytic activity (i.e. any process
that modulates the activity of enzyme), were found to be significant after Bonferroni correction
(Supplementary Data 6).

In addition to using the MAGMA data, we conducted gene-set enrichment tests using the genes prioritised
by any of the three gene mapping approaches in FUMA. Notably, prioritised genes were significantly
enriched in gene sets related to cancers (e.g. breast, prostate and glioma cancer) that do not have a
pigmentation component. As expected, gene-sets previously reported for nevus, melanoma and
pigmentation were also enriched (Supplementary Data 7). Fifteen genes (including NRAS, THBS3,
IFNA16, CCND1 and KITLG) were enriched in the critical cell survival, growth and apoptosis P13K/AKT
signaling pathway (adjusted p = 1.68×10-5). In addition, we found significant enrichments in several
immune response-related pathways. The top two pathways were the type 1 interferon receptor binding
pathway and RIG-I-like (Retinoic acid-inducible gene I) receptor signaling pathway (Supplementary Data
7). Further, 28 prioritised genes were significantly overrepresented in the gene set involved in response
to the organic cyclic compounds. Of transcriptional factor target gene sets, four gene sets showed
significant enrichment signals, of which the most significant was the target gene set regulated by
transcription factor PRDM4 (PR domain zinc finger protein 4), which is involved in cell proliferation and
tumourigenesis13,14.

Sex-specific meta-analysis GWAS for nevus count
We found differences between males and females for nevus count in QSkin I, which consists of 6,769
males and 8,577 females. To test the difference, we fitted an ordinal logistic regression model for nevus
count on a 4-point scale (none, few, moderate, more) and sex after adjusting for age, and ancestry (first
10 PCs). Sex was associated with the nevus count, and females had significantly higher odds of being in
a higher nevus count category compared with males (OR = 1.294, 95% CI = 1.216 - 1.377).

Then we performed a sex-stratified GWAS meta-analysis (N for men = 11,963; N for females = 21,766)
integrating QSkin I, QSkin II, and AGDS. Sex-specific GWAS identified nine genome-wide significant
loci in females and four in males, with three loci overlapping between the GWASs (Supplementary Fig.
20). We identified two loci (rs10873172[SYNE2] on chromosome 14 and rs35096708[FANCA] on
chromosome 16) with significantly different effect sizes by sex (P < 0.005, Bonferroni-corrected p value
for 10 tests) as indicated in Supplementary Data 8. However, genetic correlation between males and
females for nevus count calculated from LDSC regression was not significantly different to 1 (r = 0.863,
95% CI = 0.453 – 1.273, p = 3.86×10-5).

Causal association with melanoma
Mendelian randomisation (MR) showed that genetically predicted higher nevus count was significantly
associated with a higher risk of melanoma (Inverse variance weighted (IVW) OR = 4.212, 95% CI =
3.469-5.114). Sensitivity analysis yielded consistent effect estimates with overlapping confidence
intervals, reinforcing the robustness of the association (Supplementary Table 1, Supplementary Fig. 21).
The MR-Egger intercept provided only weak evidence of directional pleiotropy (intercept = 0.019, se =
0.009, p = 0.047), and its small magnitude suggests any bias is likely minimal. Although modest
heterogeneity was detected (Cochran’s Q = 35.48, p = 0.035), the agreement of the causal estimate across
methods indicates that the observed heterogeneity did substantially affect causal inference.

Polygenic risk score analysis for nevi
We built a nevus PRS (Polygenic Risk Score) using SBayesRC applied to a meta-analysis of 13 out of the
14 nevus GWASs (excluding QSkin II) (see Methods). We tested the nevus PRS on the left-out QSkin II
data (n = 6,307). The model including nevus PRS explained 5% (R2 = 0.05) of the variation and a beta
coefficient per (standard deviation) SD of the nevus PRS was 0.19 (SE = 0.0.01, p < 4.94e-324),
substantially exceeding the performance of the previous clumping PRS15 (R² = 0.0036, beta = 0.051, SE
= 0.01). Further, we applied the same nevus PRS to the Kidskin Young Adult Myopia Study (KYAMS)
cohort and the nevus PRS model performance was closely similar to that in QSKIN II where R2 was 4%
and beta was 0.20 (SE = 0.060, p = 0.0009, n = 265). Ocular nevi are a risk factor for the rare and deadly
uveal melanoma. We tested whether our PRS for cutaneous nevi was associated with iris nevi scored in
the QIMR Brisbane Twin Nevus Study (BTNS) cohort. After generating a new PRS without the BTNS
cohort, this PRS was weak but significantly associated with iris nevus count (R2 = 0.0018, beta per SD of
the PRS = 0.03, SE = 0.016, p = 0.048, n = 2,607)

Methods
GWAS cohorts
We used GWAS summary statistics for nevus count from 11 studies that were reported in a previously
published nevus GWAS meta-analysis of 52,236 individuals of European ancestry5. These studies were
combined with new data from the QSkin Sun and Health Study cohort35 Phase I (QSkin I, n = 15,346)
and Phase II (QSkin II, n = 6,608), and the Australian Genetics of Depression Study (AGDS, n = 11,775)
36. A subset of QSkin I (n = 12,930) with nevus count that was non-overlapping and unrelated (identical-
by-descent score < 0.15) to the QSkin melanoma case-control set was included in a previous pleiotropic
GWAS analysis with melanoma6. However, in the current GWAS of QSkin I, we only excluded QSkin
I melanoma cases, increasing the number of available samples for this analysis of nevus loci (n = 15,346).
Details of each cohort including genotyping, quality control, imputation, and phenotypic descriptions, are
provided in Supplementary Methods, sections 1.1 and 1.2. All study participants over 18 years provided
written informed consent; for those under 18 years consent was obtained from parents and full details of
ethics approvals (QSkin I: P1309 and P2034, QSkin II: P3434, AGDS: P2118) have been reported
elsewhere5,35,36. Our analysis only included individuals of European ancestry with no melanoma history,
similar to previously published nevus GWAS meta-analysis5.

Polygenic risk score application cohorts
KYAMS: The KYAMS study was a follow-up of participants (aged 25-30 years) in the Kidskin Study,
which began in 2015 to determine the impact of sun exposure intervention during childhood on myopia
and eye health in adulthood37. The Kidskin study was initiated in 1995, a non-randomised controlled trial
designed to reduce sun exposure in children through an educational intervention for children aged 5 to 6
years from 33 schools in Western Australia38. The KYAMS was approved by the Human Research Ethics
Committees of the University of Western Australia (RA/4/1/6807) and written informed consent was
obtained from all the participants37. For our analysis, we used the number of melanocytic nevi on the right
arm counted by a trained observer. Bio samples for 113 participants were analysed in 2016, while 183
were analysed in 2021 using Infinium Global Screening Array. For quality control of the genotyped data,
SNPs with a call rate < 0.95, HWE p <10-6, and MAF (minor allele frequency) < 0.01 were excluded.
Population outliers were identified using PCA (principal component analysis) based on participants of
known European ancestry, using the 1000 Human Genome Project reference panel. Samples that were
more than six standard deviations from the PC1 and PC2 centroids were excluded. One from each pair of
individuals with an identity-by-descent score > 0.2 was removed. Post-quality control, additional genetic
variants were imputed using the Trans-Omics for Precision Medicine (TOPMed) reference panel39. For
our analysis, high-quality imputed data (r2 > 0.5 and MAF > 0.01) and nevus count data were available
for 265 individuals.

BTNS iris nevi: The BTNS cohort has measurements for iris nevi in addition to the skin nevi used in the
nevus GWAS meta-analysis. The cohort comprises a sample of predominantly adolescent twins (age 9-
23 years) and their parents from Southeastern Queensland, Australia. The cohort was initially established
to investigate the genetic factors influencing nevus development and other risk factors for melanoma40,41.
The study was approved by the QIMR Human Research Ethics Committee and written informed consent
was obtained from all participants. Genotyping was performed on Illumina 610-Quad arrays for 2,327
individuals phenotyped to the end of 2007 and on Illumina CoreExome arrays for 934 individuals
phenotyped from 2008 to May 2013. After applying a series of quality controls, genotyped data were
imputed using 1000 Genomes version 3 as described elsewhere5. The frequency of iris nevi derived from
high-resolution digital photos of the iris was evaluated on an ordinal scale; (1) an absence of iris nevi, (2)
at least one iris nevus, (3) at least three iris nevi, (4) at least five iris nevi, based on the number of melanin
accumulations on the anterior border layer of the iris42. Although both the right and left iris were
photographed, only the score from the right iris was included in the analysis due to the high polychoric
correlation (0.91 to 0.95) between the two sides as reported in a previous study42. Post-quality-controlled
(r2 > 0.5 and MAF > 0.01) imputed dosage data and iris nevus count data were available for 2,607
individuals.

GWAS for QSkin I, QSkin II and AGDS
QSkin has two self-reported nevus measurements recorded as (1) a count on the left upper arm and (2) a
4-point categorical scale (“none”, “few”, “some”, “many”). The categorical nevus variable can be
analysed as a linear or ordinal trait. We used QSkin I data to compare the performance of a linear model
in Scalable and Accurate Implementation of the GEneralised mixed model (SAIGE)43 and the proportional
odds logistic mixed model (POLMMC) which is implemented for ordinal trait analysis44. Both methods
control sample relatedness. The performance of these GWAS models was assessed by comparing the
number of significant hits and the total SNP heritability estimated using linkage disequilibrium score
regression (LDSC)45. A linear model in SAIGE was selected (Supplementary Methods, section 1.2 and
Supplementary Fig. 1-2) and was then used to conduct GWAS in the QSkin II (Supplementary Fig. 4) and
AGDS cohorts (Supplementary Fig. 5).

Heritability and LD score intercept estimation
We estimated the proportion of variance explained by all SNPs (SNP-based heritability) for nevus count
using QSkin I through genome-based restricted maximum likelihood (REML) implemented in GCTA
v.1.2646. Following quality control, a genetic relationship matrix (GRM) was constructed using genotyped
SNPs with MAF > 0.01. The GRM was fitted in with a REML analysis with 10 PCs included as covariates;
a GRM cutoff of 0.05 was then applied to exclude related people.

LDSC regression (v1.0.0) was conducted for nevus GWAS summary results using pre-calculated LD
scores from a 1000 Genomes European reference population, to estimate the genomic inflation factor, and
to estimate the SNP-based heritability for GWAS where only summary statistics were available45.

GWAS meta-analysis
Our meta-analysis was based on 14 GWASs of nevus count, with X chromosome data available in 8 of
these GWASs (Supplementary Table 2 and Supplementary Methods, section 1.2). As discussed in the
cohort details (Supplementary Methods, section 1.2), there were slight differences in the assessment of
nevus count, individual-level genotyping, imputation and quality control criteria, and the statistical
models and analytical software employed in each GWAS. Therefore, we conducted a preliminary analysis
to demonstrate that each study was consistent for conducting a GWAS meta-analysis for the nevus count
trait. Specifically, after aligning effect alleles across all studies, we performed inverse variance weighted
Mendelian Randomisation (IVW MR)47 as implemented in the “MendelianRandomisation” R package.
The nevus count of each QSkin I, QSkin II and AGDS study was obtained on a similar 4-point scale self-
reported phenotypic assessment with moderate to high levels of reproducibility48,49 and the effect sizes of
genome-wide significant SNPs were consistent across these three GWASs (Supplementary Fig. 18).
Hence, we used the GWAS meta-analysis conducted for QSkin I, QSkin II and AGDS (Supplementary
Fig. 6) as the reference study, and its top independent genome-wide significant SNPs as the instrumental
variables to perform IVW MR. Based on the study-specific IVW MR plots (Supplementary Figs. 7-17),
all individual studies from Duffy et al. (2018) were consistent with the three new studies.

Because the measurement of nevus count varies across studies, the meta-analysis of nevus GWAS was
conducted using a fixed-effect sample size weighted Z-score method as implemented in METAL v.2020-
05-0550. We applied genomic control correction to the input GWAS summary statistics. Meta-analyses
were performed separately on autosomes and X chromosome. As no genome-wide significant associations
were detected on the X chromosome, all downstream post-GWAS analyses were restricted to autosomal
results. For post-GWAS analysis that required effect sizes and standard errors, for each SNP, we
estimated the effect size and standard error using the Z-statistics, as these were the only outputs
available51.

Before conducting any post-GWAS analyses, the following quality control steps were applied to the meta-
GWAS summary statistics. First, GWAS variants driven by a single study were dropped to ensure SNPs
were associated in multiple studies. The DENTIST (Detecting Errors iN analyses of summary
staTISTtistics) method (v.1.3.0.0) was used with default settings to obtain improved GWAS summary
statistics by reducing heterogeneity between GWAS summary data and LD reference data52. In DENTIST
analysis, the LD structure was obtained using a UK Biobank European reference panel comprising 5,000
randomly selected and unrelated individuals. In this step, 8,783,913 variants were tested by applying the
DENTIST method and 124,531 variants were filtered out.

Then, the FUMA platform v.1.6.053 was used to identify independent significant SNPs and significant
loci using pre-calculated LD structure from the 1000 Genomes Phase 3 European reference population (n
= 503). Independently significant SNPs were defined as genome-wide significant (< 5x10-8) SNPs that
were independent at LD of r2 < 0.1. Genomic loci associated with nevus count were defined as LD blocks
of independent significant SNPs within a 1 Mb window. In addition, candidate SNPs that were defined as
the SNPs having r2 ≥ 0.1 with one of the independent significant SNPs and MAF 0.01, were used for
further annotations in FUMA. Forest plots of genome-wide significant lead SNPs were observed to check
that the effect sizes were consistent across studies (Supplementary Fig. 18).

GWAS-PW for nevus count and melanoma
GWAS-PW v.0.2111 was used on nevus GWAS meta-analysis data and previously published melanoma
GWAS6 to determine whether the significant genetic loci from nevus GWAS meta-analysis were
pleiotropic for nevus count and melanoma. For both nevus count and melanoma, SNPs were assigned to
LD blocks using the recommended boundaries provided by GWAS-PW
(https://bitbucket.org/nygcresearch/ldetect-data). GWAS-PW estimates posterior probability (PPA) for
four different models: a locus (i) has an association with nevus count only; (ii) has an association with
melanoma only; (iii) has a shared association with both melanoma and nevus count through the same
genetic variant; (iv) is associated with both melanoma and nevus count through independent genetic
variants. Any genome-wide significant locus with PPA > 0.8 for hypothesis (iii) was defined as a
pleiotropic locus for melanoma and nevus count.

Functional annotation, gene-based test and gene-set analysis
Candidate SNPs were annotated for functional consequences on genes based on Ensembl genes [build 85]
annotated using ANNOVAR in FUMA53.

Then positional, expression quantitative trait loci (eQTL) and chromatin interaction mapping approaches
were conducted in FUMA to prioritise genes. First, annotated SNPs were mapped to genes, considering
the physical position of the SNP inside a gene using a 10 kb window. Further, eQTL mapping was
conducted within a 1 Mb window (cis-eQTL) to test whether GWAS SNPs are associated with the
expression of the gene at a false discovery rate (FDR) ≤ 0.05. eQTL data included were restricted to skin
tissue from TwinsUK and Genotype-Tissues Expression (GTEx) v8 (sun-exposed, not sun-exposed skin
types and fibroblasts). Moreover, chromatin interaction mapping was performed to test whether there was
a significant (FDR ≤ 1×10-6) chromatin interaction between nevus-associated regions and the promoter
of nearby genes. Promoters were defined as 250 bp upstream and 500 bp downstream of the transcription
site. Although skin-specific chromatin data were not available in FUMA, we included chromatin
interaction data for all other 14 tissues (e.g. ovary, bladder, aorta and lung) and seven cell types (e.g.
embryonic stem cell, liver, mesenchymal stem, and fibroblast). This approach enabled a comprehensive
gene mapping framework and the capture of genes overlapping with chromatin interaction regions that
may be commonly expressed across tissues. Genes were prioritised only if supported by overlapping
evidence from positional, eQTL and chromatin interaction mapping.

In addition, MAGMA (Multi-marker Analysis of GenoMic Annotation)54 gene-based and gene‐set
analyses were carried out as implemented in FUMA53. In the gene-based analyses, SNPs were mapped to
18,915 protein-coding genes obtained from Ensembl build 85. After accounting for LD, the association
test statistics for SNPs were combined into a single gene-level p value. Gene-wide significance was
defined at p = 0.05/18915 = 2.64×10-6.

Here we conducted two separate gene-set analyses using candidate genes selected from MAGMA gene-
based analysis and three gene mapping approaches. In the MAGMA gene-set analysis, MAGMA gene-
based p values were further combined into each gene set from 4,761 curated gene sets and 5,917 gene
ontology (go) terms from MsigDTB v202353. Significant gene sets were selected based on Bonferroni-
corrected p values. In addition to MAGMA gene-set analysis, genes prioritised by any of three gene
mapping techniques were tested for overrepresentation in gene sets obtained from MsigDB,
WikiPathways and the GWAS catalog if at least two prioritised genes overlapped with a testing gene set53.
P values of gene sets were adjusted using the Benjamini-Hochberg correction.

Transcriptome-wide association study and colocalisation analysis for nevus count
As an alternative approach to identifying candidate genes, a TWAS was performed for the nevus GWAS
meta-analysis. eQTL data was sourced from skin tissue in GTex v8 (sun-exposed skin; n = 508, not-sun-
exposed skin; n = 430, and fibroblasts; n = 403), whole-blood (n = 31,684) from eQTLGen55 and a set of
106 primary melanocyte samples6,56. For melanocytes and other skin tissues in GTEx, pre-computed gene
expression weights are available in a previous study6,56 and FUSION website
(http://gusevlab.org/projects/fusion/) respectively, and we performed TWAS for these tissues in FUSION
using them. For whole-blood tissue, only eQTL summary-level data are available from the eQTLGen
consortium. Since FUSION requires individual-level expression data to compute expression weights, we
performed TWAS for whole-blood tissue in SUMMIT v.1.0.2 (Summary-level Unified Method for
Modeling Integrated Transcriptome) which supports summary-level expression data57. This analysis
utilised pre-calculated gene expression imputation models in SUMMIT specifically built for eQTLGen
data (https://doi.org/10.17605/OSF.IO/7MXSA).

For each set, a Bonferroni corrected p value threshold was defined as 0.05/(no of tested genes) in all five
tissue sets (melanocytes; 3797genes, three skin tissues; (9573+9556+10661) = genes, whole blood; 11415
genes, p value is 0.05/(3797+9573+9556+10661+11415) = 1.11×10-6).

Colocalisation analysis was conducted to estimate the posterior probability of a shared putative causal
genetic variant between gene expression and nevus count (if PP4 ≥ 0.75), as a complementary analysis
for TWAS 58. For melanocytes and skin tissues, colocalisation analyses were conducted within FUSION,
which includes an interface to COLOC v.5.1.0.1 R package59. A p value threshold (--coloc_P flag) of 0.05
was used to compute COLOC statistics in FUSION. For eQTLGen data, colocalisation analysis was
performed using COLOC v.5.1.0.1 package in R v.4.2.0 with default priors (p1 = 10−4, p2 = 10−4,
p12 = 10−5; the prior probabilities that an SNP is significantly associated with the expression level of a
gene, the trait outcome, or both, respectively).

Nevus counts by sex
We tested the difference in nevus count across sex using data from QSkin I (n = 15,346, Supplementary
Table 3). Nevus count data were on a 4-point scale; therefore, the ordinal logistic regression model was
used to test the effect of sex, after adjusting for age and ancestry (top 10 PCs) in R v.4.3.1. Then, potential
sex-specific genetic effects were assessed in autosomal SNPs and X chromosome using data from QSkin
I, QSkin II and AGDS. Sex-stratified GWASs were conducted for all three datasets using SAIGE (v.1.1.4)
and then meta-analysed using the inverse-variance weighted fixed-effect method in METAL v.2020-05-
0550. To determine the statistical difference of sex-specific genome-wide significant loci in males and
females, we performed a z-test60 that compares the beta-estimates for the top lead SNPs in loci from sex-
stratified meta-analysed GWASs identified by FUMA v.1.6.0. Genetic correlation between males and
females for nevus count was assessed using LDSC regression as described under “Heritability and LD
score intercept estimation”.

Mendelian randomisation for nevus count on melanoma risk
We used MR to assess the causal association between genetically determined nevus count and melanoma
risk. Summary statistics for the outcome (melanoma) were obtained from a previous larger GWAS meta-
analysis6 and analysis was restricted to clinically-confirmed set including 30,134 cases and 81,415
controls of European ancestry. Genetic instruments were selected by clumping SNPs associated with
nevus count (p < 5×10-8) from the current GWAS meta-analysis, based on LD (r2 < 0.001 within a 10mb
window), using a reference panel from UK Biobank that included 5000 randomly selected unrelated
individuals (37 SNPs). Further, we excluded ambiguous SNPs with allele frequency above 0.42. The F-
statistic values > 10 was used to select strong instruments evaluating the relevance assumption. Cochran's
Q statistic was used to assess heterogeneity in overall SNP effect estimates; if significant heterogeneity
was detected, we excluded SNPs with SNP-level Q statistics exceeding 3.8461. We estimated the causal
association between nevus count and melanoma using the inverse variance weighted method and
performed sensitivity analyses using MR-Egger regression, weighted median, both simple and weighted
mode estimators, to provide valid estimates under different assumptions47,62. Further, the presence of
unbalanced horizontal pleiotropy was evaluated using the MR-Egger intercept. Analyses were conducted
in R v.4.3.1 using the R package “TwoSampleMR”.

Polygenic risk score for nevus count
A PRS for nevus count was derived using SBayesRC which uses functional genomic annotations to
improve polygenic prediction63. To generate PRS for testing in young adults and older adults, as a
discovery dataset, we conducted a nevus GWAS meta-analysis using 13 out of 14 cohorts, excluding
QSkin II (final n = 79,357), whereas to predict iris nevus count, a separate (cutaneous) nevus GWAS
meta-analysis was conducted, excluding BTNS (final n = 82,704) as above. In addition, LD data
calculated from UK Biobank European ancestry and functional annotation data provided by SBayesRC
(https://github.com/zhilizheng/SBayesRC) were used in SBayesRC analysis.

We used QSkin II with individuals aged above 35 years (n = 6,307), KYAMS (n = 265) and BTNS (n =
2,607) cohorts to determine the ability of nevus PRS to predict cutaneous nevus count in older adults
(average age = 64 years, SD = 11 years, range = 35-93 years), young adults (average age = 28 years, SD
= 1 year, range = 25-30 years), and then iris nevus count respectively. In addition, we used QSKIN II data
to compare the prediction performance of the new SBayesRC-derived PRS with an LD-clumping-based
PRS (collinearity r2 < 0.05 and 2Mb window) generated previously from the prior meta-analysis
GWAS5,15. To score individuals in independent cohorts, imputed allelic dosages were weighted by the
SBayesRC-derived SNP weights using PLINK v.2.0 (www.cog-genomics.org/plink/2.0/). We then
standardised the nevus PRS to have a mean of zero and SD of 1. Since related individuals were retained
in all cohorts, GCTA GREML was then used to regress inverse-normal rank-transformed cutaneous nevus
count/iris nevus count, adjusted for sex, age, age2, sex×age, sex×age2, the first 10 PCs and a GRM, onto
standardised PRS. For Kidskin and BTNS cohorts, we used a binary batch variable as an additional
covariate to control the genotype batch effect. Finally, nevus PRS performance was determined using R2
computed by applying a formula reported in64.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an important study - the largest to date of genetic variant associations with nevi. The authors have identified new loci\nof interest.\nI have just a few minor comments:\nAcronyms: FUMA is not defined.\nPRS not defined until page 20 although used earlier.\nI think for more general audience understanding, it would be helpful to be very specific about the\nSNPs used in the PRS and their weights. Maybe a supplementary table with the PRS snps and the weights used would\nmake this easier to understand and for researchers to apply to other datasets in the future. Also, how do these weights from\nSBayesRC compare with the weights previously used in these author’s other papers? Also how does the PRS prediction\ncompares to prior PRS scores. What do the extra SNPs and new method of weighting add to the prediction?", "reviewer_2": "In their paper “A large-scale genome-wide association meta-analysis for nevus count provides direct insights into the\ngenetics of melanoma”, Jayasinghe et al. performed the largest GWAS meta-analysis to date of nevus count using\npreviously published data, and adding two more cohorts, all corresponding to European genetic ancestry. They uncovered\nnovel loci associated with nevus count, which have not been associated previously with other melanoma or pigmentation\ntraits, and performed gene-based, as well as TWAS and pathway analyses to understand the biological and molecular\nprocesses involved, and their relationship with melanoma. Additionally, they created polygenic risk scores for nevus count\nand additionally for iris nevi, which are a risk factor for uveal melanoma. Overall, the authors provide important contributions\nto the field of pigmentation and melanoma genetics, using a broad range of computational approaches.\nMajor comments\n1. The authors performed GWAS in two new cohorts to perform a meta-analysis with previously published data. I\nunderstand that until recently, chromosome X in GWAS has been often neglected and rarely included, therefore I do\nnot expect the authors to re-analyze other cohorts’ data for chromosome X. However, it is now becoming fairly\ncommon, and GWAS tools such as SAIGE allow the inclusion of chromosome X in the analysis. Therefore, the\nauthors should be able to easily include the X chromosome in the two new cohorts they analyze in the present\npaper, which will strengthen the paper (even if no signals are identified), and will help pave the way for other future\nstudies to include it.\n2. The authors should report the age range of participants across cohorts used in the meta-analysis and the\nproportion of participants across the sexes (at least in the two new cohorts included). It would also be relevant to\nassess if there are phenotypic (and if so, genetic) differences between the sexes, as it is known that there are\ndifferences in the incidence of melanoma between males and females.\n3. In Supplementary Note 1.2, the authors argue that the SNP-based heritability estimates were slightly higher for\ntheir SAIGE approach, compared to POLMM, hence choosing the SAIGE approach for their downstream analyses.\nHowever, considering the SE of the estimates, they are basically the same. Therefore, the argument to choose\nSAIGE over POLMM seems a little weak. On one hand, one can argue that SAIGE is one of the standard GWAS\nmixed-models tools used by the community, whereas on the other hand, a model that considers an ordinal scale\nwould be the better choice due to the nature of the data used for the GWAS.\n4. The authors should provide more details in the Methods regarding their colocalization analyses, specifically the\npriors utilized. Additionally, in the Results section there is no mention of the follow-up results of TWAS with\ncolocalization. The authors should specify in their TWAS results that candidate genes were selected based on\nTWAS followed-up by colocalization, if that is the case (which seems to be). Finally, from the main text Results,\nthe authors should specify which TWAS analyses (i.e., which cell-type/tissue) yielded the significant results and if the\ndirection of effect was concordant when the signal was in multiple TWAS analyses (i.e., GDI2), especially given the\nless relevant whole-blood dataset used.\n5. The authors should provide the reasoning behind the inclusion of 14 different tissues for the chromatin\ninteraction mapping, given that this is a broad list of tissues/cell types that do not necessarily have an evident\nrelation to skin.\n6. The finding of a possibly antagonistic effect between nevus count and melanoma in the MC1R locus is pretty\ninteresting. As an alternative explanation, the authors mention in the Discussion that this might be a consequence\nof the difficulty in identifying nevi in individuals with red hair and very low skin pigmentation. However, I do not see\nhow specifically red (scalp) hair can possibly influence nevi count (as the authors suggest). In contrast, I could\nimagine how the presence of freckles (where MC1R is involved too) could influence nevi count.\n7. One of the main conclusions of the authors is the possibility to include nevus count genetic signals as risk\nfactors for melanoma PRS, and they have shown with pair-wise analysis the overlap of the genetic signals between\nboth traits. However, including a causal inference approach, such as Mendelian randomization, using the nevus\ncount genetic signals as exposures and melanoma GWAS as outcome, would aid in understanding the relationship\nbetween these two traits, and would potentially strengthen the conclusions of the paper.\nMinor comments\n1. Citations in the introduction are in different format (i.e., First sentence of second paragraph).\n2. In S1.2, there is a typo ‘identify-by-descent’ should be ‘identity-by-descent’.\n3. Readers would appreciate if authors could specify the genome build of the genomic coordinates used in the\nSupplementary Tables. Additionally, it would be appreciated if the column names throughout the Supplementary\nTables could be better described, as I only see an ‘Abbreviations’ tab for one table.\n4. The authors should re-organize the order in which Supplementary Tables appear. It is strange that the first\nSupplementary Table referenced in the Results section of the main text is Supplementary Table 7, whereas other\nSupplementary Tables are not referenced throughout the main text. Similarly, the Methods section for GWAS-PW\nshould be put in the same order as in the Results section.\n5. In the Functional Mapping Results section, there seems to be a word missing from the sentence: “Of 248 genes,\n26 mapped on 14 of 29 nevus-associated genomic loci that were identified by positional, eQTL and chromatin\ninteraction mapping”.\n6. For the nevus PRS, what were the age ranges selected to classify in older vs. young adults? The authors should\nspecify the age ranges in the Methods, not only the average age of each cohort."}, "author_response": "Response to Reviewers\nReferee: 1\nReviewer #1 (Remarks to the Author):\nThis is an important study - the largest to date of genetic variant associations with nevi. The authors have identified new loci\nof interest.\nI have just a few minor comments:\n1. Acronyms: FUMA is not defined.\nResponse: Corrected by defining FUMA at the first occurrence (Page 6, line 134)\n2. PRS is not defined until page 20, although used earlier.\nResponse: Corrected by defining PRS at the first occurrence (Page 11, line 241)\n3. I think for a more general audience understanding, it would be helpful to be very specific about the\nSNPs used in the PRS and their weights. Maybe a supplementary table with the PRS SNPs and the weights used would\nmake this easier to understand and for researchers to apply to other datasets in the future.\nResponse: Together with the meta-analysis GWAS summary statistics, the PRS SNPs with\nweights will be available in Zenodo, as the PRS includes 7,280,843 SNPs. The data availability statement\nhas been updated to more clearly report this.\n4. Also, how do these weights from SBayesRC compare with the weights previously used in these authors’ other\npapers?\nResponse: The weights used in the current PRS differ from those in the authors’ previous\npublications for two main reasons. First, the meta-analysis methods and thus the weight scale differ: the\nprevious PRS (e.g. used in PMID: 34778946) was derived from an inverse-variance weighted\nmeta-analysis, whereas the current PRS was based on a p-value–based meta-analysis. Second, whereas the\nprevious PRS used a limited number of independent genome-wide SNPs (8 SNPs), SBayesRC integrates\nover ~7 million SNPs and applies LD-aware modelling combined with functional annotation-informed\nBayesian shrinkage to produce more accurate SNP effect estimates, and to improve polygenic prediction\nof traits. As a result, direct comparison of the weights is not straightforward. However, they are more\nreadily compared using their predictive ability, which is answered in the reviewer’s next question.\n5. Also how does the PRS prediction compares to prior PRS scores. What do the extra SNPs and new method of\nweighting add to the prediction?\nResponse: Compared with the 2018 clumping-based PRS (R² = 0.36%, β = 0.051, SE = 0.01),\nthe SBayesRC-derived PRS shows substantially improved predictive performance (R² = 5%, β = 0.19, SE\n= 0.01). The results were added to the manuscript under the PRS analysis (Page 11, lines 245-246).\nReferee 2:\nReviewer #2 (Remarks to the Author):\nIn their paper “A large-scale genome-wide association meta-analysis for nevus count provides direct insights into the genetics of\nmelanoma”, Jayasinghe et al. performed the largest GWAS meta-analysis to date of nevus count using previously published\ndata, and adding two more cohorts, all corresponding to European genetic ancestry. They uncovered novel loci associated with\nnevus count, which have not been associated previously with other melanoma or pigmentation traits, and performed\ngene-based, as well as TWAS and pathway analyses to understand the biological and molecular processes involved, and their\nrelationship with melanoma. Additionally, they created polygenic risk scores for nevus count and additionally for iris nevi,\nwhich are a risk factor for uveal melanoma. Overall, the authors provide important contributions to the field of pigmentation\nand melanoma genetics, using a broad range of computational approaches.\nMajor comments\n1. The authors performed GWAS in two new cohorts to perform a meta-analysis with previously published data. I\nunderstand that until recently, chromosome X in GWAS has been often neglected and rarely included, therefore I do not\nexpect the authors to re-analyze other cohorts’ data for chromosome X. However, it is now becoming fairly common, and\nGWAS tools such as SAIGE allow the inclusion of chromosome X in the analysis. Therefore, the authors should be able to\neasily include the X chromosome in the two new cohorts they analyze in the present paper, which will strengthen the paper\n(even if no signals are identified), and will help pave the way for other future studies to include it.\nResponse: We thank the reviewer for the helpful suggestion. We agree that the X chromosome\nhas been under-analysed in GWAS, but that including it is increasingly feasible and will strengthen the\nfield moving forward. We have now included an analysis of the X chromosome in QSkinI, QSkin II and\nAGDS using SAIGE. In addition, data on the X chromosome were available from 5 GWASs included in\nthe previous meta-analysis (PMID: 30429480). We then integrated these GWASs into a meta-analysis of\nchromosome X. However, no genome-wide significant associations were detected in either cohort or the\nmeta-analysis. We have reported the methods and results related to this analysis in the manuscript and will\ninclude chromosome X results in the shared GWAS summary statistics (Methods: Page 21, line 474 and\nSupplementary Note S1.2 Page 29; Results: Page 5, lines 103-106)\n2. The authors should report the age range of participants across cohorts used in the meta-analysis and the proportion of\nparticipants across the sexes (at least in the two new cohorts included). It would also be relevant to assess if there are\nphenotypic (and if so, genetic) differences between the sexes, as it is known that there are differences in the incidence of\nmelanoma between males and females.\nResponse: The age ranges were available for 11 studies from the previous meta-analysis (PMID:\n30429480), and so were added to Table S1.1 in Supplementary Note (Page 25-26). Age ranges and\nproportions of participants across sex were reported for the new cohorts under the cohort description\n(Supplementary Note S1.2: Page 27-28).\nWe also thank you for your thoughtful comment regarding the addition of sex-specific analyses to the\nmanuscript. We assessed the phenotypic differences in the nevus count across sex using QSkin I data and\nobserved a higher total-body nevus count in women than in men. Therefore, we performed sex-specific\nGWAS by meta-analysing QSkin I, QSkin II, and AGDS to identify the difference in the genetic\narchitecture of nevus count between men and women. The genetic correlation was not significantly\ndifferent to 1 (rg = 0.863, 95% CI = 0.453 – 1.273), indicating that the genetic architecture of nevus count\nwas largely shared across men and women. We have updated the methods section and added the results to\nthe manuscript under “Sex-specific GWAS meta-analysis for nevus count” (Methods: Page 26, lines\n570-582; Results: Page 10, lines 213-227; Discussion: Page 15, lines 333-338).\n3. In Supplementary Note 1.2, the authors argue that the SNP-based heritability estimates were slightly higher for their\nSAIGE approach, compared to POLMM, hence choosing the SAIGE approach for their downstream analyses. However,\nconsidering the SE of the estimates, they are basically the same. Therefore, the argument to choose SAIGE over POLMM\nseems a little weak. On one hand, one can argue that SAIGE is one of the standard GWAS mixed-models tools used by\nthe community, whereas on the other hand, a model that considers an ordinal scale would be the better choice due to the\nnature of the data used for the GWAS.\nResponse: We agree that the SNP-based heritability estimates (h2) from SAIGE and POLMM are\nvery similar when considering their standard errors. While the SEs overlapped, we acknowledge that\nselecting SAIGE based on the fact that its point estimate was higher has the advantage of aligning with\nthe generally used method. We have revised the text in the supplementary note 1.2 to acknowledge these\nfacts (Supplementary Note S1.2: Page 28).\n4. The authors should provide more details in the Methods regarding their colocalization analyses, specifically the\npriors utilized. Additionally, in the Results section there is no mention of the follow-up results of TWAS with\ncolocalization. The authors should specify in their TWAS results that candidate genes were selected based on\nTWAS followed-up by colocalization,\nif that is the case (which seems to be). Finally, from the main text Results, the authors should specify which TWAS analyses\n(i.e., which cell-type/tissue) yielded the significant results and if the direction of effect was concordant when the signal was in\nmultiple TWAS analyses (i.e., GDI2), especially given the less relevant whole-blood dataset used.\nResponse: We have expanded the Methods section to detail the priors used in our colocalization\nanalyses (Page 25-26, lines 566-568). In the Results, we specified that candidate genes were identified\nthrough TWAS, followed by colocalization, along with the tissues that showed significant TWAS\nassociations for each gene. For genes like GDI2, which exhibited signals in multiple tissues, we have\nincluded the directions of effect and noted their concordance across tissues in the results (Page 8-9, lines\n184-190).\n5. The authors should provide the reasoning behind the inclusion of 14 different tissues for the chromatin interaction\nmapping, given that this is a broad list of tissues/cell types that do not necessarily have an evident relation to skin.\nResponse: Chromatin interaction mapping was performed using FUMA. Since skin-specific\nchromatin data were not available in FUMA, we included all other available 14 tissues and 7 cell types to\nprovide a comprehensive gene-mapping framework and to capture genes that overlapped with chromatin\ninteraction regions that may be commonly expressed across tissues. Importantly, genes were prioritised\nonly if supported by overlapping evidence from chromatin interaction, positional mapping, and eQTL\nmapping conducted using skin tissues. We have added this clarification to the methods (Page 24, lines\n522-528).\n6. The finding of a possibly antagonistic effect between nevus count and melanoma in the MC1R locus is pretty\ninteresting. As an alternative explanation, the authors mention in the Discussion that this might be a consequence\nof the difficulty in identifying nevi in individuals with red hair and very low skin pigmentation. However, I do not see\nhow specifically red (scalp) hair can possibly influence nevi count (as the authors suggest). In contrast, I could imagine how the presence of freckles\n(where MC1R is involved too) could influence nevi count.\nResponse: We thank the reviewer for this thoughtful comment. We agree that attributing reduced\nnevus counts specifically to red hair needs to be correct; if there is a role, it is likely to act on the skin or\nmoles directly. We appreciate the suggestion to also consider the role of freckling more directly. As noted\nin a previous work, individuals within pale-skinned, highly freckled groups often show fewer nevi;\nhowever, molecular analyses have demonstrated that MC1R RHC homozygotes exhibit lower nevus\ncounts as freckling increases across all freckling levels. This pattern indicates that the association is\nunlikely to be explained solely by difficulty in counting nevi on lightly pigmented skin. The MC1R RHC\nvariants likely reduce nevus formation by disrupting cAMP signalling, which is normally required for the\nearly melanocyte proliferation that forms a mole (PMID: 14709592), while simultaneously increasing\nmelanoma risk by impairing UV protection in melanocytes. We have revised the Discussion to reflect this\nmechanistic, rather than phenotypic-artifact, explanation (Page 13, lines 289-300).\n7. One of the main conclusions of the authors is the possibility to include nevus count genetic signals as risk factors for\nmelanoma PRS, and they have shown with pair-wise analysis the overlap of the genetic signals between both traits. However,\nincluding a causal inference approach, such as Mendelian randomization, using the nevus count genetic signals as exposures\nand melanoma GWAS as outcome, would aid in understanding the relationship between these two traits, and would\npotentially strengthen the conclusions of the paper.\nResponse: We thank the reviewer for this helpful suggestion. We have now performed Mendelian\nrandomisation using nevus count–associated variants as instruments and used a previous melanoma\nGWAS (PMID: 32341527) as the outcome. This found that genetically predicted higher nevus count was\nsignificantly associated with a higher risk of melanoma (Inverse variance weighted (IVW) OR = 4.212,\n95% CI = 3.469-5.114). We have updated the results and methods sections, integrating MR analysis into\nthe manuscript (Methods: Page 26-27, lines 584-600; Results: Page 10-11, lines 229-238; Discussion: Page\n15, 339-342).\nMinor comments\n8. Citations in the introduction are in different format (i.e., First sentence of second paragraph).\nResponse: Corrected (Page 3, line 71).\n9. In S1.2, there is a typo ‘identify-by-descent’ should be ‘identity-by-descent’.\nResponse: Corrected.\n10. Readers would appreciate if authors could specify the genome build of the genomic coordinates used in the\nSupplementary Tables. Additionally, it would be appreciated if the column names throughout the Supplementary\nTables could be better described, as I only see an ‘Abbreviations’ tab for one table.\nResponse: Genome build has been added to each table. We have also added readable column\nnames.\n11. The authors should re-organize the order in which Supplementary Tables appear. It is strange that the first\nSupplementary Table referenced in the Results section of the main text is Supplementary Table 7, whereas other\nSupplementary Tables are not referenced throughout the main text. Similarly, the Methods section for GWAS-PW\nshould be put in the same order as in the Results section.\nResponse: We have reorganised the Supplementary tables and revised the methods section\naccordingly.\nThe order of the Methods section for GWAS-PW has been changed to align with the order of the Results\nsection.\n12. In the Functional Mapping Results section, there seems to be a word missing from the sentence: “Of 248 genes,\n26 mapped on 14 of 29 nevus-associated genomic loci that were identified by positional, eQTL and chromatin\ninteraction mapping”.\nResponse: Corrected as “Of 248 genes, 26 overlapped in all three mapping approaches, and those\ngenes mapped on 14 of 29 nevus-associated genomic loci”. (Page 6 lines 139-141)\n13. For the nevus PRS, what were the age ranges selected to classify in older vs. young adults? The authors should\nspecify the age ranges in the Methods, not only the average age of each cohort.\nResponse: Thank you for this suggestion. Initially, we classified the two adult groups (younger vs.\nolder) based on the cohort-level average ages: KYAMS (mean age = 28 years, SD = 1 year, range = 25–30\nyears) and QSKIN II (mean age = 63 years, SD = 13 years, range = 18–93 years). To provide a clearer\nseparation between the two groups, we have now restricted the PRS analysis in the QSKIN II cohort to\nindividuals aged >35 years. We have specified the age ranges for both cohorts in the Methods section\n(Page 27, lines 613-614)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "All questions have been satisfactorily answered. The work is significant to the field, the conclusions are supported.\nThe methodology is sound and there is enough detail.", "reviewer_2": "The authors have made a significant improvement to their manuscript, and it reads very well.\nI just have one minor comment left, but otherwise I recommend for publication, with regards to the use of gender\nterms in the context of sex-specific GWAS. If the authors refer to genetic sex, then the correct terms to use would be\nfemale/male instead of women/men."}, "author_response": "Referee: 1\nReviewer #1 (Remarks to the Author):\nAll questions have been satisfactorily answered. The work is significant to the field, the conclusions are supported. The\nmethodology is sound and there is enough detail.\nResponse: We thank the reviewer for the positive assessment and glad that all questions have been\nsatisfactorily addressed.\nAdditionally, we would like to note that we have improved the clarity of the data access statement in\ngeneral without changing the access status.\nReferee 2:\nReviewer #2 (Remarks to the Author):\nThe authors have made a significant improvement to their manuscript, and it reads very well.\nResponse: We thank the reviewer for the positive feedback.\nI just have one minor comment left, but otherwise I recommend for publication, with regards to the use of gender terms in the\ncontext of sex-specific GWAS. If the authors refer to genetic sex, then the correct terms to use would be female/male instead\nof women/men.\nResponse: We have updated the manuscript to use female and male when referring to genetic sex,\nreplacing women and men where applicable.\nAdditionally, we would like to note that we have improved the clarity of the data access statement in\ngeneral without changing the access status."}]

Decision letter