Melanoblast transcriptome analysis reveals novel pathways promoting melanoma metastasis

Cutaneous malignant melanoma is an aggressive cancer of melanocytes with a strong propensity to metastasize. We posited that melanoma cells acquire metastatic capability by adopting an embryonic-like phenotype, and that a lineage approach would uncover novel metastatic melanoma biology. We used a genetically engineered mouse model to generate a rich melanoblast transcriptome dataset, identified melanoblast-specific genes whose expression contributed to metastatic competence, and derived a 43-gene signature that predicted patient survival. We identified a melanoblast gene, KDELR3, whose loss impaired experimental metastasis. In contrast, KDELR1 deficiency enhanced metastasis, providing the first example of different disease etiologies within the KDELR-family of retrograde transporters. We show that KDELR3 regulates the metastasis suppressor, KAI1, and report an interaction with the E3 ubiquitin-protein ligase gp78, a regulator of KAI1 degradation. Our work demonstrates that the melanoblast transcriptome can be mined to uncover novel targetable pathways for melanoma therapy.

melanoblast transcriptome. Using this approach, we identified a 43-gene embryonic melanoblast gene signature that predicts metastatic melanoma patient survival, and we highlight a new role for KDELR3 20 , distinct from other members of the KDELR family. A metastasis suppressor screen highlights KAI1/CD82 (hereafter referred to as KAI1) as a KDELR3-regulated protein. We observe that KDELR3 regulates KAI1 protein levels and post-translational modification. We demonstrate an undescribed interaction of KDELR3 with gp78, the E3 ubiquitin protein ligase known to regulate KAI1 degradation 23 . Our work shows that melanoma cells can commandeer embryonic transcriptomic programs to promote their progression to metastasis. These genes represent an untapped source of novel targetable pathways to exploit for improving melanoma treatment.

Melanoblast transcriptomic expression in melanoma metastasis
To study melanoblast genes, GFP-positive melanocytic cells were isolated from four developmental time points: Embryonic day (E) 15.5 and 17.5, and Postnatal day (P) 1 and 7 (Fig. 1b, Supplementary Fig. 1a-b). These four stages represent embryonic melanoblast development from the neural crest into differentiated quiescent melanocytes of the postnatal pup 24,25 .
Melanocytes/ melanoblasts were isolated using Fluorescence-Activated Cell Sorting (FACS) from iDct-GFP mice ( Supplementary Fig. 1c). At E15.5 and E17.5 melanoblasts are still migrating and colonizing the hair follicles within the epidermis [24][25][26] − processes that we believe are highly relevant to metastasis, particularly with respect to colonization at the metastatic site -and intrafollicular melanoblasts are still present 26 . P1 and P7 mature melanocytes were selected as a model of quiescent differentiated melanocytes; these time points are prior to the first hair follicle cycle that begins at 6 weeks post birth. Melanocytic cells were extracted from multiple litters (6-10 pups) at each developmental stage to ensure comprehensive representation of all melanoblasts/ melanocytes present. RNA was extracted for whole-transcriptome sequencing.
Genes with differential expression between embryonic melanoblasts (E15.5 and E17.5) and postnatal differentiated melanocytes (P1 and P7) were identified using DESeq2 27 with a q-value < 0.1, yielding 976 differentially expressed genes ( Supplementary Fig. 2). Of these genes, we filtered out any whose differential expression was less than 1.5 log2 fold increased in melanoblasts, as we deemed that a fold change of less than this was unlikely to be biologically meaningful. 467 melanoblast-specific genes were identified from our analyses, which we hypothesize to be putative melanoma metastasis enhancer genes (MetDev genes; Fig. 1c). To test the relevance of our melanoblast gene cohort in melanoma metastasis we interrogated this gene list in melanoma patient data. To ask if our 467-gene MetDev cohort was enriched in genes that contributed to poor progression of patients, we used a Cox proportional hazards model to associate their expression with overall survival in a training dataset of human patient samples derived from melanoma metastases (stage III and stage IV; GSE19234) 28 . We discerned a 43-gene survival risk predictor ( Fig. 1c, black/red arrows; Fig. 1d, black text, Kdelr3) that could accurately predict patient outcome in a separate testing dataset of late stage (stage III and stage IV) metastatic melanoma patient samples derived from metastases (GSE8401; Fig. 1e) 29 . These data show that not only is our MetDev cohort enriched for metastatic progression genes, but it can also predict survival in multiple independent patient datasets. Notably, gene expression levels in samples derived from early stage (stage I and stage II) primary melanoma lesions did not predict patient outcome, suggesting that MetDev genes play a key role in late-stage disease specifically (GSE8401; Fig.   1f) 29 .
To allow functional validation of our MetDev candidates in both soft agar colony forming assays and in experimental metastasis models we elected to prioritize the list of MetDev gene candidates.
To do this in an unbiased fashion we applied criteria based solely on melanoblast expression data, selecting for genes with no detectable gene expression in P7 postnatal pups. Differential expression was validated using a separate microarray expression dataset derived from our iDct-GFP model (E17.5 vs P2 and P7; q-value < 0.1) 21 . Further criteria using differences in fold-increase expression in melanoblasts vs. melanocytes and the greatest expression at embryonic stages allowed us to select 20 genes likely to be the most functionally relevant. Of these 20 we noted that 7 genes (Kdelr3, P4ha2, Gulp1, Dab2, Lum, Aspn, Mfap5) were associated with Extracellular Matrix (ECM) or trafficking. For functional analyses, we chose 4 of these 7 genes (Kdelr3, P4ha2, Gulp1, Dab2) with no established role in cutaneous melanoma metastasis (Fig. 1c, green/red arrows; Fig.   1d, red text). siRNA knockdown of our four candidate genes in B16 mouse melanoma cells inhibited both growth in soft agar colony formation assays and formation of lung metastases in experimental metastasis assays compared to non-targeting controls (Table 1). Our work demonstrates that the MetDev dataset is enriched in genes that have a functional role in melanoma metastasis. We identify four new melanoma metastasis genes and highlight ECM and trafficking as important pathways common to both melanoblast development and melanoma metastasis.
We further observed significant co-expression of three of the four functionally validated genes  Table   2). Notably, expression of Kdelr3 and P4ha2 was highly correlated throughout all datasets (Supplementary Fig. 3a-b), raising the possibility that some metastasis-associated MetDev genes may be co-regulated and serve a more coordinated role in metastasis.

KDELR3 is a Golgi-resident protein whose expression correlates with melanoblast development and melanoma progression
To understand how melanoblast genes might facilitate metastasis we chose to study one MetDev gene in depth. KDELR3 was selected as it was a positive hit in all of our analyses: KDELR3 is a trafficking protein important in the ERSR whose expression was associated with poor patient prognosis in metastatic melanomas (Fig. 1e, 43 gene signature), and KDELR3 was functionally validated in both soft agar colony formation and experimental metastasis assays ( Table 1). The KDELRs are Golgi-to-ER retrograde transporters responsible for maintaining ER localization of their protein substrates. KDELR substrates consist of protein chaperones required for protein folding and targeting unfolded proteins for degradation 20 , thereby assisting the UPR and maintaining ER quality in times of ER stress. We show that KDELR3 is localized to both the cisand trans-Golgi compartments in metastatic melanoma cells ( Supplementary Fig. 3c) and validate expression of KDELR3 in mouse melanoblasts (Fig. 2a). Moreover, within the KDELR family only KDELR3 demonstrated a melanoblast-specific expression pattern and showed consistent upregulation in melanoma cell lines ( Fig. 2b; Supplementary Fig. 3d-e). These data raise the possibility that KDELR3 plays a role in melanoma progression which is distinct from other KDELRs, despite their presumed redundancy. Analysis of human patient datasets and tumor histology microarrays confirmed an upregulation of KDELR3 expression in malignant melanoma vs. benign nevi (Fig. 2c-e).
We sought to functionally validate a role for KDELR3 in melanoma progression. We used human and mouse melanoma cells to demonstrate that small-interfering RNA (siRNA) and short-hairpin RNA (shRNA) knockdown of KDELR3 significantly reduced, and KDELR3 overexpression enhanced, anchorage-independent growth (Fig. 3a-d; Supplementary Fig. 4a-b), which cannot be attributed to a change in proliferation ( Supplementary Fig. 4c). There are two KDELR3 variants, and we selected the KDELR3-001 variant to perform rescue experiments as it is the most abundant transcript expressed in human cell lines and patient samples. We therefore performed rescue experiments via exogenous expression of KDELR3-001 Mu , whose shRNA recognition site had been mutated without altering the final protein sequence. KDELR3-001 Mu expression was restored, rescuing the anchorage-independent growth phenotype ( Fig. 3e-g; Supplementary Fig. 4d).
KDELR3 was therefore validated as a mediator of anchorage-independent growth in melanoma cells, a process required for metastasis.

KDELR3 knockdown reduces lung colonization in experimental metastasis assays
To assess the relevance of KDELR3 within the metastatic cascade, we used a tail vein experimental metastasis assay, which specifically assesses the ability of the cells to extravasate and colonize the lung, processes that are critical for metastatic capacity. Tail vein metastasis assays enable lung colonization to be assessed with greater specificity/sensitivity − biology that we suggest may be mirrored during hair follicle colonization (E17.5). Transient knockdown of KDELR3 in either mouse ( Fig. 3h-i) or human melanoma cell lines (Fig. 3j, Supplementary Fig. 5a) resulted in significantly reduced metastatic potential compared to non-targeting controls, indicating that KDELR3 expression is important for the cells' ability to extravasate/ colonize the lung, further validating that KDELR3 is a melanoblast gene that functions in metastasis (MetDev gene). Stable shRNA knockdown of KDELR3 also resulted in a reduction in lung colonization following tail vein metastasis and significantly fewer mice characterized with high metastatic burden (Supplementary Fig. 5b-f). However, no appreciable difference in cell cycle or subcutaneous in vivo tumor growth was observed ( Supplementary Fig. 5g-i), suggesting that the KDELR3-mediated metastatic phenotype cannot be attributed to a change in proliferation, and that KDELR3 is a genuine melanoma metastasis progression gene.

KDELR3 and the ER Stress Response in metastatic melanoma
To uncover how KDELR3 expression may be involved with melanoma metastasis, we asked which pathways were co-regulated with KDELR3 expression. Gene Set Enrichment Analysis (GSEA, FDR < 0.0001) of KDELR3 co-expressed genes in TCGA skin cutaneous melanoma patients (cBioPortal) 30,31 , revealed Gene Ontology (GO) term enrichment of ECM and trafficking pathways (consistent with previous data, Table 1, Supplementary Fig 2a), and pathways involved in the ERSR and response to unfolded proteins ( Supplementary Fig. 6a). Quantitative mass spectrometry was used to analyze whole cell lysates of KDELR3 knockdown compared to nontargeting controls and parental controls; GSEA analysis revealed the top-scoring, most consistent pathway using GO term enrichment showed upregulation of ER lumen proteins ( Supplementary   Fig. 6b). Enriched proteins included protein chaperones, lectins, and enzymes involved in protein folding and targeting misfolded proteins for degradation (including UGGT, ER Lectin, FKBP7, Calumenin), which is consistent with an increase in misfolded protein load in KDELR3 knockdown cells 32 We therefore asked how KDELR3's role in the ERSR response is associated with its metastasis phenotype. Metastasis is known to be linked with ER stress, activating the UPR and therefore downstream signaling events that function to alleviate this stress 17 . High doses of ER stress, or an ineffective UPR have been associated with deleterious signals and ultimately cell death. We therefore hypothesized that one role of KDELR3 in metastasis would be to alleviate ER stress-induced deleterious signaling ( Supplementary Fig. 6c). We observed in four independent mouse models of melanoma (N = 6-13 mice per model) that Perk (Eif2ak3) transcription was negatively correlated with Kdelr3 transcription (Fig. 4a), whereas Gadd34 (Ppp1r15a) transcription was positively correlated (Fig. 4b). As PERK is a protein kinase and GADD34 a protein phosphatase that both act on EIF2α 33 , we hypothesized that KDELR3 -low cells are primed to activate the PERK-EIF2α arm of the UPR. We knocked down KDELR3 (KD) in both 1205Lu and WM-46 human cell lines (shRNA knockdown; Supplementary Fig. 5b) and found that loss of KDELR3 expression resulted in increased PERK and EIF2α protein levels in untreated cells, corroborating our mouse model data (Fig. 4c). We also saw a concomitant increase in PERK and EIF2α phosphorylation, suggesting constitutive activation of the PERK-EIF2α axis in untreated KD cells (Fig. 4c). The other two branches of the UPR pathways, the IRE1-XBP1 and ATF6α axes, were inactive in untreated KDELR3 KD cells ( Supplementary Fig. 6d-e). Tunicamycin, a chemical inhibitor of N-glycosylation that induces ER stress in cells, was used as a positive control  19 . These data indicate that loss of KDELR3 expression disrupted ER homeostasis, resulting in a dysregulated UPR, which has previously been linked with ER stress-associated cell death 34. We hypothesized that KDELR3 functions to alleviate deleterious ER stress-induced signaling ( Supplementary Fig. 6c). To test this, we asked if KDELR3 knockdown sensitizes metastatic melanoma cells to ER stress-induced death. We treated cells with tunicamycin, and measured cell death through flow cytometry using Live/Dead cell stain. We observed that siRNA-mediated knockdown of KDELR3 expression resulted in a ~5-fold increase in metastatic melanoma cell death over controls (8.3%, siKDELR3; 1.6%, siControl; Fig. 4d). These data suggest that KDELR3 promotes cell survival in metastatic melanoma cells, which likely influences metastatic potential.
KDELR3-knockdown cells have an enhanced sensitivity to ER stress induction with tunicamycin (>13-fold difference in cell death: 28.4%, siKDELR3; 2.1%, siControl; Fig. 4d). These data indicate that the ability of KDELR3 to relieve ER stress is crucial for adaptation and survival in metastatic melanoma and may be instrumental to the metastatic phenotype.

KDELR3 mediates post-translational regulation of the metastasis suppressor KAI1
To further understand the role of KDELR3 in metastasis, we queried if KDELR3 knockdown would increase expression of known metastasis suppressors in melanoma. To address this, we screened protein expression of 5 melanoma metastasis suppressors following KDELR3 knockdown 35,36 . Of the 5 metastasis suppressors screened (BRMS1, Gelsolin, GAS1, NME1/NM23-H1, KAI1) only KAI1 demonstrated an increase in expression following KDELR3 knockdown (Fig. 5a), in line with our hypothesis. Moreover, we observed a change in KAI1 molecular weight distribution following KDELR3 knockdown, suggesting alterations in KAI1 post-translational modification.
KAI1 protein upregulation was independent of transcriptional changes (Fig. 5b), supporting a regulatory role for KDELR3 at the post-translational level. KAI1 has been shown to influence metastasis through multiple mechanisms, including cell-cell adhesion, cell motility, cell death and senescence, and protein trafficking in many cancer types, including melanoma 37 . To further validate the role of KDELR3 on KAI1 protein regulation, we exogenously expressed KAI1 protein in 1205Lu metastatic melanoma cells (in which endogenous KAI1 expression is relatively low) and co-expressed both KDELR3-001 and KDELR3-002. Corroborating our initial findings, we found that increased KDELR3 expression resulted in dramatically reduced KAI1 protein levels ( Fig. 5c), which could not be accounted for by KAI1 transcriptional changes ( Fig. 5d-e). KAI1 protein glycosylation pattern was impacted reciprocally by knockdown and overexpression experiments, supporting the notion that KAI1 post-translational modification pathways are regulated by KDELR3, including an upregulation of a high molecular weight band in KDELR3 knockdown cells (Fig. 5f, red arrow) that we showed corresponds to a highly glycosylated form of KAI1 (Fig. 5g). Glycosylated KAI1 has been linked to inhibition of cell motility and promotion of cell death 38 , and has been shown to influence N-cadherin clustering and bone metastasis in AML 39 .
Owing to our protein expression data, we hypothesized that KDELR3 regulates KAI1 protein degradation. We asked if KDELR3 regulates expression of the E3 ubiquitin ligase known to target KAI1, gp78/Autocrine Motility Factor Receptor 23, 40 , hereafter gp78. Although we saw no significant alterations in gp78 protein or RNA expression following KDELR3 knockdown (Fig. 5f, h), we did observe a 3-fold increase in KDELR3 transcription following gp78/AMFR knockdown, suggestive of a functional link between the two proteins (Fig. 5i). We identified a previously undescribed interaction between KDELR3 and gp78, which was supported by evidence of colocalization ( Fig. 5j-k; Supplementary Fig. 7a). Interestingly, gp78 was first identified as a motility factor associated with metastasis in several cancers 41 , including melanoma. We asked if the KDELR3-gp78 interaction impacted gp78 function. We reasoned that gp78 ubiquitin ligase substrates would be upregulated following gp78 knockdown, as these proteins would not be targeted for degradation; however, not all upregulated proteins identified will be direct gp78 substrates. Quantitative mass spectrometry was used to analyze whole cell lysates of gp78 (AMFR) knockdown or KDELR3 knockdown cells compared to non-targeting controls. We could confirm that 43-57% of upregulated proteins matched between the gp78 and KDELR3 knockdown groups.
GSEA showed that the top-scoring, upregulated pathways (FDR <0.05) for both groups using GO term enrichment were those associated with the ER (Supplementary Table 3 -4). This result suggests that both gp78 and KDELR3 act within similar cellular pathways and supports a role for KDELR3 in gp78 function, highlighting at least one mechanism through which KDELR3 can influence metastasis at the post-translational level. Since gp78 is a ubiquitin ligase known to function in ERAD, our data link KDELR3 to ERAD regulation. In summary, our work implicates KDELR3 in glycosylation of the metastasis suppressor, KAI1, and in its degradation through gp78 (and likely other ERAD effectors), thereby providing a mechanism for KDELR3's influence on the metastatic phenotype (Fig. 5l).

KDELR3 correlates with late-stage metastasis and poor prognosis in melanoma patients
To assess how KDELR3 contributes to melanoma progression in patients, we utilized multiple melanoma patient databases, The Cancer Genome Atlas 30, 42 (TCGA) and Gene Expression Omnibus (GEO; GSE8401 29 , GSE19234 28 ). We found increased expression of the KDELR3-001 transcript, but interestingly not the alternate transcript, KDELR3-002, in late-stage (stage III and IV) metastatic melanoma patients compared to early-stage (stage I and II) melanoma patients (Fig.   6a), consistent with a role for KDELR3 in melanoma progression. Metastatic melanoma patients with KDELR3 copy number amplifications demonstrated reduced survival relative to patients without such alterations ( Supplementary Fig. 7b). We next assessed melanoma patient survival using KDELR3 expression as a prognostic marker (GEO 28,29 ). High KDELR3-expressing latestage metastatic melanomas show statistically significant association with poor patient outcome, whereas KDELR3 expression levels in early-stage primary tumor samples did not (Fig. 6b-c).
Taken together these data strongly support a role for KDELR3 in the advancement of late-stage metastatic melanoma and implicate KDELR3 as a bona fide MetDev gene.

KDELR3 and KDELR1 knockdown have opposing effects on lung colonization
As KDELR3 is the only member of the KDELR family to be identified as a MetDev gene by our analyses, including embryonic-specific upregulation and consistent upregulation in melanoma cell lines, we posited that different KDELR members have different functions in melanoma etiology/progression. To address this, we asked which pathways were co-regulated with KDELR1 expression and if these are the same or different relative to KDELR3-regulated pathways. GSEA analysis (FDR < 0.0001) of KDELR1 co-expressed genes in TCGA skin cutaneous melanoma patients (cBioPortal) 30, 31 revealed a strong enrichment of mitochondrial, metabolic and protein synthesis pathways (top 10 GO term enrichment, Fig. 7a), which differed from the most enriched pathways in KDELR3 co-expressed genes that consisted predominantly of ECM, trafficking and ERSR pathways (top 10 GO term enrichment, Fig. 7b). Moreover, knockdown of KDELR3/KDELR1 did not consistently alter expression of each other, suggesting that expression of these genes is not intrinsically linked (Supplementary Fig. 7c-d). These data intimate that KDELR1 and KDELR3 play different roles in melanoma progression. To test this, we compared the behavior of KDELR3 and KDELR1 knockdown cells using experimental metastasis assays.
Notably, in contrast to KDELR3 knockdown, which predictably diminished metastasis, KDELR1 knockdown actually increased metastasis, suggesting that KDELR1 contributes in a different way to melanoma etiology and can function as a metastasis suppressor ( Fig. 7c-d). Moreover, analysis of KDELR1 expression in skin cutaneous melanoma patients (TCGA) showed, unlike KDELR3, no significant difference between early-stage melanoma patients and late-stage metastatic melanoma patients (Fig. 7e). These data demonstrate that despite assumed redundancy between KDELR family members, KDELR3 and KDELR1 must have distinct roles, at least with respect to metastatic competence.

Discussion
Here we propose that metastatic cancer cells exploit innate pathways that are hardwired within their cellular lineage to ensure proper development. These pathways, quieted in the differentiated cell, can be reactivated under pathologic conditions. The genetic/epigenetic reactivation of pathways that allow embryonic melanocytes to migrate, invade and colonize would represent an efficient strategy for melanoma cells to successfully metastasize. Here we employed a GEM model to identify, at the transcriptome level, novel genes that are required during melanocyte development and find that these are enriched in genes that are specific for progression of late-stage disease. We functionally validated 4 out of 4 genes tested, demonstrating the value of our dataset and supporting our hypothesis. We anticipate that other genes that passed our filtering criteria will ultimately prove to be functionally relevant and deserving of further analysis in future studies.
We report a mechanistic analysis of our top hit and melanoblast gene, KDELR3, a member of the KDEL receptor family. KDELR3 has neither been previously associated with cutaneous melanoma metastasis nor investigated in depth in the literature. Differences between KDELRs have been cited in the literature but the main focus has been the role of KDELR1 19,20,[43][44][45] . All three KDELR family members have been shown to mediate retrograde transport of proteins containing a Cterminus KDEL-like motif 19 . KDELRs typically reside in the cis-Golgi; however, tagged KDELRs are known to localize in both the cis-and trans-Golgi, which is consistent with our results 46 . Upon interaction with KDEL-like motif-containing proteins, KDELRs facilitate transport from the Golgi Apparatus back to the ER via COPI vesicles 47 . When this system fails, KDEL-like motifcontaining proteins have been shown to be secreted out of the cell 19 . Our data demonstrating reduced BiP protein in stable KDELR3 knockdown cells suggest that BiP is a genuine substrate for KDELR3 retrograde trafficking, and that without KDELR3 expression melanoma cells are unable to maintain normal BiP levels. KDELRs appear to differ in the substrates that they preferentially transport, suggesting they have distinct roles within the cell 19 . How preferential substrate binding of KDELRs may affect cellular biology or disease etiology is still unknown.
Our study is the first to show that distinct KDELRs mediate dramatically different experimental metastasis phenotypes. We demonstrate that the embryonic melanoblast gene, KDELR3, is a metastasis enhancer in both mouse and human melanoma cells, whereas KDELR1 suppresses metastasis, despite having extensive homology and similar retrograde trafficking functions. Our data allow a new perspective when interpreting existing KDELR literature and present a dichotomy between KDELR3 and KDELR1 metastasis phenotypes that could be leveraged in future studies to understand how these retrograde trafficking receptors function in disease. Moreover, Trychta and colleagues have reported tissue-specific KDELR expression patterns in rats implying that different KDELRs may have lineage-specific roles 19 . This is the first study to document KDELR expression in melanocyte development, and a specific role for KDELR3.
The KDEL receptors are intrinsically linked to ER stress and proteostasis. KDELR retrograde trafficking substrates include protein chaperones, protein folding chaperones, protein folding enzymes, enzymes that target proteins for degradation and glycosylation enzymes 19 . Cumulatively, these protein substrates help maintain correct protein processing, and regulate cellular response to ER stress 20 . However, the role of ER stress response in tumor progression has been much debated 48 .
The success of proteasome inhibitors in the treatment of multiple myeloma patients 49 , as well as provocative data linking ER stress pathways to vemurafenib-resistant melanoma and immunotherapy sensitization, suggest UPR/ERAD biology could be harnessed for treating metastatic melanoma [50][51][52][53][54] . Our analysis implicates both UPR and degradation pathways of the ERSR as acting downstream of KDELR3. We show that KDELR3 expression is critical for adaptation of melanoma cells to ER stress and provide evidence that PERK-EIF2α expression and activation is regulated by changes in KDELR3 expression levels. Activation of the PERK-EIF2α pathway is known to result in translational attenuation, a cellular mechanism to alleviate ER load, causing translational rewiring of cells and affecting metastasis 15,17,48,[55][56][57] , which may contribute to KDELR3's metastatic role.
We demonstrate that KDELR3 is a regulator of glycosylated KAI1, a tetraspanin glycoprotein with a well-documented metastasis suppressor role in tumors, including melanoma 23,36,37,58,59 . KAI1 functions at the cell membrane to mediate interactions between extracellular and intercellular signaling, which is key to its metastatic suppressor function. KAI1 glycosylation leads to changes in its membrane organization and therefore its ability to mediate this extracellular/intercellular signaling 38,39,60 . However, no studies have linked specific KAI1 glycosylated forms with its metastasis suppressor function in vivo. Our work notes specific glycosylated forms of KAI1 that are subject to KDELR3 regulation and associated with metastatic function. Future work would benefit from determining how critical each of these forms are to KAI1's metastatic influence in vivo. Previously KDELR1 was shown to mediate signaling and transcriptional networks 43 , and at the protein level, in the relocation of lysosomes and modulation of autophagy 61 . However, KDELR3 was shown to be inactive in these processes. Here we link KDELR3 to post-translational regulation of protein, specifically through post-translational modification (glycosylation) and degradation of the metastasis suppressor, KAI1. Our data insinuate an interaction with gp78, implicating ERAD in this process. This biology may be informative for developing therapeutics for KDELR3-high metastatic melanoma patients.
We here identify an enrichment of ECM organization and trafficking genes within our MetDev cohort, consistent with a known role for these in metastasis [62][63][64] . Further analysis of these genes/pathways may prove a rich resource to uncover novel metastasis biology. We found that two such genes, KDELR3 and P4HA2 (a collagen prolyl 4-hydroxylase involved in ECM remodeling and associated with worse clinical outcome in melanoma patients 65 ), from our 4-gene functional validation screen are tightly co-expressed in four independent mouse models and in human melanoma patients. This raises the possibility that expression of some genes within our MetDev cohort may be coordinated and/or networked to realize the complex and dynamic phenotypes exhibited by melanocytic cells during development and metastasis. Uncovering common upstream regulators of co-regulated genes could prove a powerful approach to target metastatic melanoma as multiple pathways could be targeted simultaneously.
To our knowledge this is the first example in which the mouse melanoblast transcriptome has been exploited to generate a resource of novel melanoma metastasis genes. The success of this study supports the use of developmental models to uncover innate melanoma biology that may be at the root of melanoma's propensity to metastasize 2-9, 11-13, 66 . We anticipate that further exploration of KDELR3 and other now-uncovered embryonic genes/pathways will facilitate the development of more effective treatment strategies for patients with advanced melanoma, and perhaps other tumor types. The field would further benefit from elucidation of the specific melanoblast cell characteristics/cell states that in fact contribute to metastasis. In summary, this work describes the creation of a novel resource of putative MetDev genes, enriched in genes that have functional roles in melanoma metastasis that may prove to be useful targets for designing more effective approaches to the treatment of melanoma patients.

Mouse models of melanoma
Experimental metastasis studies were performed using a filtered, single-cell suspension in PBS.

Isolation of melanoblasts and melanocytes
FVB/N iDct-GFP dams were fed doxycycline-fortified chow for the entire duration of gestation until collection of E15.5, E17.5 and P1 pups. Doxycycline was injected intraperitoneally at 80 μg/g body weight 24 hr before collection of P7 pups. A single cell suspension was generated from embryos and skin of newborn pups. Multiple litters were used for each developmental sage, and embryos/pups from each stage were pooled to ensure adequate numbers of GFP + cells. The head was removed to prevent collection of GFP positive cells in the embryonic telencephalon, and melanocytes from the inner ear or from the retinal pigmented epithelium (RPE) were discarded.
Excess tissue was also removed. The spinal cord was kept intact as some melanoblasts still remain in the neural crest area. At E17.5, P1 and P7 stages, most melanocytes have reached the dermis, thus only the skin was collected from these developmental stages. Back skin was immersed in a shallow layer of 1X PBS and subcutaneous fat was scraped off until skin appeared translucent.
E15.5 was the youngest age assessed due to the necessity to capture sufficient cells for RNAsequencing.

Preparation of single cell suspensions
Tissue was minced and incubated for 30 min at 37°C in digestion media containing RPMI 1640 (Gibco Life Technologies) with 200 units/ml Liberase TL (Roche Applied Science). Up to 1g of tissue was digested per 5 ml digestion media. Tissue was processed using a Medimachine (BD Biosciences) and sterile medicon units (BD Biosciences). Cells were extracted using 1.5-2 ml of RFD solution (24 ml RPMI media, 6 ml FBS, 300 µl 5% DNase I) through a 20 ml syringe with 18-gauge needle. Collected cells were filtered through a 50-micron filter (BD Biosciences). This process was repeated until all the tissue was processed. Cells were spun at 1200 rpm, 4°C for 5 min and resuspended twice in a solution of 1% BSA in PBS and filtered through a 30-micron filter for sorting.

Fluorescence activated cell sorting (FACS)
Embryos of the same developmental age that were heterozygous for the TRE-H2B-GFP gene but lacked the Dct-rtTA gene were used as negative controls. Cell doublets were excluded from the analysis. Cells were sorted based on GFP expression and SSC-A. Based on these reference sorts, gates were set so that background cells represented less than 10% of sorted cells.

RNA Isolation and RNA sequencing
Cells were lysed in 10-fold TRIzol reagent (w/v), phases were separated by addition of 0.2X volume of chloroform, the aqueous phase was combined with an equal volume of 70% ethanol and applied to a RNeasy Micro column (Qiagen) and processed as per the manufacturer's instructions.

Analysis of MetDev genes in patient survival
Based on the RNASeq data for the samples E15, E17, P1 and P7, we used DESeq2 to find differentially expressed genes comparing E15, E17 vs. P1, P7. We selected 467 up regulated genes with q-value < 0.1 (based on glm model) and log2FoldChange > 1.5. We then used the GEO dataset GSE19234 to perform survival analysis using Cox proportional hazards model for each gene. We to selected 43 genes that were correlated patient overall survival with the patient survival with pvalue < 0.1 and HR >1. The Figure 1c-d showed the heatmaps of the gene expression (using zscore) for the 467 genes and 43 genes respectively. The sum of the total expression of the 43 genes forms the expression signature for prognosis prediction and the signature was tested on the new dataset GSE8401. Among the late-stage patients, the patients with high expression signature had significant poor survival compared to those with low expression (P=3.486e-5, Logrank test, Figure   1e) while for the early stage patients the two groups had no difference in survival (Figure 1e-f).

Gene filtration pipeline for functional analysis
From our 467 identified melanoblast genes we first filtered for only those genes whose P7 expression level was minimal (FPKM < 2) i.e. no functional P7 gene expression to our knowledge, reasoning these would denote genes that truly had a unique role in melanoblast development compared to differentiated melanocytes. Next, we validated these by identifying the genes that are the intersect of the 467 genes with the differentially expressed genes from microarray expression data derived from our iDct-GFP model (E17.5 vs P2 and P7) 21 . The microarray differential gene expression was identified using a linear regression model with contrast to compare embryonic versus postnatal stages and selected with a q-value < 0.1. The intersect yielded 233 genes. We acknowledge that the microarray data is not as thorough a representation of melanoblast/melanocyte development as our developmental cohort and therefore we may incur false negatives, we deemed this acceptable however to shorten our list for experimental validation.
Next, we filtered the list to 81 genes with > 2.75 log2 fold increase expression in melanoblasts vs melanocytes and P-value <0.0003. Finally, we reasoned that genes with the greatest expression at embryonic stages would likely be the most functionally relevant, so selected for the top 20 greatest mean embryonic expression. Of these 20 we noted that 7 genes (Kdelr3, P4ha2, Gulp1, Dab2, Lum, Aspn, Mfap5) were all associated with Extracellular Matrix (ECM) or trafficking. Of these we chose to test the 3 least studied genes in metastasis (Kdelr3, P4ha2, Gulp1) to uncover novel metastasis biology, and the 1 gene well studied in metastasis (Dab2).

Statistical analysis of KDELR3 expression in microarray data
Mouse developing melanoblasts (E17.5, n = 3) and differentiated melanocytes (P2, n = 3) were isolated and RNA extracted for microarray analysis as previously described 21  Welch's correction was used to compare the mean expression of KDELR3 between the two developmental stages. As two probes for KDELR3 on the Illumina beadchip showed high positive correlation (r = 0.987), the average KDELR3 expression was analyzed.

Analysis of The Cancer Genome Atlas (TCGA) skin cutaneous melanoma expression
All patient samples were collected between 0-14 days after disease classification (101 patients). t-statistic test 69 was applied to test the null hypothesis both for no difference in KDELR3 expression, or for KDELR1 expression level between early and late stage melanoma patients. A P-value of 0.05 or less was considered statistically significant.

Statistics and general methods
All sample sizes were determined based on preliminary studies and prior knowledge of expected variability within assays. For animal studies, age-matched (6-8 weeks) female ATHYMIC NCrnu/nu mice and C57BL/6 mice were randomly assigned to control and test groups. Blinding was used to quantify lung metastases counts. Where blinding was not used, data was analyzed using automated image analysis software when possible. All statistical tests used were deemed appropriate and met the assumptions required. Where necessary unequal variance was corrected for, or if no correction was used variation was assumed equal based on prior knowledge of the experimental assay. All cell lines used in this paper were identified correctly as per the

Code Availability
Upon acceptance of the manuscript custom code will be made publicly available and a full code availability statement will be included here.

Melanoma cell lines
Human

Immunohistochemistry (IHC) and Immunofluorescence (IF) staining
Formalin-fixed paraffin-embedded (FFPE) immunofluorescence of iDct-GFP mouse skin sections was performed using Heat Induced Epitope Retrieval (HIER) in Target

Exogenous expression studies
For exogenous over-expression of CD82 and KDELR3 genes the following expression plasmids

Mass spectrometry
Cell lysates were extracted 4 days post siRNA knockdown of KDELR3, AMFR or non-targeting control (siGENOME) using Dharmafect #1 transfection reagent. Cell lysates (250 μg each) were digested with trypsin using the filter-aided sample preparation (FASP) protocol as previously Samples analyzing the effect of KDELR3 or AMFR siRNA knockdown were quantitated using label-free methods. Dried peptides were fractionated by high pH reversed-phase spin columns (Thermo). The peptides from each fraction were lyophilized, and dried peptides were solubilized in 4% acetonitrile and 0.5% formic acid in water for mass spectrometry analysis. Each fraction of each sample was separated on a 75 µm x 15 cm, 2 µm Acclaim PepMap reverse phase column (Thermo) using an UltiMate 3000 RSLCnano HPLC (Thermo) at a flow rate of 300 nL/min followed by online analysis by tandem mass spectrometry using a Thermo Orbitrap Fusion mass spectrometer. Peptides were eluted into the mass spectrometer using a linear gradient from 96% mobile phase A (0.1% formic acid in water) to 35% mobile phase B (0.1% formic acid in acetonitrile) over 240 minutes. Parent full-scan mass spectra were collected in the Orbitrap mass analyzer set to acquire data at 120,000 FWHM resolution; ions were then isolated in the quadrupole mass filter, fragmented within the HCD cell (HCD normalized energy 32%, stepped ± 3%), and the product ions analyzed in the ion trap.
The mass spectrometry data were analyzed and either dimethyl labeling or label-free Reactions were analyzed by NU-PAGE Nu-PAGE and immunoblotted with the rabbit monoclonal antibody to CD82 (D7G6H) was used (Cell Signaling Technology Cat# 12439S).