The Expression Quantitative Trait Loci in Immune Pathways and their Effect on Cutaneous Melanoma Prognosis

While cutaneous melanoma represents only 4% of skin cancers, it accounts for about 80% of skin cancer–related deaths. Melanoma incidence rates have been rising 1.4% on average each year over the last 10 years, while mortality rates have remained steadily high, in particular for advanced stages. In 2015, approximately 73,870 new cases of melanoma of the skin and approximately 9,940 melanoma-related deaths are estimated to occur in the United States (1), reflecting traditionally poor disease outcomes associated with more advanced stages; the 5-year melanoma survival rates for stages I–II, III, and IV are 98.3, 63.0 and 16.6%, respectively (2). The negative trends in melanoma mortality are largely attributed to difficulties in clinical prognostication especially for more advanced stages, suggesting that in addition to standard clinical predictors, there are other factors affecting the unpredictability of melanoma outcomes. Melanoma is considered to be highly immunogenic with the ability to induce an immune response that can suppress tumor growth, a phenomenon which is believed to be governed by effector T cells. Observed tumor immunogenicity modulates prognosis of cutaneous melanoma and varies greatly on the individual level (3, 4), suggesting that different capacities of the immune system control tumor growth (5, 6). The germline genetic factors emerge as possible novel, personalized markers of cancer outcomes, including melanoma (7–11), some exhibiting putative immunoregulatory capabilities with tumor impact and therefore representing plausible modulators of observed individual immune-response heterogeneity (12). In general, the germline associations with clinical outcome of complex disease traits found in small candidate studies are difficult to interpret and often need validation in larger cohorts. In addition, despite substantial efforts, evidence supporting the biologic relevance of associated germline variants remains elusive, as they map almost exclusively in noncoding, often intergenic, regions. Findings from genome-wide association studies (GWAS) estimate that 88% of disease/trait-associated germline variants are noncoding, and 12% and 34% of them overlap with transcription factor–binding regions and DNase I hypersensitive sites (i.e., markers of DNA regulatory region), respectively, thus putatively impacting gene expression in cis- or trans- configuration (13, 14). In efforts to help interpret functional consequences of germline variation, several recent studies have genome-wide mapped variants that correlate with expression levels of nearby genes known as cis expression quantitative trait loci (cis-eQTL) in several different cell types (15–19). These large analyses generated comprehensive maps of inherited genetic variation that regulate gene expression (20, 21), thus representing plausible biologic candidates for association with human traits, including common diseases. While most eQTL studies were conducted on a relatively small sample size, a recent study by Grundberg and colleagues (2013) identified cis-eQTL SNPs in three selected tissues, including lymphoblastoid cell lines (LCL) derived from a large population of 857 well-phenotyped healthy female twins of the MuTHER (Multiple Tissue Human Expression Resource) project. Besides a large sample size which gives MuTHER resource enough statistical power to detect genetic variants with meaningful eQTL properties, the advantage of the MuTHER project's twin design allows confirmation of identified eQTLs separately in each twin set (18, 19).

Capitalizing on the MuTHER resources, the aim of this study was to investigate whether biologically relevant germline polymorphisms that regulate expression levels of immune-relevant genes in cells of the immune system, (e.g., LCLs) might serve as prognostic markers of melanoma clinical outcomes. By interrogating 382 immunomodulatory genes against eQTL data from MuTHER, we have evaluated 40 expression-regulating polymorphisms and their cumulative effects for association with melanoma clinical outcomes in a large population sample of 1,221 patients with cutaneous melanoma.

This study comprises a total of 1,221 patients with cutaneous melanoma (stage I to III) of self-reported European descent who were treated at the New York University Langone Medical Center (NYUMC). Blood samples, demographic, and clinical information including age at diagnosis, gender, self-reported family history of melanoma, primary tumor characteristics: anatomic site, thickness, histologic type, 2009 AJCC stage at diagnosis, ulceration status, as well as follow up information, were obtained following criteria established by the Interdisciplinary Melanoma Cooperative Group (IMCG; ref. 12, 22, 23). All patients gave written informed consent at the time of enrollment and the study was approved by the Institutional Review Board of the NYUMC.

Candidate immunomodulatory genes were selected from exploring the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using the following search terms: somatic diversification of T-cell receptor genes, T-cell receptor V(D)J recombination, cytotoxic T-cell differentiation, T-helper cell differentiation, regulatory T-cell differentiation, T-cell costimulation, T-cell receptor complex, T-cell receptor signaling pathway, dendritic cell differentiation, cytokines-immune, related to T-cell receptor, interleukin, cytokines-immuno, T-cell cytokine production, negative regulation of regulatory T-cell differentiation, regulation of T-cell cytokine production. A total of 382 genes were selected (Supplementary Table S1).

SNPs for genotyping were selected on the basis of information from sequence-based gene expression variations in lymphoblastoid cell lines (LCL), obtained from the recently established MuTHER project, which has been extensively described elsewhere (18, 19, 24, 25). The MuTHER project generated genomic and transcriptomic data from three disease-relevant tissues, including LCLs (i.e., representing cells of the immune system), which were derived from a cohort of 856 female Caucasian twins with detailed phenotypic information from the UK Adult Twin registry (18). For the purpose of this study, a publicly available list of top cis-eQTLs per probe in LCLs was mined for all the probes representing our panel of 382 immunomodulatory genes (Supplementary Table S1). A total of 50 SNPs with most significant cis-eQTL activity (ranking with P < 4.46 × 10⁻⁸) in cells of the immune system were selected for genotyping (Supplementary Table S2). To confirm genotype-expression associations for the selected 50 probe–SNP pairs, we obtained publicly available expression data from the ArrayExpress (accession no. E-TABM-1140), while access to the genotype dataset was obtained from the Department of Twin Research (DTR), King's College London. Twins (339 twin-pairs) from the same pair were separated into two twin sets and independent eQTL analyses were performed for each twin set using Spearman Rank Correlation as previously described (19). Genotype–expression correlations were assessed in 777 participants (including 339 twin-pairs) under three genetic models of inheritance [i.e., genotypic (three genotypes were coded 1,2,3), dominant (genotypes were coded 1,2,2) and recessive (genotypes were coded 1,1,2)] using Spearman rank correlation test.

For genotyping, genomic DNA was isolated from whole blood samples using a QiaAmp kit (Qiagen). All SNPs were genotyped using MassARRAY System (Agena Bioscience Inc.) according to the manufacturer's protocol, as described in detail elsewhere (12, 23). Primer design was not successful for eight SNPs, due to highly polymorphic regions around the SNPs of interest.

Expression levels of IL19 in CD4⁺ T cells were assessed using NanoString nCounter platform. CD4+ T cells were isolated from peripheral blood mononuclear cells of a subset of 43 patients with cutaneous melanoma as described previously (12), and RNA was extracted from purified CD4⁺ cells using RNeasy RNA Isolation kit (Qiagen). Probes were designed and synthesized by NanoString nCounter technologies to probe the following sequences of target transcripts: CCACAGACATGCACCATATAGAAGAGAGTTTCCAAGAAATCAAAAGAGCCATCCAAGCTAAGGACACCTTCCCAAATGTCACTATCCTGTCCACATTGGA for IL19 (NM_013371.3), GCAAGAAGTATGCTGAGGCTGTCACTCGGGCTAAGCAGATTGTGTGGAATGGTCCTGTGGGGGTATTTGAATGGGAAGCTTTTGCCCGGGGAACCAAAGC for PGK1 (NM_000291.2) and CGGTCGTGATGTGGTCTGTGGCCAACGAGCCTGCGTCCCACCTAGAATCTGCTGGCTACTACTTGAAGATGGTGATCGCTCACACCAAATCCTTGGACCC for GUSB (NM_000181.1). NanoString nCounter analysis was performed using an input of 200 ng of total RNA from each sample and hybridization of RNA with Nanostring probes was processed according to the manufacturer's protocols (NanoString Technologies, Inc). The data were analyzed using the nCounter digital analyzer software (Version 2.5.34). Counts for target gene IL19 were subjected to a technical normalization considering the counts obtained for positive control probe sets, followed by a biologic normalization using the two housekeeping genes (PGK1 and GUSB) included in the CodeSet.

Univariate Cox proportional hazard models were used to assess associations between demographic and clinical characteristics and overall (OS) and recurrence-free survival (RFS). Time for OS and RFS was defined as the time from the date of diagnosis to the date of death (OS) and recurrence (RFS), or the date of last follow up. In addition, for the patients with no recorded recurrence event, RFS was defined from date of diagnosis to the date of death due to melanoma. The effects of individual SNPs on OS and RFS were accessed by multivariate Cox proportional hazard model adjusted for age at diagnosis (≤60 years/>60 years), gender (female/male), tumor stage (I/II/III), thickness (<1.0/1.0–2.0/2.01–4.0/>4.0), ulceration (absent/present), and histologic subtype (superficial-spreading/nodular/desmoplastic/acral-lentiginous/lentigo-maligna/other), as described previously (12). For each SNP, two genetic models (i.e., two genetic models that showed strongest expression-genotype correlations in the MuTHER data) were applied in Cox analysis, separately. To assess whether the observed associations are independent of clinical variables, the SNPs surpassing the correction for multiple tests were further analyzed using Cox models adjusted only for age and gender. SNP–SNP interaction analysis was performed by counting the number of putative unfavorable genotypes. To address the issue of multiple comparisons for Cox proportional hazard models, we used Benjamini-Hochberg procedure while restricting the false discovery rate to 0.05 (26) and assuming 80 hypotheses (40 SNPs tested under 2 genetic models of inheritance). All statistical analysis of genotype data and survival were performed using the “survival” package in R.

The study included 1,221 patients with stage I–III cutaneous melanoma recruited by the IMCG at the NYUMC. All patients were of self-reported European descent with median age of 58.1 ± 16.5 years. The median time between diagnosis and follow-up was of 52.7 months. The 5-year survival rate was 84.1% for OS and 78.6% for RFS. Demographics as well as tumor characteristics are described in Table 1. Almost 85% of all primary diagnoses were of stage I and II. The majority of tumors (75.4%) were less than 2.0 mm thick and 80.6% of all primary tumors were not ulcerated. Superficial spreading melanoma was the most common histologic subtype of melanoma (57.2%), followed by nodular melanoma (25.6%). Using Cox regression analysis, primary tumor stage, thickness, ulceration status, and histologic subtype were found significantly associated with OS and RFS (all P < 2.2 × 10⁻¹⁶), respectively. Patients' age at diagnosis (P = 5.0 × 10⁻⁸) and family history of melanoma (P < 0.001) were significantly associated with OS (Table 1).

Using MuTHER resources, a total of 398 probe–SNP pairs (i.e., eQTLs) corresponding to 268 unique immunomodulatory genes and 385 unique SNPs were identified at the P < 0.05 significance level (Supplementary Table S1). The top most significant 50 unique SNPs (ranking with P < 4.46 × 10⁻⁸) affected expression levels of 54 unique expression probes, which represent 48 unique immune-relevant genes (Supplementary Table S2). For these eQTLs we further computed expression–genotype correlations under three genetic models (genotypic, dominant, and recessive) using Spearman correlation test to identify underlying mode of inheritance of eQTL effect. Genotype–expression correlations were strongest under genotypic and/or dominant genetic model for most SNP (Table 2). To further confirm the genotype–expression associations for 50 eQTL SNPs, we separated twins into two twin sets and performed separate eQTL analyses for SNP-probe pairs of interest in each twin sets using Spearman Rank Correlation as described previously (19). All examined probe-SNP pairs were significant in both twin sets at a significance level of P < 0.05 and the direction of association with expression (i.e., negative vs. positive) remained unchanged (Table 2). The two genetic models that best correlated with expression were considered in downstream survival analysis for each eQTL.

The most significant 50 eQTL SNPs were selected for genotyping in a population cohort of 1,221 patients with melanoma. Only 40 SNPs were successfully genotyped and passed to association analyses; primer design was not successful for 8 SNPs and 2 SNPs did not pass our genotyping quality control filters (Hardy–Weinberg equilibrium P < 0.001). The associations between each individual SNP and cutaneous melanoma clinical outcomes of 1,221 patients with melanoma (OS and RFS) were evaluated using Cox proportional hazards models as detailed in Materials and Methods. The results of the single SNP analysis are summarized in Table 3. In the RFS analysis, the most significant association was observed for rs9921791, an eQTL controlling expression of MLST8 in LCLs (Fig. 1A) and ranking #16 among the top 50 eQTLs in the study (Supplementary Table S2), where carriers of at least one copy of minor T allele were associated with better outcome under the dominant model (HR, 0.52; 95% CI, 0.35–0.79; P = 0.0009). While other variants show nominal significance under different genetic models of analysis, including rs6695772 (IL19, P = 0.005), rs6695772 (BATF3, P = 0.006), rs841718 (STAT6, P = 0.015), rs11539345 (CD40, P = 0.016), and rs2276645 (ZAP70, P = 0.048), after adjusting for multiple testing, none of the variants remained statistically significant in RFS analyses.

For OS, six SNPs showed significant associations using multivariate analysis (Table 3). The most significant association with OS was observed for rs6673928 [an eQTL impacting expression of IL19 (Fig. 1B) and ranking #12 among the top 50 eQTLs in the study (Supplementary Table S2)], under the dominant model, in which the variant T allele was associated with improved OS (HR, 0.56; 95% CI, 0.41–0.77, P = 0.0002). This association remained significant after adjusting for multiple testing. Moreover, comparably significant association between OS and rs6673928 was also observed in an analysis adjusted by only age and gender (see Materials and Methods), suggesting that survival effect of rs6673928 is independent of other AJCC clinical variables (HR, 0.61; 95% CI, 0.45–0.82; P = 0.0008). Our second most significant OS association, reaching the level of multiple testing adjusted significance, was observed for rs6695772 (influencing BATF3 expression levels, Fig. 1C). The carriers of the rs6695772 minor C allele were associated with worse survival under the dominant model (HR, 1.64; 95% CI; 1.19–2.24; P = 0.0019). Polymorphisms rs6673928 and rs6695772, both exhibiting most significant association with OS, are located approximately 6 Mb apart on chromosome 1q32, previously suggested for association with melanoma outcome (12). Linkage disequilibrium (LD) analysis showed that the two SNPs are independent (r² = 0.002; D′ = 0.099), which we further confirmed by performing association between OS and rs6695772 using multivariate Cox regression analysis while also adjusting for the rs6673928 variant and found no difference in effect size and significance. Moreover, as per MuTHER data, there is no correlation between rs6673928 (eQTL for IL19) and BATF3 expression and similarly no association is observed between rs6695772 (eQTL for BATF3) and expression levels of IL19. This further suggests that both loci represent two independent eQTL effects (Supplementary Fig. S1). Other associations with OS, not reaching the adjustments for multiple testing, were observed for rs841718 (STAT6, P_{dominant model} = 0.007), rs2701652 (IRAK3, P = 0.035), rs7720838 (PTGER4, P_{recessive model} = 0.037), and rs11569345 (CD40, P_{dominant model} = 0.040). Of note, SNPs rs11569345 (CD40), rs6673928 (IL19), rs6695772 (BATF3), and rs841718 (STAT6) exhibited association with both RFS and OS, and the directionality of associations was comparable between RFS and OS for all analyzed SNPs (Table 3).

We further assessed the cumulative effect of two eQTL SNPs, rs6673928 (IL19), and rs6695772 (BATF3) on OS, as these were the most significant associations in the single SNP analysis after adjustment for multiple testing, and, interestingly, both variants, albeit genetically independent (not in LD), map in the same locus at 1q32. SNP–SNP interaction analysis was performed by counting the number of putative, unfavorable genotypes (i.e., associated with worse outcome) of both SNPs and assessing their association with OS (multivariate Cox model is presented in Supplementary Table S3). The following genotypes were considered as unfavorable: rs6673928 (wild-type) and rs6695772 (heterozygotes and variant homozygotes). Comparing with the reference subjects carrying 0 or 1 unfavorable genotypes, this analysis shows that the subjects carrying 2 unfavorable genotypes had significantly worse survival (HR, 1.92; 95% CI, 1.43–2.60; P = 1.87 × 10⁻⁵). Five-year OS rate was 86.3% for reference subjects, while 5-year survival rate decreased to 80.5% among carriers of 2 unfavorable genotypes (Fig. 2).

To further validate the genotype–expression correlation of our most significant eQTL in melanoma survival analysis, we tested the correlation of rs6679328 with IL19 expression in CD4⁺ T cells purified from a subset of 43 patients with melanoma from a population genotyped in this study. Following isolation of CD4⁺ T cells the expression level of IL19 was analyzed from total RNA. Expression levels of IL19 were compared between GG and GT genotypes for rs667328 (no homozygotes for minor T allele were observed among 43 patients). Using Spearman correlation approach, we were able to observe a positive association between T allele and expression level of IL19 (Spearman coefficient ρ = 0.32, P = 0.0347; Fig. 3), further confirming that the eQTL from UK twin studies has a true effect on IL19 expression in an independent subset of patients with melanoma tested in this study.

The discovery of biologically impactful and personalized melanoma prognostic markers complementing currently established, yet clinically less specific, histopathologic indicators is one of the key objectives of current melanoma research. Recent studies suggest that germline genetic variants may modulate cutaneous melanoma clinical endpoints, thereby representing potentially personalized and easily accessible prognostic biomarkers (7–12). However, while the majority of currently published studies on germline associations with prognosis still require independent validation in larger cohorts, the biologic and functional uncertainty of most prognostic variants further limits their clinical consideration. With the exception of few melanoma etiology related SNPs [e.g., MC1R (9), vitamin D-binding protein (11)], in general, almost exclusively, the genetic variants identified to date for associations with melanoma risk (i.e., GWAS loci) or prognosis (reviewed in ref. 8), map to noncoding regions with unknown biologic impact. In fact, from a broader point of view the general lack of knowledge on functional consequences of germline genetic variation associated with human diseases is apparent and has been a major hurdle in translating these findings into clinical practice. As such, novel strategies for the discovery of biologically impactful germline genetic surrogates associated with cancer risk (e.g., in GWAS data), prognosis, or therapy response (27–29) are highly demanded not only in regard to more predictive strategies, but importantly for improved biologic understanding ultimately leading to more targeted treatments. Considering the fact that melanomas are highly immunogenic, the genetic variants in immunomodulatory genes that are significantly and reproducibly associated with gene expression in immune cells as validated eQTLs may have a strong potential to serve as clinically actionable prognostic markers.

By exploring recent data of genome-wide eQTLs generated on 777 healthy female twins by MuTHER project (18), in this study, we have tested the hypothesis that the individual genetic variation associated with the expression of immune-related genes plays a role in modulating cutaneous melanoma outcomes, as immune response appears to be one of the key mechanisms for controlling melanoma progression. Our study provides a first in-depth analysis aimed at mapping the functionally important germline variants in immunomodulatory networks as potential markers of melanoma prognosis.

Using the unique approach of interrogating the MuTHER eQTL data, we have identified 385 SNPs significantly associated with the expression of 268 immunoregulatory genes (Supplementary Table S1). For the most significant 40 variants (P < 4.46 × 10⁻⁸; Supplementary Table S2) with comparable genotype–expression correlations in both matched healthy female twin sets, we assessed genotype–expression correlations under three genetic models (i.e., genotypic, dominant, or recessive) using Spearman rank correlation test. The best two modes of inheritance for each eQTL were examined for their association with cutaneous melanoma clinical outcomes in 1,221 patients with melanoma. Our approach for the first time identified several eQTL variants that were significantly associated with melanoma survival.

The most significant finding in our study is the association of melanoma overall survival (OS) with eQTL variant rs6673928 at 1q32.1, impacting the expression of IL19 gene [linear regression coefficient (β) = 0.12, P = 5.66 × 10⁻²³]. We observed a strong association of rs6673928 with better survival for the carriers of a minor T allele (HR, 0.56; 95% CI, 0.41–0.77; P = 0.0002), which correlated with increased expression of IL19 in LCLs from MuTHER dataset, and this correlation was also validated in CD4⁺ T cells purified from a subset of patients with melanoma from our study population (Fig. 3). To our knowledge, this is the first evidence suggesting that increased germline expression of IL19 in cells of an individual's immune system associates with better clinical outcome possibly via suppression of melanoma progression (Fig. 1B). These findings are particularly interesting in regards to sparse experimental data on IL19 function. While the putative involvement of IL19 has been reported in a range of diseases including cancer (30) or autoimmune disorders (31, 32), in these and other studies, IL19 was shown to have both suppressive and stimulatory capacities in immune regulation (33). IL19 is a member of IL10 family of cytokines and notably, both genes map within a relatively narrow region (∼26 kb apart) at 1p.32 locus. Interestingly, eQTL data from MuTHER project show that our most significant variant rs6673928 also associates with expression of IL10, albeit with less significance (P < 2.5 × 10⁻⁶) compared to the effect on IL19 expression. This is interesting, as we have recently identified strong association with cutaneous melanoma survival for a variant near IL10, rs3024493 (12). In that recent report, we showed that the association with improved OS is driven by rs3024493 heterozygotes, which secrete medium levels of IL10, as compared with low-secreting minor allele homozygotes conversely associated with worse outcome, which is consistent with directionality of the effects for IL19 in the current study. Other prior smaller scale studies also reported associations at 1q32.1 with melanoma survival for a set of three highly correlated polymorphisms in IL10 promoter: rs1800896, rs1800871, and rs1800872 (34–37). Similarly to our observations, these studies reported that the IL10 high-level expression genotypes are protective in cutaneous melanoma, while low-level expression genotypes associate with poorer disease prognosis. These consistent reports clearly suggest that the genetic variation at 1q32 may result in specific gene-expression patterns regulating several interleukin candidates in the locus with an impact on melanoma clinical outcome, likely via modulation of melanoma immune surveillance. Notably, all the variants associated with cutaneous melanoma outcomes at 1q32.1 in the current and prior reports are scattered in a relatively narrow region spanning ∼9.6 kbp, thus raising a possibility that the variants are in LD. Using the data from 1000 Genomes Pilot project, we found that rs1800896 (previous studies; ref. 35), rs3024493 [our recent study; ref. 12), and rs6673928 (current study) show little to no correlation (r²_{rs1800896-rs3024493} = 0.31, r²_{rs1800896-rs6673928} = 0.24, and r²_{rs6673928-rs3024493} = 0.07), which was also confirmed in our study population (data not shown). This indicates that the multivariant associations with cutaneous melanoma survival at 1q32.1 may be due to other mechanisms. As discussed in our recent report (12), the region of 1q32.1 appears to be in an extensive transcriptional “hot spot.” The region of associated variants shows the presence of several strong DNase I hypersensitive sites in T cells, involving multiple transcription factors spread across a 30 kbp region (Supplementary Fig. S2). It is therefore possible that 1q32.1 exerts a broader positional effect in the immune cells, affecting the expression of several genes simultaneously in this locus. Given the amount of published data showing significant expression correlation of different interleukin gene targets at 1q32.1 in immune regulation (32, 38), it is likely that genes of this locus share common gene expression–regulatory elements. The genetic variants in these regulatory elements would then produce specific expression signatures, impacting both physiologic immune response as well as disease outcome. This scenario is supported by other intriguing findings generated in our study. We found another variant on chromosome 1, rs6695772, which was associated with OS. The carriers of at least one copy of the minor C allele for rs6695772 had worse prognosis when compared to wild-type GG homozygotes (HR, 1.64; 95% CI, 1.19–2.24; P = 0.0019). The MuTHER eQTL data in LCLs associate minor C allele of rs6695772 with decreased expression levels of BATF3 gene (basic leucine zipper transcription factor, ATF-like 3) in a dose-dependent manner [linear regression coefficient (β) = −0.16, P = 6.93 × 10⁻¹⁰; Fig. 1C). BATF3 is a positive immune regulator, primarily involved in stimulation of CD8α⁺ dendritic cells (39). Our data appear to be consistent with these findings, suggesting that the decreased expression of BATF3 in immune cells predicts worse cutaneous melanoma outcome, likely due to suppressed immune surveillance of tumor progression (Fig. 1C). Interestingly, rs6695772 maps ∼6 Mb downstream from our most significantly associated variant, rs6673928, and although not in LD (r² = 0.002), the comparably significant association effects observed with melanoma survival and relatively close proximity of both loci might suggest common genetic or biologic underpinnings. Notably, overexpression of BATF3 in T cells has also been shown to stimulate Th17 cell differentiation. (40). Moreover, knocking down BATF3 in Th2 cells dramatically decreased expression of IL4 and IL10 cytokines (41, 42). This is an intriguing biologic connection as both the overexpression of IL10 and elevated levels of Th17 cells in peripheral blood were correlated with improved cancer patient survival in previous studies (43, 44). These findings therefore align with our observations that minor T allele carriers of germline eQTL rs6695772, associated with worse OS, express low levels of BATF3, which may in turn downregulate IL4 and IL10. In an attempt to support such a hypothesis, we analyzed the correlation between BATF3 expression and expression levels of IL4, IL10, IL17, and IL19, in the MuTHER data (Supplementary Fig. S3). Although we did not note the correlation of BATF3 expression with expression of IL10, nor have we confirmed positive correlation of BATF3 with IL4, it is possible that the expression correlation is T-helper subtype specific. However, we noted a significant correlation between expression of BATF3 and IL19, suggesting a possibility that there is a novel, previously unexplored molecular interaction of both proteins in immune surveillance of melanoma progression. Nevertheless, the biologic meaning and causality of these associations need to be further investigated, likely using in vivo models.

Because of functional commonalities and putative mutual interaction between the two loci (rs6673928 and rs6695772) most significantly associated with melanoma OS in our study, we explored possible cumulative effects of these two variants. We found that the joint effect of both variants on OS is substantially stronger (HR, 1.92; 95% CI, 1.43–2.60; P = 1.87 × 10⁻⁵; Fig. 2), when compared with single SNP analysis. Interestingly, the analysis adjusted for only age and gender shows similarly strong association effect (HR, 1.60; 95% CI, 1.21–2.11; P = 9.26 × 10⁻⁴) when compared with multivariate analysis adjusted for age, gender, and other established clinical markers. This clearly suggests that the observed joint effects are possibly independent from histopathologic predictors. These observations for the first time propose the germline genetic variants as independent prognostic factors and, due to the strength of their joint interaction, the clinically actionable personalized biomarkers of melanoma outcomes. However, the final replication of these associations in additional subsets of melanomas will be important to provide a definitive verdict on their consideration for a clinically valid prognostic test.

Significant association was also noted for rs9921791, located nearby MLST8 (mammalian lethal with SEC13 protein 8). The carriers of minor T allele (associated with increased expression of MLST8) recurred significantly later compared with patients with CC genotypes (Fig. 1A). MLST8 is a regulator of mTOR kinase activity (45, 46) and interestingly, mTOR signaling was recently shown to promote T-cell development (47–49). This suggests that upregulation of MLST8 expression among rs99212791-T allele carriers with cutaneous melanoma, may stimulate immune tumoral response via mTOR pathway-mediated T-cell activation, leading to suppression of cutaneous melanoma recurrence. However, despite these encouraging findings, exploring such hypothesis is premature pending the association validation and more detailed investigation into the underlying molecular role of MLST8 in melanoma clinical outcomes.

It is important to mention that the genotype–expression correlations tested in our study were based exclusively on the data from female-only MuTHER cohort. Although all our analyses were adjusted for gender, the gender-stratified tests for our most significant associations revealed a gender-specific effect for some variants, in particular rs6695772 (BATF3) association with OS observed only in males. (Supplementary Table S4). While this could likely be attributed to reduced statistical power (the samples size in each separate gender-stratified analysis was reduced approximately by half), it is also possible that the gender-specific associations are due to yet-unknown biologic underpinnings. As it is difficult to draw reliable conclusions at this stage, to confirm their potential biologic meaning, the gender-specific testing of rs6695772 (BATF3) in a larger melanoma population will be needed in subsequent efforts.

As our current report is hypothesis driven and focuses on cis-eQTL associations in relatively narrow selection of immune genes due to their biologically plausible role in melanoma, we did not assess trans-acting genotype–expression correlations on the genome-wide scale that may also be important. While this may be a potential limitation in our design, the recent MuTHER study on mapping trans-eQTLs has concluded that direct cis effects on local genes are stronger than indirect trans effects (21). Moreover, in that recent report, only a handful of trans effects at a 5% false discovery rate have been identified on genome-wide scale. This is likely attributed to the power limitations of available eQTL resources, as the interrogation of trans-acting genotype correlations involves much larger number of tests resulting in more rigorous control for false discovery rate. However, it is estimated that 65% of gene expression heritability is trans-regulated and about 52% of cis-eQTLs also have trans-acting effect (21), strongly suggesting importance of inherited trans-eQTLs as potentially important prognostic biomarkers. Therefore, with the expansion of eQTL resources, as part of ongoing and future efforts, the comprehensive assessment of trans-eQTLs, using the similar approach applied in our study, will become feasible for the identification of novel clinically relevant outcome modulators of melanoma and other cancers.

In conclusion, our unique approach of interrogating lymphocyte-specific eQTLs from healthy twins was notably powerful in identifying several immunomodulatory eQTLs, and indirectly, their gene targets, including IL19 and BATF3 at 1q32, as novel biologically relevant predictors of cutaneous melanoma prognosis. In addition, the substantially enhanced cumulative effect of these associations strongly encourages the consideration of joint screening of these variants in a prognostic clinically relevant test in the near future. Our study suggests that the eQTL-based strategy proposed here will be highly efficient in discovering novel molecular markers of outcome, risk, or therapy response in other human cancers driven by specific molecular pathways.

M. Krogsgaard reports receiving speakers bureau honoraria from LifeSci Corp and is a consultant/advisory board member for Agenus. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M. Vogelsang, J. Rendleman, A. Romanchuk, M. Krogsgaard, T. Kirchhoff

Development of methodology: M. Vogelsang, A. Romanchuk, T. Kirchhoff

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Rendleman, A. Romanchuk, R.S. Berman, M. Krogsgaard, I. Osman

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Vogelsang, C.N. Martinez, J. Rendleman, A.B. Bapodra, K. Malecek, A. Romanchuk, T. Kirchhoff

Writing, review, and/or revision of the manuscript: M. Vogelsang, C.N. Martinez, J. Rendleman, A.B. Bapodra, A. Romanchuk, E. Kazlow, R.L. Shapiro, M. Krogsgaard, I. Osman, T. Kirchhoff

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): K. Malecek, A. Romanchuk

Study supervision: R.S. Berman, M. Krogsgaard, T. Kirchhoff

The study was funded by the grants from NCI 1R21CA184924-01 (to T. Kirchhoff), 1R01CA187060-01A1 (to T. Kirchhoff) and Cancer Center Support Grant P30CA016087. The TwinsUK study was funded by the Wellcome Trust; European Community's Seventh Framework Programme (FP7/2007-2013). The study also receives support from the National Institute for Health Research (NIHR)-funded BioResource, Clinical Research Facility and Biomedical Research Centre based at Guy's and St Thomas' NHS Foundation Trust in partnership with King's College London. SNP Genotyping was performed by The Wellcome Trust Sanger Institute and National Eye Institute via NIH/CIDR.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.