Predicting Clinical Outcome through Molecular Profiling in Stage III Melanoma

Purpose: Patients with macroscopic stage III melanoma represent a heterogeneous cohort with average 5-year overall survival rates of <30%. With current algorithms, it is not possible to predict which patients will achieve longer-term survival. We hypothesized that molecular profiling could be used to identify prognostic groups within patients with stage III melanoma while also providing a greater understanding of the biological programs underpinning these differences.

Experimental Design: Lymph node sections from 29 patients with stage IIIB and IIIC melanoma, with divergent clinical outcome including 16 “poor-prognosis” and 13 “good-prognosis” patients as defined by time to tumor progression, were subjected to molecular profiling using oligonucleotide arrays as an initial training set. Twenty-one differentially expressed genes were validated using quantitative PCR and the 15 genes with strongest cross-platform correlation were used to develop two predictive scores, which were applied to two independent validation sets of 10 and 14 stage III tumor samples.

Results: Supervised analysis using differentially expressed genes was able to differentiate the prognostic groups in the training set. The developed predictive scores correlated directly with clinical outcome. When the predictive scores were applied to the two independent validation sets, clinical outcome was accurately predicted in 90% and 85% of patients, respectively.

Conclusion: We describe a gene expression profile that is capable of distinguishing clinical outcomes in a previously homogeneous group of stage III melanoma patients.

In industrial nations, the incidence of melanoma has steadily risen over the past 25 years, with the incidence in Australia being the highest in the world (1). Although the perceived “melanoma epidemic” largely represents increased detection of thin melanomas (2), melanoma affects predominantly younger age groups and results in a loss of productive life years exceeded only by childhood malignancies and testicular cancer (3, 4). Melanoma is largely unresponsive to cytotoxic chemotherapy (5), biological agents (6, 7), and various vaccination strategies (8). A small subgroup of patients seems to benefit from biological and/or cytotoxic chemotherapies, but identifying these patients a priori is currently not possible, which necessitates the exposure of many patients to substantial toxicities with a low probability of benefit.

Molecular classification of various tumor types, including breast, bowel, and lymphoma, has consistently shown subtypes of tumor with prognostic outcomes independent of those expected with standard staging procedures (9–11). In the case of breast cancer, molecular profiles have been used to subclassify cancers previously thought to be homogenous (12), allowing prediction of those most likely to benefit from neoadjuvant chemotherapy (13, 14), and to predict overall survival (15, 16). Although the application of microarray technology to melanoma is intuitively useful (17–19), few studies incorporating clinically applicable techniques have been published to date. Most research has focused on cell lines rather than ex vivo patient tumor tissue and the clinical application of these results is therefore limited (20–22).

The heterogeneity of clinical outcome in patients with locally advanced melanoma (23) suggests that factors intrinsic to the tumor or to the host's antitumor immune response may be important (24, 25). To investigate biological mechanisms within tumors that may affect clinical outcome in stage III melanoma, we did gene expression profiling on an initial training set of 29 melanoma specimens from patients with diverse clinical outcome following lymphadenectomy for stage IIIB and IIIC melanoma. We further sought to prospectively predict clinical outcome based on a molecular profile in two independent validation sets comprising 10 and 14 stage III melanoma patients. Using this molecular information, we identified cellular pathways and networks that may be differentially regulated between the two patient groups and are possible targets for therapeutic intervention.

Specimen collection and selection for microarray analysis. Melanoma tissues from 29 patients who underwent surgical lymphadenectomy for clinically palpable nodes between 1997 and 2004 at Austin Health were selected for microarray analysis. All specimens were collected under a tissue procurement protocol approved by the Austin Health Human Research Ethics Committee and with the written informed consent of each patient. Snap-frozen specimens were embedded in OCT compound and stored as tissue blocks at −80°C within the Ludwig/Austin tissue bank repository at the time of surgery. Diagnosis was confirmed by a single dedicated pathologist in all cases.

Patient samples were selected for microarray analysis based on time taken to tumor progression (TTP) from stage III to stage IV disease and included 16 “poor-prognosis” (mean TTP, 4 mo) and 13 “good-prognosis” (mean TTP, 40 mo) patients. We chose a TTP of at least 24 mo as a cutoff for good prognosis based on the tail of the survival curves for stage IIIB/C melanoma (26) and on the availability of cryopreserved tissues. Postoperative reviews in a dedicated Melanoma Unit were carried out on a 1 to 3 monthly basis for the initial 12 mo after lymphadenectomy, followed by 3 to 6 monthly reviews thereafter or according to clinical requirement until 4 y, with annual review thereafter. Staging investigations were done according to clinical suspicion or routinely every 3 to 6 mo.

Tissues were considered acceptable for this study if minimal necrosis was present and tumor cells comprised at least 60% of the total cell population. At the time of RNA extraction, two 5-μm sections were cut and stained with H&E to ensure integrity of the extracted tissue.

RNA extraction and cDNA synthesis. cDNA synthesis and hybridization with a common reference design were conducted in duplicate with tissue samples from the 29 selected patients. Total RNA was extracted from OCT-embedded tissue by immersing and homogenizing tissue sections in Tri-reagent (Molecular Research Center). Chloroform (1.5 mL) was added to the homogenate, the sample was centrifuged, and the upper phase was removed and mixed with 100% ethanol. Purification using RNeasy columns was done according to the manufacturer's instructions (Qiagen). RNA quality was confirmed based on A_{260 nm} to A_{280 nm} ratios of absorbances and integrity inspected on formaldehyde-agarose gels against rRNA standard markers. cDNA was synthesized from 20 μg RNA in the presence of oligo(dT) and aminoallyl deoxynucleotide. Cy dyes (Amersham Biosciences) were coupled to tumor cDNA and reference cDNA produced in parallel. Reference cDNA was synthesized from the pooled RNA of a variety of tumor samples and cell lines, including melanoma, as well as from normal tissues (Supplementary Table S1).

Oligonucleotide arrays and data analysis. Oligonucleotide probes (30,888), representing individual genes and internal controls, were obtained from MWG Biotech and spotted as high-density arrays using an Omnigrid robot (Gene Machines). Labeled tumor/reference cDNA was cohybridized and scanned using a GenePix 4000A microarray scanner (Axon Instruments). The matrix overlay was aligned to the scanned image and feature extraction was done using GenePix v6.0 software (Axon Instruments). The raw data were analyzed using GeneSpring v7.2 (Silicon Genetics). The data were normalized to print-tip group and then median was normalized. Briefly, a Lowess curve was fitted to the log-intensity versus log-ratio plot. Twenty percent of the data were used to calculate the Lowess fit at each point. This curve was used to adjust the control value for each measurement. Each gene was then divided by the median of its measurements in all samples.

Data for independent validation set B from the European Organization for Research and Treatment of Cancer (EORTC) melanoma study (27) were available through the ArrayExpress public data repository.⁵

The data were imported into GeneSpring v7.2 and normalized per spot, per chip, and per gene as recommended by the chip manufacturers (Agilent Technologies). In brief, the measured intensity of each gene was divided by its control channel value in each sample and then divided by the 50th percentile of all measurements in that sample. Finally, each gene was divided by the median of its measurements in all samples. Expression values for the differentially expressed genes from the training set were used to calculate a predictive score as described below.

Gene Ontology terms were used to group differentially expressed genes based on biological function (28). Using the GeneSpring software, a Mann-Whitney test was done to compare the differentially expressed genes with existing Gene Ontology–based gene lists to identify significant representation of Gene Ontology terms.

Statistical methods. Gene expression data were first subjected to a filter that excluded probes that were not present in all samples. Of the initial 30,888 probes considered, 18,807 passed this filter and were used for ANOVA, hierarchical clustering, and principal component analysis. Differentially expressed genes were discovered by doing a Mann-Whitney test with the false discovery rate controlling method of Benjamini and Hochberg (29) used to correct for multiple testing based on a P value cutoff of 0.05. Hierarchical clustering of samples was done using average linkage and Spearman correlation as the distance function.

Data availability. The raw data have been uploaded to the ArrayExpress public data repository under the accession number E-TABM-403.⁵

Quantitative real-time PCR. Quantitative real-time PCR (qPCR) was done for differentially expressed genes to confirm the array results and then in validating the predictor using validation set A. First-strand cDNA was synthesized from 2 μg of total RNA extracted for the array experiment and 1 μg of random hexamer primer (Promega). Reverse transcriptase was omitted for negative controls. Intron-spanning multiplex assays were designed for qPCR (Supplementary Table S2) using the Universal Probe Library Assay Design Centre⁶

(Roche). Reactions were carried out in duplicate for each individual patient sample using the ABI 7700 Prism Sequence Detector (Applied Biosystems). Thermal cycler conditions were as follows: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 94°C for 20 s and 60°C for 45 s. All results were normalized to 18S amplification (Applied Biosystems). We calculated relative expression using the target threshold (C_T) value for reference as a calibrator (30).

The relative expression values for individual genes were then plotted alongside the normalized log₂ ratio array values and correlation coefficients were calculated.

Clinical and pathologic features for the patients included in the training set and validation set A are listed in Table 1. All patients had information on age at initial diagnosis, gender, and number and location of positive lymph nodal metastases. All patients had at least stage IIIB melanoma with clinically palpable lymph nodes. We used the American Joint Committee on Cancer (AJCC) staging criteria because this is the most widely used clinical staging system. Once regional lymph nodes are involved with macroscopic melanoma (stage IIIB or IIIC), the thickness of the primary tumor ceases to be a significant prognostic factor. As initial diagnoses originated from several clinics, it was not possible to ascertain whether ulceration was present in the primary melanoma in some cases. Ulceration in the primary tumor is an independent prognostic factor that, if present, upstages the disease from IIIB to IIIC (31).

In the training set, the mean TTP for the good-prognosis group was 40 months compared with 4 months in the poor-prognosis group. There were no statistically significant differences in the median age and sex between the groups, although the good-prognosis group appeared younger and contained more women. A multivariate analysis using the Cox proportional hazards method showed no statistically significant differences in other known prognostic characteristics, including AJCC staging, the use of adjuvant IFN, and the presence of tumor-infiltrating lymphocytes, although this may reflect the limited number of samples available.

One patient had isolated stage IV disease confined to resected spleen, but given that this patient remained disease-free, this sample was included. Exclusion of this sample did not alter the gene expression profile.

Differentially expressed genes segregate the two prognostic groups. Unsupervised hierarchical clustering did not reveal subgroups that correlated with prognostic or other clinical information. This was expected given the similarities between the samples. To search for genes that could effectively segregate the prognostic groups, differential gene expression was investigated. Two thousand one hundred forty genes were differentially expressed between the two groups; however, the stringent application of multiple testing correction reduced this to 21 genes with highly significant differential expression (Table 2). Hierarchical clustering and principal component analysis showed the ability of the 21 genes to segregate the prognostic groups (Fig. 1). The 21 genes were further validated in the training set using qPCR and the genes with the highest correlation coefficient between the two platforms (r > 0.5; P < 0.05) were selected for further analysis (data not shown). Of these 21 genes, 17 intron-spanning qPCR assays could be designed for well-annotated genes. The four that could not be designed were from gene prediction sequences. Fifteen genes exhibited strong cross-platform correlation and these were used in the development of a predictive score.

Gene Ontology grouping of the 2,140 differentially expressed genes showed that genes involved in apoptosis, nuclear factor-κB pathways, Wnt/Frizzled signaling, immunologic signaling, and development were highly represented. Gene clusters achieving P ≤ 0.05 are reported (Supplementary Fig. S1).

Development of predictive scores. The training set was used to develop a predictor that was subsequently tested on two independent validation sets. Two predictive algorithms were developed based on the array data and the qPCR data. (a) To calculate a predictive score from the array data (aPS), the 21 most statistically significant differentially expressed genes were used. The normalized expression ratios were transformed by raising the values to the power of two. Genes down-regulated in the good prognostic group were ascribed a negative value. The final score was then calculated by the sum of values for all 21 genes. A positive score was associated with improved outcome. (b) For the qPCR data (qPS refers to the predictor), ΔΔC_T values for the 15 most correlated genes were applied to a logistic regression algorithm, which used the Akaike Information Criterion to select only those genes that contribute to class distinction (32). This selected five highly significant genes that were then used in the following equation: qPS = [1,328.15 − 187.42 (IDH)] + [137.10 (MFG8)] + [73.61 (PILRA)] + [211.22 (HLA-E)] + [143.94 (TXNDC5)]x − 1. As with the aPS, a positive score was associated with improved clinical outcome.

The predictive scores correlate with TTP and survival. As expected, both the aPS and qPS applied to the training set were capable of distinguishing the two prognostic groups. A strong correlation between individual scores and both TTP and overall survival were evident, such that the magnitude of individual scores correlated with improved outcome for both the aPS and qPS (Fig. 2). This suggests that the expression level of these differentially expressed genes is related to underlying biological mechanisms that directly influence clinical outcome, emphasizing their prognostic relevance.

Application of the predictive scores to two independent sets. We sought to apply our results to independently generated data. We applied the five-gene qPS algorithm to an independent set of 10 tumors from the Ludwig/Austin tissue bank for which qPCR data assays were done (validation set A). The predictor correctly classified all five of the good-prognosis tumors but misclassified one of the five poor samples, yielding an overall correct classification rate of 90% (Fig. 3). The incorrectly classified poor sample represented a patient in whom TTP was brief but who had a prolonged overall survival of 6 years with metastatic disease.

We identified one published gene expression data set with a subgroup of patients similar to our own. Of the 83 patients who were profiled in this study (27), 14 had stage III disease with long-term follow-up (validation set B). In this subgroup, 10 patients would have been classified as poor (mean TTP, 10 months) and 4 patients would have been classified as good (mean TTP, 62 months) using similar criteria applied in our training set. This published data were generated using a commercially available two-color oligonucleotide microarray, an independent yet comparable platform to our own. Of the 21 differentially expressed genes from our microarray experiments, only 13 probes were able to be located for this validation set. Importantly and given their predictive power, of the 13 probes, the five genes that comprise the qPS were used in this validation set. When the aPS algorithm was applied to these samples, 9 of the 10 poor patients and 3 of the 4 good patients were correctly predicted, yielding an overall correct classification rate of 85%.

We report the successful prediction of clinical outcome in an otherwise indistinguishable group of stage III melanoma patients using an expression profile derived from microarray gene expression data and qPCR. We have shown in two independent sets that two developed predictive score algorithms, based on the expression of 21 genes, can be applied to microarray and qPCR data to prospectively forecast clinical outcome in patients with stage IIIB/C melanoma.

The data were derived from clinical groups of patients who were selected based on good and poor outcomes, although all patients had similar disease stage. Several studies had previously shown more differences in gene expression between patients with similar stage disease than among autologous samples taken at different stages from the one patient (27, 33, 34), so gene expression differences were expected; however, it was unclear whether a clear dichotomy between these groups would be found that could be used to predict outcome. Indeed, the observed gene expression differences following unsupervised clustering and principal component analysis were too subtle to segregate the patients into distinct prognostic groups based on the expression of all 30,888 probes and the key genes could only be identified following supervised analysis. These key genes were not those that have previously been described to relate to melanoma progression. However, as has been shown in other microarray studies in different tumor types, the prognostic value of these genes may be quite independent from any recognized role in the pathogenesis or progression of the tumor (35). This was achieved with a high degree of statistical significance using two completely independent data sets, although both were small. First, the profile was applied prospectively to an independent patient population drawn from our own clinic (“validation set A”). The score was able to successfully predict outcome with 90% accuracy. The profile was then used to analyze patients previously reported in the EORTC study (“independent validation set B”; ref. 27). Because the available published information was limited, these subjects were not as closely controlled in terms of end points such as TTP and survival compared with validation set A. Nonetheless, the molecular profile was able to classify 85% of the patient samples accurately, although 8 of the 21 genes were not represented in the EORTC data set. Furthermore, because the scores were quantitative, it was also possible to reliably correlate the magnitude of the score with duration of the two survival end points: TTP and overall survival.

Once melanoma has metastasized to local lymph nodes, 70% of patients will die within 5 years (31). The minority of patients with prolonged survival represent a unique cohort. No current adjuvant therapies offer an overall survival benefit, and although some clinicians offer IFN-α to improve disease-free survival (36), others offer no active adjuvant treatment outside clinical trials. Predicting which patients are likely to do well could prevent needlessly toxic adjuvant therapy. Importantly, it could greatly assist in the design of clinical trials through patient stratification by prognostic group, which will reduce both type I and type II statistical errors. This will facilitate the development of superior therapeutic strategies. The 12-year update following the Eastern Cooperative Oncology Group 1684 study and other randomized studies have shown that IFN-α improves TTP but not overall survival in stage III melanoma (5, 36, 37). It has not previously been possible to control for inherent heterogeneity within populations of patients participating in adjuvant clinical trials. Such imbalances have the potential to confound such studies, particularly smaller trials. Using molecular signatures to stratify patients may allow treatments to be compared more accurately.

Previously characterized markers of progression in melanoma have predominantly been investigated using melanoma cell lines or xenografts. Because these samples have often been passaged at length outside the original host, it is not surprising that their genetic profiles showed little overlap with our own. Similarly, the overlap in genes between our study and two other similar studies [EORTC melanoma group (27) and Mandruzzato et al. (38)] was small. There are many reasons for this lack of overlap, including differences in study design, microarray platforms, and analysis methodology. The study from the EORTC melanoma group analyzed tissue predominantly from patients with early-stage melanoma but did include 14 with stage III disease in their initial cohort of 83 patients (27). Of these 14 patients, however, few had achieved long-term survival, a key inclusion criterion for this study. The study from Mandruzzato et al. predominantly examined lymph nodes; however, they also included patients with stage IV disease. These were not stratified based on survival before the array experiment and their gene list was not validated on independent samples (38). Nonetheless, there was some overlap with our own gene set.

In both our own and in the Mandruzzato study, most human leukocyte antigens (HLA) class I and II molecules, including HLA-E, HLA-F, HLA-G, HLA-DR, and HLA-DQ, were more highly expressed in the good-prognosis tumors. HLA-E belongs to the family of class I HLA molecules that are termed “nonclassic” because they are able to bind the CD94 receptor on natural killer cells and cytotoxic T lymphocytes. Cell surface expression of HLA-E is rarely found in normal tissues, and it has been postulated that its expression provides a mechanism of tumor escape from immune surveillance (39, 40). In contradiction to this, a recent study has documented loss of HLA-E expression occurring with the progression from a primary tumor to metastatic disease (41). This finding is consistent with our own, which revealed reduced expression in the poorer prognosis tumors.

In contrast, we found the interleukin-4 receptor to be more highly expressed in the poor-prognosis tumors. Among the many differences between the two studies, it is important to note that the Mandruzzato study comprised 69% of patients who were treated with IFN-α, a molecule that may alter interleukin-4 receptor expression (42), whereas only 7% of our patients had received this cytokine (37).

Several of the other molecules that make up the 21-gene list may play a role in anticancer immune responses. For example, CHST4 is a molecule that has been found to be expressed in colon cancers (43) but, more recently, has been identified as a regulator of L-selectin ligands and lymphocyte homing, especially in Peyer's patches (44). Macrophages have been shown to secrete MFGE8, a phosphatidylserine-binding protein, which tethers apoptotic cells to macrophage integrins (45). Finally, KCNIP2 has been shown to regulate T lymphocyte proliferation, thereby resulting in immune suppression (46).

Beyond the shortlist of the 21 genes used to create the predictive scores, several other molecules that have been previously shown to play a significant role in melanoma progression were identified in the larger list of 2,140 differentially expressed genes. For example, molecules involved in Wnt signaling, nuclear factor-κB, and apoptosis pathways were shown to be up-regulated in the poor-prognosis tumors. Wnt5a is a marker of a more aggressive cellular phenotype with increased expression correlating with increased invasiveness and clinical progression of disease (47). Strategies that target this molecular pathway have been shown to increase apoptosis in melanocytes (48). Furthermore, small molecules that target nuclear factor-κB and apoptosis pathways are currently being investigated as novel therapies for the treatment of melanoma (49). Our results suggest that there is real clinical relevance in developing therapeutic strategies that abrogate pathways differentially regulated between the prognostic groups.

To our knowledge, this is the first study in which a validated predictive algorithm has been developed using oligonucleotide arrays and qPCR to forecast clinical outcome prospectively in patients with stage III melanoma. The small sample size tested requires this study to be repeated on additional patients and this is a focus of ongoing work. To be able to prospectively identify which patients will not gain additional benefit from adjuvant treatments using standard laboratory tools may help limit the use of potentially toxic therapies. Perhaps more importantly, it will aid researchers who are evaluating new adjuvant strategies in stage III melanoma through the prognostic stratification of trial subjects.

No potential conflicts of interest were disclosed.

We thank Carmel Murone and Judy Browning who coordinate the Ludwig/Austin tissue bank repository, Suzanne Svobodová who provided expertise in optimizing qPCR assays, Benjamin Cowie who provided advice and assistance with statistical methods, and Andrew Holloway for advice and assistance with microarray informatics.