Fractionation and isolation of human blood cells and cell nuclei

Normal polymorphonuclear granulocytes (predominantly neutrophils) were isolated from discarded fresh white blood cells (buffy coats) from unidentified healthy donors collected at Hershey Medical Center Blood Bank (IRB protocol # HY03-136EP-A) using standard OptiPrep density centrifugation [32] and resuspended in PBS (Phosphate Buffered Saline) buffer. The AML blood samples to be used in this study were collected from patients with written informed consent at time of sample collection under guidelines and procedures approved by institutional review board of Milton Hershey Medical Center (IRB protocol 2000–186; STUDY00002518). D. F. Claxton, M.D., the attending physician, was responsible for de-identifying the samples collected from the patients. AML mononuclear cells were collected from peripheral blood of 10 patients diagnosed with AML with and without monocytic features (FAB types 1, 2, 4, and 5) at the time of diagnosis prior to treatment, isolated by Ficoll-Pacque density gradient centrifugation [33], cryopreserved in RPMI media with 10% dimethylsulfoxide and 5% human serum albumin (American Red Cross, Washington, DC, USA) and stored in liquid nitrogen. For clinical classification, the samples were subjected to flow cytometry and percent of surface expression of various antigens is given in S1 Table. Cryopreserved bone marrow CD34+ cells from normal donors were obtained from Allcells (Alameda, CA). Bone marrow CD34+ cells represent a combination of stem and progenitor cells and can be easily isolated at high cell numbers without cytokine mobilization [34]. For a scheme relating the types of AML to the process of human hematopoiesis, see Fig A in S1 File.

Human K562 cells (ATCC CCL-243) were grown in RPMI 1640 GlutaMAX medium (Invitrogen) supplemented with 10% FBS (HyClone, Logan, UT) and 1% penicillin-streptomycin (Invitrogen), in 5% CO₂ at 37^°C. K562 cells were grown for not more than 10 passages after acquiring from ATCC. K562 cells were treated with 0–1 μM UNC0638 (a G9a/GLP HMTase inhibitor; Sigma U4885), with 0.01% (vol/vol) DMSO, for 5 days; medium was replaced on day 3 with fresh growth medium containing UNC0638 or DMSO vehicle only. Cells were incubated with UNC0638 at the concentration indicated, when not specified the K562 cells were treated with 1 μM UNC0638. Three biological replicates were performed, using independently cultured K562 cells. After treatment, K562 cells were washed 2 x PBS and then resuspended in PBS for fixation.

Antibodies

For chromatin immunoprecipitation (ChIP) and western blotting analyses, we have used anti-H3K9me2 (ab1220), anti-H3K9me3 (ab8898), anti-H3 C-tail (ab8898) from Abcam; and Anti-H3K4me2 (07–030) from Upstate/Millipore. We have successfully tested these antibodies for immunofluorescence analysis, western blotting, and ChIP in previous works with myeloid [32] and other cells types [20, 35]. For quality assessment of those antibodies in ChIP applications, see [36].

PAGE and western analysis

SDS-PAGE was carried out in 15% polyacrylamide and the proteins were transferred to Immobilon-P PVDF membranes (Millipore) in Tris-glycine buffer containing 10% methanol. Membranes were blocked in 3% nonfat dry milk in TBS buffer with 0.1% Tween 20, and incubated with primary anti-H3K9me2 (ab1220), anti-H3K9me3 (ab8898), and anti-H3 C-tail (ab8898) from Abcam as previously described [32]. Semi-quantitative analysis of relative protein levels was performed using ECL Prime (Amersham, GE Healthcare) and photographic film, which was digitized and the intensity of protein bands was quantitated using ImageJ [37].

Chromatin immunoprecipitation

For formaldehyde fixation, cultured cells or freshly thawed cells were suspended in PBS at ~ 5x10⁷ cell/ml at room temp., mixed with 1/10 of the final volume of 10% formaldehyde in PBS freshly prepared from 37% stock (Fisher, ACS reagent F79-500). Cells were fixed with 1% formaldehyde for 10–15 min at room temp. Fixation was stopped by adding of 125 mM glycine. Cells were spun down at +4^°C and 4000 g for 5 min. and washed 3 times cold PBS. The fixed cells were counted to estimate the yield and stored overnight at +4^°C.

Nuclei from formaldehyde-fixed cells were isolated as described [38]. The purified nuclei were resuspended in micrococcal nuclease digestion buffer (50 mM Tris-HCl, pH = 7.6; 3 mM CaCl₂) at ~1 mg/ml DNA. The nuclear suspensions (1 ml) were digested with 60 units of Micrococcal nuclease (Roche) for 40 min. at + 37^°C to ~200 bp nucleosome sizes and digestion was stopped by adding 3 mM EDTA. After digestion, the supernatant (S1) was collected and discarded and the nuclear pellet (P) was resuspended in SDS-containing Lysis buffer (1% SDS; 10mM EDTA; 50 mM Tris HCl pH 8.0) with light sonication (Fisher model 100, 2 x 10 sec. at setting 2) to enhance solubility rather than to shear DNA. Chromatin immunoprecipitation (ChIP) was conducted essentially as described [39]. The final DNA was dissolved in 50–100 μl DNase-free water. DNA concentration was measured using Qubit fluorometer (Invitrogen).

For each sample, we conducted two parallel immunoprecipitation reactions with a) no antibody added; b) similarly prepared IgG from a none-immunized animal (e.g. rabbit IgG #2729, Cell Signaling). A typical yield from immunoprecipitation with antibodies against H3K9me2 (Abcam, ChIP Grade, ab1220) was more than 400 ng for 40 μg of the input (ChIP/Input ratio >1%). Nonspecific binding in the control samples was < 0.02%.

Library preparation and Next-generation sequencing

For SOLiD sequencing, the DNA samples (IP ChIP and IP input, ~ 100 ng each) were submitted for preparing fragment libraries and their amplifications using Applied Biosystems reagents and deep sequencing (~20x10⁶ 70 bp DNA sequence reads per sample) using the SOLiD sequencer at the Penn State Genomic Core Facility (University Park). For next-generation Illumina sequencing library preparations, we used NEBNext ChIP-Seq single-end library preparation reagent set for Illumina (cat #E6200) and NEBNext Multiplex Oligos for Illumina Index Primers Set 1 (cat #E7335). For one library, we took 10–30 ng of chromatin-immunoprecipitated or input DNA. End repair of ChIP DNA, dA-Tailing, Adaptor ligation, and cleaning of DNA by AMPure XP magnetic beads (Agencourt part # A63880), were conducted according to the NEBNext ChIP-Seq library preparation manual with exception that after adaptor ligation and AMPure cleaning (step 4), the material was eluted in 20 μl of 0.1 x TE buffer for agarose gel size selection. DNA electrophoresis was run in 2% agarose gels with SyBr Gold gel stain (Invitrogen, S11494) using TAE buffer and 100 bp DNA ladder (NEB N3231S) as size standard. We excised a band from the gel ranging in size from 175–225 bp (150 bp nucleosome plus 50 bp adaptors) using the DNA ladder as a guide. We purified the DNA using Qiagen QIAquick gel extraction kit (cat# 28704) and eluted in 50 μl of sterile water. We took 23 μl of eluted size-selected DNA for PCR enrichment using the universal primer and one of 12 index primers (from NEB cat #E7335) as described in the NEB Next ChIP-Seq library preparation manual. We used 15–17 cycles of PCR amplification depending on the amount of starting material. After PCR, we cleaned the PCR-amplified library DNA by AMPure XP magnetic beads, redissolved in 20 μl of 0.1 x TE buffer and measured the DNA concentration using Qubit fluorometer (dsDNA HS assay). Typical yield was 40–100 ng/μl. For quality controls, we assessed the DNA size distribution on a Bioanalyzer^® (Agilent, High Sensitivity DNA Assay). The DNA libraries were sent for deep sequencing using the Illumina Hi-Seq 2500 sequencer at Penn State Genomic Core Facility (University Park). We used rapid-run mode with 2 lanes, 4 samples per lane and 100 nucleotides per read. For each sample Illumina sequencing provided 32–40 million of 100 bp reads (see S2 Table for SOLiD and Illumina-sequenced samples).

Mapping to the human genome and H3K9me2 domain calling

Sequence reads were mapped to the hg19 assembly of the human genome using bowtie 0.12.9 with default parameters (see GEO database acc.# GSE71809 for the raw and processed sequence reads files associated with this study). Domains were called with the RSEG software, version 0.4.8 RC, using the provided deadzone files k41-hg19 [40]. The size domains were calculated for the RSEG predicted values of 10 kb, 20 kb and 50 kb and appeared to be almost identical (see Fig B in S1 File). The RSEG domains calculated with predicted 10 kb values were used in this study. Other parameter used was: “keep duplicates”. The chromosome X and Y mappings were duplicated for the male samples prior to calling RSEG. Sequence reads mapped to multiple genomic locations were excluded from subsequent analysis. The RSEG enriched domains were mapped into windows of sizes 10 kb, 50 kb, 100 kb, 500 kb and 1Mb, with each domain contributing to a window’s score its average read count for each of its nucleotides falling within the window. Window scores for each dataset were then normalized by the average score of all its windows and technical replicates were then averaged. Supplementary files showing enriched RSEG domains and their boundaries are included in the GEO database (acc.# GSE71809) and contain “rsegout.bed.gz” in their filename extensions.

Hierarchical cluster analysis

Cluster analysis of H3K9me2 sequence read distributions genome-wide and on specific chromosomal loci was done with open source Gene Cluster 3.0 software [41] with hierarchical clustering by complete linkage and visualized either with TreeView version 1.6 (2002) or with Java Treeview 1.1.6r2 (2011).

Correlations and fold enrichment analyses with genomic data from databases

The following human ENCODE [42] datasets were extracted using the University of California, Santa Cruz (UCSC) browser (http://genome.ucsc.edu) table tool: 1–16: chromatin state segmentation by HMM from ENCODE/Broad Institute [43]; 17, 18: European Institute of Oncology/J.C.Venter Institute CD34 DNase I EIO/JCVI CD34+/- Nuclease Accessible Sites, CD34+ mobilized UW DNase I HS hotspots and peaks [44]; 19–24: University of Washington DNaseI clusters, CD34+ mobilized DNaseI hotspots and peaks, and K562 DNaseI hypersensitive peaks [45]; 25: NKI LADs laminB1Lads [46]; 26–28: Transcription Factor Binding Sites for ZEB1, KAP1, and clustered TFBS from ENCODE/Analysis ChIP-seq (40); 29: Target Scan miRNA regulatory sites [47]; 30, 31: Hudson Alpha Institute for Biotechnology K562 Transcription Factor Binding Sites for USF1 [48], 32–39: Broad Institute K562 CTCF, POL2, histone H2AZ, histone H3K27me3, histone H3K4me2, histone H3K9me1, histone H3K9ac [43]; 40–42: Hudson Alpha Institute for Biotechnology K562 CpG methylation by Methyl 450 bead arrays and Reduced Representation Bisulfite Sequencing (RRBS) sites, [49]; 43: CpG islands [50]; 44: vista enhancers (functionally active enhancers) [51]. Additional public datasets used were, 45: H3K9me2 in T lymphocytes [52], 46: cancer SNVs [53]; 47: AML-specific SNVs [11]; 48: Long-range epigenetic silencing in cancer cells, LRES [54]. We also calculated: 49—density of all annotated genes in hg19 per 10 kb windows and 50—density of chromosomal translocation and inversions affecting specific genes per 10 kb windows. The chromosomal translocation density (see S3 Table) was derived from the list of all AML-specific chromosomal translocations and inversions in the Atlas of Genetics and Cytogenetics in Oncology and Haematology (http://atlasgeneticsoncology.org/index.html). For correlation analysis, the RSEG k41-hg19 deadzones [40] were subtracted from each of our experimental, public and computed datasets and then mapped into windows of sizes 10 kb, 50 kb, 100K, 500 kb and 1Mb for chromosomes 1–22 (see S4 Table). Pairwise Pearson correlations were calculated using R 2.15.1.

Identification of conserved and variable H3K9me2 domains

Conserved high (“Red”) and conserved low (“Green”) H3K9me2 domains were determined by calculating the standard deviation of the 10 kb window scores of the H3K9me2 granulocyte, CD34+, and AML samples. Only the male samples were used to determine the conserved domains over the Y chromosome. A threshold of 1 standard deviation was used to classify a sample’s window as ‘high’ (score > stdev) or ‘low’ (score ≤ stdev). An overall window was classified as conserved low (“Green”) H3K9me2 domains when only 2 or less samples were ‘high’ and at most one was ‘high’ from each of the granulocyte, AML, and CD34+ categories. Analogously, a window was classified as conserved high (“Red”) when only 2 or less samples were ‘low’ with at most one ‘low’ from each category. All remaining windows were classified as variable H3K9me2 domains (“Yellow”).

Identification of differential blocks of H3K9me2 (dLOCKs)

To set bases for dLOCK calculations, the 10 kb window RSEG entries were grouped into the following categories: a) granulocytes, b) CD34+ c) K562 control, d) K562-UNC0638 treated, e) AML cluster A, f) AML cluster B. For each group we counted a continuous stretch of 10 kb windows as a LOCK base if at least 50% of 10 kb windows in each row contained positive values. For each of the continuous LOCK bases, for category “c1” above (where c1 is one of a-f), a dLOCK c2>c1 10 kb window score was computed as the sum of the negative log ratio (Sum Log Ratio) of the 10 kb window score of a sample in the first listed category to that of a sample in the second for every possible first and second category sample pair. This result was then normalized by the number of pairwise terms in the sum. The X and Y chromosome dLOCK 10 kb sums only included terms of sample pairs having the same gender.

Contiguous blocks of 10 kb windows were determined to be a positive dLOCK of size spanning those windows if the total score, defined to be the sum of the corresponding dLOCK 10 kb window scores, was positive, and a negative dLOCK for a negative total score. The threshold of significant positive dLOCKs was set to the 90 percentile, and that of significant negative dLOCKs to the 10 percentile. The latter calculation was performed using the ‘quantile’ function in R, with alpha = 0.90 and 0.10 respectively excluding Y chromosome LOCKs. We thus calculated ten categories of differential LOCKs (dLOCKs): 1. CD34+ > granulocytes; 2. granulocytes > CD34+; 3. CD34+ > AML cluster A; 4. AML cluster A > CD34+; 5. granulocytes > AML cluster A; 6. AML cluster A > granulocytes; 7. CD34+ > AML cluster B; 8. AML cluster B > CD34+; 9. K562 UNC0638 > K562 control; 10. K562 control > K562 UNC0638; each category including ALL dLOCKs, and dLOCKs with alpha <0.10, >0.90, <0.05, >0.05, <0.01, >0.99 (S5 Table). For 10 categories of dLOCKs as well as for the 3 genomic zones as well as for each of the Red, Green, and Yellow genomic fractions, we calculated the enrichments of the genomic features from UCSC/ENCODE databases (S6 Table).

Pathway, gene ontology and upstream regulator analyses

For each dLOCK set described above, the genes of the top 1%, 5% and 10% scoring dLOCKs (alpha <0.01 and >0.99, <0.05 and >0.95, <0.10 and >0.90) were extracted (see Table A in S7 Table). Ingenuity® Pathways (IPA build version 212183) was used to identify upstream regulators and significantly enriched canonical pathways, networks, diseases and disorders, molecular and cellular functions, physiological system development and functions, and toxicity lists. Ingenuity Pathway was used with the following parameters: reference set to Ingenuity Knowledge Base Genes Only, relationships to include set to direct and indirect, includes endogenous chemicals, and filter set to only consider molecules or relationships for the human species and where the confidence is set to experimentally observed.

mRNA expression analysis

Total RNA was isolated from cryopreserved AML cells and human bone marrow CD34+ cells from normal donors (Allcells, Alameda, CA) using Trizol® reagent (Invitrogen, 15596–026) with 1 ml of Trizol reagent per 5–10 x 106 cells as described in the vendor’s manual. Isolated RNA was quantified using the Bioanalyzer RNA 6000 Nano Assay kit (Agilent Technologies, Santa Clara, CA). Four AML type A and two AML type B samples randomly chosen from initial 10 AML samples and 6 normal donor samples were taken for this analysis. First strand cDNA was produced from 1.0 μg of total RNA using the standard protocol from the High Capacity cDNA Reverse Transcriptase kit (Life Technologies, Grand Island, NY). cDNA concentrations were quantified by absorbance using a NanoDrop ND-1000 Spectrophotometer (Fisher Scientific, Pittsburgh, PA). cDNA expression rates in AML were assayed using Illumina HumanHT-12 array containing hybridization probes with 33,732 distinct chromosomal positions and using normal bone marrow CD34+ cells as controls (see Table D in S7 Table).

Genomic qPCR, RT-qPCR and PCR primers

All quantitative real-time PCR (qPCR) was performed using FastStart Universal SYBR Green Master (Roche) on a Roche LightCycler 480 Instrument following the manufacturer’s instructions.

For ChIP-qPCR assays, reactions were performed in triplicate, with two technical replicates for each sample. The fold enrichment of ChIP over input was calculated as 2 to the power of the quantification cycle (Cq) difference between input chromatin and ChIP samples. Genomic qPCR primers are listed in S8 Table.

For RT-qPCR, total RNA was extracted using Trizol (Invitrogen) and Direct-zol RNA MiniPrep kit (Zymo Research), and cDNA was synthesized using SuperScript II First-Strand Synthesis System and oligo(dT) (Invitrogen). Reactions for each target were performed in triplicate, with two technical replicates for each sample. Melt curve analysis was performed at the end of PCR to confirm the presence of a single, specific product. RT-qPCR primers are listed in S8 Table. The (Δ(ΔCq)) method was used to calculate target mRNA levels relative to two housekeeping genes, GAPDH and TBP, with similar results; data relative to GAPDH is shown. Fold change in gene expression was calculated by comparing DMSO control- and UNC0638-treated K562 cells. Student’s T-Test was used to evaluate the statistical significance of RT-qPCR gene expression analysis in two independent biological replicates.

Myeloid cells contain long chromosomal blocks of H3K9me2 that have low gene density and are sustained during myeloid differentiation

To determine genome-wide rearrangements of H3K9me2, we conducted ChIP with anti-H3K9me2 and several control antibodies with normal human granulocytes, normal bone marrow CD34+ hematopoietic progenitors, and K562 human myeloid leukemia cell line (see experimental flowchart in Fig A in S1 File). K562 is a blastic myeloid leukemia cell line that has been extensively characterized as a tier I model cell of the human ENCODE project [43] thus facilitating its use in genome-wide data analysis.

Most ChIP-sequencing and mapping experiments were independently repeated on two sequencing platforms (SOLiD and Illumina) each resulting in about 20–30 mln. reads per sample (S2 Table) to ensure that the results of data analysis are independent of the next generation sequencing techniques. The sequence reads (omitting duplicated sequences) were aligned to the human genome hg19 (see GEO database acc.# GSE71809 for the sequence reads files associated with this study). At the chromosomal scale, the H3K9me2 ChIP reads were distributed evenly, with exception of the haploid X and Y chromosomes in males (Fig 1A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Human myeloid cells contain conserved blocks of H3K9me2 at gene-depleted chromosomal regions.

Distribution of DNA sequence reads mapped to the human genome at 1Mb resolution (a) and at 10 kb resolution (showing a region of chromosome 5 between 0 and 5 Mb (b-g). DNA sequence reads were derived from H3K9me2 ChIP-seq of granulocytes sample 15 (a, b), bone marrow CD34+ cells sample 20 (c), K562 cells sample 3 (d); granulocyte input DNA sample 23 (e), granulocyte ChIP-seq of H3 C-tail sample 16 (f), and granulocyte ChIP-seq of H3K4me2 sample 28 (g). The reads were grouped by 10 kb windows and normalized to genome average over the total human genome. Note enrichment of H3K9me2 over a large chromosomal domain, or LOCK [18] indicated by red dashed line box and variable levels of H3K9me2 over the TERT gene forming a peak in the granulocyte and CD34+ cells but not in K562 cells indicated by vertical red arrows. Y-axes represent fold enrichment of sequence reads vs. genome average.

https://doi.org/10.1371/journal.pone.0173723.g001

At the higher resolution, the H3K9me2 ChIP-seq reveals profound variations in the distribution of H3K9me2. A map of a chromosome 5 region in granulocytes at 10 kb resolution (Fig 1B) shows an example of typical genome-wide distribution of H3K9me2 with low levels of this mark associated with a region of high gene density and an extended block of high level of H3K9me2 associated with a gene desert (underlined by a red line). Similar Large Organized Chromatin H3K9me2 domains (LOCKs) were described previously for human placenta cells [18] using ChIP-on-chip technology for selected regions of the human genome. The domain-wide organization of H3K9me2 is in sharp contrast to the even distribution of the input DNA (Fig 1E) and ChIP-seq for total histone H3 using antibodies against its unmodified C-tail (Fig 1F), as well as ChIP-seq showing sharp peaks of an active chromatin mark H3K4me2 (Fig 1G). The input DNA and the H3 C-tail levels appear to be somewhat lower at the regions with high H3K9me2 consistent with the recent observation that micrococcal nuclease accessibility and release of mononucleosomes negatively correlates with two other epigenetic markers of repressed chromatin, H3K27me3 and H3K9me3 [55] and our observation of an even stronger negative correlation of H3K9me2 with MNase sensitivity in K562 cells (see Fig 2B below). Furthermore, the ChIP-seq data normalized to input DNA and histone H3 C-tail also show the same constitutive LOCK in both cell types (Fig C in S1 File) thus attesting for unique and specific genome-wide distribution of the histone H3K9me2 mark in myeloid cells.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Genome-wide correlation analysis of H3K9me2 domains in myeloid cells.

a-c): graphs showing Pearson correlation values for normal granulocytes (a), K562 (b), and CD34+ H3K9me2 (c) ChIP HMM domains vs: input DNA, H3K4me2 ChIP HMM domains, and H3K9me2 ChIP HMM domains from the two other samples as indicted. The five columns in each sample show correlation values for 10K, 50K, 100K, 500K, and 1M windows (left to right). d-f): graphs show Pearson correlations between granulocyte (d), CD34+ (e), and K562 (f) H3K9me2 HMM domains vs. indicated biodata from human genome-wide databases (for references see Materials and Methods and S4 Table).

https://doi.org/10.1371/journal.pone.0173723.g002

The positions of the large H3K9me2 blocks were generally conserved over the chromosomes of granulocytes, CD34+ progenitors, and K562 cells (cf. Fig 1B–1D) and apparently similar to LOCKs described in human placenta [18]. However, we were intrigued by observing some gene-containing chromosomal domains (for example a domain including the TERT gene repressed during granulocyte differentiation and reactivated in many types of cancer cells [56], which showed marked changes between normal myeloid cells and K562 (Fig 1B–1D, and Fig C in S1 File; vertical red arrows). Therefore, we decided to characterize the extent and levels of H3K9me2 domains in the normal and myeloid leukemia cells in relation to other genome-wide biological data available from databases.

H3K9me2 domains in myeloid cells are enriched with lamina-associated domains (LADs) and cancer-related SNVs

To determine the genomic positions and boundaries of H3K9me2 domains and eliminate possible artifacts and repeats known as “deadzones” we used RSEG, an unbiased procedure for selecting domain boundaries based on hidden Markov model (HMM) analysis of genome-wide ChIP-seq read distribution [40]. First, we calculated RSEG domains for H3K9me2 distribution in K562 cells and analyzed co-localization between the boundaries of H3K9me2 domains with other functional genomic features from ENCODE/Broad Institute data [43]. We observed that the sharp transitions in the occurrence of transcription regulation sites (strong promoters and enhancers, and especially the peaks of CTCF, a protein known to demarcate the chromatin domain boundaries [57]) coincided with the boundaries of H3K9me2 domains showing a strong exclusion of these regulatory elements from the H3K9me2 domains. The weaker enhancers and promoters showed much more shallow transitions and hence less exclusion. Notably, H3K4me2, a marker of active genes, had a strong drop coinciding with the H3K9me2 boundaries in a sharp contrast with peaks of the two repressive marks, H3K27me3 and H3K9me3 (Fig D in S1 File).

For genome-wide correlation analysis of H3K9me2 in K562, granulocytes, and CD34+ cells and comparison with our experimental datasets and available data from human genome databases (see refs. 41–55), the RSEG domains were grouped by 10 kb, 50 kb, 100 kb, 500 kb, and 1Mb windows and Pearson correlations were calculated for all experimental and database entries (Fig 2 and S4 Table). Independent analysis of our experimental datasets showed a negative correlation of H3K9me2 with H3K4me2 in the three cell types (Fig 2A–2C) consistent with the boundary analysis. Interestingly, the two types of normal myeloid cells, granulocytes and CD34+ progenitors, both showed a higher positive correlation between themselves than with a leukemia cell line K562 (Fig 2A–2C) prompting us a more extensive analysis of H3K9me2 in human AML cell (below).

Correlation analysis with multiple genomic datasets revealed several general features consistent with the expected function of H3K9me2 in establishing repressive and/or heterochromatic chromatin states: 1) strong positive correlation of H3K9me2 in K562, granulocytes and CD34+ cells with H3K9me2 distribution in human lymphoid T-cells [52] thus showing that this histone modification is generally conserved between myeloid and non-myeloid cells; 2) positive correlations of H3K9me2 with H3K9me3 and H3K27me3 (the two repressive marks analyzed by ENCODE [43] is consistent with synergism between these three main types of repressive epigenetic marks [58, 59]; and 3) negative correlations with a number of activatory epigenetic marks, promoters, enhancers, transcriptional elongation, RNA polymerase II, gene density, and accessibility to DNase I and micrococcal nuclease (Fig 2D and 2E) thus confirming a tight association of H3K9me2 with chromatin states previously described as “repressed chromatin” and “heterochromatin” [43] in myeloid cells.

As expected from previous studies in non-myeloid cells [18, 53, 60], for all three myeloid cell types we observed strong positive correlations between H3K9me2 domains and LADs, total cancer SNV, and AML-specific SNV (Fig 2 and S4 Table). Remarkably, we observed a striking contrast between the strong positive correlation of H3K9me2 with the previously defined “repressed” chromatin state [43] and H3K27me3 in K562 and CD34+ cells and the very weak correlations in granulocytes (arrows on Fig 2). In addition, correlations with the insulator elements known to separate functional chromosomal domains [61] were notably discriminatory between CD34+ and granulocytes. H3K9me2 was strongly correlated with a constitutive heterochromatin mark H3K9me3 in K562 but not in granulocytes and CD34+ cells. No other biomarkers showed such strong differences among the myeloid cells (Fig 2 and S4 Table).

We thus concluded that the H3K9me2 domains form a set of conserved large organized chromosomal blocks (LOCKs) whose association with gene deserts, low nuclease accessibility, and regions of increased occurrence of cancer-related SNV (reflecting their relative inaccessibility to DNA repair [62]), is universal. The genome-wide distribution of H3K9me2 is clearly distinct from the other repressive heterochromatin marks (H3K9me3 and H3K27me3), and changes dramatically between the normal myeloid cells and leukemic K562 suggesting that the H3K9me2 domains may substantially rearrange during normal myeloid differentiation as well as leukemogenesis.

Histone H3K9me2 undergoes genome-wide rearrangement in AML cells

To look for possible alterations of H3K9me2 in human myeloid leukemia, we conducted ChIP-seq analysis of 10 different human AML samples, mostly cases with cytogenetically normal AML (S1 Table). For all AML H3K9me2 domains, we observed strongly positive (r > 0.5) genome-wide correlations between the individual AML cases as well as with the “normal” human samples–granulocytes and CD34+ progenitors (S4 Table). As controls, the SOLiD/Illumina platforms replicate correlations were in the range r>0.671 to 0.805 thus attesting for independence of the correlation analysis from the sequencing platforms.

Comparison of H3K9me2 domains in AML with multiple genomic biodata showed a wide variation in H3K9me2 correlations with two principle chromatin states [43]: insulators and repressed chromatin (Fig 3A and S4 Table). Seven AML samples (AML7-12, 15) displayed a positive H3K9me2 correlation with repressed chromatin and insulators thus closely resembling CD34+ progenitors. In contrast, the other 3 AML samples (AML13, 14, 16) had negative or very weak correlation with insulators and repressed chromatin, more similar to granulocytes than to CD34+ and the other AML samples.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Analysis of H3K9me2 correlations in AML cells and identification of two epigenetically different types of AML.

A: Graphs showing correlations between Granulocyte, CD34+ progenitors, 10 AML samples, and K562 H3K9me2 HMM domains vs. selected biodata with top discriminatory power (insulator, repressed chromatin, and heterochromatin). B: Hierarchical cluster analysis of genome-wide H3K9me2 distribution for 10 AML samples, CD34+, granulocytes, and K562 cells. The insert on top shows the color key and histogram for H3K9me2 enrichment/depletion score.

https://doi.org/10.1371/journal.pone.0173723.g003

Furthermore, a non-supervised cluster analysis of the H3K9me2 HMM domains showed that the AML samples clearly segregate into two distinct clusters (Fig 3B). One cluster designated here as AML cluster A (high level of correlation between H3K9me2 vs. repressed chromatin on Fig 3A) contained all 3 CD34-positive cases of AML with features of less differentiated cells. The other cluster, designated AML cluster B (strong correlation with heterochromatin) showed features of more differentiated AML M4 and M5 and was characterized by high monocytic cell counts (S1 Table).

Identification of large variable blocks of H3K9me2 (dLOCKs) associated with myeloid differentiation and leukemogenesis

In addition to the changes in genome-wide correlations, mapping of H3K9me2 revealed prominent changes in large blocks of H3K9me2 often containing multiple genes. A striking example of such blocks on chromosome 19 showing the strongest H3K9me2 changes between normal and AML cells (Fig 4A) contains dozens of genes encoding zinc finger (ZNF) proteins. An unsupervised cluster analysis of the chromosome 19 (50–59 Mb) locus (Fig 4B) showed the same two distinct clusters as above (Fig 3B), with AML cluster A showing strongly reduced H3K9me2 and AML cluster B showing distributions much closer to that of K562 cells and granulocytes. A 10 kb map of H3K9me2 HMM domains confirmed the differences between dLOCKs of the AML cluster A and cluster B (cf. panels c-h, k vs. i, j, l on Fig E in S1 File) and also showed that the differences revealed by ChIP reflected actual histone H3K9me2 level variations and not the genomic copy number variations since the input DNA maps of normal and AML cells were almost identical (Fig F in S1 File, panels 5–6). In addition, western blotting analysis showed no significant difference between the AML samples in general H3K9me2 levels (Fig G in S1 File, panel a). Furthermore, parallel ChIP-qPCR experiments with gene regions from selected chromatin regions in granulocyte and four AML type A and two AML type B samples randomly chosen from initial 10 AML samples were consistent with the results of ChIP-seq analysis (Fig G in S1 File, panels b-h) though qPCR data exhibited stronger individual variations than ChIP-seq. Therefore, in subsequent work, we focused our analysis of ChIP-seq data derived from AML cluster A, which is better represented in the current dataset.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Mapping of H3K9me2 blocks shows consistent alterations associated with myeloid differentiation and AML.

A: 50K resolution maps of chromosome 19 showing gene density (1) and H3K9me2 ChIP-sequence reads (Bowtie) of human granulocytes (2), normal CD34+ cells (3), and AML myeloblasts (4). Sequence reads were normalized to average and plotted in Log₂ scale. Red arrows show difference in H3K9me2 levels at three loci at chromosome 19. B: Hierarchical cluster analysis of H3K9me2 distribution over the chromosome 19 (50–59.12 Mb) locus for 10 AML samples, CD34+, granulocytes, and K562 cells.

https://doi.org/10.1371/journal.pone.0173723.g004

The observation of large variable H3K9me2 blocks prompted us to conduct genome-wide comparative analysis of large blocks of H3K9me2. Therefore, we calculated differential LOCKs (dLOCKs) of H3K9me2 for six comparison categories: two representing epigenetic changes during normal myeloid differentiation (CD34+ > granulocytes and granulocytes > CD34+) and four representing epigenetic changes between AML and normal cells (CD34+ > AML type A; AML type A > CD34+; granulocytes > AML type A; AML type A > granulocytes). The dLOCKs are ranked by their Sum Log Ratios (see materials and methods) reflecting the extent of the developmentally-regulated and/or AML-related changes in H3K9me2 levels (Table 1). The threshold of significant dLOCKs was set to the 95 and 5 percentile using the ‘quantile’ function in R, with alpha = 0.95 and 0.05 respectively, excluding Y chromosome LOCKs. The list of all H3K9me2 dLOCKs with changes in Sum Log Ratio (alpha < 0.10 and alpha >0.90) is shown in S5 Table.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Main characteristics of transient H3K9me2 blocks (dLOCKs) with strongest variations among myeloid cells.

https://doi.org/10.1371/journal.pone.0173723.t001

Our analysis revealed extended chromosomal domains with average size range between approximately 200 and 900 Kb. Despite the overall genome-wide correlation of H3K9me2 domains with low gene density (Fig 2), the dLOCKs did not show a substantial depletion in the gene density (Table 1 and Fig 5). Analysis of dLOCK gene intersections between the normal and AML cells helped us to answer the question whether the genome-wide changes in H3K9me2 merely reflect an interrupted differentiation process in AML or they include some genes showing gains or loss of H3K9me2 specific for the disease state. We found that, contrary to the former hypothesis that would predict complete intersection of the dLOCKs changed in the disease with those changed during myeloid differentiation, just ~38% of genes from the dLOCK^CD34+>AMLA did intersect with dLOCK^{CD34+>granulocyte} while ~51% of genes from dLOCK^CD34+>AMLA intersected with dLOCK^{granulocyte>AML}. Reciprocally, only ~32% of genes from dLOCK^AMLA>CD34+ intersected with dLOCK^{granulocyte>CD34+} while ~40% of genes from dLOCK^AMLA>CD34+ intersected with dLOCK^{AML>granulocyte} (Fig H in S1 File). We thus conclude that AML cells, exhibit large-scale alterations (both gains and losses) in the levels of H3K9me2 that are either specific to the disease state or are transiently forming during some early pre-leukemia differentiation stages and preserved in AML.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Major genomic features associated with genes located within H3K9me2 domains altered during myeloid differentiation and AML.

Associations (fold enrichment) of selected genomic features within differential H3K9me2 domains altered in myeloid cells. The top panel shows: domains of constant high (red) variable (yellow) and constant low (green) H3K9me2 in all myeloid cells; The bottom panel shows: H3K9me2 dLOCKs maximal and minimal in CD34+ vs. AML cluster A. Data in the top panel are based on dLOCKs at alpha>0.90, 0.95, 0.99 and dLOCKs at alpha<0.10, 0.05, 0.01). Data in the bottom panel are based on dLOCKs at alpha >0.95 and <0.05. For references (as numbered) see Materials and Methods and S6 Table).

https://doi.org/10.1371/journal.pone.0173723.g005

Large blocks of H3K9me2 show AML-specific variations in genomic features related to gene regulation and genomic stability

The observation of the conserved and variable blocks of H3K9me2 (Fig 2) prompted us to analyze their selective enrichment with specific genomic features. First, we designated the conserved high (“Red”), transient (“Yellow”), and conserved low (“Green”) H3K9me2 zones by calculating the standard deviation of the 10 kb window scores of H3K9me2 in the all myeloid samples combined (see Materials and Methods). The Red, Yellow, and Green zones occupy 12.2%, 66.4% and 21.4% of the human genome respectively. For these genomic zones as well as for each of the six categories of dLOCKs, we calculated the enrichments of the genomic features from UCSC/ENCODE databases (Fig 5 and S6 Table).

Consistent with their repressed state, the Red H3K9me2 domains are depleted in all “active chromatin” traits and are notably enriched (≥1.5 fold) in just a few genomic features: repetitive copy number variations (CNV), repressed chromatin and histone H3K27me3 (Fig 5). Among the dLOCKs, the strongest difference between the CD34+ progenitors and AML is found for the “repressed” chromatin and H3K27me3 (Fig 5), in agreement with the genome-wide correlations (Fig 3). These epigenetic marks are strongly (more than 2-fold) enriched in dLOCK^{CD34+>AML A} and depleted in dLOCK^{AML A>CD34+} (Table 1 and Fig 5). We also observed a notable enrichment of dLOCK^CD34+>AMLA for repetitive copy number variations (CNVs). Since the last result could potentially reflect changes in DNA copy numbers rather than ChIP, we additionally analyzed the DNA annotated as repetitive CNV among our samples and observed no significant variation among all primary cells with exception of increased copy numbers in K562 (Fig I in S1 File).

The prominent specific features of the H3K9me2 blocks enriched in AML, as compared to those enriched in CD34+, are the depletion of DNA methylation sites and density of CpG islands as well as the enrichment with vista enhancers (Fig 5). These findings suggest that the loss of histone H3K9me2 dLOCK^CD34+>AMLA could cause de-repression of genes controlled by DNA methylation in AML. In contrast gain of H3K9me2 at vista enhancers within dLOCK^AMLA>CD34+ would cause repression of corresponding genes in AML.

The most universally conserved genomic features enriched in all six types of the H3K9me2 dLOCKs are the AML SNV densities confirming the association of H3K9me2 with regions of increased mutagenesis, apparently due to inhibited DNA repair in heterochromatin [62]. In contrast, the enrichment of genes involved in chromosomal translocations known to cause AML (S3 Table) was prominent only in dLOCK^AMLA>CD34+ (Table 1 and Fig 5) suggesting that, in addition to transcriptional regulation, the rearrangement of large blocks of H3K9me2 in acute myeloid leukemia is linked to chromosomal translocations and genomic instability.

AML-specific blocks of H3K9me2 contain repeated and coordinately regulated sets of genes sharing common upstream transcriptional regulators

Gene ontology analysis using Ingenuity® pathway search showed that the differentiation-related and the AML-related dLOCKs were moderately (1.7–2.7 fold) but very significantly (p<10⁻⁶ –p<10⁻¹⁵) enriched with genes listed under disease or function annotation “acute myeloid leukemia” (Fig H in S1 File). The most prominent among the H3K9me2 blocks depleted in AML A (dLOCK^CD34+>AMLA) are gene-rich clusters containing multiple repeats of ZNF genes (KRAB-ZFP, zinc finger proteins containing a Kruppel-associated or KRAB domain [63]) genes on chromosome 19 and elsewhere in the genome (Fig 4 and the top lanes in Tables C and E in S5 Table). Other examples of large repeated gene families within the dLOCKs depleted with H3K9me2 in AML include repeats of small nucleolar RNA SNORD116 and SNORD115 on chromosome 15 and clusters of protocadherin genes on chromosome 5 (lanes 12 and 32 in supp. Table E in S5 Table).

In contrast to the AML-depleted dLOCKs, the AML—gained dLOCKs, (dLOCK^AMLA>CD34+ and dLOCK^AMLA>Gran) are not enriched with repeated genes and have a smaller magnitude of H3K9me2 variations (Tables D and F in S5 Table). Remarkably, they contain unique genes associated with hematological development and AML such as RUNX1, ETV6, ERG, MECOM and ETS2 among other genes (Table B in S7 Table).

We have used Ingenuity® pathway analysis to identify upstream transcription regulatory molecules associated with coordinated changes in the dLOCKs (alpha < 0.05 and alpha >0.95). We found the regulatory molecule ERG standing at the top within the H3K9me2 blocks enriched in AML (dLOCK^AMLA>CD34+, p<10⁻¹¹). In contrast, in the H3K9me2 blocks depleted in AML the most significant upstream regulators (p<10⁻¹⁰) were CTCF and RAD21 (cohesin), the transcriptional factors that also regulate chromosome domain boundaries. TRIM28 (KAP-1), a common co-repressor of KRAB-ZFP (ZNF) genes was prominent in the H3K9me2 blocks enriched in granulocytes (Fig 6A). Thus, our data show that analysis of genes within differential H3K9me2 domains reveals specific regulatory pathways affected by large-scale epigenetic changes in AML.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Genes within the transient H3K9me2 blocks are controlled by specific upstream regulators and undergo massive transcriptional changes in AML.

A: Top upstream transcriptional regulators associated with total genes within the four types of dLOCKs (a>0.95; a<0.05) shown at the top. B, C: Gene set enrichment analysis [64] with AML and CD34+ control cell populations reveals contrasting association of the H3K9me2 dLOCK^CD34+>AML (B) and AML-repressed genes within dLOCK^AML>CD34+ (C) with genes up- and down-regulated in AML. D, E: Top upstream regulators associated with the AML-expressed genes within dLOCK^CD34+>AML (D) and AML-repressed genes within dLOCK^AML>CD34+ (E).

https://doi.org/10.1371/journal.pone.0173723.g006

AML-enriched blocks of H3K9me2 are associated with massive gene silencing

To examine whether the genome-wide changes in histone H3K9me2 domains correspond to altered expression in the underlying genes, we conducted gene expression array experiments in AML cluster A cells and normal CD34+ bone marrow stem cells (chosen as the normal cells most similar to AML cluster A in their genome-wide H3K9me2 distribution (Fig 3B)). mRNA expression was assayed using Illumina HumanHT-12 array containing 33,732 hybridization probes with distinct expressed sequence tags. Of these, 7346 transcription units showed significant changes (p<0.05) between AML and control (CD34+) samples (Table C in S7 Table).

We used Broad Institute gene set enrichment analysis, GSEA, [64] to determine associations between the genes from the CD34+/AML gene expression arrays and genes within the transient H3K9me2 blocks in AML and CD34+ cells. GSEA revealed a dramatic difference between the enrichment of AML-activated genes in dLOCK^CD34+>AML and enrichment of AML-repressed genes in dLOCK^AML>CD34+ (Fig 6B and 6C and Table D in S7 Table), suggesting that transcriptional repression of multiple genes in AML is associated with the AML-enriched H3K9me2 blocks. For the top-scoring genes associated with dLOCK^CD34+>AML and upregulation in AML, Ingenuity® pathway analysis of upstream transcriptional activators revealed several transcriptional regulators including CTCF and RAF21 (Fig 6D) but not TRIM28 apparently because, although its target KRAB-ZFP genes are enriched in H3K9me2, they are known to resist silencing by G9a/GLP [17].

Remarkably, among the top-scoring transcriptional regulators associated with dLOCK^AMLA>CD34+ and downregulation in AML (Fig 6E) we found ERG as the top transcriptional regulator (P<10⁻⁵) just as it was prominent in the total dLOCK^AMLA>CD34+ genes (Fig 6A). ERG expression and ERG-regulated pathways are associated with stem cell characteristics and poor prognosis of AML [65, 66]. We thus concluded that the ERG-regulated genes associated with dLOCKs and silenced in AML are the most likely to be regulated by H3K9me2 and to fulfill functions contributing to AML progression.

Reversing of the H3K9me2 blocks and associated gene silencing by G9a inhibition

Since G9a/GLP, the histone methyltransferase specific for H3K9me2, is involved in repression of several genes associated with AML [24, 25] and contributes to AML progression [26], we asked if we could reverse formation of AML-specific H3K9me2 blocks in human myeloid cells using a pharmacological G9a inhibitor UNC0638 [67]. First, we treated K562 cells with increasing concentrations of UNC0638 and observed a substantial (>2-fold) decrease in total H3K9me2 level by western blot analysis (Fig 7A). We then conducted H3K9me2 ChIP-seq analysis of UNC0638-treated and control K562 cells as we did for the other myeloid cells. Remarkably, analysis of genome-wide correlations showed a dramatic decrease in correlation with insulators, repressed chromatin, and H3K27me3 induced by UNC0638 (Fig 7B) and recapitulating the trend observed during terminal myeloid differentiation (Fig 3A). Accordingly, analysis of upstream transcriptional regulators revealed TRIM28 and ERG among the top transcriptional regulators in dLOCK^UNC0638>K562 and dLOCK^K562>UNC0638 respectively (Fig 7C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Inhibition of G9a/GLP in K562 decreases H3K9me2 levels and activates genes within AML-enriched transient H3K9me2 blocks.

A: Western blotting probed with antibodies against H3K9me2 and unmodified histone H3 show decreasing levels of H3K9me2 in K562 cells treated with G9a/GLP inhibitor UNC0638. B: Graphs showing correlations between H3K9me2 domains in control and 1 μM UNC0638-treated K562 cells vs. selected biodata with top discriminatory power (insulator, repressed chromatin, H3K27me3, H3K9me3, and heterochromatin). C: Ingenuity pathway analysis showing upstream regulators associated with the dLOCK^{UNC0638>K562control} and dLOCK^{K562control>UNC0638}. D: H3K9me2 levels determined by ChIP-qPCR in control and UNC0638-treated K562 cells for selected gene probes representing dLOCK^CD34+>AML (ZNF274, ZNF544), constitutively high H3K9me2 (NLRP11), constitutively low H3K9me2 (TRIM28), dLOCK^AML>CD34+ (MECOM, ETS1, ERG) and dLOCK^AML>Gran (CDH1). E: RT-PCR analysis of gene expression level in control and UNC0638-treated K562 cells for selected gene probes representing dLOCK^CD34+>AML (ZNF274, ZNF544), constitutively low H3K9me2 (TRIM28), dLOCK^AML>CD34+ (MECOM, ETS1, ERG) and dLOCK^AML>Gran (CDH1). p–values represent Student’s t-test for 2 tailed, unpaired equal variance.

https://doi.org/10.1371/journal.pone.0173723.g007

Next, we tested expression of selected genes residing within the transient H3K9me2 blocks by RT-qPCR in UNC0638-treated and control K562 cells. We chose genes belonging to the two regulatory pathways: TRIM28-regulated KRAB-ZFP (ZNF) genes most represented in dLOCKs depleted in AML but not reduced by G9a inhibition and ERG-regulated genes (ERG, ETS1, MECOM) residing in the dLOCKs gained in AML and strongly reduced by G9a inhibition in K562 cells (Fig 7D). We found that transcription of KRAB-ZFP genes was not affected by G9a inhibition consistent with previous findings [17]. Furthermore, MECOM gene was also not affected despite a significant loss in H3K9me2 thus showing that reduction of H3K9me2 per se is not sufficient to reverse repression of all silent genes in dLOCKs. Probably, this is due to some additional levels of epigenetic regulation, such as DNA methylation, could be established upon histone methylation [68] to prevent gene reactivation. In contrast, ERG-regulation pathway genes, ERG and ETS1, were significantly (p < 10⁻⁸ and <10⁻⁶ respectively) activated by UNC0638 treatment (Fig 7E). Since both of these genes are involved in AML leukemogenesis [69–71], we concluded that G9a/GLP histone methyltransferases could have an important role in restricting expression of these oncogenes by forming the AML-specific H3K9me2 blocks. Furthermore, since ERG is the top upstream regulator of the genes within the AML-enriched blocks (Fig 6A) and is itself downregulated by G9a (Fig 7E), our finding indicates that ERG may act as a part of a regulatory circuit in which its repression by G9a/GLP may lead to a decreased expression of downstream genes positively regulated by ERG and negatively by G9a/GLP-mediated H3K9me2.