Whole Genome Expression Microarray Analysis of Highly Versus Poorly Tumorigenic Murine Melanoma Cell Lines Provides Insights into Factors That Regulate Tumor Growth, Metastasis, and Immunogenicity.

Introduction

Melanoma, a cancer derived from pigment-producing melanocytes of the skin, is a disease of major public health significance. Although it is the least common form of skin cancer, melanoma is by far the most lethal due to its propensity to metastasize to a number of vital organs, including the brain, lungs, liver, and other visceral organs (1). According to the American Cancer Society’s 2015 Cancer Facts & Figures statistics, melanoma accounts for <2% of all skin cancer cases but is responsible for nearly 75% of skin cancer-related deaths, and it is predicted that melanoma will cause 9,940 deaths in the U.S. this year (2). Worldwide estimates of mortality for melanoma are that this disease is responsible for over 65,000 deaths annually (3). Moreover, recent data collected by the Surveillance, Epidemiology, and End Results (SEER) Program show that melanoma incidence rates have continued to increase over the previous 40 years (4). In the U.S. alone, the current annual costs for treatment and productivity losses associated with melanoma are near $3.3 billion (5). These numbers are even more staggering when considering a recent study by the International Agency for Research on Cancer, which ranked North America only fourth in melanoma incidence worldwide (6), thus highlighting the need to address melanoma as a global public health concern. With these statistics in mind, it is vital that we increase our understanding of the progression of melanoma so that we may improve current – and develop new – therapies for the treatment of this cancer.

The B16 murine melanoma cell line is one of the most widely used models for studying the progression and treatment of melanoma in pre-clinical settings (7, 8). Characterized by its rapid outgrowth as both primary tumors and metastatic lesions, the B16-F1 subline has also been shown to elicit dysfunctional CD8+ T cell responses (9) that mimic those frequently observed in many melanoma patients (10–12), and soluble factors derived from this tumor have been shown to alter the function of dendritic cells (DC) as well (13–15). We have recently described a chemically mutated variant of B16 melanoma, D5.1G4, that is significantly less tumorigenic and immunosuppressive than its wild-type counterpart (13, 15), though the basis for these differences between B16-F1 and D5.1G4 are poorly understood. Therefore, we sought to gain insights into factors that regulate melanoma tumorigenicity by comparing whole genome expression profiles of these murine melanomas. The microarray dataset obtained from these experiments and described herein offers a great resource for investigators wishing to study melanoma growth and progression and will likely drive the design of future experiments aimed at understanding the role of both individual gene products and entire pathways in melanoma tumorigenesis.

Materials and Methods

Cell lines

D5.1G4 murine melanoma cells were a generous gift of Dr. Jerry Neiderkorn (University of Texas Southwestern Medical School). B16-F1 murine melanoma cells were a generous gift of Dr. Victor Engelhard (University of Virginia). These tumor cell lines were grown in RPMI-1640 medium (Thermo Scientific, Hudson, NH, USA) supplemented with 10% fetal bovine serum (Premium Select, Atlanta Biologicals, Norcross, GA, USA), 2 mM L-glutamine, 2 grams/liter glucose, 2 grams/liter sodium bicarbonate, 100 U/ml penicillin (American Type Culture Collection, Manassas, VA, USA), and 100 μg/ml streptomycin (ATCC). All cells were grown at 37°C in 5% CO2 and were passaged at 80–90% confluence.

Tumor challenge

To evaluate tumor outgrowth, 4 × 105 B16-F1 or D5.1G4 melanoma cells in 0.2 ml of endotoxin-free 1× PBS were injected either subcutaneously in the nape of the neck or intravenously via the lateral tail vein. Mice injected subcutaneously were monitored every day for tumor formation. Following tumor formation, tumor area was determined every 2–3 days using digital calipers to take perpendicular diameter measurements of the tumor. Mice were euthanized once tumor area reached >300 mm2 or tumors became necrotic. Mice challenged intravenously were euthanized 17 days post-tumor challenge, and lungs were harvested for counting metastatic lesions under a dissecting microscope. Mice were used between 8 and 12 weeks of age, and all experiments were performed in accordance with regulatory standards and guidelines approved by the Hampden-Sydney College Animal Care and Use Committee.

RNA isolation

Prior to RNA isolation, tumor cells (1 × 106 cells/well) were plated in 6-well tissue culture plates and incubated at 37°C in 5% CO2. After 24 h of culture, B16-F1 and D5.1G4 melanoma cells were collected by cell scraping and centrifuged at 1,500 rpm for 5 min. Cell pellets were lysed with Buffer RLT containing β-ME, and RNA was extracted using an RNeasy Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s recommendations. Following RNA extraction, samples were treated with Amplification Grade DNase I (Life Technologies, Grand Island, NY, USA) to digest any residual DNA, and RNA was quantified using an Epoch Spectrophotometer (BioTek, Winooski, VT, USA). A260/280 ratios for all samples were >2.0.

Whole genome expression microarray

RNA samples (replicates of 3 for each melanoma cell line) were diluted to 100 ng/μl, and ~5 μg of RNA was shipped overnight on dry ice to Arraystar, Inc. (Rockville, MD, USA) for analysis using the company’s gene expression service and Agilent Mouse 4 × 44K Gene Expression Microarray v2 (Agilent Design ID 026655). RNA integrity was assessed by standard denaturing agarose gel electrophoresis. Sample labeling was performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology). Briefly, total RNA from each sample was linearly amplified and labeled with Cy3-UTP. The labeled cRNAs were purified using the RNeasy Mini Kit (Qiagen). The concentration and specific activity of the labeled cRNAs (pmol Cy3/μg cRNA) were measured using a NanoDrop ND-1000. One microgram of each labeled cRNA was fragmented by adding 11 μl 10× blocking agent and 2.2 μl of 25× fragmentation buffer, then heated at 60°C for 30 min, and finally 55 μl 2× GE hybridization buffer was added to dilute the labeled cRNA. One hundred microliters of hybridization solution were dispensed into the gasket slide and assembled to the gene expression microarray slide [Agilent Mouse 4 × 44K Gene Expression Microarray v2 (Agilent Design ID 026655)]. The slides were incubated for 17 h at 65°C in an Agilent hybridization oven. The hybridized arrays were washed and fixed prior to scanning. Slides were scanned on the Agilent DNA microarray scanner (G2505C) using the one-color scan setting for 1 × 44K array slides (scan area = 61 mm × 21.6 mm, scan resolution = 5 μm, dye channel set to Green, and Green PMT set to 100%).

Data processing and analysis

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software (Agilent Technologies). After quantile normalization of the raw data, genes for which at least three out of six samples had flags in Detected (“All Targets Value”) were chosen for further data analysis. Differentially expressed genes with statistical significance were identified through Volcano Plot filtering. D5.1G4 melanoma served as the reference sample, and genes expressed at greater than or equal to twofold higher or lower levels in B16-F1 (with a p value ≤0.05) were defined as genes exhibiting differential expression between the two tumors with statistical significance. Pathway analysis was performed using the latest Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the standard enrichment computation method.

Data deposition

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus (16) and are accessible through GEO Series accession number GSE69908 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69908). Data are freely available from this repository under file name GSE69908_RAW.tar, which contains .txt files of the normalized signal intensities for each Agilent Probe Name for all replicates of each sample.

Results

B16-F1 is a well-characterized murine melanoma cell line that models the growth, metastasis, and disruption of anti-tumor immunity often exhibited by aggressive melanomas in many patients with this disease. D5.1G4 is a chemically mutated variant of B16 that is significantly less tumorigenic than its wild-type counterpart, both with respect to its growth as primary tumors in the skin (Figure 1A) and with regard to its growth as metastatic lesions in the lungs (Figure 1B). In order to gain insight into factors that regulate melanoma growth and progression, RNA was isolated from these two melanoma cell lines for whole genome expression microarray analysis with the Agilent Mouse 4 × 44K Gene Expression Microarray v2 (Agilent Design ID 026655) that profiles the expression of 39,430 genes. Volcano Plot filtering was used to identify differentially expressed genes exhibiting greater than or equal to twofold upregulation or downregulation in B16-F1 melanoma as compared to D5.1G4 melanoma (Figure 1C). The dataset obtained from this analysis, accessible through the NCBI GEO data repository (16) at http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE69908, revealed statistically significant upregulation of 1,462 genes and statistically significant downregulation of 1,935 genes in B16-F1 melanoma. In order to better understand how these differences in gene expression might impact the biological properties of these melanoma cell lines, KEGG Pathway analysis was performed to identify biological pathways in which there was a significant enrichment of differentially expressed genes. In all, 71 KEGG pathways exhibited a statistically significant enrichment (p < 0.05) in upregulated genes in B16-F1, and 15 KEGG pathways exhibited a statistically significant enrichment in downregulated genes in this melanoma. The pathways with the 10 most significant enrichments scores for both upregulated and downregulated genes in B16-F1 are shown in Figure 1D, and the pathway-associated genes exhibiting a greater than or equal to twofold change in expression in this tumor versus the D5.1G4 melanoma are listed in Table 1.

Figure 1

Whole genome microarray analysis of B16-F1 and D5.1G4 murine melanomas. C57Bl/6 mice were challenged with 4e5 B16-F1 or D5.1G4 melanoma cells either subcutaneously (A) or intravenously (B). Subcutaneous tumor outgrowth was monitored over time, and lung metastatic lesions were enumerated 17 days post-tumor challenge. Data represent five mice per group and are graphed as the average with error bars designating standard error of the mean. Data shown are representative of three independent experiments. (C) RNA isolated from each respective melanoma cell line was used for whole genome microarray expression analysis as described in Materials and Methods. Volcano plot filtering was performed to identify genes that were differentially expressed with statistical significance between the B16-F1 and D5.1G4 melanomas. (D) KEGG pathway analysis was performed to identify biological pathways significantly enriched with differentially upregulated and downregulated genes in B16-F1.

Table 1

Differentially expressed genes associated with pathways with highest enrichment scores.

Pathway ID* (Mus musculus)	Differentially expressed genes associated with indicated KEGG pathway
B16-F1 upregulated pathways
Melanogenesis	Adcy7, Adcy9, Calm1, Creb3l1, Dvl2, Fzd2, Fzd6, Fzd7, Fzd9, Kit, Map2k1, Mapk1, Mc1r, Nras, Tcf7, Tcf7l1, Tyr, Wnt3a, Wnt9b
Rap1 signaling pathway	Adcy7, Adcy9, Adora2b, Akt2, Angpt2, Arap2, Calm1, Efna1, Efna4, Fgfr1, Itgb1, Kit, Lpar2, Map2k1, Mapk1, Mapk11, Nras, P2ry1, Pard6a, Pfn4, Pik3cd, Plcg1, Rac1, Rapgef3, Rasgrp2, Rasgrp3, Rgs14, Rhoa, Sipa1l1, Thbs1
Proteoglycans in cancer	Akt2, Col1a1, Ctsl, Ddx5, Eif4b, Fgfr1, Fzd2, Fzd6, Fzd7, Fzd9, Grb2, Hpse, Itgav, Itgb1, Kit, Lpar2, Map2k1, Mapk1, Mapk11, Nras, Nudt16l1, Pik3cd, Plcg1, Plcg2, Rac1, Rhoa, Tgfb1, Thbs1, Tlr2, Wnt3a, Wnt9b
Pathways in cancer	Akt2, Axin1, Cdk2, Col4a1, Col4a4, Dvl2, Egln1, Epas1, Fgfr1, Fos, Fzd2, Fzd6, Fzd7, Fzd9, Grb2, Itgav, Itgb1, Kit, Lama5, Lamb1, Map2k1, Mapk1, Msh3, Nras, Pik3cd, Plcg1, Plcg2, Pten, Rac1, Rhoa, Rxrg, Tcf7, Tcf7l1, Tgfb1, Traf4, Wnt3a, Wnt9b, Xiap
PI3K–Akt signaling pathway	Akt2, Angpt2, Atf2, Cdk2, Col1a1, Col4a1, Col4a4, Creb3l1, Efna1, Efna4, Eif4b, Fgfr1, Gnb4, Grb2, Itga1, Itga7, Itgav, Itgb1, Kit, Lama5, Lamb1, Lpar2, Map2k1, Mapk1, Nr4a1, Nras, Pck2, Pik3cd, Ppp2r1a, Ppp2r2b, Ppp2r3d, Prakaa1, Pten, Rac1, Thbs1, Them4, Tlr2, Ywhaz
N-glycan biosynthesis	Alg14, Alg5, Alg8, B4galt2, B4galt3, Fut8, Man1c1, Man2a1, Rft1, St6gal1, Stt3a
Estrogen signaling pathway	Adcy7, Adcy9, Akt2, Atf2, Calm1, Creb3l1, Fos, Grb2, Hspa8, Itpr3, Map2k1, Mapk1, Nras, Pik3cd, Prkcd, Shc4
Inositol phosphate metabolism	Aldh6a1, Inpp5b, Ippk, Isyna1, Minpp1, Pi4k2b, Pik3cd, Pikfyve, Pip4k2b, Plcg1, Plcg2, Pten
Colorectal cancer	Akt2, Axin1, Fos, Map2k1, Mapk1, Msh3, Pik3cd, Rac1, Rhoa, Tcf7, Tcf7l1, Tgfb1
VEGF signaling pathway	Akt2, Map2k1, Mapk1, Mapk11, Nras, Pik3cd, Plcg1, Plcg2, Ppp3cc, Rac1, Sphk1
B16-F1 downregulated pathways
Axon guidance	Ablim1, Dcc, Epha1, Epha2, Epha3, Epha4, Epha7, Ephb6, Gnai1, Plxna4, Plxnb3, Ptk2, Rgs3, Robo1, Sema3c, Sema3f, Sema4d, Sema4f, Sema6a, Sema6c, Slit2, Srgap1
Protein digestion and absorption	Coll11a2, Col17a1, Col18a1, Col1a2, Col2a1, Col4a2, Col4a4, Col9a1, Mme, Prcp, Slc7a9, Try4, Try5
Alanine, aspartate, and glutamate metabolism	Abat, Ass1, Ddo, Gfpt1, Gfpt2, Gls, Glul, Il4i1
Proteoglycans in cancer	Ank3, Araf, Camk2b, Camk2d, Cav2, Cd44, Col1a2, Fas, Fgf21, Fgf9, Fn1, Fzd8, Gpc1, Hbegf, Hspg2, Igf1r, Mapk13, Mras, Mtor, Ppp1r12a, Ptk2, Sdc4, Slc9a1, Wnt10a, Wnt5b, Wnt6
SNARE interactions in vesicular transport	Stx11, Stx16, Stx17, Stx3, Vamp1, Vamp5
NOD-like receptor signaling	Ccl2, Ccl5, Map3k7, Mapk13, Naip7, Nlrp1a, Nod1, Tnfaip3
Cocaine addiction	Cdk5r1, Creb3, Gnai1, Gnas, Grm2, Maoa, Slc18a1
ECM-receptor interaction	Cd44, Col11a2, Col1a2, Col2a1, Col4a2, Col4a4, Fn1, Gp9, Hspg2, Sdc4
Gap junction	Gja1, Gnai1, Gnas, Gucy1a3, Gucy1b2, Pdgfrb, Plcb4, Tuba1a, Tubb2b, Tubb3
Wnt signaling pathway	Apc, Camk2b, Camk2d, Fzd8, Map3k7, Nfatc1, Nkd1, Nkd2, Plcb4, Porcn, Vangl1, Wnt10a, Wnt5b, Wnt6

*Pathways with the 10 highest enrichment scores for differentially upregulated and downregulated genes in B16-F1 as compared to D5.1G4. Genes listed are those within the indicated pathways that exhibit a greater than or equal to twofold change in expression in B16-F1 with a .

Discussion

The gene expression microarray data described herein and available for further inspection and analysis at the publicly accessible NCBI GEO repository are a valuable tool for investigators wishing to study the roles played by individual genes and biological pathways in regulating melanoma tumorigenicity. These data are likely to be useful for those wishing to gain insights into melanoma-associated genes that influence various facets of tumor growth and progression, including cell cycle regulation, survival, cell adhesion and motility, tissue invasion and metastasis, and immune evasion. Our laboratory has been particularly interested in this latter phenomenon, as CD8+ T cell responses to well-established B16 melanoma have been shown to exhibit the type of dysfunction often associated with anti-tumor T cell responses in melanoma patients (9–12). Of note, our microarray data revealed a 16.77-fold upregulation (p < 0.0001) of Tgfb1 gene expression in B16-F1 as compared to its poorly tumorigenic D5.1G4 counterpart. It is therefore possible that melanoma-derived TGFβ1 suppresses CD8+ T cell effector function, leading to tumor immune escape and enhanced outgrowth. We have also recently shown that B16-F1-derived TGFβ1 alters that maturation and activation of tissue-resident DC and promotes their acquisition of a phenotype that is likely to favor tumor progression (15). Tumor-associated glycans have also been reported to promote immunosuppressive activity in DC (17), and KEGG pathways relating to glycan biosynthesis and cancer-associated proteoglycans were enriched with several genes exhibiting significant upregulation in B16-F1 (Table 1). It is interesting to speculate that overactivation of these and related immunoregulatory genes in B16-F1 may suppress cytotoxic T cell function or preclude efficient stimulation of CD8+ T cells by tumor-associated DC, thereby promoting tumor immune escape and enhanced outgrowth.

In addition to the role that overexpressed genes might play in promoting immune escape by B16-F1 melanoma, it is likely that downregulation of certain genes by B16-F1 may also contribute to escape from anti-tumor immune responses, even those that are robustly activated and exhibit significant effector function. For instance, our microarray data revealed significant downregulation of Cd72 (9.93-fold, p < 0.001) and Fas (9.91-fold, p < 0.01) gene expression in B16-F1 melanoma. CD72 functions as a cell adhesion and costimulatory ligand for the CD5 receptor on T cells (18), and its downregulation may limit high-affinity interactions between melanoma cells and tumor antigen-specific T cells. Likewise, melanoma-associated downregulation of Fas may promote tumor escape from activated tumor-specific T cells expressing FasL. Such findings might explain B16-F1 melanoma outgrowth in the face of the functional anti-tumor CD8+ T cell responses that are elicited against this tumor in early stages of its growth (9), and subsequent tumor progression could then be compounded by the overexpression of genes with active immunosuppressive functions as described above. It is also interesting that there is a significant enrichment of downregulated genes associated with the NOD-like receptor signaling pathway in B16-F1 melanoma. Downregulation of the Ccl2 and Ccl5 genes associated with this pathway may preclude efficient recruitment of various immune effectors into tumors or tumor-bearing tissue. Downregulation of other genes associated with this pathway, such as the Nod1 gene that promotes inflammation and Nlrp1 gene that functions to induce apoptosis, would also be expected to contribute to immune escape, enhanced survival, and the overall protumorigenic nature of this melanoma.

In conclusion, this comparative analysis of the whole genome expression profiles for the highly tumorigenic B16-F1 and poorly tumorigenic D5.1G4 melanomas offers a resource of data that provides significant insights into factors that regulate melanoma tumorigenicity. Those genes and pathways highlighted in this report are just a sample of the many that are likely to be of interest to tumor immunologists and cancer biologists, and it is hoped that the data obtained from this study will drive future investigation into the dysregulation of individual genes and biological pathways in melanoma that will enhance our understanding of this tumor’s progression and inform the design of new and improved therapies for the treatment of this cancer.

Author Contributions

KH was responsible for all aspects of this article and the work related to it.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This research was supported by funding from Virginia’s Commonwealth Health Research Board (Grant #375-01-14) and a Virginia Academy of Science Jeffress Research Grant from the Thomas F. Jeffress and Kate Miller Jeffress Memorial Trust (Grant #J-1018).