Abstract
1α,25-Dihydroxyvitamin D3 signals via the vitamin D receptor (VDR). Higher serum vitamin D is associated with thinner primary melanoma and better outcome, although a causal mechanism has not been established. As patients with melanoma commonly avoid sun exposure, and consequent vitamin D deficiency might worsen outcomes, we interrogated 703 primary melanoma transcriptomes to understand the role of vitamin D–VDR signaling and replicated the findings in The Cancer Genome Atlas metastases. VDR expression was independently protective for melanoma-related death in both primary and metastatic disease. High tumor VDR expression was associated with upregulation of pathways mediating antitumor immunity and corresponding with higher imputed immune cell scores and histologically detected tumor-infiltrating lymphocytes. High VDR–expressing tumors had downregulation of proliferative pathways, notably Wnt/β-catenin signaling. Deleterious low VDR levels resulted from promoter methylation and gene deletion in metastases. Vitamin D deficiency (<25 nmol/L ∼ 10 ng/mL) shortened survival in primary melanoma in a VDR-dependent manner. In vitro functional validation studies showed that elevated vitamin D–VDR signaling inhibited Wnt/β-catenin signaling genes. Murine melanoma cells overexpressing VDR produced fewer pulmonary metastases than controls in tail-vein metastasis assays. In summary, vitamin D–VDR signaling contributes to controlling pro-proliferative/immunosuppressive Wnt/β-catenin signaling in melanoma and this is associated with less metastatic disease and stronger host immune responses. This is evidence of a causal relationship between vitamin D–VDR signaling and melanoma survival, which should be explored as a therapeutic target in primary resistance to checkpoint blockade.
VDR expression could potentially be used as a biomarker to stratify patients with melanoma that may respond better to immunotherapy.
Introduction
1α,25-Dihydroxyvitamin D3 is the ligand for the dimeric vitamin D receptor (VDR) and retinoid X receptor (RXR): ligand–receptor binding facilitates transcription of target genes containing the vitamin D response element (1). The physiologic effect of the vitamin D–VDR signaling axis is often target tissue specific (2).
The association of low serum 25-hydroxyvitamin D2/3 levels (henceforth referred to as vitamin D) with higher cancer incidence has been reported (3), but the significance of the association has been debated (4). Extensive in vitro evidence however indicates an antiproliferative role of vitamin D, with 1,25(OH)2 vitamin D3 treatment shown to induce expression of proapoptotic genes and antiproliferative genes in prostate (5), breast (6), colon (7), squamous cell carcinoma, and leukaemia cells (8, 9).
We have previously reported that higher serum vitamin D levels at recruitment were associated with lower American Joint Committee on Cancer (AJCC) stage and better melanoma-specific survival (MSS) in the Leeds Melanoma Cohort (LMC; ref. 10). This was subsequently verified in five additional studies (11–15), which collectively indicate a significant role for vitamin D–VDR signaling in melanoma progression.
The unique dataset used in this study was derived from 703 formalin-fixed, paraffin-embedded (FFPE) primary melanomas, from the LMC, a population-based and extensively annotated cohort with a long follow-up (10). Tumor-derived transcriptomic data, clinical, histopathologic, and whole-genome copy number alteration (CNA) data were jointly analyzed to assess the pan-genome effects of vitamin D–VDR signaling and to determine the processes most associated with this pathway. Importantly, these findings were replicated in The Cancer Genome Atlas (TCGA) metastatic melanomas and then functionally validated using in vitro and in vivo experiments.
Materials and Methods
The LMC tumor transcriptome
LMC tumor transcriptome processing was described previously (accession no. EGAS00001002922; ref. 16). Briefly, tumor samples were taken from FFPE primary melanomas and RNA was extracted to generate whole-genome gene expression data (Illumina DASL HT12.4 array). Background correction and quantile normalization were applied; singular value decomposition was used to assess the confounding factors that were subsequently adjusted out. Participants in the LMC gave written informed consent; the study was conducted in accordance to international ethical guidelines (Declaration of Helsinki) and was approved by the national ethics committee (MREC 1/03/57 and PIAG3-09(d)/2003).
Measurement of serum vitamin D at diagnosis
For 554 of the 703 participants with transcriptomic data, 25-hydroxyvitamin vitamin D2 and D3 (nmol/L) was measured as described previously (10) and adjusted for season (see Supplementary Materials and Methods).
TCGA melanoma dataset
TCGA metastatic melanoma data (n = 365 samples) such as transcriptomic (RNA-seq), clinical, methylation, and copy number data were downloaded from cBioPortal (http://www.cbioportal.org/). The same statistical tests and software/packages were used to analyze Leeds and TCGA data.
Statistical analyses
Association of VDR expression with clinical variables.
Univariate and multivariate linear regression were used to test the association between VDR and AJCC stage, mitotic number, tumor site, age, and sex.
Association of VDR expression with MSS.
Univariate and multivariate Cox proportional hazards models were used to test the prognostic effect of VDR expression level after adjusting for AJCC stage (7th Edition), tumor mitotic number, tumor site, and tumor-infiltrating lymphocytes (TIL). For this analysis VDR expression was on the continuous log2 scale, meaning that HR per VDR unit corresponds to change in hazard when expression is doubled.
VDR copy number changes in the LMC tumors
Copy number profiles were estimated from next-generation sequencing output from sequenced DNA samples from a subset of LMC tumors (n = 276, 39%) as described previously (17). Gistic2.0 was used to identify VDR copy number estimates (see Supplementary Materials and Methods).
Whole-transcriptome correlation with VDR
Linear regression analysis and Benjamini–Hochberg multiple testing correction (FDR) was used to test the correlation between the expression of each gene and VDR in the tumors. Genes with FDR < 0.05 and |Reg Coef| > 0.2 (regression coefficient; Reg Coef) were plotted in a volcano plot (R function: “plot”) and used for functional enrichment analysis.
Because VDR is expressed by keratinocytes (18), a FLG2-adjusted whole-transcriptome correlation with VDR was additionally performed [Fillagrin family member 2 (FLG2) being a maker of keratinocyte differentiation). This sensitivity analysis was conducted to account for any bias in VDR expression, which might have originated from keratin-rich melanoma subsets, which have been reported previously (19–21).
Functional enrichment analyses
The “Gene set/mutation analysis” feature of ReactomeFIViz (22) was used to identify significantly enriched pathways (Benjamini–Hochberg FDR < 0.05, from hypergeometric test) from a given input gene list.
VDR-binding regions
The genomic regions identified as having VDR-binding peaks across six tissue types (23) were downloaded as BED files. In addition, genomic regions known to contain the VDR-binding motif were downloaded from Motifmap (as BED files). In both cases, genes associated with genomic regions were identified using GREAT 3.0.0 (24): “basal plus extension” approach was used with Human GRCh37 assembly, whole genome as background regions, and the gene regulatory domain set to ±20 Kb upstream and ±400 Kb distal. The genes that mapped to these regions (“region–gene associations”) were exported and their overlap with VDR-correlated genes in the LMC (at FDR < 0.05) was assessed.
VDR expression across Lund and TCGA molecular phenotypes
Nearest centroid method (21) was used to classify the LMC tumors into TCGA and Lund molecular melanoma phenotypes (see Supplementary Materials and Methods). Differential VDR expression across these subgroups was tested with the Mann–Whitney test.
Imputed immune scores
As we described previously (25), briefly, 28 immune cell scores were calculated as mean expression of genes pertaining to an immune cell type, after deducting genes identified as potentially nonimmune cell specific, from the initial immune gene signature described by Angelova and colleagues (26). Correlation analysis of each immune cell score with VDR expression was conducted.
Expression of VDR and response to checkpoint blockade
VDR expression was compared between responders and nonresponders in two published albeit small transcriptomic data sets: (i) 38 patients treated with PD-1 blockade (27) and (ii) 40 patients treated with CTLA4 blockade (28). Both studies used pretreatment biopsies. Fold change was computed as ratio of mean VDR expression in responders to nonresponders.
Vitamin D-VDR subgroup analysis
X-tile (29) was used to identify patient subgroups with the most contrasted survival profiles (melanoma specific) based on their tumor VDR expression in the LMC. This approach was trained in randomly selected two-thirds of the samples and validated in the remaining third. The cut-off points defining patient subgroups in terms of percentiles were applied to the TCGA metastatic melanoma dataset for replication.
The identified VDR groups were further stratified on the basis of the participants' serum vitamin D levels at recruitment (season adjusted) in the LMC. Clinically defined, vitamin D deficiency is generally considered to be <25 nmol/L (30), with recent evidence defining vitamin D deficiency able to compromise bone health as ≤30 nmol/L (31). Thus, in the LMC dataset, serum vitamin D levels < 25 nmol/L and ≥25 nmol/L were classified as “deficient” and “sufficient.”
Vitamin D treatment of human melanoma cells
SK-MEL-28 and MeWo cells (courtesy of Professor Alan Melcher, Leeds CRUK cell line bank) were authenticated using the PowerPlex 16 System (Promega), cultured, treated with 1,25 (OH)2 vitamin D3 (Supplementary Materials and Methods), collected at different time points, RNA extracted, and used to generate gene expression profiles (HG-U133 plus 2.0 array, Affymetrix). Genes differentially expressed in vitamin D–treated versus control cells were identified at each timepoint (FDR < 0.10) and were used for functional enrichment using Reactome FIViz.
Generation of VDR-B16BL6 and control-B16BL6 cells
B16BL6 cells (chosen owing to lowest endogenous VDR expression among B16 strains, estimated from transcriptomic data by Dr. Martin del Castillo, Wellcome Sanger Institute, Hinxton, United Kingdom, personal communication) were purchased from the M.D. Anderson Cancer Center Cell Line Core facility. Cells were screened for the presence of Mycoplasma and mouse pathogens (at Charles River Laboratories) before culturing. Early passage B16BL6 cells were cultured and transfected with murine VDR cDNA (synthesized by GeneArt) or empty backbone plasmid to generate the VDR-B16BL6 (clones V1 and V2) and control-B16BL6 (clones C1 and C2) cells, respectively (Supplementary Materials and Methods), verified by Western blot analysis and qRT-PCR.
In vivo tail-vein metastasis assay
Although subcutaneous B16 models have been used to demonstrate the role of natural killer (NK) cells in host responses to melanoma, we adopted the tail-vein model based on successful demonstration by the group of a role for effector T cells and NK cells in B16 pulmonary metastases studies (32) and other reports of a role for T cells in B16 melanoma cells (33–35). The care and use of all mice in this study were in accordance with the United Kingdom Animals in Science Regulation Unit's Code of Practice for the Housing and Care of Animals Bred, or Used for Scientific Purposes, the Animals (Scientific Procedures) Act 1986 Amendment Regulations 2012, and all procedures were performed under a United Kingdom Home Office Project license, which was reviewed and approved by the Sanger Institute's Animal Welfare and Ethical Review Body. Housing and husbandry conditions were as described previously (32) and mice were maintained on Mouse Breeders' Diet (Laboratory Diets, 5021-3) throughout the study. V1, V2, C1, and C2 cells (detailed above) were tail-vein administered to 6- to 10-week-old sex-matched wild-type (C57BL/6NTac) mice (104 cells in 0.1 mL PBS). After 21 days, mice were humanely sacrificed and their lungs macroscopically examined to determine the number of metastatic deposits in all five lobes (“met count”). Lungs were formalin-fixed, processed, and embedded in paraffin wax blocks, from which consecutive 5-μm sections were cut and used for hematoxylin and eosin (H&E) and CD3 staining. H&E sections were digitally scanned (Leica Aperio AT2) and total metastatic area (“met area” μm2) was calculated as the total area of all metastatic deposits, across all five lung lobes using Aperio Imagescope (Leica Biosystems). CD3+ lymphocytes were estimated as described in Supplementary Materials and Methods.
Wnt/β-catenin signaling in VDR-B16BL6 and control-B16BL6 cells
cDNA from V1, V2, C1, and C2 cells were analyzed using a RT-PCR array of 84 mouse Wnt/β-catenin pathway genes (see Supplementary Materials and Methods). Relative expression was calculated using the ΔΔCt method, normalized to average Ct of the five housekeeping genes provided in the array. Fold change (FC) of the VDR-B16BL6 clones relative to control-B16BL6 clones was calculated as follows: {\rm F{C_{V1( {or} )V2}} = 2^{ - \Delta \Delta {C_{\rm{t}}}_{_{{\rm{V1(or)V2}}}}}}, where {\rm \Delta \Delta {C_{{{\rm{t}}_{{\rm{V1}}( {{\rm{or}}} ){\rm{V2}}}}}}{\rm{ }} = {\rm{ }}\Delta {C_{{{\rm{t}}_{{\rm{V1}}( {{\rm{or}}} ){\rm{V2}}}}{\rm{ }}}} - {\rm{ }}\Delta {C_{{{\rm{t}}_{{\rm{avg}}( {{\rm{C1 }} \,&\, {\rm{ C2}}} )}}}}.}
Results
VDR expression is independently protective for melanoma death in the LMC primary and TCGA metastatic melanomas
VDR expression was significantly lower in tumors of higher AJCC stage, higher mitotic rate, and tumors on the trunk and sun protected sites (compared with tumors on the head; Univariate analysis; Table 1) in the 703 LMC primaries. In multivariate analysis, lower VDR expression was independently associated with higher mitotic rate (P = 0.001) and tumor site (P = 0.001 for tumors of trunk compared with those on head), with borderline significance for higher stage (P < 0.06; Table 1).
. | Univariate . | Multivariate . | ||||
---|---|---|---|---|---|---|
Correlation of tumor VDR expression with . | Regression coefficient . | SE . | P . | Regression coefficient . | SE . | P . |
Age (per year) | −0.005 | 0.002 | 0.04 | −0.003 | 0.002 | 0.14 |
Sex | ||||||
Femalesa | ||||||
Males | −0.18 | 0.06 | 0.003 | −0.09 | 0.06 | 0.16 |
AJCC stage | ||||||
Stage Ia | ||||||
Stage II | −0.17 | 0.07 | 0.012 | −0.12 | 0.07 | 0.08 |
Stage III | −0.25 | 0.09 | 0.009 | −0.18 | 0.09 | 0.06 |
Mitotic rate | ||||||
<1 mitoses/mm2 tumora | ||||||
≥1 mitoses/mm2 tumor | −0.23 | 0.06 | 0.0004 | −0.20 | 0.06 | 0.001 |
Tumor site | ||||||
Heada | ||||||
Limbs | −0.04 | 0.10 | 0.63 | −0.11 | 0.10 | 0.28 |
Trunk | −0.35 | 0.10 | 0.001 | −0.36 | 0.10 | 0.001 |
Rare (sun-protected sites) | −0.44 | 0.12 | 0.001 | −0.38 | 0.13 | 0.003 |
. | Univariate . | Multivariate . | ||||
---|---|---|---|---|---|---|
Correlation of tumor VDR expression with . | Regression coefficient . | SE . | P . | Regression coefficient . | SE . | P . |
Age (per year) | −0.005 | 0.002 | 0.04 | −0.003 | 0.002 | 0.14 |
Sex | ||||||
Femalesa | ||||||
Males | −0.18 | 0.06 | 0.003 | −0.09 | 0.06 | 0.16 |
AJCC stage | ||||||
Stage Ia | ||||||
Stage II | −0.17 | 0.07 | 0.012 | −0.12 | 0.07 | 0.08 |
Stage III | −0.25 | 0.09 | 0.009 | −0.18 | 0.09 | 0.06 |
Mitotic rate | ||||||
<1 mitoses/mm2 tumora | ||||||
≥1 mitoses/mm2 tumor | −0.23 | 0.06 | 0.0004 | −0.20 | 0.06 | 0.001 |
Tumor site | ||||||
Heada | ||||||
Limbs | −0.04 | 0.10 | 0.63 | −0.11 | 0.10 | 0.28 |
Trunk | −0.35 | 0.10 | 0.001 | −0.36 | 0.10 | 0.001 |
Rare (sun-protected sites) | −0.44 | 0.12 | 0.001 | −0.38 | 0.13 | 0.003 |
. | Univariate . | Multivariateb . | ||||
---|---|---|---|---|---|---|
. | HR . | SE . | P . | HR . | SE . | P . |
Effect of tumor VDR expression on MSS (per unit expression) | 0.75 | 0.05 | 0.0001 | 0.80 | 0.06 | 0.008 |
. | Univariate . | Multivariateb . | ||||
---|---|---|---|---|---|---|
. | HR . | SE . | P . | HR . | SE . | P . |
Effect of tumor VDR expression on MSS (per unit expression) | 0.75 | 0.05 | 0.0001 | 0.80 | 0.06 | 0.008 |
aBaseline group used for comparison with relevant groups.
bMultivariate survival analysis was adjusted for AJCC stage, mitotic rate, tumor site, and tumor immune infiltrate. Regression coefficient and P value from linear regression.
Higher tumor VDR expression was protective for melanoma-related death independent of stage, tumor mitotic rate, tumor site, and histologically reported TILs (HR for melanoma-related death = 0.80; P = 0.008; Table 1). Although VDR expression correlated significantly with other members of the NR1lL family such as LXRB, FXR1, FXR2, and PXR expressions (Supplementary Table S1), it remained significantly prognostic after adjusting for the expression of those genes (adjusted HR for melanoma-related death = 0.77; P = 0.001).
Because VDR forms heterodimers with RXR and RXRγ signaling has been reported to drive epithelial–mesenchymal transition and invasion in melanoma, we tested the prognostic effect of the expression of genes coding for RXR receptors. In the LMC, VDR did not correlate significantly with RXRγ (Supplementary Table S1) and the prognostic significance of VDR expression was independent of the expression of RXRα, RXRβ, and RXRγ (adjusted HR for melanoma-related death = 0.75; P = 0.0001).
IHC of LMC primary melanomas revealed that VDR expression was predominantly in tumor (rather than immune) cells (representative images in Supplementary Fig. S1A–S1C, quantified in Supplementary Table S2).
In the TCGA metastatic melanomas (n = 353), higher VDR expression was protective for death (overall survival HR = 0.82; P = 0.03).
Tumor VDR expression is associated variably with CNAs and promoter methylation in primary and metastatic melanomas
In the LMC primary melanomas, VDR expression showed a weak positive correlation with VDR copy number, which failed to reach statistical significance (P = 0.12; Fig. 1A). However, in the TCGA metastatic melanomas, VDR expression was lowest in tumors with hemizygous deletion compared with those with “no changes/neutral” (P = 0.01) and was significantly higher in tumors with “gain” (P = 0.007) and “high-level amplifications” (P < 0.00001; Fig. 1B). VDR copy number was more frequently reduced in distant metastases compared with regional lymph node metastases (P = 0.015; Fig. 1C). Concomitantly, the TCGA data also showed progressively reduced VDR expression from primary to lymph node than distant metastases (Fig. 1D). VDR expression was significantly and inversely correlated with VDR promoter methylation (P = 0.0001; Fig. 1E) in TCGA metastases.
VDR correlates positively with genes enriched for immune-related pathways and negatively with proliferation-related pathways in LMC primary and TCGA metastases
In the LMC, an agnostic whole-transcriptome correlation analysis identified genes positively (n = 2,025) and negatively (n = 3,408) correlated with VDR (Fig. 2). The negatively correlated genes were enriched for proliferation-related pathways such as mitotic prometaphase, Wnt signaling, mitochondrial translation, tricarboxylic acid cycle, and cadherin signaling (Fig. 2; Supplementary Table S3). In contrast, the positively correlated genes were enriched for immune-related pathways such as cytokine–cytokine receptor interaction, TNF, IFNγ, IL12-mediated, NF-κB, and chemokine signaling (Fig. 2; Supplementary Table S4). VDR-correlated genes remained largely unchanged after FLG2 adjustment (Materials and Methods; Supplementary Table S5), indicating that confounding from epidermal sampling is unlikely.
We assessed whether the VDR-correlated genes in the LMC were known to have a VDR-binding site. Tuoresmaki and colleagues previously reported 54 nonoverlapping genomic VDR-binding regions recurrent in six tissue types, based on meta-analysis of VDR chromatin immunoprecipitation-sequencing (ChIP-seq) data (23). We mapped the 54 genomic-binding regions to be associated with 73 genes (GREAT, Supplementary Materials and Methods), of which, 43 genes (58%) were among the significant VDR-correlated genes in the LMC. Alternatively, 60% of genes mapped to genomic regions containing the VDR-binding motif (identified by Motifmap, see Materials and Methods) also correlated significantly with VDR in the LMC.
In the TCGA metastatic melanomas, VDR correlated negatively with genes enriched for: ECM organization, cadherin signaling, eukaryotic translation initiation, TGFβ, and VEGFR1 signaling; and positively with: NF-κB, TNF, IFNα/β, IFNγ, IL12-mediated, T-cell receptor, and chemokine signaling in naïve CD4 T-cells pathways (Supplementary Tables S6 and S7). The majority overlapped with those observed in LMC primaries.
Tumor VDR expression is associated with reported melanoma phenotypes, imputed immune cell scores, and reduced Wnt/β-catenin signaling
The LMC primaries were classified on the basis of previously described melanoma molecular phenotypes (Supplementary Materials and Methods; refs. 20, 36). The TCGA signature (20) classified the 703 LMC melanomas into immune (n = 192), keratin (n = 247), and MITF-low (n = 150) subtypes. VDR expression was significantly higher in the keratin and immune subtypes compared with the MITF-low subtype (P = 1.1 × 10−6; Fig. 3A). The Lund signature (36) classified tumors into high-immune (n = 173), normal-like (n = 198), pigmentation (n = 222), and proliferative (n = 83) subtypes. VDR expression was significantly higher in high-immune subtype compared with the poorer prognosis proliferative (P = 7.5 × 10−8) and pigmentation subtypes (P = 6 × 10−13; Fig. 3B).
VDR expression correlated positively with 25 of the 26 immune cell scores (Supplementary Table S8), of which, the strongest correlation (correlation coefficient > 0.30) was with dendritic cells, myeloid-derived suppressor cells, neutrophils, central memory CD4, NK, Th1, Th2, and T cells (Fig. 3C). Concordantly, VDR expression was significantly lower in tumors with “absent” immune infiltrate compared with tumors with “non-brisk” (P = 0.02) and “brisk” immune infiltrate (P = 0.004; Fig. 3D), according to histopathologic scores that were available for 601 (86%) of the LMC melanomas.
We previously reported an immunologically “cold” tumor subtype in the LMC (Consensus Immunome Cluster 4) with increased β-catenin signaling, reduced imputed immune scores for cytotoxic, T cell, and activated dendritic cells (aDC), and expression of genes coding for checkpoint molecules (16). VDR expression was lowest in that tumor subtype (Fig. 3E), which was concordant with the agnostic correlation analysis identifying Wnt/β-catenin as the pathway most strongly associated with low immune signals in both LMC primary and TCGA metastatic melanomas.
To compare VDR expression with response to immunotherapy, we used previously published therapy-response datasets (27, 28). VDR expression did not vary significantly between responders (n = 15) and nonresponders (n = 13) to anti-PD-1 therapy (P = 0.27; Supplementary Fig. S2B), nor to anti-CTLA4 treatment (P = 0.12) in a dataset of responders and (n = 13) nonresponders (n = 22; Supplementary Fig. S2A), although in the latter, VDR expression was 1.35-times higher in responders compared with nonresponders (FC = 1.35).
Deficient levels of serum vitamin D are associated with more melanoma-related deaths within the context of VDR expression
The 703 LMC primary melanomas were stratified into three groups using a survival-based stratification approach (X-tile, see Materials and Methods). VDR expression thresholds that best predicted differential survival identified the following groups: 17% with the lowest VDR expression (low-VDR group, n = 119), 17% of tumors with highest VDR expression (high-VDR group, n = 119), and middle 66% (intermediate-VDR group, n = 465) having the worst, best, and intermediate survival, respectively (P = 5.2 × 10−8; Fig. 4A and B). This was replicated in the TCGA metastatic melanomas (n = 353) using the same VDR expression percentiles (P = 0.03; Fig. 4C and D).
Among the three VDR groups: deficient serum vitamin D levels (<25 nmol/L) were associated with poorer prognosis compared with sufficient vitamin D levels (≥25 nmol/L) in the intermediate-VDR group (HR = 1.73; P = 0.02), but not in the low-VDR (P = 0.66) or high-VDR (P = 0.55) groups (Fig. 4E). The deleterious association with vitamin D deficiency was therefore apparent within the context of VDR expression. Intermediate-VDR group participants with deficient vitamin D were characterized by higher Breslow thickness (P = 0.02), higher frequency of pathologist-reported vascular invasion (P = 0.01), and AJCC stage II tumors (compared with stage I: P = 0.01), when compared with those with sufficient vitamin D. An agnostic whole-transcriptome analysis identified no gene significantly differentially expressed (at FDR < 0.10) between participants with deficient or sufficient vitamin D in the intermediate-VDR group. However, among the pathways that correlated significantly with VDR (Fig. 2), the following were significantly underexpressed in participants with vitamin D sufficiency: NK-cell–mediated cell killing (P = 0.02), IL12 (P = 0.03), and T-cell receptor signaling on naïve CD4 and CD8 signaling (P = 0.05).
Vitamin D treatment and increased VDR expression inhibit Wnt/β-catenin pathway and melanoma cell growth in vitro and in vivo
We functionally validated our findings from the LMC primary and TCGA metastatic melanomas, where high VDR-expressing tumors (with active vitamin D–VDR signaling) had better survival and reduced expression of proliferative pathways, in particular the Wnt/β-catenin signaling. The human melanoma cell lines, SK-MEL-28 and MeWo, treated with 1,25 (OH)2 vitamin D3 showed reduced proliferation posttreatment (Fig. 5A). At 24 and 48 hours posttreatment, VDR expression was significantly upregulated, while Wnt signaling and ECM organization genes were among the top downregulated pathways in both cell lines. The pathways upregulated in both cell lines at both time points were: MAPK, IFNα/β, TGFβ, and TLR signaling (Fig. 5B).
In a second functional validation model, the in vivo metastatic potential of murine melanoma B16BL6 cells overexpressing VDR (“VDR-B16BL6 cells”) was compared with control cells expressing low/no VDR (“control-B16BL6 cells”; Supplementary Fig. S3A and S3B). Mice injected with VDR-B16BL6 cells developed significantly lower pulmonary metastatic load (metastatic area and metastases count) compared with those dosed with control-B16BL6 cells (P < 0.04), with comparable results from both clones in two independent experiments (pooled analyses from two experiments represented in Fig. 5C). Differential expression of Wnt/β-catenin pathway genes between VDR-B16BL6 and control-B16BL6 was compared using a preformatted qRT-PCR–based array. Of the 84 Wnt/β-catenin genes tested, 62 genes had lower expression (FC < 1, of which, 26 had FC ≤ 0.5) in both VDR clones compared with control clones (Fig. 5D). Twelve genes had increased expression (FC > 1, none with FC ≥ 2) in both VDR clones.
In comparing the tumor immune infiltrate, although the number of CD3+ cells/100 mm2 met area was not statistically significantly higher in mice injected with VDR-B16BL6 cells (P = 0.11, compared with control-B16BL6 cells), there was a trend for increased CD3+ immune infiltrate (Supplementary Fig. S4).
Collectively, both models provide causal evidence that elevated vitamin D–VDR signaling in melanomas cells inhibits Wnt/β-catenin signaling and tumor proliferation.
Discussion
Melanoma is one of the most immunogenic malignancies, increased lymphocytic infiltration in both primaries (37, 38) and metastases (39) is associated with improved outcomes (16, 21) and melanoma responses to immune checkpoint blockade are high (40, 41). However, dampened immune responses and therapeutic resistance attributed to oncogenic pathways such as Wnt/β-catenin signaling (16, 42) mean that only 58% of patients with stage IV melanoma have significant benefit. Identification of factors that boost antitumor immunity is required to improve outcomes.
Low vitamin D levels are associated with thicker, poorer-prognosis primary melanomas (11), and lower VDR expression is associated with melanoma progression (43, 44), but neither causality nor the mechanistic basis has been established. In this article, we demonstrate that vitamin D–VDR signaling is protective for melanoma-related death at least, in part, through inhibition of Wnt/β-catenin signaling, impacting on melanoma proliferation and antitumor immune response.
We report that VDR expression was significantly lower in advanced tumors in the LMC primary melanomas. Importantly, high VDR expression was independently protective of melanoma-related death after adjusting for AJCC stage, mitotic rate, and TILs. Survival benefit was replicated in the TCGA metastatic tumors, highlighting the significance of vitamin D–VDR signaling in both primary and metastatic melanoma progression. This protective effect of VDR was independent of the expression of other NR1L family genes, despite reports of integrated activity between nuclear receptors (45). In assessing factors that could control VDR expression, VDR copy number was not significantly associated with expression in LMC primaries. However, distant metastases (which have worse prognosis) had lower VDR copy number compared with regional lymph node metastases in the TCGA, suggesting a progression-associated genomic loss of VDR. Low VDR–expressing metastatic tumors from TCGA were also more likely to be hypermethylated, consistent with previous reports of the epigenetic control of VDR expression (46). Despite the progressive deletion of the VDR locus with tumor progression in metastatic disease, the data suggest that therapeutic manipulation of vitamin D–VDR signaling could have adjuvant therapeutic benefit in primary disease where we saw little evidence of VDR deletion.
An agnostic correlation analysis revealed that VDR expression was strongly positively correlated with immune-related pathways and negatively correlated with proliferation-associated pathways in LMC primaries and TCGA metastases. An additional sensitivity analysis adjusting for FLG2 expression produced no significant changes in the correlated pathways, indicating that artefact from keratinocyte-rich tumor populations was unlikely. VDR expression was higher in “brisk” immune-infiltrated primaries compared with tumors with no immune infiltrate, which is in agreement with a previous report of a smaller IHC study (43) and is an independent validation of the transcriptomic imputation of immune infiltration. Fifty-eight percent of the genes within reported genomic VDR-binding regions (identified from VDR ChIP-seq data as well as VDR-binding motif) were found to correlate with tumor VDR expression in the LMC. This is consistent with direct transcriptional control by the VDR transcription factor for a proportion of differentially expressed genes.
Validation using published prognostic melanoma molecular phenotypes (16, 36) revealed that VDR expression was significantly higher in high-immune subtypes compared with proliferative subtypes, consistent with the view that the prognostic significance of VDR is associated with increased immune and decreased proliferative signaling. We also assessed whether VDR preferentially correlated with a particular immune cell type, which was previously uncharacterized in primary melanomas. However, VDR was strongly positively correlated with all imputed immune cell scores, and we have previously reported simultaneous upregulation of adaptive and innate immunity in good prognosis primary melanomas (16). Furthermore, the proimmune effect of vitamin D–VDR signaling was supported by strong positive correlation of VDR with genes involved in pathways such as extracellular matrix organization, TNFα, NF-κB, IFNγ, and IL12-mediated signaling.
In comparing pretreatment gene expression between responders and nonresponders to immunotherapy, VDR expression was higher in responders to anti-CTLA4 (FC = 1.35), albeit not statistically significant (P = 0.12) in this very small dataset. Although data from these immunotherapy studies did not support VDR expression as a biomarker of response, we posit that the datasets available to explore VDR as a biomarker were insufficient to properly explore this possibility.
Wnt/β-catenin signaling was among the top negatively correlated pathways with VDR in the LMC and TCGA melanomas, concordant with reports of vitamin D-VDR–mediated inhibition of Wnt/β-catenin signaling in colon cancer (47). We further explored this relationship and report that VDR expression was lowest in previously reported subset of tumors characterized by high Wnt/β-catenin expression, reduced immune infiltrate, and high mortality (16). Collectively, our findings support the hypothesis that, as in colon cancer, some of the effects of vitamin D–VDR signaling in melanoma are mediated by inhibition of Wnt/β-catenin signaling.
Vitamin D deficiency (≤25 nmol/L, ∼10 ng/mL) was associated with worse prognosis only in participants with intermediate-VDR expression (albeit the majority). Although an agnostic analysis identified no significant transcriptomic differences associated with vitamin D in this subgroup, there was some evidence for paradoxically reduced expression of immune-associated pathways with higher vitamin D levels in a candidate gene expression analysis. The reported associations between vitamin D and the immune system are numerous and complex. The findings of this study are therefore not inconsistent with the view that vitamin D deficiency should be avoided but that high levels would not necessarily be beneficial to all patients.
The lack of a protective effect of higher vitamin D levels in participants with low-VDR tumors was not unexpected as low receptor expression could preclude effective signaling despite ligand sufficiency. A lack of benefit in the high-VDR tumors was more surprising and we postulate that this could reflect receptor saturation as reported in other NHR family receptors (48) or a ligand-independent effect of VDR, which has been described in other cancers (49, 50).
In functional validation, the observed vitamin D-treatment–induced reduction in cell proliferation is concordant with previous reports (51). However, the pan-transcriptome–based findings that vitamin D treatment of human melanoma cell lines promotes VDR expression and inhibition of protumor pathways including the Wnt/β-catenin pathway, is novel. Furthermore, elevated VDR expression in murine melanoma cells decreased their in vivo metastatic potential, the expression of key Wnt/β-catenin genes, some of which (Dkk1 and Sfrp2) have previously been shown to be inhibited by VDR (52), but not in melanomas. Interestingly, the “classic” noncanonical Wnt ligands Wnt5a, Wnt5b, Wnt10a, Wnt7, and Wnt11 were also downregulated by VDR. This finding is of significance because Wnt5a (and some other noncanonical Wnt ligands) affect cell motility and invasion and is implicated in worse melanoma prognosis (53). Thus, the findings functionally validate the transcriptome-derived evidence for the inverse association between VDR and Wnt/β-catenin signaling.
We report a trend toward increased T-cell infiltration of lung metastases produced by VDR-expressing B16BL6 melanoma cells, although this did not reach statistical significance. Rejection of B16-derivative melanoma cell lines in vivo requires both NK and T cells, and in the pulmonary metastasis model, NK-cell–derived IFNγ is important but T cells also exert important antitumor activity (33–35, 54–57). The lack of statistically significant reduction in T cell numbers does not exclude the possibility that T cells are more active and/or that NK-cell recruitment and activity is also modified. To this effect, further analyses of T-cell and NK-cell recruitment and activity following VDR manipulation is warranted using both the subcutaneous and metastatic B16 models.
This study reports that vitamin D–VDR signaling bestows a prognostic benefit for patients with melanoma by inhibiting Wnt/β-catenin signaling and increasing immune cell infiltration. These findings also suggest that activating vitamin D–VDR signaling has the potential to enhance antitumor immunity in an adjuvant setting. Notably, our findings suggest that vitamin D deficiency (<25 nmol/L) is deleterious for melanoma survival rather than that high levels are protective. As melanoma is causally related to intense sun burn (58), sun avoidance is frequently recommended to patients in follow-up. Our data suggest a causal relationship between reduced vitamin D–VDR signaling and therefore, as sun exposure is the dominant vitamin D source in most populations, simultaneous avoidance of vitamin D deficiency is important health advice.
Disclosure of Potential Conflicts of Interest
S.J. O'Shea has received speakers bureau honoraria from L'Oreal Ltd UK, is an unpaid consultant/advisory board member for British Association of Dermatologists Research Subcommittee, and has provided expert testimony for Janssen Ltd (sponsorship for travel and subsistence to attend the EADV meeting in 2018, the BAD meeting in 2019, and EADV meeting in 2018) and Pierre Fabre Ltd (sponsorship for travel and subsistence to attend the Avene hydrotherapy centre in 2018). D.J. Adams is a consultant (paid consultant) for Microbiotica and reports receiving a commercial research grant from BMS. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: S. Muralidhar, J. Reichrath, D.T. Bishop, J. Newton-Bishop
Development of methodology: S. Muralidhar, D.T. Bishop, J. Newton-Bishop
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Filia, J. Nsengimana, S.J. O'Shea, M. Harland, L. van der Weyden, D.J. Adams, D.T. Bishop, J. Newton-Bishop
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Muralidhar, A. Filia, J. Nsengimana, J. Poźniak, J.M. Diaz, M. Harland, L. van der Weyden, D.J. Adams, D.T. Bishop, J. Newton-Bishop
Writing, review, and/or revision of the manuscript: S. Muralidhar, A. Filia, J. Nsengimana, J. Poźniak, S.J. O'Shea, J.M. Diaz, M. Harland, J.A. Randerson-Moor, J. Reichrath, J.P. Laye, L. van der Weyden, D.J. Adams, D.T. Bishop, J. Newton-Bishop
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.A. Randerson-Moor, D.T. Bishop, J. Newton-Bishop
Study supervision: J.P. Laye, D.T. Bishop, J. Newton-Bishop
Acknowledgments
This work was funded by Cancer Research UK C588/A19167, C8216/A6129, and C588/A10721 and NIH CA83115. J. Poźniak, J.M. Diaz, and S. Muralidhar were funded by Horizon 2020 Research and Innovation Programme no. 641458 (MELGEN). Our grateful thanks to the participants who gave of their time and their blood samples. Also to the research nurses and technicians who collected the data over many years: Susan Leake, Susan Haynes, Birute Karpavicius, Paul Affleck, Kairen Kukalizch, Linda Whitaker, Sharon Jackson, Edwina Gerry, Elaine Fitzgibbon, Clarissa Nolan, Saila Waseem, Yvonne Taylor, Pauline Brunyee, Paul King, Tracy Lee, Samira Lobo, and Minttu Polso. To Jo Gascoyne and May Chan who provided critical support in managing the complex dataset. We thank Dr. Robert Salmond, University of Leeds (Leeds, UK), for his guidance and access to cell culture facilities and Professor Graham Cook for immunologic insights. We acknowledge the advice and assistance of Kenneth MacLennan, Consultant Professor of Tumour Pathology at Leeds Teaching Hospitals NHS Trust, in reviewing IHC sections. We thank Dr. James Hewinson (Experimental Cancer Genetics group, Wellcome Sanger Institute, Cambridge, UK) for his valuable assistance to the in vitro transfection experiments. We thank Dr. Martin del Castillo for aiding in the choice of appropriate cell line for the in vitro and in vivo experiments.
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.