Genomic and signalling pathway characterization of the NZM panel of melanoma cell lines: A valuable model for studying the impact of genetic diversity in melanoma.

Melanoma is a disease associated with a very high mutation burden and thus the possibility of a diverse range of oncogenic mechanisms that allow it to evade therapeutic interventions and the immune system. Here, we describe the characterization of a panel of 102 cell lines from metastatic melanomas (the NZM lines), including using whole-exome and RNA sequencing to analyse genetic variants and gene expression changes in a subset of this panel. Lines possessing all major melanoma genotypes were identified, and hierarchical clustering of gene expression profiles revealed four broad subgroups of cell lines. Immunogenotyping identified a range of HLA haplotypes as well as expression of neoantigens and cancer-testis antigens in the lines. Together, these characteristics make the NZM panel a valuable resource for cell-based, immunological and xenograft studies to better understand the diversity of melanoma biology and the responses of melanoma to therapeutic interventions.

The discovery of mutually exclusive mutations in BRAF (Akbani et al., 2015) and NRAS (Alexandrov et al., 2013) genes and recent findings of inactivating alterations in NF1 (Al‐Khadairi & Decock, 2019; Andersen et al., 1993) gene in melanoma together suggested the RAF‐MEK‐ERK pathway was a major driver of melanoma tumorigenicity. This resulted first in the development of drugs that specifically target V600E mutant forms of BRAF, such as vemurafenib (Busam, Hedvat, Pulitzer, von Deimling, & Jungbluth, 2013) and dabrafenib (Chatterjee et al., 2018). However, these drugs are only effective in some, but not all, tumours with BRAF mutations and can cause paradoxical MAPK pathway activation in tumours with wild‐type BRAF (Colombino, 2012; Davies et al., 2002; Dreno et al., 2018). The initial clinical success of these targeted therapies in treating BRAF‐mutant tumours has been tempered by the fact that a number of tumours display inherent resistance and that most of those that respond initially go on to develop adaptive or acquired resistance quite rapidly (Hatzivassiliou et al., 2010). The efficacy of BRAF inhibitor treatment is improved by the addition of MEK inhibitors, but again not all tumours respond and resistance still develops (Hatzivassiliou et al., 2010). Furthermore, a majority of advanced melanoma patients failed to respond to immune checkpoint blockade therapies. Therefore, understanding why some tumours do not respond to these therapies while others do has been the subject of intensive research (Hauschild et al., 2012; Heidorn et al., 2010). Particularly, the high levels of mutation burden, genetic instability and marked heterogeneity of gene expression patterns in melanoma are likely to be contributing factors to both initial drug response and intrinsic or acquired resistance (Hélias‐Rodzewicz et al., 2015). Currently, cultured cell models of melanoma are widely used to study mechanisms and drug responses. However, due to the genetic heterogeneity of melanoma cell lines, it would be advantageous to have access to large panels of well‐characterized cell lines that are representative of the multiplicity of driver mechanisms and varied patterns of gene expression seen in the disease. We have established a panel of NZM cell lines cultured from melanoma patient samples (Hélias‐Rodzewicz et al., 2017; Henare et al., 2012), some of which have been used previously to study melanoma function (Huang et al., 2013; Improta et al., 2013; Jeffs et al., 2009; Kakadia et al., 2018). Here, we present a more comprehensive genomic and gene expression characterization of the cell panel using Sequenom MassARRAY, whole‐exome and RNA sequencing to allow wider use of these lines as a resource for studying the impact of genetic diversity in melanoma on therapeutic responses.

In total, 102 NZM melanoma cell lines were prepared from biopsies of metastatic melanoma samples from patients presenting at clinics in Auckland and Palmerston North, New Zealand. All patients gave appropriate written informed consent, as previously described (Improta et al., 2013). At the time of sample collection, the standard of care for melanoma in New Zealand was surgery and chemotherapy. Therefore, most samples collected were from drug‐naïve patients. The preparation and use of the cell lines were conducted in accordance with protocols approved by the Northern Region Health and Disability Ethics Committee (AKL/2000/184/AM04). Four of the lines were repeated biopsies from later surgeries, so lines NZM1 and 2; NZM41&47; NZM42&48; NZM62&64 are pairs of lines from the same subject. The cells were passaged in α‐modified minimal essential medium (MEM‐α) supplemented with penicillin (100 U/ml), streptomycin (100 µg/ml), amphotericin B (0.25 µg/ml; GIBCO Life Technologies), insulin (5 μg/ml), transferrin (5 μg/ml), sodium selenite (5 ng/ml; Roche Diagnostics GmbH) and 5% foetal bovine serum (FBS). Short tandem repeats (STR) analysis was performed for cell line authentication (Table S1). Importantly, cells were maintained in a 5% oxygen environment at all times to minimize changes in cell characteristics caused by growth in non‐physiological atmospheric conditions used in standard cell culture.

First, the NZM cell lines were genotyped using Agena Bioscience Oncocarta V1.0 and MelaCarta MassARRAY panels, which, respectively, interrogated 238 mutations in 19 genes and 72 mutations in 20 genes (see Supplementary Methods). The results of this genotyping are summarized in Figure 1 and Table S2. Next, whole‐exome sequencing was performed on a sub‐panel of 52 NZM cell lines using an Ion Proton next‐generation sequencing platform (Thermo Fisher) and the manufacturer's protocol (see Supplementary Methods). Raw exome sequence data can be downloaded using the NCBI accession numbers in Table S3. Raw data analysis, alignment and variant calling were performed using a Torrent Suite Software and Ion Reporter Software v5.6 (Thermo Fisher). A summary of the driving mutations of BRAF, NRAS and NF1 identified in the NZM lines is shown in Figure 2 and Table S4. All but one of the variants identified in mass array analysis were validated by sequencing except in NZM3, where the mass array indicated BRAF V600K while whole‐exome sequencing indicated V600E. This analysis also revealed potentially functional NF1 variants occur in 7 of 52 lines (Figure 2 and Table S4). Furthermore, whole‐exome sequencing data also allowed a more detailed understanding of genetic variations in tumour suppressor genes. Reduction in copy numbers of tumour suppressor genes was common in the NZM lines (Figure 3a and Table S5). In particular, homozygous deletion of CDKN2A was found in all major genotype groups, while PTEN deletion was mainly found in the BRAF‐mutant group (Figure 3b). In contrast, TP53 mutations and deletions are rare in the NZM cells as is the case in melanoma tumours. Furthermore, expression patterns of key signalling molecules were also profiled by Western blotting (Figure S1).

Figure 1

Mutational landscape of the panel of 102 NZM cell lines as divided by oncogenic mutations of BRAF and NRAS. DNA extracted from 102 NZM cell lines were genotyped by Sequenom analysis for hotspot mutations commonly found in melanomas. Lists of those mutations were presented in Supplementary Methods. The top row indicates major genotype groups based on BRAF and NRAS driver mutations. The second row indicates the number of hotspot mutations per cell line. Middle rows indicate colour‐coded individual mutations found in 102 NZM cell lines. Side panels indicate the percentage of cell lines with hotspot mutations per gene. The bottom row indicates cell line names

Figure 2

Mutational landscape of NZM cell lines as divided by mutations of BRAF and NRAS and NF1. Whole‐exome sequencing was performed in 52 NZM cell lines. The top row indicates major genotype groups based on BRAF and NRAS driver mutations, of which the double‐negative group was further divided into the NF1‐mutant group and the triple‐negative group. The second row indicates the number of hotspot mutations per cell line. Middle rows indicate colour‐coded individual mutations found in 52 NZM cell lines. Side panels indicate the percentage of cell lines with mutations per gene. The bottom row indicates cell line names

Figure 3

Deletion of tumour suppression genes in NZM lines. The deletion of tumour suppression genes (TSGs) was identified using whole‐exome sequencing data of 52 NZM cell lines. (a) Percentages of the NZM lines with TSGs deletion. (b) Distribution of homozygous deletion of PTEN and CDKN2A in the mutational landscape of NZM lines. The top row indicates the number of mutations of BRAF, NRAS, PTEN and CDKN2A genes per cell line. Middle rows indicate colour‐coded individual mutations of BRAF, NRAS, PTEN and CDKN2A genes found in 52 NZM cell lines. Side panels indicate the percentage of cell lines with mutations per gene. The bottom row indicates cell line names

The overall frequency of major melanoma‐associated genetic alterations in the NZM cell lines studied is summarized in Figure 4a. This analysis revealed the ratio of major known melanoma mutations in the NZM cell line panel is very similar to that described in primary melanoma tumours (Lim, Menzies, & Rizos, 2017; Maertens et al., 2013). A total of 44 lines of the panel (43%) had mutations in BRAF, of which 38 lines were V600E mutant, 3 lines were V600K mutant and one each of V600R, G466A and G469A mutations. A further 23 lines (23%) had NRAS mutations, of which 21 lines were Q61 mutant. A further 2 lines had known KRAS‐activating mutations. As expected, RAS and BRAF mutations were mutually exclusive in the lines. Notably for the major driver mutations, the V600 mutations in BRAF and Q61 mutation in NRAS were on average heterozygous. In contrast, other variants in these genes tended to be homozygous (Figure 4b,c). These patterns were similar to the zygosity status of BRAF V600 and NRAS Q61 observed in melanoma tumours in previous studies (Marshall et al., 1993; Marshall et al., 1992; Mitsiades et al., 2011; Poulikakos, Zhang, Bollag, Shokat, & Rosen, 2010). Other cell lines representing rarer melanoma genotypes were also observed. For example, two of the lines had GNA mutations, but not NRAS or BRAF mutations, and so were presumably of uveal origin (Reifenberger et al., 2004). Activating mutations in PK3CA (4 lines), PDGFRA (2 lines), CTNNB1 (3 lines) and MET (1 line) were also observed. In addition, the mass array analysis revealed 60/102 of the cell lines contained C228T or C250T mutations in the TERT promoter (Figure 5a), which is a similar ratio to that previously described in melanoma tumours (Stagni et al., 2018). Interestingly, we also note that the C228T variant was more commonly seen in heterozygous form compared to the C250T variant (Figure 5b).

Figure 4

Comparison of % of cell lines or tumour samples containing various genetic alterations. (a) Comparisons are with published data from TCGA (Cancer Genome Atlas, 2015) except for TERT promoter mutation frequency incidence (Huang et al., 2013). (b) Comparison of allele frequency of BRAF and NRAS mutations. Genotype and allele frequency were assessed using Sequenom analysis as described in Supplementary Methods. Statistical difference was analysed by a two‐tailed t test with ** and *** indicates p < .01 and p < .001, respectively

Figure 5

TERT promoter mutations in NZM lines. DNA extracted from NZM cell lines was genotyped by Sequenom analysis for hotspot mutations, including the 2 TERT promoter mutations, C228T and C250T. (a) Oncoprint plot of TERT promoter mutations in alignment with BRAF and NRAS driver mutations in NZM cell lines. (b) Frequency of mutant alleles of TERT C228T and TERT C250T mutations. Statistical difference was analysed by a two‐tailed t test with **** indicates p < .0001

Next, we performed RNA sequencing analysis for a panel of 28 NZM lines (see Supplementary Methods). Raw RNA sequence data can be downloaded using the NCBI accession numbers in the Table S6. Variant allele frequencies were highly correlated between the DNA and RNA sequence data (Figure S2). Hierarchical clustering of gene expression profiles revealed four broad subgroups of cell lines (Figure 6), of which differences in the major pathways corresponded to interleukins (C1), extracellular matrix proteins (ECM; C2), cell cycle (C3) and stress response (C4; Figure 6 and Figure S3). We noted that BRAF and NRAS mutations randomly distributed between the four subgroups of cell lines (Figure 6), consistent with similar descriptions of discordance between gene expression and mutation status noted in previous reports (Jeffs et al., 2009; Sweetlove et al., 2015; Tsai et al., 2008). We further characterized the expression patterns of specific genes of interest across the subgroups (MITF, SOX10, SOX9, SMAD3, CTNNB1, AXL, NGFR, EGFR, ERBB3, YBX1, EBF3 and PXDN; Figure S4). Overall, except for cluster C1, clusters C2 to C4 corresponded broadly in gene expression patterns to neural crest‐like (low MITF expression signature), melanocytic (high MITF expression signature) and transitory clusters, respectively, identified in earlier studies (Jeffs et al., 2009; Tsoi et al., 2018).

Figure 6

Transcriptionally defined subgroups of NZM lines. RNA extracted from 28 NZM lines was sequenced, and the data were used for hierarchical clustering analysis. Four major clusters annotated C1 to C4 were identified. Pathway analysis was performed on the top 500 expressed ENSEMBL genes (RPKM) of each cluster against Reactome database, which revealed differences in the major pathways corresponding to interleukins (C1), extracellular matrix proteins (ECM; C2), cell cycle (C3) and stress response (C4; more details about pathway analysis can be found in Figure S3)

Finally, we performed immunogenotyping in the NZM lines to facilitate the use of these cell lines for immuno‐oncology research (see Supplementary Methods). Firstly, for each NZM line, we identified up to two class I HLA haplotypes in both DNA and RNA sequence data with generally good concordance between RNA‐based and DNA‐based HLA haplotype identification (Table S6). Then, neoantigens were predicted by combining: (a) DNA and RNA sequence variants, (b) the ability of the cells to present peptides encoded by these variants to the immune system, given their HLA haplotypes, and (c) expression of the variants (Figure 7a,b). Cancer–testis antigens (CTA) are not generally expressed in somatic cells but can be recognized by the immune system when expressed in tumours (Welsh, Rizos, Scolyer, & Long, 2016); some CTAs such as MAGE‐A3 have been used in immunotherapeutic vaccine trials in melanoma (Wilmott et al., 2013). The expression of RNAs encoding CTA was quantified in these cell lines to facilitate their future use in CTA research. For example, NZM104 and NZM3 co‐expressed detectable levels of multiple CAGE and MAGE family RNAs, respectively (Figure S5).

Figure 7

Neoantigen identification in NZM lines. Neoantigens were predicted by combining DNA and RNA sequence variants, the ability of the cells to present peptides encoded by these variants to the immune system and the expression of the variants (more details can be found in the Supplementary Methods). (a) Number of expressed neoantigens in NZM lines. (b) Correlation between variant number and expressed neoantigen number

In summary, the NZM early passage primary melanoma cell panel has already proved valuable in melanoma research. Here, we have shown that these lines represent a wide range of melanoma genotypes, signalling pathway subsets inferred from RNA expression profiles and immunogenic features. Given that all lines in the panel are directly comparable to one another since they were derived and maintained identically, they provide a uniquely useful resource in which to study the impact of the genetic, gene expression and immunogenic diversity found in melanoma. These cell lines can be made available upon reasonable request for non‐commercial research to investigators able to cover the handling and shipping costs involved.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Fig S4

Click here for additional data file.

Fig S5

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.

Table S6

Click here for additional data file.

Supplementary Methods

Click here for additional data file.