1 INTRODUCTION

Cutaneous melanoma is still increasing in incidence and mortality in many countries. Although immunotherapy has improved outcomes for patients with advanced disease, 40% do not benefit, and the need to better understand the biology of the tumour and host variation in response remains crucial. Microscopic ulceration of primary melanoma is a strong independent predictor of melanoma death, but its biology remains unclear. We have reported evidence previously that ulceration may at least in part, be driven by host systemic inflammation in that obesity, diabetes, vitamin D deficiency and smoking were associated with ulceration, in the Leeds Melanoma Cohort (LMC) (Newton-Bishop et al., 2015). This hypothesis was subsequently supported by an Australian study in which ulceration in thick tumours was associated with diabetes, and that statin use (which is reported to reduce IL-6 levels (Sepehri et al., 2016)) was protective for ulceration (von Schuckmann et al., 2017). IL-6 and other molecules resulting from IL-1β signalling mediate systemic inflammation associated with an increased risk of cardiovascular disease, and suppression of this pathway has been reported to reduce lung cancer deaths (Ridker et al., 2011). Drugging systemic inflammation therefore remains of interest as a potential adjuvant therapy for melanoma.

Genomic studies of primary melanomas are relatively few as a result of the tiny volume of many tumours (Mar et al., 2014) with little evidence of a genomic basis for ulceration (Arbiser, 2014; Mascaro et al., 1984). In a recent study by Koelblinger et al. (2019), ulcerated melanomas had a significantly higher proportion of tumour cells expressing the immunomodulatory protein PD-L1 (programmed death ligand 1), suggesting that the biological processes driving ulceration could impact immune responses. The study also reported increased expression of CD11c (a marker of dendritic cells) and CD68+ and CD163+ macrophages (Koelblinger et al., 2019) in ulcerated tumours. The hypothesis of our study was that tumoural genomic variation plays a significant role in driving the ulceration phenotype in primary melanomas, modified by potentially modifiable environmental exposures. A systems biology approach to its exploration was taken as a model of approaches to understanding complex interactions between host and tumoural variation.

2 METHODS 2.1 Human data and samples

The LMC is a prospective cohort of 2184 primary melanoma patients recruited from a geographically defined area of the UK in the period 2000–2012, see Supplementary Information. The Vitamin D and Immunity Study (VDI) REC reference 13/YH/0237 is a second cohort of 393 primary melanoma patients recruited at diagnosis, making their first visit to the Leeds Melanoma Multidisciplinary Clinic. Here, we report the analysis of ulceration status in association with clinicopathological variables and serum vitamin D levels at diagnosis. The VDI histopathological data were extracted from the summary of the Leeds multidisciplinary team review, which was carried out for all participants according to protocol.

2.2 Statistical and bioinformatic analyses

See Supplementary Methods.

2.3 Association of tumour and host variables with ulceration and melanoma-specific survival

The analysis of clinicopathological factors used data collected from the entire LMC cohort of 2184 people and that from 393 participants in the VDI study. Tumour analyses were carried out on the subset of LMC 703 tumours which were large enough to sample yet leave sufficient tumour in the formalin-fixed tumour blocks for subsequent clinical testing if needed. The quantity of DNA and RNA extracted from tumour tissue was limited by tumour size, and therefore, transcriptomic and sequencing data did not overlap completely (Figure S1).

2.4 Measurement of circulating inflammatory markers and serum vitamin D from blood samples taken at recruitment to the study

Sera stored at −80℃ since recruitment (duration of storage 6–18 years) were used to measure C-reactive protein (CRP) and albumin in LMC participants whose tumour genomic data were available (performed at the University of Glasgow) (Salim et al., 2016). This approach was chosen as these measures are established to be stable in cryopreserved serum over time. We computed the Glasgow Prognostic Score (GPS) (Ohmura et al., 2017) derived from CRP/albumin levels, as this is reported to better reflect systemic inflammation and to predict poorer outcomes from colorectal cancer (Nozoe et al., 2014). The individual values and the score were tested as associated with ulceration status and melanoma-specific survival (MSS). The vast majority (93%) were scored 0 (i.e. CRP ≤10 mg/L and albumin ≥35 g/L) and a minority scored 1 (i.e. CRP >10 mg/L or albumin <35 g/L). Patients who had both a serum elevation of CRP and hypoalbuminaemia were allocated a GPS of 2. Vitamin D levels were measured as reported previously (Newton-Bishop et al., 2009).

In the VDI study, serum vitamin D levels were measured at diagnosis as reported for the LMC (Newton-Bishop et al., 2009) and high sensitivity CRP was measured by the Leeds NHS laboratory in fresh samples.

2.5 Tumour genomics and ulceration in the LMC 2.5.1 Genomic data generation

Formalin-fixed paraffin-embedded (FFPE) primary melanomas were sampled using a 0.6 mm diameter tissue micro array (TMA) needle inserted horizontally consistently through the least inflamed part of the invasive tumour. DNA/RNA were extracted from 820 tumour cores (703 unique patients and 117 duplicates) as described previously using Qiagen AllPrep DNA/RNA FFPE kits (Jewell et al., 2015; Nsengimana et al., 2018). Gene expression was quantified using the Illumina DASL Human HT12 v4 array). DNA samples were used to generate copy number data from 303 of the samples as reported (Filia et al., 2019). Tumour DNA was also processed at the Wellcome Sanger Institute to generate somatic variant calls (mutations) from 524 of the primaries compared with matched blood DNA (465 cases) or a high depth control DNA sample filtered using data from ExAc (http://exac.broadinstitute.org/). Samples were selected only on the basis of there being sufficient DNA available. Targeted capture was performed on 554 genes using Agilent SureSelectXT probes as described by Chen et al. (Birkealv et al., under review, Nature Communications, December 2021). The baits included 164 melanoma-associated genes, 245 known to be associated with other solid cancers, common melanoma promoter mutations, and 101 genes from the interferon signalling pathway (Gao et al., 2016) (Supplementary Excel File S1). Human leucoyte antigen (HLA) regions were also screened. Sequencing designed to generate copy number data was carried out using the Illumina HiSeq4000 platform, using 75 bp paired-end reads and data were mapped to GRCh37d5 with BWA-MEM 0.7.15 and somatic variants were called using the Caveman algorithm (v.1.11.2). The mutation load was analysed for its association with ulceration, defined as summated non-synonymous mutations. As the melanoma primary samples are small, there was insufficient material to process all samples using the three different platforms and the overlaps are illustrated in Figure S1.

2.5.2 Analysis of genomic data

Associations between individual mutations detected and ulceration status were tested using logistic regression (univariate and adjusted for thickness, age and sex). A candidate gene approach was also taken to the analysis of expression of genes coding for cytokines and their association with ulceration status. Differential gene expression by ulceration status was used to identify biological pathways associated with ulceration in MetaCoreTM (Table 2). The transcriptomic data were also subjected to a bioinformatic inference of the presence of specific immune cell subgroups in tissues (Angelova et al., 2015), in association with ulceration status. Gene level data associated with ulceration status were generated from segmented copy number using GISTIC (Software.broadinstitute.org). Finally then, the copy number and transcriptomic data were combined in a neural network analysis (Abdel-Fatah et al., 2016). The resultant analysis identified ‘influencer’ genes and ‘influenced’ genes associated with ulceration, which were then subjected to pathway and network analysis using MetaCore to explore significantly associated biological systems. We did not analyse mutation data in this combined analysis as the samples in common between all three genomic data sets were significantly smaller (Figure S1).

2.5.3 Utilisation of artificial neural network-based network inference to determine molecular influence and perturbation of molecular interaction networks in merged transcriptomic and copy number data

Here, we implemented a machine learning approach based on an artificial neural network (ANN) combined with a concordance analysis conducted across multiple Monte Carlo data splits (where four levels of cross validation are conducted). We previously showed this to be highly effective at eliminating false discovery and overfitting, while maximising generalisation of the biomarkers identified (Abdel-Fatah et al., 2016) (Figure S2). The concordant set of markers determined was fed into an ANN-based network inference algorithm (ANNi) (Tong et al., 2014). This approach uses each transcript or gene identified with a copy number change, from the enriched set as a network output and all remaining features as network inputs. This approach was repeated for each data point, creating a group of interaction models. The trained, optimised interaction models were then collectively analysed to define an interaction map for the genomic data points (Figure S4), using the summed weights leading from a given input to a given output. The visualisation of interactome network maps of concordant genes was undertaken in Cytoscape (version 2.8). The Pearson correlation coefficient r with a cut-off value of 0.7 was implemented in the algorithm for removing the least significant interaction scores.

3 RESULTS 3.1 Clinicopathological data analyses

There were 187 ulcerated tumours in the LMC (of 576 total examined, 32%) and 153 ulcerated in the VDI study (out of 391 total 39%). Ulceration was more common with increasing age in both the LMC and VDI data sets, and this was independent of other variables in the larger LMC data set (Table 1). Ulceration was more common in men but this was not independent of other factors. Vitamin D deficiency and smoking were independently associated with ulceration in the LMC as previously reported (Newton-Bishop et al., 2015). Increased body mass index (BMI) was associated with ulceration odds ratio (OR) 1.03 per BMI unit (1.00–1.05), p = .02 but, after taking into account vitamin D status, the effect size was less significant (OR 1.02 per unit [0.99–1.04] p = .2).

Although the VDI study was very much smaller (having only 18 type II diabetics and 26 current smokers) than the LMC therefore reducing statistical power, the VDI data were consistent with the previous observation of a role for vitamin D. The OR for ulceration in vitamin D deficiency was = 1.18 (0.42, 3.26), p = .8 and for MSS: hazard ratio (HR) = 1.48 (0.46, 4.81), p = .5, adjusted for age and sex. There was no significant association with sex or smoking and there was a borderline association between BMI and ulceration status in this data set (OR 1.03 [1.00, 1.04], p = .07) after adjusting for age and sex.

TABLE 1. Association of clinicopathological factors with ulceration in LMC and VDI data sets Leeds Melanoma Cohort (LMC) Vitamin D and immunity study (VDI) Adjusted for age and sex Multivariable (n = 1007) Adjusted for age and sex Multivariable (n = 368) N OR (95% CI) p OR (95% CI) p N OR (95% CI) p OR (95% CI) p Age (pa) 2147 1.03 (1.02, 1.04) <.001 1.02 (1.00, 1.04) .01 384 1.02 (1.00, 1.03) .03 1.01 (0.99, 1.02) .4 Sex 2147 384 Female 1223 1 - 1 - 187 1 - 1 - Male 924 1.24 (1.00, 1.5) .05 0.89 (0.60, 1.31) .6 197 0.99 (0.65, 1.49) 1 0.88 (0.53, 1.44) .6 TILs 1039 370 Absent 251 1 - 1 - 52 1 - 1 - Non-brisk 604 0.82 (0.57, 1.17) .3 1.07 (0.70, 1.65) .8 261 1.17 (0.63, 2.17) .6 1.58 (0.78, 3.18) .2 Brisk 184 0.55 (0.34, 0.92) .02 1.18 (0.65, 2.14) .6 57 0.68 (0.30, 1.52) .3 1.22 (0.49, 3.07) .7 Breslow, mm 2128 382 <1.0 605 0.28 (0.17, 0.44) <.001 0.29 (0.11, 0.75) .01 44 0.22 (0.06, 0.79) .02 0.31 (0.09, 1.11) .07 1–0.99 814 1 - 1 - 101 1 - 1 - 2–3.99 474 3.41 (2.57, 4.51) <.001 3.83 (2.51, 5.86) <.001 ** 131 2.04 (1.14, 3.63) .02 1.75 (0.96, 3.21) .07 >4 235 8.65 (6.19, 12.1) <.001 6.17 (3.64, 10.5) <.001 *** 106 5.65 (3.05, 10.5) <.001 4.51 (2.27, 8.94) <.001 Site 2147 376 Trunk 765 1 - 1 - 138 1 - 1 - Limbs 210 0.74 (0.56, 0.97) .03 0.73 (0.47, 1.13) .2 56 1.31 (0.79, 2.15) .3 1.38 (0.79, 2.41) .3 Head/neck 959 1.00 (0.68, 1.46) 1 0.76 (0.41, 1.43) .4 155 1.16 (0.60, 2.24) .7 0.92 (0.44, 1.91) .8 Sun-protected sites 213 3.98 (2.84, 5.57) <.001 3.96 (2.27, 6.90) <.001 * 27 2.36 (1.01, 5.54) .05 2.11 (0.79, 5.62) .1 Mitoses 1844 376 <5 mm−2 1331 1 - 1 - 198 1 - 1 - >5 mm 513 4.9 (3.87, 6.28) <.001 2.69 (1.85, 3.89) <.001 ** 178 3.52 (2.26, 5.48) <.001 2.17 (1.31, 3.59) .003 Univariable analysis Log coding mutation count per unit 498 0.84 (0.72, 0.99) .04 Uni-variable analysis**** Note Data are presented as univariable and multivariable analyses for each data set, for which logistic regression models were fitted with ulceration as the outcome variable. Missing data, for example where mitotic rates had not been reported by the clinical pathologists, explain the differing numbers in the univariable analyses. The total mutation count was available only for 498 individuals in the LMC data set and the statistical power was considerably reduced if this variable was included in the multivariable analysis (274 individuals); therefore, we report only a univariable analysis for this variable. The stars (as below) show the level of significance, however, for variables adjusted for mutation count. Associations significant at the 5% level or higher are highlighted in bold. CI, confidence interval; OR, odds ratio. 3.1.1 Clinicopathological factors and mutation load associated with ulceration

Ulceration was independently more common in thicker tumours and was associated with a higher mitotic rate (in both LMC and VDI studies) (Table 1). Ulceration was more common in tumours arising in sun-protected sites such as acral and genital tumours in the LMC data set independent of thickness and age: OR for ulceration was 3.96 95% confidence interval (CI) 2.27, 6.90, p < .001 in the LMC multivariable analysis with support for this in the smaller VDI data set. Brisk tumour-infiltrating lymphocytes (TILs) (as reported in clinical reports) were less common in ulcerated tumours (OR 0.55 95% CI 0.34, 0.92, p < .02) in the LMC in a univariable analysis, but this was not significant in the multivariable analysis nor in the VDI study. Total tumour mutation load data were available for 498 participants in whom blood test results were also available, and in a univariable analysis, there was evidence that ulceration was less frequent in tumours with a higher mutation load. This did not, however, persist in the multivariable analysis, which was then a relatively small data set and therefore inadequately powered.

Vascular invasion was seen more frequently in the ulcerated tumours (p = .0003) in the LMC data set, and this was independent of tumour thickness p < .001. In the smaller VDI data set the difference did not reach statistical significance (p = .3) (Table S1). Tumour regression was not consistently commented on in the clinical histopathology reports but was assessed in a single observer analysis of the LMC tumours sampled for genomic studies. Ulcerated tumours were more likely to have regression when the whole slide was considered, compared to non-ulcerated tumours (chi-squared test, p = .001) and the depth of regression was also more likely to be greater within ulcerated tumours (whole tumour measure: none, <0.71 mm, 0.71 to <1.14 mm, 1.14 to <1.63 mm, ≥1.63 mm, Fisher's exact test, p = .00001).

3.1.2 Examination of the cored tumour site by single observer

There was no difference in the presence or absence of an immune infiltrate detected histologically in the region of the tumour sampled for genomic studies (Fisher's exact test, p = .4) comparing ulcerated with non-ulcerated tumours. Ulcerated tumours were, however, more likely to have a core immune infiltrate classified as ‘none or barely perceptible’ compared to non-ulcerated tumours (Fisher's exact test, p = .001). In terms of specific subtypes of immune cells, ulcerated tumours were less likely to have lymphocytes (Chi-squared test, p = .0005) or macrophages (chi-squared test, p = .0003) but more likely to contain neutrophils within the core, although the absolute numbers were very small (Fisher's exact test, p = 3 × 10−6). There were only 12 cases with neutrophils in the core and 7 of these were ulcerated, representing 4% ulcerated cases, compared to 0.9% non-ulcerated cases.

3.1.3 Clinicopathological factors and MSS in the LMC

Table S2 shows the analysis of clinicopathological factors predictive of death from melanoma in the mature LMC and as was expected, age, tumour thickness, tumour site, mitotic rate and TILs were independently predictive. Ulceration was significantly predictive only in the univariable analysis. The data are presented in the table as univariable and multivariable analyses except for mutation load as this factor was measured only in a subset of tumours. In a univariable analysis, a higher mutation load was significantly associated with fewer melanoma deaths. We will report detailed analysis of associations with MSS in the VDI data set when the survival data are more mature.

3.1.4 Ulceration and death from causes other than melanoma in the LMC

We looked at ulceration and death from causes of death other than melanoma as we had hypothesised that ulceration may also predict deaths related to diseases mediated by systemic inflammation (Table S3). There were 153 non-melanoma deaths in 2018 participants in this data set, and the majority occurred in individuals with diseases associated with systemic inflammation. Of these, 2 had autoimmune disease, 56 cardiovascular events (including cerebrovascular accidents), 10 had chronic lung disease, 11 chronic neurological disease including Alzheimer's, 56 died of other cancers and there were 14 additional deaths caused by rare events such as trauma or the cause of death was unclear. Death from causes other than melanoma in participants with ulcerated tumours was higher in the LMC overall HR 1.47 95% CI (0.98–2.21) p = .06 although this was not statistically significant when the analysis adjusted for age and sex HR = 1.14, 95% CI (0.76–1.72), p = .5. This was a little stronger in males HR 2.03 95% CI (1.25–3.31), p = .004, although this also lost significance when adjusted for age and sex HR = 1.57, 95% CI (0.96–2.57), p = .1.

3.2 Circulating inflammatory markers

When analysed on a continuous scale, blood CRP and albumin in 695 samples from the LMC showed no significant association with ulceration, TILs, AJCC stage (Figure 1a) or MSS (Figure 1a–d). However, when dichotomised on the median, 11% of participants whose albumin level was lower than the median, had brisk TILs compared to 21% of those with higher levels (p = .02, Figure 1b). Ninety-three percent were scored 0 using the Glasgow Predictive Score (GPS) (i.e. CRP <10 mg/L and albumin >35 g/L), 7% scored 1 (i.e. CRP >10 mg/L or Albumin <35 g/L), and none scored 2. We tested the log2(ratio) of CRP/albumin against MSS, generating a survival curve for the top quartile compared with the others (Figure 1e). All 45 patients with GPS = 1 were in the top quartile of the ratio, with an additional 128 with GPS = 0. Neither the GPS nor the ratio predicted ulceration status: values for the GPS score (p = .64), AJCC stage (p = .58) or TILs (p = .16), but GPS predicted MSS (HR = 1.6, p = .05) as did the log2ratio (Figure 1e). Taken individually, CRP and albumin were negatively correlated with each other (Spearman rho = −.27) and had significant opposite and independent associations with MSS: HR = 1.1 (p = .02) per CRP unit (log2 scale) and HR = 0.39 (p = .02) per albumin unit (log2 scale). These effects were independent of the AJCC stage, although they fell short of significance when dichotomised on the median (Figure 1).

Associations between inflammatory markers CRP and albumin in the LMC participants whose tumours were sampled. Analysis of circulating inflammatory markers in 695 participants in the LMC whose tumours were analysed. (a) The p values for tests of association between CRP, albumin levels and the GPS and key variables, tumour ulceration, TILs and AJCC stage. GPS = 0 if CRP ≤10 mg/L and albumin ≥35 mg/L. GPS = 1 if either CRP was higher than range or albumin lower and GPS = 2 if both CRP and albumin were higher than threshold. (b) Association between albumin and CRP levels and the presence of brisk TILs in the primary tumours: here, the analysis was stratified at the median. (c) Kaplan–Meier curve for MSS by CRP level. (d) Kaplan–Meier curve for MSS by albumin and (e) Kaplan–Meier curve for MSS by log ratio2 ratio CRP/albumin. CRP, C-reactive protein; GPS, Glasgow Prognostic Score; LMC, Leeds Melanoma Cohort; MSS, melanoma-specific survival; TIL, tumour-infiltrating lymphocyte

In the VDI study, 372 samples had measurable CRP and 386 samples measurable fibrinogen. No significant associations were seen with ulceration for either CRP (adjusted for age, sex OR = 1.00 95% CI 0.37, 2.66, p = 1) or fibrinogen (adjusted for age, sex OR = 1.58 95% CI 0.85, 2.91, p = .1). Nor were associations seen with presence of brisk TILs for either CRP (adjusted for age, sex OR = 0.35 95% CI 0.05, 2.74, p = .3) or fibrinogen (adjusted for age, sex OR = 1.37 95% CI 0.59, 3.19, p = .5). No significant associations were seen for either CRP or fibrinogen with TILs (categorised as absent, non-brisk, brisk).

3.3 Tumour genomics and ulceration in the LMC 3.3.1 Mutations and ulceration status

In the LMC data, the commonest mutations detected were TERT promoter mutations along with coding mutations in BRAF, CDKN2A and NRAS. These, and all the other coding mutations identified were tested individually for association with ulceration. A number of them showed a marginal association with ulceration in a univariate analysis and after adjusting for tumour thickness, mitotic rate and age of patient at diagnosis, although none remained significant when the data were corrected for multiple testing. Figure 2a shows the mutations significantly associated with ulceration (unadjusted p < .05 from at least 10 tumours with one or more mutations). TP53 and APC were the most significantly associated with ulceration with opposing associations: mutations in TP53 increased the odds of ulceration 2.5-fold (unadjusted p = .004) while APC mutations decreased the odds of ulceration 10-fold (unadjusted p = .001, see Figure 2a,b). These associations remained significant after adjusting for thickness, age and sex. APC gene expression was borderline significantly lower in the tumours with APC mutations (p = .06) as was also true for TP53 (p = .03). Of the most common driver mutations, those in BRAF and NRAS mutations were neutral in terms of MSS while NF1 was protective for death (adjusted OR = 0.36, p = .02, see Figure 2a). The total coding mutation load in the 554 genes sequenced was not associated with ulceration status (Figure 2c) but a higher load predicted an improved MSS (Figure 2d). There was no association of mutation load with TILs (p = .21) or AJCC stage (p = .51).

Mutations associated with ulceration status. Association between non-synonymous mutations detected in 455 tumours subjected to targeted sequencing, and ulceration or MSS. (a) The most significant associations between individual gene mutations (presence/absence) and ulceration status, adjusted and unadjusted for covariates (not corrected for multiple testing). (b) Heterozygous versus homozygous mutations in the two genes most strongly associated with tumour ulceration in 166 ulcerated and 289 non-ulcerated tumours. (c) Lack of association between ulceration and the load of non-synonymous mutations. (d) Association between the load of non-synonymous mutations and MSS. Comparisons were made between samples with mutation loads recorded at above or below the median. MSS, melanoma-specific survival

3.3.2 Transcriptomic data

We then sought to identify biological pathways associated with ulceration in the 703 LMC primary melanomas for which transcriptomic data were available. A whole-genome univariate test (Mann–Whitney) revealed that 4660 genes varied in expression significantly with ulceration (Benjamini–Hochberg false discovery rate 5%). Genes whose expression was significantly higher in ulcerated tumours (n = 2681) were enriched for pathways related to cell proliferation, such as mitotic prometaphase, signalling by Rho-GTPases, cell cycle checkpoints and mitochondrial translation. Conversely, the genes which had significantly lower expression in ulcerated tumours (n = 1979) were enriched for pathways related to the tumour microenvironment including immune signalling such as extracellular matrix (ECM) organisation, ECM–receptor interaction, cytokine–cytokine receptor interaction, interferon gamma signalling and cell adhesion pathways such as β1 integrin cell surface interactions and cell adhesion molecules (Table 2). The inferred immune cell subgroups (Table 3) showed evidence for increased signals of eosinophils and activated CD4 cells in ulcerated tumours, but the data suggested that overall signals from immune cell subgroups were largely reduced. The subgroups reduced in ulcerated tumours included natural killer and T cells including Th1 and Th2 cells.

TABLE 2. Pathway analysis of the tumour genes differentially expressed in ulcerated compared with non-ulcerated tumours Pathway associated Direction of association of genes up- or downregulated in ulcerated tumours Significance of the association (p-value) Mitotic prometaphase(R)