Comprehensive molecular profiling of UV-induced metastatic melanoma in Nme1/Nme2-deficient mice reveals novel markers of survival in human patients.

Hepatocyte growth factor-overexpressing mice that harbor a deletion of the Ink4a/p16 locus (HP mice) form melanomas with low metastatic potential in response to UV irradiation. Here we report that these tumors become highly metastatic following hemizygous deletion of the Nme1 and Nme2 metastasis suppressor genes (HPN mice). Whole-genome sequencing of melanomas from HPN mice revealed a striking increase in lung metastatic activity that is associated with missense mutations in eight signature genes (Arhgap35, Atp8b4, Brca1, Ift172, Kif21b, Nckap5, Pcdha2, and Zfp869). RNA-seq analysis of transcriptomes from HP and HPN primary melanomas identified a 32-gene signature (HPN lung metastasis signature) for which decreased expression is strongly associated with lung metastatic potential. Analysis of transcriptome data from The Cancer Genome Atlas revealed expression profiles of these genes that predict improved survival of patients with cutaneous or uveal melanoma. Silencing of three representative HPN lung metastasis signature genes (ARRDC3, NYNRIN, RND3) in human melanoma cells resulted in increased invasive activity, consistent with roles for these genes as mediators of the metastasis suppressor function of NME1 and NME2. In conclusion, our studies have identified a family of genes that mediate suppression of melanoma lung metastasis, and which may serve as prognostic markers and/or therapeutic targets for clinical management of metastatic melanoma.

Introduction

Metastatic spread is the primary cause of death in melanoma patients[57]. Thus, new therapeutic strategies and prognostic markers are needed that characterize the metastatic potential of primary tumors. Metastasis suppressor genes inhibit metastatic activity of cancer cell lines in vitro and in vivo[58]. NME1, the first metastasis suppressor described[63], attenuates metastasis in melanoma[75], breast carcinoma[34] and other cancers[20]. Ablation of mouse Nme1 or Nme2 increases metastasis in chemical-induced hepatocarcinogenesis[9] and UV-induced melanoma[45]. NME1 has been proposed to inhibit cell motility-driving signaling pathways[39] and induces dramatic transcriptomic alteration in breast carcinoma[23] and melanoma[32] cells that is mediated by target genes such as EDG2[24] and ITGB3[33]. NME1 may also suppress cancer progression by enhancing repair of ultraviolet light (UV) radiation-induced mutations[27, 72] and regulating double-strand break repair[27, 56, 71]. Using a genetically engineered mouse strain that is susceptible to UV-induced melanoma due to overexpression of hepatocyte growth factor and a homozygous deletion of the 16/Ink4α locus (HP mice)[51], we recently localized metastasis suppressor activity to the Nme1 and Nme2 loci[45]. While melanomas in HP mice display histological features of human melanoma, their metastatic potential is low[44, 51]. In contrast, ablation of Nme1 or Nme2 in HP mice results in robust metastasis to the lymph nodes and lung, suggesting this model is ideal for studying changes in gene mutation and expression associated with melanoma metastasis. Unlike other models of metastatic melanoma, tumors in HP mice are initiated by UV, a well-recognized human melanoma risk factor[18]. Importantly, the HP model often a yields single melanoma per mouse, enabling profiling in matched primary and metastatic tumors.

We employed whole genome sequencing and RNA-seq transcriptomic analysis to obtain molecular profiles comparing primary and metastatic tumors in HP mice harboring hemizygous deletions of both the Nme1 and Nme2 loci[50]. Here we identify novel mutations and gene expression changes associated with metastatic potential of UV-induced melanoma in mice, and show that forced expression of a representative set of these genes in human melanoma cells suppresses their invasive activity. Of particular relevance to human melanoma, patterns of expression of these genes predict survival in patients with cutaneous and uveal melanoma indexed in The Cancer Genome Atlas (TCGA).

Results

Hemizygous Ablation of the Nme1/Nme2 Locus Confers Strong Metastatic Activity in a Mouse Model of UV-Induced Melanoma

A hemizygous-null deletion spanning the Nme1/Nme2 (Nme1/2 +/Δ) locus (Fig. 1a) was introduced into the HP strain (Fig. 1b, Supplementary Fig. 1a–c) to create HPN mice. Melanoma was initiated in HP and HPN mice by exposure to UV (9 kJ/M2) at postnatal day 4 (Fig. 1c). Melanoma incidence was high in both HP (91%) and HPN (87%) mice at day 220 (Supplementary Table 1a), while the median number of melanomas per mouse was low in both groups (Supplementary Table 1b). The melanomas were strongly pigmented, with those achieving critical volume (500 mm3) acquiring a domed shape and central necrosis (Fig. 1d, Supplementary Fig. 1d). Melanomas were predominantly located on dorsal, UV-exposed surfaces, but occasionally on limbs and tail. Melanotic melanoma cells of both epithelioid and spindle/dendritic morphology were observed in the dermal compartment and they often invaded the muscle and other subcutaneous tissues (Fig. 1e, Supplementary Table 1c). Median time of melanoma onset after UV exposure was not significantly different between HP (136 days) and HPN (114d) mice (Supplementary Table 1b) and HP and HPN mice exhibited no differences across five different indices of tumor growth (Supplementary Fig. 2, Supplementary Table 1b). Melanoma incidence was higher in males than females (Supplementary Table 1a) but sex had no significant effect on melanoma tumor growth (Supplementary Table 1d).

Fig. 1

Hemizygous deletion of the Nme1/Nme2 locus increases metastatic activity of UVR-induced melanomas in the HP mouse model.

a Genomic organization of the mouse Nme1 and Nme2 loci and their ablation in transgenic mouse strains. Exons are represented by boxes (blue for Nme1, green for Nme2) and introns/flanking sequences by a black horizontal line. Above, locus of wild-type Nme1 and Nme2 (Nme1/2). Below, Nme1/Nme2 locus harboring a single large deletion (Δ) of all coding regions for both genes as shown (Δ Nme1/2). Deleted sequence is highlighted within a grey box. b Diagram depicts the two strains of mice employed in the study. A mouse strain engineered for overexpression of HGF (HGF) and homozygous deletion of the Ink4a/p16 locus (p16) designated “HP” was used to create a novel HP hybrid strain harboring the ΔNme1/2 deletion in a hemizygous-null condition (“HPN”). c Protocol for UVR initiation and monitoring of melanomas in HP mouse strains. d Representative images of primary melanomas, lungs and lymph nodes from HP and HPN mice. Primary and metastatic melanoma tissues were stored at −80°C following placement in either RNA stabilization solution (RNAlater, QIAGEN) or O.C.T. compound, or at room temperature in 10% formalin-PBS for histopathological analysis. Numbers following HP/HPN genotype indicate mouse identifiers (Supplementary Table 1c). To assess metastatic infiltration of lymph nodes, left and right nodes (cervical, brachial, axillary and inguinal) of melanoma-positive mice were dissected and their volumes derived (π/6 × L × W2). Node volumes were scored in decades of increasing size (1–10 mm3 = 1, 11–20 mm3 = 2, 21– 30 mm3 = 3, etc.). A composite lymph node enlargement score was derived for each mouse as the sum of scores for all lymph nodes. Metastatic growth within lung and liver were quantified in terms of number and size of visible pigmented colonies. Metastases were counted and measured, sorted into four bins of increasing diameter, and the respective bin counts weighted accordingly (< 0.5 mm, 1X; 0.6– 1mm, 2X; 1– 5 mm, 3X; > 5 mm, 4X). e Summary of lung metastasis and lymph node enlargement scores for all melanoma-positive HP and HPN mice. f Summary of lung metastasis and lymph node enlargement scores for all HP and HPN mice subjected to molecular profiling. g Microscopic images of hematoxylin/eosin-stained tissue sections from representative primary melanomas and metastatic tissues of HPN mice. Scale bars are provided in the lower right corner of each image.

In contrast, the incidence of lung metastasis was markedly higher in HPN (67%) versus HP mice (25%) (Fig. 1f, Supplementary Table 1a). Incidence of lung metastasis was significantly higher in females than males, although HPN males exhibited a much higher incidence of lung metastasis than did HP males. Incidence of lymph node enlargement did not differ between sexes. Lung metastasis and lymph node enlargement scores were correlated with each other across all HP and HPN mice (Supplementary Table 1e) and were significantly higher in HPN mice (Fig. 1f, Supplementary Table 1b). Lung metastasis scores were higher in females than males, while lymph node scores were not different (Supplementary Table 1d). Together, these findings demonstrate a robust increase in UV induced metastatic melanoma in HPN mice.

An eight-gene missense mutation signature in HPN melanomas is correlated with lung metastasis

Primary and metastatic (lung, lymph node) melanoma tissues from HP (4 male, 5 female) and HPN (5 male, 3 female) mice (Supplementary Table 1c) were selected for whole genome sequencing. The median lung metastasis score was significantly higher in this subset of HPN mice (P = 0.04), while the median lymph node score was marginally higher (P = 0.07) (Fig. 1g, and Supplementary Table 1f). No effect of sex on metastatic activity was observed (Supplementary Table 1g). Lung metastasis and lymph node scores were correlated within the subset chosen for profiling (Supplementary Table 1h). Whole genome sequencing was performed only on metastatic lung and lymph node specimens of the HPN genotype (Supplementary Table 1c) due to higher availability of matched primary and metastatic melanoma tissue in this strain. Normal lung and/or liver tissue from melanoma-positive HP mice and HPN mice provided germline sequence information concordant with that of the parental C57BL/6 strain.

Primary tumors from both HP and HPN mice harbored non-silent, mostly missense, mutations (Fig. 2a) which were present at an equal rate in lung metastatic lung tumors and the corresponding primary tumor. We detected fewer mutations in lymph nodes as compared to the corresponding primary tumor (Fig. 2b), likely because of reduced sequencing depth due to stromal contamination. Eight genes were mutated in 38% of primary melanomas of HPN mice that were never mutated in HP mouse tumors (Fig. 2c). The number of these genes that were mutated per mouse (mutation score) correlated directly with lung metastasis score. HP and HPN primary melanomas exhibited equal numbers of mutations consistent with COSMIC signatures of UV-induced damage, APOBEC deficiency, alkylating agent exposure, BRCA1/BRCA2 deficiency and POLH deficiency (Supplementary Table 2).

Fig. 2

Summary of nonsilent mutations identified in primary melanomas of HP and HPN mice.

Genomic sequencing data from mouse melanoma specimens was deposited in SRA (Accession number PRJNA749368). a Number of nonsilent mutations (missense, nonsense and splice site) in primary melanomas of individual HP and HPN mice. b, Number of nonsilent mutations in primary melanomas and paired metastatic tissues from HPN mice. c, Missense mutations in eight signature genes in primary melanomas of HP and HPN mice. Red boxes indicate mutated genes, with amino acid substitutions noted. Listed at the right of each row is the P value for the Chi-Square test comparison of mutation frequency between HP and HPN mice for each gene. Listed below each column are the number of genes harboring nonsilent mutations (HPN signature mutation score) and the lung metastasis scores for each mouse. The Spearman rank test was used to analyze the correlation between mutation and metastasis scores. d Missense mutations in five G-alpha protein genes in primary melanomas of HP and HPN mice, analyzed as described in panel c. e Missense mutations in HPN signature and G-alpha genes in primary melanomas and paired metastatic tissues of HPN mice. Listed at the right of each row is the concordance of mutations between paired primary and metastatic tissues for each gene.

Many HP and HPN melanomas (8 of 17) harbor activating mutations in the Q209-encoding codon of Gnaq or its Gα subunit paralogue Gna11 (1 of 17) (Fig. 2d), as recently reported in UV-induced melanomas with the HP mouse model[48]. Mutations in the Gα-encoding genes Gna14, Gna15 and Gnas were also identified in the current study. However, incidence of Gα protein mutations was not correlated with lung metastasis. Mutations in GNAQ and related Gα genes have been described for human uveal melanomas and blue nevi[68,69]. The HPN mutation signature and Gα mutations found in primary melanomas were detected in 100% of paired lung metastases and 35% of paired lymph nodes (Fig. 2e). Fifty-three other genes were mutated only in metastatic tissues and not in paired primary tumors, although most (44/53) occurred in only one mouse (Supplementary Fig. 3). Thus, these studies identify an eight-gene missense HPN mutation signature that is associated with increased lung metastasis in UV-induced melanoma.

Identification of copy number variation and gene fusion events

Copy number variations were detected by significance testing using STAC for chromosomal-level analysis[13] and GISTIC (The Genomic Identification of Significant Targets in Cancer) for focal-level analysis[8] (FDR/Q <5%; Methods). Regions of copy number gain were seen in melanomas of multiple HPN mice in chromosomes 1 and 6, as well as copy number loss in chromosomes 8, 10–13, 16 and 19 (Supplementary Fig. 4a). Copy number gain was also seen at chromosome 6 in the majority of HP melanomas, with modest gains also seen in chromosomes 11, 15 and to a lesser extent in other chromosomes. HPN melanomas displayed modest focal copy number gains in chromosomes 1, 5, 10, 12 and 14, along with losses at chromosomes 5 and 6 (Supplementary Fig. 4b). Modest copy number loss was detected in HP melanomas at chromosomes 5, 6, 12 and 13. Copy number variation in 59 chromosomal regions was significantly different between HP and HPN melanomas (Supplementary Fig. 5). Copy number variation event scores within a given melanoma were correlated with lung metastasis scores for 14% of the chromosomal regions (Supplementary Fig. 5). Copy number variation events in primary melanomas were detected in 63% of paired lung metastases and 43% of paired lymph nodes (Supplementary Fig. 6).

RNA-seq analysis detected numerous gene fusion events in primary melanomas from 7 HP and 7 HPN mice, and in 6 lymph node and 2 lung metastatic lesions from HPN mice. Eighty-four fusions were identified in primary melanomas of HP and/or HPN mice, with most occurring between gene homologs (Supplementary Table 3). Five fusion events occurred more frequently in HPN versus HP tumors (Supplementary Fig. 7a), with none occurring more frequently in HP tumor. Nine fusions occurred in at least one-third of primary melanomas, regardless of genotype. Fusion events found in primary melanomas were present in 100% of paired lung metastases and 58% of paired lymph nodes (Supplementary Fig. 7b). Six fusion events occurred only in metastatic lesions (Supplementary Fig. 7c). Frequency of individual fusion events was not correlated with lung metastasis scores.

Identification of a 32-gene expression signature associated with the HPN genotype and lung metastatic activity

RNA-seq analyses detected 18,169 transcripts in primary melanomas derived from HP and HPN mice. HPN melanomas exhibited reduced expression of NME1 (2.4-fold, P = 0.038) and NME2 (2.8-fold, P = 0.005) mRNA. 106 differentially expressed genes (DEGs; >2.5-fold change, P <0.05) exhibited significantly lower mean expression in HPN versus HP primary melanomas (Fig. 3a and Supplementary Table 4), while none were expressed at higher levels. 179 genes had significantly lower expression in HPN metastases compared to paired primary melanomas and 741 genes had higher expression (Fig. 3b, Supplementary Table 4). Principal components analysis of the 106 DEGs from primary melanomas identified 32 genes with the highest levels of covariance for PC1 (eigenvectors >0.03) (Fig. 3c and Supplementary Table 5). PC1 and PC2 scores for HPN melanomas were tightly clustered near the origin, indicating expression of the 32 genes was coordinately reduced (Fig. 3d), but these scores were not clustered in HP melanomas, reflecting a diverse pattern of gene expression. Importantly, low PC1 scores were associated with high lung metastasis scores (Fig. 3e). We refer to these 32 genes as the HPN lung metastasis signature (HPN-LMS).

Fig. 3

Identification of a 32-gene expression signature associated with lung metastasis.

RNA-seq data from mouse melanoma specimens was deposited in GEO (Accession number GSE181207). a Volcano plot based on P values (-log10) and fold-changes (log2) derived from comparisons of means for all detectable transcripts for primary melanomas of HP and HPN mice, as measured by RNA-seq. Genes expressed at a significantly lower level (p < 0.05, fold-change > 2.5) in primary melanomas of HPN versus HP mice (Supplementary Table 5) are highlighted in green. b Volcano plot as in panel a, comparing transcripts with differential expression between HPN primary melanomas and paired metastases. Genes expressed at significantly lower levels in metastatic specimens are highlighted in green; genes expressed at higher levels are highlighted in red. c Principal components analysis (PCA) of the 106 DEGs. Shown are PC1 and PC2 eigenvectors of the covariance matrix for 32 genes with PC1 eigenvectors > 0.03. This list of 32 genes was denoted the “HPN lung metastasis signature” (HPN-LMS). d Plot of PC1 and PC2 scores for each primary melanoma. Scores represent sums of products of the eigenvector and the relative expression of each of the 32 DEGs for each melanoma specimen. e P-values for Mann-Whitney rank sum tests of differences in median PC1 and PC2 scores between primary melanomas of HP and HPN mice. Also shown are Spearman rank correlation coefficients and P-values for the comparison between PC scores and lung metastasis scores across all HP and HPN mice. f Heatmap displaying relative expression of HPN-LMS genes in the indicated primary melanoma specimens from HP and HPN mice. Color coding for individual specimens denotes the fold difference from the mean expression for each gene. Asterisks denote genes whose expression was individually correlated with lung metastasis scores. HPN-LMS scores were calculated as the average relative expression of all HPN-LMS genes for each specimen. Shown below the heatmap is the P-value for a Mann-Whitney rank sum test of the difference in median LMS scores between primary melanomas from HP and HPN mice. Also shown are the Spearman rank correlation coefficient and P-value for the comparison between HPN-LMS and lung metastasis scores across all HP and HPN mice. g Plot of PC1 versus PC2 scores for primary and metastatic melanomas of HPN mice, as derived from a PCA conducted with all 179 genes expressed at lower levels in metastatic versus paired primary melanoma specimens (see Supplementary Table 4). The PCA identified 19 genes whose eigenvectors for PC1 were greater than 0.03 (see Supplementary Table 5). These eigenvectors were used to calculate PC1 and PC2 scores for each primary and metastatic melanoma specimen as described in (d). h Plot of PC1 versus PC2 scores for primary and metastatic melanomas of HPN mice, as derived from a PCA conducted with all 741 genes expressed at higher levels in metastatic versus paired primary melanoma specimens (see Supplementary Table 4). The PCA identified 95 genes whose eigenvectors for PC1 were greater than 0.03 (see Supplementary Table 5). These eigenvectors were used to calculate PC1 and PC2 scores for each primary and metastatic melanoma specimen as described in d. Summarized below the graphs in g and h are P-values for paired t-tests of differences in mean PC1 and PC2 scores between primary and metastatic melanomas of HPN mice.

A heatmap reveals reduced expression of signature genes in primary HPN tumors (Fig. 3f). Expression of 19 signature genes was negatively correlated with lung metastasis scores. A composite HPN-LMS score was calculated for each primary melanoma by averaging the relative expression of all 32 HPN-LMS genes. The median HPN-LMS score was significantly lower in HPN versus HP primary melanoma, and low composite scores were correlated with high lung metastasis scores across all mice. No pattern emerged when comparing expression of individual HPN-LMS genes between HPN primary melanomas and paired metastases (Supplementary Fig. 8). Likewise, composite HPN-LMS scores for all 32 genes were not significantly different between HPN primary melanomas and paired metastases. Thus, HPN-LMS genes comprise a signature for primary melanomas, and their expression is not altered significantly in metastatic tissues.

In contrast to the HPN-LMS, other genes that were differentially regulated in metastatic tissues exhibited distinctive patterns of covariance (Fig. 3g, 3h and Supplementary Table 5). PCA of 179 genes with lower mean expression in metastases yielded tightly clustered PC scores for lung and lymph node specimens near the origin, suggesting coordinate down-regulation of these genes. Principal component analysis of 741 genes with higher mean expression in metastases led to tight clustering of scores for metastatic lung tissues with primary melanomas, while scores for metastatic lymph node tissues were widely scattered. These findings strongly suggest that the increase in expression of these genes is due stromal contamination in lymph node specimens.

HPN-LMS genes were used to identify associations with functional pathways by Ingenuity© Pathway Analysis (IPA, QIAGEN) (Supplementary Table 6). Genes from two IPA pathways, “Migration of Cells” and “Recruitment of Macrophages” were collectively expressed at lower levels in HPN mouse primary melanomas and negatively correlated with lung metastasis scores (Fig. 4a, 4b). Interestingly, expression of mRNAs encoding other important regulators of melanoma phenotype[3] including Mitf, Axl, Wnt5a/5b and Egfr were not different in HP and HPN melanomas, nor was expression of 15 housekeeping genes (Fig. 4c). In conclusion, these analyses describe a 32-gene signature that is strongly associated with lung metastasis in UV-induced melanoma.

Fig. 4

Expression of genes associated with four melanoma metastasis-relevant pathways is negatively correlated with metastatic activity of primary melanomas of HP and HPN mice.

Heatmaps display relative expression levels of genes associated with two IPA pathways: a “Migration of Cells” and b “Recruitment of Macrophages”; subsets of genes with an HP/HPN ratio of 2-fold or more were included in the heatmap. c Relative expression levels of a panel of housekeeping genes. Individual specimens are color-coded based on the fold-difference in gene expression relative to the mean of all specimens for that gene. The expression score for a given mouse is the mean of the relative expression for all genes in that pathway. Lung metastasis scores for each mouse are listed below the average heat scores. Mann-Whitney rank sum tests were used to compare expression scores between HP and HPN primary melanomas. Spearman rank correlation tests were used to compare expression scores and lung metastasis scores across all HP and HPN mice.

Gene mutations predict mutational targets in human melanoma

Mutation frequencies of HPN signature genes and GNAQ (Fig. 2) were determined in 12 human databases classifying cutaneous, acral, desmoplastic and metastatic melanoma (Fig. 5 and Supplementary Fig. 9). All eight HPN signature genes were mutated at high rates in one or more databases, with NCKAP5, KIF21B, ATP8B4, PCDHA2 and BRCA1 displaying the highest frequencies (15–30%). Interestingly, NCKAP5, KIF21B, BRCA1, ARHGAP35 and IFT172 were mutated more frequently in databases enriched for metastatic melanoma[25, 67]. KIF21B, ATP8B4, PCDHA2 and BRCA1 mutations were present in desmoplastic melanomas, while signature genes were rarely mutated in uveal melanomas. No CNVs or gene fusions occurring at high frequency in HPN melanomas (Supplementary Figs. 4–7) were detected in the human melanoma databases.

Fig. 5

Genes enriched for missense mutations in metastatic melanomas of HP and HPN melanomas are frequently mutated in human melanomas.

Displayed is the incidence of mutations in genes comprising the HPN mutation signature and GNAQ in human melanoma. Frequency of point mutations, gene amplifications, deep deletions and other alterations is displayed in the bar graphs, with color-coding for the different types of alterations indicated at the top of the figure. Graphs were generated using the Cancer Types Summary function in CBioPortal (https://www.cbioportal.org/). Human melanoma databases queried for the analysis included the following, as identified in each panel: TCGA/SKCM (Firehose Legacy/TCGA/Skin Cutaneous Melanoma, 479 samples), SKCM-Yale[31] (147 samples), SKCM-Broad1[22] (121 samples), SKCM-Broad2[66] (78 samples), Melanoma-MSK[59] (4 samples), Melanoma-Broad[7] (Broad/DFCI, 26 samples), Metastatic-UCLA[25] (38 samples), Metastatic-DFCI[67] (110 samples), Acral[35] (38 samples), Desmoplastic[55] (20 samples), Uveal-TCGA (Firehose Legacy, 80 samples) and Uveal[26] (QIMR, 28 samples).

Expression of HPN-LMS genes predicts survival in cutaneous and uveal melanoma patients

To determine whether the HPN-LMS identified in mice reflected profiles of gene expression in human melanoma patients, we analyzed expression of HPN-LMS genes in the skin cutaneous melanoma cohort of the Cancer Genome Atlas (TCGA-SKCM). Five patient clusters were identified with distinct profiles of gene expression (Fig. 6a). Overall survival was significantly different among the five clusters (Fig. 6b), with cluster 1 exhibiting the longest median survival time (168m), followed by clusters 2 (96m), 3 (66m), 4 (58m), and 5 (49m). HP and HPN melanomas displayed a high rate of activating mutations in Gα subunit genes, a hallmark of uveal melanoma in humans[68]. Accordingly, we evaluated utility of the HPN-LMS in predicting survival in a uveal melanoma cohort of 80 patients (TCGA-UVM, Firehose Legacy). Five patient clusters were identified that exhibited distinct profiles of HPN-LMS gene expression (Fig. 7a) and large differences in overall survival (Fig. 7b). Mean overall survival varied significantly among patients with different primary diagnoses of UVM: spindle (54 months), epithelioid (25 months) and mixed spindle/epithelioid (39 months) (Fig. 7c). A highly significant interaction was observed between HPN-LMS gene expression and primary diagnosis (Fig. 7d). Spindle UVM was more frequent in HPN-LMS clusters with higher survival, whereas epithelioid UVM was more frequent in HPN-LMS clusters with lower survival (P = 8.55 × 10−36).

Fig. 6

The HPN-LMS predicts survival in cutaneous melanoma patients (TCGA-SKCM).

a Heatmap of HPN-LMS gene expression in TCGA-SKCM subjects. Access to the dataset “Disease-Specific (Melanoma) in Melanoma Genome Sequencing Project” was approved by the National Institues of Health (project #24178, David Kaetzel, P.I.). Human RNAseq data in raw count format for the TCGA-SKCM and TCGA-UVM cohorts were transformed into log2CPM values using R package “edgeR”, then normalized for each gene by dividing the value for each individual sample by the mean for all samples. ADH7 and PLCXD2 were excluded from this heatmap because mean RNA expression levels of these genes were below detection. Hierarchical cluster analyses were performed using ComplexHeatmap (v1.99.7). Small clusters from the dendrogram were combined into five larger clusters that exhibited differences in mean survival time, then numbered from best to worst survival. b Survival log rank analysis of the five clusters of TCGA-SKCM subjects using metadata from the TCGA-SKCM database.

Fig. 7

The HPN-LMS predicts survival in three subtypes of uveal melanoma patients (TCGA-UVM).

a Heatmap of HPN-LMS gene expression in TCGA-UVM subjects. Access to the dataset “Disease-Specific (Uveal Melanoma, NPU) in Genome-wide Discovery of Novel Cancer Predisposing Mutations in Melanoma” was approved by the National Institues of Health (project #24178, David Kaetzel, P.I.). RNA-seq counts were transformed as log2 and normalized to the mean of all subjects for each gene. ADH7, FLRT3 and PLCXD2 were excluded from this heatmap because mean RNA expression levels of these genes were below detection. Hierarchical cluster analyses were performed as described in Figure 5. b Survival log rank analysis of TCGA-UVM subjects, grouped in five hierarchical clusters as indicated in panel a, using metadata from the TCGA-UVM database. Chi-square analysis demonstrated that the clusters had significantly different survival times (P = 0.00003). c Survival log rank analysis of TCGA-UVM subjects grouped by primary diagnosis. d Grouping of TCGA-UVM subjects from each primary diagnosis category into HPN-LMS clusters (see panels a and b). Chi-square analysis was used to compare the distribution of subjects within each primary diagnosis category across the five HPN-LMS RNA clusters.

Patterns of gene expression within HPN-LMS clusters in the SKCM and UVM cohorts revealed associations between individual genes and patient survival (Fig. 8a). Higher survival was associated with increased expression of NYNRIN and RND3 in both SKCM and UVM subjects, with FLRT3, GALNT6, S100A4 and TFPI in SKCM subjects, and with ARRDC3, MAOA, MXRA8, PDE3B and TMOD2 in UVM subjects (blue circles in Venn diagram). By contrast, lower survival was associated with increased expression of ELOVL7, SHTN1 and SLC22A23 in both SKCM and UVM subjects, and with CASP1, CLRX, ITGA2, P2RY1, S100A4 in UVM subjects (red circles). Other genes, while part of the multifactorial pattern of expression in the HPN-LMS clusters, were not individually associated with survival in the SKCM and UVM cohorts (yellow circle).

Fig. 8

Four HPN-LMS genes (ARRDC3, FLRT3, NYNRIN and RND3) associated with higher overall survival in human melanoma patients exhibit metastasis-suppressor functions in human melanoma cells.

a Venn diagram identifying individual genes associated with survival in TCGA-SKCM and TCGA-UVM subjects. Genes correlated with higher survival (blue circles) had median expression in the top quartile (75% or higher) in subjects within Cluster 1 of TCGA-SKCM (see Fig. 6) or Clusters 1 and 2 of TCGA-UVM (see Fig. 7) compared to median expression in subjects within other clusters. Genes correlated with lower survival (red circles) had median expression in the top quartile in subjects within Cluster 5 of TCGA-SKCM or TCGA-UVM compared to median expression in subjects within other clusters. Expression of genes in the central yellow circle was not significantly correlated with survival. b siRNA-mediated knockdown of gene expression in WM9 human metastatic melanoma cells (gift from M. Herlyn, Wistar Inst.; authenticated by STR profiling prior to experiments) robustly suppresses expression of four HPN-LMS genes. WM9 cells were subjected to two cycles of electroporation with the indicated siRNAs as described[17], total cellular RNA was collected 24h later, then levels of the indicated RNAs were quantified by qRT-PCR. Data are expressed as mean + SEM (n = 3). Asterisks indicate that the RNA level was significantly lower in cells treated with targeted siRNA compared to control siRNA (t-test; P < 0.001). The level of each protein was detected by immunoblot, with β-actin used as a loading control. Sources of antibodies: anti-ARRDC3 (ab134289, Abcam, Cambridge, UK); anti-FLRT3 (NBP2–36566, Novus Biological, Littleton, CO, USA); anti-NYNRIN (orb475599, Biorbyt, St. Louis, MO, USA); anti-RND3 (ab171799, Abcam); anti-β-actin (A5441, Sigma-Aldrich, St. Louis, MO, USA); ECL-anti-mouse IgG (NA931, GE Healthcare, Little Chalfort, UK); ECL-anti-rabbit IgG (NA934, GE Healthcare). c Knockdown of three HPN-LMS genes enhances invasion activity of WM9 human metastatic melanoma cells. WM9 cells were treated with the indicated siRNAs as in panel b, then 24h later cells invading through a Matrigel matrix were counted as described[17, 28]. On the left of the bar graph are representative images of DAPI-stained nuclei. Data in the bar graph were pooled from two independent experiments and expressed as mean + SEM (n = 6). Significant differences in invasion between cells treated with gene-specific vs. control siRNA were determined by 2-way ANOVA and Kruskal-Wallis post-hoc tests. P-Values are indicated above the bar for each siRNA. d Knockdown of two HPN-LMS genes enhances proliferative activity of WM9 human metastatic melanoma cells. WM9 cells were treated with the indicated siRNAs as in panel b, then 24h and 48h later cells were counted as described[17, 28]. To the left of the bar graph are representative images of cell cultures 24h after plating of siRNA-treated cells. The bars indicate 50 microns. Data in the bar graph are expressed as mean + SEM (n = 3). Significant differences in invasion between cells treated with gene-specific vs. control siRNA were determined by 2-way ANOVA and Kruskal-Wallis post-hoc tests. Asterisks indicate that proliferation was significantly higher in cells treated with the indicated siRNA vs. control siRNA after 24h (green) or 48 h (red) (P < 0.001). The red hashmark indicates that proliferation was significantly lower in cells treated with the indicated siRNA vs. control siRNA after 48 h (P < 0.001).

To obtain biological evidence that HPN-LMS genes regulate cancer-relevant endpoints in melanoma cancer cells, we examined the impact of ARRDC3, FLRT3, NYNRIN and RND3 knockdown on invasion and proliferation of WM9 human metastatic melanoma cells. These genes were selected based on the strength of their association with improved survival (Fig. 8a) and their relative novelty with respect to human melanoma. siRNA-dependent knockdown was highly efficient, with robust decreases in steady-state RNA expression (> 95%) and protein expression achieved for all four target genes (Fig. 8b). Knockdown of ARRDC3, RND3 and NYNRIN resulted in highly significant increases in invasion of WM9 cells through Matrigel, suggesting that these gene products suppress invasion (Fig. 8c), while FLRT3 knockdown had no effect. The effects of ARRDC3 and RND3 knockdown on invasion were replicated in a second human melanoma cell line, WM1158 (Supplementary Fig. 10). In addition, knockdown of FLRT3 and RND3 elicited modest increases in cell proliferation rates, suggesting proliferation-suppressing activities of these genes, while ARRDC3 knockdown resulted in a small decrease in proliferation after 48h of growth (Fig. 8d). NYNRIN knockdown had no effect on cell proliferation. The small increase in proliferation elicited by RND3 knockdown at 24h (28%) was insufficient to explain the dramatic increase in invasion activity (88%) observed during that time period. Taken together, these studies suggest ARRDC3, RND3 and NYNRIN are key effectors of the metastasis suppressor functions of NME1 and NME2, and support the potential value of HPN-LMS genes as therapeutic targets and markers of survival in human melanoma patients.

Discussion

NME1 often exhibits reduced expression in melanoma or other cancers but the gene is rarely mutated or deleted, strongly suggesting that transcriptional and/or post-transcriptional mechanisms are at play[10]. Moreover, we previously identified rare cell subpopulations in human melanoma cell lines that exhibit coordinate downregulation of NME1 and NME2 expression, and these cells exhibit elevated metastatic activity[60]. Our HPN melanoma model phenocopies these findings and reveals genomic and transcriptomic changes resulting from Nme1/Nme2 deficiency that correlate with enhanced lung metastasis. HPN tumors harbor non-silent mutations in eight genes (Arhgap35, Atp8b4, Brca1, Ift172, Kif21b, Nckap5, Pcdha2 and Zfp869) and these mutations are associated with enhanced lung metastasis. Screening of The Cancer Genome Atlas shows that these genes were also missense-mutated in human melanoma suggesting a role in disease progression. Brca1 mutations, that compromise double-strand break repair, are associated with enhanced tumor initiation, progression and metastasis[2,49]. Our studies show that the Brca1 mutations in HPN mice are associated with increased lung metastasis and additional analysis shows that these mutations occur at a high rate in melanoma patients with metastatic disease (Fig. 5)[25, 67]. This suggests that impaired Brca1 results in insufficient DNA repair at the later stages of melanoma progression. ARHGAP35 is a RhoA-specific GTPase implicated as a suppressor of the cancer phenotype in melanoma[42] and colon cancers[11]. Arhgap35 mutations in HPN tumors were located within the RhoGAP and pseudo-GTPase domains of the protein, suggesting they disrupt its GTPase and metastasis-suppressing functions. In addition, both HP and HPN melanomas exhibited a high frequency of coding mutations (e.g. Q209) in Gnaq and other Gα homologs (i.e. Gna11, Gna15, Gnas). The Q209 mutation disrupts the intrinsic ATPase function of human GNAQ, resulting in constitutive G protein and oncogenic activity[40, 69]. Interestingly, these mutations occur at very high frequency (> 30–50%) in human uveal melanomas and blue nevi, but at lower rates in subcutaneous and other melanoma subtypes (Fig. 5). Through direct interactions with DNA, NME1 plays important roles in maintenance of genomic stability in human melanoma cell lines, with evidence strongly suggesting participation by its NDPK and 3’−5’ exonuclease activities[27]. A potential mechanism is suggested by participation of NME1 and its cognate enzymatic activities in the transcription-coupled repair mode of NER, which is restricted to actively transcribed genes[62].

Thirty-two genes whose expression was significantly lower in HPN melanomas and was inversely correlated with lung metastasis scores were collectively termed the HPN-LMS. Cancer-relevant activities have been ascribed to most of these genes including Arhgef5[21], Casp1[41], Itga2[19], Notch1[5, 74], S100A4[16], Rnd3/RhoE[6, 54], Slc22a23[12, 15] and ElovL7[64]. The pattern of expression of HPN-LMS genes in the SKCM and UVM cohorts of the TCGA predicted reduced overall survival. Interestingly, higher expression of NYNRIN and RND3 was associated with longer survival in SKCM and UVM patients, while higher expression of ELOVL7, SHTN1 and SLC22A23 was associated with shorter survival in these patients. The similarity in expression patterns of these genes between SKCM and UVM patients strongly suggests these genes regulate metastatic activity and survival.

We previously reported gene signatures associated with metastatic activity of human melanoma cell lines with variable expression of NME1 and/or NME2[32,60]. Genes induced by forced NME1 expression[32] did not overlap with the LMS identified herein, which could be secondary to the established melanoma cell line used for the analyses, the application of forced NME1 expression, and/or the selective modulation of NME1 but NME2 expression. The LMS was as efficacious a predictor of survival in patients of the TCGA-SKCM cohort as the signature elicited by forced NME1 expression. Interestingly, reduced expression of the LMS gene Itga was also a component of the signature associated with elevated metastatic activity in a rare subpopulation of melanoma cells with coordinately reduced expression of NME1 and NME2[60].

We identified ARRDC3, NYNRIN and RND3 (Fig. 8) as potential effectors of the metastasis suppressor function of the NME1/NME2 locus. ARRDC3 (α-arrestin domain-containing protein-3) is implicated in ubiquitination[4] but has yet to be studied in the context of melanoma. It exhibits antiproliferative and anti-invasive activities in cultured triple-negative breast cancer (TNBC) cells[14, 61] but not in human TNBC patients. Tumor suppressor activities have been attributed to ARRDC3 in carcinomas of the prostate[76] and kidney[70] but metastasis suppressor function was not addressed in either study. Our study points to ARRDC3 as being an important metastasis suppressor in human melanoma. Interestingly, increased ARRDC3 expression is associated with improved survival in UVM patients but not in SKCM patients. NYNRIN (NYN domain and retroviral integrase-containing) is an RNA processing factor which has received only limited attention in the context of human cancer. NYNRIN is part of a 10-gene expression signature associated with poor survival in acute myelocytic leukemia (AML)[73], suggesting a tumor-promoting function. In contrast, NYNRIN-inactivating mutations are enriched in Wilms tumor[37], suggesting a tumor suppressor activity. However, neither study directly tested the impact of NYNRIN on tumor cell phenotype, and its putative AML-promoting activity is opposite from the metastasis suppressing function we have demonstrated in human melanoma cells.

RND3 (Rho family GTPase 3) is an atypical Rho GTP exchange factor that exhibits both tumor-promoting and suppressor functions, depending on cancer type[47]. Silencing of RND3 expression inhibits invasive activity in some melanoma cell lines (e.g. WM793 and WM115)[30], consistent with an invasion-driving function. In contrast, our findings show reduced RND3 expression correlates with increased metastatic activity in HPN mice, an association with increased survival in SKCM and UVM patients, and an ability to suppress invasion of WM9 cells. This difference is likely due to context-dependent effects of RND3 on metastasis-relevant activities.

In conclusion, our novel findings reinforce the concept that metastasis is controlled by gene networks rather than by altered activity of one or a few genes. In addition, our findings suggest that the mutational and transcriptomic signatures we identified may provide novel markers for disease prognosis as well as novel therapeutic targets. Our data indicate a necessity to target multiple genes within the network to achieve maximum prognostic and/or therapeutic potential. As larger and better-annotated cohorts of melanoma patients become available, stronger associations of the HPN signature with disease progression and patient survival should be possible. We propose that the HPN mouse model holds considerable promise as a preclinical platform for identification of novel biomarkers and therapeutic targets for clinical management of melanoma patients. This model is particularly attractive for uveal melanoma, by virtue of its high incidence of activating mutations in Gα subunit genes and the efficacy of HPN-LMS gene expression profiles for predicting survival in uveal melanoma patients.

Materials and methods

Generation and surveillance of melanoma in mice

Generation of mice is summarized in Supplementary Fig. 1a. The HP and Nme1/2+/Δ mouse strains, including methodologies for genotype screening, have been described[50,51]. The Nme1+/− / Nme2+/− strain harbors a deletion that eliminates all coding sequence within the respective genes[50]. Because mice of different genotypes were generated randomly during the breeding process, no further randomization was required for the HP and HPN experimental groups. Melanoma was initiated in a blinded fashion at postnatal day four with a single erythematous dose of UV radiation[45] (Fig. 1C). Incidence and dimensions of pigmented skin lesions were monitored weekly[65]. Protocols for murine experiments were compliant with ethical regulations and approved by the Institutional Care and Use Committee at the University of Maryland-Baltimore (Protocols 0612013 and 0515008; D. Kaetzel, P.I.).

Sequencing of whole genomes and transcriptomes, and analyses of copy number, differential gene expression and functional pathways

The presence of melanoma cells (>80%) was verified by histopathological staining in all primary melanoma specimens subjected to whole genome sequencing and RNA-seq analyses. Genomic DNA and RNA were prepared using the AllPrep Mini Kit (QIAGEN, Germantown, MD). Paired tumor/normal whole tumor genome libraries were prepared using the KAPA HyperPrep Kit (Roche, Basel, Switzerland). RNA sequencing libraries were generated using the NEBNext® Ultra™ II Directional RNA Library Prep Kit (E7765, NEB, Ipswich, MA). 150bp paired-end sequencing was performed using the NovaSeq 6000 v1 system (Illumina, San Diego, CA). FASTQ files were converted using BCL2FASTQ v1.8.4. Whole genome sequencing data were aligned using mm9_build37.fa reference, and variant calling was performed using GATK best practices and validation with internal pipelines as previously described[1, 38]. Copy number analysis was performed utilizing Nexus Copy Number v10 (BioDiscovery, El Segundo, CA). RNAseq data was aligned using STAR v2.5.3a, DESeq2 was utilized for differential expression analysis, and Salmon v0.7.2 package to quantify isoforms[36, 43, 46, 52, 53]. TopHat Fusion[29], FusionCatcher (DOI: https://10.1101/011650), and STAR-Fusion (DOI: 10.1101/120295) were used to detect potential oncogenic fusion events. Only fusions that were present in at least 2 of the callers were selected. The R v3.4.0 package was utilized for statistical and graphical processes.

Statistical methods

Sample size for mouse experiments was based on the estimate that 10 melanoma-positive mice/group are sufficient to detect a 20% difference in metastasis incidence (p < 0.05, power of 0.8)[45]. Parametric analyses (Student’s t-test, ANOVA, Pearson correlation) were used to compare normally distributed variables, while nonparametric analyses (Mann-Whitney rank sum test, ANOVA on ranks, Spearman rank correlation) were used to compare non-normally distributed variables. Paired t-tests were used for comparisons between variables from HPN primary tumors and paired metastatic tissues. Statistical tests for equal or non-equal variance were employed as appropriate for each comparison. Simple linear regression was used to calculate rates of growth of primary melanomas. Chi-square analysis was used to compare frequencies of categorical variables between groups. To compare human subjects in TCGA-SKCM and TCGA-UVM cohorts, Kaplan-Meier log rank survival plots were used to assess differences in overall survival based on a single categorical variable, i.e., subjects grouped by patterns of RNA expression or a single metadata annotation. The Cox regression stratified model was used to assess the combined effects on overall survival of a single categorical variable and one or more covariates. For invasion and proliferation assays, number of cell culture replicates required to achieve statistical significance (P < 0.05, power of 0.8) was based on prior observations of variability in those systems[17].