Intratumoral genetic heterogeneity in metastatic melanoma is accompanied by variation in malignant behaviors

Regional DNA copy number heterogeneity in metastatic melanoma

Eight regions of lymph node metastasis Tumor 1 were assessed for the presence of chromosomal amplifications and deletions. DNA extracted from cores taken from three separate FFPE tissue blocks was analyzed using the Affymetrix Oncoscan 2.0 platform. H&E staining was used to identify regions composed primarily of tumor cells prior to coring, with sections taken from immediately below analyzed fragments to control for contaminating normal tissue (Figure 1A and Additional file 1: Figure S1). Hierarchical clustering of DNA copy number profiles separated the samples into two groups, with visual inspection of the heatmap showing that cores taken from the same tissue block often demonstrated very different patterns of amplifications and deletions (Figure 1B). Statistically significant regions of amplification and deletion were next defined using a segmentation algorithm, and the occurrence of specific aberrations compared across the tumor regions. The sampled tumor regions harbored between 44 and 133 significant regions of copy number changes (Figure 2A), encompassing between 23 and 59 percent of the genome (Figure 2B). The greatest proportion of changes was present in all regions; however, many aberrations were present in only one or two cores (Figure 2C). Heterogeneity was observed in genomic regions containing genes with demonstrated potential to impact melanoma biology, such as the high level amplification (greater than 5 copies) of chromosome band 1q21 observed in Core 2 from Block 1–2. This region encompasses the gene for histone methyltransferase SETDB1, recently identified as an oncogene [12] and a candidate susceptibility gene [13] in melanoma. Detailed probe level and segmentation results from Chromosome 1 and Chromosome 17 are shown in Figure 3 and Additional file 2: Figure S2 respectively.

The Oncoscan 2.0 platform also includes probes that test for 541 individual mutations in 62 well known cancer genes, with the presence of a mutation indicated by high probe intensity. Despite the presence of mutation-specific probes with highly heterogeneous intensity, no heterogeneous mutations could be validated by capillary sequencing. Similar results were obtained following Ion Torrent sequencing of amplicons covering a similar panel of cancer genes, which was performed in regions from two additional melanomas (Tumors 2 and 3). These experiments are described in the Additional file 3, in Additional file 4: Table S1, and in Additional file 5: Figures S3 and Additional file 6: Figure S4.

Early passage metastatic melanoma cell lines retain genetic heterogeneity

The copy number profiles of a cell line derived from Tumor 1 (LM-MEL-62), ten single cell clones derived from LM-MEL-62, and a fresh frozen tumor fragment of Tumor 1 were assessed using Illumina SNP microarrays. All cell line-derived samples clearly shared many common copy number changes with the original tumor, such as gain of 6p (Figure 4A). The clones harbored between 55 and 69 copy number changes, which represents significantly less variation than was observed in the different regions of Tumor 1 (Figure 4B). Nevertheless, the cell clones were also heterogeneous for many of the detected aberrations (Figure 4C), and the LM-MEL-62 clones displayed heterogeneity at chromosome regions similarly affected in the archival FFPE tumor material (Figure 5A and B). This suggests that early passage melanoma cell lines are polyclonal, and that they can retain or recapitulate genetic heterogeneity representative of that found in the patient’s tumor.

We generated copy number profiles for clones from two other cell lines derived from Tumor 2 and Tumor 3 (Figure 6A & B respectively). Clones from both cell lines had copy number changes that were identified in all clones, and some that were present in some clones but not others (Figure 6C & D). The clones from the three cell lines assessed contained significantly different quantities of copy number changes (one way ANOVA, p < 0.0001; Figure 6E), however they did not differ in the proportion of the genome affected by copy number abnormalities (one way ANOVA, p = 0.165; Figure 6F).

The copy number profile of the parental cell lines represent an average signal derived from all cells in the line. Therefore, the clones with profile closest to the parental line as determined by clustering (such as Clones 7, 8, and 9 in LM-MEL-62) represent the most prevalent or dominant clone(s), with the others representing relatively minor populations. This information provided the opportunity to assess functional differences between dominant and minor cell populations.

Single cell clones from a metastatic melanoma cell line are functionally heterogeneous

In order to assess whether the genetic heterogeneity observed in cell lines was accompanied by functional variation, we compared clones from LM-MEL-62 to the parental line. Clear differences were seen in their sensitivity to cytotoxic drugs paclitaxel (Figure 7A) and 5-fluorouracil (5FU) (Figure 7B), and in their ability to form adherent colonies from single cells (Figure 7C), and in soft-agar (Figure 7D). Clones with copy number profiles that did not cluster with the parental line (presumably less prevalent in the pooled population) were those that demonstrated the strongest behaviors. For example, Clone 3 was the most clonogenic, and Clones 1 and 10 were the least sensitive to cytotoxic drugs. Clone 2 was not assessed for drug sensitivity, and Clone 10 was not assessed in soft-agar assays, as they ceased proliferating before the assays could be performed. This indicates that the clones also differed in their long-term replicative potential following isolation.

Single cell clones from metastatic melanoma cell lines display evidence of differential pathway activation based on gene expression profiling

The gene expression profiles of LM-MEL-62 and derived clones were assessed using Illumina HT-12 microarrays. Unsupervised hierarchical clustering placed some clones into different clusters than observed based on copy number data, but Clones 7, 8, and 9 still clustered with the parental cell line (Figure 8A). Again as the profile of the parental line represents an average of all cells/clones, this supports Clones 7, 8, & 9 as the most prevalent cell types.

Notably, the expression of interferon-inducible genes has been associated with decreased sensitivity to paclitaxel [14, 15]. Consistent with this, the clone expressing the highest levels of inflammation and interferon-inducible genes (Clone 1; Figure 8B) was significantly less sensitive to paclitaxel than the parental cell line (Figure 7A). Treating parental LM-MEL-62 with increasing concentrations of paclitaxel likewise enriched for a cell population that expressed higher levels of several of interferon inducible genes from the enriched gene sets (Figure 8C). This demonstrates that a phenotype that was present in only a minority of cells has the potential to play an important role in survival and re-growth during and after drug treatment.

Tumor material and cell lines

The FFPE blocks from an individual tumor represent contiguous transverse slices; however, tissue orientation was not recorded during embedding. Hematoxylin and eosin (H&E) staining was used to determine regions of viable tumor cells in each tissue block, and 3mm cores were then removed from target regions using a tissue microarrayer. The tip (~2 mm) of the cylindrical core was removed with a sterile scalpel blade and used for DNA extraction. For the cores from Tumor 1 used for DNA copy number analysis, the remaining core was embedded into a recipient block of paraffin so that the upper surface could be sectioned and the proportion of tumor cells analysed. Cores from Tumors 1 and 2 were used only for sequencing-based mutation profiling.

The melanoma cell line establishment, culture methods, and RF10 growth media formulation used by our laboratory have previously been reported [20]. Single cell-derived clonal sublines (clones) were isolated through low-density plating and colony isolation using 5 mm plastic cylinders.

Microarray analysis

Nucleic acids were extracted from cell line pellets and fresh frozen tumor pieces using the AllPrep Mini Kit (Qiagen). All extractions for cell line clones were performed before the clones had been passaged five times in culture. DNA extraction from FFPE samples used the Arcturus Picopure DNA Extraction kit (Applied Biosystems) with a 24 hour Proteinase K incubation, followed by further purification using DNeasy columns (Qiagen).

DNA from cell lines, patient blood, and fresh frozen materials was analyzed on Illumina Human610-Quad genotyping arrays at the Memorial Sloan-Kettering Cancer Center Genomics Core Laboratory and imported into Partek Genomics Suite (PGS). Data from blood samples was used to create paired copy number data for each patient. Segmentation algorithm settings were: minimum markers 15, p-value 0.0001, expected range 0.5, amplification signal-to-noise ratio 0.4, deletion signal-to-noise ratio 0.8.

DNA from FFPE samples was analysed using the Oncoscan Express 2.0 service from Affymetrix, which employs arrays containing 334 000 copy number probes, and 541 probes specific for somatic cancer mutations. Oncoscan copy number data were processed and normalized by Affymetrix according to previously published methods [21], and copy number is calculated in reference to Oncoscan data from an Affymetrix panel of normal reference samples. Segmentation algorithm settings were: minimum markers 10, p-value 0.0001, expected range 0.5, amplification signal-to-noise ratio 0.4, deletion signal-to-noise ratio 1.0.

Hierarchical clustering of copy number data used Euclidian distance and average linkage.

Illumina HT-12 gene expression arrays were processed at the Australian Genome Research Facility. Using R and the Bioconductor package Limma [22], raw data were background corrected using the normexp function, log-transformed, and quantile normalized. Hierarchical clustering in Partek GS used Euclidian distance and average linkage. Gene Set Enrichment Analysis (GSEA) [23] employed gene set permutation, a 5% FDR cutoff, and gene sets from the MSigDB database v3.0, category C2.

The Illumina HT-12 gene expression array and Illumina 610-Quad SNP array data are available in the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-1753.

Mutation profiling by amplicon sequencing

Somatic mutation profiling using the Ion AmpliSeq Cancer Panel, Ion Torrent sequencing, and the Ion Variant Caller (Life Technologies) was performed by AIT Biotech (Singapore). Libraries were barcoded using the Ion Xpress Barcode Kit and eight samples were sequencing together in a single Ion Torrent sequencing run using Ion 318 chips (all Life Technologies). Default Ion Variant Caller settings were employed, with the exception of a SNP QV minimum of 14. This value was chosen as it eliminated false negatives in the identification of BRAF mutations in Tumor 1, which was previously confirmed via capillary sequencing. Known SNPs and synonymous changes were removed, as were low-confidence insertions and deletions associated with homopolymer runs were also disregarded due to known issues with false positives using this sequencing platform [24]. Variants were analyzed using the Ensemble Variant Effect Predictor [25], which included SIFT [26] and PolyPhen [27] predictions of the consequence of a sequence variant on protein function. The amplicon sequencing reads and variant prediction results are available in the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-1794.

Validation of high frequency variants was performed using standard capillary sequencing. PCR products covering mutations indentified by Ion Torrent sequencing were amplified using GoTaq mastermix (Promega). PCR products were sequenced at AGRF, using a 3730xl sequencer and BigDye Terminator v3.1 chemistry (both Applied Biosystems). Primers (all shown 5′ to 3′) used to validate the BRAF and NRAS mutations in Tumor 1 were as follows: BRAF G649E Forward Primer – TACCATGCCACTTTCCCTTG, Reverse Primer TTTTCTGTTTGGCTTGACTTGA; NRAS G12S Forward Primer – GGTTTCCAACAGGTTCTTGC, Reverse Primer CTCACCTCTATGGTGGGATCA.

Validation for low frequency variants was performed by high-resolution melt (HRM) analysis on a Rotor-Gene Q (Qiagen), before and after uracil-DNA glycosylase (UDG) treatment to control for FFPE sequence artifacts as previously described [28]. The HRM primer sequences for FGFR3 exon 6 variant analysis were F 5′-CAGTGGCGGTGGTGGTGAGG-3′ and R 5′-ACCTTGCAGTGGAACTCCACGTC-3′. PCR cycling and HRM was performed on the Rotor-Gene Q (Qiagen, Hilden, Germany). The reaction mixture in a final volume of 20 μL was prepared as follows; 1 × PCR buffer, 2.5 mM MgCl2, 400 nM of each primer, 10 ng of FFPE DNA, 200 μM of dNTPs, 5 μM of SYTO9 (Invitrogen), and 0.5 U of HotStarTaq polymerase (Qiagen). The PCR cycling and melting conditions were as follows; an initial incubation at 95°C for 15 min, followed by 55 cycles of 96°C for 15 s, 70°C for 20 s, 72°C for 30 s; one cycle of 97°C for 1 min and a melt from 70°C to 95°C rising 0.2°C per second. All samples were tested in duplicate. For UDG treatment, 0.5 × UDG buffer and 0.5 U of uracil-DNA glycosylase (New England Biolabs) were added to the PCR master mix. The same PCR conditions were used except an addition of initial incubation at 37°C for 30 min before the activation of HotStarTaq polymerase.

A PIK3CA E549D mutation detected at a low frequency (4%) in the Core 2 of Tumor 2 by AmpliSeq was analysed by limited copy number (LCN)-HRM [29]. The reaction mixture and cycling conditions were those used for the FGFR3 HRM assay with UDG treatment, except with an annealing temperature of 60°C and 100 pg of template were used. The Core 2 sample was tested in 100 replicates, but no variants were identified. The primer sequences for PIK3CA exon9 analysis were F 5′- AAGAACAGCTCAAAGCAATTTCTACAC-3′ and R 5′-AATCTCCATTTTAGCACTTACCTGTGAC-3′.

Functional assays

Single cell colony formation assays involved seeding single cells in 96 well plates and counting the wells with colonies of greater than 50 cells after four weeks.

Soft agar colony formation assays were performed in 24-well plates, with 2000 cells cultured in 500 μL of RF10 in 0.9% agarose, on a base layer 1.4% agarose in RF10. 500 μL of RF10 was layered over top of the agarose, and plates were incubated for four weeks at 37°C. Colonies were visualized and counted follow staining with MTT (3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide; Sigma Aldrich) for 60 minutes.

Paclitaxel and 5FU (Sigma Aldrich) were resuspended in DMSO and sterile water respectively. Drug treatments were applied to 7500 cells/well in 100 μL RF10 in 96-well plates, with final concentrations of 20 nM paclitaxel and 200 μM 5FU. After 72 hours the amount of viable cells remaining was quantified using the CellTiter AQueous One Solution Cell Proliferation Assay (Promega). Values from drug treatments were compared to vehicle-only controls. A paclitaxel-resistant derivative of LM-MEL-62 was created by exposing the parental cell line to drug for 72 hours, then allowing the cells to recover to confluence before re-treating. Three rounds were performed, with doses doubled each time, with a final treatment of 40 nM. All functional assays were performed on cell line clones before they had reached ten passages in culture.

Quantitative PCR (QPCR)

cDNA was produced from RNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). QPCR was performed using the QuantiFast SYBR Green PCR Kit (Qiagen), with reactions run on a Stratagene MX3005P thermocycler. Expression levels were normalized to beta-actin. Primer sequences were as follows: IFI6 F 5′- AGCAGCAGGTAGCACAAGAA-3′ and R 5′- GGGCTGAAGATTGCTTCTCTT -3′; IFI27 F 5′-GCCACAACTCCTCCAATCAC-3′ and R 5′- ATCAGCAGTGACCAGTGTGG -3′; MX1 F 5′- GATGATCAAAGGGATGTGGC -3′ and R 5′- AGCTCGGCAACAGACTCTTC -3′.