Cell Cycle Gene Networks Are Associated with Melanoma Prognosis

BACKGROUND: Our understanding of the molecular pathways that underlie melanoma remains incomplete. Although several published microarray studies of clinical melanomas have provided valuable information, we found only limited concordance between these studies. Therefore, we took an in vitro functional genomics approach to understand melanoma molecular pathways. METHODOLOGY/PRINCIPAL FINDINGS: Affymetrix microarray data were generated from A375 melanoma cells treated in vitro with siRNAs against 45 transcription factors and signaling molecules. Analysis of this data using unsupervised hierarchical clustering and Bayesian gene networks identified proliferation-association RNA clusters, which were co-ordinately expressed across the A375 cells and also across melanomas from patients. The abundance in metastatic melanomas of these cellular proliferation clusters and their putative upstream regulators was significantly associated with patient prognosis. An 8-gene classifier derived from gene network hub genes correctly classified the prognosis of 23/26 metastatic melanoma patients in a cross-validation study. Unlike the RNA clusters associated with cellular proliferation described above, co-ordinately expressed RNA clusters associated with immune response were clearly identified across melanoma tumours from patients but not across the siRNA-treated A375 cells, in which immune responses are not active. Three uncharacterised genes, which the gene networks predicted to be upstream of apoptosis- or cellular proliferation-associated RNAs, were found to significantly alter apoptosis and cell number when over-expressed in vitro. CONCLUSIONS/SIGNIFICANCE: This analysis identified co-expression of RNAs that encode functionally-related proteins, in particular, proliferation-associated RNA clusters that are linked to melanoma patient prognosis. Our analysis suggests that A375 cells in vitro may be valid models in which to study the gene expression modules that underlie some melanoma biological processes (e.g., proliferation) but not others (e.g., immune response). The gene expression modules identified here, and the RNAs predicted by Bayesian network inference to be upstream of these modules, are potential prognostic biomarkers and drug targets

The edges immediately downstream of gene network hubs that were enriched for cell cycle-associated children.

<p>Each red dot represents a probe set of either a gene network hub that has significant cell cycle enrichment of children, or one of the children. Blue lines represent gene network edges.</p

Laboratory investigation of gene network hubs.

<p>A–C, general cell biological effects of plasmid over-expression. Human 293T embryonic kidney cells (A and B) and human Mel501 melanoma cells (C) were transfected with <i>control</i> plasmids encoding lacZ and with plasmids encoding <i>UBE2S</i> (A), <i>ELMOD1</i> (B and C) and <i>TMCO1</i> (B and C). At 0, 2 and 4 days after the transfection, the number of viable cells was assessed using MTT assays. X-axes represent time in days while y-axis represent the OD570 absorbance (indicating viable cell number). Error bars represent standard deviation of the mean from four replicate wells. All graphs are representative of at least three independent experiments. D, Cell cycle analysis. 48 hours after transfection of plasmids into Mel501 cells, the cells were analysed by flow cytometry to identify the % cells in different phases of the cell cycle. Numbers show the percentage of hypodiploid cells. E, Fluorescent microscopy suggests that GFP-tagged over-expressed Elmod1 and TMC01 proteins have a punctate cytoplasmic distribution. F Western blotting indicates PARP cleavage in cells transfected by <i>Elmod1</i> and <i>TMC01</i> plasmids. 48 h after transfection of Mel501 melanoma cells with Elmod1 and TMC01 plasmids, protein lysates were analysed by Western blot using anti-β-actin and anti-PARP (a Caspase target degraded during apoptosis) antibodies. G–I, survival analysis in metastatic melanoma. Graphs compare survival of patients whose metastatic melanomas had above (green) or below (red) the 50^th percentile of the particular RNA expression in the Bogunovic 2009 data series <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0034247#pone.0034247-Bogunovic1" target="_blank">[57]</a>. All experimental data shown in panels A–F of this figure are representative of at least three independent experiments.</p

RNAs that are correlated in an E2F1-associated A375 cell-derived cluster are also correlated in datasets of primary and metastatic melanomas.

<p>A and B show heatmaps of Spearman's correlation coefficients between the E2F1 cluster probe sets shown in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0034247#pone-0034247-g003" target="_blank">Figure 3</a>, across primary and metastatic melanoma data. C and D show heatmaps of Spearman's correlation coefficients between members of a random list of probe sets the same size as the E2F cluster used in panels A and B, across primary and metastatic melanoma datasets. E and F show in red kernel density plots of the correlations shown in A and B, respectively. In grey they show Spearman's correlations across primary and metastatic melanoma datasets, respectively, between 10 random list of probe sets the same size as the E2F cluster. The colour key at the top left of each heatmap maps Spearman's correlation coefficient between probe sets to colour. The deepest red colour represents the Spearman's correlation coefficient of −1 and the deepest green colour Spearman's correlation coefficient of 1. Note that, in order to illustrate the broad range of correlations observed, the scale used here is different from that used in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0034247#pone-0034247-g003" target="_blank">Figure 3</a>.</p

Correlations within functionally-related gene sets in A375 cells and melanomas.

<p>All possible gene-gene correlations within gene sets defined by TRANSFAC and GO were calculated. The proportion of the gene pairs within each gene set that had Spearman's |ρ|≥0.5 was calculated separately for A375 siRNA microarray data (blue bars), a composite Affymetrix dataset from three primary melanoma studies (red bars), and a composite Affymetrix dataset from four metastatic melanoma studies (green bars). The x-axis represents TRANSFAC or GO gene sets. The y-axis represents the proportion of gene pairs from each gene set that had Spearman's |ρ|≥0.5 for the A375 data, primary tumour data, and metastatic tumour data, respectively.</p

Gene signatures from multiple microarray studies of melanomas from patients.

<p>The table summarises melanoma prognostic, progression and metastatic gene signatures that have been generated from a set of high-quality published microarray studies.</p>+<p>indicates that the array study generated the corresponding type of signature.</p

Frequency distribution of % ‘knock-down’ of target mRNA abundance.

<p>X-axis represents the percentage knock down in mRNA target abundance that was induced by the specific siRNA, as reported by microarrays (calculated based on mRNA abundance in the siRNA-targeted array/median mRNA abundance in all other arrays). The left y-axis represents the frequency of knockdown for each of the x-axis bins (blue bars) and the right y-axis represents the cumulative frequency (pink line).</p

Enrichment of transcription factor targets in A375 clusters.

<p>The first column is the cluster identifier. The second column is the enriched TRANSFAC Pro v8.2 transcription factor motif in the promoter of the genes in the cluster. The third column is the numbers of RNAs with the transcription factor target annotation in the cluster. The fourth and fifth columns are from the GATHER web tool – indicating the Bayes factor and the permutation <i>p</i>-value for the Bayes factor (indicating how often ≥ this Bayes Factor may be expected due to chance), respectively.</p

siRNA targets.

<p>‘OGS’ designates the official gene symbol of the target mRNA, ‘Degree of knockdown’ is the fold reduction in expression of the target RNA after siRNA incubation relative to median expression of the target RNA in all microarrays, ‘% reduction in expression’ is the % that the target RNA expression is reduced relative to median expression of the target RNA in all other microarrays, and ‘Batch’ is the experimental batch in which the siRNA was used.</p

Association between melanoma patient survival and the abundance of RNAs identified as gene network hubs that have their gene network children enriched for specific GO annotation categories.

<p>The x-axes represent GO annotations. GO annotations with related functions are coloured and grouped together: red (cell cycle related); sky blue (DNA repair related); green (RNA related); purple (transcription related); dark blue (metabolism related); yellow (protein related) and grey (miscellaneous). The left-most box (white) in each panel represents <u>all</u> RNAs on the microarray. The y-axes represents minus log_base2 p value, based on a Cox proportional hazards model, to indicate the strength of association between RNA abundance and patient survival. The horizontal line (y = 4.32) represents a threshold for a significant association with survival (equivalent to p = 0.05, above this line is a significant association). Within each box the dark horizontal line represents the median, the coloured area of the box the interquartile range. Panels A and B include all RNAs on the microarrays that have the specified GO annotations. Panels C and D include the Bayesian network hub RNAs, for which the function of downstream gene network children is significantly enriched for the specified GO annotations. Panels A and C represent primary tumours <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0034247#pone.0034247-Winnepenninckx1" target="_blank">[58]</a> while panels B and D represent metastatic tumours <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0034247#pone.0034247-Bogunovic1" target="_blank">[57]</a>.</p