Dynamic Gene Expression Analysis Links Melanocyte Growth Arrest with Nevogenesis

Expression profiles often capture states rather than processes. Supervised techniques partitioning two states often recover a panel of expression differences that define a “signature” and that impute biological function. However, time series experiments that track dynamic gene expression may more effectively uncover expression “programs” that correlate with or dictate state transitions. This common approach has been undertaken in the genomic analysis of melanocytic tumors (Fig. 1A). Several microarray analyses have attempted to isolate expression signatures for subclasses of melanoma (M′; reviewed in ref. 1) and specific profiles for melanoma relative to benign nevi (M′ versus N′; refs. 2, 3). The M′ and N′ signatures represent a static view of a petrified sample, and thus very little can be deduced about the specific programs that led to the clinical lesions. As shown in Fig. 1B, a dynamic view of cultivated normal human melanocytes, or any primary human cells, readily reveals that a state transition occurs between an early, proliferative cell population (m₁) and a later growth-arrested cell population (m₂). The transcriptional program(s) that correlate with this melanocyte growth transition (m′) is distinct from classic clustering approaches for static lesions because it necessarily incorporates a sequential time domain. Moreover, fragments of the m′ program may still be captured in the M′ and N′ signatures thereby providing possible clues to the preclinical pathophysiology leading to melanomas and benign nevi. We thus set out in an unsupervised fashion to identify genes that may be captured by m′ programs and that may be potentially responsible for melanocyte growth arrest.

To discover programs of similar expression behavior over the course of melanocyte cultivation, we used a Bayesian method for model-based clustering of gene expression dynamics termed Cluster Analysis of Gene Expression Dynamics (CAGED; ref. 4). This method represents gene expression dynamics as autoregressive equations to search for the most probable set of clusters given the available data and is a specialized version of a more general class of methods called Bayesian clustering by dynamics (5). The main novelty of this approach is the concept of similarity in that two time series are comparable when they are generated by the same stochastic process. With this in mind, the Bayesian approach consists of searching the “most probable” set of processes generating the observed time series. Besides its ability to account for dynamic expression profiles, this method automatically identifies the number of clusters and partitions the gene expression time series in different groups on the basis of the principled measure of the posterior probability of the clustering model. Using CAGED, we were able to recover multiple clusters including several that paralleled the decelerating proliferation rate. We defined these clusters in aggregate as the melanocyte growth arrest program (MGAP) and show that it is still partially embedded in the signature of benign nevi, and that collagenous elements of the MGAP can still be shown in human melanocytic tumors.

For microarray analysis, we obtained normal human melanocytes from Dr. Mark Pittelkow (Mayo Clinic, Rochester, MN) as previously described (6). Melanocytes were all lightly pigmented and derived from foreskins of white neonates. Cells were grown in Medium 154 supplemented with human melanocyte growth supplement (Cascade Biologics) and penicillin/streptomycin. Melanocytes were grown for a total of 640 d (Fig. 1C) until a terminal doubling time of ∼4,994 d per doubling.

Total RNA from cells collected at days 6, 66, 122, 178, 234, 290, 346, 402, 458, 514, and 640 was isolated using the Qiagen RNAeasy kit and then submitted to the Harvard/Partners Center for Genetics and Genomics Microarray facility for analysis with the Affymetrix Human Genome U133 Plus 2.0 Chips. This microarray includes 1,300,000 unique oligonucleotides covering more than 47,000 transcripts, which in turn represent ∼39,000 of the best-characterized human genes. RNA purity was determined initially by 260/280 = 1.85–2.01 and finally by scanning with an Agilent 2100 Bioanalyzer using the RNA 6000 Nano LabChip (18S/28S ratio was ∼2.0).

The steps of cDNA synthesis, cDNA cleanup, and in vitro transcription and labeling of cRNA were carried out using standard Affymetrix kits. Hybridization was carried out according to the Affymetrix GeneChip Manual. Twenty micrograms of labeled, fragmented in vitro transcription material was the nominal amount used on the GeneChip arrays. A model 320 hybridization oven was used to incubate the arrays at 45°C overnight. Preparation of microarrays for scanning was carried out with Affymetrix wash protocols on a Model 450 Fluidics station under the control of Affymetrix GeneChip Operating Software.

Scanning was carried out on an Affymetrix Model 3000 scanner; percent present calls (ranges between 30% and 55%), background (<100), and spike-in control measurements were made using the GeneChip Operating Software algorithms and represented the fraction of detected transcripts of the total number of transcripts present on the array. Image scanning and data preprocessing to obtain transcript measurements were done using Affymetrix MAS 5.0 following the protocol recommended by the manufacturer.⁶

The analytic algorithm is outlined in ref. 4. The complete database consisted of 11 expression measurements of 54,673 genes. Data were normalized as ratios of the first value in each series. Only those genes with at least one observation of >3-fold increase or <3-fold decrease were considered in the analysis, for a total of 32,089 genes. The recorded expression values range between 1.33929E−4 and 274.1, with an average expression value of 1.952.

The analytic task is to partition the genes into clusters sharing a similar behavior over the 11 different conditions. A Bayesian approach assumes that the observed data were generated by a set of unknown processes. We do not know how many processes were responsible for the observed data, and the clustering function aims at identifying these processes. We defined two gene profiles as similar if they were generated by the same process. In a Bayesian framework, we can define this similarity measure as the posterior probability of a clustering model (i.e., a way of grouping individuals into groups) given the observed data. We started the clustering method adopted herein by assuming that all clustering models were equally probable (i.e., assuming a uniform prior distribution over these models), computed the posterior probability of each model given the data, and selected the most probable one. This approach spared us the effort of defining an arbitrary similarity threshold to decide whether two individuals are similar and had the further advantage of taking simultaneously all the data into account, rather than looking at single pairwise measures of similarity. Time series were merged into a cluster if the probability of the model with such a merge was at least seven times higher than the probability of the model without it.

We assumed the observations to be generated by a polynomial model of order 4, and set the prior precision and the γ value to 1 and 0, respectively. Polynomial models capture the pattern of dynamics along time. The prior precision was the size of the sample on which the prior distribution is built, whereas the γ value was the rate to zero of the prior precision, with 0 representing the case of perfect ignorance. The method requires a similarity measure to guide the search procedure and a Euclidean distance measure between gene profiles was adopted. Note that this measure was used only for search purposes and was not involved in the actual decision to merge a set of gene profiles together. Goodness of fit of the resulting model was assessed by checking the normality of the standardized residuals of each cluster. Data were natural logarithmic transformed.

Selected genes exhibiting differential expression in early and late samples were validated using TaqMan Gene Expression Assays (Applied Biosystems) as previously described (6). All target gene probes were labeled with FAM as the fluorescent reporter, whereas the endogenous control, GUSB, was labeled with VIC as the fluorescent reporter. Real-time PCR was done using the MX4000 (Stratagene), cycling for 2 min at 50°C, 10 min at 95°C, 1 min at 58°C to 66°C for 40 cycles. The normalized, relative levels of the genes between the early and late samples were expressed as the −log 2 ratio (i.e., 10 = 2⁻¹⁰ = 1,024-fold reduction in expression in early versus late sample).

The online interactive software is available at the company's website.⁷

The public data set used in this study (GDS1375) are available at the Gene Expression Omnibus (GEO) website⁸

For defining nevus-enriched genes (NEG) and melanoma-enriched genes (MEG), we first identified all genes that were differentially expressed between the 18 nevi and 45 melanoma specimens at the P = 0.05 level by a two-tailed test. Genes that were found to be expressed, on average, more than 2-fold in nevi were designated NEGs and vice versa for the MEGs.

Immunohistochemical staining was done on formalin-fixed, paraffin-embedded tissue with antihuman collagen III (BioGenex clone HWD1.1) at 1:100 following heat antigen retrieval in sodium citrate for 20 min at 95°C. Labeling was done with the Vectastain ABC-Kit/DAB (Vector Labs) and sections were counterstained with hematoxylin QS.

For this study, we used the Reticulum II staining Kit (Ventana Medical Systems). A total of eight cases (four nevi and four melanomas) were retrieved from the archives. Two of these melanomas have arisen in a background of benign nevus. Four-micron sections from formalin-fixed paraffin-embedded tissue were stained using the modified Gomori silver method. The samples were then analyzed by two of the authors using an Olympus BX51 microscope. Pictures were obtained using an Olympus DP25 digital camera.

Melanocytes were cultivated in culture over the course of 640 days. As shown in Fig. 1C, this represents ∼39 population doublings (gray boxes and line). Also depicted is the number of doublings per day, which is a population-based estimate of the proliferation rate (open circles, solid black line). As expected, the proliferation rate declined over the 2-year period as the cumulative doublings leveled off. At about 250 days, there was a slight inflection whereby an apparent surge in the proliferation rate occurred after a monotonic decline. The source of this “pause and restart” is unclear but does not seem to be related to extrinsic conditions such as media or supplemental growth factors. Fortuitously, as will be described, this biphasic growth kinetic served as an effective landmark during our analysis of expression dynamics.

The cluster analysis described in Materials and Methods yielded 10 clusters. Supplementary Fig. S1 illustrates the expression profiles of all members, and Supplementary Table S1 enumerates the identity of each member within the 10 clusters. Figure 2 shows the average expression profile of each of the 10 identified clusters; a cluster profile is the prototypical behavior of its members computed as point-wise average of the gene profiles within each cluster. The number of members of each cluster is also shown in parentheses. Supplementary Table S2 provides the descriptive statistics associated with each cluster group.

We can first use these expression levels to probe the composition and stability of the melanocyte culture over time. Direct visualization at the start of the experiment showed a pure population of stellate and bipolar cells consistent with melanocytes (Supplementary Fig. S2A). Moreover, the melanocyte-specific transcript, tyrosinase-related protein-1 (TYRP1), was the ninth most abundant transcript in the earliest sample after several housekeeping genes, suggesting that pigment cells likely represent the predominant cell type; TYRP1 remained the eighth most abundant transcript in the final sample. Overall, melanocyte-specific transcripts for TYRP1, tyrosinase, and Melan-A, in aggregate, accounted for 0.33% of the entire microarray signal intensity (compared 0.12% for β-actin) in the earliest sample and 0.31% of the microarray signal intensity (compared with 0.11% for β-actin) in the final sample. As shown in Supplementary Fig. S2B, transcript levels of TYRP1, tyrosinase, Melan-A, and vimentin all remained relatively constant throughout the 2-year period, suggesting that other nonmelanocytic cells neither dropped out nor emerged during the experiment.

Of the 32,089 transcript probes that were eligible for cluster assignment, 22,544 (70.2%) transcript probes partitioned into clusters 2 and 5, which showed minimal deviation across the entire period. This suggests that most genes do not significantly fluctuate with growth arrest physiology; furthermore, the stability of the expression levels argues against the presence of any individual outlier samples. In contrast, clusters 6, 8, 9, and 10 exhibited a striking parallel to the biphasic and diminishing proliferation rate (Figs. 1C and 2) although changes in gene expression levels preceded the physiologic changes (data not shown). We thus designated these genes as the MGAP (Fig. 3A). In total, these four clusters are composed of 659 uniquely annotated genes, of which 519 had single transcript representation whereas 140 had multiple transcript representations. Despite a wealth of biological information extant in the remaining clusters in Supplementary Table S1, we prioritized our pursuit of the MGAP genes based on levels of change, apparent physiologic relevance, and resource limitations.

To enrich for genes that may be more technically robust, we focused on genes that were represented by multiple probes within the MGAP clusters (MGAP* genes; Table 1). We characterized these genes first using standardized gene ontology (GO) designation. Of the 140 MGAP* genes, 133 mapped to at least one known GO feature (Supplementary Table S3). Within molecular GO functions, transcription factors (i.e., DNA binding) were the single largest class of proteins (ARNTL2, GLIS3, HMGA2, HOXA11, MKX, MSX2, PRRX1, RUNX1T1, SHOX2, SIM1, SIX4, SOX11, SOX9, TBX3, TBX5, TWIST2, ZEB1), whereas growth factors (BMP2, FGF1, FGF5, HGF, PDGFC, WNT5A, WNT5B) and receptors (EDNRA, EGFR, ITGA4, ITGBL1, LIFR, NRP1, PDGFRA, PLA2R1, TNFRSF11B, TNFRSF21) in aggregate were also frequently represented. Within cellular GO compartments, many of the MGAP* genes were identified as extracellular by its designation as a secreted factor or extracellular matrix (ECM) protein (43 of 133, 32.3%; ADAMTSL1, ASPN, BGN, BMP2, C1S, COL1A1, COL1A2, COL3A1, COL5A1, COL6A1, COL6A2, CPE, CYR61, DCN, DPT, EFEMP1, FBLN1, FGF1, FGF5, FN1, GLIPR1, GPC4, GREM1, GREM2, HGF, LOX, LTBP1, LTBP2, PAPPA, PCSK5, PDGFC, PLAU, POSTN, SERPINE1, SFRP2, STC1, STC2, TFPI2, THBS1, TNFRSF11B, VCAN, WNT5A, WNT5B). Taken together, the GO analysis suggests that, concomitant with proliferative slowdown, there is a reduction in the trophic milieu as modulated by declining autocrine and paracrine factor levels, growth receptor levels, and ECM constituents—all of which may be important in creating the microenvironment with tumors. To determine if any of our melanocyte program genes can still be recovered within the signature of melanocytic tumors, we performed an independent in silico survey of published nevus and melanoma profiles (3).

With the availability of large data sets on the GEO, it is now possible to interrogate biologically and clinically meaningful parameters to both generate and test hypotheses. One of the data sets available through GEO—GDS1375 (3)—used an Affymetrix U133 platform to profile 18 nevi and 45 melanoma specimens. We selected this panel given the similarities in microarray platform and the reasonable sample size. The in silico validation was done in three stages. First, of the 22,283 available GDS1375 elements, we identified all genes that were differentially expressed between nevi and melanoma at the P < 0.05 level (two-tailed test); this captured 10,681 total probes. Second, we mapped these probes to individual loci and defined two large clusters—a NEG set (mean nevi/melanoma expression ratio >2.0) and a MEG set (mean melanoma/nevi expression ratio >2.0); the NEGs and MEGs are shown in Supplementary Table S4. Third, we determined if MGAP* genes were overrepresented in either NEG or MEG signatures. Because the MGAP* gene expression levels fluctuate over time but the GEO-GDS1375 nevus and melanoma signatures are derived from a single static time point, we did not compare expression levels in this analysis but examined only representation.

Overall, there were 1,435 and 775 genes represented in the NEG and MEG sets, respectively. As shown in Supplementary Table S5, 40 of the 140 MGAP* genes (28.5%) were among the NEG set, whereas only 10 of 140 (7.1%) were represented in the MEG set (P = 0.02). When we considered all 659 MGAP genes, 144 (21.0%) overlapped with the NEGs, whereas only 33 (5.0%) with the MEGs (Fig. 3B); this association between the MGAP genes and the NEGs was highly significant (P < 0.001, Pearson's χ² test). Thus, fragments of the MGAP that were isolated from cultivated primary melanocytes seem to be retained preferentially in the transcriptional signature of benign nevi compared with the signature of malignant melanomas in vivo. This strongly suggests that MGAP genes may provide functional clues to melanocytic nevogenesis.

To determine if MGAP components coalesce into interactive pathways, we next subjected the MGAP* genes to the Ingenuity pathway analysis software. As shown in Table 1, there were statistically significant mappings to 13 canonical pathways using our panel of genes; the “hepatic fibrosis/hepatic stellate cell activation” pathway (Fig. 4A) showed the greatest significance (P = 5.129E-07). This pathway is complex and is modulated by multiple MGAP growth factor systems (e.g., TGF-β, PDGF, EGF, and FGF). To validate this mapping, we performed quantitative PCR and confirmed the dramatic decrease of many of the growth factor signaling and fibrosis pathway genes (Fig. 4B) as the melanocytes underwent growth arrest; in all cases, quantitative PCR paralleled the microarray studies. The collection of collagen and ECM genes found among the MGAP genes raises the possibility that this fibrosis pathway may play a role in pigment cell biology or melanocytic tumor formation. Because GDS1375 provides no pathologic detail, we turned to a set of human specimens for in vivo validation.

To recover evidence for a fibrogenic physiology among benign nevi and melanoma, we examined the protein levels of one of the MGAP genes, COL3A1 (type III collagen), in a collection of 21 melanocytic tumors (9 benign nevi, 3 dysplastic nevi, and 9 primary melanomas; Fig. 5). Overall, 9 of 9 benign nevi and 3 of 3 dysplastic nevi showed type III collagen either around the nests or around single cells; however, in stark contrast, only 3 of 9 melanoma samples showed any staining (P = 0.0015, two-tailed Fisher's exact test). The combined staining pattern is lost and becomes considerably weaker as lesions progress to dysplastic nevi and, ultimately, melanoma, consistent with the notion that type III collagen expression is part of a growth arrest program within melanocytes in vivo. Reticulum staining for an independent set of benign nevi and melanoma showed essentially the same results (Supplementary Fig. S3), suggesting that collagen III and other argyrophilic MGAP species are concentrated selectively around nevomelanocytes and not melanoma cells.

In this report, we (a) identified a MGAP from cultivated melanocytes, (b) showed that fragments of this transcriptional program can be preferentially recovered from benign nevi over melanoma, and (c) uncovered a fibrogenic pathway that partially defines this program and showed its existence in benign human nevi. In the parlance of Fig. 1A and B, the Bayesian CAGED method was successfully deployed to recover m′, which then revealed greater relevance to N′ and novel insights into nevogenesis. Although there were many clusters that contribute to distinct programs, we focused on the MGAP genes because it most closely mirrored the proliferation rate contour and thus most likely contained genes that relate to melanocyte growth. There were other clusters that we did not pursue in this study. For instance, cluster 7 (Supplementary Table S1: 114 probes) increased in levels as the proliferation rate decreased. Supplementary Table S6 shows that cluster 7 did not significantly map to any single pathway although cluster 7 components did map to several networks (Supplementary Table S7). Furthermore, although there were functionally intriguing genes within this group, no single gene was multiply represented like many of the elements in the MGAP genes.

Except for a brief recovery in the proliferation rate, melanocyte growth deceleration started early in culture and trailed off during the monitored period; the trigger for this attenuation is uncertain and under investigation. Many growth signaling molecules diminished over time (EDNRA, BMP2, FGF1, FGF5, WNT5A, WNT5B, SOX11, SOX9), suggesting the possibility that a trophic milieu may be locally shaped by young melanocytes in an autocrine manner and that depletion of these factors may in fact be a major component of growth cessation. It is also notable that many of these MGAP elements recapitulate key elements of pigment cell ontogeny. For instance, both endothelins (7) and Bmp2 (8) have been implicated in neural crest or melanocyte precursor development. Moreover, in model organisms, FGFs and Wnt proteins have been shown to regulate melanocyte and neural crest precursor development through Sox transcription factors (9). These developmentally vital factors may confer a more undifferentiated and proliferative phenotype, and as these stimuli decline, the cells enter a more differentiated and less divisive state.

Our results also point to a novel connection between melanocyte growth arrest, nevus formation, and fibrogenic signaling. Direct comparison of MGAP genes with NEGs and MEGs indicates that components of the MGAP can be preferentially detected within the nevus signature; thus, it is possible that the nevogenic history is partially fossilized in clinical specimens. This is biologically plausible because both nevi and primary melanocytes in vitro undergo a state transition from proliferation to stasis, whereas melanomas are in theory immortal. Moreover, nevomelanocytes are much less genetically deranged than melanoma cells and are thus more like their normal counterparts (10). Our findings also highlight an unexpected link between fibroplasia and nevogenesis. The MGAP panel contained a high number of collagen genes (COL10A1, COL13A1, COL14A1, COL1A1, COL1A2, COL20A1, COL21A1, COL3A1, COL5A1, COL6A1, COL6A2, COL6A3, COL7A1) that seem to be synthesized preferentially by early melanocytes in culture. If, in fact, the biology of these early melanocytes in vitro resembles that of early nevomelanocytes but not that of melanoma cells in vivo, then one would expect to find collagen fibers encapsulating individual nevus cells but not melanoma cells. Our type III collagen and silver stain results directly strengthen this biological association.

Earlier studies have reported that nevus cells produce basement membrane material (11, 12) including elastin (13) and type III collagen (14), which was clearly evident in our staining of benign and dysplastic nevi. Although the melanomas in our collection express lower levels of type collagen III and argyrophilic fibers, the role of these collagenous molecules in desmoplastic variants of melanoma remains unexplored. Taken together, it is conceivable that as nevomelanocytes synthesize basement membrane and collagenous proteins, the physical constraint of the dermis leads to local pericellular deposition and an encasement of material around the cells. The eventual growth arrest that defines the benign nevus is then accompanied by a drop in the production of these structural molecules.

So what are the potential regulators of this fibrogenic signal? The programmatic decline in the levels of TGFB1, LTBP1, and LTBP2 and a host of transcriptional targets downstream of TGF-β/Smad signaling points to an essential role for TGF-β in melanocyte aging. A well-known and important aspect of TGF-β function is regulation of fibrosis (15). Although it has not been established that all of the MGAP collagen genes are under the control of TGF-β, this function of TGF-β explains, in part, the mapping of the MGAP genes to the hepatic fibrosis pathway. It is conceivable that the drop in TGF-β synthesis heralds the decline in its target collagen genes. In humans, TGF-β1 is uniformly expressed in both benign nevi and melanoma (16) although its expression is thought to be lower in epidermal and cultured melanocytes (16, 17). Our data suggest that one of the reasons for this inconsistency may be the age of the melanocyte; younger, more proliferative melanocytes express higher levels of TGF-β1. In silico analysis of GDS1375 confirms the published data and also shows no difference in TGFB1 levels between benign nevi and melanoma (mean intensity, 434.7 versus 547.3; P = 0.167); however, a significant difference exists between normal melanocytes and benign nevi/melanoma as a group (mean intensity, 183.1 versus 515.1; P = 0.0039). Paradoxically, TGF-β1, derived from melanocytes, is known to inhibit proliferation (17) and to induce apoptosis (18–20)—an effect that is abrogated by FGF-2. Because many FGF isoforms, including FGF-2, are also members of the MGAP genes, it is conceivable that the balance of TGF-β and FGF dictates the growth phenotype of the target melanocytes and that, in fact, part of the growth decay may be due to autosuppression. Interestingly, melanoma cells are resistant to TGF-β-mediated suppression thereby explaining the high constitutive levels of TGF-β in actively dividing melanoma cells (19).

In summary, we have isolated a melanocyte arrest program that is partially recapitulated in the transcriptional signature of benign nevi. Bioinformatic analysis of genes within this program showed a large number of extracellular structural proteins, which includes collagens and which can be shown in human nevi specimens. These studies also show the power of dynamic expression analysis as an instructive means of recovering hidden programs in cellular physiology that may otherwise be overlooked in more static approaches to expression data.

No potential conflicts of interest were disclosed.

Grant support: Financial help from the generous philanthropic donors at MGH.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.