Immunomodulatory germline variation associated with the development of multiple primary melanoma (MPM)

Multiple primary melanoma (MPM) has been associated with a higher 10-year mortality risk compared to patients with single primary melanoma (SPM). Given that 3–8% of patients with SPM develop additional primary melanomas, new markers predictive of MPM risk are needed. Based on the evidence that the immune system may regulate melanoma progression, we explored whether germline genetic variants controlling the expression of 41 immunomodulatory genes modulate the risk of MPM compared to patients with SPM or healthy controls. By genotyping these 41 variants in 977 melanoma patients, we found that rs2071304, linked to the expression of SPI1, was strongly associated with MPM risk reduction (OR = 0.60; 95% CI = 0.45–0.81; p = 0.0007) when compared to patients with SPM. Furthermore, we showed that rs6695772, a variant affecting expression of BATF3, is also associated with MPM-specific survival (HR = 3.42; 95% CI = 1.57–7.42; p = 0.0019). These findings provide evidence that the genetic variation in immunomodulatory pathways may contribute to the development of secondary primary melanomas and also associates with MPM survival. The study suggests that inherited host immunity may play an important role in MPM development.

Recent genome-wide association studies (GWAS) [10][11][12][13][14] and candidate pathway studies 15 have identified germline genetic variants associated with cutaneous melanoma risk [10][11][12][13] ; however, the literature on MPM-specific genetic risk markers remains scarce. While most melanoma GWAS analyses were performed on patient populations with predominantly single primary diagnoses, few GWAS studies 12,13 investigated whether the presence of additional primaries modulated differential risk for specific variants. Besides GWAS, the candidate scans have been performed on known melanoma-related risk pathways and associations were found between polymorphisms in XPD 16 , TYR 17 , and MTAP 17 and MPM risk. Other studies have examined associations between MPM risk and the pigmentation pathway by assessing variations in MC1R [18][19][20] and MITF 21 , as well as vitamin D synthesis among individuals exposed to high levels of sunlight in recreational activities 22 , although these results were of borderline significance. However, besides these GWAS or candidate scans focused on genetic risk variation derived from predominantly SPM populations, the potentially novel MPM-specific risk loci have not been investigated. Also, to our knowledge, the genetic variation potentially impacting host immunity and its role in the development of additional primary tumors has not been explored. Recently, large datasets, such as the Multiple Tissue Human Expression Resource (MuTHER) 23 , linked the genetic variation with the gene expression in immune lymphoblastoid cells, allowing for opportunities to investigate the inherited basis of host immunity related to cancer outcomes. We have recently reported that the germline genetic variants linked with the expression of immunomodulatory genes (immunomodulatory expression quantitative trait loci, or ieQTLs) in lymphoblastoid cell lines (LCLs) are associated with melanoma recurrence 24 and survival 25 . Here we hypothesize that MPM may be in part attributed to impaired immune regulation during melanomagenesis in patients with SPM, and this altered host immunity is determined by inherited genetic variation. In this study, using a cohort of 977 melanoma patients, including 147 MPM patients, we sought to assess the potential of 41 ieQTLs in LCLs from MuTHER as putative biomarkers of host immunity associated with the development of additional melanoma primaries. We have also tested whether ieQTLs differentially modulate survival in patients with MPM.
Materials and Methods study population. The study population included 977 patients of self-reported non-Hispanic European descent, with no reported family history of melanoma, treated between 2002 and 2013 at New York University Langone Health (NYULH) with stage I to III cutaneous melanoma, prospectively enrolled in the NYU Interdisciplinary Melanoma Cooperative Group's Institutional (IMCG) Review Board-approved clinicopathological database and biorepository and actively followed up on a protocol-driven schedule 24,26,27 . Blood samples, demographic, follow-up, and clinical information including age at pathological diagnosis, gender, self-reported family history of melanoma, and 2009 AJCC staging at diagnosis, were collected as part of IMCG protocol, and all patients (parent or legal guardian for any patient under 18) provided written informed consent at time of enrollment. Subjects with a reported family history of melanoma were excluded from the final patient cohort. Control study participants consisted 1047 cancer-free non-Hispanic subjects from a cutaneous melanoma case-control study at MD Anderson Cancer Center from March 1998 to August 2008 (accession number phs000187.v1.p1) who were acquaintances of patients that presented at other clinics 14 . The methods in this study were carried out in accordance with guidelines set forth by the IRB at NYULH relating to the use of patient samples in genetic studies, and all experimental protocols were approved by NYULH.
Selection of immunomodulatory eQTLs and genotyping. Previously, we utilized publicly available databases to identify 382 immunomodulatory genes and these loci were used to isolate candidate SNPs of interest for genotyping 25 . Employing resources from the MuTHER project 28,29 , 50 SNPs with the most significant cis-ieQTL activity on probes in LCLs from these immunomodulatory loci were selected for genotyping. While skin and adipose tissue data are also available in the MuTHER project, the cis-ieQTLs have been assessed from the LCL expression data, as the scope of the study is focused on the host immune cell component involved in MPM development. The selection procedure was described in detail elsewhere 25 .
Genomic DNA from all 977 cases was isolated from whole blood samples using a QiaAmp kit (Qiagen). All SNP genotyping was completed using the MassARRAY System (Agena Bioscience Inc.) according to the manufacturer's protocol. To ensure high-quality genotyping, quality control filters were used to remove SNPs with call rate <90%, samples with call rates <90%, and SNPs with a significant departure from Hardy-Weinberg equilibrium (p < 1E-04) resulting in 41 ieQTLs for analysis. Additionally, rare SNPs (defined by minor allele frequency <0.05) were removed prior to logistic regression analysis. statistical analysis. Multivariate logistic regression genetic association analyses were conducted, adjusting for age at pathological diagnosis, sex (male vs. female), and Ashkenazi Jewish (AJ) status (yes vs. no), as a fraction of the patient cohort are of AJ ancestry (n = 281, 13.9%) (see Table 1). The logistic regression models were conducted comparing patients with MPM and no family history versus cancer-free controls, and those with MPM and no family history versus SPM and no family history. Because the GWAS controls did not provide AJ ancestry information, we performed Principal Component Analysis (PCA) on the control samples in order to identify any control samples of potential AJ ancestry, and then utilized this information to adjust all models for possible confounding due to AJ status. Additionally, we used whole genome sequencing (WGS) information obtained in our parallel ongoing studies for 49 of our AJ samples in this cohort to confirm the accuracy of self-reported AJ status in our clinical dataset. The PCA was done by first identifying overlapping variants among previously published AJ reference population 30 , AJ patients with WGS information, the GWAS control SNPs 14 , and genetic information from 1000 Genomes Phase 3, which served as a European population control. The 243,847 overlapping variants were pruned to a reduced set in linkage equilibrium using PLINK v1.90beta with parameters set to remove variants showing correlation r 2 > 0.2 in 500-variant windows, and finally, the non-autosomal variants were removed. The remaining 97,000 variants were used for PCA, and putative AJ status was determined by samples with principal www.nature.com/scientificreports www.nature.com/scientificreports/ Multivariate Cox proportional hazards (Cox PH) models were utilized to assess the association between SNPs and overall survival (OS) among either patients with MPM and no family history or patients with SPM and no family history of melanoma. Cox PH models were calculated from the age of diagnosis of the primary tumor until last follow-up, death, or treatment with immunotherapies and adjusted for age at pathological diagnosis, AJ status, sex (male vs. female), primary tumor histologic subtype (superficial-spreading vs. nodular vs. desmoplastic vs. acral-lentiginous vs. lentigo-maligna vs. other), and stage (I, II, III). MPM patients in this analysis were defined as those with 1 or more additional primary melanoma observed at initial diagnosis or developed subsequently after primary tumor diagnosis. Kaplan-Meier plots were produced for SNPs that reached significance after adjustment for multiple comparisons. Additive models of disease were used for the main effect logistic regression analysis, and both dominant and recessive models were selected for the survival analysis. Participants who received immunotherapy (N = 79) were censored at date of treatment. The adjustment for multiple testing was performed by Bonferroni procedure, which was employed in both logistic and survival models. Proportionality of hazards was assessed for all Cox models. All statistical analyses were conducted using PLINK 31 and the "Survival" package in R. ethical approval and informed consent. Written informed consents for the use of the blood specimens and clinical information were obtained at the time of enrollment from all participants, and the study was approved by the Institutional Review Board (IRB) at all NYULC (IRB#10362).
Results study population. The study population stems from a cohort of patients with stage I-III cutaneous melanoma who were ascertained as part of The Interdisciplinary Melanoma Collaboration Group (IMCG) at NYU Langone Health, as previously described [24][25][26] . The patients with family history of melanoma were specifically excluded from the analysis to reduce any confounding effects of known germline genetic factors associated with MPM or SPM development. 49.3% of patients were over 60 years of age and all patients were self-reported to be of European descent. The overall 5-year survival was 84.1% and the median time between diagnosis and follow up was 52.7 months. The majority of primary diagnoses were stage I and II (83%), and the most common histological subtype was superficial spreading melanoma (46.5%) followed by nodular melanoma (23.1%) and other subtypes (30.4%). The study population had no reported family history of melanoma and consists of 147 patients with multiple primary melanomas and 830 patients with single primary melanomas. For the case/control analyses, we have utilized data from 1047 control samples ascertained as part of a Genome Wide Association Study (GWAS) (phs000187.v1.p1) 14 . The clinical and demographic data for the patient population used in the current analyses are summarized in Table 1.

The stratification of Ashkenazi Jewish (AJ) ancestry in case/control populations. A fraction of
the patient population used in this analysis is of self-reported Ashkenazi Jewish (AJ) descent (n = 243, 24.9%); therefore, we ensured the stratification of AJ ancestry was appropriately adjusted. Also, the AJ ancestry status was not available for the GWAS control data. To validate the accuracy of self-reported AJ status information in IMCG cohort, we have exploited the whole genome sequencing (WGS) data available for a subset of 49 melanoma cases, sequenced as part of our parallel efforts and also used in the current analysis. Using WGS data from 49 self-reported AJ participants, the GWAS data of n = 128 AJ references samples genotyped as part of our prior study 30 , the GWAS control data 14 , and the genotyping data from the 1000 Genomes Project Phase 3, we identified 243,847 SNPs overlapping with the cancer-free GWAS control dataset. After linkage disequilibrium pruning, we performed principal component analysis (PCA) on 97,000 SNPs and plotted the two major principal components.  www.nature.com/scientificreports www.nature.com/scientificreports/ 49 samples self-reported as AJ and 38 cancer-free GWAS control samples clustered around the AJ reference set. The remaining case and control samples spread among the different European sub-populations ( Supplementary  Fig. 1) as expected. For all further analyses in this study, we considered the controls clustering in the AJ reference set to be of AJ ancestry along with those of self-reported AJ ancestry in our clinical IMCG dataset.
The association of ieQTLs with MPM risk. To explore MPM-specific associations, we performed an association analysis comparing MPM and SPM cases under an additive logistic regression model adjusting for age, sex and AJ status. The most significant association, surpassing the adjustment for multiple testing (Bonferroni threshold p < 0.001), was found for rs2071304 (OR = 0.60; 95% CI = 0.45-0.81; p = 0.0007), ( Table 2,  full results Supplementary Table 1). In order to test whether the most significant associations from this analysis were MPM-specific, we analyzed rs2071304 and rs665241 using two association tests comparing 147 MPM or 830 SPM, respectively, with 1047 healthy controls from publicly available imputed GWAS data (phs000187.v1.p1) 14 , adjusting for age, sex, and AJ status (for controls detected by PCA as described in Materials and Methods) ( Table 2, full results Supplementary Tables 2 and 3). While no significant associations have been found in SPM case/control comparison, the associations were confirmed, albeit with reduced significance, for rs2071304 in MPM case/control analysis (OR = 0.59; 95% CI = 0.41-0.83; p = 0.0025). Again, the results showed that the alternate allele G of rs2071304, which associates with decreased expression of SPI1 in LCLs (Fig. 1) in the MuTHER data 23 , confers a protective effect for MPM. It is possible that the results may be confounded by time of follow-up for SPM patients, biasing potential future diagnosis of MPM. To address this, for the 2 most significant variants from the initial analysis (rs2071304 and rs665241) we have also performed a subset analysis comparing MPM patients to those with SPM with at least 8 years of follow-up, a threshold established in prior studies 32 Table 4).
ieQTLs associated with OS of patients with multiple melanomas. We have previously reported that ieQTLs may be clinically actionable prognostic biomarkers predictive of OS among melanoma patients in the general population regardless of their MPM status 25 . Based on this suggestive evidence, we sought to test whether ieQTLs may also predict survival among patients with MPM. We performed Cox PH analysis on 147 MPM cases genotyped for 41 ieQTLs in this study. The results of the additive and dominant Cox proportional hazards model analyses are summarized in Table 3 (full results Supplementary Table 5). The most significant variant associated with MPM survival was rs6695772 (additive model: HR = 3.42; 95% CI = 1.57-7.42; p = 0.0019; dominant model: HR = 18.69; 95% CI = 3.34-104.55; p = 0.0009) (Fig. 2), linked to the expression of BATF3, which was the same ieQTL predictive of melanoma survival in our recent study 25 , with the alternate allele of rs6695772 linked to the decreased expression of BATF3 in LCLs (Fig. 1)

Discussion
A role of host immunity in controlling melanoma progression has been recently demonstrated by the successes of melanoma immune checkpoint blockade (ICB) therapies [33][34][35] . The ICB significantly improves survival of metastatic melanoma patients and there is a clear promise of ICB benefits at earlier stages as evidenced from ongoing adjuvant trials in melanoma 36 . These groundbreaking developments represent a paradigm showing that melanoma, an immunogenic tumor, attracts the attention of the host immune response which may be a modulating www.nature.com/scientificreports www.nature.com/scientificreports/ factor in tumor burden reduction. While this link has been clearly established in the context of advanced disease, it is highly plausible that the host immunity will also play a critical role in melanoma development. This is clearly supported by the studies showing that the risk of melanoma 37 , or tumor recurrence 9 are increased among individuals who are immunosuppressed, through causes such as HIV infection 38,39 or organ transplants. In this study, we tested a hypothesis that the genetic variants impacting host immune regulation may contribute to the development of multiple primary melanoma. We postulate that the immune surveillance of early stage primary tumors may be impaired by genetic factors, which in turn can be direct surrogate markers of additional primary tumor development. As we suggest, this may be a convenient clinical complement to the lengthy process of follow-up of patients with single primary melanomas, often undergoing many skin biopsies in order to reduce the likelihood of second melanoma, which is financially and emotionally burdensome. Thus, given the lifetime risk of additional primary melanomas ranging between 1-12% 40 , the availability of novel personalized biomarkers predictive of multiple primary tumors would be of significant clinical importance.
To address these assumptions, we analyzed 41 ieQTLs that associate with the expression of immune responsive genes in immune lymphopbasloid cells, from a healthy twin cohort in the MuTHER consortia 23,25 , in order to identify variants predictive of risk of MPM development. We have previously reported that such ieQTLs impact melanoma survival 25 . Here, we report significant associations of these 41 ieQTLs with the development of MPM. These findings suggest that immune modulation controlled by inherited genetic variants may be a contributing factor affecting the development of additional primary tumors in patients with single primary melanoma. The most significant result was found for rs2071304 (OR = 0.60; 95% CI = 0.45-0.81; p = 0.0007) comparing MPM patients with single primary cases. We found that patients who carried the alternate allele (G) in the study population are 40% less likely to develop MPM. In this analysis, the consideration of time to second primary from single primary diagnosis is critical, and as such we have also tested whether this association remained significant after accounting for sufficient follow-up time after initial SPM diagnosis. Importantly, rs2071304 remained among the top 2 most significant results after considering only SPM patients with at least 8 years of follow up after single primary diagnosis, an interval suggested in prior studies 32 . In addition, when we compared MPM samples with healthy controls, we found that rs2071304 showed borderline significance (p = 0.0025), which was not observed in the analysis of SPM patients versus healthy controls, further suggesting that the association is MPM-specific. While these findings suggest that rs2071304 may modulate risk of MPM, the data also raise a possibility that the risk alleles affecting the development of second primary melanomas from single primary diagnoses may   www.nature.com/scientificreports www.nature.com/scientificreports/ not be necessarily identifiable by MPM case/control analysis. In controls, the risk alleles are not under selection pressure and are comparably distributed, as they do not associate with melanoma risk per se, but with the increased susceptibility to MPM, given the prior SPM diagnosis. This may explain attenuated significance of rs2071304 associations in comparison of MPM cases with healthy controls. Hence, the MPM-specific alleles may be fully identifiable only by comparison of MPM versus SPM patients, an analysis in which the observed associations reached the strongest significance in our study. The minor allele of rs2071304 was shown to correlate with decreased expression of SPI1 in the MuTHER dataset (Fig. 1). SPI1 has previously been shown to be involved in the development of several different types of immune lineage precursor cells, including, T-cells, B-cells, dendritic cells (DC) and monocytes [41][42][43][44] . Lowered levels of PU.1, the protein product of SPI1, have also been shown to result in preferential development of B1 B cells 45 , which are involved in innate immunity and are often autoreactive 46 , and other studies suggest that decreased levels of PU.1 associate with autoimmune conditions 47 . Given this evidence and the results of this study, the patients with the "low-expressing" SPI1 allele could potentially have increased sensitivity to self-antigens, what may allow for improved clearing of "micro" melanomas before they developed into detectable MPMs 48 . In addition to its role in the immune system, SPI1 has been shown to upregulate Mcl-1 transcription in melanoma cells 49 and this upregulation prevents melanoma cells from undergoing endoplasmic reticulum stress-induced apoptosis 50 . As rs2071304 is strongly associated with expression of SPI1 in skin cells (p < 1E-15) in the MuTHER data, with the same directionality as in LCLs, this is another potential mechanism explaining how the alternate allele of rs2071304 is protective against MPM.
Lastly, we tested if ieQTLs were able to predict survival in MPM patients. The most significant association with OS was found for rs6695772, showing that the minor allele (C) associates with significantly worse survival in MPM patients (additive model: HR = 3.42; 95% CI = 1.57-7.42; p = 0.0019, dominant model: HR = 18.69; 95% CI 3.34-104.55; p = 0.0009), although the relatively wide confidence interval suggests that this result needs to be validated in a cohort with larger statistical power. Interestingly, this association with OS was reported in our recent study in a general melanoma population 25 , in which we suggested that rs6695772 is a putative inherited risk marker predictive of melanoma survival. While only borderline association was observed in SPM patients in this study, the pooled analysis of SPM and MPM patients showed association effect of rs6695772 with survival, comparable with MPM analysis alone. Hence, while further validation in MPM and SPM patients will be needed to confirm MPM specific association of this variant with survival, the findings from the current study suggest that the effect of rs6695772 on survival is more pronounced in MPM patients. Noteworthy, the association observed in univariate analysis remained comparably significant, strongly suggesting that rs6695772 associates with OS independently of other clinical and prognostic predictors. We have also attempted to test this association with melanoma-specific survival. However, this comparison was underpowered due to the lack of melanoma-specific death information for most of the patients. Nevertheless, by restricting the analysis to melanoma-specific survival we have found that rs6695772 was still the most significant association with MPM, albeit with reduced significance level (HR = 11.27; 95% CI = 2.05-61.90; p = 0.0053). As we discussed extensively in our prior report 25 , this ieQTL associates with decreased expression of BATF3 in LCLs in MuTHER dataset (Fig. 1). It has been documented that the loss of BATF3 decreases the ability of dendritic cells in presenting cell-associated antigens via the MHC-I complex, thereby impairing baseline antitumor response 51 . As such the effect of "low expressing allele" at this locus associated with MPM survival may have significant implications for impaired immunity affecting melanoma progression. Given the findings in this study and our prior observations 25 , this region should be examined further in larger studies for its prognostic value in improving the current clinical outcome assessments in melanoma progression.
In summary, our study has identified several ieQTLs that associate with the development of multiple primary melanomas following the single primary diagnosis. In particular, rs2071304 showed the strongest association with MPMs compared to patients with SPM. The results of our study suggest that the potential modifications to www.nature.com/scientificreports www.nature.com/scientificreports/ the host immune response competency via an interplay of germline genetic factors, may be a crucial trigger in altering the likelihood of patients with primary melanomas developing subsequent primary tumors. Given that patients who develop additional primary melanomas are at a greater 10-year mortality risk compared to those with single primary tumors 3 , the genetic markers predictive of MPM development reported here, if validated, can potentially provide new opportunities for improving screening and clinical management of a high-risk population in which early identification has been challenging. While requiring further validation in large collaborative efforts, our results suggest that besides other pathological and clinical factors of progression from single primary tumors to multiple primary melanomas (tumor genetics, microenvironment, baseline immunity, etc.), additional research should also strongly consider germline genetic underpinning. As we emphasize, addition of germline genetic biomarkers identified in such efforts could have substantial clinical benefit for MPM high-risk patients as well as those in general melanoma population.