Crosstalk with keratinocytes causes GNAQ oncogene specificity in melanoma

Different melanoma subtypes exhibit specific and non-overlapping sets of oncogene and tumor suppressor mutations, despite a common cell of origin in melanocytes. For example, activation of the Gαq/11 signaling pathway is a characteristic initiating event in primary melanomas that arise in the dermis, uveal tract, or central nervous system. It is rare in melanomas arising in the epidermis. The mechanism for this specificity is unknown. Here, we present evidence that in the mouse, crosstalk with the epidermal microenvironment actively impairs the survival of melanocytes expressing the GNAQQ209L oncogene. We found that GNAQQ209L, in combination with signaling from the interfollicular epidermis (IFE), stimulates dendrite extension, leads to actin cytoskeleton disorganization, inhibits proliferation, and promotes apoptosis in melanocytes. The effect was reversible and paracrine. In contrast, the epidermal environment increased the survival of wildtype and BrafV600E expressing melanocytes. Hence, our studies reveal the flip side of Gαq/11 signaling, which was hitherto unsuspected. In the future, the identification of the epidermal signals that restrain the GNAQQ209L oncogene could suggest novel therapies for GNAQ and GNA11 mutant melanomas.


INTRODUCTION
To identify the pathways that change in response to GNAQ Q209L expression in IFE melanocytes, 259 we performed RNA sequencing (RNA-seq) immediately after sorting WT and GNAQ Q209L 260 melanocytes from tail IFE (n=6 mice and n=3 libraries per genotype). There were 14,461 genes 261 with an average FPKM > 0.1 in WT and/or GNAQ Q209L melanocytes (Supplementary Table 1). 262 Seven genes previously shown to be directly involved in pigment production and melanosome 263 biology were on the list of the top 20 most highly expressed genes in WT melanocytes, 264 validating our melanocyte isolation protocol (Supplementary Table 2). These genes were Pmel, 265 Dct, Tyrp1, Mlana, Cd63, Slc24a5, and Gpnmb. We analyzed differential gene expression using 266 DEseq2. 1,745 genes were found to be differentially expressed (DE), of which 729 genes were 267 down-regulated and 1016 genes were up-regulated in GNAQ Q209L melanocytes (q-value < 0.05) 268 (Supplementary Tables 3 and 4).  Table 5). Overlapping genes supported a second signature for the terms: axons, axon guidance 277 and nervous system development (Supplementary Table 6 of DE genes (with no cut-off), we also found that two Rho GEFs that interact with Gα q/11 , Trio 283 and Kalirin (Kalrn), were up-regulated (LFC 0.5 and 1.8). RhoGEFs catalyze GDP to GTP 284 exchange on small Rho guanine nucleotide-binding proteins, regulating the actin cytoskeleton. 285 We also noticed that several genes that play an essential role in lamellipodia formation 21 286 including RhoB, RhoC, Cdc42, and Rock were all up-regulated by LFC > 1.0. Increased adhesion 287 and cell-ECM interactions could underlie the changes in cell shape and increased cellular area in 288 GNAQ Q209L melanocytes co-cultured with IFE (Figure 3). 289 290 To examine the actin cytoskeleton in GNAQ Q209L and WT melanocytes, we measured the overall 291 level of actin alignment, termed cell fibrousness, using an in vitro quantitative approach 292 previously described by Haage et al 22 . First, WT and GNAQ Q209L tdTomato-positive 293 melanocytes co-cultured with IFE were stained for f-actin using phalloidin ( Figure 6B). Then, 294 confocal z-projections of individuals melanocytes were obtained and processed using 295 quantitative image analysis to determine cell fibrousness. GNAQ Q209L expressing melanocytes 296 exhibited less cell fibrousness than WT melanocytes (p=1.8x10 -12 , Figure 6C) Using the RNAseq data, Gene set enrichment analysis (GSEA) revealed that GNAQ Q209L IFE 328 melanocytes express a pattern of genes related to cellular stress and apoptosis. GSEA showed 329 that GNAQ Q209L melanocytes are enriched in the hallmarks for "apoptosis", "p53 pathway", and 330 "hypoxia" (Figure 6D). Among the top-ranked DE genes was Stanniocalcin (Stc1) which 331 encodes a glycoprotein hormone involved in calcium/phosphate homeostasis. Stc1 was up-332 regulated by a LFC of 7.1 in GNAQ Q209L melanocytes. Significant up-regulation of Stc1 has been 333 reported in tumors under hypoxic or oxidative stress 27-31 . Furthermore, Stc1 expression down-334 regulates pro-survival ERK1/2 signaling and reduces the survival of MEFs under conditions of 335 oxidative stress 28 . Cdkn2a and Ccnd1 were up-regulated in GNAQ Q209L expressing melanocytes 336 by a LFC of 5 and 1.7, respectively, similar to the intracellular response to oxidative stress 337 previously described in melanocytes in the disease, vitiligo, which results in the loss of 338 melanocytes from the epidermis 32-34 . Moreover, the S100A8, S100A9 and BNIP3 genes, involved 339 in cell death and autophagy via ROS-mediated crosstalk between mitochondria and lysosomes, 340 were all up-regulated in GNAQ Q209L melanocytes by a LFC of 0.7 to 1.7 35 . 341

342
In terms of apoptosis, S100b, a negative regulator of p53 in melanocytes, was down-regulated by 343 (LFC -2.6), while Annexin V (Anxa5), a marker of early apoptosis, was up-regulated (LFC 0.6). 344 We also found up-regulation of the pro-apoptotic Bok (LFC 1.3), an effector of mitochondrial 345 outer membrane permeabilization (MOMP) that triggers membrane permeabilization and 346 apoptosis in the absence of Bax and Bak 36 . Cystatin C (Cst3), a major inhibitor of cathepsins 37 , 347 was up-regulated (LFC 2.3) in GNAQ Q209L melanocytes. In response to DNA damage, activated 348 p53 can induce Cystatin C expression by binding to regulatory sequences in the first Cst3 intron. 349 Hence, we hypothesized that cellular stress was stimulating increased apoptosis of GNAQ Q209L -350 expressing melanocytes. To assess the amount of cell death in vivo, we dissociated IFE from 351 mouse tails and stained the cells for Annexin-V, which labels cells in early apoptosis. Using 352 FACS, we quantified the percentage of tdTomato and Annexin-V double-positive cells ( Figure  353   6E). There was a greater percentage of Annexin-V positive melanocytes from GNAQ Q209L tail 354 epidermis (4.2% versus 1.2%). This supported the GSEA data that suggests that the loss of 355 GNAQ Q209L expressing melanocytes from the IFE involves cellular stress and apoptosis. 356 357

Low frequency of GNAQ and GNA11 oncogenic mutations in human cutaneous melanoma 358
Our results suggested that oncogenic mutations in GNAQ and its closely related homolog, 359 GNA11, are rarely found in cutaneous melanoma because the epidermal microenvironment 360 causes Gα q/11 signaling to inhibit melanocytes. To determine the frequency of these mutations, 361 we searched the COSMIC database to identify all cutaneous melanoma cases with a GNAQ or 362 cases carried one or more mutations in genes also found in uveal 38, 39 or mucosal 40, 41 melanoma, 368 but not in cutaneous melanoma: BAP1, SF3B1 or EIF1AX and hence were similarly suspect as 369 misclassified. One case was from a cell line and in another case, the frequency of an odd 370 GNAQ Q209H mutation in the tumor was less than 5%. This left eight cases that we think could 371 have arisen in the epidermis, for an incidence of 0.15% (4/2753) for GNAQ and 0.17% (4/2295) 372 for GNA11. Details on two of these cases can be found in Supplementary Table 9. 373 374

Mutations in PLCB4 in human cutaneous melanoma 375
Phospholipase C beta 4 (PLCB4) encodes a protein that binds to Gα q and Gα 11 and serves as the 376 primary effector of signaling to downstream components for q class heterotrimeric G proteins. 377  Table 9). The 407 14 | P a g e NRAS Q61 and BRAF V600 oncogenic mutations occur in an exact mutually exclusive pattern in the 408 TCGA-SKCM dataset (116 and 205 cases, respectively). Therefore, we began by comparing the 409 frequency of PLCB4, PTEN, TP53 and NF1 mutations in these subsets. PLCB4 was mutant in 410 24% of the NRAS Q61 mutant cases, 17% of the BRAF V600 mutant cases and 23% of the remaining 411 cases (hereafter referred to as 'other'). TP53 and NF1 were also less frequent in BRAF V600 mutant 412 cases and more frequent in the NRAS Q61 and other cases ( Figure 7A). NF1 mutations were 413 highly enriched in the 'other' subset, as expected from the literature 8, 40 . PTEN mutations, in 414 contrast, were more likely to be found in the BRAF V600 mutant subset compared to the NRAS Q61 415 mutant subset or other cases. We calculated the average total number of all gene associated small 416 somatic mutations in each of the three subsets and found that the 'other' subset had 1.7-fold 417 greater mutation load than the NRAS Q61 mutant subset. Hence, TP53 and PLCB4 mutations in the 418 NRAS Q61 subset could be enriched when considering the overall mutation burden. Lastly, 32% of 419 the PLCB4 mutant cases had more than one mutation in PLCB4 (i.e. potential biallelic mutation), 420 compared to 38% for NF1, 14% for TP53 and 4% for PTEN ( Figure 7B). 421

422
We then identified all cases that had a mutation in more than one of these four genes of interest 423 (PTEN, NF1, PLCB4 and TP53) and counted the number of times a double hit for any of the two 424 gene combinations was found ( Figure 7C). In agreement with the frequencies in Figure 8A, 425 PTEN mutations were less often found with other members of this group, particularly NF1. 426 PLCB4 mutations were most often seen with TP53 and NF1, which were found less frequently 427 with each other than with PLCB4. genes in each bin. The mutation rate of BRAF V600 and NRAS Q61 was highest in the bin with the 436 smallest mutation load. As mutation load increased, the mutation rate of BRAF V600 and NRAS Q61 437 15 | P a g e consistently decreased. PTEN also generally followed this pattern. In the 2000 to 3000 mutation 438 load/case range, the mutation rate of BRAF V600 , NRAS Q61 and PTEN was low, while NF1 and 439 PLCB4 reached their maximum. The rate of NF1 and PLCB4 mutations declined in the 3000+ 440 mutation load/case bin. TP53 showed no particular pattern across the bins. 441

442
We then organized all of the current TCGA-SKCM cases into Figure 8, reflecting on how the 443 overall mutation load relates to mutations in these six genes of interest. At the lowest mutation 444 load (average of 35 total mutations per case), we note that there were five cases of GNAQ or 445 GNA11 oncogenic mutations plus a mutation in either BAP1 or SF3B1. We suspect that these 446 cases could be tumors from uveal or mucosal melanomas, although they were not described as 447 such. Next, there were 74 cases with an average of 431 mutations per case, with no mutations in 448 any of the genes of interest described here. Cases with just BRAF V600 or NRAS Q61 closely 449 followed, with an average of 453 and 587 mutations per case, respectively. As the mutation load 450 continued to increase, hits in multiple genes of interest accumulated. First, there were cases with 451 just one tumor suppressor mutation plus BRAF V600 or NRAS Q61 (at 900 mutations/case). The 452 tumor suppressor mutation was more likely to be PTEN in the BRAF V600 mutant subset. The overall pattern is that gain of function drivers are more frequently found when mutation load 466 is lower and loss of function drivers predominate when mutation load is higher. Perhaps when 467 mutagen exposure is low, a stronger acting BRAF or NRAS mutation can drive the melanoma 468 because an inactivating hit to the same allele is less likely to occur. However, when mutagen 469 exposure is high, it alternatively allows for the accumulation of loss of function mutations in 470 multiple genes, which has a positive synergistic effect when they occur in tumor suppressors. Previous studies have shown that stable transfection of GNAQ Q209L induced anchorage-495 independent growth in soft agar of hTERT/CDK4 R24C /p53 DD mouse melanocytes in a TPA-496 independent manner, a feature associated with cellular transformation 3, 51 . In addition, 497 GNAQ Q209L signaling produced an abnormal cell morphology with irregular contours. A similar 498 observation was described in GNAQ Q209L -transformed human melanocytes and melanocytes of a 499 Gnaq Q209L zebrafish model 52 . Our study complements these previous results to show that 500 melanocytes directly taken from the mouse IFE and grown in primary culture without any other 501 cues survive better if they express GNAQ Q209L . We also found that Braf V600E had a stronger effect 502 than GNAQ Q209L in this same situation, although the promoters driving the expression of these 503 two oncogenes were not the same. In addition, both WT and GNAQ Q209L expressing melanocytes 504 survived better in culture when grown on MEFs. Therefore, the attrition of GNAQ Q209L 505 expressing IFE melanocytes from the mouse tail epidermis is not due to an inherent characteristic 506 18 | P a g e in IFE melanocytes because the phenotype is reversible by changing the microenvironment. This 507 is exciting because it suggests that there might be a way to shut the GNAQ Q209L oncogene off 508 through some external cues. 509

510
The IFE microenvironment inhibits melanocytes expressing GNAQ Q209L 511 It is well known that keratinocytes, which make up the vast majority of the epidermis, secrete 512 various growth factors and cytokines in a paracrine manner that regulate growth, survival, 513 adhesion, migration, and differentiation of melanocytes 53 . In WT melanocytes, we found that the 514 presence of dissociated IFE in co-cultures provided a significant boost. In contrast, co-culturing 515 GNAQ Q209L expressing melanocytes with IFE reduced their survival capacity and inhibited cell 516 division. Moreover, the effects of the IFE on melanocyte survival did not require direct cell-cell 517 contact and could be replicated in a transwell culture system. Our studies of IFE melanocytes 518 suggests that GNAQ Q209L expression causes increased cell adhesion, a disorganized actin 519 cytoskeleton and the promotion of long dendrite extensions, which frequently break, possibly 520 related to changes in semaphorin signaling in pathways shared with neurons. We conclude that 521 the microenvironment surrounding melanocytes is the main factor controlling whether 522 GNAQ Q209L signaling is oncogenic or inhibits growth. IFE melanocytes have the innate capacity 523 to be transformed by GNAQ Q209L , but whether or not they are permitted to do so is dictated by 524 the microenvironment in which they reside. substitutions. In addition, 32% of PLCB4 mutant cases have more than one mutation in the gene 572 (potential biallelic mutation), which was similar to NF1 at 38% of cases. The pattern of PLCB4 573 mutation across the set was not random; instead, PLCB4 mutations were less frequent in 574 BRAF V600 mutant cases and more frequent in cases with a higher overall mutation load, alongside 575 NF1 and TP53. Our finding that GΝΑQ Q209L signaling restrains melanocyte growth in the 576 epidermis provides an explanation for the frequent occurrence of PLCB4 mutations in cutaneous 577 melanoma. It will be of interest to study the relationship between PLCB4 loss and frequency and 578 location of metastases, since metastasis involves melanoma cells leaving the epidermal 579 microenvironment and interacting with mesodermal fibroblasts and stromas. 580 581 In conclusion, we have shown that the GNAQ Q209L oncogene switches from promoting to 582 inhibiting melanocyte growth due to paracrine crosstalk with the IFE microenvironment. Our 583 studies revealed evidence for multiple potential mechanisms, including oxidative stress, 584 apoptosis, inhibition of cell division, changes to cell adhesion and cell fragmentation. Hence, the 585 Gα q/11 signaling pathway is not only important for uveal melanoma; it also has implications for 586 both vitiligo and cutaneous melanoma. The critical paracrine signal(s) remain to be determined 587 through further biochemical and cell culture studies. It is also possible that a better understanding 588 of the ability of GNAQ Q209L to switch from acting as an oncogene to an inhibitor of melanocyte 589 growth could lead to new therapies for uveal melanoma, which lacks effective treatment options 590 for metastatic disease. culture with MEFs, 6000 to 8000 FACS sorted melanocytes were plated into 96-well plates 652 23 | P a g e previously seeded with MEFs that had formed a confluent monolayer. For direct co-culture 653 experiments with IFE, the IFE was dissociated and plated at 100,000 cells per well (i.e. with no 654 FACS sorting for the melanocytes first). For transwell co-culture with IFE, 6000 to 8000 FACS 655 sorted melanocytes were plated in the well suspended above 100,000 dissociated IFE cells plated 656 below. All cultures were incubated at 37 o C in 5% CO 2 . Media was changed every third day by 657 exchanging 1/3 of the existing volume. Images were taken at 5x, 10x, and 20x magnification. MatLab. Images were first blurred with a Gaussian filter, and then an edge detection algorithm 673 was applied to identify cell borders. The resultant binary images were refined through successive 674 dilations and erosions to yield the final cell contour. These contours were used to measure cell 675 area, and circularity (4π*area/perimeter 2 ), protrusions. The number and length of protrusions 676 were quantified automatically in MatLab. First, cell contours were identified as outlined above. 677 Next, we obtained the contour coordinates of the convex hull of the binary image representing 678 the cell area. At each point along the cell contour, we computed the minimum distance between 679 the convex hull and the actual cell contour. Based on these distances, minima corresponding to 680 protrusions could be identified. To be counted as protrusions, minima had to be at least 10 pixels 681 peaks, the width, height, and aspect ratio of protrusions could be computed. 683 684 Actin fiber organization 685 First, we identified cell contours as described above. Next, the cell was subdivided into 32x32 686 pixel windows overlapping by 50%. We then computed the two-dimensional Fourier transform 687 of each window. If a window contains no fibers, the Fourier transform will be a central, diffuse 688 point of bright pixels. However, if a window contains aligned fibers, the Fourier transform will 689 consist of an elongated accumulation of bright pixels at a 90 degree angle to the original fibers. 690 Based on the aspect ratio and orientation of the Fourier transform, we determined fibrousness 691 and fiber orientation in a given window. The data for individual windows could then be 692 compared across the entire cell to estimate the cell fibrousness, defined here as the percentage of 693 cell area (% of windows) with aspect ratio greater than a cut-off value. 694

RNA-seq data and bioinformatics analysis 696
Total RNA from sorted cells was extracted using Trizol (Life Technologies) following the 697 manufacturer's protocol. For each library preparation, sorted melanocytes from two mice of the 698 same litter were pooled to increase the number of cells for RNA extraction. We generated 3 Wt 699 libraries and 3 GNAQ libraries from a total of 12 mice. Sample quality control was performed 700 using the Agilent 2100 Bioanalyzer. Qualifying samples (RNA Integrity Number ≥ 9) were then 701 prepped following the standard protocol for the NEB next Ultra ii Stranded mRNA (New 702