Introduction
DNA methylation is a key epigenetic mechanism whereby a methyl group is added to the 5′ position of a cytosine pyrimidine ring at promoter regions of a gene (1). In humans, DNA methylation occurs mainly in the dinucleotide CpG sites (2). Aberrant DNA methylation has been identified as one of the important factors that contributes to cancer due to silencing of the tumour suppressor genes (3). The mechanisms underlying the silencing effects of DNA methylation are associated with chromatin configuration (4). When the CpG sites remain unmethylated, DNA binding proteins and transcription factors can easily access to the promoter site and active gene expression will occur (5,6). Methylation of CpG sites results in changes in the chromatin structure that interfere with the binding capacity of the standard machinery for transcription factors resulting in gene silencing (7,8). In general, there are two types of DNA methylation including promoter hypermethylation and global hypomethylation.
Promoter hypermethylation is associated with low expression or silencing of tumour suppressor genes, for example, the AT-rich interactive domain 1A (ARID1A) gene in breast cancer (9). Hypermethylation of the ARID1A promoter has resulted in transcription inactivation despite the changes in the gene copy number, mutations and histone modifications (10). Previous studies reported several hypermethylated genes as potential DNA methylation markers for breast cancer including TUSC5, DOK7, KLF11, SIM1, NT5E and OTP (11–15). However, there is a need to further explore the reliability, sensitivity and specificity of these biomarkers.
Focal hypermethylation can affect not only discrete genes but also a large region of the chromosome that leads to long-range epigenetic silencing (LRES) (16). In LRES, methylation of a particular gene may also suppress expression of neighboring unmethylated genes (17). This phenomenon was observed in two specific gene clusters that were related to breast cancer (16,18). In one of these examples, aberrant methylation in the CpG islands of HOXA1, HOXA7, HOXA5 and HOXA9 in chromosome 7 caused silencing of HOXA1 to HOXA10 genes. In the second report, hypermethylation of CpG islands in protocadherin (PCDH) in chromosome 5 resulted in low expression of the PCDH gene cluster (16,18). Global hypomethylation commonly occurs at repetitive elements such as satellite DNA sequences (19). An example of this mechanism was documented in a breast cancer study where there was hypomethylation of the satellite sequence of SATR-1 gene (20). A study using deep sequencing identified long-range hypomethylation with focal hypermethylation of nuclear-associated domains and these were likely to cause gene silencing in colorectal cancer (21).
Recent publications have used multiple genomic datasets to get a clearer picture of the genetic events in cancers including breast cancer (22–24). This is because a single dataset analysis gives limited information and does not fully reflect the actual events in the cells. In the present study, we performed an integrative analysis combining DNA methylation and gene expression profiling datasets to identify the significant genes and their biological pathways that can served as potential therapeutic targets for breast cancer.
Materials and methods
Clinical samples
Approval for the study was obtained from the Ethics Committee of Universiti Kebangsaan Malaysia (ref. UKM 1.5.3.5/244/SPP/UMBI-003-2012). Subject recruitment was carried out at the Universiti Kebangsaan Malaysia Medical Centre and Hospital Kuala Lumpur, Malaysia. Primary breast tumour and adjacent non-cancerous breast tissue samples were collected from 87 female patients after a written informed consent. The collected tissues were kept in liquid nitrogen and stored at −80°C before further analysis. All tissues were sectioned at 5–7 μm thickness using a cryostat (Microtome Cryostat HM550; Microm International GmbH, Walldorf, Germany) and stained with haematoxylin and eosin (H&E) before being viewed and confirmed by a histopathologist. Tissues that contained > 80% of malignant cells were included. We only used non-cancerous tissues which were free from malignant or inflammatory cells and contained mainly ductal and lobular cells. Laser capture microdissection (Arcturus Engineering, Mountain View, CA, USA) was used for tissues that contained <80% of cancer or non-cancerous cells, and this followed the previously described procedure (25).
DNA methylation profiling
Genomic DNA was isolated using two different kits, the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) and the QIAamp® DNA Micro kit (Qiagen), according to the protocol of the manufacturers. DNA integrity was assessed by 1.0% agarose gel electrophoresis under 60 V for 1 h. We quantified concentration and purity of the extracted DNA using NanoDrop (Thermo Fisher Scientific, Leicester, UK).
Bisulphite conversion of genomic DNA was carried out using EZ DNA Methylation-Gold kit (Zymo Research, Irvine, CA, USA). DNA methylation analysis of 76 tumours and 25 adjacent normal breast samples was performed using the Illumina’s Infinium HumanMethylation27 Beadchip kit (Illumina, Inc., San Diego, CA, USA) based on the manufacturer’s protocol and followed the previously described procedure (26). The data generated by GenomeStudio was further analysed using Partek Genomics Suite 6.6 (Partek Inc., St. Louis, MO, USA). A 4-way analysis of variance (ANOVA) was used to compare CpG loci methylation data across tumour and normal groups. The different sources of variation in the entire data in this ANOVA model included sample group (tumour/normal), batch effect, status of oestrogen receptor (positive/negative) and sources of tissue type (frozen section/laser capture microdissection). Adjustment for the different sources of variation such as batch effect, status of oestrogen receptor and sources of tissue type was done. We further used the filtering characteristics of fold-change -2 to 2 and a false discovery rate (FDR) at P<0.05 to identify the differentially methylated genes.
Gene expression profiling
Total RNA was extracted using the RNeasy Plus Mini kit (Qiagen) according to the manufacturer’s protocol. RNA was quantified using NanoDrop (Thermo Fisher Scientific) and Agilent RNA 6000 Nano kit (Agilent Technologies GmbH, Waldbronn, Germany). RNA samples with OD 260/280 of 1.8–2.2 and RNA integrity number ≥6.5 were included in the present study. Gene expression profiling of 15 tumours and 5 adjacent normal breast tissues was performed using GeneChip® Human Gene 1.0 ST array (Affymetrix, Santa Clara, CA, USA). Data were extracted using the Affymetrix® Genotyping Console™ (Affymetrix) and were further analysed using Partek Genomics Suite 6.6 (Partek Inc.). Data were normalised using quantile normalisation and robust multi-array analysis (RMA) background correction. Filtering characteristics of fold-change -1.5 to 1.5 and a FDR at P<0.05 were used in identifying the differentially expressed genes.
Integrative genomic and epigenomic analysis
In order to identify overlapping genes across different datasets, gene symbols from each of the datasets were used. Filtered datasets in either in CSV (comma separated values) or tab delimited format with significant cut-off values from each of the datasets were import into MySQL relational database for downstream data analysis. Each of the datasets was compared in pairwise (gene expression vs. methylation). Unique gene symbol found between the overlapping comparisons were used as a gene list for downstream analysis. These overlapping genes were then analysed using KEGG pathway and DAVID v6.7 for functional annotation, classification and enrichment analysis. Functional classification and signalling pathway that showed P≤0.05 was considered significantly enriched. Expression values and methylation values were also extracted from the datasets for circular map generation.
Methylation-specific multiplex ligation-dependent probe amplification (MS-MLPA)
MS-MLPA was performed to confirm the methylation microarray results. The SALSA MLPA Kit P200-A1 Reference-1 (MRC-Holland, Amsterdam, The Netherlands) was used to detect aberrant methylation for genes GPX7, SPARC, TIMP3, BTG4 and SFRP2 in 70 tumours and 23 normal samples. Briefly, 5 μl of 50 ng/μl DNA was denatured at 95°C for 30 min followed by cooling down to 25°C and adding of the probes mix. The sample was then incubated at 60°C for 18 h. The sample was then divided into 2 tubes, one used as a control without the HhaI enzyme while another tube contained the HhaI enzyme (Promega, Madison, WI, USA) for the digestion. PCR was performed on all samples using the thermal cycler (Applied Biosystems, Foster City, CA, USA). The amplified products were further subjected to fragment analysis using the ABI 3500 Genetic Analyzer (Applied Biosystems). Data analysis was performed using the Coffalyser version 1.0.0.43 (MRC-Holland). Quantification of methylation status for each gene was obtained by comparing the probes relative peak area ratio from the digested samples with those obtained from the undigested samples. Digested samples with probes of relative peak area ratio ≥0.25 were considered as methylated.
Methylation-specific quantitative polymerase chain reaction (MS-qPCR)
MS-qPCR was performed to further confirm the methylation microarray results. Bisulphite conversion of genomic DNA was carried out using BisulFlash DNA modification kit (Epigentek Group Inc., New York, NY, USA) followed the manufacturer’s protocol. Methylation specific qPCR Fast kit (Epigentek Group Inc.) was used to detect aberrant methylation for genes SFRP1 and NRG1 in 50 samples. Basically, there were two types of primers in this assay, namely methylated primer and unmethylated primer. The methylated primer sequences for SFRP1 were 5′-GTTTTTAGTCGGATATCGGTTC-3′ (forward) and 3′-CACGTTATAACACAACCGCA-5′ (reverse) while the unmethylated primer sequences were 5′-GTGAGTTTTTAGTTGGATATTGGTTT-3′ (forward) and 3′-CCCACATTATAACACAACCACA-5′ (reverse). For gene NRG1, the methylated primer sequences were 5′-CGGATTGGGGTAAAATAAGTTC-3′ (forward) and 3′-ACAATAATAACAACAACGACAACGA-5′ (reverse) while the unmethylated primer sequences were 5′-AGAGTTGGGTAGAGTTTGAATTGA-3′ (forward) and 3′-CAACAATAATAACAACAACAACAACAAC-5′ (reverse). We used β-actin as the positive control in this assay. The MS-qPCR was carried out using Applied Biosystem 7500 Fast Real-Time PCR system (Applied Biosystems) followed the manufacturer’s protocol. Data analysis for gene SFRP1 was based on 43 tumour and 7 normal breast tissues. For gene NRG1, data analysis was performed based on 47 tumour and 2 normal breast tissues as one sample had been identified as an outlier and excluded from the data analysis. Data generated by MS-qPCR were further analysed by using Microsoft Excel. Ct value of each sample was normalised with the Ct value of the β-actin. Percentage of methylation level for each sample was calculated based on a previous study (27). Unpaired t-test was used to test the significance of the results.
Results
Novel gene clusters are hypermethylated in breast cancer
DNA methylation analysis was conducted on 87 breast cancer patients with a mean age of 55.90±11.17 years. The epidemiological data of the patients are shown in Table I. Principal component analysis (PCA) showed the tumour and normal samples were clustered distinctly (Fig. 1). Supervised hierarchical clustering distinctly separated the tumour group from the adjacent normal breast tissues. Filtering using a fold-change −2 to 2 and a FDR P<0.05 generated 1,411 CpG sites which cover 1,049 genes. These significantly methylated CpG sites were further grouped into 1,389 sites with high methylation level or hypermethylation and 22 with low methylation level or hypomethylation (Fig. 1). We excluded 133 false CpG sites from the pathway analysis and 8 CpG sites which were located on the X chromosome. There were 9 clusters of genes that were hypermethylated including two which (HOXA and PCDH) were previously reported in breast cancer (Table II) (16,18).
Gene expression pattern in breast cancer through genome-wide expression microarray
Gene expression microarray analysis was performed on 15 tumour and 5 normal methylation-matched samples. Filtering using a fold-change of −1.5 to 1.5 and a FDR P-value of <0.05 identified 867 differentially expressed genes. Principal component analysis (PCA) showed the tumour and normal samples were clustered distinctly (Fig. 2). The supervised hierarchical clustering revealed 404 upregulated and 463 downregulated genes in cancer compared to non-cancerous samples (Fig. 2). The top 10 upregulated genes were CASC5, CENPF, KIF23, DTL, MK167, TPX2, NUF2, KIF4A, NUSAP1 and BUB1B while the top 10 downregulated genes were PAK3, B3GALT1, CX3CL1, EDN3, KCNMB1, HOXA5, NRG1, KLHL13, TSHZ2 and IL17RD. Gene Ontology enrichment analysis revealed that most of the genes were enriched in cell proliferation, viral reproduction, pigmentation, growth, rhythmic process, cell killing and metabolic process under the biological process. For the molecular function, most of the genes were enriched in the chemoattractant, structural molecule, translation regular, enzyme regulator, transporter and binding activity. For the cellular component, most of the genes were active in extracellular region and synapse. For the gene set enrichment analysis (GSEA) with P<0.05, most of the genes are involved in cell proliferation, spindle, M phase of mitotic cell cycle, microtubule-based movement, microtubule motor activity, condensed chromosome kinetochore and microtube.
Integrative analysis showed 64 driver genes involved in breast cancer
The integrative analysis revealed 64 overlapping significant genes between DNA methylation and gene expression analysis. Notably, all of the overlapped genes were hypermethylated with 51 showing negative association and 13 positive association (Table III). The 64 overlapping genes were further mapped to the KEGG database as shown in Table IV.
For the pathway analysis, two of the enriched pathways identified were the focal adhesion and extracellular matrix-receptor interaction. Seven genes (PAK7, COL4A1, CCND2, COL1A2, RELN, COL1A1 and THBS2) are involved in the focal adhesion while five genes (COL4A1, COL1A2, RELN, COL1A1 and THBS2) are involved in extracellular matrix-receptor interaction.
Gene Ontology (GO) enrichment analysis was carried out on the overlapping genes. The enrichment analysis showed that for the biological process, most of the genes were enriched in the skeletal system development. For the molecular function, most of the genes were enriched in the platelet-derived growth factor binding while for the cellular component, most of the genes were involved in collagen type I (Fig. 3). The circular map generated showed the overall visualisation of overlapping genes according to chromosomes (Fig. 4).
MS-MLPA and MS-qPCR confirmed the methylation profiling results
A total of 7 genes were selected for the validation of methylation microarray results. The genes were selected based on their association of breast cancer and the pathways of these genes involved. Five genes included GPX7, SPARC, TIMP3, BTG4 and SFRP2 were validated with MS-MPLA while another two genes, NRG1 and SFRP1, were validated using MS-qPCR. Both MS-MLPA and MS-qPRC results showed that the breast cancer samples have higher percentages of methylation compared to normal samples for all the selected genes (Figs. 5 and 6). The highest percentage of methylation was seen in TIMP3 followed by SPARC, SFRP2, BTG4 and GPX7 in MS-MLPA while MS-qPCR showed that the average methylation percentage in tumour samples for NRG1 and SFRP1 were 53.5 and 27.2%, respectively, with P-value <0.001.
Discussion
DNA methylation mechanism is a reversible process (28). Reversibility of DNA methylation can result in gene re-expression that leads to normal gene regulation (29). Based on this process, it can serve as a potential therapeutic target in cancer. Up to date, numerous DNA methylation studies have been carried out in breast cancer. However, most of the studies provided results that were obtained from single dataset analysis and gave limited information in terms of associating methylation to the transcriptome datasets. We have performed an integrative genomic analysis to identify the significant genes that can serve as important prognostic indicators and potential therapeutic targets for breast cancer. We found 9 gene clusters in breast cancer from the DNA methylation profiling analysis. Two of the gene clusters namely HOXA and PCDH were previously reported in breast cancer (16,18). These two gene clusters were involved in LRES which can cause multiple genes from the same cluster to become silenced. It is highly possible that the same mechanism could appear in the other 7 gene clusters that may lead to silencing of their nearby genes.
Integrative analysis identified 64 hypermethylated genes with 51 showing negative association and 13 with positive association. Out of the 51 hypermethylated genes with low gene expression, 14 genes were associated with breast cancer, including DAB2IP, NDRG2, AGTR1, CXCL2, CCND2, DKK3, FGF2, KLK10, NRG1, PTN, PTPRZ1, SFRP1, RELN and KIT (30–43). Five of these genes, DAB2IP, DKK3, KLK10, NGR1 and SFRP1, were tumour suppressor genes (41,44–47). Similar findings were documented in previous breast cancer studies that confirmed the reliability of our results (30,35,37,38,41,48–50).
DAB2IP is a Ras GTPase-activating tumour suppressor protein which plays an important role in maintaining cell homeostasis (51,52). Hypermethylation of this gene might disrupt the cell homeostasis and lead to breast carcinogenesis. The DKK3 encodes a protein that plays a role in inhibiting planar cell polarity pathway which regulates cell adhesion, motility and polarity (35,53). Therefore, inactivation of this gene by hypermethylation might activate planar cell polarity pathway and cause metastasis in breast cancer. The SFRP1 is the modulator of Wnt signaling pathway which plays a significant role in embryonic development, cell differentiation and proliferation (41,54). SFRP1 protein can inhibit the Wnt signalling pathway by binding to WNT1 molecules (50). It has been proposed that inactivation of SFRP1 is an early event in breast cancer and downregulation of this gene has also been shown to be associated with poor prognosis in breast cancer (41,49).
Examples of non-tumour suppressor genes which were hypermethylated with low expression and have crucial roles in breast carcinogenesis included NDRG2, CCND2 and FGF2. NDRG2 plays a role in stress responses, cell proliferation and differentiation (55). A previous study showed that this gene can regulate CD24 expression to decrease metastasis in breast cancer (56). Thus, hypermethylation that leads to silencing of NDRG2 might contribute to metastasis of breast cancer cells. CCND2 is involved in cell cycle regulation as it regulates the transition from G1 to S phase during cell cycle (34). This gene was found to be inactivated and hypermethylated in other studies (34,57). In addition, FGF2 is a growth factor that regulates epithelial cell proliferation, migration and angiogenesis (58). Loss of expression in FGF2 might probably lead to uncontrolled cell proliferation and metastasis of breast cancer.
Our integrative analysis approach showed that PTPRZ1, AGTR1, PTN and CXCL2 were hypermethylated and showed low expression. This observation contradicts with other studies on gene expression (32,33,39,40). However, those studies have no information on the methylation status of these genes.
We discovered 13 hypermethylated genes with high expression. Three genes, namely COL1A2, FBN1 and COL4A1, were previously reported in breast cancer (59–61). Previous studies showed that COL1A2 and FBN1 were hypermethylated and down regulated in the breast cancer cell lines (61,62). However, the present study showed that both of these genes were hypermethylated and upregulated in breast tumour tissues. The difference between the results might be due to the nature of the biological samples which gave distinct gene expression patterns (63). For gene COL4A1, its protein was reported to be elevated in the serum of primary breast cancer patients (60). This finding somewhat supported our results as COL4A1 was also upregulated. It has been suggested that the positive association between hypermethylation and gene expression could probably be explained by the presence of long non-coding RNA (64).
We further compared our results with the genes listed in the MammaPrint assay (Agendia, Amsterdam, The Netherlands), Oncotype DX assay (Genomic Health, Redwood City, CA, USA) and Ion Ampliseq™ Comprehensive Cancer Panel (Life Technologies, Carlsbad, CA, USA) to check for overlapping genes. There were only 5 genes (CASC5, CCND2, COL1A1, EPHB1 and KIT) found to overlap with the gene list in Ion Ampliseq™ Comprehensive Cancer Panel whereas none of our genes were in the gene list of the MammaPrint and Oncotype DX assays.
For the pathway analysis, two of the enriched pathways which have been identified are the focal adhesion and extracellular matrix-receptor interaction. There were 7 genes (PAK7, COL4A1, CCND2, COL1A2, RELN, COL1A1 and THBS2) involved in focal adhesion while 5 genes (COL4A1, COL1A2, RELN, COL1A1 and THBS2) involved in extracellular matrix-receptor interactions. PAK7 belongs to the PAK family of Ser/Thr protein kinases which are known to regulate cytoskeleton dynamics, proliferation and cell survival signalling (65) while RELN encodes a glycoprotein that acts as a regulator for neuronal migration (66). Previous study showed that RELN was inactivated by hypermethylation and silencing of RELN was associated with poor prognosis in breast cancer (43). Finally, THBS2 is involved in inhibiting angiogenesis (67). Aberrant methylation of these genes might disrupt the pathways involved and further contributes to breast cancer.
In conclusion, we have successfully performed the integrative genomic analysis from DNA methylation and gene expression profiling datasets and revealed a focused list of key genes in breast cancer. This list will further be used to study in more detail the pathogenesis of breast cancer using the interactive genomics data.