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

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.

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. 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.
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.

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 Braf V600E 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).
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, Braf V600E -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 (

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 BRAFtransformation from quiescent fibroblasts ('Kuilman', 21 genes [Kuilman et al., 2008]  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.
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 nevusassociated 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 log 2 -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 1-2, 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 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 log 2 -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 proliferationassociated genes in non-nevus, hair follicle melanocytes (Mel 3) when compared with nevus melanocytes, was consistent between Braf V600E -expressing and control mice (Figure 3-figure supplement 2), suggesting that tissue context plays a role in whether Braf V600E -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 BRAF V600E -transduced and normal human melanocytes in culture, under conditions in  Figure 2A and B were averaged by cluster and, for each gene, expressed as a ratio to the average expression level of that gene in the entire skin sample ('Whole Skin') or just the melanocyte clusters ('Melanocytes'). The heat map displays the logarithm of that ratio, with blue representing upregulation and red downregulation, relative to the appropriate average. Two proposed signatures for genes upregulated in senescence ('Classical' and 'Universal Up') are shown (For other signatures, see   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 cellautonomous 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).
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 mm (corresponding to an area of 0.019 mm 2 , 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/2power 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 Source data 1. Raw data used to generate histogram in Figure 4E.  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-b 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).
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 agentbased 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 TGFb 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). (E-F) Using spatial coordinates and areas of 122 nests in seven individual fields at P21, nests were modeled as disks of equivalent area, and mean sizes of neighboring nests falling within successively larger annuli around each target nest were determined. Distributions of mean neighbor sizes up to 45 mm away from (E) large (radius >20 mm) and (F) small (radius <20 mm) nests (histograms) are compared with a 'null distribution' derived by random permutation (blue fields). Arrows show deviations in neighbor size distribution greater than expected at random, and substantially different for large versus small nests. The online version of this article includes the following source data and figure supplement(s) for figure 5: Source data 1. Raw data used to generate histograms and permutation tests in Figure 5E-F. 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 mm) with the equivalent distributions around the 70 largest nests.
As shown in Figure 5E-F, within annuli extending 45 mm away from the perimeters of target nests, we saw fewer examples of large neighbors (radius >30 mm) 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 5 th and 95 th 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 mm, 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 mm, 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 mm ( 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 b-galactosidase, Ki67, p16 INK4a (CDKN2A), g-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 oncogenemediated 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.
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 semisynchronous 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. Mouse treatment for nevus development Braf V600E , 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'-CGGATCCGCCGCA TAACCAGTG À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
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
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
Formalin fixed paraffin embedded skins were sectioned 8 mm 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
Braf WT , Tyr-CreER or Braf V600E , 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 mm 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 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
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 mm; 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 mm sections of back skin from Albino Braf V600E 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 mm, and counted approximately 14.4 cells per 10 4 mm 3 of nest. In pigmented animals, we measured by MPM an average nest cross sectional area of 1385 mm 2 , an average nest volume of 38792 mm 3 , 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
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 opensource 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 pixel 2 , 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 pixel 2 , 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
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 5 th and 95 th percentiles for the permuted data, at each target nest size, and compare the observed data with the bounds of that envelope.