Multi-platform whole-genome microarray analyses refine the epigenetic signature of breast cancer metastasis with gene expression and copy number

Abstract

Background: We have previously identified genome-wide DNA methylation changes in a cell line model of breast cancer metastasis. These complex epigenetic changes that we observed, along with concurrent karyotype analyses, have led us to hypothesize that complex genomic alterations in cancer cells (deletions, translocations and ploidy) are superimposed over promoter-specific methylation events that are responsible for gene-specific expression changes observed in breast cancer metastasis.

Methodology/principal findings: We undertook simultaneous high-resolution, whole-genome analyses of MDA-MB-468GFP and MDA-MB-468GFP-LN human breast cancer cell lines (an isogenic, paired lymphatic metastasis cell line model) using Affymetrix gene expression (U133), promoter (1.0R), and SNP/CNV (SNP 6.0) microarray platforms to correlate data from gene expression, epigenetic (DNA methylation), and combination copy number variant/single nucleotide polymorphism microarrays. Using Partek Software and Ingenuity Pathway Analysis we integrated datasets from these three platforms and detected multiple hypomethylation and hypermethylation events. Many of these epigenetic alterations correlated with gene expression changes. In addition, gene dosage events correlated with the karyotypic differences observed between the cell lines and were reflected in specific promoter methylation patterns. Gene subsets were identified that correlated hyper (and hypo) methylation with the loss (or gain) of gene expression and in parallel, with gene dosage losses and gains, respectively. Individual gene targets from these subsets were also validated for their methylation, expression and copy number status, and susceptible gene pathways were identified that may indicate how selective advantage drives the processes of tumourigenesis and metastasis.

Conclusions/significance: Our approach allows more precisely profiling of functionally relevant epigenetic signatures that are associated with cancer progression and metastasis.

Figure 1. Experimental design for gene expression, promoter methylation and copy number analysis, and data integration.

Individual microarrays in replicates (red, light blue or gray boxes for expression, promoter methylation, or copy number variation analysis respectively) were imported into Partek Genomics Suite (PGS) software and background corrected using the RMA algorithm. Genes significantly altered in expression, promoter methylation, or in copy number were then compared using Venn analysis in PGS, and further validated. GFP: MDA-MB-468GFP cells; LN: MDA-MB-468GFP-LN cells; RMA: Robust Multichip Averaging.

Figure 2. Mapping of genomic DNA copy number variation to individual chromosomes.

The thicker vertical bar in the centre of each scan represent the normal diploid number, and the points represent smoothed averages of the probes on the array. Points falling to the right (or left) of the bar indicate regions of copy number gain (or loss, respectively) in 468LN vs 468GFP cells.

Figure 3. Chromosomal mapping of regions of copy number alteration.

Regions appearing increased in copy number are shown in red, and those decreasing in copy number in blue. A: 468GFP samples vs reference Yoruba population (YRI), B: 468LN samples vs YRI reference population and C: 468LN vs 468GFP.

Figure 4. Chromosomal mapping of copy number variations detected in Human SNP 6.0 arrays.

As a reference population, we used a subset (60 Yoruba females; YRI) of the 270 samples from the International HapMap Project run on the Affymetrix Human SNP 6.0 array . Copy numbers were normalized to 2 for ease of comparison. Chromosomal locations of regions of significant copy number alteration are shown: 468GFP vs the Yoruba reference population, 468LN vs Yoruba, and 468LN vs 468GFP.

Figure 5. Representative multi-array alignment of data from the SNP/copy number variation (top), gene expression (middle) and promoter methylation (bottom) are shown for chromosome 6.

In the upper panel 1 indicates relative copy number for 468GFP:468LN. For the U133 array data, plus represents significantly upregulated genes, and minus represents downregulated genes. For the promoter array data, plus represents significant regions of hypermethylation, and minus represents hypomethylated regions. Also, the promoter array data is presented on a log₂ scale, while for expression data, the height of the bars representing individual genes is proportional to the expression fold change. Similar multi-array data specific for all other chromosomes are provided in Figure S1.

Figure 6. Proportional Venn analysis of significantly changed gene regions as determined by multiarray analyses.

A: Venn analysis of genes predicted to be hypermethylated, decreased in expression, and showing a loss in copy number; specific regions of functional overlap are indicated (1 or 2). B: Venn analysis of genes predicted to be hypomethylated, increased in expression, and showing a gain in copy number; specific regions of functional overlap are indicated (3 or 4). The diameter of each circle is proportional to the number of genes identified by that specific array analysis.

Figure 7. Ingenuity Pathway analyses.

A: Top functional categories and B: canonical pathways from our data set based on significance.

Figure 8. Network analysis was performed to provide a graphical representation of genes having known biological relationships.

The EGFR and Mapk networks presented are shown in duplicate, with (A and C), displaying genes comprising the four insecting subregions shown in Figure 6 (regions 1,2,3 and 4) and (B and D) displaying genes comprising the two intersecting subregions (regions 1 and 3) that have a methylation/expression status that is independent of copy number. Green icons indicate downregulated genes and red indicates upregulated genes. The arrows indicate selected genes that have a variable methylation status that is dependent on copy number status.

Figure 9. Sodium bisulfite sequencing, gene expression and copy number.

Sodium bisulfite sequencing of representative genes detected with aberrant methylation with (or without) a concomitant change in copy number. Each square represents a CpG (open square: unmethylated; closed square: methylated). Each row of squares one cloned PCR sequence across the gene promoter (5–20 clones were sequenced per gene). Percentages indicate degree of methylation at each gene locus.

Figure 10. Expression and copy number analyses.

A: Quantitative real time RT-PCR expression data for each of these genes, including EGFR. Scale of the y-axis is log₁₀ of the fold change. B: Determination of copy number by quantitative real time PCR (qRT-PCR). Primers spanning an exon of the gene of interest were designed using ExonPrimer software, and qRT-PCR performed using 15 ng genomic DNA from 3 biological replicates each of LN and GFP DNA as template. Shown is the LN/GFP ratio, ± standard error of the mean (SEM). Data were normalized to β-globin as a reference gene. A p-value <0.05 indicates that mean normalized LN copy numbers in LN triplicates were significantly different from those in GFP.

Figure 11. Exposure to 5-aza-2-deoxycytidine and Trichostatin A.

A–C: 468LN cells were cultured for 7 days in A: the absence or B: the presence of 5-aza-2-deoxycytidine followed by C: an additional 16 h exposure to the histone deacetylase inhibitor Trichostatin A (TSA). D: Induction of epigenetically down-regulated genes. 468LN cells were cultured in the absence (Control) or presence of 5-aza-2-deoxycytidine (5AZA) for 72 hours (5AZA 72 hr), 88 hours (5AZA 88 hr), or for 72 hours followed by the addition of Trichostatin A (TSA) for 16 hours (5AZA+TSA). Total RNA was extracted and qRT-PCR performed as described in the text. Significant group differences were determined using ANOVA followed by the Student-Newman-Keuls multiple comparison procedure. a: significantly different (p<0.05) vs Control. b: significantly different vs 5AZA 72 hr. c: significantly different (p<0.05) vs 5AZA 88 hr.

Figure 12. Four scenarios (A–D) are presented to explain the relationship between copy number and DNA methylation.

Heat map analyses (Figure 3) revealed regional changes in copy number that reflected chromosomal aberrations in comparisons between the 468GFP and 468LN cells and the reference YRI population. These changes can be interpreted in the context of net gains (A,B) or losses (C,D) in DNA methylation. Specific chromosomal regions displaying these additive events are presented in the final column of this Figure.