Integrated Profiling Reveals a Global Correlation between Epigenetic and Genetic Alterations in Mesothelioma

Pleural mesothelioma is a rapidly fatal asbestos-associated malignancy, with ∼80% of pleural mesothelioma cases reporting a known history of exposure (1). Asbestos fibers have been shown to be both cytotoxic and clastogenic in vitro (2, 3), and interference with the mitotic machinery by asbestos can lead to abnormal chromosome segregation and hence deletion events (4). Although the number of phenotypically important point mutations in pleural mesotheliomas has been shown to be relatively low (5, 6), reports of extensive gene copy number (CN) alterations in this disease are numerous, and common regions of allele loss include 1p, 3p21, 6q, 9p21, 15q11-15, and 22q (7–14). There have also been several reports of promoter hypermethylation at tumor-suppressor gene loci in mesothelioma (6, 15–21). Gene silencing by DNA methylation at CpG dinucleotides in promoter regions is a well-recognized mechanism of gene inactivation and known to contribute to tumorigenesis (22). We recently reported that increased asbestos exposure is significantly associated with an increased prevalence of gene promoter methylation at cell cycle control genes (23). Other recent work from our group confirmed the relationship between asbestos exposure and DNA methylation in an investigation of more than 750 cancer-related genes, and defined the distinct profiles of coordinated DNA methylation in mesotheliomas relative to nontumorigenic pleural tissues (24).

The spectrum of somatic alterations in mesothelioma and all human cancers includes both genetic and epigenetic events that act in concert to drive tumorigenesis and promote progression to a malignant phenotype. To date, there are only a few examples of studies that aim to perform an integrative analysis of multiple types of somatic alterations. The Cancer Genome Atlas Research Network examined mutations, CN alterations, and DNA methylation in an integrative manner and highlighted alterations of core pathways in glioblastoma (25). Further, a subset of colorectal cancers is known to possess a DNA methylator phenotype, and an association between tumors harboring this phenotype and a lower degree of chromosomal aberrations strongly suggests that there are independent mechanisms of tumor progression in this disease (26). Although integrative genomics is still in its infancy, these types of approaches have tremendous potential to affect cancer biology and translational research as they may provide insight into more generalized mechanisms at work in carcinogenesis and thus targets for therapeutic intervention that could have broader impact.

In pleural mesotheliomas, we investigated both genetic and epigenetic alterations using high-throughput, array-based approaches and established robust analytic strategies with the objective of gaining a more complete understanding of how genetic and epigenetic alterations may interact to drive tumorigenesis.

Mesotheliomas (n = 158) and grossly nontumorigenic parietal pleura (n = 18) were obtained following surgical resection at Brigham and Women's Hospital through the International Mesothelioma Program from a pilot study conducted in 2002 (n = 70) and an incident case series beginning in 2005 (n = 88) as previously described (27). All patients provided informed consent under the approval of the appropriate institutional review boards. Clinical information, including histologic diagnosis, was obtained from pathology reports. The study pathologist confirmed the histologic diagnoses and further assessed the percent tumor from resected specimens (mean >60%).

DNA from fresh frozen tissue and matched whole blood was isolated with QIAamp DNA mini kit (Qiagen) following the manufacturer's protocol. Tumor DNA was modified with sodium bisulfite using the EZ DNA Methylation Kit (Zymo Research). Illumina GoldenGate methylation bead arrays interrogated 1,505 CpG loci associated with 803 cancer-related genes processed at the University of California San Francisco Institute for Human Genetics, Genomics Core Facility, using methods described in ref. (28). The GoldenGate methylation data used in the analysis has been previously described (24). Methylation array data are publicly available through the Gene Expression Omnibus (GSE20989).

From the total study population, 23 tumors from the incident case series were chosen for CN alteration profiling by hybridizing 5 μL containing ≥50 ng/μL of tumor or matched peripheral blood DNA to each of the two GeneChips that comprise the Human Mapping 500K single-nucleotide polymorphism (SNP) array set (Affymetrix; ref. 29), following manufacturer protocols and standard operating procedures at Harvard Partners Microarray Core servicing facility. Probe intensities at each locus were determined in the GCOS software, and genotype calls were generated using the Genotyping Analysis Software (Affymetrix; ref. 30). Probe signals were normalized to their matched reference peripheral blood sample data using the Copy Number Analysis Tool v4.0.1 software (CNAT; ref. 31; Affymetrix) with median scaling and default tuning parameters, and CN states were inferred by hidden Markov model analysis. SNP array data are publicly available through the Gene Expression Omnibus (GSE21057), and the combined set of array data (methylation and SNP data) can be accessed through the Gene Expression Omnibus superseries GSE21058.

Illumina BeadStudio Methylation software was used for methylation dataset assembly. Fluorescent signals for methylated (Cy5) and unmethylated (Cy3) alleles give methylation level β = [max(Cy5, 0)] / (|Cy3| + |Cy5| + 100), with ∼30 replicate bead measurements per locus. Detection P-values determined poorly performing CpG loci (n = 8), which were removed from analysis. X chromosome loci were also removed, leaving 1,413 CpG loci associated with 773 genes. All 23 tumors analyzed for CN alteration passed quality control.

Subsequent analyses were carried out in the R statistical software package (32) and were restricted to autosomal chromosomes. Kaplan-Meier survival strata and log-rank tests were used for univariate survival analyses, and Cox proportional hazards models were used to control for potential confounders. For inference, data were clustered using a recursively partitioned mixture model (RPMM), and the number of classes was determined by recursively splitting the data with two-class models, with Bayesian information criterion used at each potential split to decide whether the split was to be maintained or abandoned as described by Houseman and colleagues (33) and as used in refs. (24, 34).

GoldenGate CpG loci were matched to Affymetrix SNPs in the manner described in ref. (35), typically within 1 kb. The local relationship between methylation and CN alterations at 1,413 loci was tested by calculating the Pearson product-moment correlation coefficients. P values were calculated through permutation test (5,000 permutations), and the qvalue package in R was used to correct for multiple comparisons (Q values). To assess correlation at the gene (rather than locus) level, equivalent statistics were used on CpGs matched to genes by chromosomal position, assigned promoter status if they were upstream of the transcription start site, and averaged by gene. Similarly, CN calls were averaged together for all SNPs associated with a gene.

To investigate molecular alterations individually at each locus, two-sided, two-sample t tests assuming unequal variance were used to compare methylation between the 23 tumors and 18 nondiseased tissues, and a one-sample t test was used with the normal log₂ ratio assumed to be 0, for CN alterations. CN log₂ ratios were used in this analysis due to the possibility that discretization imposed by CNAT might bias associations toward zero. Mean alteration differences were considered statistically significant where P < 0.05. The association between DNA CN and methylation class membership was tested using the mean value of |CNS-2|, where CNS is the CN state of each of 500,446 loci, as defined by the CNAT. A permutation test with 5,000 iterations using the Kruskal-Wallis test statistic was performed.

A representative series of mesotheliomas and their matched peripheral blood–derived DNA as referents (n = 23) from the total study population were subjected to 500K SNP mapping arrays to assess CN alteration (Supplementary Table S1). A summary of the CN alterations identified is provided in Supplementary Table S2. Consistent with previous reports (36–38), we observed prevalent allele loss at 1p36 (35%), 1p21.3 (30%), 4q22 (30%), 4q31-32 (35%), 3p21.3 (44%), 6q25 (39%), 9p21 (39%), and 22q (44%). In addition, we observed prevalent gains at 1q23 (35%), 5p (22%), 7p (22%), and 8q24 (22%). Across tumor samples, the overall prevalence of CN losses and CN gains varied widely, with a mean prevalence of CN loss events of 8% (SD 9%) and a mean prevalence of CN gains of 5% (SD 6%; Supplementary Table S3). Tests for association between extent of CN alteration and age, gender, histology, and asbestos exposure burden were not statistically significant.

Using prior data generated from interrogating DNA methylation in 158 mesotheliomas (24), we applied RPMM to infer methylation profile classes wherein the 23 tumors profiled for CN alterations were included. This approach built classes of samples based on profiles of methylation with data from all measured autosomal loci using a mixture of β distributions to recursively split the tumors into parsimoniously differentiated classes (39–41). Because methylation profile classes from the established model of all 158 mesotheliomas represent a more stable classification, we did not model the subset of 23 tumors with CN data separately. Instead, original class membership data were used in our analyses, and among the seven RPMM methylation classes, the 23 tumors with CN data had a distribution of class membership that included all but methylation class 5. Tumors in methylation class 2 were previously shown to be from patients with significantly lower asbestos burden than other classes. Moreover, relative to these class 2 cases, classes 4 and 7 had significantly poorer prognosis (24).

Figure 1 directly compares CN alterations and DNA methylation state in side-by-side, chromosome-specific plots of these alterations arranged by methylation class membership. Remaining chromosomes are illustrated in Supplementary Fig. S1. Qualitatively, when focusing on single loci or specific genes, there was no obvious or consistent pattern of altered methylation when CN alterations were present. However, we next quantitatively evaluated whether—consistent with Knudson's two-hit paradigm of gene inactivation—coordinate dysregulation of CN state and methylation was present. No loci were found to have a significant correlation (Q < 0.05) between CN and methylation state across the examined tumors (Supplementary Table S4).

By averaging CN and methylation state within gene-specific regions, we expanded this approach to a gene-level query. Only ITGB1 and FYN had significantly correlated (Q < 0.05) CN and methylation among 663 examined genes (Supplementary Table S4). However, Fig. 2 shows that these significant correlations were not due to coordinate allele loss and hypermethylation. Other genes with highly correlated CN and methylation were also plotted. CDKN2A had widespread methylation of one of its target CpGs, although this site is normally methylated in nontumor pleura (Fig. 2). Only two genes, TGFB2 and GDF10, had visually detectable coordinate methylation and CN loss, and only one tumor had both alterations in each case (Fig. 2). Although not among the genes in Supplementary Table S4, RASSF1A is commonly inactivated in mesothelioma; moreover, in our samples, seven tumors (30%) had allele loss, two (22%) were hypermethylated (Fig. 2), and, consistent with results above, only one tumor had coordinate methylation and CN loss.

To assess whether global rather than discrete trends for methylation and/or CN alteration might be observed, we plotted the magnitude and direction of CN and methylation alterations (log₂ ratios versus a ratio of 0 for unaltered CN and tumor versus nontumor methylation at each gene) to genes by their associated P value (Fig. 3A). The volcano plot of CN alterations was skewed to the left, indicating genome-wide trends for gene-level CN losses (15,790 genes). Similarly, the volcano plot of methylation alterations was skewed to the left, indicating a genome-wide trend for loss of methylation relative to nontumor pleural samples (663 genes). To determine whether CN and methylation alterations differ based on methylation profile, we split samples by the first partition of RPMM methylation classes. More specifically, the unsupervised nature of the RPMM model allows stratification of samples by branches of its associated dendrogram; because the model recursively partitions the methylation data, class membership can be retraced to the original parent methylation classes—the “left branch” and the “right branch” classes. Here, the left-branch samples are those belonging to methylation classes 1 or 2 (n = 10), the daughter classes of the initial left branch class, whereas the right-branch samples are those belonging to methylation classes 3 to 7 (n = 13), the daughter classes of the initial right-branch RPMM class. The volcano plots for CN alterations were strikingly different between methylation class branches: tumors in the left-branch classes of the RPMM had far more loci with significantly altered CN compared with tumors in the right branch of RPMM (Fig. 3B). Tumors in the left branch had a wider range of methylation alteration compared with tumors in the right branch, indicating that the trend toward an overall increase or decrease in the degree of methylation was associated with allele CN loss. These results prompted further investigation of the global association between CN alterations and methylation state.

The evidence for a global relationship between CN alterations and methylation alterations could suggest that CN alterations of master epigenetic regulatory genes influence overall tumor methylation. Maintenance methyltransferase DNMT1 had CN loss in seven tumors (30%; Supplementary Fig. S1, chromosome 19), and de novo methyltransferase DNMT3B had no CN alterations (Supplementary Fig. S1, chromosome 20). Plotting the average CpG methylation for tumors with DNMT1 allele loss versus tumors without loss, we observed a significant trend for increased methylation among tumors with no allele loss at DNMT1 compared with tumors with allele loss (Fig. 4; P = 0.05). Further, DNMT1 allele loss was associated with significantly reduced patient survival in a Cox proportional hazards model controlling for age, gender, and tumor histology (HR, 5.07; 95% C.I. 1.23–20.9; Supplementary Table S5). In similar investigations of alterations of critical DNA double-strand break repair genes (e.g., ATM, XRCC4), prevalent CN alterations or methylation silencing events were not observed.

In an effort to better visualize and test the global trends of association between methylation and CN alterations, we plotted the CN alteration profiles of tumors based on RPMM methylation class membership (Fig. 5). This plot illustrates the significant difference in the extent of CN alterations among methylation classes (permutation test P < 0.02). Tumors in methylation classes 2, 6, and 7 exhibit a greater extent of CN alteration than tumors in methylation classes 1, 3 or 4 (Fig. 5). The mean and SEM values of class-specific percents of SNPs with CN alterations are also illustrated in Fig. 5.

The phenotype of any individual tumor arises as the result of a constellation of somatic alterations, and only recently has technology afforded the opportunity to profile and examine multiple types of alterations on a genome-wide scale. In an effort to further our understanding of the relationship between two forms of somatic alterations known to be important in the genesis of mesothelioma, we profiled CN alterations in a subset of tumors for which we had available DNA methylation profiling data. Investigating single loci and specific genes, we found no significant correlations between CN and methylation, indicating that two-hit gene inactivation is not commonly achieved by coordinate hypermethylation and allele loss in mesothelioma. However, when looking genome-wide, we found that the extent of CN alterations was significantly related to the DNA methylation profiles of these tumors, suggesting a strong link between genetic and epigenetic dysregulation in mesothelioma. The association between asbestos and mesothelioma carcinogenesis is well established: asbestos fibers are clastogenic, and patient asbestos burden has been associated with methylation alterations in mesothelioma (3, 24). Combined with the known, low mutagenic action of asbestos (e.g., mesotheliomas almost never harbor TP53 mutations and generally have few mutations; ref. 5), these tumors are an excellent candidate for integrative genomic profiling of CN and epigenetic alterations.

It is possible that the observed global association between CN alteration and methylation could be influenced by a bias in methylation measurement at regions with CN alteration. If methylation measurements were biased by CN state, the bias could contribute to the classification of DNA methylation profiles modeled by RPMM. In fact, previous work has indicated that certain microarray platforms will be affected by just such a bias (42). However, we recently performed an extensive investigation including multiple tumor types and showed that the Illumina GoldenGate bead-array-based platform for methylation measurement is not subject to significant bias from CN alterations (35).

Despite generally extensive CN and methylation alterations, we observed very few instances of locally coordinate allele loss and methylation in these tumors. Although two genes had significant correlations between CN loss and methylation, neither case had high enough methylation levels to be consistent with gene silencing. On the other hand, among other genes with highly ranked correlation between methylation and CN alteration, there were two single instances of coordinate CN loss and hypermethylation, one tumor at TGFB2 and one tumor at GDF10. The specific functions of TGFB2 are dependent on cell type and the presence of both TGF-β type II and type I receptors (43). In mesothelioma, it is thought that TGFB2 may play complex roles in growth and regulation of immune responses through cytokine production (44). Interestingly, the other target of coordinate methylation and CN loss was GDF10, which encodes bone morphogenic protein 3b; the bone morphogenic proteins are part of the transforming growth factor-β superfamily, and a low prevalence of GDF10 methylation has previously been shown in mesotheliomas (45).

Despite a lack of locally coordinated CN and methylation alterations, we observed strong evidence for an association between genome-wide extent of CN alteration and methylation alterations. In fact, a global propensity for both CN and methylation losses was observed. The trend for methylation losses in these tumors was previously reported in the total study population (n = 158): relative to nontumor pleura, among 1,413 autosomal CpG loci, nearly 1,000 CpGs had significantly altered methylation between nontumor pleura and mesotheliomas (24), and more than 75% were losses of methylation. Hence, one potential explanation for the lack of coordinate CN loss and hypermethylation is a low prevalence of hypermethylation events in these tumors. In the context of a propensity for methylation loss, further investigation of the global coordination between CN alterations and methylation led us to the methyltransferases themselves. Not only was there prevalent CN loss of the maintenance methyltransferase, but also, tumors without DNMT1 loss had higher average methylation across all CpGs than those with loss. It is reasonable to suggest that loss of DNMT1 may also result in hypomethylation of repeat elements and thereby further contribute to genomic instability and increased CN alterations, although additional investigations are necessary. Furthermore, we found that loss of DNMT1 was associated with a significantly increased risk of death. We have previously reported associations between methylation profile class membership and survival in mesothelioma (46). It is possible that DNMT1 loss is also associated with methylation class membership (potentially those classes with significantly poorer survival), although additional investigation is necessary. Nonetheless, when studying human tumors, it may not be possible to estimate the timing of critical events that contribute to dysregulation of genetic and/or epigenetic maintenance.

We have shown that fulfillment of the two-hit paradigm is generally not achieved by coordinate methylation and CN alteration mechanisms in mesothelioma, but that there is a strong association between CN alteration and methylation state in these tumors. Initial exploration of the cause(s) behind overall CN and methylation associations indicate promise for the investigation of master regulatory genes (such as the DNA methyltransferases) in additional studies with integrative genomic approaches. Other forms of gene inactivation mechanisms not measured here (e.g., miRNA/mRNA alterations, gene mutation) may still act coordinately with methylation or CN alteration to fulfill the two-hit paradigm at a gene-specific level, may also be globally associated with one another, and should be investigated in future studies. Certainly, further studies including higher-resolution methylation arrays, additional forms of somatic alteration, and/or different tumor types with distinct etiologies are all warranted. In fact, in solid tumors of the head and neck, our group has now reported correlation between overall CN alteration extent and methylation profiles in the absence of gene-specific coordinate inactivation events, strong evidence that our findings may extend to many other forms of cancer (47). Collectively, these data further show the feasibility and promising utility of integrative genomic approaches for the advancement of cancer biology and translational research.

No potential conflicts of interest were disclosed.

Grant Support: R01CA126939, R01CA100679, International Mesothelioma Program at Brigham and Women's Hospital (research grant), and Mesothelioma Applied Research Foundation (research grant).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.