# Dynamics of nevus development implicate cell cooperation in the growth arrest of transformed melanocytes

1. Center for Complex Biological Systems, University of California, Irvine, United States
2. Department of Developmental and Cell Biology, University of California, Irvine, United States
3. Department of Dermatology, University of California, Irvine, United States
4. Beckman Laser Institute, University of California, Irvine, United States
5. Department of Mathematics, University of California, Irvine, United States
6. Department of Biological Chemistry, University of California, Irvine, United States
Research Article

## Abstract

Mutational activation of the BRAF proto-oncogene in melanocytes reliably produces benign nevi (pigmented ‘moles’), yet the same change is the most common driver mutation in melanoma. The reason nevi stop growing, and do not progress to melanoma, is widely attributed to a cell-autonomous process of ‘oncogene-induced senescence’. Using a mouse model of Braf-driven nevus formation, analyzing both proliferative dynamics and single-cell gene expression, we found no evidence that nevus cells are senescent, either compared with other skin cells, or other melanocytes. We also found that nevus size distributions could not be fit by any simple cell-autonomous model of growth arrest, yet were easily fit by models based on collective cell behavior, for example in which arresting cells release an arrest-promoting factor. We suggest that nevus growth arrest is more likely related to the cell interactions that mediate size control in normal tissues, than to any cell-autonomous, ‘oncogene-induced’ program of senescence.

## eLife digest

Melanocytes are pigment-producing cells found throughout the skin. Mutations that activate a gene called BRAF cause these cells to divide and produce melanocytic nevi, also known as “moles”. These mutations are oncogenic, meaning they can cause cancer. Indeed, BRAF is the most commonly mutated gene in melanoma, a deadly skin cancer that arises from melanocytes. Yet, moles hardly ever progress to melanoma.

A proposed explanation for this behavior is that, once activated, BRAF initiates a process called “oncogene-induced senescence” in each melanocyte. This process, likened to premature aging, is thought to be what causes cells in a mole to quit dividing. Although this hypothesis is widely accepted, it has proved difficult to test directly.

To investigate this notion, Ruiz-Vega et al. studied mice with hundreds of moles created by the same BRAF mutation found in human moles. Analyzing the activity of genes in individual cells revealed that nevus melanocytes that have stopped growing are no more senescent than other skin cells, including non-mole melanocytes.

Ruiz-Vega et al. then analyzed the sizes at which moles stopped growing, estimating the number of cells in each mole. The data were then compared with the results of a simulation and mathematical modeling. This revealed that any model based on the idea of cells independently shutting down after a number of random events could not reproduce the distribution of mole sizes that had been experimentally observed. On the other hand, models based on melanocytes acting collectively to shut down each other’s growth fit the observed data much better.

These findings suggest that moles do not stop growing as a direct result of the activation of BRAF, but because they sense and respond to their own overgrowth. The same kind of collective sensing is observed in normal tissues that maintain a constant size. Discovering that melanocytes do this not only sheds light on why moles stop growing, it could also help researchers devise new ways to prevent melanomas from forming.

## Introduction

Activating BRAF mutations (e.g. BRAFV600E) are the most common oncogenic mutations in melanoma, seen in about 66% of cases (Davies et al., 2002). Curiously, the same mutation is found in 89% of melanocytic nevi (Pollock et al., 2003)—the benign, pigmented ‘moles’ found on the skin of most individuals. In animal studies, melanocyte-specific expression of BRAFV600E efficiently produces nevi, but only very rarely melanoma (Dankort et al., 2009; Dhomen et al., 2009; Patton et al., 2005). The widely-accepted explanation is that transformed melanocytes undergo oncogene-induced senescence (OIS), arresting proliferation before additional oncogenic events can occur (e.g. Bennett, 2003; Huang et al., 2017; Kaplon et al., 2014; Michaloglou et al., 2005).

Nevus melanocytes are indeed growth-arrested, but the assumption that OIS is the cause remains untested, in part because of a lack of criteria to rigorously define OIS in vivo (Damsky and Bosenberg, 2017). Initially studied as a consequence of forced expression of oncogenes in cell cultures (Serrano et al., 1997), OIS has come to be seen as a distinctive cellular stress response characterized by a phenotype of growth arrest, morphological and metabolic changes, chromatin alterations, and secretion of growth factors, chemokines, cytokines and proteases (Campisi and d'Adda di Fagagna, 2007; Gorgoulis et al., 2019; Ito et al., 2017; Kuilman et al., 2010).

Given an abundance of ‘hallmarks’ of senescence, one might think that recognizing this cell state in vivo should be straightforward. Yet no single hallmark distinguishes senescence from other growth-arrested cell states. Phenotypes once thought to be ‘gold standards’, such as expression of lysosomal beta-galactosidase, cyclin-dependent kinase inhibitors, or p53, commonly mark only subsets of senescent cells (Wiley et al., 2017), as well as non-senescent cells (Tran et al., 2012). Moreover, observations of supposedly senescent cells resuming proliferation (e.g. Beauséjour et al., 2003), imply that permanent cell cycle exit cannot be used as a distinguishing feature. In vivo senescence, as a result, is currently somewhat of a Gestalt diagnosis, that?is assessed by a collection of traits, no subset of which is necessary or sufficient. Yet there is no clear consensus on which traits are best to assess, and recent meta-analyses of gene expression suggest that some of the most commonly assessed features are not ‘core’ to senescence at all (Hernandez-Segura et al., 2017).

The reason it is important to clarify how BRAF-transformed nevus melanocytes stop growing is that it shapes how we think about the origins of melanoma. OIS is usually portrayed in cell-intrinsic terms: oncogene expression within a transformed cell produces a stress within that cell, which triggers it to senesce. Even those who acknowledge a possible role for paracrine signals (Acosta et al., 2013; Elzi et al., 2012; Ito et al., 2017; Wajapeyee et al., 2008) still portray the process as something initiated and orchestrated by cell-autonomous responses to oncogenes. This naturally leads to an approach to melanoma prevention and treatment that focuses on understanding how oncogenes derange intracellular processes; how those derangements elicit stress responses; and what might enable cancer cells to circumvent those responses (e.g. Bennett, 2003; Damsky and Bosenberg, 2017; Vredeveld et al., 2012; Yu et al., 2018). In contrast, as we argue below, it is possible that the growth arrest displayed by nevus melanocytes has little to do with oncogene-induced stress, and may have more to do with networks of cell–cell communication that are characteristic of melanocytes, independent of whether they are transformed. In this case, the most effective path to understanding how to prevent or treat melanoma could be to better elucidate the normal physiology of melanocytes in their environment.

Here, we investigate the details of nevus growth arrest in a model in which melanocyte-specific Braf activation generates hundreds of nevi on the skin of mice (Dankort et al., 2009). By examining both single-cell transcriptomes and the dynamics of growth arrest in nevus-associated melanocytes, we make two key observations: First, patterns of gene expression in arrested nevus melanocytes fail to identify them as any more senescent than other skin cells or normal melanocytes, arguing against a primary role for any form of senescence in their arrest. Second, the timing and statistics of nevus formation effectively argue against any relatively simple cell-autonomous process as being the cause of growth arrest. Ultimately, we propose a model in which arrest is driven not by oncogene stress, but by feedback mechanisms similar to those commonly involved in normal tissue homeostasis.

## Results

### Dynamics of nevus growth

Characterizing the dynamics of nevus growth and arrest requires observing nevi that started growing at known times. We took advantage of a mouse model in which Cre-mediated recombination introduces the activating V600E mutation into the endogenous Braf locus. When crossed onto a background carrying a Tyr-CreER transgene, the mice acquire the BrafV600E mutation only in cells of the melanocytic lineage, and only after Cre activation by 4-hydroxytamoxifen (4-OHT), applied either systemically or through painting on the skin.

As shown previously (Dankort et al., 2009), 4-OHT treatment of these mice leads to development of numerous pigmented nevi. Visualization of nevi is hindered, however, by the strong pigmentation in hair follicles which, except at microscopic resolution, can be difficult to distinguish from nevi. One way to circumvent this difficulty is to observe nevi only during the telogen phase of the hair cycle, when follicle-associated pigment is not present (conveniently, synchronization of hair cycles may be maintained on a large patch of skin through depilation).

As shown in Figure 1, in mice whose back skin was treated with topical 4-OHT at postnatal day 2 (P2), P3 and P4, nevi were apparent macroscopically at telogen (P50; Figure 1A). Live?imaging, using multi-photon microscopy (MPM; Saager et al., 2015), revealed that, like human nevi, mouse nevi consist of scattered nests of pigment-containing cells (Figure 1B). Nevi could also be visualized post-mortem, using a dissecting microscope, on the undersurface of pieces of telogen-stage back skin (Figure 1C).

Figure 1 with 1 supplement see all

An alternate approach to visualization that did not require hair synchronization was to generate nevi by painting 4-OHT on glabrous (hairless) skin, such as the ventral surface of the paw, permitting tracking of individual nevi on a daily basis. As shown in Figure 1D, when forepaws were treated with 4-OHT from P2 through P4, tiny nevi could be detected as early as P6. Serial observation indicated that most nevi reach a maximum size somewhere between P16 and P21 (Figure 1D ,?Figure 1—figure supplement 1A). This suggests that, in the mouse, BrafV600E-transformed melanocytes arrest within 2–3 weeks. To confirm this, we used BrdU labeling to monitor DNA synthesis. Because melanin readily obscures immunohistochemical signals, these experiments were done in an albino (unpigmented) genetic background, using premelanosome protein (Pmel) staining to identify melanocytes. As shown in Figure 1F, albino mice generate nests similar to those seen in pigmented mice. In such animals, BrdU readily incorporated into hair follicle melanocytes (Figure 1E, Figure 1—figure supplement 1B), whereas by p21 nests within nevi were uniformly negative for BrdU, implying growth arrest (Figure 1F, Figure 1—figure supplement 1C).

The conclusion that Braf-induced nevi are already growth-arrested by P21 agrees with the reports of others (Damsky and Bosenberg, 2017), and is lent further support by time course measurements of nest size by MPM (Figure 1C), which show that nest size distributions change insignificantly between P21 and P50 (Figure 4—figure supplement 1A-B).

### Do nevi undergo OIS?

As discussed above, senescence is usually accompanied by distinctive gene expression. Various gene expression ‘signatures’ have been developed to help investigators identify senescent cells and distinguish them from cells that have become growth-arrested by other processes. We considered several of these (Source data 1):

1. A set of genes encoding the most commonly considered ‘hallmarks’ of senescence, that?is p53, Rb, lysosomal beta-galactosidase, H2AX, and three cyclin-dependent-kinase inhibitors (‘Classical’, seven genes [Campisi and d'Adda di Fagagna, 2007; Collado and Serrano, 2006]);

2. A gene signature that distinguishes cultured human fibroblasts growth-arrested by BRAF-transformation from quiescent fibroblasts (‘Kuilman’, 21 genes [Kuilman et al., 2008])

3. Results of a meta-analysis (Hernandez-Segura et al., 2017) of publications on fibroblasts, melanocytes, and astrocytes, comparing senescence (induced by multiple different stresses) with quiescence, yielding ‘universal’ signatures of genes that are up- and downregulated specifically in senescence (‘Universal Up’, 31 genes, and ‘Universal Down’, 23 genes) as well as signatures of genes specifically up- or down-regulated by senescence induced in melanocytes (‘Melanocyte Up’, 397 genes and ‘Melanocyte Down’, 135 genes).

4. The most statistically significant genes in a recent meta‐analysis (Chatsirisupachai et al., 2019) of 20 replicative senescence microarray datasets from the Gene Expression Omnibus (‘Chatsirisupachai Up’, 237 genes and ‘Chatsirisupachai Down’, 244 genes).

5. A list of genes characteristic of the ‘Senescence-Associated Secretory Phenotype’ (‘SASP’, 81 genes), compiled from 38 literature references (for citations see Source data 1).

To determine whether any of these proposed signatures fits nevus melanocytes, we performed single-cell RNA?sequencing on dissociated cells from the back skin of nevus-bearing mice at both P30 and P50 (i.e. after nevi have stopped growing), using wildtype skin as a control. Using known cell-type marker genes (Figure 2—figure supplement 1A-B), we identified 14 different cell types in the skin, including melanocytes (Figure 2A). Unsupervised clustering further sub-divided the melanocytes into four groups (Figure 2B): Two of them, Mel 0 and Mel 1, were composed of cells found only in nevus-bearing, and not wildtype, skin (Figure 2C); they are highly similar in gene expression, primarily differing in having a slightly lower level of pigment gene expression in Mel 1 versus Mel 0 (Figure 2D). We identify them as the ‘nevus melanocytes’, because they are seen only when nevi are present, and are by far the predominant melanocyte population in such animals.

Figure 2 with 1 supplement see all

Mel 2 cells express the lowest levels of pigmentation genes (Figure 2D), and are seen in both genotypes (Figure 2C) at all stages (although expanded in number in nevus-bearing animals (Figure 2F)). Their pattern of gene expression bears a strong resemblance to one recently published for melanocyte stem cells isolated from telogen-stage hair follicles (Zhang et al., 2020). In particular, they express Cd34, which has been proposed to be a marker for bulge-associated melanocyte stem cells (Joshi et al., 2019).

Finally, cells of cluster Mel 3, which express the highest levels of pigment genes (Figure 2D), are found in both mutant and wildtype mice, but only at the P30 time point (Figure 2E–F). We thus identify them as mature hair follicle melanocytes, as such cells are present exclusively during anagen phase of the hair cycle (P30), and disappear during telogen (P50).

Because gene signatures are based on the idea of up- and downregulation of expression relative to some baseline state, to test whether nevus melanocytes fit a known signature it is necessary to have comparison transcriptomes. We made two types of comparisons: nevus melanocytes versus every other cell type in the skin (which, with the possible exception of mature keratinocytes, we would not expect to be senescent); and the four melanocyte subclusters (two of which are nevus-associated and two of which are not) versus each other. In each case we computed average expression for each gene in every cell type or cluster, together with a standard error of the mean as a measure of dispersion. Expression values were then normalized to average expression across all of the cells being compared (i.e. all skin cells, or all melanocytes, depending on which comparison was being done) and log2-transformed, so that positive values signify upregulation (relative to the average for that gene), and negative downregulation. Gene expression was then visualized using heat maps (Figure 3A and Figure 3—figure supplements 12, with positive values in blue and negative in red). Because gene expression levels inferred from single-cell RNA sequencing tend to be noisy, particularly for genes with low expression, we ranked all genes by their minimum level of noise (i.e. normalized standard error of the mean in the least noisy cell type), and used this value (‘n-SEM’, which is also presented graphically as a bar to the right of each heatmap) to sort gene lists, so that maps vary from most to least reliable as one goes from top to bottom.

Figure 3 with 3 supplements see all

Figure 3A shows the results for the ‘Classical’ and ‘Universal Up’ signatures (heat maps for the other signatures are shown in Figure 3—figure supplement 1). Here we see no strong enrichment of blue over red signals in nevus melanocytes, nor in most other cells. When compared with whole skin, using the ‘Classical’ signature, only Cdkn2a stands out as strongly upregulated in nevus melanocytes, but it is similarly upregulated in skin fibroblasts (it also has the noisiest data among genes in the signature). With the ‘Universal Up’ signature, more genes are downregulated than upregulated in nevus melanocytes. To quantify such impressions, we summed the log2-transformed data in each column in every heat map, producing the bar graphs in Figure 3B. We reasoned that summation of log-transformed data would emphasize consistent trends in the data while suppressing effects of noise (random positive and negative variation would tend to cancel out). The results suggest that skin fibroblasts better fit the ‘Classical’ senescent signature than any skin cell type, including nevus melanocytes, or indeed melanocytes of any cluster. As a control—to demonstrate the ability of this approach to correctly associate cell types with gene signatures— we analyzed the same data using a signature of cell proliferation, ‘meta-PCNA’, that represents 129 human genes most positively correlated with proliferation marker PCNA in a compendium of normal tissues (Venet et al., 2011). As shown in Figure 3B (also see Figure 3—figure supplement 1), this signature (122 genes of which had unambiguous mouse orthologs; Source data 1) identified two keratinocyte populations (‘IFE-cycling’ and ‘Outer Bulge 1’) as highly proliferative (in agreement with Joost et al., 2020), and mature (postmitotic) keratinocytes as non-proliferative. Importantly, it also correctly identified nevus melanocytes as non-proliferative—and other melanocytes as proliferative—in agreement with Figure 1E–F. Interestingly, the relatively high expression of proliferation-associated genes in non-nevus, hair follicle melanocytes (Mel 3) when compared with nevus melanocytes, was consistent between BrafV600E-expressing and control mice (Figure 3—figure supplement 2), suggesting that tissue context plays a role in whether BrafV600E-expressing cells even arrest growth.

Figure 3C and Figure 3—figure supplement 1 extend this analysis to the remaining eight potential signatures of senescence (five consisting of genes that are upregulated; three of genes that are downregulated). In five of the eight cases, nevus melanocytes rank as less senescent than the average skin cell; in two of the cases nevus melanocytes are about average. In only one case (‘Chatsirisupachai Down’) does nevus melanocyte gene expression go in the predicted direction for senescence. However, the Chatsirisupachai signatures had not been curated to remove genes associated simply with proliferation/quiescence (Chatsirisupachai et al., 2019), and inspection of the ‘Chatsirisupachai Down’ gene list shows that 61 of its 250 genes are shared with the 129-gene meta-PCNA signature; that?is it is more likely a signature of proliferation than ‘non-senescence’ (note the strong similarity between the ‘Chatsirisupachai Down’ bar graph in Figure 3C and the meta-PCNA graph in Figure 3B).

To confirm that the senescence-associated gene expression signatures used here truly could identify melanocytes that had become senescent, we also analyzed published data comparing gene expression in BRAFV600E-transduced and normal human melanocytes in culture, under conditions in which the former developed definitive morphological characteristics of senescence (Pawlikowski et al., 2013; see Figure 3—figure supplement 3). Of the 23 ‘Universal Down’ signature genes, 15 were significantly differentially expressed, and 100% of these were decreased in expression. Of the 31 ‘Universal Up’ genes, 19 were significantly differentially expressed, and 84% of these were elevated in expression.

Together these data do not support the view that any sort of senescence—oncogene-induced or otherwise—is characteristic of nevus melanocytes and therefore a possible cause of their growth arrest.

### Does a cell-autonomous process arrest nevi?

As discussed above, OIS is usually presented as a cell-autonomous process (e.g. Dankort et al., 2009; Dhomen et al., 2009; Michaloglou et al., 2005; Serrano et al., 1997). The simplest cell-autonomous process that one might imagine is a probabilistic switch: Once oncogene activation commences, cells arrest with a fixed probability (per time or per cell cycle). Regardless of the molecular details, such a model makes distinctive predictions about clonal dynamics.

Consider the clonal descendants of a single oncogene-transformed founder cell. For any value of the per-cell-cycle senescence probability (which we denote here as ‘s’), how many cells should we expect that clone to contain at any given time? How many cell cycles should it take before all of the cells in most clones should have arrested? Such questions are well studied in mathematics (Athreya and Ney, 1972), and easily solved by computer simulation. For this particular problem, there are two key results (Figure 4).

Figure 4 with 1 supplement see all

First, the time after which one can expect clones to have stopped growing (e.g. when all cells will have arrested in, say, 95% or 99% of clones) is a steep function of s. If s?<?0.5, (i.e. less than a 50% chance of arrest per cell cycle), then some clones will never stop growing. If s is, say, 0.53, all clones will eventually stop growing, but one must wait 51 cell cycles before 99% of them do so (Figure 4A). Given typical lengths of postnatal mammalian cell cycles, and the fact that we observe cessation of mouse nevus growth in about 2–3 weeks, we may consider 30 to be a generous estimate for the maximum number of cell cycles by which nevi stop growing. To achieve 99% clonal arrest by 30 cell cycles, s must be around 0.56 or higher; to achieve arrest in 95% of clones, s must be greater than 0.52 (Figure 4A).

From the same simulations one may also calculate predicted distributions of clone sizes. There is a clear reciprocal relationship between mean clone size and the fraction of clones that arrest by 30 cell cycles of time (Figure 4B). For example, a value of s that enables 95% of clones to arrest produces a mean clone size of only 18.5 cells. For comparison, we estimate cell numbers per nevus to be in the range of 100–1000 cells (see Materials?and?methods).

The explanation for the small mean clone sizes produced by simulations can be appreciated by examining the full size distributions. As shown in Figure 4C–D, such distributions are extremely heavy-tailed, with a very large number of very small clones and a small number of very large clones (the histograms in Figure 4C–D are plotted with logarithmic abscissa to facilitate display of all clone sizes).

Qualitatively, this is very different from what we observed for nevi on the backs of p21 mice. Nevi displayed a mean radius of 76.8 μm (corresponding to an area of 0.019 mm2, in excellent agreement with the results of Damsky et al., 2015) and, when plotted on a logarithmic scale, individual radii displayed a Gaussian-like distribution (Figure 4E; a Gaussian shape on a logarithmic axis is usually referred to as ‘log-normal’). Interestingly, nest sizes (quantified by MPM) also seem to be log-normally distributed (Figure 4—figure supplement 1) whether at P21 (panel A) or P50 (panel B). So are the nests within human melanocytic nevi, despite the latter being an order of magnitude larger than those in mice (Figure 4—figure supplement 1C). It should be noted that using different units to represent simulation results (cell numbers) and empirical data (linear dimension) in Figure 4 and its supplement does not confound comparing the shapes of distributions, thanks to the logarithmic abscissa: as long as cell number scales as some power of linear dimension, values associated with log-transformed bins are simply scaled by a constant factor.

The above comparison of observed distributions with the results of computer simulation is not entirely fair, however: Simulations track all clones, no matter how small, whereas empirical distributions undoubtedly omit nevi with sizes below some threshold of observability. To correct for this, one can truncate simulated distributions to remove clones smaller than some threshold size. With no a priori way to know what threshold to use, we examined the entire range of plausible truncations (up to the largest clone sizes). As it turns out, the relative shapes of simulated distributions were about the same regardless of where they were truncated. The reason for this behavior can be understood by displaying simulated distributions (with bin sizes of one cell) on a log-log plot, and observing that they fall, over most clone size ranges, on a straight line (see Supplemental Material). This implies an approximately ‘power law’ relationship which, by definition, is scale-free, that?is has the same relative shape over any range of observations. In fact, for s reasonably close to 0.5, the approximate probability of observing any clone size can be shown analytically to vary inversely with the 3/2-power of size (for derivation, see Appendix 1).

These data imply that observed nevus size distributions cannot be generated by any cell-autonomous, random, time independent, one-step process. But they do not speak to whether a more complicated random process, for example, one with several steps, might suffice. To address this, one can simulate clonal evolution under multi-step models. Again, dynamic predictions can be made. First, to achieve clonal stopping times within 30 cell cycles, the minimum per-step transition probability increases with the number of steps (Figure 4F). For example, if it takes three random events to arrest growth in 99% of clones, the average probability of each event needs to be at least 0.66 per cell cycle; with five random events that number is 0.74.

Second, although distributions generated by such models still tend to be heavy-tailed (Figure 4G), they become less so as the number of steps increases (Figure 4H), gradually approaching something that looks log-normal. This makes mathematical sense: as per-cell-cycle probabilities approach unity, the system approaches a clock that simply ticks off a fixed number of cell cycles before stopping. A scenario in which all clones stop at roughly the same time, plus or minus some variation, necessarily produces a log-normal distribution, since the logarithm of cell size will be proportional to the number of cell divisions.

To determine how many independent steps would be required for a random cell-autonomous process to produce distributions that fit those we observed for nevi, we simulated up to eight random stages, over a variety of transition probabilities (Figure 4I). We subjected the results to a range of possible truncations, from 0 to 1600 cells, to mimic any observability cutoffs in the empirical data, and recorded the median clone sizes produced under each of these scenarios.

As described below (see Materials?and?methods), we estimate that the average nevus has about 500–1000 cells, but given possible errors in the estimate, we consider here a range of values between 100 and 3000 (gray-shaded area in Figure 4I). Subject to the constraint that enough clones must arrest within 30 cell cycles, and that observation thresholds cannot be so high that the observed median is less than twice the threshold, we find that, to produce nevi of even 200 cells requires 4–5 independent events (stages), depending on whether one requires 95% or 99% clonal arrest; to reach 500 cells requires 6–7 events. To reach even larger numbers—as would be found in human nevi, or in other mouse models (Chai et al., 2014)—would require even more stages.

### Does a collective process arrest nevi?

The above results indicate that, to generate in vivo-like distributions of nevi, a process something like a clock is needed, with cells either counting elapsed divisions (or time) since oncogene activation, or progressing through a sufficiently large sequence of random processes, with tightly controlled probabilities, so that the net outcome is clock-like.

Cell-autonomous counting of cell cycles (up to about 12) can occur in early, cleavage-stage embryos (Tadros and Lipshitz, 2009), but no mechanism has been described to enable growing (as opposed to merely cleaving) cells to track more than a small handful of divisions (or the equivalent amount of time). Erosion of telomeres can mark the passage of large amounts of time in some cells, but this does not seem to occur to any significant degree in nevus melanocytes (Michaloglou et al., 2005).

In contrast, if growth arrest is not cell-autonomous, but driven by cell–cell communication, then clock-like behavior is easily achieved, without any sort of intrinsic cell memory: Consider a simple communication circuit in which every cell’s arrest probability is simply a monotonically increasing function of the number of cells around it that have already arrested (Figure 5A). This mechanism describes a dynamically well-understood feedback process that normal tissues use to control size (Lander, 2011; Lander et al., 2009). Termed ‘renewal control’ (Buzi et al., 2015), because differentiated cells control the probability that progenitor cells self-renew, the process is often mediated by secreted TGF-β superfamily members such as myostatin, activin and GDF11 (Gokoffski et al., 2011; Lee et al., 2005). Because it implements the engineering principle of ‘integral negative feedback’, renewal control produces highly robust final population sizes that are independent of parameters such as cell cycle speed or the starting numbers of cells (Buzi et al., 2015; Lander, 2011; Lander et al., 2009).

Figure 5 with 1 supplement see all

When growth arrest due to renewal control is simulated as a probabilistic process (Figure 5B), the observed size distributions of clones are very close to log-normal. This is because renewal control effectively enforces cell cooperation, so that once a small fraction of a clone has arrested, the entire clone stops soon thereafter. The resulting narrow distribution of stopping times produces size distributions that are approximately log-normal, that?is that emulate a clock.

This behavior is a generic outcome of feedback control, and does not depend strongly on the details of how feedback is implemented. Similar distributions are obtained whether we model nevi as progressing reversibly or irreversibly through more than one proliferative stage, or use agent-based simulations in which renewal is controlled by the concentration of a secreted molecule that accumulates according to the laws of diffusion and local uptake (Figure 5C–D). The point of these simulations is not to argue in favor of a specific feedback mechanism, but rather to show that, where cell-autonomous mechanisms of arrest struggle to fit nevus dynamics, almost any sort of (collective) feedback does so easily.

Although nevus size distributions alone cannot shed light on the molecular details of how feedback might be implemented in nevi, it is interesting that those cells that we identify as nevus melanocytes (Figure 2) express multiple genes encoding ligands with known or suspected growth inhibitory activities, together with the receptors for those ligands. These include TGFβ superfamily members Gdf11, Gdf15, Tgfb1, Tgfb2, and Tgfb3, as well as other genes associated with growth inhibition, such as Angptl2, Angptl4, Il6, Sema3a, Sema3b, and Sema3f (Attisano and Wrana, 1996; Neufeld et al., 2016; Santulli, 2014).

The evaluation of such candidates (as well as other genes expressed at too low a level to be reliably detected by single-cell RNA sequencing) will no doubt require further study. In the meantime, we reasoned that any feedback mechanism based on secreted, diffusible factors should induce spatial correlations among clones. In particular, when clones (or subclones) get close enough to each other, they should inhibit each other’s growth, leading to a smaller final size. The distance over which such effects could occur should reflect the spatial decay length of diffusible molecules in the skin, which is thought to be on the order of no more than a few hundred microns (Chen et al., 2015). Although our data on macroscopic nevi (Figure 1D), which had been collected in a manner that included spatial coordinates, did not contain enough examples of nevus spacings in this range to test this hypothesis, our data on the nests within individual nevi did, as the median spacing between nests at postnatal day 21 was approximately 79 microns.

To assess whether nests were significantly smaller when located near other nests (especially large ones), we extracted the coordinates and nest areas from seven separate fields at P21 (representing 122 individual nests) and, modeling each nest as a disk of equivalent area, calculated the mean sizes of neighboring nests falling within successively larger annuli around each target nest. We compared the distributions of mean neighbor sizes near the 52 smallest nests (radius?<20 μm) with the equivalent distributions around the 70 largest nests.

As shown in Figure 5E–F, within annuli extending 45 μm away from the perimeters of target nests, we saw fewer examples of large neighbors (radius?>30 μm) around large nests (panel E) than around small ones (panel F). To determine whether the difference was statistically significant, we used a permutation test in which we randomly swapped nest areas (but not locations) within each field 5000 times, and recalculated the distributions. This allowed us to plot the envelope enclosing the 5th and 95th percentiles for permuted data, onto which we overlaid the observed data. Unlike the observed data, the envelopes of the permuted data (blue zones in Figure 5E–F) look similar whether target nests are large or small. Moreover, the observed data extended outside of the envelopes only for median neighbor sizes?>?30 μm, with the data for small target nevi extending well above the relevant envelope and the data for large target nevi lying at to the bottom of the envelope (Figure 5E–F, arrows). These results argue that proximity is associated with a small, but significant decrease in nest size, supporting the view that nests inhibit each other’s growth. Interestingly, when we repeated the same analysis using annulus sizes of 150 μm, differences in the sizes of neighbors of small and large nests were not seen, consistent with the view that whatever is promoting coordination among nests has a spatial range of?<150 μm (Figure 5—figure supplement 1).

## Discussion

Studies in man, mouse and fish establish that most melanocytic nevi form by mutational activation of BRAF, which triggers proliferation followed by growth arrest (Damsky and Bosenberg, 2017; Dankort et al., 2009; Dhomen et al., 2009; Kaufman et al., 2016; Michaloglou et al., 2005; Patton et al., 2005; Shain and Bastian, 2016). Nevus growth is often considered a paradigmatic example of OIS, but here we question two of the major tenets of the OIS hypothesis: that nevus melanocytes are actually senescent; and that growth arrest is a direct effect of oncogene action on the individual cell.

To assess whether nevi are senescent, we used single-cell RNA sequencing in a mouse model of Braf-driven nevus formation, comparing gene expression of nevus melanocytes with that of other cell types. Across a wide variety of gene expression signatures, especially those developed to distinguish senescence from other growth-arrested cell states, we failed to find any evidence in support of the OIS hypothesis. By gene expression criteria, nevus melanocytes were less senescent than many other normal skin cells, including non-nevus melanocytes (Figure 3, and its supplements). These results support earlier work that also questioned, based on immunohistochemical staining of human nevi for markers including lysosomal β-galactosidase, Ki67, p16INK4a (CDKN2A), γ-H2AX and p53, whether nevus melanocytes should be considered senescent (Tran et al., 2012). In agreement with other studies, we do find that Cdkn2A is highly expressed in nevus cells; it is in fact the only ‘classical’ senescence marker that clearly distinguishes nevus melanocytes from other melanocytes (Figure 3A). Yet, as others have shown, Cdkn2A is neither necessary nor sufficient for oncogene-mediated melanocyte growth arrest (Haferkamp et al., 2009; Zeng et al., 2018). Thus, to the extent that nevus melanocytes do execute even part of a common senescence program, there is little to support the view that this why they stop proliferating.

As for the question of exactly how Braf-induced nevus growth arrest occurs, Figure 6 presents a continuum of models: In model A, oncogene action elicits a cell-autonomous stress response which, after some time lag, triggers senescence that shuts proliferation down. This is the form in which the OIS hypothesis is most frequently presented.

Figure 6

In model C, growth?arrest is not a direct effect of oncogene action, but rather a consequence of growth itself. This type of feedback is commonly used by adult tissues to maintain constant size, and also enables developing tissues to produce precise numbers of differentiated cells (Kunche et al., 2016; Lander, 2011; Lander et al., 2009). Because of the collective nature of this mechanism—cells that have stopped dividing tell other cells in their vicinity to do likewise—it naturally produces semi-synchronous arrest of spatially-coherent cell clones, and the distinctive log-normal clone size distribution that comes along with that. In contrast, a purely cell-autonomous mechanism (panel A) has great difficulty producing such distributions (Figure 4), either necessitating the operation of some kind of multi-cell-cycle clock, or requiring cells to complete a long sequence of independent probabilistic events prior to arresting (Figure 4).

One can, of course, build a model in between these two extremes (model B), in which oncogenes induce growth arrest directly, but paracrine signals (i.e. SASP factors) help maintain it. If the paracrine role is important enough, this mechanism might also produce clone size distributions that are approximately log-normal, so we cannot categorically rule this model out. However, our gene expression data do not support any versions of it that have been explicitly proposed for nevi. So far, several groups (working predominantly from in vitro observations) have claimed a critical role for SASP factors in melanocyte OIS: Wajapeyee et al., 2008 argued that IGFBP7 plays a necessary role in the establishment of BRAF-V600E-induced melanocyte senescence (a conclusion disputed by some [Scurr et al., 2010]); Feuerer et al., 2019 proposed that MIA (melanoma inhibitory activity) secreted by senescent melanocytes is required to maintain senescence; and Damsky and Bosenberg have proposed that IL1, IL6, IL8 (encoded in mouse by Cxcl15), and type 1 interferons produced by nevus cells play a role in their arrest (Damsky and Bosenberg, 2017).

Our in vivo data do not support any of these hypotheses. For example, we observed that the vast majority of Igfbp7 transcripts are produced by fibroblasts and endothelial cells and that, among melanocytes, nevus melanocytes express lower levels of Igfbp7 than non-nevus melanocytes. We did not detect any Mia transcripts in nevus melanocytes, although it was expressed at detectable levels in non-nevus melanocytes and various other skin cells. Likewise, of Il1 family members, only Il1a transcripts were detected in nevus melanocytes and they were at levels lower than in many other skin cell types. Il6 was also only weakly expressed in nevus melanocytes, especially when compared with other cells. Transcripts for type one interferons were not detected in any melanocytic cells, and Cxcl15 transcripts were not detected in any skin cells at all.

Of course, the accuracy of single-cell RNA sequencing can be limited for weakly expressed genes, so we cannot completely eliminate the possibility that these factors play some role in nevus growth arrest. But given these results, and the evidence that nevus melanocytes are not senescent, we strongly favor the renewal-feedback model (model C). Adopting this model also makes it easier to accommodate long-standing evidence that nevus growth arrest is not permanent (Shain and Bastian, 2016). For example, it is known that nevi may exhibit a low level of mitoses (Glatz et al., 2010); that they can grow in response to stimuli such as UV light (Rudolph et al., 1998) or immunosuppression (Shain and Bastian, 2016) and, perhaps most tellingly, they can re-grow after incomplete surgical resection—stopping again when they reach a typical nevus size (Vilain et al., 2016). The latter result is inherently problematic for any non-feedback model, but is precisely what renewal feedback predicts (Lander, 2011; Lander et al., 2009).

Because feedback control of renewal implements a generic strategy (integral negative feedback [Lander, 2011]), it places no constraints on the molecular details of feedback, short of the fact that whatever is mediating it needs to rise with the number of cells already arrested. One possibility is that nevi are responding to some of the same signals that are used in melanocyte homeostasis. For instance, during anagen phase of the hair cycle, melanocyte stem cells produce progeny that migrate out of the hair follicle bulge as they differentiate, leaving functional stem cells behind for future cycles. A variety of experimental and pathological circumstances that allow small numbers of melanocytes to differentiate within the bulge cause differentiation and loss of the entire stem cell pool (with concomitant hair graying [Nishimura et al., 2005]). This sort of behavior—where differentiated cells drive the differentiation of their progenitors—is exactly the sort of behavior that drives feedback models of renewal (Buzi et al., 2015; Lander, 2011; Lander et al., 2009).

Nevi are but one of many types of benign, clonal, proliferative lesions that arise due to the activation of oncogenes, but rarely if ever progress to malignancy (Adashek et al., 2020). Notwithstanding the disruptive influence that oncogenes can have on cell physiology, the existence of such lesions suggest that homeostatic mechanisms persist and function at many stages along the road to cancer. New avenues for cancer prevention and treatment are likely to follow from the detailed elucidation of such mechanisms.

## Materials and methods

Key resources table

### Mouse treatment for nevus development

Request a detailed protocol

BrafV600E, Tyr-CreER (C56BL/6) mice (RRID:MGI:5902125) were genotyped by PCR as previously described (Bosenberg et al., 2006; Dankort et al., 2007). The primers used in this study are: Braf forward 5’-TGAGTATTTTTGTGGCAACTGC ?3’, Braf reverse 5’-CTCTGCTGGGAAAGCGCC ?3’, Cre forward 5’- GGTGTCCAATTTACTGACCGTACA-3’ and Cre reverse 5’- CGGATCCGCCGCATAACCAGTG ?3’. Topical administration of 4-hydroxytamoxifen (4-OHT; 25 mg/mL or 75 mg/mL in DMSO; 98% Z-isomer, Sigma-Aldrich) was administered to pups on their back and/or paws at ages P2, P3, and P4. Images of nevi on back and paw skin were taken with a digital camera at the indicated ages. Nevi from the underside of the skin were imaged using a dissection microscope. All mouse procedures were approved by UCI’s IACUC.

### Live imaging of the skin by MPM

Request a detailed protocol

Mice were sedated, shaved, and depilated with wax strips at the indicated ages (during a telogen phase) and the dorsal skin was imaged to capture the intrinsic fluorescent signal from keratin, melanin, as well as the second-harmonic-generation signal from collagen, using the LSM 510 NLO Zeiss system. Excitation was achieved with a femtosecond Titanium: Sapphire (Chameleon-Ultra, Coherent) laser at 900 nm. Emission was detected at 390–465 nm for second harmonic generation (blue) and 500–550 nm (green) and 565–650 (red) for fluorescence.

### In vivo labeling with BrdU

Request a detailed protocol

BrdU was prepared in sterile PBS at 10 mg/mL and injected intraperitoneally into mice that were 20 days old at 100 mg/kg of body weight. 24 hr later the mouse was shaved, depilated with wax strips and the skin was removed and fixed in 10% formalin for 16 hr.

### Immunofluorescence

Request a detailed protocol

Formalin fixed paraffin embedded skins were sectioned 8 μm thick, deparaffinized with Xylene, and dehydrated in a series of increasing concentration of ethanol washes. Antigen retrieval was performed with 10 mM citric acid buffer at pH 6.0 for 10 min in a steamer. Samples were washed with PBS, incubated with TrueBlack for 30 s to reduce autofluorescence, and washed again with PBS. All antibodies were diluted at a 1:500 and incubated overnight at 4°C. Samples were washed and incubated with the appropriate secondary antibody. Melanocytes were identified with a Pmel antibody (EP4863(2); ab137078, Abcam; RRID:AB_2732921). Cells that incorporated BrdU were visualized with a BrdU antibody (ab6326, Abcam; RRID:AB_305426).

### Cell isolation for single-cell RNA sequencing

Request a detailed protocol

BrafWT, Tyr-CreER or BrafV600E, Tyr-CreER mice were euthanized at either P30 (n?=?2 of each genotype) or P50 (n?=?3 of each genotype), shaved, and depilated. A 2 × 3 cm section of the dorsal skin was removed, and the fat scraped off from the underside. The piece was then diced into smaller pieces and suspended in dissociation buffer (RPMI, liberase 0.25 mg/mL, Hepes 23.2 mM, Sodium Pyruvate 2.32 mM, Collagenase:Dispase 1 mg/mL) for 50 min at 37°C with gentle agitation. After incubation, DNaseI (232U) was added for 10 min and then inactivated with fetal bovine serum and EDTA (1 mM). The tissue suspension was further dissociated mechanically with the GentleMACS using the setting m_imptumor_04.01, which runs for 37 s at various speeds. Single-cell suspensions were filtered twice through a 70 μm strainer and dead cells removed by centrifugation at 300 x g for 15 min. The live cells were washed with 0.04% UltraPure BSA:PBS buffer, gently re-suspended in the same buffer, and counted using trypan blue.

### Library preparation for single-cell RNA sequencing and analysis

Request a detailed protocol

Libraries were prepared using the Chromium Single Cell 3’ v2 protocol (10X Genomics). Briefly, individual cells and gel beads were encapsulated in oil droplets where cells were lysed and mRNA was reverse transcribed to 10X barcoded cDNA. Adapters were ligated to the cDNA followed by the addition of the sample index. Prepared libraries were sequenced using paired end 100 cycles chemistry for the Illumina HiSeq 4000. FASTQ files were generated from Illumina’s binary base call raw output with Cell Ranger’s (v2.1.0; RRID:SCR_017344) `cellranger mkfastq` command and the count matrix for each sample was produced with `cellranger count`. All ten samples (four samples from P30 [two control (wild type) and two mutant] and six samples from P50 [three control and three mutant]) were aggregated together with the `cellranger aggr` command to produce one count matrix that includes all samples. Data analysis was performed with Scanpy [v1.3.6; RRID:SCR_018139] (Wolf et al., 2018). Cells with fewer than 200 detected genes, and genes detected in less than three cells, were discarded. We calculated the percent mitochondrial gene expression and kept cells with less than 13% mitochondrial gene expression, and cells with fewer than 4000 genes/cell (35,141 cells). Each cell was normalized to total counts over all genes. In the final preprocessing step, we regressed out cell-cell variation driven by mitochondrial gene expression and the number of detected UMI. To identify clusters, we first performed principal component analysis on log-transformed data, using highly variable genes, Louvain clustering (Levine et al., 2015), and visualization with t-distributed stochastic neighbor embedding (tSNE).

### Quantification of nevus and nest size and cell content

Request a detailed protocol

To quantify the sizes of nevi in mice, dorsal skin was excised and the underside visualized using a dissecting microscope. Nevi were traced, and area calculated using ImageJ. Nest sizes were quantified in live mice by MPM. Sizes of human nests were measured from histological samples (n?=?5) obtained from the UCI Department of Dermatology. Samples were stained with hematoxylin and eosin and imaged with a microscope. A dermatologist manually identified the nests on each slide, and nest area was quantified using ImageJ. Human studies were performed under IRB protocol HS# 2019–5054.

Estimates of cell numbers for mouse nevi were obtained in two different ways: First, we used estimates from?Chai et al., 2014 for melanocytic nuclei per square area of mouse nevus, together with our observed median nevus radius of 76.8 μm; this approach led to an estimate of 897 cells/nevus. As the data of Chai et al., 2014?come from a different genetic model, we also estimated cell number as follows: Using 8 μm sections of back skin from Albino BrafV600E mice, we used fluorescence microscopy to measure the sizes of 194 Pmel-stained melanocytes within the nests of nevi, obtaining an average cell diameter of 5.68 μm, and counted approximately 14.4 cells per 104 μm3 of nest. In pigmented animals, we measured by MPM an average nest cross sectional area of 1385 μm2, an average nest volume of 38792 μm3, and an average number of nests per nevi of approximately 12, yielding an estimate of 672 cells/nevus. Given uncertainties in these measures, analyses in the manuscript take into account the possibility of an average that falls anywhere between 100 and 1000 cells.

### Simulations and agent-based modeling

Request a detailed protocol

Stochastic, non-spatial simulations of renewal control were obtained by Monte Carlo simulation, in which cells duplicated every cell cycle, and then chose randomly whether to differentiate or continue dividing according to a probability modified by feedback from non-dividing cells. A Hill function, with Hill coefficient?=?1, was used to represent the feedback.

To model feedback in a spatial context, we used CompuCell3D (RRID:SCR_003052), an open-source platform for Cellular Potts modeling (Swat et al., 2012). In CPM, every generalized object or ‘cell’ is associated with a list of attributes such as cell type, surface area, volume, etc. These enter into the calculation of an effective energy, which can be summarized as the sum of the contact energy between neighboring cells and the effective energy due to volume constraints.

Simulations were initialized by seeding a single cell, with a size of 25 pixel2, in the center of a 300 × 300 pixel lattice, which grew and divided according to rules and parameters summarized in Source data 2. To add variability to cell growth, cells randomly chose one of three different growth rates after every cell division. To add variability to cell division times, cells randomly chose a target area, between 72 and 80 pixel2, at which to divide. Growth rates were chosen to be sufficiently slow that the mean time between cell divisions came out to approximately 382 time steps. At division, each cell was divided in half by a randomly-oriented division plane.

Upon division, a cell either remained dividing or became permanently arrested. All cells had a minimum 1% probability of arrest per cell division. Once an arrested cell was generated, it began continuous secretion of a signaling molecule that diffuses and promotes the transition from dividing to non-dividing (Kunche et al., 2016). Diffusion and decay of the feedback factor was modeled deterministically, with parameters chosen to produce a steady state decay length of 15 pixels. The concentration of this factor at the center of mass of each cell then augmented the arrest probability of that cell by an amount determined by a Hill function (see Source data 2).

### Statistical analysis

Request a detailed protocol

Statistical analyses for single-cell RNA sequencing were performed using Scanpy (RRID:SCR_018139). Other statistical testing was done using Mathematica (RRID:SCR_014448). For the spatial analysis in Figure 5E–F, nest areas in each field were randomly swapped, with positions held constant, 5000 times, and the distributions of neighboring nest locations and sizes recalculated each time. This allowed us to generate an envelope enclosing the 5th and 95th percentiles for the permuted data, at each target nest size, and compare the observed data with the bounds of that envelope.

## Appendix 1

### Modeling Oncogene-Induced Senescence as a Cell Autonomous Process—Mathematical Results

Consider the clonal descendants of an oncogene-transformed cell. The simplest model of a cell-autonomous oncogene-induced arrest process is one in which cells replicate at a constant rate, and undergo senescence with a constant probability s per cell cycle. Once all the cells in a clone have undergone senescence, we refer to the clone as ‘extinct’. For any given s, we wish to find the probability that extinction has occurred by any given time, as well as the distribution of clone sizes that should be expected at that time.

If the decision to senesce is independent in each daughter cell at each cell division, then this scenario describes a branching process in which each cell produces two senescent cells with probability s2, two dividing cells with probability (1?s)2, and one dividing and one senescent cell with probability 2s(1?s). The theory of branching processes may then be used to obtain the probability generating function (PGF) for the number of dividing cells after n cell cycles. The offspring distribution of a single cell as described above may be written as:

$\mathit{\text{f(z)}}={\mathit{\text{s}}}^{2}+2\mathit{\text{s}}\left(1?\mathit{\text{s}}\right)\mathit{\text{z}}+\left(1?\mathit{\text{s}}{\right)}^{2}{\mathit{\text{z}}}^{2}$

where z is a dummy variable whose exponent indicates the number of dividing (non-senescent) cells produced by an event, and the coefficient in front of each zn is the probability of that event. It can be shown that the PGF for the number of dividing cells after n cell cycles is equal to the n-th composition of the offspring distribution f unto itself (Hawkins and Ulam, 1944). For example, the PGF for the number of dividing cells after 2 cell cycles is:

$\begin{array}{ll}\mathit{\text{F}}\left(2,\mathit{\text{z}}\right)& =f\left(f\left(z\right)\right)={s}^{2}+2s\left(1?s\right)\left({s}^{2}+2s\left(1?s\right)z+\left(1?{s}^{2}\right){z}^{2}\right)\\ & +\left(1?s{\right)}^{2}\left({s}^{2}+2s\left(1?s\right)z+\left(1?s{\right)}^{2}{z}^{2}{\right)}^{2}\end{array}$

And in general:

$\mathit{\text{F}}\left(\mathit{\text{n}},\mathit{\text{z}}\right)=\mathit{\text{F}}\left(1,\mathit{\text{F}}\left(\mathit{\text{n}}?1,\mathit{\text{z}}\right)\right)=\mathit{\text{f}}\left(\mathit{\text{F}}\left(\mathit{\text{n}}?1,\mathit{\text{z}}\right)\right)$

Thus F(n, 0) is the cumulative probability that there are 0 dividing cells—that?is a clone has extinguished—after n cell cycles.

The theory of branching processes also tells us that the probability that a clone eventually goes extinct (over the long run) is just the smallest non-negative root of z?=?f(z). We solve the equation:

z = s2 + 2 s(1?s)z + (1?s)2z2 and obtain:?z = s2/(1?s)2, 0?≤ s?<?0.5 and z?=?1, s?≥ 0.5. This confirms that only if s?≥ 0.5 is eventual extinction of all clones guaranteed.

### Monte Carlo Simulations

The behavior of the above branching processes may be easily observed using Monte Carlo Simulation, seeding each clone with a single dividing cell and specifying the value of s. In the simulations show in Figure 4, cells divided synchronously and s was used to determine the fate of each daughter cell after each division. For each simulation it was recorded (1) whether the clone went extinct, and if so, after how many cell cycles it did so and (2) the number of cells at the time of extinction (or the end of the simulation).

### Distribution of clone sizes at extinction

To derive the distribution of clone sizes at extinction it is helpful to think about clonal development not in terms of the number n of cell cycle times that have elapsed, but in terms of the cumulative number of cell division events that have occurred up to any given time within a clone, which we will represent as η. Assuming clones begin from one cell, the number of cells at extinction (Tc) will simply be Tc?=?1 + η;

To find the expected distribution of values of Tc, we begin by computing the probability that a clone goes extinct at a given value of η, which we will call pE(η). Recalling that two choices are made at every division (one per daughter cell), it may be seen that to extinguish at η requires 2 η choices, exactly η +1 of which are choices to become senescent, and η ?1 events are choices to remain proliferative. This implies that:

$\begin{array}{ll}{p}_{E}\left(1\right)& ={\mathit{\text{A}}}_{1}{\mathit{\text{s}}}^{2}\\ {p}_{E}\left(2\right)& ={\mathit{\text{A}}}_{2}\left(1?\mathit{\text{s}}\right){\mathit{\text{s}}}^{3}\\ {p}_{E}\left(3\right)& ={\mathit{\text{A}}}_{3}\left(1?\mathit{\text{s}}{\right)}^{2}{\mathit{\text{s}}}^{4}\\ ...\\ {p}_{E}\left(\eta \right)& ={\mathit{\text{A}}}_{\eta }\left(1?\mathit{\text{s}}{\right)}^{\eta ?1}{\mathit{\text{s}}}^{\eta +1}\end{array}$

where the coefficients Aη are constants that capture the number of different ways that each combination of choices of s and 1?s can happen. Aη can be thought of as the number of unique full binary trees that end with η + 1 senescent cells. As senescent cells are the dead ends in those trees, they are referred to as ‘leaves’ of the tree. The sequence that counts the number of unique full binary trees is called the Catalan numbers. Ck, the kth Catalan number is the number of unique binary trees with k+1 leaves. It is defined by:

${C}_{k}=\frac{1}{k+1}\left(\begin{array}{c}2k\\ k\end{array}\right)$

Accordingly,

${p}_{E}\left(\eta \right)=\frac{1}{\eta +1}\left(\begin{array}{c}2\eta \\ \eta \end{array}\right){\left(1-s\right)}^{\eta -1}{s}^{\eta +1}$

Consequently, the probability of Tc?=?m cells at extinction is simply

${p}_{{T}_{c}}\left(m\right)={p}_{E}\left(m-1\right)$

In Appendix 1—figure 1, panel A, we plot the value of pE(η) as a function of η for different values of s, using logarithmic axes. Notice that for s sufficiently close to 0.5, the relationship approximates a line of slope –3/2. Thus, the approximate probability of finding a clone of size m varies inversely with the 3/2-power of m. In panel B, the analytical results for s?=?0.56 are superimposed on results obtained by Monte Carlo simulation of 500,000 cases.

Appendix 1—figure 1

## References

1. 1
2. 2
3. 3
Branching Processes
(1972)
Dover.
4. 4
5. 5
6. 6
7. 7
8. 8
9. 9
Cellular senescence: when bad things happen to good cells (2007)
Nature Reviews Molecular Cell Biology 8:729–740.
https://doi.org/10.1038/nrm2233
10. 10
11. 11
12. 12
13. 13
14. 14
15. 15
16. 16
17. 17
18. 18
19. 19
20. 20
21. 21
22. 22
23. 23
24. 24
25. 25
26. 26
Theory of Multiplicative Processes I
(1944)
27. 27
28. 28
29. 29
30. 30
31. 31
32. 32
33. 33
34. 34
35. 35
36. 36
The essence of senescence (2010)
Genes & Development 24:2463–2479.
37. 37
38. 38
39. 39
40. 40
41. 41
42. 42
43. 43
44. 44
45. 45
46. 46
47. 47
48. 48
49. 49
50. 50
51. 51
52. 52
53. 53
54. 54
From melanocytes to melanomas (2016)
Nature Reviews Cancer 16:345–358.
https://doi.org/10.1038/nrc.2016.37
55. 55
56. 56
57. 57
58. 58
59. 59
60. 60
61. 61
62. 62
63. 63
64. 64
65. 65
66. 66
67. 67

## Decision letter

1. Richard M White
Senior Editor; Memorial Sloan Kettering Cancer Center, United States
2. C Daniela Robles-Espinoza
Reviewing Editor; International Laboratory for Human Genome Research, Mexico
Reviewer
4. Mariana Gómez-Schiavon
Reviewer; UCSF, United States
5. Heinz Arnheiter
Reviewer; NINDS, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

In this elegant manuscript, Ruiz-Vega and collaborators test the generally accepted hypothesis that nevi are growth-arrested due to oncogene-induced senescence, which is a cell-autonomous process thought to be triggered by the oncogenic mutation BRAF V600E, present in the majority of nevi. They analyse single-cell transcriptome data from mice nevi and do not find any evidence of these cells being senescent. Rather, and showing great creativity, they find through the use of mathematical models that the distribution of nevi sizes is better explained by cooperative mechanisms, of the kind used in size control of normal tissues.

Decision letter after peer review:

Thank you for submitting your article "Dynamics of nevus development implicate cell cooperation in the growth arrest of transformed melanocytes" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Richard White as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Peter Adams (Reviewer #1); Mariana Gómez-Schiavon (Reviewer #2); Heinz Arnheiter (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (http://www.asadoresteatre.com/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

In this manuscript, the authors present evidence that benign nevi do not simply reflect oncogene-induced senescence. The authors propose, based on single cell RNA-seq of BRAF oncogene-expressing mouse melanocytes from mouse skin, that nevus growth arrest is more likely related to the cell interactions that mediate size control in normal tissues, than to any cell-autonomous, "oncogene-induced" program of senescence. Then, the authors propose a basic integral feedback control as the underlying mechanism behind nevi growth regulation. The editors and reviewers consider that this manuscript is exquisitely done, clear and convincing, and an important contribution to an important topic. As the authors note, the prevailing dogma is that nevus melanocytes are senescent.

All issues raised by the reviewers can be considered minor, and are listed below.

1) About the concept behind the paper:

The reviewers highlight two conceptual questions that would be in need of an answer:

a) what are the suggested mechanisms leading to growth arrest of the first cells in a clone? If this were due to some probabilistic event that, however, still depended on oncogene activation, would that not favor the hybrid model B (at least for these initial cells) rather than the authors' preferred model C?

b) While the model is clearly exclusively or nearly exclusively cell non-autonomous, it is still nevus-autonomous. The authors argue that the OIS model is mostly based on in vitro observations, but should a nevus-autonomous model not also be reproducible in vitro? If it were not, then one would have to invoke an influence of the nevus microenvironment, which would leave room for, among other factors, a senescence-associated program generated by the surrounding cells.

What would be the authors' opinion on these two points?

2) Regarding the single-cell RNAseq data:

a) The authors show data to indicate that BRAF-expressing nevi are no-proliferating, compared to WT hair follicle melanocytes. What about BRAF-expressing melanocytes in hair follicles? Does the system permit an analysis?

b) The authors compare their single cell RNA-seq data to many different senescence datasets and signatures – but most of these are from fibroblasts, including Hernandez-Segura et al., 2017. However, Pawlikowski et al., 2013, reported gene expression profiling specifically in primary human melanocytes expressing activated BRAF. Can a comparison to the dataset in this paper be made, as it is judged to be better suited?

3) Regarding the mathematical model:

a) When the authors are exploring the hypothesis of the multiple stages for a cell-autonomous growth arrest (on Figure 4 legend, and in the Results), the authors use "observation thresholds cannot be so high that the observed median is less than twice the threshold" as an argument to reject some models with few stages. Can you please explain what the logic is behind this restriction? What are the implications?

b) On Figure 4C-D, a reviewer found it difficult to compare the two models. Can you please explain what the intended message is? It might be clearer if the two distributions are overlapped.

c) A reviewer thought that Figure 4I and its legend are confusing. Particularly, the sentence "Curves become dashed at the point where the observability threshold exceeds 50% of the median cell number" is unclear. The main text does a better job explaining it, but the figure panel and legend can be improved. Can this be clarified, please?

d) The authors refer to the number of "Monte Carlo time steps" in Figure 5 and supplementary information. Can the authors please explain how this is meaningful? Why is this number relevant? A reviewer writes, "The actual kinetic parameter values can have implications on the system behavior when stochasticity is considered, and these parameters being congruent with the expected cell cycle length is relevant. Nevertheless, the relevance of the Monte Carlo time steps on the simulations isn't clear to me."

e) The authors say "repeated the same analysis using annulus sizes of 150 μm, differences in the sizes of neighbors of small and large nests were not seen." Can you please include these results in the supplementary material?

https://doi.org/10.7554/eLife.61026.sa1

## Author response

[…] All issues raised by the reviewers can be considered minor, and are listed below.

1) About the concept behind the paper:

The reviewers highlight two conceptual questions that would be in need of an answer:

a) What are the suggested mechanisms leading to growth arrest of the first cells in a clone? If this were due to some probabilistic event that, however, still depended on oncogene activation, would that not favor the hybrid model B (at least for these initial cells) rather than the authors' preferred model C?

It is true that, for a mechanism in which growth arrest is promoted by a signal from already-arrested cells, something must account for the arrest of the very first cell (or cells). However, that event need not be related to oncogene activation. In our model, as in most models of tissue growth control, dividing cells are simply assumed to have some pre-existing, small probability of spontaneously arresting. Over a fairly wide range of that probability there is little effect on the final size at which clones achieve complete arrest (Kunche et al., 2016)

b) While the model is clearly exclusively or nearly exclusively cell non-autonomous, it is still nevus-autonomous. The authors argue that the OIS model is mostly based on in vitro observations, but should a nevus-autonomous model not also be reproducible in vitro? If it were not, then one would have to invoke an influence of the nevus microenvironment, which would leave room for, among other factors, a senescence-associated program generated by the surrounding cells.

Autonomous collective nevus growth arrest should, in principle, be reproducible in vitro, but in tissue culture many conditions are usually very different from what occurs in vivo. For example, the amount of extracellular volume into which secreted factors are diluted in vitro is typically orders of magnitude more than what a cell might encounter in vivo. Moreover, it is well known that dimensionality (2D versus 3D) has a very large effect on cell behaviors. These factors would make it difficult to interpret any negative in vitro result. Regardless, we strongly agree with the reviewer that the mechanism of growth arrest need not be nevus-autonomous. Cells in the microenvironment could very well have an essential role in relaying or amplifying signals that originate with melanocytes. We do not believe that statements in this manuscript conflict with that possibility. Yet, we would still argue against terming such a mechanism “senescence-associated”, given the evidence we present that arrested nevus cells are not senescent.

What would be the authors' opinion on these two points?

2) Regarding the single-cell RNAseq data:

a) The authors show data to indicate that BRAF-expressing nevi are no-proliferating, compared to WT hair follicle melanocytes. What about BRAF-expressing melanocytes in hair follicles? Does the system permit an analysis?

Indeed, the cluster we identify as “hair follicle melanocytes” (cluster 3) contains both cells from wildtype (untreated) and mutant (Braf-activated) mice, so a comparison of gene expression is possible. We have added this (using the metaPCNA signature to evaluate proliferation) as a new supplementary figure (Figure 3—figure supplement 2). With the caveat that, as we split up cell groups into smaller subsets we lose statistical power, we failed to observe any statistically significant difference in gene expression between wildtype and mutant cells of this cluster. There was a small trend toward lower proliferative gene expression in the mutant group, so perhaps with more data such an effect might become significant. Regardless, the data show that Braf activation does not drive growth arrest in every context, consistent with the view that cell-cell interactions, rather than Braf itself, are the direct cause of growth arrest.

b) The authors compare their single cell RNA-seq data to many different senescence datasets and signatures – but most of these are from fibroblasts, including Hernandez-Segura et al., 2017. However, Pawlikowski et al., 2013, reported gene expression profiling specifically in primary human melanocytes expressing activated BRAF. Can a comparison to the dataset in this paper be made, as it is judged to be better suited?

Pawlikowski et al., compared gene expression (by microarray) between normal and BRAFV600E-transduced human melanocytes in vitro. They showed that cultured, BRAF-activated melanocytes display classical morphological features of senescence. This gave us the opportunity to compare how these cells differ from normal cultured melanocytes, and how our nevus cells differ from normal in vivo melanocytes. To the extent that we observe shared differences, we could ask whether those differences are related or unrelated to genes thought to be associated with senescence.

A new supplementary figure (Figure 3—figure supplement 3) addresses that question. In it, we use scatter plots in which the statistically significant gene expression changes seen by Pawlikowski (there are 5063 of them that have unambiguous mouse orthologues) are contrasted with the differences we saw between nevus melanocytes and wildtype cells from either cluster 2 (putative melanocyte stem cells) or cluster 3 (hair follicle melanocytes).

Overall (Figure 3—figure supplement 3A-B), we see only a weak correlation in fold-changes when we compare the Pawlikowski et al., data with our nevus-versus-cluster 3 fold-changes (panel B), a correlation that improves slightly if we restrict the data from the present study to genes that passed a significance cutoff of adjusted p<0.05. That correlation mainly derives from a common subset of genes that is strongly downregulated in both the Pawlikowski samples and ours. When we examine these genes closely, we find many are proliferation-associated (panel D), and others are Mitf-target genes (panel E; interestingly, in the Pawlikowski data, Mitf and all of its target genes are coordinately downregulated, whereas in our case Mitf itself is not, and only about half of its targets are). The remaining genes are an interesting set (panel C), but do not overlap in any obvious way with gene sets that have been associated with senescence. On the other hand, when we look specifically at genes that others (specifically, Hernandez-Segura et al., 2017) have associated with senescence (panel F), we see that they very nicely characterize the changes (in both directions) in the Pawlikowski data—but not our data—confirming that the gene expression signatures we chose to use in this manuscript really do apply to melanocytes when they are actually senescent. We thank the reviewer for suggesting this informative analysis.

3) Regarding the mathematical model:

a) When the authors are exploring the hypothesis of the multiple stages for a cell-autonomous growth arrest (on Figure 4 legend, and in the Results), the authors use "observation thresholds cannot be so high that the observed median is less than twice the threshold" as an argument to reject some models with few stages. Can you please explain what the logic is behind this restriction? What are the implications?

We apologize for the lack of clarity in the legend to Figure 4, and we have revised it to explain this issue better. What we mean is this: whenever making measurements of a single dispersed quantity many times, one notices an observed dynamic range, i.e. how big the average measurements are compared with the smallest. This number contains empirical information about limits of observability. In modeling, we realized we should stop considering theoretical observation thresholds once they required the observed dynamic range to be implausibly narrow. As shown in Figure 4E, in our in vivo measurements, we see nevi as large as 150 μm in radius, a median of 76.8 μm, and can detect them as small as 40 μm in radius. Assuming a quadratic relationship between radius and cell number, that would represent a dynamic range of observation of 3.7-fold between the median and the smallest. In the interests of being conservative, we chose to consider models implausible once they required a dynamic range of less than two-fold. The dashed lines are used to indicate when modeling parameters produce regimes that should be considered implausible.

b) On Figure 4C-D, a reviewer found it difficult to compare the two models. Can you please explain what the intended message is? It might be clearer if the two distributions are overlapped.

As with the above comment, we were trying to make the legend brief, and ended up being unclear. In panel C, the arrest probability was set at the level required to ensure 99% arrest, whereas in D it was set at the level required to ensure 95% arrest (see panel A). The figure legend has been revised to make this clear. The point is to show that the two distributions are very similar, i.e. the observed heavy-tailed shape is general, and not an artifact of choosing a particular value for the arrest probability. Since we are trying to show that the two distributions are very similar, we find that an overlapped graph is not very visually appealing.

A reviewer thought that Figure 4I and its legend are confusing. Particularly, the sentence "Curves become dashed at the point where the observability threshold exceeds 50% of the median cell number" is unclear. The main text does a better job explaining it, but the figure panel and legend can be improved. Can this be clarified, please?

This is essentially the same as comment 3a. Please see the response above.

d) The authors refer to the number of "Monte Carlo time steps" in Figure 5 and supplementary information. Can the authors please explain how this is meaningful? Why is this number relevant? A reviewer writes, "The actual kinetic parameter values can have implications on the system behavior when stochasticity is considered, and these parameters being congruent with the expected cell cycle length is relevant. Nevertheless, the relevance of the Monte Carlo time steps on the simulations isn't clear to me."

In cellular Potts modeling, cells grow and are displaced over time, and everything about them gets repeatedly recalculated –these are the Monte Carlo time steps. Consistent with reality, cell division in our models is probabilistic—cells experience a nonlinearly increasing likelihood of dividing over time and with increasing size. As a result, there is no exact relationship between Monte Carlo step and cell cycle length, only an empirically-observed, approximate one (which we report as 382 time steps per cell cycle). The importance of reporting this value is to show that times associated with cell growth and division are very much longer than the time intervals over which we re-calculate. As in differential equation modeling, the key to guarding against numerical artifacts is to ensure that calculation time steps are short relative to the characteristic times associated with the events being modeled.

e) The authors say "repeated the same analysis using annulus sizes of 150 μm, differences in the sizes of neighbors of small and large nests were not seen." Can you please include these results in the supplementary material?

We are happy to oblige. The data are now presented as Figure 5—figure supplement 1.

https://doi.org/10.7554/eLife.61026.sa2

## Article and author information

### Author details

1. #### Rolando Ruiz-Vega

1. Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
2. Department of Developmental and Cell Biology, University of California, Irvine, Irvine, United States
##### Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Writing - original draft, Writing - review and editing
##### Competing interests
No competing interests declared
2. #### Chi-Fen Chen

Department of Dermatology, University of California, Irvine, Irvine, United States
##### Contribution
Data curation, Formal analysis, Investigation
##### Competing interests
No competing interests declared

Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
##### Contribution
Data curation, Investigation
##### Competing interests
No competing interests declared
4. #### Priya Vasudeva

Department of Dermatology, University of California, Irvine, Irvine, United States
##### Contribution
Data curation, Investigation
##### Competing interests
No competing interests declared
5. #### Tatiana B Krasieva

Beckman Laser Institute, University of California, Irvine, Irvine, United States
##### Contribution
Data curation, Investigation
##### Competing interests
No competing interests declared
6. #### Jessica Shiu

Department of Dermatology, University of California, Irvine, Irvine, United States
##### Contribution
Data curation, Investigation
##### Competing interests
No competing interests declared
7. #### Michael G Caldwell

Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
##### Contribution
Formal analysis, Investigation
##### Competing interests
No competing interests declared
8. #### Huaming Yan

Department of Mathematics, University of California, Irvine, Irvine, United States
##### Contribution
Formal analysis
##### Competing interests
No competing interests declared
9. #### John Lowengrub

1. Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
2. Department of Mathematics, University of California, Irvine, Irvine, United States
Methodology
##### Competing interests
No competing interests declared
10. #### Anand K Ganesan

1. Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
2. Department of Dermatology, University of California, Irvine, Irvine, United States
##### Contribution
Conceptualization, Supervision, Funding acquisition, Writing - original draft, Writing - review and editing
##### Competing interests
No competing interests declared
11. #### Arthur D Lander

1. Center for Complex Biological Systems, University of California, Irvine, Irvine, United States
2. Department of Developmental and Cell Biology, University of California, Irvine, Irvine, United States
3. Department of Biological Chemistry, University of California, Irvine, Irvine, United States
##### Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing
##### Competing interests
No competing interests declared

### Funding

#### National Cancer Institute (U54CA217378)

• John Lowengrub
• Anand K Ganesan
• Arthur D Lander

#### National Institute of Arthritis and Musculoskeletal and Skin Diseases (P30AR075047)

• Arthur D Lander
• Anand K Ganesan
• John Lowengrub

#### National Institute of Biomedical Imaging and Bioengineering (T32EB009418)

• Michael G Caldwell

#### University of California (Postdoc Fellowship)

• Rolando Ruiz-Vega

#### National Academies of Sciences, Engineering, and Medicine (Postdoc Fellowship)

• Rolando Ruiz-Vega

#### National Cancer Institute (T32CA009054)

• Rolando Ruiz-Vega

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

### Acknowledgements

The authors are grateful to Eric Mjolsness, Allon Klein, Randy Heiland and Paul Macklin for helpful discussions. Funding: This work was supported by National Institute of Health (NIH) grants U54CA217378?and?P30AR075047 (ADL, AKG, JL). MGC was supported by NiH-T32EB009418. RRV was supported by the UC Presidents fellowship, FORD Foundation Fellowship and NIH-T32CA009054. Author Declaration of interests: No conflict of interest declared by any author. Data Availability: Raw sequence data are available at GEO (GSE154679).

### Ethics

Animal experimentation: This study performed in strict accordance with the recommendation from University Laboratory Animal Resources (ULAR). All the animals were handled according to approved institutional animal care and use committee (IACUC) protocol (#AUP-17-230) at the University of California Irvine.

### Senior Editor

1. Richard M White, Memorial Sloan Kettering Cancer Center, United States

### Reviewing Editor

1. C Daniela Robles-Espinoza, International Laboratory for Human Genome Research, Mexico

### Reviewers

2. Mariana Gómez-Schiavon, UCSF, United States
3. Heinz Arnheiter, NINDS, United States

### Publication history

2. Accepted: September 11, 2020
3. Version of Record published: October 13, 2020 (version 1)

? 2020, Ruiz-Vega et al.

## Metrics

• 1,045
Page views
• 119
• 0
Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

A two-part list of links to download the article, or parts of the article, in various formats.

### Open citations (links to open the citations from this article in various online reference manager services)

1. Cancer Biology

# eIF4E S209 phosphorylation licenses myc- and stress-driven oncogenesis

Hang Ruan et al.
Research Article Updated

To better understand a role of eIF4E S209 in oncogenic translation, we generated EIF4ES209A/+ heterozygous knockin (4EKI) HCT 116 human colorectal cancer (CRC) cells. 4EKI had little impact on total eIF4E levels, cap binding or global translation, but markedly reduced HCT 116 cell growth in spheroids and mice, and CRC organoid growth. 4EKI strongly inhibited Myc and ATF4 translation, the integrated stress response (ISR)-dependent glutamine metabolic signature, AKT activation and proliferation in vivo. 4EKI inhibited polyposis in ApcMin/+ mice by suppressing Myc protein and AKT activation. Furthermore, p-eIF4E was highly elevated in CRC precursor lesions in mouse and human. p-eIF4E cooperated with mutant KRAS to promote Myc and ISR-dependent glutamine addiction in various CRC cell lines, characterized by increased cell death, transcriptomic heterogeneity and immune suppression upon deprivation. These findings demonstrate a critical role of eIF4E S209-dependent translation in Myc and stress-driven oncogenesis and as a potential therapeutic vulnerability.

1. Cancer Biology

# Inhibiting IRE1α-endonuclease activity decreases tumor burden in a mouse model for hepatocellular carcinoma

Nata?a Pavlovi? et al.
Research Article Updated

Hepatocellular carcinoma (HCC) is a liver tumor that usually arises in patients with cirrhosis. Hepatic stellate cells are key players in the progression of HCC, as they create a fibrotic micro-environment and produce growth factors and cytokines that enhance tumor cell proliferation and migration. We assessed the role of endoplasmic reticulum (ER) stress in the cross-talk between stellate cells and HCC cells. Mice with a fibrotic HCC were treated with the IRE1α-inhibitor 4μ8C, which reduced tumor burden and collagen deposition. By co-culturing HCC-cells with stellate cells, we found that HCC-cells activate IREα in stellate cells, thereby contributing to their activation. Inhibiting IRE1α blocked stellate cell activation, which then decreased proliferation and migration of tumor cells in different in vitro 2D and 3D co-cultures. In addition, we also observed cell-line-specific direct effects of inhibiting IRE1α in tumor cells.