We performed cytosine methylation sequencing on genetically diverse patients with acute myeloid leukemia (AML) and found leukemic DNA methylation patterning is primarily driven by nonpromoter regulatory elements and CpG shores. Enhancers displayed stronger differential methylation than promoters, consisting predominantly of hypomethylation. AMLs with dominant hypermethylation featured greater epigenetic disruption of promoters, whereas those with dominant hypomethylation displayed greater disruption of distal and intronic regions. Mutations in IDH and DNMT3A had opposing and mutually exclusive effects on the epigenome. Notably, co-occurrence of both mutations resulted in epigenetic antagonism, with most CpGs affected by either mutation alone no longer affected in double-mutant AMLs. Importantly, this epigenetic antagonism precedes malignant transformation and can be observed in preleukemic LSK cells from Idh2R140Q or Dnmt3aR882H single-mutant and Idh2R140Q/Dnmt3aR882H double-mutant mice. Notably, IDH/DNMT3A double-mutant AMLs manifested upregulation of a RAS signaling signature and displayed unique sensitivity to MEK inhibition ex vivo as compared with AMLs with either single mutation.
Significance: AML is biologically heterogeneous with subtypes characterized by specific genetic and epigenetic abnormalities. Comprehensive DNA methylation profiling revealed that differential methylation of nonpromoter regulatory elements is a driver of epigenetic identity, that gene mutations can be context-dependent, and that co-occurrence of mutations in epigenetic modifiers can result in epigenetic antagonism. Cancer Discov; 7(8); 868–83. ©2017 AACR.
This article is highlighted in the In This Issue feature, p. 783
Acute myeloid leukemia (AML) is a genetically and biologically heterogeneous disease that features a relative paucity of somatic mutations compared with other tumor types (1). In contrast, patients with AML exhibit widespread disruption of cytosine methylation patterning (1, 2). These specific DNA methylation profiles have been exploited to classify AML into 16 different epigenetically defined subtypes that were partially linked to genetic lesions (3, 4) and enabled discovery of novel leukemia driving mechanisms (2–4). All of these studies focused primarily on cytosine methylation patterning at promoter areas and CpG islands. Yet, more recent studies in cancer epigenetics have highlighted potential roles for other types of genomic elements, such as enhancers and CpG shores, in contributing to perturbed epigenetic patterning in cancer (5–7). Enhancers are distal regulatory elements that often become active or inactive in a tissue-specific manner to mediate particular developmental and functional processes (8). Enhancer activity status depends on the presence of specific histone modifications and other epigenetic marks, including cytosine modifications. Enhancers are increasingly shown to play critical roles in myelopoiesis and leukemogenesis (9–11). Recent studies implicate epigenetic enhancer dysfunction through aberrant cytosine methylation as a contributing factor for development of AML in mouse models, for example, in the case of EVI1 leukemias (9, 10) or mutant TET2 in combination with FLT3 internal tandem duplications (ITD; ref. 12). Hence, it is of great interest to explore whether these elements are indeed a prominent and perhaps even defining feature of epigenetic patterning in human primary AML.
One potential driving force for epigenetic reprogramming in AML is the occurrence of somatic mutations affecting proteins that regulate cytosine methylation. DNMT3A, which encodes a de novo DNA methyltransferase, is affected by somatic mutations in approximately 25% of AML cases (1). Evidence to date would appear to indicate that certain DNMT3A mutants could function in a dominant-negative manner to disrupt active DNMT3A tetramers (13–15). DNMT3A mutations occur during the preleukemic phase of AML pathogenesis (16, 17), attesting to a potential critical role for aberrant epigenetic reprogramming in malignant transformation of hematopoietic cells. The IDH1 and IDH2 loci encode metabolic enzymes that are also mutated in 15%–20% of AML cases (1, 18). Mutant IDH1 and IDH2 proteins generate the oncometabolite 2-hydroxyglutarate. This inhibits dioxygenases, such as the TET family enzymes, which are involved in DNA demethylation (19), resulting in a DNA hypermethylation phenotype in AMLs carrying these mutations (4). Murine models of DNMT3A and IDH mutations have revealed that aberrant methylation driven by these lesions is not restricted to promoter regions but instead targets a significant proportion of intergenic and intronic regions (20, 21).
To better understand the nature of aberrant epigenetic programming in AML, we utilized enhanced reduced representation bisulfite sequencing (ERRBS) on a panel of primary AML cases with detailed cytogenetic, molecular, and epigenetic information (2). The expanded coverage and base-pair resolution of ERRBS allowed us to evaluate cytosine methylation changes at higher resolution and with greater breadth than previously possible. We evaluate the contribution of individual gene mutations to these epigenetic patterns and specifically focus on the effect of co-occurring epigenetic mutations.
Comprehensive DNA Methylation Profiling Extending beyond Promoter Regions Robustly Captures AML Biology
AML can be classified into biologically distinct subtypes using gene promoter methylation patterns (2, 3). However, increasing evidence points toward a key regulatory role for DNA methylation beyond promoter regions (6, 7). We hypothesized that higher-resolution methylation profiling using next-generation bisulfite sequencing would better define the basis for aberrant epigenetic programming in AML. We therefore performed ERRBS (22) on 119 adult patients with AML selected from a cohort that had been previously studied using the promoter-centric HELP microarray assay and classified into 16 epigenetically defined clusters (2). Samples were randomly selected from within each of the original 16 epigenetically defined clusters and curated to include examples of the key somatic mutations currently understood to be associated with AML pathogenesis (summarized in Table 1).
|.||N (%) .||.|
|Gender Male||66 (55.46%)|
|M4||21 (17.65%)||1 M4E|
|NA||5 (4.2%)||Includes 4 “AML non-M3,” 1 RAEB MDS|
|tri 8||2 (1.68%)|
|Poor||13 (10.92%)||Includes 1 favorable + unfavorable|
|NA||6 (5.04%)||Unavailable karyotype|
|Double mutation||12 (10.08%)|
|Single mutation||1 (0.84%)|
|EVI1 overexpression||8 (6.72%)|
|.||N (%) .||.|
|Gender Male||66 (55.46%)|
|M4||21 (17.65%)||1 M4E|
|NA||5 (4.2%)||Includes 4 “AML non-M3,” 1 RAEB MDS|
|tri 8||2 (1.68%)|
|Poor||13 (10.92%)||Includes 1 favorable + unfavorable|
|NA||6 (5.04%)||Unavailable karyotype|
|Double mutation||12 (10.08%)|
|Single mutation||1 (0.84%)|
|EVI1 overexpression||8 (6.72%)|
Abbreviations: FAB, French-American-British acute myeloid leukemia classification; M4E, M4 with eosinophilia; MDS, myelodysplastic syndrome; RAEB, refractory anemia with excess blasts; tri 8, Trisomy 8.
Unsupervised analysis by hierarchical clustering using Euclidean distance and based on the most variable CpGs across all patients (n = 142,643) yielded 14 clusters at a distance threshold that minimizes the number of clusters without merging samples with distinct cytogenetic or mutation profiles (Fig. 1A). Clusters characterized by inv16, t(8;21), t(15;17), t(v;11q23), and CEBPA silencing remained identical to those specified by HELP. In contrast, cases that originally fell into a single cluster harboring patients with CEBPA double mutation (CEBPA-dm) were instead split into two distinct groups by ERRBS, one that was almost entirely cytogenetically normal (cluster 10), and another that included cases with both normal karyotype and cytogenetic abnormalities (cluster 9). Similarly, IDH1/2-mutant cases, which were previously mainly grouped into two clusters including both IDH1 and IDH2 mutant cases (4), were now grouped by ERRBS profiling into a predominantly IDH1-mutant cluster (cluster 1) in which all mutant cases also harbored DNMT3A mutations, and a second cluster exclusively carrying IDH mutations without co-occurring DNMT3A mutations, almost all of which were IDH2 mutations (Supplementary Table S1). The remaining clusters were enriched for a combination of DNMT3A, NPM1, and FLT3 abnormalities. These data indicate that comprehensive methylome sequencing more precisely links genetic lesions to perturbations in cytosine methylation patterning than previous, promoter-centric approaches (1–3).
To further and more comprehensively assess the link between cytosine methylation signatures and somatic mutations, we performed additional targeted resequencing of 49 genes implicated in myeloid malignancies (Supplementary Table S2). Using a series of multinomial linear regression models (see details in the Methods section), we then quantified the ability of individual mutations, genetic rearrangements, or known epigenetic abnormalities to specify epigenetic clustering. In addition to inv16, t(8;21), t(15;17), and t(v;11q23), only CEBPA-dm, NPM1, DNMT3A, and IDH2 mutations, or CEBPA gene silencing provided additional significant predictive ability (Supplementary Fig. S1A). No additional mutation or pair of mutations (including those affecting FLT3) contributed additional predictive power (Supplementary Fig. S1B). Notably, the presence of an IDH1 mutation, which was cluster-defining only in combination with DNMT3A mutation, did not emerge as an independent predictor of clustering. Therefore, even when performing more comprehensive methylome analysis, only a subset of somatic mutations can be linked to specific cytosine methylation patterning effects. TET2 mutation status also did not emerge as a cluster-defining lesion, instead being distributed among several clusters.
Differential Methylation at Non–Promoter Gene Regulatory Regions Is the Primary Driver of AML Epigenetic Subtypes
The genomic coverage of ERRBS extends beyond promoters to cover intergenic and intragenic CpGs, as well as beyond CpG islands to also cover CpG shores and distal intergenic regions (22). In order to determine which genomic elements have the greatest influence on determining cluster-defining epigenetic patterning, we first annotated CpGs covered by ERRBS into the following categories: (i) “promoter”—CpG +/− 2 kb from gene transcription start sites (TSS); (ii) “gene neighborhood”—CpGs between 2 kb and 50 kb away from the TSS or transcription end site (TES); (iii) “downstream”—CpG +/− 2 kb from transcription end sites; (iv) “intergenic”—CpGs greater than 50 kb from the nearest gene; (v) “exons”—CpGs within any annotated exon; (vi) “introns”—CpGs between exons (Fig. 1B). A total of 950,953 common CpGs with coverage ranging between 10× and 400× were uniformly captured in all 119 samples. Among these common CpGs, 44.7% were located within promoters, 16.1% within gene neighborhoods, 7.9% within intergenic regions, 1.5% within the downstream window, 9.2% within non-first exons, and 20.6% within non-first introns (Supplementary Fig. S2A). From a CpG island–centric point of view, 42.6% of CpGs were at CpG islands, whereas 9.7% were at CpG shores, and 47.6% were located elsewhere (Supplementary Fig. S2B). Taken together, these CpGs captured 82.6% of all RefSeq promoters, 91.6% of gene neighborhoods, 88.9% of intergenic regions, 16.7% of non-first exons, 33.9% of non-first introns, as well as 93.2% of all CpG islands, and 72.6% of CpG island shores (Supplementary Fig. S2C). Overall, these CpGs represent an extensive sampling of the epigenome.
To determine which of these genomic regions contributes the most to epigenetic patterning in AML, we used the adjusted Rand index to compare clustering using only CpGs located at the various genomic regions to that of the overall clustering. We found that CpGs located in gene neighborhoods were best able to accurately recapitulate epigenetic clusters despite representing only 15% of the CpGs used for the unsupervised clustering analysis. When compared with the same number of CpGs selected from the genome at random, we found this prediction to be significantly better than chance (P < 0.01, z-test). Other genomic subsets also predicted clustering better than chance, but with poorer cluster specification as indicated by their lower adjusted Rand indices (Fig. 1C). Although promoters classified patients poorly, the subset of CpGs in promoters within CpG shores performed far better than CpGs within CpG islands, underlining that CpG shores and not islands may be the principal sites for contributing to epigenetic identity in AML (Fig. 1D). Next, we sought to determine whether this general observation was equally true across all clusters or whether certain AML subtypes were best captured by specific compartments. A cluster-by-cluster performance analysis for each genomic compartment revealed that gene neighborhoods were indeed most accurate in predicting cluster assignments in most cases (10 cases out of 14). However, they were only partially able to capture the IDH1/2 clusters, one of the NPM1 clusters, and the EVI1 cluster. This may reflect the fact that differential methylation in IDH mutant cases tends to occur relatively nearer to TSSs (22). Overall, the remaining genomic compartments showed weaker clustering accuracy compared with gene neighborhoods (Fig. 1E). The only AML subtypes captured uniformly by all genomic regions were those with t(15;17), inv(16), t(v;11q23), and those with CEBPA silencing. These findings suggest that oncogenic drivers primarily disrupt cytosine methylation at gene neighborhood regions, which are enriched in regulatory elements such as enhancers and transcription factor binding sites.
Active Enhancers Exhibit Focal Enrichment of Differential Methylation in AML
To determine the contribution of enhancer CpGs to epigenetic cluster identity, we determined the location of hematopoietic stem and progenitor cell (HSPC) enhancers using chromatin immunoprecipitation sequencing (ChIP-seq) datasets generated in human-mobilized CD34+ cells (23) based on H3K4me1 and H3K27ac peaks in the absence of a H3K4me3 peak or an annotated promoter. Active enhancers were distinguished from poised enhancers by the presence of an H3K27ac peak. Enhancer-associated CpGs were predominantly within introns and gene neighborhoods, with 19.4% of active enhancers captured located at gene neighborhoods and 48% within introns (Supplementary Fig. S2D). Gene neighborhoods were notably enriched for both active and poised enhancer CpGs (P < 0.01 for each, Fisher exact test).
To determine if enhancers are generally disrupted by leukemic cytosine methylation patterning, we specifically interrogated methylation change within active HSPC enhancers. We found that these loci exhibited focal methylation changes in AMLs when compared with normal CD34+ cells, and that the average enhancer methylation change at all covered CpGs was greater than surrounding sequences (Fig. 2A; Supplementary Fig. S3A and S3B). This pattern was in sharp contrast to promoters, which exhibited little methylation change relative to surrounding sequences (Fig. 2B; Supplementary Fig. S3C and S3D). These changes were not equal in nature in all clusters and included both some with increases in average methylation and some with decreases in average methylation (Supplementary Fig. S4A). In absolute terms, the changes at promoters were both less marked and more uniform in nature across the different clusters. In addition, the change consisted almost exclusively of gains in methylation in areas surrounding the TSS with little to no change at the TSS itself (Supplementary Fig. S4B).
At differentially methylated CpGs, the directionality of methylation change within these two compartments also differed. Among CpGs exhibiting significant differential methylation in AML (≥25% methylation difference, β-binomial q ≤ 0.01) relative to normal HSPC controls, the majority of clusters showed a predominance of hypomethylation at active enhancers with the exception of CEBPA-silenced and IDH2-mutant clusters. In contrast, promoters were generally hypermethylated, with the exception of DNMT3A and MLL fusion leukemias (Fig. 2C and C). Therefore, aberrant cytosine methylation of HSPC enhancers is a defining feature of AML epigenetic patterns, with directionality of methylation dependent on driver oncogenes, whereas gene promoters are involved to a lesser extent.
Epigenetically Defined AML Patient Clusters Exhibit Distinct Regional Distribution Patterns and Directionality
To characterize the extent to which cytosine methylation is perturbed in each epigenetically defined AML subtype, we performed supervised analyses comparing the DNA methylation profiles of the AML cases in each cluster to ERRBS performed in healthy human bone marrow CD34+ HSPC controls (n = 22). In all but one cluster, differentially methylated CpGs (DMC; ≥ 25% methylation change, β-binomial q ≤ 0.01) represented less than 10% of CpGs covered by ERRBS (Fig. 3A; Supplementary Fig. S4C). The CEBPA-silenced cluster was the only outlier, with nearly 19% of covered CpGs hypermethylated relative to CD34+ controls. Notably, the cluster defined by coexistent DNMT3A and IDH1 mutations contained the fewest total DMCs (Fig. 3A; Supplementary Fig. S4D).
Next, we evaluated the distribution of differential methylation in each cluster of patients with AML (Supplementary Table S3). Gene neighborhoods were the only genomic region enriched for DMCs in all clusters (P < 0.01 for all, Fisher exact test). However, DMCs within clusters characterized by a predominantly hypermethylated signature, such as the CEBPA-silenced cluster, the CEBPA-dm clusters, and del (5/7q), were preferentially enriched at promoter regions (P < 0.01 for each, Fisher exact test). In contrast, DMCs within clusters characterized by a predominantly hypomethylated signature, such as the DNMT3A mutant and t(v;11q23) clusters, were preferentially enriched at distal intergenic and intronic regions (P < 0.01 for each, Fisher exact test; Fig. 3B and C). Interestingly, less than half of the clusters showed enrichment of differential methylation at promoter regions (Fig. 3C; Supplementary Fig. S3D), but all showed significant enrichment at HSPC enhancers, regardless of genomic compartment (Fig. 3D). Altogether, these data underline the markedly diverse nature of epigenetic programs occurring in AML and again point to enhancers and other nonpromoter regulatory regions as the critical disease-defining sites for epigenetic dysregulation.
AMLs with Co-occurring DNMT3A and IDH Mutations Exhibit Epigenetic Antagonism
Clusters 4 and 11 capture the hypermethylated and hypomethylated profiles for IDH2 and DNMT3A mutations, respectively, consistent with the known molecular mechanisms of these two mutant proteins (4, 13, 15, 21, 24). However, patients with AML with dual IDH and DNMT3A mutations exhibited a markedly distinctive cytosine methylation pattern from patients in either cluster 4 or 11. Among these patients, 8/11 (73%) exhibited dual IDH1/DNMT3A mutations, so this cohort will be referred to as IDH1/DNMT3A. Indeed, IDH1/DNMT3A patients manifested the smallest degree of differential methylation of any cluster, with only 31,882 DMCs relative to CD34+ controls. Of these, 47.3% of DMCs (15,072 CpGs) were hypomethylated and 52.7% (16,810 CpGs) were hypermethylated (Fig. 3A). In striking contrast, DNMT3A mutations alone (cluster 11) resulted in 105,939 DMCs, 84.8% of which were hypomethylated, whereas isolated IDH2 mutations (cluster 4) led to a predominantly hypermethylated profile with 148,961 DMCs, 93.9% of which were hypermethylated.
To better understand the nature of these differences, we compared and contrasted the regional distribution of differential methylation in DNMT3A, IDH2, or IDH1/DNMT3A patients (Supplementary Table S4). Whereas IDH2 mutants exhibited an average DMC methylation increase >40% above normal HSPC controls at neighborhoods, promoters, and gene bodies, DNMT3A-mutant cases manifested ∼20% loss of methylation at gene neighborhoods and bodies, and little effect on promoters. In marked contrast, IDH1/DNMT3A double mutants completely lost the characteristic differential methylation of either IDH2 or DNTM3A at gene neighborhoods. IDH1/DNMT3A also manifested attenuated IDH effects at promoters and exons and attenuation of the DNMT3A effect at introns and intergenic regions (Fig. 4A).
We next examined active enhancers in more detail because, as shown above, these regions represent a dominant site of epigenetic variability in AML. This analysis revealed that IDH2-mutant cases manifested overall enhancer hypermethylation when all covered CpGs were considered whereas DNMT3A-mutant cases exhibited predominant enhancer hypomethylation. However, there was very little deregulation of enhancer cytosine methylation in IDH1/DNMT3A cases (Fig. 4B). At promoter regions, differential methylation at all covered CpGs exhibited focal restriction centered at the transcription start site in all three subtypes. Regions flanking the TSS once again manifested overall hypermethylation at all covered CpGs in the IDH2 AMLs, whereas there was a mild degree of hypomethylation in the DNMT3A mutants. In contrast, the IDH1/DNMT3A had very little effect on promoters with almost complete lack of the IDH2 promoter hypermethylation effect (Fig. 4C).
This analysis demonstrates that the overall gain in regional cytosine methylation seen in IDH2 AMLs and the general loss of cytosine methylation observed in DNMT3A AMLs were profoundly attenuated in the IDH1/DNMT3A double mutants. Specifically, these two mutations have opposing effects on gene neighborhoods and enhancers, yet both of these effects are mostly lost in IDH1/DNMT3A double-mutant AML. Enhancers serve as platforms for the action of transcription factors. Using a bioinformatic approach, we observed that differentially methylated HSPC enhancers in IDH2 AMLs were enriched for SMAD3, RXRA, and USF2, all of which contain strong “CG” motifs (Supplementary Fig. S5A), whereas DNMT3A-mutant cases were enriched for SMAD3 and REL NF-κB motifs (Supplementary Fig. S5B). Even though IDH1/DNMT3A double-mutant cases manifest much less enhancer perturbation, there was still some enrichment for NR2C2 motifs, which were not present in IDH- or DNMT3A-mutant cases (Supplementary Fig. S5C). All of these transcription factors are expressed in these leukemias, suggesting that they could be in some way involved or affected by aberrant cytosine methylation (Supplementary Fig. S5D). As a complementary approach, we examined available ChIP-seq data from leukemia and hematopoietic cells (Supplementary Fig. S5E). All three sets of patients featured enrichment for TAL1, GATA1, RUNX1, and others at differentially methylated enhancers. However, GATA2, CEBPA, CEBPB, PRAME, NFYA, and TCF7L2, although uniquely enriched in both single-mutant AMLs, were missing from the double-mutant leukemias, which in turn did not manifest enrichment for any unique set of ChIP-seq transcription factor binding sites.
Mouse hematopoietic stem cells have been observed to contain large domains of unmethylated cytosines containing binding sites of transcription factors involved in hematopoiesis (termed “canyons”; ref. 25). To determine whether these regulators of hematopoiesis are perturbed in AML, we identified 207 of 1,104 of these regions with analogous sites in the human genome covered by our ERRBS data. IDH2-mutant AMLs exhibited relative hypomethylation at these canyons when compared to neighboring regions, although these canyons were still hypermethylated as compared with CD34+ controls. Notably, both DNMT3A only and IDH1/DNMT3A cases also showed mild hypermethylation at these canyons relative to normal HSPC controls, though this was much more attenuated than that seen in IDH2 cases. (Fig. 4D). Hence, although IDH and DNMT3A mutations also exert distinct effects on cytosine methylation at canyons, these seem to trend toward hypermethylation in AMLs relative to normal controls.
Given the opposing effects of these two mutations on methylation profiles, we sought to determine whether these mutations targeted the same CpG sites with opposing results, or whether each targets a distinct set of CpGs. We found that when compared with control HSPCs, AMLs with mutant IDH2 and DNMT3A for the most part target distinct CpG sites (Fig. 4E), and the CpGs that become hypermethylated in IDH2 AMLs are different from those that become hypomethylated in DNMT3A-mutant AMLs. Strikingly, even though these CpG sites do not overlap, a majority of the uniquely hypermethylated cytosines in IDH2 and uniquely hypomethylated cytosines in DNMT3A did not manifest differential methylation compared with control HSPCs in IDH1/DNMT3A double-mutant AMLs (Fig. 4E). These findings indicate that even though IDH2 and DNMT3A individually target a different set of CpG sites, the pattern of methylation change in IDH1/DNMT3A mutants largely reflects antagonistic actions of IDH1 and DNMT3A mutations at most CpGs, while at the same time maintaining a subset of the respective hypermethylation and hypomethylation effects of each allele. The majority of DNMT3A mutations in both DNMT3A alone and the IDH1/DNMT3A cases were at the R882 position (Supplementary Table S1); thus, the expected impact of mutant DNMT3A on the epigenome in these two clusters is similar in both groups.
IDH2 and DNMT3A Directly Induce Epigenetic Antagonism in Preleukemic Hematopoietic Stem Cells
We next sought to determine whether DNMT3A and IDH2 are directly responsible for the epigenetic antagonism observed in human patients, or whether these changes are simply a late consequence of malignant transformation. For this purpose, we collected LSK cells from conditional Mx1-Cre;Dnmt3aR878H knock-in mice, which corresponds to the human DNMT3A R882 mutation (26), 6 months after inducing recombination. Notably, DNA methylation analysis by ERRBS revealed a similar hypomethylated distribution as seen in human DNMT3A-mutant blasts, even though at this time point Dnmt3aR878H knock-in mice do not manifest any hematologic phenotype. Similarly, DNA methylation profiling of LSK cells from Idh2R140Q knock-in mice (27) collected at 6 months after recombination presented a similar aberrantly hypermethylated profile as seen in human AMLs carrying IDH1/2 mutations. Again, these mice do not manifest hematologic phenotypes at this time point. Finally, in order to determine whether the co-occurrence of both mutations results in epigenetic antagonism, we generated a double knock-in with Dnmt3aR878H/Idh2R140Q and performed ERRBS profiling on LSK cells at 6 months after recombination. Importantly, at this time point, these Dnmt3aR878H/Idh2R140Q mice also did not manifest a hematologic phenotype. The observed cytosine methylation profile in this genetically engineered model accurately reflected the epigenetic antagonism in the human double-mutant AMLs, with Dnmt3aR878H/Idh2R140Q LSKs exhibiting far fewer differentially methylated CpGs than either Dnmt3aR878H or Idh2R140Q (Fig. 5A; Supplementary Table S5). The changes observed both at specific genomic regions as well as at individual CpGs for each allele paralleled the behavior seen in the human AML counterparts (Fig. 5B as compared with Fig. 3B). Active enhancers exhibited focal enrichment of methylation variability with Idh2R140Q LSKs exhibiting gain of methylation and Dnmt3aR878H LSKs exhibiting focal loss of methylation (Fig. 5C). At promoters, Idh2R140Q and Dnmt3aR878H LSKs again exhibited hypermethylation and hypomethylation respectively with a restriction in methylation variability at the TSS (Fig. 5D). Within canyons, there was a relative restriction in methylation variability, with immediate canyon flanking sequences exhibiting hypermethylation and hypomethylation in Idh2R140Q and Dnmt3aR878H respectively (Fig. 5E). Although the absolute number of perturbed CpGs is minimal in the Idh2R140Q/Dnmt3aR878H double knock-in mice, there was a tendency for these few cytosines to be less methylated (Fig. 5C–E). Similar to the pattern observed in their AML counterparts, differential methylation in the Idh2R140Q and Dnmt3aR878H mice occurred in largely different genomic compartments, and the majority of CpGs in the double-mutant mice remained unchanged relative to normal controls (Fig. 5F as compared with Fig. 4E). These findings confirm that the epigenetic antagonism observed in human IDH/DNMT3A double-mutant AMLs is directly linked to these alleles. The data also underscore that these epigenetic changes occur early on in preleukemic cells as a consequence of the expression of the mutant alleles and preceding the phenotypic changes associated with complete transformation of the cells.
Transcriptome Analysis Reveals Unique Expression Profile of IDH/DNMT3A Double- Mutant AMLs, Suggesting Increased Sensitivity to MEK Inhibition
Although antagonistic in their effects on cytosine methylation, DNMT3A and IDH mutations must still cooperate to drive the leukemic phenotype. To determine how this might occur, we next evaluated the gene expression profiles of these patients with AML. When compared with normal HSPC controls, IDH1/DNMT3A-mutant AMLs contained 2,524 differentially expressed genes (≥1.5-fold change, adjusted P ≤ 0.05). Of these, only 535 genes were uniquely differentially expressed in the IDH1/DNMT3A cluster, with the remaining genes being shared with the DNMT3A and IDH cluster patients (Fig. 6A). Gene set enrichment analysis (GSEA) using the Hallmark MSigDB gene sets was performed to identify sets of genes unique to IDH1/DNMT3A cases. The gene sets uniquely enriched are suggestive of competing pro- and antiproliferative stimuli. KRAS signaling (NES = 1.36, FDR q = 0.048), which favors cell growth, appears to be opposed by the downregulation of certain MYC targets (V2: NES = −2.03, FDR q < 0.01) as well as upregulation of apoptosis-associated genes including multiple caspases (NES = 1.43, FDR q −0.034). Normal hematopoietic processes such as heme metabolism are ultimately downregulated (NES = −1.8, FDR q < 0.01; Fig. 6B). Given the unique upregulation of the RAS signaling gene set in the double mutant AMLs, we hypothesized that these leukemias might be uniquely sensitive to inhibition of this signaling cascade. We tested this using a liquid culture assay in which primary AMLs carrying either single IDH mutations, single DNMT3A mutations, double IDH/DNMT3A mutations, or AMLs with neither of these two mutations, were cultured in the presence of cytokines for 7 days and the MEK inhibitor MEK162 was added to the culture every 48 hours. Cells with co-occurrence of IDH/DNMT3A mutations showed significantly greater sensitivity to MEK inhibition than any of the other genotypes (P < 0.01; Fig. 6C). This observation was further validated in an orthogonal assay on methylcelullose in which a number of total colonies formed by either primary AMLs or normal mononuclear cells from cord blood were scored after cells were cultured for 14 days in the presence of the inhibitor (Fig. 6D). Although IDH single mutants did not form colonies in this assay, once again IDH1/DNMT3A cases were more sensitive to RAS signaling inhibition via MEK inhibition than either DNMT3A, IDH2, or other types of AML (Fig. 6C and D). Taken together, these data suggest that IDH and DNMT3A mutations each result in mixed pro- and antiproliferative gene expression signatures, which is further perturbed in IDH1/DNMT3A cases in a manner that is at least partly mediated by RAS signaling and not dependent on aberrant cytosine methylation.
Epigenetic dysregulation has emerged as one of the hallmarks of AML. Broad epigenetic changes occur during hematopoietic differentiation, and disruption of this process accompanied by proliferation of immature myeloid forms underlies the development of AML. The diversity of epigenetic signatures occurring in AML is indicative of multiple underlying mechanisms leading to the development and maintenance of the leukemic state. Herein, through the use of base pair resolution DNA methylation sequencing that represents the genome far more thoroughly than previous microarray studies, we reveal the basic nature of the specific and distinct epigenetic programs that occur in AML. Although microarray profiling focused on promoters was sufficient to identify biological subtypes with distinct responses to therapy, significant heterogeneity remained, and several of the resulting categories did not have characteristic mutation profiles (2). We show here that when methylation profiling is extended into gene neighborhoods and intergenic regions, we are better able to classify AMLs and clarify their linkage to mutation profiles.
A major finding of this study is that methylation changes in the gene neighborhood region are better able to predict epigenetic subtype than changes within genes, promoters, or even CpG islands or CpG island shores. In addition, we found that differential methylation appears to particularly affect enhancers that are active in CD34+ HSPCs. In contrast, promoters show a general restriction in differential methylation. Thus, in spite of a relative paucity of somatic mutations, AMLs manifest profound epigenetic disruptions of enhancers that likely reflect differences in the mechanisms that drive malignant transformation. Given that many of the genetic lesions in AML directly or indirectly perturb hematopoietic master regulatory transcription factors, and that these factors largely function through enhancer binding, it is appealing to speculate that the observed distinctive enhancer profiles that distinguish AMLs from each other reflect perturbation of different combinations of transcription factors. Hence, the epigenetic pattern of these AMLs may serve as an imprint of events that occur as hematopoietic stem/progenitor cells undergo stepwise transformation. It is also worth noting that epigenetic patterning is more tightly linked to differential methylation of CpG shores and not CpG islands, suggesting that alterations of CpG islands are more stochastic in AML and not as reflective of how specific somatic mutations perturb the epigenome. Finally, it is not clear whether the observed enhancer changes are simply the reflection of the transformation process, or whether their differential methylation contributes in any way to maintaining the malignant phenotype. Our data point to the need to carefully map epigenetic marks in experimental models of malignant transformation to explore such questions.
Another striking result of this study was the unexpected epigenetic antagonism between co-occurring IDH1 and DNMT3A mutations. In addition to our cohort, these mutations exhibit co-occurrence in other AML cohorts (The Cancer Genome Atlas: log odds 1.64, P < 0.01; Papaemmanuil and colleagues: FWER < 0.05; refs. 1, 28). Our study allows the direct comparison of the methylomes of IDH-mutant and DNMT3A-mutant AMLs, and we show that these cases contain opposing cytosine methylation profiles with differential methylation occurring at largely distinct CpGs in different genomic compartments. IDH2 mutations result in hypermethylated DMCs that uniformly extend throughout the gene and gene neighborhood compartments. In contrast, DNMT3A mutations result in hypomethylation within the gene body, gene neighborhood, and intergenic compartments at a largely distinct set of CpGs. Co-occurrence of both IDH1 and DNMT3A mutations results in a broad “loss” of differential methylation from both compartments. Given that AML develops in a stepwise manner, these data do not exclude that whichever mutation occurred first did not affect and reprogram the epigenome as per the effect of each mutation. Indeed, in mouse models of IDH or DNMT3A mutation we observe similar effects to what we see in humans in which IDH and DNMT3A are independent of each other (4, 20, 21). It is also important to note that even if cytosine methylation does not appear perturbed, other chromatin marks such as 5′ hydroxymethylation or histone modifications might be affected.
Regardless, our data still point toward functional cooperation between mutant IDH1 and DNMT3A as indicated by gene expression signatures implicating upregulation of KRAS signaling, downregulation of certain MYC targets, and increased in vitro sensitivity to MEK inhibition. This pattern is not the result of simple KRAS mutation in IDH1/DNMT3A cases, which do not co-occur in this cohort. In addition, recent murine models have shown that DNMT3A R882 mutations specifically cooperate with NRAS to promote leukemogenesis via concomitant specific DNA hypomethylation and gain of active histone modifications (29). In this context, we speculate that this finding might indicate a vulnerability of these leukemias to therapeutic targeting of signaling pathways linked to RAS. In summary, this first large-scale study of genome-wide cytosine methylation patterning underlines the power of the epigenome to reveal insights into the mechanism of action of oncogenes and point to putative unexpected actionable targets in specific sets of patients with AML.
Interestingly, the presence of a TET2 mutation did not emerge as a cluster-defining lesion as did IDH2 and IDH1, despite the fact that TET2 and IDH mutations both exert their effect on the epigenome through deregulation of the DNA demethylation pathway (4, 21, 24). This is consistent with prior studies (4, 30) and may be due to compensation by other TET enzymes. Although IDH1 and IDH2 mutations inhibit all TET enzymes, the functional effect of TET2 deficiency alone might produce a more modest phenotype.
The AML patient samples analyzed in this study were collected at the Erasmus University Medical Center (Rotterdam, the Netherlands) from 1990 to 2008 and processed as previously described (31, 32). All samples consisted of >90% blasts. The samples used for analysis were chosen from our prior study using HELP analysis of DNA methylation (2), with 5–10 patient samples selected per epigenetic cluster based on sample availability. Twenty-two normal CD34+ bone marrow specimens were obtained from AllCells. In addition, 6 AML samples were obtained from the Memorial Sloan Kettering Hematologic Oncology Tissue Bank for colony-forming assays. Two cord blood controls were obtained from the New York Blood Center. This research was approved by the institutional review boards at Weill Cornell Medical College, University of Michigan Medical School, Erasmus University Medical Center, and Memorial Sloan Kettering, and written donor informed consent was obtained in accordance with the Declaration of Helsinki.
High-Throughput Methylation Profiling
ERRBS was performed as described previously (22, 33). Briefly, DNA was isolated using standard phenol-chloroform extraction and digested with MspI, followed by end repair, A-tailing, and ligation of methylated Illumina adapters. Fragment lengths between 250 and 400 bp were gel-isolated, treated with sodium bisulfite, and PCR amplified. Libraries were then sequenced on an Illumina HiSeq2000 per manufacturer's recommended protocol for 50 bp single-end read runs. Image capture, analysis, and base calling were performed using Illumina's CASAVA 1.7. Alignment was performed using Bowtie version 0.12.7 and Bismark version 0.5.4 (34) after adapter trimming using FAR version 2.15 as part of the Weill Cornell Epigenomics bisseq pipeline (https://pbtech-vc.med.cornell.edu/git/thk2008/bisseq). Raw aligned reads and methylated base calls for CpGs with 10 to 400× coverage were deposited in GEO (GSE98350).
Gene Expression Analysis
Gene expression analysis for this cohort using Affymetrix Human Genome 133 Plus2.0 GeneChips has been published previously (ref. 32; GEO accession number GSE6891). For this study, raw data were reprocessed using the HGU133plus2.0 BrainArray annotation version 17.0.0 (35), and differential expression was evaluated using limma version 3.22.7. Loci with ≥ 1.5 fold change relative to CD34+ controls and an adjusted P value of ≤ 0.05 were considered significant. Multiple probes mapping to the same locus were summarized using the median probe expression.
Identification and Analysis of Genetic Lesions
Targeted amplicon gene enrichment was performed on all samples and 12 normal controls with a ThunderStorm instrument (RainDance Technologies) using a custom myeloid malignancy gene panel sequenced by the manufacturer. The 260-bp paired end read sequencing was performed using Illumina chemistry (MiSeq v3) in 3 batches with 48 samples multiplexed per lane. Demultiplexed reads were mapped to the human reference genome (hs37d5) using BWA MEM. Post-alignment primer trimming was performed using AmpTools followed by mutation calling using FreeBayes. Variants were annotated using SnpEff. Amplicons displaying a nonzero median allele frequency in controls were dropped, and mutations were called for lesions displaying >5% variant allele frequencies with >10 reads representing the variant allele and having >100× coverage at the locus.
Identification of mutations in NPM1, FLT3, IDH1, IDH2, DNMT3A, CEBPA, and FLT3, overexpression of EVI1, and silencing of CEBPA were performed as described in prior studies (2, 4, 36–41).
To quantify the predictive ability of each mutation relative to translocation data and all other individual mutations, we built multinomial logistic regression models containing (1) each pair of features, (2) the set of translocations alone, and (3) each pair of features adjusted for translocation status. To measure the goodness of fit of each model, we used the Akaike Information Criterion (42). To compare models in (3) to model (2), we calculated the probability of minimizing information loss using exp((AICtranslocation – AICmutation+translocation)/2). This measurement allows the calculation of the probability that the model including mutation data minimizes information loss relative to one containing translocation data alone.
Analysis of Differential Methylation
The base-resolution methylation state of AML and CD34+ control samples relative to the human genome (hg19) was quantified using Bismark (version 0.5.4) after which CpGs with 10–400× coverage were selected. Differential methylation between epigenetic AML clusters and CD34+ controls was identified using a corrected beta-binomial test performed by methylSig (version 0.1.12; ref. 43) with a minimum of 3 CpGs per group required to call differential methylation. DMCs were called for any CpG with a q value ≤ 0.01 and absolute methylation change ≥ 25%. Individual CpGs within genes were associated with the overlapping gene(s) or nearest gene, with promoter > exon > intron > neighborhood > intergenic priority using ChIPseeqer 2.1 (44) and a set of custom PERL and Python scripts (see Supplementary Methods). CpGs falling within the TSS ± 2 kb window were associated with promoters. CpGs falling within the TES ± 2 kb window were labeled as “downstream.” CpGs overlapping any annotated exon were associated with exons. CpGs within introns and not associated with either exons or promoters were associated with introns. Exons and introns were also separated into first exon/first intron and other exon/other intron categories. CpGs falling within the region extending from 2 kb to 50 kb from either the TSS or TES were termed “neighborhoods.” Regions ≥50 kb from the nearest TSS or TES were termed “intergenic.” A 1:1 annotation scheme was used in which every CpG was assigned a single genomic compartment.
Unsupervised Analysis and Generation of Clusters
Hierarchical clustering of the samples was performed using a Euclidean distance measurement applied to the subset of CpGs covered in all AML samples with standard deviation values in the top 15% of those surveyed (n = 142,643). Dendrograms were constructed using Ward's method, and clusters were produced by choosing the highest possible cutoff that did not merge well-characterized translocation-based AML samples with other samples. We then labeled the clusters according to their dominant distinct molecular characteristics, requiring that > 50% of samples in a cluster exhibit the cluster-defining feature.
Analysis of CpG Predictive Ability
The predictive ability of CpGs within a genomic subset was assessed by comparing the clusters produced by CpGs in a given subset to those produced when all CpGs were used. The same 15% most variable CpGs genome-wide were used for each subset and overall. This comparison was quantified using the adjusted rand index (45), and the significance of each annotation's prediction was compared with 500 random shufflings of the CpGs reported by ERRBS in this study using a Z test.
Feature-Centered DMC Analysis
Feature-centered DMC enrichment analysis was performed using a custom set of PERL scripts (see Supplementary Methods). For enhancer and other feature-centered enrichment analyses, a set of virtual bins was created with the central third containing CpGs within the feature and the outer two thirds containing flanking CpGs. For all analyses in this study, 30 bins were used with flanking regions ranging from 50 bp to 10 kb. Significance testing was performed by comparing the number of DMCs falling within active enhancers to those in the 10 kb flanking regions using a Fisher exact test. For analysis of fixed genomic loci, such as TSSs, an analogous approach was used with the 30-bin window centered on the feature position.
IDH1/2 and DNMT3A Cohort Analysis
In order to analyze the interaction of IDH and DNMT3A mutations, patients with IDH1, IDH2, or DNMT3A mutations were identified irrespective of their cluster assignment. These patients were divided into groups with an IDH mutation alone, DNMT3A mutation alone, or an IDH mutation co-occurring with a DNMT3A mutation. Gene expression was analyzed as above by comparing each group to CD34+ normal controls using limma. Identification of differentially methylated CpGs was also performed as above by comparing each group to CD34+ controls using methylSig. Gene set enrichment analysis was performed using GSEA version 184.108.40.206 via comparison to the hallmark gene sets deposited in the Molecular Signatures Database version 5.2 (46).
To identify de novo motifs contained within differentially methylated active enhancers, HOMER version 4.8.3 was used to analyze active enhancers with differentially methylated CpGs. Enrichment was assessed relative to active enhancers with nondifferentially methylated CpGs and filtered for motifs containing strong CG signatures. To identify transcription factor binding site enrichment using experimentally generated ChIP-seq data, we used a subset of the ReMap database (47) generated from myeloid-derived cell lines (K562, CD34 derived, SKNO1, monocyte derived, macrophage derived, U937, APL derived, erythroid derived, hematopoietic stem cell derived, Kasumi1, embryonic stem cell derived, AML derived, hematopoietic progenitor derived, NB4, THP1). We pooled data by transcription factor, and then assessed for enrichment in differentially methylated active enhancers compared with those with nondifferentially methylated CpGs. Enrichment was assessed using the Fisher exact test.
Murine Model of IDH2/DNMT3A Antagonism
All mice were housed at the Memorial Sloan Kettering Cancer Center animal facility with all procedures being approved by the Institutional Animal Care and Use Committee. Rosa26:CreERT2 mice were obtained from The Jackson Laboratory and were crossed to mice harboring a Cre-activatable knock-in Idh2R140Q mutation (27) as well as mice harboring a Cre-activatable knock-in Dnmt3aR878H mutation (26). These crosses generated Idh2 single mutants, Dnmt3a single mutants, and Idh2-Dnmt3a double mutants. Bone marrow was isolated from 8- to 12-week-old mice and transplanted into lethally irradiated CD45.1 hosts as previously described (26). Briefly, mice were euthanized and bone marrow was flushed from femurs and tibia followed by red blood cell lysis with ammonium–chloride–potassium (ACK) buffer. Mononuclear cells (1 × 106) were injected via tail vein into lethally irradiated (950 cGγ) CD45.1 hosts. Engraftment was monitored by peripheral blood chimerism of CD45.2 every 4 weeks. After engraftment, mice were administered tamoxifen via oral gavage (8 mg).
Mice were euthanized at indicated time points for bone marrow isolation from femur and tibia. Cells were ACK-lysed and then lineage depleted (Stem Cell 19756). Cells were incubated for 1 hour at 4°C with the following antibodies (CD34 RAM34 FITC, CD45.2 104 PE, cKIT 2B8 APC, CD16/32 93 AlexFluor700, Sca-1 D7 PE-Cy7, and Lineage cocktail APC-Cy7). Cells were washed with PBS+2%BSA and death was monitored with DAPI. Flow-cytometric analysis was performed on a BD LSR Fortessa and sorting was performed on a BD FACS Aria III at the MSKCC Flow Cytometry Core Facility. Single live, Cd45.2+, Lineage−, Sca-1+, and cKIT+ (LSK) cells were sorted into lysis buffer from the Gentra Pure Blood Kit (Qiagen), and DNA was isolated for ERRBS. DNA methylation differences were assessed using methylSig as above in which the Idh2, Dnmt3a, and Idh2/Dnmt3a were compared to normal LSK controls.
In Vitro Colony-Forming Assays of Primary IDH1/DNMT3A AML Samples
Human cord blood mononuclear cells (CB-MNC) and bone marrow mononuclear cells (BM-MNC) from non–RAS-mutant IDH2, DNMT3A, and IDH1/DNMT3A patients in biological duplicate and technical triplicate (8 total) were thawed and were seeded at 2.0 × 105 cells per well in triplicate into methylcellulose medium (Methocult, H4230, STEMCELL Technologies) with GM-CSF 10 ng/mL and DMSO only, 0.1 μmol/L MEK162, and 1 μmol/L MEK162, obtained from Chemietek. Plates were placed into an incubator at 37°C and 5% CO2 for 14 days after which colony counts were obtained.
Liquid Culture of Primary IDH/DNMT3A AML Samples
Primary human acute myeloid leukemia cells were derived from bone marrow or peripheral blood and blasts and mononuclear cells from three patients with IDH2 mutations without concurrent DNMT3A mutation, three with DNMT3A mutations without concurrent IDH1/2 mutation, two with IDH1/DNMT3A dual mutations, and one with IDH2/DNMT3A dual mutation. Two of the double-mutant samples had IDH1 mutations and one had an IDH2 mutation. Cells were purified by Ficoll–Hypaque (Nygaard) centrifugation and cryopreserved. After thawing, cells were left to recover overnight in Iscove's Modified Dulbecco's Medium (IMDM; Life Technologies) supplemented with 20% BIT 9500 serum substitute (STEMCELL Technologies), 16.7 μg/mL human low density lipoproteins (EMD Millipore), 55 μmol/L beta-mercaptoethanol (Gibco), and a cocktail of recombinant human cytokines (G-CSF 20 ng/mL, GM-CSF 20 ng/mL, IL3 20 ng/mL, IL6 20 ng/mL, FLT3 ligand 50 ng/mL, and SCF 50 ng/mL (STEMCELL Technologies). The next day, 6 × 104 cells were plated in 200 ul of fresh medium in 96-well plates in the presence of either MEK162 (1 μmol/L), MEK162 (0.1 μmol/L), or DMSO vehicle control, all in technical triplicate. On days 2 and 4, fresh MEK162 or vehicle was added. On day 6, viability and cell count were assessed by flow cytometry using staining with 7AAD (BD Biosciences). Data were acquired on a FACS Canto flow cytometer (BD Biosciences) and analyzed using FlowJo software v10.0.8 (FlowJo).
Disclosure of Potential Conflicts of Interest
C.B. Thompson is on the Board of Directors of Merck and Charles River Laboratories; has ownership interest (including patents) in Merck, Agios, Charles River Laboratories, University of Michigan, University of Pennsylvania, and University of Chicago; and is a consultant/advisory board member for Agios. A. Melnick reports receiving commercial research grants from GSK, Roche, Eli Lilly, and Janssen and is a consultant/advisory board member for Epyzyme and Roche. No potential conflicts of interest were disclosed by the other authors.
Conception and design: J.L. Glass, H. Kunimoto, R. Levine, R. Delwel, A. Melnick, M.E. Figueroa
Development of methodology: J.L. Glass, D. Hassane, O.A. Guryanova, M. Fall, A. Alonso, M.E. Figueroa
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.L. Glass, B.J. Wouters, H. Kunimoto, F.E. Garrett-Bakelman, O.A. Guryanova, R. Bowman, S. Redlich, A.M. Intlekofer, M.L. Guzman, P.J.M. Valk, C.B. Thompson, R. Levine, R. Delwel, A. Melnick, M.E. Figueroa
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.L. Glass, D. Hassane, B.J. Wouters, H. Kunimoto, R. Avellino, F.E. Garrett-Bakelman, C. Meydan, T. Qin, P.J.M. Valk, R. Levine, O. Elemento, A. Melnick, M.E. Figueroa
Writing, review, and/or revision of the manuscript: J.L. Glass, R. Avellino, F.E. Garrett-Bakelman, O.A. Guryanova, C.B. Thompson, R. Delwel, A. Melnick, M.E. Figueroa
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.L. Glass, A. Melnick, M.E. Figueroa
Study supervision: J.L. Glass, R. Levine, O. Elemento, A. Melnick, M.E. Figueroa
We would like to acknowledge the Weill Cornell Medicine Epigenomics Core and the Weill Cornell Medicine Applied Bioinformatics Core, which performed sample processing, ERRBS, and initial bioinformatic sample processing. We would like to thank Caroline Sheridan, Tak Lee, and Yaseswini Neelamraju for technical support. J.L. Glass would like to acknowledge Matt Teater for his advice and feedback during the preparation of this manuscript.
J.L. Glass is funded by T32CA009207. O.A. Guryanova is funded by K99CA178191. A. Melnick is funded by LLSSCOR 7006-13 and R01 CA198089. M.E. Figueroa is funded by 1R01HL126947-01 and LLS QFC-3154-14. R. Avellino received funding from the Dutch Cancer Society “koningin Wilhelmina Fonds” and from the Lady Tata Foundation. F.E. Garrett-Bakelman is funded by NCIK08CA169055 and the American Society of Hematology ASHAMFDP-20121 award under the ASH-AMFDP partnership with the Robert Wood Johnson Foundation.
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.