|Home | About | Journals | Submit | Contact Us | Français|
Epigenomic and genomic changes affect gene expression and contribute to tumor development. The histone modifications trimethylated histone H3 lysine 4 (H3K4me3) and lysine 27 (H3K27me3) are epigenetic regulators associated to active and silenced genes, respectively and alterations of these modifications have been observed in cancer. Furthermore, genomic aberrations such as DNA copy number changes are common events in tumors. Pheochromocytoma is a rare endocrine tumor of the adrenal gland that mostly occurs sporadic with unknown epigenetic/genetic cause. The majority of cases are benign. Here we aimed to combine the genome-wide profiling of H3K4me3 and H3K27me3, obtained by the ChIP-chip methodology, and DNA copy number data with global gene expression examination in a malignant pheochromocytoma sample. The integrated analysis of the tumor expression levels, in relation to normal adrenal medulla, indicated that either histone modifications or chromosomal alterations, or both, have great impact on the expression of a substantial fraction of the genes in the investigated sample. Candidate tumor suppressor genes identified with decreased expression, a H3K27me3 mark and/or in regions of deletion were for instance TGIF1, DSC3, TNFRSF10B, RASSF2, HOXA9, PTPRE and CDH11. More genes were found with increased expression, a H3K4me3 mark, and/or in regions of gain. Potential oncogenes detected among those were GNAS, INSM1, DOK5, ETV1, RET, NTRK1, IGF2, and the H3K27 trimethylase gene EZH2. Our approach to associate histone methylations and DNA copy number changes to gene expression revealed apparent impact on global gene transcription, and enabled the identification of candidate tumor genes for further exploration.
Pheochromocytoma is a rare catecholamine-producing tumor that arises from chromaffin cells of the adrenal medulla or at extra-adrenal location (paraganglioma). The tumors mostly occur sporadically although the hereditary syndromes multiple endocrine neoplasia (MEN) 2A and 2B, von Hippel Lindau disease (VHL), neurofibromatosis type 1 (NF1) and paraganglioma syndromes (PGL) 1, 2 and 3 are all associated with the development of pheochromocytoma/paraganglioma. Mutations in the proto-oncogene RET or the tumor suppressor genes VHL, NF1, SDHD, SDHC, and SDHB predisposes to tumor development in these disorders (Xu et al., 1992; Latif et al., 1993; Mulligan et al., 1993; Baysal et al., 2000; Niemann and Muller 2000; Astuti et al., 2001a, 2001b). Approximately 10-25% of apparently sporadic tumors harbor mutation in one of these predisposing genes (Neumann et al., 2002; Bryant et al., 2003). However, the genetic/epigenetic cause(s) of the majority of apparently sporadic tumors is still unknown. Although most pheochromocytomas are benign, the prevalence of malignancy in sporadic adrenal pheochromocytomas is up to 10% (Bravo and Tagle 2003; Elder et al., 2005). The presence of metastases or invasive behaviour is the only criterium to determine malignancy. The diagnostic distinction between benign and malignant disease may however be difficult in the absence of metastases. Novel molecular markers with prognostic relevance are needed.
Loss of heterozygocity and metaphase/array-comparative genomic hybridization (CGH) analyses of both sporadic and familial pheochromocytomas have shown frequent losses on particularly chromosome arms 1p, 3p/q, 11p/q, 17p and 22q and gains on chromosome arms 1q, 12q, 17q, 19p/q and 20q (Bender et al., 2000; August et al., 2004; Cascon et al., 2005; Aarts et al., 2006; Jarbo et al., 2006; van Nederveen et al., 2009; Sandgren et al., 2010). DNA copy number changes in tumors may reflect an important cause of altered gene expression and by locating these changes one may identify individual genes that have implications in tumor development and progression. Global gene expression profiling of pheochromocytomas has mainly focused on, and revealed significant differences between familial and subsets of sporadic forms (Eisenhofer et al., 2004; Dahia et al., 2005) as well as between malignant and benign cases (Brouwers et al., 2006; Thouennon et al., 2007).
Epigenetic changes are important in cancer biology since they may affect gene transcription. Epigenetic profiling of pheochromocytoma has been restricted so far to investigations of promoter CpG-methylation of selected genes (Geli et al., 2007, 2008; Margetts et al., 2008). Thus, it is of considerable interest to characterize genome-wide the epigenetic mechanisms involved in benign and malignant pheochromocytoma development.
Posttranslational modifications of histones are additional major epigenetic transcriptional control mechanisms that establish domains of active and inactive chromatin. These modifications include methylation of lysine residues that can serve as either repressive or active marks. Specifically, trimethylation on lysine 27 (H3K27me3) and lysine 4 (H3K4me3) on histone H3 are associated with silent and active gene expression, respectively (Kirmizis et al., 2004; Bernstein et al., 2005; Barski et al., 2007; Mikkelsen et al., 2007; Pan et al., 2007). Several recent studies have reported global distribution patterns of these histone modifications in various cell types and their association to gene expression, mostly in embryonic stem (ES) cells of human and mouse origin, but also in human and mouse T cells and prostate cells (Barski et al., 2007; Mikkelsen et al., 2007; Pan et al., 2007; Zhao et al., 2007; Araki et al., 2009; Ke et al., 2009; Wei et al., 2009). However, global profiling of H3K27me3 and H3K4me3 and their suggested effects on gene expression has not been extensively studied in tumor tissue. H3K27me3-mediated tumor suppressor gene silencing as well as H3K4me3 and H3K27me3 enrichment associated to differential gene expression have been observed in prostate cancer cell lines (Kondo et al., 2008; Ke et al., 2009). The importance of this modification in tumor biology is further supported by reports that have demonstrated pre-marking by H3K27me3 of de novo-DNA methylated genes in a colon cancer cell line (Schlesinger et al., 2007) and frequent observations of polycomb protein EZH2 overexpression in several cancers (Simon and Lange 2008; Ke et al., 2009). The latter finding may be associated to the polycomb repression signature of metastatic prostate cells that was related to cancer outcome (Yu et al., 2007). EZH2 is part of the polycomb repressive complex 2 that specifically trimethylates H3K27 on target genes (Kirmizis et al., 2004). Furthermore, overexpression of the H3K4-specific methyltransferase SMYD3 has been observed in colorectal, hepatic and prostate cancer cells (Hamamoto et al., 2004; Ke et al., 2009), indicative of H3K4 hypermethylation at promoters of oncogenes. A role for H3K4me3 in tumorigenesis was also suggested by the finding that knockdown of SMYD3 by RNA interference inhibited cervical carcinoma cell growth and invasion in vitro (Wang et al., 2008a).
Here we have applied chromatin immunoprecipitation coupled to oligonucleotide arrays with 35 bp resolution (ChIP-chip), to genome-wide profile H3K4me3 and H3K27me3 occupancy in a malignant pheochromocytoma. Gene expression profiling in relation to normal adrenal medulla was also performed as well as array-CGH analysis to detect tumor-associated DNA copy number changes. Comprehensive genome wide array-CGH results for a large set of pheochromocytomas including the present sample were recently described (Sandgren et al., 2010). Epigenetic profiling in tumor tissues is preferred to tumor cell lines (Meissner et al., 2008) and so far global H3K4me3/H3K27me3 profiling has been done only in a limited number of cases using fresh or frozen tissues (Yu et al., 2007; Acevedo et al., 2008; Rada-Iglesias et al., 2009; Zhang et al., 2009). This study provides an integrated and detailed view of the genomic and epigenomic status of a malignant pheochromocytoma, and investigates the impact of copy number changes and histone marks on gene expression in the search for candidate tumor genes.
Using high-resolution whole-genome Affymetrix arrays we comprehensively examined the epigenetic signatures of the histone modifications H3K4me3 and H3K27me3 in the entire repeat free part of the human genome. Using stringent criteria, we identified 16096 and 27694 enriched regions of H3K4me3 and H3K27me3, respectively. We identified 2806 sites of overlapping H3K4me3 and H3K27me3 enrichment, i.e. bivalent regions. It should be noted that we could not, using the current method, distinguish between bi- or mono-allelic enrichment. Hence, some of the regions called bivalent could be the result of allele-specific histone modifications, i.e. one allele is enriched with H3K4me3 and the other with H3K27me3.
Consistent with previous reports (Barski et al., 2007; Guenther et al., 2007; Pan et al., 2007; Wei et al., 2009), the majority of H3K4me3 regions were located in close proximity of transcription start sites (TSSs), while H3K27me3 regions were more evenly distributed along the 20 kb region surrounding the TSS (Figure 1A). Moreover, the individual H3K27me3 and bivalent enrichments did, on average, occupy larger regions of the genome than H3K4me3 (Supplemental Data Table S1). Based on the observations of preferential locations of both marks and on the results of previous reports, we associated the histone marks to a conservative set of transcript predictions, UCSC Known Genes (KGs), if located within 2.5kb, upstream or downstream, of the corresponding TSS. As expected, a higher fraction of the H3K4me3 regions (~40%) were associated to KGs than H3K27me3 (~10%) or bivalently enriched regions (~15%) (Figure 1B). A number of these associations were confirmed by quantitative PCR (Supplemental Data Figure S1). By also considering potential TSSs as given by the locations of CAGE tag clusters (cap-analysis of gene expression) (Kawaji et al., 2006) the percentage of enriched regions for the two histone modifications at TSSs increased (Figure 1B), supporting the role of at least H3K4me3 as a frequent mark in the promoter regions of both known and novel transcripts. As expected, the average H3K27me3 signal around the TSSs covered a broader region and was also present at intergenic and intragenic locations further down in the gene body. Footprints of the H3K4me3 and H3K27me3 signals around associated KGs are shown in Figures 1C-1D.
KGs that were associated with non-overlapping enrichment of both H3K4me3 and H3K27me3, were termed semi-bivalent. Approximately 3.5 times more KGs were associated with H3K4me3 than with H3K27me3 (Table 1), excluding the semi-bivalent and bivalent KGs. As expected, we also found that fewer KGs were enriched with bivalent marks than what was observed in embryonic stem cells (Mikkelsen et al., 2007; Pan et al., 2007; Zhao et al., 2007). Gene Ontology (GO) enrichment analysis revealed that H3K4me3-associated KGs were significantly over-represented (P < 10-3; Fisher's exact test, Bonferroni adjusted) with housekeeping and metabolic processes (Figure 1E) in accordance with previous reports (Pan et al., 2007). H3K27me3-associated KGs were, on the other hand, significantly over-represented with developmental process terms (Figure 1F). KGs associated with bivalent marks were over-represented with multicellular organismal and anatomical structure development, both of which are shared with the H3K27me3 over-represented GO terms.
The tumor-specific DNA copy number alterations of the malignant pheochromocytoma analyzed here has been determined using a 32K BAC array. The genomic profiling revealed a complex pattern of chromosomal alterations including gains on chromosomes 5, 7, 12, 13, 15, 19 and 20 and losses on chromosome arms 11q and 16q as well as whole chromosome 18 (see below, Figures 3C and 3D). The total fraction of altered genomic sequence represented 29.4% of the genome. Chromosomal regions with gain were more common, affecting in total 769.7 Mb in contrast to 170.8 Mb with observed losses. The results of a comprehensive copy number analysis of pheochromocytomas in detail are presented elsewhere (Sandgren et al., 2010).
Examination of the autosomal located enriched regions of histone modifications in relation to the copy number data showed a similar distribution for H3K4me3, H3K27me3 and bivalency (Supplemental Data Figure S2). Relative to the fraction of altered genomic sequence, regions with increased copy number contained a high number of enriched regions. Coupling KGs to the histone modifications showed a similar pattern and indicated that the altered regions of copy number were more gene-dense (Supplemental Data Figure S3), thus explaining the lightly skewed distribution of histone modifications in relation to copy number.
To further elucidate the characteristics of the malignant pheochromocytoma, we measured autosomal gene expression levels in the tumor sample and in a normal adrenal medulla sample using Affymetrix U133A 2.0 plus arrays. The tumor mRNA expression levels of KGs marked with H3K4me3 and/or H3K27me3 and the expression levels of KGs located in chromosomal regions with increased copy number (gained), decreased copy number (deleted) or no copy number change (normal) are shown in Figure 2A. A clear trend of higher expression levels was evident for the group of KGs associated with H3K4me3-enrichment or gained regions in comparison with the expression values for KGs marked with H3K27me3 solely, bivalently associated, or in deleted regions. The average log2-transformed expression values for all KGs were similar to the expression values for KGs with no detected histone mark, i.e. neither H3K27me3 nor H3K4me3, and the expression values for KGs in diploid regions. For high-expressed KGs, we observed a significant over-representation of H3K4me3- associations (P < 10-177, Fisher's exact test) and gained regions (P < 10-60) while a significant under-representation of H3K27me3 (P < 10-50) and deleted regions (P < 10-25) when compared to all KGs measured on the array. For low-expressed KGs, the opposite was observed for H3K4me3 (under-representation, P < 10-172), gained regions (under-representation, P < 10-28), H3K27me3 (over-representation, P < 10-86) and deleted regions (over-representation, P < 10-27). Furthermore, nearly 70% of the high-expressed KGs were either associated with H3K4me3-enrichment or located within gained regions, indicating that these events are largely associated with high transcription levels. These results further strengthen the previously reported relationship of these two histone marks and their apparently opposite effect on transcription (Barski et al., 2007; Mikkelsen et al., 2007), as well as the impact of copy number changes on the expression level of individual genes.
When the expression levels of the malignant pheochromocytoma was related to the normal adrenal medulla sample, 3080 and 4290 differentially expressed KGs were identified (associated with 1043 and 1477 HGNC genes) with increased and decreased expression, respectively. The fold change of expression for all KGs with measured expression level and DNA copy number status grouped by their inferred associations are shown in Figure 2B. Evidently, the H3K27me3-marked KGs displayed lower expression in the tumor compared to the normal tissue, suggesting a transition of active histone modifications to the silent marker in some of these KGs. The H3K4me3-enriched KGs showed generally a higher fold change than H3K27me3-associated KGs. As expected, we also observed that KGs in deleted regions presented lower fold change and that KGs mapping to regions of gain showed higher fold change when compared to the median value of all KGs. Both observations are suggestive of an important impact on gene expression due to copy number changes. However, when we compared the fraction of KGs annotated with a certain histone modification in the increased and decreased groups relative to the fraction in the group of all measured KGs (Figures 2C-2H), we only found significant changes in the group of decreased KGs (Figures 2G and 2H). H3K4me3-enrichment was significantly under-represented (P < 10-115; Fisher's exact test, Bonferroni adjusted) while H3K27me3 was significantly over-represented (P < 10-74) in this group. A significant over-representation of bivalency (P < 10-8) was also detected in the group of decreased KGs. DNA copy number status seemed nevertheless to have a considerable impact on both groups. Gained regions were significantly over-represented in the group of increased KGs (Figure 2E; P < 10-90) and under-represented among KGs with decreased expression (Figure 2H; P < 10-35), while deleted regions were under-represented among those with increased expression (Figure 2F; P < 10-5) and over-represented among KGs with decreased expression (Figure 2G; P < 10-32).
The observed frequencies of KGs with altered DNA copy number status among the increased and decreased groups suggest that chromosomal changes may have great impact on gene expression levels, as reported in several studies of cancer genomes (Hughes et al., 2000; Pollack et al., 2002; Stransky et al., 2006). Moreover, the presence or absence of TSS-proximal epigenetic marks seems to be important for gene silencing.
The results from the above analyses showed that either histone modification or chromosomal alteration, or both, constituted the underlying cause responsible for the altered expression of a substantial fraction of genes in the investigated malignant pheochromocytoma. Their relationships are illustrated in Figures 3A-3B and the KG-associated genes with increased and decreased expression that were located in regions of altered DNA copy number and/or marked with H3K27me3 or H3K4me3 are listed in Supplemental Data Tables S2 and S3. We aimed to systematically search for genes related to the tumor by examining the groups of genes with altered expression, DNA copy number change and association to the H3K27me3 or H3K4me3 histone marks.
Firstly, genes with decreased expression that were associated with H3K27me3 and located in heterozygously deleted regions (Supplemental Data Table S2) constitute candidate tumor suppressors and are mechanistically interesting since it is possible that both their alleles are inactivated due to heterozygous deletion of one allele and a repressive histone modification, e.g. H3K27me3, at the other. In this group of 9 genes (associated with 17 KGs), we focus on 6 genes with potential tumor-related functions: ZNF423, TGIF1, AMICA1, PHLPP, ZBTZ7C and DSC3. The ZNF423 (16q12) gene product (zinc finger protein 423) is a transcription factor that can activate and repress target gene expression. It is involved in olfactorory neurogenesis, controls proliferation and differentiation of neural precursors in cerebellar vermins formation. ZNF423 is also a critical component required for retinoic acid induced differentiation. Low expression of ZNF423 has been associated with poor disease outcome in neuroblastoma patients (Huang et al., 2009). The TGIF1 gene (18p11.31) encodes a negative regulator of TGF-β signaling. AMICA1 (11q23.3) codes for an adhesion molecule and PHLPP (18q21.3) for a phosphatase that negatively regulates Akt and interestingly, functions as a tumor suppressor in colon cancer (Liu et al., 2009). ZBTB7C also maps to 18q 21.1, a locus commonly found deleted in many tumors. Reduced expression of ZBTB7C has been reported in cervical carcinoma cell lines (Reuter et al., 1998). The desmosome related gene DSC3 (18q12.1) is of interest, because lack of desmosomal adhesion has been implicated in conversion of early stage tumors to invasive cancers (Chidgey and Dawson 2007).
Secondly, a group of 174 genes were associated with decreased expression and marked with H3K27me3 (Supplemental Data Table S2). Several interesting candidate genes with possible tumor suppressor function were found, such as the TNFRSF10B (8p21.3) and HOXA9 (7p15.2) that both have been found to be CpG methylated in pheochromocytomas (Margetts et al., 2005, 2008). Promoter CpG methylation of HOXA9 has been reported in lung cancer (Rauch et al., 2007) and in neuroblastoma (Alaminos et al., 2004). The TNFRSF10B gene seems to play a role as a dosage-dependent tumor suppressor gene whose mono-allelic deletion can impair TRAIL-induced apoptosis in B-cell lymphoma (Rubio-Moscardo et al., 2005). It also induces both a caspase-dependent apoptotic pathway and the FADD pathway. RASSF2 (20p13) is an interesting candidate since it acts as a KRAS-specific effector-protein and may promote apoptosis/cell cycle arrest, and has been found hypermethylated in colorectal adenomas and gastric cancers (Hesson et al., 2005; Maruyama et al., 2008). Moreover, the PTPRE gene (10q26.2) is worth attention since it is a member of the protein tyrosine phosphatase family, which is known to regulate cellular processes such as cell growth and oncogenic transformation.
Thirdly, potential tumor suppressors were also found in the group of genes that showed decreased expression and were situated in regions of loss (124 genes, Supplemental Data Table S2), for instance the distinctly down-regulated gene WFDC1 (16q24)(fold change: -4.6) and the CDH11 gene (16q22.1). The WFDC1 gene is interesting since loss at this locus has been observed frequently in other types of tumors, including prostate, breast and hepatocellular carcinoma (Larsen et al., 2000; Watson et al., 2004). The CDH11 gene is a tumor suppressor candidate with decreased expression in osteosarcoma metastases (Nakajima et al., 2008). Moreover, LOH has been reported at this locus in retinoblastoma (Marchong et al., 2004).
Noteworthy, the pheochromocytoma-associated tumor suppressor gene VHL also showed decreased expression, however, no H3K27me3 enrichment or allelic loss were found at this genomic region. In addition, the CASP8 (2q33-2q34) and SERPINH1 (11q13.5) genes, that have been found CpG methylated in a subset of pheochromocytoma samples (Margetts et al., 2005, 2008), were down-regulated. Interestingly, the chromosome arm 11q is often deleted in pheochromocytomas, including the present sample (Figure 3D).
A relatively large group of 265 genes that displayed increased expression in the tumor were marked with H3K4me3 and were also located in regions of gain (Supplemental Data Table S3). Among these genes, we found several candidates that may promote tumorigenesis, for instance the GNAS gene (20q13.32) that is activated by mutations in kidney cancer (Kalfa et al., 2006) and the ETV1 gene (7p21.2) that is frequently overexpressed in prostate cancer (Cai et al., 2007). Furthermore, the INSM1 gene (20p11.23), which is overexpressed in endocrine tumors and small cell lung cancer (Goto et al., 1992; Lan et al., 1993), and the DOK5 gene (20q13.3) which product enhances Ret-dependent activation of the MAPK pathway (Grimm, et al., 2001) were identified in this group. RET (10p11.21) was also overexpressed in the malignant pheochromocytoma and marked with H3K4me3.
In the group of 562 genes that showed both increased expression and copy number (Supplemental Data Table S3), we found the genes NTRK1 (1q22) and PTTG1 (5q35.1). NTRK1 is a membrane-bound kinase receptor that, upon neurotrophin binding, phosphorylates itself and members of the MAPK pathway. Overexpression of NTRK1 has been observed in the related tumor form neuroblastoma (Nakagawara et al., 1992) and, for instance, in pancreatic cancer (Liu et al., 2007). It is located on chromosome arm 1q, which is often gained in pheochromocytomas (Sandgren et al., 2010), suggesting that this gene may be important in tumor development. The PTTG1 protein plays a central role in chromosome stability and in negative regulation of the p53 pathway.
Furthermore, among the 431 genes that had increased expression in the tumor and were marked with H3K4me3 (Supplemental data Table S3) we found IGF2 (11p15.5) and PLXNB1 (3p21.31), both with possible growth promoting functions. The IGF2 gene is imprinted and the regulated pattern of expression is maintained epigenetically in almost all tissues, but is frequently seen altered in cancers, where also subsequent overexpression is a recurrent theme (Chao and D'Amore 2008). PLXNB1 was evidently up-regulated in this malignant pheochromocytoma (fold change: 4.8) and has a role in axon guidance, invasive growth and cell migration, which makes it an interesting candidate oncogene.
Finally, among the genes that showed increased expression, and also in regions of gain and semi-bivalently marked, we found the polycomb associated trimethylase gene EZH2 (7q35-q36). Overexpression of EZH2 mRNA has been observed in several other cancers (Kondo et al., 2008; Simon and Lange 2008; Ke et al., 2009) suggestive of an important role for H3K27me3 in tumor development and as a potential biomarker. Increased expression of known H3K4 tri-methyltransferases (SMYD3, MLL3, SET) was not observed.
To investigate the biological relationship between genes with aberrant expression in the malignant pheochromocytoma, we analyzed gene interaction networks. We also aimed to examine the possible contribution of histone modifications and DNA copy number changes on the interacting genes.
The most enriched network for genes with increased expression in the tumor was related to cell morphology and cell transport (Figure 4A, Supplemental Data Table S5). This network included the oncogene RET, the RET- and MAPK-pathway associated gene DOK5 and the neuroendocrine tumor marker gene INSM1. To examine the possible cause of the increased expression of these genes, we overlaid the information about copy number status (Figure 4B) and associated histone modifications (Figure 4C). Along with several other genes in this network, both INSM1 and DOK5 are marked with H3K4me3 and in regions of gains, as also mentioned above. Noteworthy, many genes with increased expression also had direct or indirect interaction with the well-known cell signaling promoting gene ERK.
Interestingly, a gene network related to apoptosis was identified for the genes with decreased expression in the malignant pheochromocytoma (Figure 5A, Supplemental Data Table S6). This gene network indicated interactions between TNFRSF10B, FAD, CASP4, CASP8 and several other genes involved in apoptosis. When we overlaid the information about copy number changes (Figure 5B) and associated histone modifications (Figure 5C) to the genes in this network, it is illustrated that for instance, CASP4 was situated in a region of loss, whereas TNFRSF10B was marked with H3K27me3, as stated above.
In this study we present an integrative analysis of genome-wide histone modification status of H3K4me3 and H3K27me3, DNA copy number aberrations and gene expression changes in a malignant pheochromocytoma. By this combined approach we were able to pinpoint tumor biologically important genes, and also possible mechanisms behind the observed decrease or increase in expression for a substantial fraction of the genes.
Studies of genome-wide patterns of H3K4me3 and H3K27me3 in various cell lines have been performed recently (Pan et al., 2007; Zhao et al., 2007; Wang et al., 2008b; Araki et al., 2009; Ke et al., 2009; Wei et al., 2009), but the global, or locus-specific, patterns of H3K4me3- and H3K27me3-enrichment in pheochromocytoma have not been investigated before. In fact this has, to the best of our knowledge, rarely been done in any tumor or normal tissues with a method of comparable resolution (Acevedo et al., 2008; Rada-Iglesias et al., 2009; Zhang et al., 2009).
Our genome-wide profiling of H3K4me3 and H3K27me3 in the malignant pheochromocytoma sample showed a higher fraction of H3K4me3-regions associated with genes compared to H3K27me3 regions. This is in agreement with what has been reported in ES cells, T cells, and prostate cancer cell lines (Pan et al., 2007; Zhao et al., 2007; Araki et al., 2009; Ke et al., 2009; Wei et al., 2009), reinforcing the strong link between H3K4me3 and promoters (Barski et al., 2007; Guenther et al., 2007). Our observation of more and longer regions enriched for H3K27me3 than for H3K4me3 (Figure 1A) contrasts to human ES cells, where H3K4me3-regions are more frequent (Pan et al., 2007; Zhao et al., 2007). The higher frequency of H3K27me3-regions may be due to loss of H3K4me3 at several loci during differentiation, which has been observed for bivalent regions in differentied cells of mouse origin compared to ES cells (Mikkelsen et al., 2007). Also, intergenic regions may gain H3K27me3 modifications during differentiation in contrast to H3K4me3 enrichment, which is supported by the genome-wide report on substantially more H3K27me3 islands in six different T cell populations compared to H3K4me3 islands (Wei et al., 2009). Moreover, and important from a tumor perspective, prostate cancer cells have been found to contain more H3K27me3- than H3K4me3-enriched regions when compared to normal cells (Ke et al., 2009). In line with this and potentially interesting is the observed increase of EZH2 expression in the investigated tumor sample when compared to normal adrenal medulla. Up-regulation of this polycomb-associated methyltransferase, which specifically methylates H3K27, may be of significance in pheochromocytoma tumorigenesis, since it could contribute to silencing of important tumor suppressor genes. Overexpression of EZH2 has been observed in various cancers including the aforementioned prostate cancer cells (Simon and Lange 2008; Ke et al., 2009).
The finding from GO analysis that H3K27me3-associated KGs were significantly over-represented with developmental process terms have also been observed in ES cells, even though a major fraction of those H3K27me3-associated genes were also marked with H3K4me3, i.e. bivalent (Pan et al., 2007; Zhao et al., 2007). The enrichment of developmental related GO terms for H3K27me3-marked genes in differentiated cells is ambiguous. Our results are consistent with those of some cell types, i.e. colon mucosa and prostate cancer cells, but not with others, such as colon cancer cells and prostate cells, (Squazzo et al., 2006; Ke et al., 2009; Rada-Iglesias et al., 2009) in which various different GO terms were over-represented for genes with this repressive epigenetic mark.
Clearly, data on the genome-wide histone modification status of H3K4me3 and H3K27me3 in adrenal medulla from a healthy individual would have been beneficial for the study. Nevertheless, the combined analysis of global H3K4me3/H3K27me3 profiling, differential gene expression and copy number changes resulted in valuable new information regarding genes possibly important for pheochromocytoma development.
Firstly, as expected we observed a strong relationship between DNA copy number alterations, i.e. gains and losses, and mRNA expression levels for many genes. Specifically for the group of genes that showed increased expression, we observed a significant over-representation of gained genes and a significant under-representation of deleted genes. The contrary was observed for genes that displayed decreased expression, i.e. under-representation of gained genes and over-representation of deleted genes.
Secondly, in accordance with previous reports (Barski et al., 2007; Mikkelsen et al., 2007; Araki et al., 2009; Ke et al., 2009), we observed generally that high-expressed genes were associated with H3K4me3-enrichment while low-expressed genes were associated with H3K27me3. Moreover, we found that genes with decreased expression in relation to normal adrenal medulla were significantly related to the presence of H3K27me3- and absence of H3K4me3-enrichment. This was not the case for genes with relatively increased expression. Hence, several genes that were down-regulated in the tumor and associated with a H3K27me3 mark might have obtained this histone modification and, possibly, lost an active H3K4me3 mark thus enabling for reduced transcription. On the contrary for genes with increased expression in the tumor, acquisition of the H3K4me3 mark and possibly loosing H3K27me3 seemed less common.
The relationships between gain and loss of chromosomal material to increased and decreased gene expression, respectively, were recently reported in a liver cancer study. Also, comparisons between normal and tumor tissues regarding the global distribution of the repressive histone modifications H3K9me3 and H3K27me3, as well as RNA polymerase II occupancy were performed (Acevedo et al., 2008). It was concluded that copy number alterations were responsible for the majority of changes in gene expression in the tumor. We also observed that a larger fraction of genes seemed to be affected by DNA copy number than the absence or presence of H3K4me3 and H3K27me3. Presence or absence of these epigenetic marks and their strong relation to up- and down-regulated genes have recently been shown in prostate cancer cells and CD8+ T cells (Araki et al., 2009; Ke et al., 2009). Also in prostate cancer, H3K27me3-mediated silencing of tumor suppressor genes and associations of a polycomb epigenetic signature with adverse clinical outcomes has been reported (Yu et al., 2007; Kondo et al., 2008).
Our approach to focus on differentially expressed genes with characteristic histone modifications and/or DNA copy number status narrows down the list of tumor gene candidates. Relatively many genes were identified in regions of copy number gain with increased expression and marked with H3K4me3. Remarkably fewer genes showed decreased expression, copy number loss and H3K27me3 enrichment. Obviously, many other mechanisms can be responsible for decreased and increased expression that may explain the fraction of differentially expressed genes with no obvious genomic copy number imbalance and/or association to the histone modifications H3K27me3 and H3K4me3 (Figures 3A and 3B). Possible mechanisms include the presence or absence of other active epigenetic marks such as H3K36me, H3K79me, H3K9ac, H3K14ac, repressive modifications such as H3K9me and H4K20me (Wang et al., 2008b) and CpG-methylation. DNA point mutations and balanced structural rearrangements, for instance, may also affect gene transcription.
Nevertheless, of the nine heterozygously deleted and down-regulated genes with H3K27me3, we found interesting tumor associated functions for the majority (see above). However, of these, TGIF and DSC3 are worth more attention. TGIF inhibits TGF-β signaling that plays a fundamental role in cell regulation; it negatively controls cell proliferation, induces apoptosis and suppresses immune function. Disruption of this pathway occurs in several cancers (Liu, 2008). Loss of expression of adhesion molecule DSC3 is often observed in combination with reduced β-catenin expression in several cancers and is related to invasive behaviors (Wang et al., 2007). DSC3 has been reported also to be down-regulated and marked with H3K27me3 in prostate cancer cells as well as epigenetically silenced in breast cancer due to CpG methylation (Oshiro et al., 2005; Ke et al., 2009).
HOXA9 and TNFRSF10B, additional candidate genes with possible tumor suppressor functions, showed decreased expression and were marked with H3K27me3. Both genes have been reported to be CpG hypermethylated in pheochromocytomas (Margetts et al., 2005, 2008). Interestingly, in several reports EZH2 was shown to have functional connection to DNA methylation (Simon and Lange, 2008), and a link between pre-markings of H3K27me3 at de novo CpG-methylated promoters was for instance shown in colon cancer cells (Schlesinger et al., 2007). Several HOX genes, for instance HOXB13 and HOXA5, have been found DNA methylated in human cancers (Okuda et al., 2006; Piotrowski et al., 2006). HOXA9 promoter methylation has been identified in lung cancer (Rauch et al., 2007) and in neuroblastoma (Alaminos et al., 2004). However, the only HOX gene with decreased expression in the malignant pheochromocytoma was HOXA9. On the contrary, nine of 24 HOX genes showed increased expression. De-regulated expression of HOX genes is commonly observed in other solid tumors and cell lines (Abate-Shen, 2002).
The gene-interaction network tool identified several genes related to apoptosis that showed reduced tumor expression, such as the H3K27me3-marked gene TNFRSF10B, CASP8, CASP4, and several others. Noteworthy, CASP4 is located in the deleted region on chromosome 11q in the studied sample, a chromosome arm that frequently presents with losses in pheochromocytomas (Sandgren et al., 2010).
We found several candidate oncogenes among the genes located in regions of copy number gain that showed increased expression and H3K4me3-association. A few were included in the highly enriched gene network related to cell morphology and cell transport, for instance the INSM1 and DOK5 genes. Both genes are situated on chromosome arm 20q, a region commonly associated with gain in malignant pheochromocytomas (Sandgren et al., 2010). The former is overexpressed in endocrine tumors and small cell lung cancer (Goto et al., 1992; Lan et al., 1993). The DOK5 gene is of interest since the encoded protein is a substrate for the Ret receptor tyrosine kinase and enhances Ret-dependent activation of the MAPK pathway. RET was also found overexpressed in this tumor, and located in the same identified gene network. Mutations in the this well-known proto-oncogene predispose to pheochromocytoma development in the familial cancer syndrome MEN type 2, but mutations rarely occur in apparently sporadic cases, but overexpression has been observed (Dannenberg et al., 2003). The MAPK pathway is more prevalently activated in adrenergic sporadic and MEN-associated pheochromocytomas in contrast to VHL-associated and noradrenergic sporadic cases (Eisenhofer et al., 2004; Dahia et al., 2005). RET was also enriched with an H3K4me3 mark, along with a large fraction of the genes in the same gene network. A number of genes in this network were also associated with copy number gain, indicating again that both copy number changes and histone modifications of individual genes may affect the expression of functionally connected genes.
The integrated analysis of genomic and epigenomic statuses in tumor tissue is a promising approach in the quest for new insights of cancer development and progression. The association of histone modifications and DNA copy number changes to differentially expressed genes revealed apparent impact on global gene transcription and enabled the identification of interesting candidate tumor genes. Their contributions to tumor development and possibly as prognostic markers will be the focus of future studies.
The tumor tissue from adrenal medulla was obtained from a female patient diagnosed with malignant pheochromocytoma sample 39 in (Sandgren et al., 2010). The normal medulla sample used for expression analysis was obtained from the surrounding healthy adrenal medulla tissue from a patient suffering from a benign pheochromocytoma. The samples were surgically removed at Uppsala University Hospital, snap frozen and stored at -70. Examination and confirmation of tumor and normal cell content in the studied samples were carried out by an experienced pathologist. Samples were collected with informed consent from the patients and the local ethical committees approved the study.
The extracted DNA from the tumor sample was hybridized on a whole genome BAC array consisting of more than 32,000 BAC clones. Array construction is outlined in a previous report (Diaz de Stahl et al., 2008). DNA labeling, hybridization, washing and scanning were performed as previously described with minor modifications (Buckley et al., 2002). In short, 1,5 µg of genomic test and reference DNA were labeled using random primers (BioPrime Array CHG Genomic Labeling kit, Invitrogen, Carlsbad, CA), CY3CTP (PA53021, GE HealthCare) and CY5 (PA55021, GE HealthCare). Purified probe DNA was mixed with 150 µg human COT-1 DNA, vacuum evaporated and resuspended in 50 µl of hybridization solution (50% formamide, 5% SDS; 10% Dextran Sulfate, M.W. 500000, 2 × SSC). The mixture was then denaturated for 5 min at 95 and incubated for 2 h at 45 before being applied to the array and covered by a lifterslip. The hybridization at 45 lasted 20 h in slide chambers (Corning Inc. New York,). Washing of the slides was performed in 25% formamide, 2 × SSC, 0.1% SDS at 45 for 17 min followed by 10 min in 1 × PBS at room temperature and finally rinsed in ddH2O for a few seconds before immediately drying them using pressurized dust free air.
Image capture and detection of Cy5 and Cy3 intensities for each spot on the arrays were performed using GenPix4000B scanner (Axon Instruments Inc, union City, CA) and acquired data were analyzed with GenePixPro v6 image analysis software (Axon Instruments). Raw data files were uploaded to a LIMS database for storage, hosted by the Linnaeus centre for Bioinformatics (LCB, base.lcb.uu.se). Filtering, normalization and copy number profiling of the data were performed within the LCB Data Ware House (LCB-DWH, http://dw.lcb.uu.se) (Ameur et al., 2006). Oversaturated spots, spots with low Signal to Noise Ratio (SNR<3) in either channel and spots manually or automatically flagged as bad or absent in GenePixPro were filtered out. Print-tip loess normalization was applied to adjust for dye bias and/or spatial effects. The classification of clones as gained, deleted or balanced level was done using SMAP (Andersson et al., 2008), an open source software package available from Bioconductor (Gentleman, et al., 2004). Each individual clone was subsequently assigned a preliminary copy number class (CNC): gain, deleted or balanced. Single clones having a different CNC than neighboring clones were re-assigned to the CNC of the neighbors. Subsequent analyses were performed using R (R Development Core Team 2008).
Chromatin immunoprecipitation was performed as previously described (Rada-Iglesias et al., 2005) with some modifications. Briefly, approximately 40-60 mg of tissues were cut in small pieces, cross-linked with 0.37% formaldehyde for 10 min and resuspended in cell lysis buffer for 10 min on ice. Nuclei were resuspended in 1× RIPA buffer and kept on ice for another 10 min. Chromatin was sonicated to a size of 150-250 bp and pre-cleared by incubating with Protein G-agarose (Roche) for at least 1 h at 48 with slow rotation. At this step, a fraction of the pre-cleared chromatin was kept as input DNA and the rest was incubated with 10 µg antibody at 48 overnight, and 100 µl of protein G-agarose was used for each ChIP reaction. Protein G-agarose was washed four times with 1 RIPA buffer, once with ChIP washing buffer and once with 1 TE buffer. DNA-protein complexes were eluted, treated with RNase A (Amersham Biosciences) and incubated at 65 for 6 h in order to reverse cross-links. Proteins were degraded by Proteinase K (Amersham Biosciences), and DNA was extracted by phenol/chloroform/isoamyl alcohol extraction, purified and resuspended in water.
The antibodies against H3K4me3 (product number ab8580) and H3K27me3 (product number 07-449) were purchased from Abcam and Millipore, respectively. Two biological independent replicates were performed for each histone modification together with the corresponding input DNA, used as total genomic DNA reference.
DNA amplification, labeling and fragmentation were performed according to Affymetrix recommendations. Quantitative PCR, before and after amplification, were performed to detect, if present, introduced bias. The labeled samples were hybridized on Affymetrix GeneChip Human Tiling 2.0R Array sets, (seven arrays) with an average resolution of 35 bp. The recommendations of array hybridization, washing and staining conditions by Affymetrix were followed. Raw array data are available from ArrayExpress, accession: E-TABM-812.
Affymetrix Tiling Analysis SDK, revision 4, was used to quantile-normalize each replicate of ChIP measurements together with the input measurements. Probes were smoothed using a sliding window of 301 bp (bandwidth = 150) assigning the median value of the window. The resulting replicates of ChIP measurements over input measurements log2-ratios were quantile-normalized and the averaged. The latter and subsequent analyses were performed using R (R Development Core Team 2008).
Enriched regions were defined using a Z-score based sliding window approach. For each probe on the array, a Z-score was calculated from the average of log2-ratios of probes within a window of 150 base pairs centred on the probe against the distribution of log2-ratios on the array. Enriched regions were defined from enriched probes (Z-score >= 6) using the following criteria; (i) at least 2 enriched probes within the region, (ii) each enriched probe must be flanked on both sides by existing measurements within 110 base pairs, (iii) the maximum gap between enriched probes within a region is 110 base pairs. Regions were prolonged with probes positioned within 75 base pairs of the probes defining the tails of the regions. Finally, regions were merged if their distance was shorter than 1,000 bp. Bivalent regions were defined as regions sharing at least three probes between H3K4me3 and H3K27me3 regions out of which at least two were enriched.
27 regions were selected for confirmation of enrichment/no enrichment of H3K4me3 and H3K27me3 in tumor tissue detected by array analysis (see Supplemental Data Table S4 for primer sequences). SYBR Green detection chemistry and a standard curve method were used. A standard curve was generated from a dilution series of the input DNA to produce four quantities. The primer efficiency was calculated using the formula E% = (10(-1/slope)-1) × 100. The relative fold enrichment for each H3K4me3 positive region tested, in total 13, were compared to all H3K4me3 negative regions tested, in total 14 (Supplemental Data Figure S1A). The H3K4me3 positive and negative regions were selected for being H3K27me3 negative and positive regions, respectively. The same procedure was followed for validating the H3K27me3 positive regions (Supplemental Data Figure S1B).
Total RNA was extracted from tissue sections of the tumor and normal tissue according to Trizol protocol (Invitrogen) and was followed by a cleanup step using Qiagen RNaesy kit. RNA from the two samples was then hybridized on Affymetrix U133 Plus 2.0 arrays according to Affymetrix recommendations. Raw array data are available from ArrayExpress, accession: E-TABM-811.
Gene expression data were normalized using the MAS5 procedure of the R (R Development Core Team 2008) package simpleaffy, available from Bioconductor (Gentleman et al., 2004). Only probes that were called "Present" in either of the two arrays were considered reliable and kept for further analysis. From the resulting set of genes in the tumor sample, we defined the highly expressed ones as those with expression level above the average level plus one standard deviation. An equal amount of genes were selected among the lowest ranked ones and were considered to be of low expression. We determined genes with increased and decreased expression using a fold change (log2) cutoff of 1.5 and -1.5, respectively.
Regions enriched with H3K4me3, H3K27me3 or both were annotated to UCSC (Karolchik, et al., 2008) known genes (KGs) in a sequential manner. Firstly, all KGs enriched with a histone modification within 2.5 kb downstream of transcription start site were associated with that enrichment. Secondly, KGs were associated with enrichment if located within 2.5 kb upstream of the transcription start site if that enrichment were not located intragenic in any other KG. Using this procedure, KGs became associated with H3K4me3, H3K27me3, bivalency, semi-bivalency or no enrichment. The semi-bivalent KGs were considered those enriched with both H3K4me3 and H3K27me3 but not in a bivalent manner.
KGs were associated with DNA copy number statuses gain, loss or normal if according to overlap with the inferred copy number regions (see above). When we compared the distribution of histone modification and copy number associations in the various expression groups, we considered only non-conflicting cases. For instance, KGs overlapping with more than one copy number region were discarded.
Affymetrix probe sets were mapped to KGs using UCSC annotations. KGs were mapped to Entrez Genes and HGNC gene symbols using UCSC annotations.
The web based Gene Ontology (GO) analysis tool DAVID (Dennis, et al., 2003; Huang da, et al., 2009) was used for functional annotation and to determine statistical over-representation of GO terms among genes (Entrez Genes) with enrichment of histone modifications. Only GO Biological Process terms of level 2 were considered.
The gene networks were generated using Ingenuity Pathways Analysis (Ingenuity® Systems, www.ingenuity.com) on selected genes (Entrez genes). A fold change cutoff of 1.5 was set to identify genes whose expression was differentially regulated. These selected genes, called focus genes, were overlaid onto a global molecular network extracted from the Ingenuity Pathways Knowledge Base. Networks of these focus genes were then generated based on their connectivity. Along with the gene expression, we used information about copy number status and histone modification associations of the genes to examine possible causes of differential expression.
Supplemental Data include six tables and three figures and can be found with this article online at http://e-emm.or.kr/article/article_files/SP-42-7-02.pdf.
We thank Mehdi Motallebipour for technical advice and Peyman Björklund for support. G.W. was supported by the Swedish Research Council and G.Å. by the Swedish Cancer Society. R.A., S.E., A.R.-I. and J.K. were partially supported by the Knut and Alice Wallenberg Foundation and by the Swedish Foundation for Strategic Research. C.W. was supported by the Swedish Research Council, Science and Technology and by the Swedish Cancer Society.