Ligand-activated BMP signaling inhibits cell differentiation and death to promote melanoma

The identification of new cancer-promoting genes has yielded mechanistic insights into tumor progression and led to the development of several targeted therapies. In cutaneous melanoma, the finding of common BRAF mutations highlighted the importance of ERK pathway activation in tumor initiation and maintenance. These studies also triggered the design of vemurafenib and other MAPK pathway inhibitors, which were the first drugs to extend the survival of patients with advanced disease (1–4). Immunotherapies, such as the CTLA4 inhibitor ipilimumab and the PD1 inhibitors nivolumab and pembrolizumab (5, 6), have also had a major impact on melanoma therapy, as they have dramatically improved the long-term survival rates of patients with advanced-stage disease (7, 8). Despite this progress, many patients do not respond to certain therapies, whereas others suffer from drug toxicity, therapy resistance, or disease relapse (9–11), underscoring the need to identify additional targets for therapeutic intervention.

Along with identifying BRAF and other recurrently mutated cancer genes, sequencing strategies have also defined genomic intervals subject to recurrent copy number variations (CNVs). However, cancer-promoting genes in CNVs have been difficult to identify, because (a) they are often not affected by point mutations, and (b) they are typically present in large CNVs alongside several bystander genes that have no effect on tumor progression. Analysis of CNVs has the potential to uncover several new cancer-promoting genes in solid tumors such as melanomas, in which a large percentage of the genome is subject to recurrent CNV (12).

Oncogenomic studies have also revealed expression profiles that reflect broad changes in cellular identity that distinguish cancer cells from normal tissue (13). In many cancers, tumor cells adopt cellular and molecular identities that overlap with their lineally related embryonic cells. Adopting these identities can endow tumor cells with properties, such as the ability to proliferate or migrate, not found in their differentiated counterparts (14–16). Reawakening of neural crest character in nascent melanomas, as exemplified by the expression of the neural crest specification factors SNAI2 (SLUG) and BRN3A (POU4F1) has been shown to enable protumorigenic properties like migration and survival, respectively (17, 18). Subsequent studies have noted additional gene expression and functional relationships between melanoma and neural crest cells (19, 20). Whereas similarities between tumor and embryonic cells in melanomas and other cancers have been recognized, the factors that establish and maintain an embryonic identity in tumor cells are poorly understood. Specifically, it is not clear whether embryonic genes are regulated separately to reconstitute an embryonic identity or whether a programmatic change that simultaneously regulates many genes is involved. In addition, the consequences of abrogating embryonic identity in melanoma and other cancers have not been thoroughly investigated.

In this study, we report the identification of the growth differentiation factor 6 (GDF6) oncogene, a BMP factor that is recurrently amplified and specifically expressed in melanomas. GDF6, which is expressed in the embryonic neural crest, regulates a trunk neural crest gene signature in melanomas. Loss of GDF6 results in the differentiation and death of melanoma cells, indicating that GDF6 and the BMP signaling pathway are required for tumor maintenance and are thus potentially important targets in melanoma therapy.

Comparative oncogenomics and expression analyses identify GDF6 dysregulation in melanoma. We hypothesized that a cross-species comparative approach with zebrafish would aid in the identification of cancer genes in regions of broad CNVs. Humans and zebrafish are diverged by 420 million years (21), and the genomic reorganization that has occurred over time has been predicted to frequently place orthologous driver genes next to different neighboring genes in each species. Consequently, orthologous driver genes would be altered in both species, but changes to neighboring passenger genes would be limited to a single species. To test this hypothesis, we sought to compare genes that are recurrently amplified in human melanomas, roughly 10% of the genome (22), to genes recurrently amplified in zebrafish melanomas. Using melanomas that arose autochthonously in a Tg(mitfa:BRAFV600E);p53(lf) zebrafish strain (23), we performed array comparative genomic hybridization (aCGH) to generate CNV profiles (Figure 1A and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI92513DS1). The Tg(mitfa:BRAFV600E);p53(lf) strain combines a transgene that drives expression of human BRAFV600E in the melanocyte lineage with a p53 loss-of-function mutation, and every animal of this genotype ultimately develops 1 or more melanomas. aCGH values were analyzed with the JISTIC algorithm (24) to define recurrently varied intervals, which largely overlapped with intervals obtained in an independent study of zebrafish melanomas (25). Recurrently amplified genes from JISTIC intervals were compared with their human orthologs. The degree of overlap between orthologs amplified in both species (374 genes) is greater than would be expected by chance (247 genes) (Figure 1B and Supplemental Table 2), suggesting that amplification of similar driver genes mechanistically underlies tumor formation in both species. As a further indication of mechanistic conservation, known melanoma drivers, including TERT, MYC, and SETDB1, were recurrently amplified in both species,(26–28).

Figure 1

GDF6 is recurrently amplified and specifically expressed in melanomas. (A) Circos plot displaying gene copy number gains and losses of zebrafish melanomas across 25 chromosomes. JISTIC G-scores are displayed as pale red shading (amplifications [minimum = 0; maximum = 1,550]) and blue shading (deletions [minimum = 0; maximum = 2,150]). –log₁₀-transformed JISTIC Q-values with a cutoff of 0.6 (corresponding to an untransformed Q-value of 0.25) are shown as bold red lines (amplifications [minimum = 0; maximum = 11]) and bold blue (deletions [minimum = 0; maximum = 11]). Dotted circles represent the –log₁₀-transformed Q-value of 0 (center) and 11 (outer: amplification; inner: deletion). (B) Venn diagram of orthologous genes significantly amplified in human and zebrafish melanomas from a total of 10,380 human-zebrafish gene pairs (hypergeometric test, P value: 2.0 × 10^–15). (C) Genes significantly upregulated in zebrafish melanomas as compared with melanocytes (microarray data set) are plotted in order of their fold change. Only genes with a fold change of greater than 2 and an adjusted P value of less than 0.05 are plotted. Recurrently amplified genes with amplified human orthologs are indicated in red. gdf6b (large red dot) and gdf6a (large black dot) are indicated. Dashed horizontal line represents a fold change of 2. (D) Immunostaining of Tg(mitfa:BRAFV600E);p53(lf) zebrafish scales bearing melanoma cells or normal melanocytes. DAPI (blue), Gdf6b (green), Mitfa (red), and a merged image of all channels are shown. Mitfa antibody specificity is shown in Supplemental Figure 2B. Scale bars: 10 μm.

Expression analyses were used to further winnow the list of candidates. Since copy number–amplified driver genes predominantly act by upregulation of WT transcripts, we obtained transcriptional profiles of zebrafish melanomas and normal melanocytes. Briefly, unpigmented EGFP-positive melanocytes and melanoma cells were sorted from the scales and tumors, respectively, of Tg(mitfa:BRAFV600E);p53(lf);Tg(mitfa:EGFP;alb(fl) fish. RNA was prepared from each cell population and subjected to both microarray analysis and massively parallel RNA sequencing (RNA-seq). Genes recurrently amplified and transcriptionally upregulated in both species (120 genes; fold change >2, adjusted P value < 0.05) included the BMP factor GDF6 (also known as BMP13) and its zebrafish ortholog gdf6b (Figure 1, B and C, Supplemental Figure 1, A–C and Supplemental Table 3). To determine whether Gdf6b protein was similarly enriched in melanomas, we generated an antibody that specifically recognizes Gdf6b (Supplemental Figure 2A). Whereas Gdf6b protein was readily detected in tumor cells from Tg(mitfa:BRAFV600E);p53(lf) fish, we did not detect Gdf6b in normal melanocytes from these same animals (Figure 1D). Human genes often have 2 zebrafish orthologs because of a partial genome duplication in the teleost lineage. The second zebrafish ortholog of human GDF6, gdf6a, was not recurrently amplified but was among the most transcriptionally upregulated genes in melanomas compared with that detected in melanocytes (Figure 1C and Supplemental Figure 1B). Together, these data highlight the recurrent amplification and tumor-specific expression of GDF6 genes in human and zebrafish melanomas.

GDF6 orthologs were particularly interesting, because their expression pattern in zebrafish embryos suggested that they may regulate melanocyte development. Zebrafish orthologs of GDF6 are expressed during neurulation and development of the neural crest, the embryonic tissue that gives rise to melanocytes (29–31). Using ISH, we confirmed the expression of gdf6a and gdf6b in the neural tube and neural crest, respectively. Later in development, we found that the expression of these genes was absent from differentiating embryonic melanocytes (Supplemental Figure 3). Factors involved in neurulation and neural crest signaling have been previously implicated in promoting melanoma progression (15, 17, 32), and we were intrigued by the notion that a developmental role for GDF6 genes could be reiterated to promote melanomagenesis.

GDF6 modulation alters melanoma onset in zebrafish. To assess whether GDF6 genes could promote melanoma, we first examined how their elevated expression affected tumor onset in zebrafish. We expressed gdf6a and gdf6b in the melanocyte lineage of Tg(mitfa:BRAFV600E);p53(lf);mitfa(lf) zebrafish using the miniCoopR system, as previously described (28). The mitfa(lf) mutation prevents melanocyte development and melanoma formation in Tg(mitfa:BRAFV600E);p53(lf) fish. When single-cell embryos from the Tg(mitfa:BRAFV600E);p53(lf);mitfa(lf) strain were injected with miniCoopR-EGFP, 21% of these animals later developed chimeric rescue of melanocytes. However, in embryos injected with miniCoopR-gdf6a or miniCoopR-gdf6b, only,0.6% and 7% of injected animals had melanocyte rescue, respectively (Supplemental Figure 4A). We did not observe this lack of melanocyte rescue in control embryos that were injected with miniCoopR vectors that had premature stop codons engineered into the gdf6 genes (Supplemental Figure 4A). In addition to the low percentages of miniCoopR-gdf6a– or miniCoopR-gdf6b–injected embryos that showed melanocyte rescue, the embryos that were rescued had significantly lower numbers of melanocytes (Supplemental Figure 4B). Because of this poor rescue, we were only able to perform melanomagenesis assays with miniCoopR-gdf6a–injected animals. When allowed to develop to adulthood, the fish with melanocyte-driven gdf6b expression had accelerated melanoma onset (median onset = 13 weeks) as compared with EGFP-expressing controls (median onset = 17 weeks) (Figure 2A). Accelerated onset was dependent on BRAFV600E and loss of p53, as expression of gdf6b in the Tg(mitfa:BRAFV600E) transgene (n = 33) or p53(lf) (n = 24) backgrounds alone did not produce tumors.

Figure 2

GDF6 modulation alters melanoma growth. (A) Melanoma-free survival curves for Tg(mitfa:BRAFV600E);p53(lf);mitfa(lf) zebrafish injected with miniCoopR-gdf6b or miniCoopR-EGFP. Statistical analysis was performed with a Wilcoxon rank-sum test. (B) Melanoma-free survival curves for Tg(mitfa:BRAFV600E);p53(lf) and Tg(mitfa:BRAFV600E);p53(lf);gdf6a(lf) zebrafish. Statistical analysis was performed with a Wilcoxon rank-sum test. (C) Immunoblot showing expression and quantification of GDF6 protein levels (relative to GDF6 protein in A375 melanoma cells) in melanoma cell lines. GAPDH was used as a loading control. Copy number values of the GDF6 locus in the different melanoma cell lines obtained from the COSMIC database (http://cancer.sanger.ac.uk/cosmic) are shown. (D) Immunoblots showing expression of GDF6 and GAPDH in A375 melanoma cells and M14 melanoma cells overexpressing GDF6. (E) Tumor formation in mice injected with A375 cells and M14 cells overexpressing GDF6 or empty vector control (1 × 10⁶ cells were injected per mouse). Error bars represent the mean ± SEM. n = 3. (F) Immunoblots showing expression of GDF6 in A375 cells (top) expressing an shRNA targeting EGFP or 2 independent GDF6-targeted shRNAs and M14 cells (bottom) expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. (G) Colony formation assay with A375 cells (left) expressing an shRNA targeting EGFP or 2 independent GDF6-targeted shRNAs and M14 cells (right) expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. Error bars represent the mean ± SEM. n = 3. (H) Tumor formation in mice injected with A375 cells expressing an shRNA targeting EGFP or 2 independent GDF6-targeted shRNAs and M14 cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1 (1 × 10⁷ cells were injected per mouse). Error bars represent the mean ± SEM. n = 3. *P < 0.05, **P < 0.01, and ***P < 0.001, by 2-tailed Student’s t test (E and G [right] and H [right]) or 1-way ANOVA with Dunnett’s test (G [left] and H [left]).

We next assayed the consequences of GDF6 loss in vivo using a zebrafish melanoma model. Since gdf6a loss-of-function animals were previously established (33), we used these mutants to test the role of GDF6 loss in melanoma onset. We bred a gdf6a loss-of-function, mutation into tumor-prone Tg(mitfa:BRAFV600E);p53(lf) zebrafish and found that the resulting Tg(mitfa:BRAFV600E);p53(lf);gdf6a(lf)–mutant animals had substantially delayed melanoma onset as compared with their control siblings (Figure 2B). Together, these results indicate that GDF6 ortholog activity is important for melanoma initiation, as overexpression of gdf6b caused accelerated tumor onset, whereas loss of gdf6a caused a delay in disease onset.

We also found that gdf6a-mutant zebrafish had a profoundly altered pigmentation pattern. We observed an increase in melanocytes in the flanks of gdf6a-mutant homozygotes, whereas WT siblings had normal pigmentation (Supplemental Figure 4C). Given these observations, in combination with our finding that increased expression of GDF6 orthologs caused a reduction in the number of embryonic melanocytes, we speculate that these genes have a role in inhibiting melanocyte development.

GDF6 modulation affects the tumorigenicity of human melanoma cells. We next wanted to test whether GDF6 modulation altered the growth of human melanoma cells. We first identified cell lines with GDF6 amplification (A375, SK-MEL-28, and SK-MEL-5) and others without amplification (M14 and C32) (34). Those with amplification had higher levels of GDF6 protein as compared with levels in the nonamplified lines (Figure 2C; see complete unedited blots in the supplemental material). We overexpressed GDF6 in A375, SK-MEL-28, and M14 melanoma cells (Figure 2D and Supplemental Figure 5A), followed by their transplantation into immunocompromised mice. In each case, elevation of GDF6 expression increased the tumor-forming potential of these melanoma cells as compared with the empty vector controls (Figure 2E and Supplemental Figure 5, B and C). GDF6 overexpression also caused an increase in the colony-forming potential of A375 cells, indicating that it can also enhance the tumorigenic capacity in vitro (Supplemental Figure 5D). These data indicate that increased GDF6, regardless of its endogenous levels, can promote the tumorigenicity of melanoma cells.

To determine the effects of GDF6 loss, we knocked down endogenous GDF6 in melanoma cells with amplifications or normal copy number of the GDF6 locus. Knockdown in amplified A375, SK-MEL-28, and SK-MEL-5 cells led to a growth disadvantage in vitro, as measured by anchorage-dependent colony formation,(Figure 2, F and G, and Supplemental Figure 5, E and F). In xenografts of A375 and SK-MEL-28 cells with GDF6 knockdown, we observed a substantial reduction in the tumorigenic potential compared with that of control cells with EGFP knockdown (Figure 2H and Supplemental Figure 5G). However, knockdown of GDF6 in nonamplified M14 and C32 cells caused little change in anchorage-dependent colony formation,(Figure 2, F and G, and Supplemental Figure 5, E and F). In addition, in xenografts, we found no change in the tumor-forming potential of M14 cells with GDF6 knockdown compared with that of control M14 cells with EGFP knockdown (Figure 2H). Therefore, whereas all cells examined showed increased tumorigenic potential upon GDF6 overexpression, only cells with amplification and higher levels of GDF6 protein were sensitive to GDF6 knockdown. These results suggest that GDF6 does not serve a housekeeping function, but rather that cells are dependent on its amplification and high expression for their tumorigenic potential.

BMP signaling is active in melanomas and is driven by GDF6 ligand. GDF6 genes, which encode BMP ligands, are predicted to act through SMAD1/5/8 transcription factors. For this reason, we investigated whether SMAD-dependent signaling was activated in melanomas. Using an antibody that specifically recognizes phosphorylated SMAD1/5/8 (p-SMAD1/5/8) proteins to monitor BMP pathway activity, we found robust p-SMAD1/5/8 nuclear staining in zebrafish melanomas (Figure 3A). Furthermore, transcriptome analyses of zebrafish melanomas indicated an upregulation of genes that support BMP signaling, including the BMPR1A and BMPR2 receptor subunits through which GDF6 is known to act (35) (Figure 3B). gdf6a and gdf6b were the only BMP ligands found to be upregulated in zebrafish melanomas, leading us to hypothesize that BMP signaling in melanomas is largely dependent on GDF6 (Figure 3C). To address this hypothesis, we modulated GDF6 activity in human melanoma cells. GDF6 knockdown caused a profound reduction in p-SMAD1/5/8 levels (Figure 3D), whereas GDF6 overexpression led to an increase in p-SMAD1/5/8 (Supplemental Figure 5H). p-SMAD1/5/8 proteins translocate to the nucleus, where, in complexes with SMAD4 and/or other regulatory proteins, they bind DNA and modulate the transcription of target genes. To determine whether GDF6 regulates SMAD1/5/8 DNA-binding activity, we performed ChIP and massively parallel sequencing (ChIP-seq) of p-SMAD1/5/8 in control and GDF6-knockdown A375 melanoma cells. Binding of p-SMAD1/5/8 to promoter regions of target genes was markedly reduced upon GDF6 knockdown (Figure 3E). Likewise, in a broader consideration of all p-SMAD1/5/8–bound regions, we found that knockdown of GDF6 caused a general reduction in binding (Supplemental Figure 6A). The reduction in binding was, in many cases, accompanied by transcriptional changes; for example, reduced binding and transcriptional downregulation co-occurred at the well-established p-SMAD1/5/8 target genes ID1 and ID3 (Supplemental Figure 6, B and C). These results indicate that BMP signaling is active in melanomas, and much of this activity is driven by GDF6.

Figure 3

GDF6-dependent BMP activity in melanomas. (A) Transverse sections of a Tg(mitfa:BRAFV600E);p53(lf) zebrafish bearing an invasive melanoma in the dorsal musculature. Top: H&E staining; bottom: p-SMAD1/5/8 staining. Scale bars: 500 μm (left) and 50 μm (enlarged insets on right). For p-SMAD1/5/8 staining, normal muscle (top) and a tumor region (bottom) are shown. Note that normal scale tissue (running vertically through the middle of the image) in the tumor region is p-SMAD1/5/8 negative. T, tumor; N, normal. (B) Heatmap of expression of BMP pathway genes (Reactome gene set R-HSA-201451.3; https://reactome.org/PathwayBrowser/) in zebrafish melanomas as compared with expression in melanocytes. Human orthologs of zebrafish genes are shown. (C) log₂-transformed fold change of gene expression in zebrafish melanomas as compared with expression in melanocytes (y axis). Expression of BMP ligands in microarray analysis and RNA-seq analysis. Only BMP ligands with a significant dysregulation (adjusted P < 0.05) are shown. (D) Immunoblots of p-SMAD1/5/8 and total SMAD1/5/8 in A375 melanoma cells expressing an shRNA targeting EGFP or 2 independent GDF6-targeted shRNAs. (E) Aggregation plot of p-SMAD1/5/8 ChIP-seq enrichment at annotated transcriptional start sites (TSSs) in A375 melanoma cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. The P value was calculated by 2-sample Kolmogorov-Smirnov (K-S) test after summing TSS proximal reads (–2kb to 2kb) for each gene (n = 49,344 TSSs). RPM, reads per million. (F) Tumor formation in mice injected with A375 cells (1 × 10⁶ cells injected per mouse) treated with vehicle control or 25 mg/kg DMH1 every other day. Error bars represent the mean ± SEM. n = 8. (G) Tumor formation in mice injected with A375 empty or A375-SMAD1DVD cells expressing 2 independent GDF6-targeted shRNAs. Each mouse was injected with 1 × 10⁷ cells. Error bars represent the mean ± SEM. n = 3. **P < 0.01 and ***P < 0.001, by 2-tailed Student’s t test (F) or 1-way ANOVA with Bonferroni’s test (G).

GDF6 acts via the BMP/SMAD pathway to promote tumor progression. To test whether the role of GDF6 in promoting melanoma progression is mediated by the SMAD1/5/8 axis of BMP signaling, we modulated the pathway activity in A375 cells. Knockdown of SMAD1 led to defects in cell growth and tumorigenic potential (Supplemental Figure 7, A–C), as was observed in GDF6-knockdown cells. We also used a small-molecule inhibitor of BMP signaling, DMH1, to block pathway activity. DMH1 can suppress the growth of BMP-dependent ovarian and lung cancer cells (36, 37), but its efficacy in melanoma has not been reported. DMH1 inhibits the kinase activity of ALK2 and BMPR1A (ALK3) receptors but not of BMPR1B (ALK6), thereby abrogating phosphorylation and activation of the SMAD1/5/8 transcriptional cascade (38). GDF6 has been shown to act through BMPR1A (ALK3) (35), and we found that treatment with DMH1 reduced p-SMAD1/5/8 levels and decreased cell growth and tumorigenicity (Figure 3F and Supplemental Figure 7, D and E). To test the relationship of GDF6 to SMAD1 in genetic epistasis analyses, we expressed a phosphomimetic variant of SMAD1, SMAD1DVD, in A375 cells (Supplemental Figure 7F) (39). Whereas GDF6 knockdown abrogated cell growth and the tumorigenic potential of A375 control cells, A375-SMAD1DVD cells subjected to GDF6 knockdown were rescued, showing robust growth in colony formation and xenotranplantation assays (Figure 3G and Supplemental Figure 7, G and H). The growth defects caused by treatment with DMH1 were also reversed by SMAD1DVD (Supplemental Figure 7I). Together, these data indicate that GDF6 signals via SMAD-dependent BMP signaling to promote tumorigenesis, and inhibition of this signaling achieves a reduction in tumor growth.

GDF6-dependent BMP signaling maintains a trunk neural crest gene signature. Since GDF6 acts through SMAD transcription factors, we were interested in identifying gene expression changes that could illuminate how this signaling axis regulates tumorigenesis. To do this, we modulated GDF6 and SMAD1 and sought genes commonly regulated by both. Given our genetic epistasis results, we predicted that the expression of important genes would change upon GDF6 knockdown but that such changes would be reversed when GDF6 knockdown was rescued by SMAD1DVD. Using RNA-seq, we defined the set of genes that was differentially regulated upon GDF6 knockdown and showed reciprocal differential regulation in SMAD1DVD-expressing cells that were subjected to GDF6 knockdown (Figure 4A and Supplemental Table 4). Pathway analysis showed that this gene set overlapped most significantly with genes involved in ossification and neural crest development (Figure 4B). GDF6-SMAD regulation of neural crest genes is intriguing, since melanocytes initially develop from this embryonic tissue. Several genes upregulated by GDF6 and SMAD1DVD — SOX10, TFAP2B, FOXD3, and SNAI2 — are neural crest “specifier” genes that are expressed broadly in the neural crest and help to maintain neural crest identity (Figure 4C) (40). SOX9 is initially broadly expressed in the neural crest, but as development proceeds, its expression becomes excluded from the trunk neural crest and limited to the cranial neural crest, from which mesenchymal tissues such as craniofacial cartilage are derived. Conversely, SOX10 expression becomes restricted to the trunk neural crest from which nonmesenchymal cells, including melanocytes, develop. Since SOX10 is upregulated and SOX9 is downregulated by GDF6/SMAD signaling, the pattern of gene regulation most closely resembles trunk neural crest tissue. Adopting a neural crest–like identity can contribute to the aggressive nature of melanoma cells (16, 41). For these reasons, we hypothesized that GDF6-driven BMP signaling (GDF6/SMAD), by promoting a trunk neural crest gene signature, enables melanoma cells to adopt and maintain an undifferentiated, protumorigenic state.

Figure 4

GDF6 and SMAD1 regulate a neural crest gene signature in melanomas. (A) Genes differentially regulated in A375 melanoma cells upon GDF6 knockdown (KD) (purple circle) and genes reciprocally regulated in SMAD1DVD-expressing A375 cells upon GDF6 knockdown (green circle). (B) Pathway analysis with the 605 reciprocally regulated genes (minimum overlap ≥10 genes; adjusted P < 0.01). Wikipathway analysis was done using Webgestalt toolkit (http://www.webgestalt.org/webgestalt_2013/). (C) Heatmap of neural crest genes identified in the pathway analysis.

GDF6 inhibits melanoma cell differentiation by repressing MITF. To determine how GDF6-driven BMP signaling could regulate the differentiation of melanoma cells, we considered target genes that were bound by p-SMAD1/5/8 and transcriptionally regulated in a GDF6-dependent manner. Among these genes was MITF, the master regulator of melanocyte differentiation. In control A375 cells, we detected p-SMAD1/5/8 binding in the MITF locus in a region that is intronic for longer MITF isoforms and upstream of the smaller MITF-M isoform, the predominant species of MITF in melanocytes and melanomas (Figure 5A) (42). Binding to this region was abrogated in GDF6-knockdown cells. This loss of binding in GDF6-knockdown melanoma cells was coupled with a transcriptional increase in MITF (Figure 5B). MITF itself is a transcriptional factor that orchestrates differentiation, in part, by activating the expression of melanin biosynthesis genes. The increase in MITF upon GDF6 knockdown was accompanied by an increase in,the melanin biosynthesis gene TRP1 (Figure 5C), indicating that melanocyte differentiation was invoked upon GDF6 loss. Upregulation of MITF and TRP1 was less pronounced when GDF6 knockdown was performed in SMAD1DVD-expressing cells (Figure 5, B and C). Depending on the cofactors involved, p-SMAD1/5/8 can promote transcription, as with ID1 and ID3, or repress transcription, as we propose for MITF (43). To determine whether GDF6 regulates the differentiation of melanomas in vivo, we examined tumors that arose in Tg(mitfa:BRAFV600E);p53(lf);gdf6a(lf) zebrafish. Melanomas from these animals had increased transcript levels of mitfa and trp1b, which are orthologs of the human MITF and TRP1 genes (Figure 5D). Furthermore, gdf6a-mutant tumors showed a profound increase in melanization as compared with control tumors (Figure 5E). These results show that GDF6 maintains melanomas in an undifferentiated state, and we speculate that preventing differentiation helps melanoma cells retain a neural crest identity.

Figure 5

GDF6-induced BMP signaling blocks melanoma cell differentiation. (A) p-SMAD1/5/8 binding to the MITF locus in A375 melanoma cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. Traces of 2 independent biological replicates are shown. (B) qRT-PCR analysis showing expression of MITF in A375 empty or A375-SMAD1DVD cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. (C) qRT-PCR analysis of TRP1 under the same conditions. Left brackets: MITF or TRP1 expression was upregulated upon GDF6 knockdown. Right brackets, MITF or TRP1 expression was less upregulated in SMAD1DVD-expressing cells upon GDF6 knockdown. Error bars represent the mean ± SEM. n = 3. (D) qRT-PCR analysis showing expression of mitfa and trp1b in control and gdf6a(lf) zebrafish melanomas. Error bars represent the mean ± SEM. n = 3. (E) H&E staining of transverse sections from Tg(mitfa:BRAFV600E);p53(lf) and Tg(mitfa:BRAFV600E);p53(lf);gdf6a(lf) zebrafish melanomas invading the dorsal musculature. Scale bars: 100 μm (left) and 25 μm (right). *P < 0.05 and ***P < 0.001, by 1-way ANOVA with Bonferroni’s test (B and C) or 2-tailed Student’s t test (D).

GDF6 represses SOX9 to promote melanoma cell survival. Knockdown of GDF6 and SMAD1 caused defects in melanoma cell growth. We sought to understand this defect and whether the trunk neural crest signature was involved. A variety of analyses suggested that the growth defect was linked to the regulation of apoptosis. Specifically, gene set enrichment analysis (GSEA) revealed that GDF6 expression was negatively correlated with expression of apoptotic pathway genes in cells with GDF6 modulation as well as in patients’ samples (Figure 6A and Supplemental Figure 8, A and B). In direct assessments, we found that GDF6 loss increased apoptotic cell death in A375 cells as well as in vivo in zebrafish melanomas and mouse xenografts (Figure 6, B and C, and Supplemental Figure 8, C–E). By contrast, GDF6 overexpression xenografts showed reduced apoptotic cell death as compared with basal levels of cell death in control xenografts (Supplemental Figure 8, F and G). GDF6 overexpression xenografts had a slightly increased Ki67-proliferative index, suggesting that the reduction in cell death was not caused by a failure to generate new cells with the potential to die (Supplemental Figure 8H). Finally, the cell death caused by GDF6 knockdown was rescued by SMAD1DVD expression, indicating that GDF6 acts via BMP signaling to promote melanoma cell survival (Figure 6D and Supplemental Figure 9). The involvement of GDF6 in cell death is consistent with findings in which loss of GDF6 orthologs in fish and Xenopus caused a substantial increase in apoptosis during eye and neural development (44–46). We speculate that GDF6 knockdown causes terminal differentiation of melanoma cells, leading to cell-cycle exit followed by cell death.

Figure 6

GDF6 and BMP signaling repress SOX9 to promote melanoma cell survival. (A) GSEA showed that expression of an apoptotic gene set (MSigDB M10169) was negatively enriched in GDF6-overexpressing A375 cells. (B) Caspase-3/7 activity measured as relative luciferase units (RLU) in A375 cells upon GDF6 knockdown. Error bars represent the mean ± SEM. n = 3. (C) Fluorescent TUNEL staining of Tg(mitfa:BRAFV600E);p53(lf) (top) or Tg(mitfa:BRAFV600E);p53(lf);gdf6a(lf) (bottom) zebrafish melanoma sections. TUNEL (green), DAPI (blue), and a merged image of both channels are shown. Scale bars: 25 μm. Error bars represent the mean ± SEM. n = 100 fields. (D) TUNEL staining of mouse xenografts of A375 cells expressing SMAD1DVD upon GDF6 knockdown. Scale bar: 25 μm. Error bars represent the mean ± SEM. n = 100 fields. (E) Immunoblots showing expression of SOX9 and GAPDH in A375 empty or A375-SMAD1DVD cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. (F) Caspase-3/7 activity measured as relative luciferase units (RLU) in A375-nonsilencing or A375-shSOX9 cells expressing an shRNA targeting EGFP or the GDF6-targeted shRNA GDF6.1. Error bars represent the mean ± SEM. n = 3. (G) Tumor formation in mice injected with A375-nonsilencing or A375-shSOX9 cells expressing 2 independent GDF6-targeted shRNAs. Each mouse was injected with 1 × 10⁶ cells. Error bars represent the mean ± SEM. n = 3. ***P < 0.001, by 1-way ANOVA with Dunnett’s test (B), 2-tailed Student’s t test (C), or 1-way ANOVA with Bonferroni’s test (D, F, and G). ^###P< 0.001, by 1-way ANOVA (G).

To determine how GDF6 and BMP signaling regulate melanoma cell survival, we focused on the reciprocally regulated genes defined previously. In particular, SNAI2 and SOX9 were assessed because of their importance in specifying neural crest and regulating cell survival (47–50). Whereas modulation of SNAI2 had no effect on GDF6-driven survival (Supplemental Figure 10), SOX9 was intimately involved. In studying SOX9, we first confirmed that changes in GDF6 and BMP signaling affected SOX9 expression. A375 cells with GDF6 knockdown had increased SOX9 RNA and protein levels, and this increase was much less pronounced in SMAD1DVD cells that were subjected to GDF6 knockdown (Figure 6E and Supplemental Figure 11A). Additionally, gdf6a-mutant zebrafish melanomas showed elevated sox9b levels (Supplemental Figure 11B). To determine whether GDF6 mediates cell survival by regulating SOX9, we measured whether knockdown of SOX9 (Supplemental Figure 11, C and D) could suppress the growth defects and cell death resulting from GDF6 knockdown. In colony formation assays, cells with combined knockdown of SOX9 and GDF6 grew much better than did GDF6 single-knockdown cells (Supplemental Figure 11E). Similarly, cell death, as measured by caspase-3 cleavage and annexin V positivity, was greatly reduced in GDF6/SOX9 double-knockdown cells (Figure 6F and Supplemental Figure 11F). Last, we measured the ability of SOX9 knockdown to rescue the tumor-forming capacity of GDF6-knockdown cells. When performed together with GDF6 knockdown, SOX9 knockdown enabled cells to engraft and rapidly grow (Figure 6G and Supplemental Figure 11G). These data indicate that a major function of GDF6 is to repress SOX9 expression, thereby inhibiting cell death and promoting tumor growth.

Clinical significance of GDF6 expression and signaling. We next assessed the expression of GDF6 in human melanomas and examined the potential clinical implications. IHC on an initial cohort of patients’ samples detected high levels of GDF6 protein in melanomas; however, normal melanocytes from adjacent skin (Figure 7A) or tumor-infiltrating cells (Supplemental Figure 12) rarely expressed GDF6. In the same cohort, we found high BMP pathway activity, as measured by nuclear p-SMAD1/5/8 staining, in tumor cells but not in normal tissue (Figure 7A and Supplemental Figure 12). To determine whether expression of GDF6 correlated with the clinical characteristics of melanoma patients, we performed IHC on a microarray with 104 melanoma tissue cores (78 primary and 26 metastatic melanomas). Consistent with the initial cohort, we observed robust GDF6 expression in a majority of melanomas (80% of total; n = 104 cases). Importantly, in analyzing clinical aspects of the melanoma tissue microarray, we found that patients whose tumors at diagnosis expressed high amounts of GDF6 had a lower survival probability than did those whose tumors expressed no or low GDF6 (Figure 7B and Supplemental Table 5). This association was mainly driven by patients with primary melanomas (Supplemental Figure 13A). Additionally, GDF6 expression in primary melanomas significantly correlated with lymph node metastasis (Figure 7C). In these primary melanomas, expression of GDF6 tended to be higher than in metastatic lesions (Supplemental Figure 13B). Together, these data indicate that GDF6 is a negative prognostic marker for early-stage melanomas.

Figure 7

Clinical impact of GDF6 expression and BMP pathway inhibition. (A) H&E, cytoplasmic GDF6, and nuclear p-SMAD1/5/8 immunostaining of adjacent normal skin and melanoma tissue from the same section. Melanocytes in normal skin sections are indicated by arrowheads. Images of individual cells are shown immediately to the right. Scale bars: 25 μm. Original magnification: ×63. Graphs indicate the percentage of patients’ samples with no or low expression or high expression of these proteins in normal melanocytes and melanomas. (B) Left, percentage of patients’ samples with no or low or high GDF6 expression in the melanoma tissue microarray. Graph of Kaplan-Meier analysis for the melanoma tissue microarray samples shows the overall survival of patients with no or low GDF6 expression (blue line) versus those with high GDF6 expression (red line). Statistical analysis was performed with a Mantel-Cox log-rank test. (C) GDF6 staining score in patients with primary melanomas with (n = 61) or without (n = 19) lymph node (LN) metastasis. **P < 0.01, by 2-tailed Welch’s test. (D) Mice bearing A375 xenografts were treated with vehicle, DMH1, dabrafenib plus trametinib, or a combination of all 3 drugs. Normalized tumor volumes following the beginning of drug treatments are shown. Error bars represent the mean ± SEM. n ≥ 8 animals. (E) Model for GDF6 activation and function in melanomas. **P < 0.01, by 1-way ANOVA with Bonferroni’s test (D); ^##P < 0.001 by 1-way ANOVA (D).

Next, we wanted to test whether targeting GDF6-driven BMP signaling could be combined with existing therapies. We treated established A375 melanoma xenografts with DMH1 or dabrafenib plus trametinib or a combination of all 3 drugs. While treatment with dabrafenib plus trametinib caused substantial regression, the combination of DMH1 with dabrafenib and trametinib showed even further regression (Figure 7D). Treatment with DMH1 alone had little effect, although we presume that it does not fully inhibit GDF6-driven BMP signaling in vivo, as indicated by our previous experiments (compare Figure 2H and Figure 3F). These results indicate that targeting GDF6-driven BMP signaling, in combination with current therapeutic strategies, may have profound clinical benefits.

The genomic and functional analyses performed in this study identified GDF6 as a melanoma oncogene. The results of these and additional mechanistic studies allowed us to develop a model for the role of GDF6 in melanoma tumorigenesis (Figure 7E). First, during the course of melanomagenesis, the GDF6 locus is transcriptionally activated, in some cases stemming from copy number amplification. The production of GDF6 protein then leads to either autocrine or paracrine activation of BMP pathway signaling in nascent melanoma cells. Next, BMP signaling triggers SMAD1/5/8-dependent downregulation of MITF and SOX9, which inhibits melanoma cell differentiation and death, respectively. GDF6-dependent binding of SMAD1/5/8 to the MITF locus suggests that the downregulation in this case is direct. There are no SMAD1/5/8-bound regions in the SOX9 locus, so its downregulation is either indirect or mediated by long-range interactions with SMAD1/5/8 that is bound, in a GDF6-dependent fashion, to sites in neighboring loci. The outcome of this signaling is to promote a neural crest gene signature within melanoma cells, which enables tumor progression.

GDF6 and BMP signaling in melanoma maintenance and initiation. BMP pathway activity has previously been implicated in melanoma progression. BMP4 and BMP7 expression has been shown to promote tumor cell invasion and migration in an autocrine fashion (51), whereas BMP2 and BMP4 promote angiogenesis through paracrine signaling (52). Independent studies indicate that BMP7 can block melanoma cell growth and inhibit metastasis (53, 54), suggesting that the effects of BMP7 on melanoma progression may be complex and potentially cell line specific. More relevant to our studies, BMP signaling has been proposed to promote melanoma survival by inducing the antiapoptotic factor DIDO1 (55). We tested whether DIDO1 is involved in the antiapoptotic role of GDF6 in melanoma cells. Upon GDF6 knockdown, we did not observe changes in DIDO1 transcript levels. Furthermore, we found that while p-SMAD1/5/8 bound to a site within the DIDO1 gene, this site was not differentially bound upon GDF6 knockdown. Thus GDF6, via SOX9, likely interfaces with the cell death machinery through other factors.

Along with preventing cell death, our studies indicate that GDF6 and BMP signaling promote melanoma maintenance through additional mechanisms. GDF6 downregulates the expression of MITF, which is likely accomplished by induction of p-SMAD1/5/8 and its direct binding to the MITF gene. MITF is a master regulator of the melanocyte lineage, and its regulation is critically important for the behaviors of cells within this lineage. MITF expression is governed by a variety of signals that ultimately produce an expression level that dictates cellular phenotype. In this way, MITF is proposed to act as a rheostat (41, 56, 57). Specifically, high levels of MITF promote cell-cycle arrest and terminal differentiation, whereas lower expression levels stimulate proliferation. Even lower levels endow cells with stem cell–like and invasive properties. We speculate that by repressing MITF, GDF6 and BMP signaling keep MITF expression in a range that not only inhibits terminal differentiation but also protects against cell death and endows cells with properties important for tumor maintenance.

Our data also implicate GDF6 and BMP signaling in melanoma initiation. In zebrafish, cells in the melanocyte lineage that expressed gdf6b gave rise to melanomas more rapidly than did control cells. Furthermore, melanoma onset was delayed in gdf6a-mutant zebrafish. Recently, it was shown that cells of origin for zebrafish melanomas adopt neural crest characteristics that distinguish them from normal melanocytes (32). These neural crest characteristics are proposed to be instrumental in melanoma initiation. Given that GDF6 and its orthologs are melanoma oncogenes that promote a neural crest gene signature, we speculate that GDF6 and BMP signaling are important at the earliest stages of melanoma formation.

Reiteration of embryonic GDF6 activities in melanomas. GDF6 is expressed during embryogenesis, and its functions during development mirror those upon its reactivation in melanomas. In Xenopus and zebrafish embryos, GDF6 is expressed at the edges of the neural plate and in the eye fields of the anterior neural plate. Upon neural tube closure, expression is prominent in the dorsal neural tube and neural crest (29, 31, 45). Loss-of-function studies in both species suggest that GDF6 promotes ectodermal cell survival (45, 58). In Xenopus and zebrafish, knockdown of GDF6 and gdf6a, respectively, caused a reduction in eye size that resulted from inappropriate death of retinal neuron progenitor cells (44, 46, 59). GDF6 also acts during embryogenesis to regulate cell differentiation. In the mouse, Gdf6 inhibits differentiation of the mesenchymal progenitors that develop into the coronal suture, and precocious differentiation of these cells results in fusion (60, 61). Thus, in certain tissues, GDF6 can promote cell survival during development as well as regulate terminal differentiation.

We have found that GDF6 promotes a neural crest signature in melanoma cells. In particular, GDF6 maintains expression of the trunk neural crest factor SOX10, while repressing the cranial neural crest factor SOX9. Some data support a similar function for GDF6 in embryogenesis. As noted above, GDF6 orthologs in Xenopus and zebrafish are expressed in the neural crest and adjacent dorsal neural tube. In zebrafish, knockdown of gdf6a reduces expression of sox10 in neural crest cells (31), consistent with a function in maintaining a trunk neural crest gene signature. As development proceeds, gdf6 paralog expression is progressively lost in a rostrocaudal direction, which occurs concomitantly with rostrocaudal differentiation of crest cells, including melanocytes. Perhaps the loss of gdf6 expression removes a barrier to melanocyte differentiation. Support for such a barrier has been shown in avian embryos, where BMP4-driven pathway activity inhibits melanocyte fate specification during neural crest lineage segregation (62). We speculate that GDF6 may activate BMP signaling in the neural crest to repress MITF. With MITF repressed, neural crest cells can remain undifferentiated and responsive to proliferative signals as well as those that specify alternative cell fates. In melanomas, BMP signaling likewise could promote a neural crest cell identity that is not terminally differentiated and therefore conducive to proliferation.

Targeting GDF6 and BMP signaling in melanoma. Our patient cohort analyses show that a major fraction of melanomas express GDF6 and have an active BMP pathway. Since GDF6 expression is correlated with poor patient survival, and inhibition of GDF6 leads to cell death, targeting this gene or the BMP pathway may prove to be an effective therapeutic intervention for melanomas. GDF6 itself is an attractive target, since its expression is very low or undetectable in most adult tissues (63). The knockout of mouse Gdf6 indicates that it is a nonessential gene, although it is necessary for the development of certain joints, ligaments, and cartilage (61). These developmental activities of GDF6 would likely not complicate the treatment of adult patients with anti-GDF6 therapy, although it is not known whether GDF6 inhibition would affect the limited repair of connective tissues that occurs in adults. On the other hand, BMP signaling regulates the function or repair of some adult tissues, including muscle, bone, and lung (64–66). Thus, GDF6 inhibition could be more specific to melanoma tissue as compared with a broad BMP pathway inhibitor such as DMH1, which would block BMP signaling in normal tissues. As a secreted molecule, GDF6 inhibition could be accomplished by a variety of means, including cell-impermeable therapies. Such therapies could be used in different molecular subtypes of melanoma, as cell death caused by GDF6 inhibition does not depend on underlying BRAF, NRAS, or other driver mutations. Our data do suggest that tumors expressing high levels of GDF6 would be particularly sensitive to GDF6 inhibition. This elevated expression correlates with GDF6 gene amplification, which is present in 38% of melanomas (67). However, 80% of the melanomas we analyzed showed high GDF6 expression, suggesting that mechanisms other than amplification can lead to higher expression. Ultimately, it would be useful to determine whether the level of GDF6 expression is predictive of therapeutic response to a GDF6 or BMP inhibitor. As indicated by the potent activity against xenografts of DMH1 with dabrafenib plus trametinib, such an inhibitor treatment could be used in conjunction with BRAF inhibitors and other MAPK pathway inhibitors to treat this lethal disease.

Additional methods are described in the Supplemental Methods.

Cell lines and cell culture. A375, M14, C32, and HEK293T cells (ATCC) were maintained in DMEM, and SK-MEL-28 and SK-MEL-5 cells (ATCC) were maintained in RPMI 1640 media supplemented with 10% FBS and 2 μg/ml penicillin-streptomycin (Gibco, Thermo Fisher Scientific) at 37°C and 5% CO₂. Cells cultured at the same time were pooled, counted, and then seeded in a 10-cm dish. Wells and dishes were then subjected to treatment with lentiviral vectors.

Lentiviral infection. Lentiviral infections were performed as described previously (28). For stable gene knockdowns, we used pLKO-1 lentiviral vectors to deliver shRNAs (obtained from the RNAi Consortium [TRC]/Broad Institute through the UMMS RNAi core facility) specific for GDF6 (GDF6.1, TRCN0000141818, target sequence: GCCAAGTGTTACATTGAGCTT; GDF6.2, TRCN0000140097, target sequence: GTGTCCATGCTCTCAGACAAA) or SMAD1 (SMAD1.1, TRCN0000021781, target sequence: CGGTTGCTTATGAGGAACCAA; SMAD1.2, TRCN0000021782, target sequence: GCCGATGGACACAAACATGAT) or EGFP (TRCN0000072181, target sequence: ACAACAGCCACAACGTCTATA). Virus was made using a second-generation lentiviral packaging system in HEK293T cells and quantified using a p24 ELISA Kit (Clontech). Cells were infected with virus at a MOI of 2.5, with 8 μg/ml polybrene followed by puromycin selection (2 μg/ml) for 2 days in appropriate media. For genetic epistasis experiments with SOX9, we used the pGIPZ lentiviral vectors (obtained from Thermo Fisher Scientific through the UMMS RNAi Core facility) to deliver shRNAs specific for SOX9 (V3LHS_396212, target sequence: AGTCGTACTGTGAGCGGGT) or used the nonsilencing control (target sequence: CTTACTCTCGCCCAAGCGAGAG). A375 melanoma cells expressing an shRNA targeting GDF6 or EGFP were treated with virus delivering SOX9 or nonsilencing shRNA. The viral dosage was determined such that 100% of the cells were EGFP positive and therefore contained the pGIPZ vector expressing either SOX9 or nonsilencing shRNA. For transgene expression, we used Gateway cloning (Life Technologies, Thermo Fisher Scientific) to insert the GDF6 or SNAI2 ORF (GE Life Sciences) or the SMAD1DVD ORF (provided by Takenobu Katagiri, Saitama Medical University, Saitama, Japan) into the pLenti CMV Hygro DEST (w117-1) vector (provided by Paul Kaufman, UMMS, Worcester, Massachusetts, USA). Infection and monitoring was performed as previously described (28), except that selection was done with 300 ug/ml hygromycin for 10 days.

Mouse xenografts. A375 cells stably expressing an EGFP, GDF6, or SMAD1 shRNA and/or an empty vector or a GDF6-expressing or SMAD1DVD-overexpressing vector were s.c. injected into the flanks of 6- to 8-week-old BALB/c nu/nu female mice (Taconic Farms) to produce orthotopic primary tumors. Primary tumor growth was monitored every 3 days with calipers, and tumor volume was calculated as described previously (68). For GDF6 knockdown, SMAD1 knockdown, and epistasis experiments with SMAD1DVD overexpression, 1 × 10⁷ live cells were injected. For GDF6 overexpression experiments and epistasis experiments with SOX9 knockdown, 1 × 10⁶ live cells were injected. For GDF6 overexpression and GDF6 knockdowns, a representative result of 2 independent experiments (n = 3 animals per experiment) is shown. For the DMH1 drug experiments shown in Figure 3F, 1 × 10⁶ live A375 cells were s.c. injected into the flanks of BALB/c nu/nu female mice. Beginning on the day cells were injected, mice were injected i.p. with vehicle (12.5% 2-hydroxypropyl-β-cyclodextrin) or 25 mg/kg DMH1 in vehicle every other day. This experiment was repeated twice, and the weighted average of both experiments (n = 8 animals total) is represented. For the drug experiments depicted in Figure 7D, 1 × 10⁶ A375 cells were s.c. injected into the flanks of BALB/c nu/nu female mice. Once the tumor volumes reached 75 mm³, the groups of mice were treated with vehicle, DMH1, dabrafenib plus trametinib, or a combination of all 3 drugs. DMH1 was administered i.p. in 7-day cycles, with 25 mg/kg given twice a day for 5 days and once a day for 2 days. Dabrafenib (300 mg/kg) plus trametinib (0.3 mg/kg) were administered by oral gavage once per day. The control group received the vehicle 12.5% 2-hydroxypropyl-β-cyclodextrin i.p. in 7-day cycles, with 25 mg/kg given twice a day for 5 days and once a day for 2 days. Control animals also received 0.5% hydroxypropyl methylcellulose once a day by oral gavage. These treatment regimens were continued until the end of the study or until tumors reached 400 mm³ in size. Eight or nine animals were used in each group.

Zebrafish stocks and husbandry. Zebrafish were maintained at 28.5°C with a 14-hour on/10-hour off light cycle. Tg(mitfa:BRAF^V600E) (23), p53(lf) (69), mitfa(lf) (70), and gdf6a(lf) (33) zebrafish strains were used. gdf6a(lf) mutants were provided by Herwig Baier (Max Plank Institute of Neurobiology, Martinsried, Germany). AB was used as the WT strain.

miniCoopR assay. The miniCoopR assay measuring the effect of gdf6b on melanoma onset in zebrafish was performed as previously described (28). For miniCoopR-EGFP experiments, a weighted average of 2 independent experiments is represented, and for miniCoopR-gdf6b experiments, a weighted average of 4 independent experiments is represented. For the embryonic melanocyte rescue analysis, embryos were treated with epinephrine (1 mg/ml) 4 days after fertilization to contract pigment to distinguish any overlapping cells, and melanocytes were counted manually.

cDNA synthesis and qRT-PCR. For adult zebrafish, total RNA was extracted from melanoma cells and from normal scale-associated melanocytes of Tg(mitfa:BRAFV600E);p53(lf);alb(lf);Tg(mitfa:EGFP) zebrafish. For isolation of melanoma cells, melanomas were dissected, dissociated using Liberase TH treatment (MilliporeSigma), and subjected to FACS to isolate EGFP-positive cells. The same protocol was used for normal melanocytes, except dorsal scales from zebrafish were plucked to isolate melanocytes. Total RNA from zebrafish melanomas and melanocytes was isolated using TRIzol-chloroform extraction, followed by RNA clean-up using an RNeasy Kit (QIAGEN). Total RNA was reverse transcribed using the Superscript 2 Reverse Transcriptase Kit (Invitrogen, Thermo Fisher Scientific). Quantitative reverse transcription PCR (qRT-PCR) with SYBR Green Master Mix (Bio-Rad) was performed using the following primers: gdf6a, forward, CTGAGAAACTGGGGCTCAAT; gdf6a, reverse, CGACCAGCTCCTCTTTGTCT; gdf6b, forward, CGTCTAAAGCAGCAAACACC; gdf6b, reverse, CCAAAGTGGAGAGTTCAAATGG; actb1, forward, CGAGCAGGAGATGGGAACC; actb1, reverse, CAACGGAAACGCTCATTGC. qRT-PCR was performed using the following primers: mitfa, forward, CTGGACCATGTGGCAAGTTT, mitfa, reverse, GAGGTTGTGGTTGTCCTTCT; tyrp1b, forward, CGACAACCTGGGATACACCT, tyrp1b, reverse, AACCAGCACCACTGCAACTA.

For A375 human melanoma cells with GDF6 and/or SMAD1DVD modulation, total RNA was prepared in the same manner, and qRT-PCR was performed using the following primers: ID1, forward, CCAACGCGCCTCGCCGGATC; ID1, reverse, CTCCTCGCCAGTGCCTCAG; ID3, forward, CTGGACGACATGAACCACTG; ID3, reverse, GTAGTCGATGACGCGCTGTA; SNAI2, forward, TGTTGCAGTGAGGGCAAGAA; SNAI2, reverse, GACCCTGGTTGCTTCAAGGA; SOX9, forward, GTACCCGCACTTGCACAAC; SOX9, reverse, TCTCGCTCTCGTTCAGAAGTC; MITF, forward, AAACCCCACCAAGTACCACA; MITF, reverse, ACATGGCAAGCTCAGGAC; TRP1, forward, GTAACAGCACCGAGGATGG; TRP1, reverse, TCCAAGCACTGAGCGACAT; GAPDH, forward, TGCACCACCAACTGCTTAGC; GAPDH, reverse, GGCATGGACTGTGGTCATGAG.

GSEA and pathway analysis. For GSEA, the enrichment score (ES), normalized enrichment score (NES), and FWER were calculated on the basis of a running metric, which increased when a gene (vertical line in the graphical representation) in the gene set was encountered and decreased when one was not. For GSEA of the apoptotic pathway gene signature (Molecular Signatures Database [MSigDB] M10169) (71), a rank-ordered gene list was made with fragments per kilobase per million mapped reads (FPKM) values from GDF6-overexpressing A375 melanoma cells as compared with empty vector control cells or from A375 cells expressing an shRNA targeting EGFP as compared with GDF6.1 shRNA-expressing cells. Default parameters of GSEA were used, and the Student’s t test was used to calculate significance. For GSEA based on TCGA samples, a rank-ordered gene list was derived from the expression profiles of 385 melanoma samples, using the GDF6 expression level as a continuous variable. Default parameters of GSEA were used, and Pearson’s correlation was used to calculate significance. Pathway analysis was performed using the WEB-based Gene SeT AnaLysis Toolkit (WebGestalt) (72). Default parameters were used, except the minimum number of genes for a category was set at 10.

Accession information. The RNA-seq, ChIP-seq, and gene expression microarray data sets reported in this article were deposited in the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE83400).

Statistics. Data are presented as the mean ± SEM. Differences between groups were assessed by a 2-tailed Student’s t test, except in Figure 2, G and H (left), Figure 3G, Figure 5, B and C, Figure 6, B, D, F, and G, Figure 7D, Supplemental Figure 4, A and B, Supplemental Figure 6C, Supplemental Figure 7, B, G, and I, Supplemental Figure 8C, Supplemental Figure 9, A, B, and C, Supplemental Figure 10, A–D, and Supplemental Figure 11, A, C, E, and F, where 1-way ANOVA was used (GraphPad Prism 7; GraphPad Software). In Figure 7C, Supplemental Figure 1C, and Supplemental Figure 13B, a 2-tailed Welch’s t test was used (GraphPad Prism 7). In Figure 3E, a 2-sample Kolmogorov-Smirnov test was used. In Figure 1C, Figure 3C, Figure 4A, Figure 6A, and Supplemental Figure 8, A and B, a FWER P value was used to account for multiple comparisons. For multiple comparisons, a Dunnett’s or Bonferroni’s test was used as appropriate. Kaplan-Meier curves were generated using GraphPad Prism 7 (GraphPad Software), and the differences between groups were assessed by Wilcoxon rank-sum analysis for Figure 2, A and H, and by Mantel-Cox log-rank analysis for Figure 7B and Supplemental Figure 12. A P value of less than 0.05 was considered significant.

Study approval. All zebrafish and mouse study protocols were approved by the IACUC at UMMS. Approval for the tissue microarray study was granted by the Brigham and Women’s Partners Human Research Committee. Informed consent was not necessary, as all tissue samples were discardable and deidentified.