Abstract
Aging is associated with functional decline of hematopoietic stem cells (HSC) as well as an increased risk of myeloid malignancies. We performed an integrative characterization of epigenomic and transcriptomic changes, including single-cell RNA sequencing, during normal human aging. Lineage−CD34+CD38− cells [HSC-enriched (HSCe)] undergo age-associated epigenetic reprogramming consisting of redistribution of DNA methylation and reductions in H3K27ac, H3K4me1, and H3K4me3. This reprogramming of aged HSCe globally targets developmental and cancer pathways that are comparably altered in acute myeloid leukemia (AML) of all ages, encompassing loss of 4,646 active enhancers, 3,091 bivalent promoters, and deregulation of several epigenetic modifiers and key hematopoietic transcription factors, such as KLF6, BCL6, and RUNX3. Notably, in vitro downregulation of KLF6 results in impaired differentiation, increased colony-forming potential, and changes in expression that recapitulate aging and leukemia signatures. Thus, age-associated epigenetic reprogramming may form a predisposing condition for the development of age-related AML.
AML, which is more frequent in the elderly, is characterized by epigenetic deregulation. We demonstrate that epigenetic reprogramming of human HSCs occurs with age, affecting cancer and developmental pathways. Downregulation of genes epigenetically altered with age leads to impairment in differentiation and partially recapitulates aging phenotypes.
This article is highlighted in the In This Issue feature, p. 983
Introduction
The world's population is aging. By 2050, more than 1.5 billion people, or roughly double the current percentage of the population, will be 65 years of age or older (1). This change in demographics will require a better understanding of aging-associated features and medical challenges, including the higher incidence of cancer. At the hematopoietic level, aging is associated with an increased rate of isolated cytopenias of unknown significance, loss of adaptive immunity, and an increased predisposition to develop myeloid malignancies such as myelodysplastic syndromes (MDS) and acute myeloid leukemia (AML; refs. 2–5). In addition, aging correlates with an increased incidence of clonal hematopoiesis of indeterminate potential (CHIP) with frequent mutations in the epigenetic modifiers DNMT3A, TET2, and ASXL1, as well as in proteins involved in alternative splicing such as SF3B1 and SRSF2. However, although patients carrying these mutations experience an overall increased likelihood of developing hematologic malignancies, this occurs at a rate of only 0.5% to 1% per year and does not affect every carrier (6–10). Moreover, not all age-related myeloid malignancies occur in the context of CHIP mutations. Thus, additional events must be necessary to initiate malignant transformation in healthy individuals.
Aged hematopoietic stem cells (HSC) display decreased self-renewal, reduced homing ability, and an ultimate myeloid differentiation bias due to reduced lymphoid potential (11–15). This loss of function is due in part to altered DNA damage response and accumulation of DNA damage, increased reactive oxygen species (ROS), and modifications in WNT signaling (16–21). Notably, reprogramming of HSCs through a pluripotent intermediate partially rejuvenates aged murine HSCs, suggesting that epigenetic changes with age contribute to impaired HSC function (22).
In line with this, changes in the epigenome have been reported with normal aging. Aged murine HSCs display regional gains in DNA methylation as well as an increase in the number and breadth of H3K4me3 peaks, suggesting that epigenetic shifts with age may contribute to HSC loss of function (23, 24). In human hematopoietic progenitors, global changes in histone modification levels have been reported (25), but detailed studies into chromatin profiles of aged human HSCs are still lacking. Likewise, although DNA methylation has been profiled in human specimens, these studies utilized either peripheral blood, bulk mononuclear cells (MNC), or unfractionated CD34+ bone marrow cells (26, 27). Thus, it is as yet unknown if the age-related epigenetic differences observed in aged murine HSCs are representative of those that occur with human HSC aging, or whether these epigenetic changes contribute to the increased risk of malignant transformation.
In recent years, there has been increasing evidence for the role of deregulated enhancers and other regulatory elements in hematologic malignancies (28–32). Thus, given the epigenetic changes observed with mouse HSC aging, we hypothesized that aging of human HSCs is associated with changes in the epigenetic landscape that likely affect enhancers and other regulatory elements, and contribute to aged HSC loss of function and increased risk of malignant transformation.
Results
Aged HSCs Display Systematic Changes in Histone and Cytosine Modifications
We hypothesized that epigenetic changes in human HSCs with aging may underlie the age-related functional decline observed in these cells. Therefore, to perform a comprehensive characterization of the epigenetic and transcriptional landscape of these cells, bone marrow MNCs were collected from a total of 41 young (18–30 years old) and 55 aged (65–75 years old) healthy donors. We first evaluated the frequencies of the HSC-enriched lineage-negative (Lin−) CD34+ CD38− (herein HSCe) fraction, granulocyte–monocyte progenitors (GMP; Lin−CD34+CD38+CD45.RA+CD123+), common myeloid progenitors (CMP; Lin−CD34+CD38+CD45.RA−CD123+), and megakaryocyte–erythrocyte progenitors (MEP; Lin−CD34+CD38+ CD45.RA−CD123−) in our cohort. As shown previously, we observed an increase in HSCe frequency with age (P = 1.59e-5) as well as a decrease in the frequency of GMP with an increase in MEP among the Lin−CD34+CD38+ fraction (P < 0.001; refs. 33, 34; Supplementary Fig. S1A–S1C).
To understand the incidence of clonal hematopoiesis in our aging cohort and how this may affect our epigenomic studies, we performed targeted sequencing of a panel of 128 genes. From a total of 59 donors (n = 27 young and 32 aged donors) used to complete all the genome-wide epigenetic and transcriptome profiling, 20 had sufficient material available to also perform mutational profiling. Only one of these donors presented a mutation, which corresponded to a DNMT3A mutation with variant allele frequency (VAF) of 0.12. Thus, due to the low VAF and the overall paucity of mutations, we concluded that it would be unlikely that any epigenetic or gene-expression changes in the bulk HSCe population would be driven by the presence of CHIP-related mutations in this or other donors (Supplementary Table S1).
To capture dynamic age-associated epigenetic changes in human HSCe, we characterized histone and cytosine modifications in young and aged HSCe isolated from healthy donors. We first sought to determine how histone modifications change with normal HSCe aging. Using a low-input protocol for chromatin immunoprecipitation followed by massively parallel sequencing (micro ChIP-seq), we determined the genome-wide distribution of H3K4me1, H3K27ac, H3K4me3, and H3K27me3, on 4 to 7 independent donors per age, per mark. Differential peak calling analysis for each of these four histone marks across the two age groups revealed that aging is associated with significant reductions in H3K4me1, H3K27ac, and H3K4me3, affecting 14,077 (11% of all young peaks), 10,512 (24%), and 18,195 (31%) peaks, respectively (log10 likelihood ratio > 3, absolute fold-change ≥ 1.5). Notably, relatively few regions had an increase in signal intensity of H3K4me1 (n = 229 peaks), H3K4me3 (n = 17 peaks), or H3K27ac (n = 35 peaks) with age (Fig. 1A; Supplementary Fig. S2A and S2B). Importantly, although age-related changes in H3K27ac, H3K4me3, and H3K4me1 were distributed across the genome, these changes were not random, but were instead focalized at specific areas of the genome that were reproducibly affected across the different donors. In contrast, H3K27me3 displayed only 1,564 H3K27me3 peaks changed with aging, consisting of both gains (25% of differential peaks) and losses (75% of differential peaks) of H3K27me3. However, despite fewer H3K27me3 peaks being affected with aging, these peaks displayed the greatest magnitude of change in signal intensity (Supplementary Fig. S2A and S2B; Supplementary Table S2).
Focal histone and DNA methylation alterations with HSCe aging. A, Heat map representation of regions with either loss or gain (log10 likelihood ratio > 3, absolute fold change >1.5) of H3K4me1, H3K27ac, H3K4me3, or H3K27me3 signal in aged HSCe compared with young. The log2(IP/Input) signal is plotted for each replicate, centered on the differential peak ± 5 kb. Each column is representative of an individual donor. B, Functional annotation using ChIP-enrich Gene Ontology Biological Processes for genes annotated to peaks that have reduced H3K4me1, H3K27ac, H3K4me3, or H3K27me3 signal in aged HSCe compared with young. Select significant (FDR < 0.05) annotations are shown. C, Row-scaled heat map of the percent methylation for the 529 differentially methylated regions (DMR; FDR < 0.05, absolute methylation difference ≥ 20%) in aged HSCe versus young HSCe. Each row corresponds to a unique region with differential methylation and each column corresponds to one donor. D, Bar plot representation of significant (FDR <0.05) pathways associated with genes that are differentially methylated in aged HSCe.
Focal histone and DNA methylation alterations with HSCe aging. A, Heat map representation of regions with either loss or gain (log10 likelihood ratio > 3, absolute fold change >1.5) of H3K4me1, H3K27ac, H3K4me3, or H3K27me3 signal in aged HSCe compared with young. The log2(IP/Input) signal is plotted for each replicate, centered on the differential peak ± 5 kb. Each column is representative of an individual donor. B, Functional annotation using ChIP-enrich Gene Ontology Biological Processes for genes annotated to peaks that have reduced H3K4me1, H3K27ac, H3K4me3, or H3K27me3 signal in aged HSCe compared with young. Select significant (FDR < 0.05) annotations are shown. C, Row-scaled heat map of the percent methylation for the 529 differentially methylated regions (DMR; FDR < 0.05, absolute methylation difference ≥ 20%) in aged HSCe versus young HSCe. Each row corresponds to a unique region with differential methylation and each column corresponds to one donor. D, Bar plot representation of significant (FDR <0.05) pathways associated with genes that are differentially methylated in aged HSCe.
To confirm the robustness of these observations and determine the false discovery rate (FDR) for our approach, we performed 100-fold permutation analysis of young and aged ChIP-seq data for each histone modification and determined the FDR to be < 0.05 for all marks except for H3K4me3, for which it was 0.14. Moreover, both Western blot and ChIP-seq analyses of histone 3 (H3) confirm that age-related differences in H3 modifications are not due to a reduction in total H3 levels with age, which is in line with the findings of Cheung and colleagues by mass cytometry (ref. 25; Supplementary Fig. S3A–S3D). Finally, to control for any potential gender effects, differential peak calling was also performed in a sex-dependent manner and showed that even when restricted to only male or female donors, the above observations were still present with significant reductions in H3K4me1, H3K4me3, and H3K27ac, and only minimal changes in H3K27me3, suggesting that age-related reductions in these marks are not sex-dependent (Supplementary Fig. S3E and S3F).
Functional annotation of age-related differences in histone profiles revealed that sites with decreased H3K4me1 were linked to genes involved in myeloid and erythroid differentiation and functions, whereas loss of H3K27ac was linked to genes associated with leukocyte activation, apoptotic signaling, and histone modifications. Changes in H3K4me3 and H3K27me3 were associated with genes involved in developmental processes, including cell fate commitment and the WNT receptor signaling pathway (Fig. 1B; Supplementary Table S2). These findings further indicate that age-related changes in histone modifications are not governed by random changes, but instead represent a systematic change in the epigenetic program of HSCe.
Next, to determine whether cytosine modifications are also reprogrammed during aging, we analyzed DNA methylation landscapes in young and aged HSCe by an enhanced reduced representation bisulfite sequencing (ERRBS) approach that captures approximately 3 million CpGs across the genome (ref. 35; n = 5–7 per age group). We identified 529 differentially methylated regions (DMR; beta binomial test, Benjamini–Hochberg corrected P < 0.05 and absolute difference ≥ 20%),encompassing 2,249 differentially methylated cytosines (DMC) annotated to 748 unique genes (Fig. 1C; Supplementary Table S3). Notably, despite the reduced representation aspect of the ERRBS assay as well as the focal nature of the changes observed, we were still able to capture biologically relevant DMRs associated with aging, which were significantly enriched for pathways relevant to HSC biology such as WNT, cadherin, and cell-adhesion pathways (Fig. 1D).
Enhancer Deregulation Is a Feature of Aged HSCe
To obtain a more comprehensive picture of the epigenetic changes observed and identify groups of genes with similar age-related changes in human HSCe, we performed an integrative analysis across all regions that displayed changes in the histone marks studied. Using k-means clustering, we were able to identify 12 different genomic clusters based on their age-associated epigenetic changes (Fig. 2A; Supplementary Table S4). Genomic annotation of these clusters revealed that clusters A through F were associated with epigenetic changes at promoter regions, whereas clusters G through L were associated with changes within gene bodies and distal intergenic regions.
Decreased H3K27ac at immune and cancer-associated enhancers with age. A, Heat map depicting the 12 clusters identified using k-means clustering of regions with significant (LLR > 3 and absolute fold change > 1.5) changes of H3K4me1, H3K4me3, H3K27me3, or H3K27ac with age (n = 37,058 peaks). The fold change (aged/young) signal is plotted for each histone modification for each peak. Annotation to active and poised enhancers as well as bivalent promoters identified in young HSCe is also shown. B, Heat map of H3K4me1 and H3K27ac signal at the enhancer enriched clusters J–L. The log2(Pooled IP/Pooled Input) signal is plotted for each age group, centered on the differential peak ± 5 kb. Select genes within each cluster are denoted on the right of the heat map. C, Bubble plot representation of select KEGG pathways that are enriched in clusters J–L. The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. D, Example of transcription factor binding DNA motifs enriched in clusters J–L. Motifs with significant enrichment (q-value < 0.05) are denoted with a white circle. E, UCSC genome browser track examples of active enhancers within clusters J–L that overlap transcription factor ChIP-seq peaks. Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding input. Lines below tracks represent transcription factor peaks annotated in publicly available datasets of CD34+ cells. Y, young; A, aged.
Decreased H3K27ac at immune and cancer-associated enhancers with age. A, Heat map depicting the 12 clusters identified using k-means clustering of regions with significant (LLR > 3 and absolute fold change > 1.5) changes of H3K4me1, H3K4me3, H3K27me3, or H3K27ac with age (n = 37,058 peaks). The fold change (aged/young) signal is plotted for each histone modification for each peak. Annotation to active and poised enhancers as well as bivalent promoters identified in young HSCe is also shown. B, Heat map of H3K4me1 and H3K27ac signal at the enhancer enriched clusters J–L. The log2(Pooled IP/Pooled Input) signal is plotted for each age group, centered on the differential peak ± 5 kb. Select genes within each cluster are denoted on the right of the heat map. C, Bubble plot representation of select KEGG pathways that are enriched in clusters J–L. The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. D, Example of transcription factor binding DNA motifs enriched in clusters J–L. Motifs with significant enrichment (q-value < 0.05) are denoted with a white circle. E, UCSC genome browser track examples of active enhancers within clusters J–L that overlap transcription factor ChIP-seq peaks. Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding input. Lines below tracks represent transcription factor peaks annotated in publicly available datasets of CD34+ cells. Y, young; A, aged.
To further characterize these age-related changes in the context of specific regulatory regions, we used the histone profiles from young HSCe to establish the normal landscape of poised enhancers, active enhancers, and bivalent promoters in HSCe. Enhancers were defined as nonpromoter regions marked by both (H3K4me1 > H3K4me3) and (H3K27me3 present) or by both (H3K4me1 > H3K4me3) and (H3K27ac present) for poised and active enhancers, respectively; 2,066 poised enhancers and 11,569 active enhancers were identified in young HSCe. Bivalent promoters were defined as promoter regions with overlapping enrichment for H3K4me3 and H3K27me3, and 3,604 genes with such promoters were identified in young HSCe (Supplementary Table S5). Annotation of the 12 clusters of differentially regulated regions in aged HSCe to these genomic regulatory regions revealed five categories of genomic regions that are altered with age: active transcription start site (TSS) I, active TSS II, bivalent, non-TSS, and enhancer. Clusters J, K, and L, which were defined by nonpromoter regions with loss of H3K27ac in the absence of any H3K4me3 changes, were enriched for regions overlapping with active enhancers (n = 4,124 enhancer peaks; Fig. 2B; Supplementary Table S4). Notably, fewer poised enhancers (n = 225 peaks) were annotated to any of the clusters that showed differential histone enrichment with aging, indicating that this type of regulatory element is less affected by age than active enhancers. Genes annotated to enhancers displaying age-related deregulation were involved in pathways associated with lymphoid and immune signaling as well as myeloid leukemias (Fig. 2C; Supplementary Table S4). Among the genes associated with loss of active enhancers were genes encoding several transcription factors such as RUNX1/2/3, HIF1A, IKZF1, MEIS1, and KLF6, as well as ETV6 and GFI1, which have known tumor suppressor functions. In addition to transcription factors, enhancers affected with aging were also annotated to genes encoding several epigenetic modifiers such as KAT6B, KDM8, CBX7, and DOT1L (Fig. 2B).
DNA motif analysis of the enhancer clusters identified strong enrichment for FLI1, ERG, ETV1, IRF2, and RUNX1 motifs, indicating that the establishment of new epigenetic programs in aged HSCe may be mediated through changes in the function of these transcription factors. We next took advantage of publicly available ChIP-seq datasets to validate binding of these transcription factors to these regions in human CD34+ cells (36) and found that, indeed, enhancers lost with age overlapped with FLI1, ERG, and RUNX1 binding (Fig. 2D and E). Together, these results suggest that enhancer deregulation of hematopoietic transcription factors and immune pathways may be one mechanism contributing to HSC decline with age.
Bivalent Promoters in HSCe Display a Switch to Repression during Aging
Among the three categories of age-associated epigenetic changes that were enriched for promoter regions, two (active TSS I and active TSS II) were defined by reduction of H3K4me3 at non–bivalently marked promoters whereas the third one was annotated to bivalent promoters (Figs. 2A and 3A and B). Active TSS I (clusters A and B) included active promoters marked by the presence of both H3K4me3 and H3K27ac, and predominantly lost the acetylation mark with aging, with a milder loss of H3K4me3. These promoters are enriched in pathways associated with B- and T-cell receptor signaling, apoptosis, cell cycle, mTOR and MAPK signaling, as well as several cancer-related pathways. In contrast, peaks within the active TSS II category (clusters C and D) were marked by the presence of H3K4me3, which was significantly lost with aging. These promoters were enriched almost exclusively in cancer-related pathways (Fig. 3B–D).
Loss of activating histone modifications at promoter regions with age. A, Heat map of H3K4me1, H3K27ac, H3K4me3, and H3K27me3 signals at age-associated clusters that are enriched for active promoters (top and middle) and bivalent promoters (bottom). The log2(Pooled IP/Pooled Input) signal is plotted for each age group, centered on the differential peak ± 5 kb. B, Density plots of the log2(Pooled IP/Pooled Input) signal for the characteristic histone marks for the peaks within the active TSS I (top), active TSS II (middle), and bivalent promoter (bottom) categories. C, UCSC genome browser tracks of genes with altered promoters from the active TSS I cluster (GFI1 and NFATC4) and active TSS II cluster (CARM1 and PER1). Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding Input for ChIP-seq. Light green bars below tracks represent the differential promoter regions identified from the cluster analysis. D, Bubble plot representation of select KEGG pathways that are enriched in the active promoter clusters (clusters A–D). The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. E, Bubble plot representation of select KEGG pathways that are enriched in clusters E and F. The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. F, UCSC genome browser tracks of the HOXC cluster, which contains bivalent promoters within clusters E and F that are altered with age. Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding Input for ChIP-seq. Teal bars below tracks represent the differential bivalent promoter regions identified from the cluster analysis. G, Example of transcription factor–binding DNA motifs enriched in clusters E and F. Motifs with significant enrichment (q-value < 0.05) are denoted with a white circle.
Loss of activating histone modifications at promoter regions with age. A, Heat map of H3K4me1, H3K27ac, H3K4me3, and H3K27me3 signals at age-associated clusters that are enriched for active promoters (top and middle) and bivalent promoters (bottom). The log2(Pooled IP/Pooled Input) signal is plotted for each age group, centered on the differential peak ± 5 kb. B, Density plots of the log2(Pooled IP/Pooled Input) signal for the characteristic histone marks for the peaks within the active TSS I (top), active TSS II (middle), and bivalent promoter (bottom) categories. C, UCSC genome browser tracks of genes with altered promoters from the active TSS I cluster (GFI1 and NFATC4) and active TSS II cluster (CARM1 and PER1). Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding Input for ChIP-seq. Light green bars below tracks represent the differential promoter regions identified from the cluster analysis. D, Bubble plot representation of select KEGG pathways that are enriched in the active promoter clusters (clusters A–D). The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. E, Bubble plot representation of select KEGG pathways that are enriched in clusters E and F. The size of each bubble corresponds to the number of genes within the cluster that are within the given gene set. Highly significant (FDR < 0.05) categories are colored in yellow. F, UCSC genome browser tracks of the HOXC cluster, which contains bivalent promoters within clusters E and F that are altered with age. Tracks are of pooled replicates for each age group, normalized to reads per million and to the corresponding Input for ChIP-seq. Teal bars below tracks represent the differential bivalent promoter regions identified from the cluster analysis. G, Example of transcription factor–binding DNA motifs enriched in clusters E and F. Motifs with significant enrichment (q-value < 0.05) are denoted with a white circle.
Like the active promoter clusters, clusters E and F showed reduction of H3K4me3 with age. However, these clusters were highly enriched for bivalent promoters; 79% (n = 2,171) of genes within clusters E and F had bivalent promoters that were altered with age (Fig. 3A and B; Supplementary Table S4). These clusters were annotated to developmental pathways such as WNT and Hedgehog signaling, as well as cancer-related pathways (Fig. 3E and F), and were enriched for TCF4 and homeobox motifs, including HOXA9 (Fig. 3G). Notably, aging was associated with a marked loss of H3K4me3 at bivalent promoters, with only minimal changes in H3K27me3, thus resulting in an effective switch from bivalency to repression, an epigenetic change comparable to that seen in bivalent genes during malignant transformation (37, 38). Such a switch, although having only minimal impact on the steady-state expression of those genes that are normally expressed at very low levels, changes instead their potential for activation in response to signaling cues.
Human HSCe Experience Age-Related Downregulation of Epigenetic Modifiers and Hematopoietic Transcription Factors
To identify the genes that may contribute to the marked histone and methylation changes we observed in aged HSCe, we performed a supervised analysis of RNA-sequencing (RNA-seq) profiles of young and aged HSCe (n = 10 per age group) and determined that 1,133 genes were differentially expressed with age (Benjamini–Hochberg adjusted P < 0.05 and absolute fold change ≥ 1.5) with 517 genes upregulated and 616 genes downregulated with aging (Fig. 4A; Supplementary Table S6). Age-related expression changes in HSCe were strongly associated with epigenetic reprogramming at these loci. Seven hundred thirty-two genes presented changes at regions within the active promoter clusters A–C, and 290 at regions within the enhancer-enriched clusters J–L (Fig. 4B and C).
Differential gene expression of transcription factors and epigenetic modifiers in aged HSCe. A, Volcano plot of the log2 fold-change (aged/young) gene expression versus the −log10(Padj). Significant downregulated and upregulated genes in aged HSCe compared with young are colored in blue and red, respectively (FDR < 0.05 and fold change < −1.5 or fold change > 1.5, respectively). Select differentially expressed genes are labeled. B, Bar plot depicting the number of age-associated differential genes that also epigenetically deregulated with age at the different genomic regions identified in Fig. 2. The numbers above the bars denote the number of genes that are both differentially expressed and within the epigenomic category. DEG, differentially expressed genes. C, UCSC tracks of pooled replicates for young and aged HSCe at the KLF6 and LMNA loci. Colored bars below the tracks denote differential epigenetic regions as identified in Fig. 2A. D, Bar plot representation of the normalized enrichment score (NES) for select gene set enrichment pathways that are up or downregulated with HSCe aging (n = 10 donors per age group; FDR < 0.05).
Differential gene expression of transcription factors and epigenetic modifiers in aged HSCe. A, Volcano plot of the log2 fold-change (aged/young) gene expression versus the −log10(Padj). Significant downregulated and upregulated genes in aged HSCe compared with young are colored in blue and red, respectively (FDR < 0.05 and fold change < −1.5 or fold change > 1.5, respectively). Select differentially expressed genes are labeled. B, Bar plot depicting the number of age-associated differential genes that also epigenetically deregulated with age at the different genomic regions identified in Fig. 2. The numbers above the bars denote the number of genes that are both differentially expressed and within the epigenomic category. DEG, differentially expressed genes. C, UCSC tracks of pooled replicates for young and aged HSCe at the KLF6 and LMNA loci. Colored bars below the tracks denote differential epigenetic regions as identified in Fig. 2A. D, Bar plot representation of the normalized enrichment score (NES) for select gene set enrichment pathways that are up or downregulated with HSCe aging (n = 10 donors per age group; FDR < 0.05).
Among the genes that were differentially expressed were genes previously reported as deregulated in aged HSCe, such as CDC42, JUN, and FOSL (24, 39). In addition, we identified downregulation of certain splicing factors (U2AF1 and SREK1), as well as age-related changes in key hematopoietic transcription factors, including upregulation of GATA2, GFI1b, and EGR1, and downregulation of HIF1A, HSF1, CBFB, BCL6, and the KLF factors 3, 6, 7, and 10. Finally, we also detected age-related changes in the expression levels of several epigenetic modifiers and corepressors (upregulation of CITED2 and HDAC11, with downregulation of KAT7, KAT8, KDM3A, SETD6, SETD8, and SETD1A). To verify that these results were not biased by donor sex, we restricted the differential gene expression analysis to only female donors (n = 3 young and n = 9 aged) and found that 88% of the age-associated differentially expressed genes were still detected, confirming that the aging gene signature reflects differences in donor age rather than sex (Supplementary Fig. S4A and S4B).
Gene set enrichment analysis (GSEA) revealed downregulation of genes involved in DNA damage response and lymphoid signaling, with concomitant upregulation of inflammatory response and IFN-related gene sets in the aging HSCe fraction. In addition, enrichment for the previously reported aged human HSCe signature from Crews and colleagues (40) was also confirmed in our cohort (Fig. 4D; Supplementary Table S7). To confirm the HSC-specific nature of the aging signature identified, we compared it with a publicly available dataset of tissue-specific aging profiles and found very little overlap with any of the reported signatures (41). This finding indicates that the observed expression changes are not a general signature of aging, but rather they are highly specific for HSCe and clearly distinguishable not only from other tissues, but even from differentiated cells in the peripheral blood (Supplementary Table S8).
Age-Related Epigenetic Changes Are the Result of HSC Reprogramming
To determine whether the observed changes in the HSCe epigenome and transcriptome with age were due to the expansion over time of a previously existing subpopulation within the young HSCe fraction, or whether, instead, these changes stemmed from true reprogramming, we performed single-cell RNA-seq (scRNA-seq) of HSCe from young and aged donors (n = 5 donors < 40 years old and n = 4 donors > 60 years old). As the HSCe compartment contains both long-term HSCs (LT-HSC) and more mature progenitors, we analyzed the composition changes of this population with age. Classification of the HSCe compartment (n = 338 young HSCe and n = 310 aged HSCe) using cellHarmony (42) and gene expression data from 18 reference human bone marrow progenitor cell types from the Human Cell Atlas (HCA; ref. 43) revealed that the majority of young and aged HSCe (78.4% and 87.4%, respectively) corresponded to noncycling HSC (nc-HSC). Within both the young and aged HSCe populations, we identified low frequencies of MEP (3.6% and 3.5%, respectively), as well as multilineage progenitors (0.3% and 0.6%, respectively). Although we observed a decrease in cycling HSCs (12% in young vs. 7% in aged) and lymphoid-primed multipotent progenitors (LMPP; 5.9% in young and 1.6% in aged) with age, these small shifts in percent composition are not sufficient to explain the magnitude of the epigenetic changes detected in bulk HSCe with aging (Supplementary Fig. S5A; Supplementary Table S9). Because the majority of the cells profiled represented nc-HSCs and there was an insufficient number of the other types, we focused the rest of our analysis on nc-HSCs. Differential analysis of these cells identified 364 differentially expressed genes (absolute fold change > 1.5; eBayes t test P < 0.05; FDR < 0.05) with age (Supplementary Fig. S5B; Supplementary Table S9). These included genes that were also differentially expressed at the bulk RNA-seq level, such as EGR1, KLF6, HIF1A, and JUN. HOPACH clustering of all genes differentially expressed in aged nc-HSCs compared with young nc-HSCs identified three clusters, or gene expression programs, that largely separated young and aged nc-HSCs (Fig. 5A; Supplementary Fig. S5B). To determine the frequency of these different gene expression programs in young and aged nc-HSCs and determine whether a detectable subpopulation of young nc-HSCs exists that presents an expression profile comparable with aged nc-HSCs, we calculated expression centroids for the young (n = 1) and aged (n = 2) nc-HSC clusters and then applied cellHarmony to reclassify each nc-HSC individually, using these references. Notably, although nc-HSCs isolated from aged donors were frequently classified as young based on their expression profile (26/271; 9.6%), nc-HSCs isolated from young donors were only exceptionally classified as aged (3/265; 1.1%) with only a subset of young donors showing this phenotype (Fig. 5B). Thus, a defined aged-like population could not be reproducibly identified within the sequenced young nc-HSCs and, although we cannot definitively rule out that the rare aged-like HSCs detected in some of the young donors may contribute to some degree of population shift over time, these findings are compatible with a scenario in which the observed epigenetic and transcriptional changes are largely the result of reprogramming during aging.
Age-associated changes in epigenetic landscapes are driven by reprogramming. A, Schematic of the analysis used to identify and classify aged HSC using scRNA-seq data from young (ages 24–37 years) and aged (ages 64–71 years) HSCe. B, Heat map of young and aged nc-HSCs. Expression centroids were calculated from the three clusters (C1, young in green; C2, aged in purple; C3, aged in red) identified using the genes differentially expressed at the single-cell level (shown in Supplementary Fig. S5B) and used to cluster nc-HSCs. Each column corresponds to an individual cell, and age group of the donor is depicted in the bottom bar.
Age-associated changes in epigenetic landscapes are driven by reprogramming. A, Schematic of the analysis used to identify and classify aged HSC using scRNA-seq data from young (ages 24–37 years) and aged (ages 64–71 years) HSCe. B, Heat map of young and aged nc-HSCs. Expression centroids were calculated from the three clusters (C1, young in green; C2, aged in purple; C3, aged in red) identified using the genes differentially expressed at the single-cell level (shown in Supplementary Fig. S5B) and used to cluster nc-HSCs. Each column corresponds to an individual cell, and age group of the donor is depicted in the bottom bar.
Age-Associated Epigenetic Changes May Predispose for AML
Given that alterations in DNA methylation and deregulation of enhancer elements are observed in age-associated myeloid malignancies such as AML (32, 44–46), we next sought to identify changes in DNA methylation or histone modifications with age that may predispose for this malignancy. Using our published ERRBS methylation data from 119 patients with AML (ages 15–77 years; ref. 32) and our young and aged HSCe samples, we sought to identify aging DMRs that were also differentially methylated in AML, regardless of patient age. We stipulated that if an age-related DMR was aberrantly methylated only in aged patients with AML, it could simply be reflective of a normal aging change, like that of our aged HSCe donors. However, if the DMR is likewise aberrantly methylated in young and middle-aged patients, it would indicate that the DMR is a universal feature of AML and may therefore be part of the process of malignant transformation such that AML in the young also display this change. Thus, starting with our age-related DMRs, we first used k-means clustering to identify groups of DMRs that behave in a similar manner. This resulted in 11 different DMR clusters, 8 of which corresponded to DMRs with hypermethylation with aging and 3 with age-related hypomethylation. A total of 7 DMR clusters (5 with age-related hypermethylation and 2 with hypomethylation (P < 0.05, Mann–Whitney rank-sum test, corrected for multiple testing) were identified that consistently presented the same methylation changes across all AML age groups as detected in aged HSCe, whereas the remaining 5 clusters were either variable in their behavior or consistently displayed methylation changes in the opposite direction (Fig. 6A and B; Supplementary Table S10). These AML-like DMRs were enriched for genes associated with cell-cycle regulation as well as transcription factors such as KLF6, RUNX1, and HOXC4/6 in the hypermethylated clusters, whereas hypomethylated DMRs included RXRA, SKI, and SOCS1 (Supplementary Table S10).
Age-associated epigenetic changes may predispose for AML. A and B, Box plots of the percent methylation for k-means clusters of regions that are hypermethylated (k = 8; A) or hypomethylated (k = 3; B) in aged HSCe compared with young. Percent methylation was calculated for each DMR within the clusters generated by k-means clustering. Arrows denote clusters of DMRs that may be predisposing to AML (P < 0.05, Mann–Whitney rank sum test, and corrected for multiple testing). C, Heat map of H3K27ac signal at peaks within the enhancer category that is altered with HSCe aging. Rows are ordered using k-means clustering (k = 10) performed using young and aged HSCe and AML blasts (n = 71 patients). Bars above heat map denote donor age and cell type. D, Box plots of the log2(H3K27ac/Input) enrichment for peaks (n = 4,931) within clusters 1–7 (C) that showed consistent changes in H3K27ac in aged HSCe and AML compared with young HSCe (P < 0.05, Mann–Whitney rank-sum test, and corrected for multiple testing). E, Select significant (FDR < 0.05) gene ontology biological processes associated with peaks in clusters 1–7 (C). F, UCSC genome browser tracks of the read-normalized fold enrichment (H3K27ac/Input) signal at MEIS1, a gene predicted by k-means clustering to have reduced H3K27ac at its intergenic enhancer. Each track is derived from 1 donor, and 2 biological replicates are shown for each condition.
Age-associated epigenetic changes may predispose for AML. A and B, Box plots of the percent methylation for k-means clusters of regions that are hypermethylated (k = 8; A) or hypomethylated (k = 3; B) in aged HSCe compared with young. Percent methylation was calculated for each DMR within the clusters generated by k-means clustering. Arrows denote clusters of DMRs that may be predisposing to AML (P < 0.05, Mann–Whitney rank sum test, and corrected for multiple testing). C, Heat map of H3K27ac signal at peaks within the enhancer category that is altered with HSCe aging. Rows are ordered using k-means clustering (k = 10) performed using young and aged HSCe and AML blasts (n = 71 patients). Bars above heat map denote donor age and cell type. D, Box plots of the log2(H3K27ac/Input) enrichment for peaks (n = 4,931) within clusters 1–7 (C) that showed consistent changes in H3K27ac in aged HSCe and AML compared with young HSCe (P < 0.05, Mann–Whitney rank-sum test, and corrected for multiple testing). E, Select significant (FDR < 0.05) gene ontology biological processes associated with peaks in clusters 1–7 (C). F, UCSC genome browser tracks of the read-normalized fold enrichment (H3K27ac/Input) signal at MEIS1, a gene predicted by k-means clustering to have reduced H3K27ac at its intergenic enhancer. Each track is derived from 1 donor, and 2 biological replicates are shown for each condition.
Similarly, we took advantage of two publicly available cohorts of AMLs that had been profiled for H3K27ac alone (n = 52; ref. 46), or H3K27ac in conjunction with H3K4me1, H3K4me3, and H3K27me3 (n = 19 for H3K27ac, n = 16 for H3K4me1, n = 19 for H3K4me3, and n = 17 for H3K27me3; ref. 47). Using the categories of regulatory elements that we had identified as being altered with HSCe aging, we examined the enrichment or depletion of the defining histone marks for these categories in AML. Similar to our analysis of the age-related DNA methylation changes, we identified 6 clusters of H3K27ac peaks (n = 4,582 peaks) and 4 clusters of H3K4me1 peaks (n = 2,559 peaks) that had decreased enrichment of these marks at enhancers (Padj < 0.005, Mann–Whitney rank sum test) in normal aged HSCe and AML blasts of all ages compared with young HSCe. In addition, 1 cluster of H3K27ac and 1 of H3K4me1 with enrichment of these marks at enhancers showed a similar pattern of mimicking the age-related increases seen in HSCe in AMLs of all ages (Fig. 6C and D; Supplementary Fig. S6A; Supplementary Table S10). These enhancers were associated with RAS signaling, cell adhesion, immune cell activation, and myeloid differentiation, and included genes such as BCL6B, NFATC2, PML, RXRA, GADD45A, HES1, MEIS1, and DNMT3B (Fig. 6E and F; Supplementary Fig. S6B). In addition, we also observed chromatin alterations at active TSSs in aged HSCe that were also present in AML, with 7 clusters (n = 6,028 peaks) having reduced H3K27ac enrichment at active TSS I, and 4 clusters (n = 2,661 peaks) showing loss of H3K4me3 at active TSS I and II in aged HSCe and across AMLs of all ages. These active TSSs were associated with VEGF, RAS, and p53 signaling, and developmental processes, respectively. Finally, we identified 9 clusters that showed loss of H3K4me3 and 3 clusters with loss of H3K27me3 (n = 4,780 and n = 596 peaks, respectively) at bivalent promoters in aged HSCe and AMLs, which were associated with cell fate and specification (Supplementary Fig. S6A and S6B).
Downregulation of KLF6 Impairs Differentiation and Recapitulates Aging and Malignant Expression Profiles
Because KLF6 is among the genes most strongly downregulated with HSCe aging (Fig. 4A) and the DNA methylation changes seen at the locus are also detected in AMLs across all age groups, we next evaluated the impact of KLF6 downregulation on normal CD34+ cells. KLF6 encodes for the Kruppel-like factor 6, a transcription factor essential for yolk sac hematopoiesis and involved in inflammation and myeloid differentiation (48, 49). We first tested the effect of CRISPR/Cas9-mediated knockout of KLF6 on the phenotype of CD34+ cells (Supplementary Table S11). Knockout of KLF6 resulted in an increase in both total colony numbers and granulocyte–monocyte colonies when plated on methylcellulose, demonstrating that downregulation of this gene in CD34+ induces an increase in myeloid colony–forming potential (P < 0.05; Fig. 7A; Supplementary Fig. S7A and S7B). We further tested the impact of KLF6 downregulation on myeloid and erythroid differentiation potential in liquid culture. After 7 days in myeloid- or erythroid-promoting conditions, CD34+ cells transfected with sgKLF6 had decreased expression of CD11b and CD71, respectively, with increased expression of CD34 compared with control single-guide RNA (sgRNA; n = 5; P < 0.05; Fig. 7B). Persistent robust knockdown of KLF6 at this time point was confirmed using intracellular flow cytometry staining (Supplementary Fig. S7C).
Loss of KLF6 impairs differentiation and leads to expression changes reminiscent of leukemia. A, Colony-forming unit (CFU) assays of CD34+ cells with CRISPR/Cas9 knockout of KLF6. Normalized colony numbers per 500 CD34+ cells plated are plotted for total colony number, granulocyte–macrophage (GM), granulocyte–erythrocyte–macrophage–megakaryocyte (GEMM), and burst-forming unit erythroid (BFU-E). Colony numbers for each biological replicate (n = 5) are normalized to the total colony number for that replicate. B, Representative flow cytometry histograms for 2 donors and contour plots from 1 donor for CD34+ cells transfected with sgCTRL or sgKLF6 and cultured in myeloid (top and middle) or erythroid (bottom) promoting conditions for 7 days. Dot plot representation of the percentage of CD34+ cells (top) CD34− CD11b+ (middle), and CD235a− CD71+ (bottom) cells is also shown (n = 5). Cells transfected with sgCTRL were gated on KLF6+, whereas cells transfected with sgKLF6 were gated on KLF6−. C, GSEA leading edge plots showing the enrichment of the gene sets for genes upregulated or downregulated with HSCe aging in CD34+ cells with sgKLF6 knockout (n = 4 replicates). GSEA was run using a list preranked by the Wald statistic (DESeq2), with the weighted enrichment score. D, Bar plot representation of the normalized enrichment score (NES) for the topmost gene-set enrichment pathways that are upregulated or downregulated with sgKLF6 knockout (FDR < 0.05). P values from paired one-tailed t test, calculated with Prism. ns, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01.
Loss of KLF6 impairs differentiation and leads to expression changes reminiscent of leukemia. A, Colony-forming unit (CFU) assays of CD34+ cells with CRISPR/Cas9 knockout of KLF6. Normalized colony numbers per 500 CD34+ cells plated are plotted for total colony number, granulocyte–macrophage (GM), granulocyte–erythrocyte–macrophage–megakaryocyte (GEMM), and burst-forming unit erythroid (BFU-E). Colony numbers for each biological replicate (n = 5) are normalized to the total colony number for that replicate. B, Representative flow cytometry histograms for 2 donors and contour plots from 1 donor for CD34+ cells transfected with sgCTRL or sgKLF6 and cultured in myeloid (top and middle) or erythroid (bottom) promoting conditions for 7 days. Dot plot representation of the percentage of CD34+ cells (top) CD34− CD11b+ (middle), and CD235a− CD71+ (bottom) cells is also shown (n = 5). Cells transfected with sgCTRL were gated on KLF6+, whereas cells transfected with sgKLF6 were gated on KLF6−. C, GSEA leading edge plots showing the enrichment of the gene sets for genes upregulated or downregulated with HSCe aging in CD34+ cells with sgKLF6 knockout (n = 4 replicates). GSEA was run using a list preranked by the Wald statistic (DESeq2), with the weighted enrichment score. D, Bar plot representation of the normalized enrichment score (NES) for the topmost gene-set enrichment pathways that are upregulated or downregulated with sgKLF6 knockout (FDR < 0.05). P values from paired one-tailed t test, calculated with Prism. ns, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01.
RNA-seq profiling of CD34+ cells after downregulation of KLF6 confirmed a partial recapitulation of the aging HSCe signature with enrichment of both upregulated and downregulated age-related expression changes (Fig. 7C). Furthermore, GSEA of RNA-seq of CD34+ with KLF6 knockout showed strong enrichment of several AML- and cancer-associated gene sets (FDR < 0.05; Fig. 7D; Supplementary Table S6). Together, these findings indicate that age-related loss of KLF6 results in an aberrant differentiation phenotype that may predispose to the development of myeloid neoplasms.
Of note, although many other genes change in a manner similar to KLF6, not all age-related changes may predispose to malignancy in a manner similar to that induced by KLF6 knockout. LMNA, which encodes for the nuclear Lamin A/C protein and is one of the genes mutated in Hutchinson–Gilford progeria syndrome, was also one of the most downregulated genes in aged HSCe, and age-related H3K27ac changes at this locus parallel those seen in AMLs of all ages. Downregulation of LMNA by shRNA in CD34+ cells recapitulated features of aging, which consisted of partial recapitulation of the aging signature and increased colony-forming potential, and variable impact on the differentiation potential, with CD34+ cells from different donors experiencing different degrees of myeloid or erythroid differentiation impairment (Supplementary Fig. S7D–S7G). Thus, we show that both LMNA and KLF6 are epigenetically and transcriptionally deregulated with age, and as expected only age-associated alteration of select genes, such as KLF6, may predispose to myeloid malignancy.
Discussion
In the hematopoietic system, aging is associated with an increase in the prevalence of CHIP and the risk of developing myeloid leukemias (7, 9, 50). Although the contribution of the epigenome to aged HSC dysfunction has been explored in murine models, little has been described in human HSCs. Given that the aging HSC phenotype varies among mouse strains, and that genomic elements such as enhancers are not always conserved among species, studying the epigenome in human cells is of particular importance (51, 52). Here we present, to our knowledge, the first map of regulatory elements, including active and repressed promoters, enhancers, and bivalent promoters in normal human HSCe. In addition, we show that aging results in reproducible genome-wide changes at the epigenomic level, affecting both cytosine and histone modifications and targeting key developmental and cancer pathways, including several known tumor suppressor genes in AML and MDS. This shift in epigenetic program of HSCe with aging was particularly evident at active enhancers, indicating that, similar to what has been observed for cancer (28, 31, 53, 54), age-related dysfunction of human HSCe may also be mediated through the dysfunction of these regulatory elements. As such, many of the putative active enhancers that lose H3K27ac with age were associated with genes involved in lymphoid and immune signaling, opening the possibility that, through a change in their activation potential, HSCe enhancer dysfunction may be one of the mechanisms through which immune function is compromised with aging. Importantly, despite a certain degree of interpersonal variability, age-related epigenetic changes targeted a core set of regions across all donors, which was reproducibly detected despite the relatively small number of donors analyzed, indicating that these changes are unlikely to be stochastic. Moreover, scRNA-seq analysis further confirmed that the observed changes in the epigenetic program are likely not due to the expansion of a preexisting subpopulation within the bone marrow, but rather are the consequence of true epigenetic reprogramming. Thus, the epigenetic characterization within this study will not only allow for the identification of those epigenetic regulatory elements that are vital for maintaining normal human HSC function, but also those with the potential to be targeted to rejuvenate aged HSCs.
Although additional work will be required to determine the mechanism driving age-associated epigenetic reprogramming, a possible mechanism that may explain some of the observed epigenetic changes with age is the altered gene expression or splicing of epigenetic modifiers with age, such as the downregulation of the histone 3 lysine 4 methyltransferase SETD1A, which may explain the reduced number of H3K4me1 and H3K4me3 peaks observed in aged HSCe. Because we also observed gene-expression alterations in KAT5/7/8, KDM2A/3A, SETD6/8, and HDAC11, it is possible that HSCe aging may also present focal alterations in other histone modifications, as has been observed at the global level by others (25). In addition, as alterations in RNA splicing have been observed with HSC aging (40), and genes encoding the U2AF1 and SREK1 splicing factors become downregulated with age, it is possible that aberrant alternative splicing of specific epigenetic modifiers with aging may be contributing to age-associated epigenetic reprogramming. However, whether differential expression or usage of epigenetic modifiers leads to alterations in the epigenome or instead the reverse is true is yet to be determined.
In addition, genes encoding multiple transcription factors become epigenetically deregulated with aging, including FLI1, RUNX3, BCL6, MEF2C, and KLF6. Furthermore, regions with age-related epigenetic changes are enriched for motifs of many of these same transcription factors. Therefore, it is possible that changes in the binding of these transcription factors drives epigenetic reprogramming with age by altering the recruitment of epigenetic modifiers to these target sites.
Notably, a subset of the regions that show chromatin changes with normal aging are also altered in both elderly and young AML. These regions include enhancers that putatively regulate immune and cell adhesion–associated genes, as well as active promoters annotated to genes involved in RAS and VEGF signaling, and bivalent promoters associated with developmental genes. Likewise, a subset of age-associated DMRs involved in cell-cycle regulation, apoptosis, and WNT signaling also appear to be aberrantly methylated in all AMLs. These findings would suggest that alterations in DNA methylation and enhancer regulation with age may also contribute to leukemogenesis. Among these age-associated DMRs, we identified KLF6, which not only is among the most strongly downregulated transcription factors with aging, but is also aberrantly methylated in AMLs. Notably, downregulation of KLF6 in normal CD34+ cells was sufficient to change the differentiation phenotype of the cells and recapitulate expression profiles seen in AMLs, indicating that deregulation of these transcription factors with aging is likely to contribute to both aging HSCe dysfunction and the increased risk of malignant transformation. However, given that not all people develop myeloid malignancies with age, it is likely that age-associated epigenetic deregulation alone is not sufficient for transformation, but rather other cooperating events are also required. Future functional screens in the context of known clonal hematopoiesis or leukemic drivers will also be helpful for identifying which age-associated epigenetic changes are most vital for increasing risk of myeloid malignancy with age.
Although our study focused on intrinsic characteristics of HSCe aging, extrinsic signals from the bone marrow niche should not be ignored. It is possible that the changes in the epigenetic program of aged HSCe are the result of microenvironment cues from an aged bone marrow niche, which becomes more adipogenic, with less bone cell formation, and several studies have demonstrated the influence of the aged niche on HSC function (55–57). However, the precise mechanism by which the bone marrow niche may contribute to HSC epigenetic reprogramming with aging remains to be elucidated.
In summary, human HSCe experience dynamic age-related changes in multiple layers of the epigenome that result in a shift in the epigenetic program of aged HSCe likely to affect their functional capabilities. Most notably, scRNA-seq of HSCe points to a true reprogramming mechanism underlying these age-related changes, as opposed to an expansion of preexisting subpopulations already present in the young bone marrow compartment. However, why cells would undergo age-related reprogramming and whether this reprogramming inevitably results in detrimental consequences for the cells or instead forms part of a protective mechanism to help compensate for other deleterious changes in the cells and the microenvironment is yet to be answered. Moreover, whether remodeling of the epigenome with age favors the development of age-related clonal hematopoiesis (ARCH)/CHIP is yet to be determined. Given the low incidence of ARCH in our cohort and that reproducible age-associated epigenetic changes were observed regardless of the presence of ARCH-accompanying mutations, it is possible that the remodeling of the epigenome with age may precede the development of clonal hematopoiesis. Further functional studies utilizing models of clonal hematopoiesis and single-cell experiments will be required in the future to help establish this as well as to discern how any deleterious changes contribute to driving the aged HSCe phenotype and predisposition to myeloid malignancies.
Methods
Human Bone Marrow Samples
All aged (ages 65–75 years) samples were derived from femurs from individuals undergoing hip replacement surgery at the University of Michigan Medical School (Ann Arbor, MI) or the University of Miami Miller School of Medicine (Miami, FL). The study was approved by the Institutional Review Boards of the University of Michigan Medical School and University of Miami Miller School of Medicine. All material included in this study was deidentified, leftover material to be discarded after pathologic review of the hip specimens, and was therefore exempt from human subject research regulations. Donors had no known history of hematologic malignancy. Bone marrow (BM) MNCs were isolated using Ficoll-based density centrifugation and cryopreserved. Additional BM MNCs were purchased from Allcells, Hemacare, StemCell Technologies, University of Pennsylvania Stem Cell and Xenograft Core, and Cincinnati Children's Hospital Processing Core for young (ages 18–30 years) donors. Supplementary Table S1 describes donor information in detail.
FACS Isolation of HSCe
Miltenyi MACS magnetic bead purification was used to enrich BM MNCs for CD34+ cells. The CD34+ fraction was then stained with CD2-PE-Cy5 (eBioscience, clone RPA-2.10), CD3-PE-Cy5 (eBioscience, clone UCHT1), CD4-TC PE-Cy5 (Invitrogen, clone S3.5), CD7-TC PE-Cy5 (Invitrogen, clone CD7–6B7), CD8-TC PE-Cy5 (Invitrogen, clone 3B5), CD10-PE-Cy5 (eBioscience, clone eBioCB-CALLA), CD11b-PE-Cy5 (eBioscience, clone ICRF44), CD14-TC PE-Cy5 (Life Technologies, clone TuK4), CD19-PE-Cy5 (eBioscience clone HIB19), CD20-PE-Cy5 (eBioscience, clone 2H7), CD56-TC (Invitrogen, clone MEM-188), Glycophorin A-PE-Cy5 (BD, clone GA-R2), and Thy1-biotin (eBioscience, clone 5E10, followed by CD34-APC (BD Biosciences, clone 8G12), CD38-PE-Cy7 (eBioscience, clone HIT2), CD123-PE (BD Biosciences, clone 9F5), streptavidin-APC-Cy7 (Life Technologies), CD45RA-FITC (Invitrogen, clone MEM-56), and DAPI. HSCe (DAPI−, Lin−, CD34+, CD38−) were cell sorted using either a BD FACSAria I with a 70-µm nozzle or a BD FACS SORP Aria-Ilu.
ChIP-seq
HSCe were sorted into 1 mL Iscove's Modified Dulbecco's Medium (IMDM) 20% FBS. For H3K4me1, H3K4me3, and H3K27me3, 14,000 to 20,000 HSCe were used per immunoprecipitation; for H3K27ac, 15,000 to 35,000 cells were used per immunoprecipitation, and for total H3, 10,000 to 15,000 cells were used per immunoprecipitation. ChIP-seq was then performed using the True MicroChIP (Diagenode, #C01010130) kit and antibodies that had been validated for specificity and reactivity using the MODified Histone Peptide Array (Active Motif, #13001). The manufacturer's protocol was followed using the following modifications. After quenching with glycine and washing with PBS, samples were suspended in 100 µL undiluted Lysis buffer with 1× Diagenode protease inhibitor cocktail and 5 mmol/L sodium butyrate per 10,000 cells. Samples used for H3 were sonicated in 1.5 mL TPX tubes in a Bioruptor Pico for 6 cycles of 30 seconds on and 30 seconds off. All other samples were sonicated in a Bioruptor XL for 55 cycles of 30 seconds on and 30 seconds off. Chromatin was immunoprecipitated for 12 hours at 4°C using 1 µg H3K27me3 (Millipore 07-449, lot #21494165), 1 µg H3K4me3 (Abcam ab8580, lot #GR164207-1), 0.5 µg H3K4me1 (Diagenode C15410194, lot #A1862D), 0.5 µg H3K27ac (Abcam ab4729, lot #GR155970-2), or 0.5 µg H3 (Abcam ab10799, lot #GR275925-1). After reverse cross-linking, DNA was purified using the minElute PCR Purification Kit (Qiagen, #28004) and eluted in 16 µL of Tris-HCl pH 8.0. Enrichment was verified using qPCR with the primers GAGAGTCCTGGTCTTTGTCA and ACAGTGCCTAGGAAGGGTAT for H3K4me1 and H3, AGGGAGGGAATTAATCTGAG and ACAGTGCCTAGGAAGGGTAT for H3K4me3 and H3, TACTTGGTTTCTGCATCCTT and TCACTAAAGAAACCGTTCGT for H3K27me3 and H3, and GAGCAGAGGTGGGAGTTAG and TCAGACCCTTTCCTCACC for H3K27ac. The remaining DNA was then used for library preparation with the V1 MicroPlex Library Preparation kit (Diagenode, #C05010011). For the PCR amplification, a total of 16 amplification cycles was used. Libraries were purified using a 1:1 Ampure bead cleanup and eluted in 16 µL of Tris-HCl pH 8.0. Fold enrichment over input was then verified using qPCR (primer sequences in Supplementary Table S11). Multiplexed libraries were sequenced either on an Illumina NextSeq 500 or a HiSeq-2500 sequencer.
ChIP-seq Alignment
FastQC was used to evaluate library quality (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Illumina adapters were trimmed using Cutadapt (version 1.12; ref. 58). Reads were aligned to hg19 using Bowtie (version 2.0.5; ref. 59). All samples had at least 28,000,000 aligned reads. Only reads mapping to unique locations were retained for downstream analysis.
ChIP-seq Analysis and Differential Peak Calling
Peak calling for individual IPs was performed using the callpeak function from the model-based analysis of ChIP-seq 2 (MACS2 v.2.0.10.20131216) program (60). For visual representation and differential peak calling, aligned reads for each mark for each age group were pooled using samtools v1.3.1 (61). Peaks were re-called using the pooled IP and the corresponding pooled Input with the nomodel option. For H3K4me1, H3K4me3, and H3K27ac, a q-value cutoff of 0.0001 with the narrow peak option was used. Broad peak calling with a q-value cutoff (–broad-cutoff) of 0.1 was used for H3K27me3. The IP and Input bedgraph files produced by macs2 peak calling were used for differential peak calling with the macs2 bdgdiff function, which adjusts for sequencing depth. To determine the fold change for significant (log10 likelihood ratio > 3) differential peaks that were reduced with age, the fold enrichment of (Pooled IP/Pooled Input) was calculated for each histone modification and age group and normalized to reads per million. The deepTools2 multiBigwigSummary function (62) was then used to determine the number of counts for each differential peak, for each pooled sample. The fold change of Aged/Young was then calculated with gtools (63). Significant peaks for each age group and differential peaks (FDR < 0.05, log10 likelihood ratio > 3, and absolute fold change > 1.5, respectively) were annotated to the nearest TSS of Refseq hg19 using the R packages ChIP-enrich, and Broad-enrich for H3K27me3 (v1.8.0; ref. 64). For the permutation analysis, random pools composed of young and aged samples were generated using R for H3K4me1, H3K4me3, H3K27ac, and H3K27me3, with n = 5, 6, 4, and 6 samples per pool, respectively. Peak calling and differential peak calling were performed as described above (n = 100 pair-wise comparisons per modification, respectively). The FDR was calculated for each histone modification using the significance analysis of microarrays (SAM) method (65). Briefly, for each modification, the median number of differential peaks across all permutations was divided by the total number of peaks identified in young and aged HSCe by macs2 bdgdiff.
Clustering of Aging Histone Changes
k-means clustering was performed using all peaks that had significant (LLR > 3, absolute fold change > 1.5) H3K4me1, H3K4me3, H3K27ac, or H3K27me3 changes with age (n = 37,058 peaks) and R version 3.5.1. For the input matrix, scores for each histone modification, for each age group, were calculated using deeptools2 multiBigwigSummary (62) with bigwig files that were read normalized and contained the fold enrichment of the Pooled IP/Pooled Input. The fold change of aged/young was calculated for each modification for each peak using gtools (63). Peaks with low counts for a given histone modification (score aged <3 and score young <3) were filtered by converting the fold change to zero for that respective mark. K-means clustering of the fold change values was then performed using base R, with set.seed(123), centers = 12, and nstart = 50. Genomic annotation of each peak was performed using the R-package Genomation (v1.2.2; ref. 66), with promoter regions being defined as ± 3,000 bp of the TSS. Genomation was also used to annotate peaks to bivalent promoters and active enhancers identified in young HSCe. The heat map was plotted using ComplexHeatmap (67), splitting between each cluster.
Comparison to AML ChIP-seq
ChIP-seq data from AML blasts was downloaded from the Sequence Read Archive (SRP103200; ref. 46) and the European Genome–Phenome Archive (EGAD00001002340 and EGAD00001002418; ref. 47) and processed as described above. In total, there were 71 patients for H3K27ac, 16 for H3K4me1, 19 for H3K4me3, and 17 for H3K27me3 available for download and used in the analysis. Only patients for whom donor age was annotated were kept for downstream analysis. Deeptools2 bamCompare (62) was then used to create log2-normalized bigwig files for each sample and for young and aged HSCe samples with the following parameters: bamCompare –bamfile1 ip –bamfile2 input –outFileName name_log2.bw –ratio log2 -p 30 –verbose –outFileFormat bigwig –scaleFactorsMethod readCount. For each category of interest (enhancer, active TSS I and II, and bivalent promoter), a matrix containing the counts for each region within the category was generated using the bigwig files for the relevant histone modification and deeptools2 multiBigwigSummary. H3K27ac samples from both Yi and colleagues and McKeown and colleagues were used together. The optimal number of k-means clusters was calculated using the gap statistic method in factoextra (68). K-means clustering was performed using R, and the heat map was plotted using ComplexHeatmap (67), splitting between each cluster, and using light gray for NA values. Boxplots were plotted using ggplot2 (69). Clusters that had significantly different enrichment for the relative histone mark with age and AML were identified using the one-sided Mann–Whitney rank-sum test (implemented in python 3). Only those clusters where young HSCe was significantly (P < 0.05 and Bonferroni Padj < 0.005) different from all other conditions were considered to be predisposing to AML.
Enhanced Reduced Representation Bisulfite Sequencing
FACS-isolated HSCe (n = 5–7 per age group) were sorted into Qiagen RLT+. DNA was extracted using the Allprep DNA/RNA Micro Kit (Qiagen, #80204) and 2 to 10 ng of DNA per sample was used for library preparation. ERRBS was performed as described previously, with the following modifications to accommodate the low input amount of DNA (35). Prior to adapter ligation, the methylated adapters were diluted to 150 nmol/L, and for the gel size selection, fragments of 150 to 450 bp were excised. Libraries were sequenced on an Illumina HiSeq 2500.
Differentially Methylated Regions Analysis
Reads were aligned against a bisulfite-converted human genome (hg19) using Bowtie and Bismark (versions 0.12.8 and 0.4.1, respectively; ref. 70). After filtering and normalizing by coverage, a methylDiff object containing the methylation differences and locations for cytosines that were present in at least three samples per age group (meth.min = 3) was identified using MethylKit (version 0.9.4; ref. 71) and R statistical software (version 3.2.1, R-Core Team, http://www.R-project.org/). DMCs were then input into the R package eDMR (v.0.6.3.1; ref. 72), which was used to identify significant DMR (regions with at least 3 CpG, 2 of which are DMC, with absolute mean methylation difference ≥ 20, and DMR q-value ≤ 0.05). DMRs were annotated to hg19 using the following parameters: DMRs overlapping a gene were annotated to that gene, intergenic DMRs were annotated to all neighboring genes within 50 kb, and if no gene was present within a 50-kb window, DMRs were annotated to the nearest TSS (73). Functional annotation was performed using DAVID version 6.7 (74, 75). For the heat map, the percent methylation for each DMR was calculated for each sample. The row-scaled values were than plotted with heatmap.2, using the complete clustering method with the Euclidian distance, and using light gray for NA values.
Comparison to AML DNA Methylation
The percent methylation for each age-associated hyper- or hypo-DMR was calculated for each young and aged HSCe donor and patient with AML (32). NA values were removed prior to clustering. The optimal number of k-means clusters was calculated using the gap statistic method in factoextra (68). K-means clustering was performed using R. Boxplots plotted using ggplot2 (69). Clusters that had significantly different methylation with age and AML were identified using the one-sided Mann–Whitney rank-sum test (implemented in python 3). Only those clusters where young HSCe methylation was significantly (Bonferroni Padj < 0.00625 for hyper and Padj < 0.017 for hypo) different from all other conditions were considered to be predisposing to AML.
RNA-seq
FACS-isolated HSCe were sorted into Qiagen RLT+ (n = 10 biological replicates per age group). RNA was immediately isolated using the Allprep DNA/RNA Micro Kit (Qiagen, 80204) using the manufacturer's protocol to extract total RNA. Samples with a RIN > 8.5 were used for library preparation. Ribosomal RNA was removed using RiboGone (Clontech, #634846). Stranded libraries were prepared by the University of Michigan Sequencing Core using the SMARTer Stranded RNA-seq Kit (Clontech, #634836). Libraries were sequenced on the HiSeq-2500 with 50 bp paired-end sequencing.
RNA-seq Alignment
Using Cutadapt (version 1.12), all reads were trimmed to 48 base pairs, and adapters were removed (58). Reads were aligned to the hg19 gencode v19 reference genome using the STAR aligner (version 2.5.2b), specifying the following parameters: outFilterType = BySJout, outFilterMultimapNmax = 20, alignSJoverhangMin = 8, alignSJDBoverhangMin = 1, outFilterMismatchNmax = 999, alignIntronMin = 20, alignIntronMax = 1000000, alignMatesGapMax = 1000000, alignEndsType = EndToEnd (76).
Differential Gene Expression Analysis
Gene counts were calculated using QoRTs (v1.0.7; ref. 77). QoRTs was run in second stranded mode using the hg19 gencode annotation file without entries for ribosomal RNA. Differential gene expression analysis was performed using DESeq2 v1.10.1 (78). A multifactor design was used to control for sex of the donor as well as any batch effect during library preparation. Dispersions were calculated using samples from both age groups and then contrasts were established for pair-wise comparisons. Significant genes were defined as having a fold change > 1.5 and Padj < 0.05.
GSEA
GSEA was performed using a list of genes preranked by the Wald statistic (stat column from DESeq2 output; ref. 79). A weighted enrichment score was used and gene set size was limited to 15 to 500 genes. To test enrichment for the Crews and colleagues aging signature, the published list of genes upregulated in aged HSCe (FPKM > 1,P < 0.05, L2FC > 1) was used as a gene set in GSEA (40).
Cell Isolation for scRNA-seq
Mononuclear cells were isolated from 5 young (ages 24–37 years) donors and 4 aged (ages 64–71 years) donors with Ficoll-Paque (GE Healthcare), then enriched using CD34 MicroBeads on an Automacs Pro separator (Miltenyi Biotec). CD34+ cells were stained with lineage antibodies conjugated to PerCP-Cy5.5: CD2 (Becton, Dickinson and Company, clone RPA-2.10), CD3 (Becton, Dickinson and Company, clone UCHT1), CD4 (Becton, Dickinson and Company, clone RPA-T4), CD7 (Becton, Dickinson and Company, clone M-T701), CD8 (Becton, Dickinson and Company, clone RPA-T8), CD10 (Becton, Dickinson and Company, clone HI10a), CD11b (BioLegend, clone ICRF44), CD14 (Becton, Dickinson and Company, clone M5E2), CD19 (eBioscience, clone HIB19), CD20 (eBioscience, clone 2H7), CD24 (BioLegend, clone ML5), CD56 (Becton, Dickinson and Company, clone B159), CD66b (Becton, Dickinson and Company, clone G10F5), Glycophorin A (BioLegend, clone HIR2), and CD90-biotin (eBioscience, clone 5E10). To isolate HSCe, lineage-stained cells were stained with streptavidin APC-Cy7 (Becton, Dickinson and Company), CD123-PE (eBioscience, clone 6H6), CD34-APC (eBioscience, clone 4H11), CD38-PE-Cy7 (eBioscience, clone HIT2), CD45RA-FITC (Invitrogen, clone MEM-56). Cell sorting was performed on a BD FACSAria II with a 100-µm nozzle.
scRNA-seq
Single-cell Lin−CD34+CD38–isolated cells were captured using the C1 Single-Cell Auto Prep System (Fluidigm), according to the manufacturer's instructions. In short, flow-sorted cells were counted and resuspended at a concentration of 35,000 cells per 100 µL PBS, then loaded onto a primed C1 Single-Cell Auto Prep Integrated Fluidic Chip for mRNA-Seq (5–10 µm). After the fluidic step, cell separation was visually scored, between 50–84 single cells were normally captured. Cells were lysed on ChIP and reverse transcription was performed using Clontech SMARTer Kit using the mRNA-Seq: RT + Amp (1771 ×) according to the manual. After the reverse transcriptase step, cDNAs were transferred to a 96-well plate and diluted with 5 µL C1 DNA Dilution Reagent. cDNAs were quantified using Quant-iT PicoGreen dsDNA Assay Kit (Life Technologies) and Agilent High Sensitivity DNA Kit (Agilent Technologies). Libraries were prepared using Nextera XT DNA Library Preparation Kit (Illumina Inc.) on cDNAs with an initial concentration > 180 pg/µL that were then diluted to 100 pg/µL. In each single-cell library preparation, a total of 125 pg cDNA was tagmented at 55°C for 20 minutes. Libraries were pooled and purified using AMPure bead-based magnetic separation before a final quality control using Qubit dsDNA HS Assay Kit (Life Technologies) and Agilent High Sensitivity DNA Kit. We required the majority of cDNA fragments to be between 375 and 425 bp to qualify for sequencing. Single-cell libraries were subjected to paired-end 75 bp RNA-seq on a HiSeq 2500 (Illumina Inc.). Ninety-six scRNA-Seq libraries were sequenced per HiSeq 2500 cell (∼300 million bp/cell).
scRNA-seq Analysis
The Fluidigm C1 RNA-seq data were aligned to the human transcriptome (hg19) and quantified using the software RSEM (80). Cells with a minimum of 100,000 aligned reads were retained for further analysis. To assign cell types to these profiles, the software cellHarmony (42) was used to classify individual cells according to 18 cell populations from the bone marrow single-cell HCA (43). Differential expression analyses were performed in the software AltAnalyze (81), using an empirical Bayes-moderated t test, following FDR correction where indicated. Clustering for the heat map of noncycling young versus aged differentially expressed genes was performed using HOPACH in AltAnalyze. To determine whether any young HSC had an aged gene signature profile, expression centroids from the three clusters identified by HOPACH (n = 1 young cluster and n = 2 aged clusters) were used to reclassify young and aged HSC.
Enhancer Analysis
To identify young poised enhancers [(H3K4me1>H3K4me3), H3K27me3+, H3K27ac−, and not within the promoter region], peak calling of young pooled H3K4me1 was performed as described above. The enrichment of H3K4me1 and H3K4me3 at these regions was calculated using deeptools2 multiBigwigSummary (62) with bigwig files that were read normalized and contained the fold enrichment of the Pooled IP/Pooled Input. Peaks where H3K4me1 enrichment was not greater than twice that of H3K4me3 were discarded. Genomic annotation of the remaining peaks was performed using the R-package Genomation (v1.2.2; ref. 66), and peaks within 3 kb of a TSS were removed from further analysis. Peaks overlapping (by at least 1 base pair) H3K27me3, but not H3K27ac peaks, were found with bedtools intersect. Active enhancers were defined as H3K27ac+, (H3K4me1 > H3K4me3), H3K27me3−, and not within the promoter region. To identify young active enhancers, peak calling of H3K27ac was performed as described above. Genomic annotation of the resulting peak files was then performed, and peaks within 3 kb of a TSS were discarded. Bedtools intersect was used to find sites that overlapped with regions where H3K4me1 fold enrichment/H3K4me3 fold enrichment > 2, but did not overlap with H3K27me3 peaks. Active and poised enhancers were annotated to the nearest TSS using the R package ChIP-enrich with hg19 RefSeq (64).
Bivalent Promoter Analysis
Peak calling of pooled young HSCe H3K4me3 and H3K27me3 was performed as described above. Promoter peaks (± 3,000 bp of TSS) were identified using the R-package Genomation (v1.2.2; ref. 66). The bedtools intersect function was used to identify regions that were bound by both H3K4me3 and H3K27me3 in young HSCe (overlap >25%). Peaks overlapping bivalent promoters (n = 4,662) were annotated to the nearest TSS using hg19 RefSeq and using ChIP-enrich to identify the total number of genes with bivalent promoters (n = 3,604; ref. 66).
Western Blot Analysis
Western blotting was performed on 5,000 HSCe per donor (n = 3 per age group) and 10,000 HSCe per donor with 1–2 technical replicates for an additional 5 young and 4 aged donors. HSCe from a 53-year-old donor were used to determine the linear range of the assay. Protein was extracted using a modified trichloroacetic acid (TCA; MP Biomedicals #196057) precipitation protocol (82). Briefly, HSCe were sorted into 10% TCA and precipitated at 4°C. Extracts were spun at 13,000 rpm for 10 minutes at 4°C. Supernatant was discarded and samples were washed twice with acetone. Pellets were resuspended in TCA lysis buffer (9 mol/L Urea, 2% Triton X-100, 1% DTT) and lithium dodecyl sulphate (Life Technologies NuPage, #NP0007). For samples with 5,000 HSCe, 10 mmol/L Tris-HCl pH 8.0 was also added. Samples were denatured at 70°C for 10 minutes. Protein from samples with 10,000 HSCe were loaded into a 4% to 12% Bis-Tris gel (Thermo Fisher Scientific, #NP0322BOX), 20 µL/well and then transferred to a polyvinylidene difluoride (PVDF) membrane. Blots were imaged on a Bio-Rad ChemiDoc. Densitometry was performed using the Bio-Rad Image Lab software. For samples with 5,000 HSCe, and the standard curve, protein was loaded into a 15% SDS-PAGE gel, 13 µL/well, and then transferred to a Immobilon FL PVDF membrane (Millipore, #IPFL10100). Blot was imaged on a LI-COR Odyssey CLx. Densitometry was performed using the LI-COR Image Studio software. For all samples, intensity of H3 was normalized to that of a loading control, GAPDH. GAPDH-normalized values were then normalized to the average young intensity. P value was calculated using the parametric unpaired t test, which was performed using Prism software. The following antibodies were used: rabbit αGAPDH (Santa Cruz Biotechnology, #SC25778, lot K0615), mouse αH3 (Abcam, #ab10799), αMouse IgG, HRP-linked (Cell Signaling Technology, #7076S, lot 32), and αRabbit IgG, HRP-linked (Cell Signaling Technology, #7074S, lot 26), αRabbit IR-Dye 800 CW (LI-COR, #92632213 lot C61012-02), and αMouse IR-Dye 680 RD (LI-COR, #925-68070, lot C70613-11).
Transcription Factor Binding Motif Analysis
The Homer findMotifs function was used to detect significant (q-value < 0.05) enrichment of transcription factor binding motifs in clusters altered with age (83). For comparison of HSCe ChIP-seq to transcription factor binding sites in CD34+ cells, bedtools intersect was used to find the overlap between published transcription factor ChIP-seq peaks in CD34+ cells and the peaks within the enhancer clusters J–L (36).
Peak Visualization
Peaks were visualized using the UCSC genome browser. The macs2 SPMR and bdgcmp functions were used to generate fold enrichment ChIP-seq tracks that are normalized by read count and to the IP's corresponding Input. Destranded RNA-seq tracks of the individual and pooled replicates were created using the STAR aligner parameters: –outWigType bedGraph –outWigNorm RPM –outWigStrand Unstranded.
Gene Ontology Analysis
Functional annotation of ChIP-seq peaks, age-associated clusters, and regions predisposing to AML was performed using the R package ChIP-enrich (64). Regions were annotated to RefSeq hg19 Gene Ontology Biological Processes data sets of less than 500 genes, using the method = chipenrich function for H3K4me1, H3K4me3, H3K27ac, and clusters. H3K27me3 was annotated using method = broadenrich. All ChIP-seq peaks and age-associated clusters were annotated to the nearest TSS, and DMRs predisposing to AML were annotated to the nearest gene. Results with FDR < 0.05 were considered highly significant.
Heat Maps and Density Plots of ChIP-seq
Heat maps of differential ChIP-seq peaks, enhancers, and bivalent promoters were created using deepTools2 (v2.5.0.1; ref. 62). Individual and pooled ChIP-seq IP replicates were normalized to their corresponding inputs using the deepTools2 bamCompare function and the options –ratio log2 –outFileFormat bigwig –scaleFactorsMethod readCount (62). The resulting bigwig files were then used in conjunction with deepTools2 computeMatrix and plotHeatmap to calculate and plot enrichment at regions of interest. Heat maps are peak centered on the differential regions.
Box Plots of ChIP-seq Data
For box plots of histone modifications lost with age, pooled ChIP-seq IP replicates were log2 normalized as described above. Enrichment of each IP was then calculated for H3K4me1, H3K4me3, H3K27me3, and H3K27ac peaks lost with age using the deepTools2 computeMatrix function. For box plots of genome-wide H3 enrichment, individual replicates (n = 2 young and n = 2 aged) were log2 normalized to their corresponding inputs as described above. The deepTools2 multiBigwigSummary bins function was then used to calculate H3 enrichment over 150-bp bins covering the whole genome. Box plots were then plotted using base R.
Targeted Genomic Sequencing
Genomic DNA was used to prepare libraries for hybrid capture as per the manufacturer's protocol (Agilent). Libraries were quantified and pooled up to 24 samples per reaction in equimolar amounts totaling 500 mg of DNA. Agilent Custom SureSelect In Solution Hybrid Capture RNA baits were used to hybridize libraries, targeting 443 kbp of exonic DNA with 16,890 probes as described in Lindsley and colleagues (84). Each capture reaction was washed, amplified, and sequenced on two lanes of an Illumina HiSeq 2000 100-bp paired-end run. Subsequent analysis of the target region was restricted to the regions corresponding to the following genes: ANKRD26, CEBPA, ETV6, KIT, PRPF40B, SF3A1, TERC, ASXL1, CREBBP, EZH2, KRAS, PRPF8, SF3B1, TERT, ATRX, CSF1R, FLT3, LUC7L2, PTEN, SH2B3, TET2, BCOR, CSF3R, GATA2, MPL, PTPN11, SMC1A, TP53, BCORL1, CSNK1A1, GNAS, NF1, RAD21, SMC3, U2AF1, BRAF, CTCF, GNB1, NPM1, RAD51C, SRSF2, U2AF2, BRCC3, CUX1, IDH1, NRAS, RUNX1, STAG1, WT1, CALR, DDX41, IDH2, PHF6, SBDS, STAG2, ZRSR2, CBL, DNMT3A, JAK2, PIGA, SETBP1, STAT3, CBLB, EP300, JAK3, PPM1D, SF1, and STAT5B.
Variant Calling
Fastq files were aligned to the hg19 version of the human genome with The Burrows Wheeler Aligner (BWA v0.7.12) MEM module for paired-end reads (85). Duplicate reads were flagged and removed using Picard tools (V1.91). GATK v3.2 was used for base recalibration prior to variant calling, and also for local realignments for insertions/deletions (indels) using the reference variant databases (86). Somatic variants were called using LoFreq v2.1.1 for all variants at ≥ 1% variant frequency (87). Additional variant filtering after variant calling was used with the following parameters: variant frequency < 0.05, read depth at variant site <20, GQ and/or QUAL scores <30, IndelRepeatFilter > 8. Variants with excessive strand bias and indels with VAF < 10% adjacent to homopolymer repeats were excluded by manual curation. Filtered called variants were first annotated using ANNOVAR (88). Variants predicted to alter splicing were assessed as described in Jian and colleagues (89). Variants located outside protein-coding regions or splice sites of the genes listed above, and synonymous variants that were not predicted to alter splicing were filtered out. To remove common polymorphisms, variants with population frequencies of ≥1% in either 1,000 genomes (90) or the Exome Aggregation Consortium (ExAC v.3.1) were similarly excluded unless they were also listed as confirmed somatic mutations in the Catalogue of Somatic Mutations in Cancer (COSMIC; ref. 91). Remaining variants were manually reviewed and those considered likely somatic-coding variants were included in further analyses. Variant calling and interpretation were performed blinded to sample identifiers and associated phenotype information.
Lentiviral Production and Transduction of CD34+ Cells
The 293FT cells were maintained according to supplier instructions (Thermo Fisher Scientific). A custom pLKO.1 containing turboGFP vector was used to express shRNA targeting LMNA (5′-GATGATCCCTTGCTGACTTAC-3′; Millipore Sigma). Lentiviruses were produced by cotransfection with packaging plasmids psPAX2 and pMD2.G using polyethylenimine transfection reagent (Polysciences). Lentivirus-containing supernatant was collected 48 to 72 hours post-transfection, filtered through a 0.45-µm syringe filter, and concentrated using PEG-it virus precipitation solution (System Biosciences, #LV825A-1). Primary human CD34+ cells were freshly isolated from mobilized peripheral blood obtained from the University of Michigan Comprehensive Cancer Center Tissue Procurement Service or purchased from Fred Hutchinson Hematopoietic Cell Procurement Services. CD34+ cells were isolated using magnetic bead purification as described for bone marrow (Miltenyi Biotec). Isolated CD34+ cells were cultured in 20% FBS containing IMDM or SFEM II (StemCell Technologies, #09605) and prestimulated with recombinant human SCF (100 ng/mL), FLT3-L (10 ng/mL), IL6 (20 ng/mL), and TPO (100 ng/mL; PeproTech). Lentiviral transduction of CD34+ cells was performed in the presence of 8 µg/mL polybrene (Millipore Sigma, TR-1003-G). Four days post-transduction, cells were sorted for CD34 and GFP double-positive cells on a FACSAria IIu (BD Biosciences). Empty vector was used as control comparison.
Colony-Forming Unit Assay
Sorted CD34+ were seeded in methylcellulose, MethoCult H4435 (StemCell Technologies, #04435), in duplicate onto a 6-well SmartDish (StemCell Technologies, #27302) at varying densities of 500 to 1,500 cells. Colonies were scored on a STEMvision after 14 days of incubation (StemCell Technologies; n = 7 and n = 5 for shLMNA and sgKLF6, respectively).
sgRNA Synthesis
Protospacer sequences were identified using the CRISPRdesign algorithm (http://crispr.mit.edu) or CRISPRScan (www.crisprscan.org; ref. 92). DNA templates for sgRNAs contain a T7 promoter, the protospacer sequence, and the sgRNA scaffold sequence (93). They were produced by PCR using custom forward primers and a reverse primer that amplifies the sgRNA scaffold of the plasmid pKLV-U6gRNA-PGKpuro2ABFP (Addgene #62348). The PCR products were purified and in vitro transcribed with the HiScribe T7 High Yield RNA Synthesis Kit (NEB #E2040S) following manufacturer instructions. In vitro transcribed sgRNA products were purified using RNA Clean & Concentrator Kit (Zymo Research #R1015). See also Supplementary Table S11.
CRISPR/Cas9 Targeting of CD34+ Cells
Peripheral blood CD34+ cells were thawed and cultured in stem cell–promoting media [IMDM + 20% Human Serum AB + hTPO (100 ng/mL), hSCF (100 ng/mL), hIL6 (20 ng/mL), hFLT-3 (10 ng/mL), and nonessential amino acids] under low oxygen conditions to maximize maintenance and expansion of LT-HSCs (94–96). After 72 hours in culture, cells were then transfected with equimolar amounts of Cas9 protein (PNA Bio; 1 µg/µL in PBS) and sgRNAs using the Neon Transfection System (Thermo Fisher Scientific). Electroporation conditions used were 1,600 volts, 10 milliseconds, and 3 pulses. A nontargeting sgRNA was used as control. After transfection, cells were cultured in stem cell–promoting media for 72 hours under low oxygen conditions (94–96) before sorting for CD34+ cells. After an additional 72 hours in culture, cells were used for RNA extraction, colony-forming unit assay, or plated in myeloid or erythroid differentiation media (see conditions below). Viability was assessed by Trypan blue exclusion and was generally above 55% for 48 to 72 hours post-transfection. Cutting was validated using the T7E1 assay as follows: PCR amplicon spanning the cleavage site was gradually hybridized in a thermal cycler. Hybridized fragments were then digested with 1.25 U of T7 endonuclease I (NEB) for up to 30 minutes at 37°C. Cleavage at heteroduplex mismatch sites were detected by agarose gel electrophoresis. Knockout of KLF6 was also verified using flow cytometry with an anti-KLF6 antibody (Millipore, #MABN119) and goat-anti-mouse IgG-AF594 (Thermo Fisher Scientific, #A-11032) or goat-anti-mouse IgG-AF488 (Thermo Fisher Scientific, #A-11059).
RNA-seq of LMNA Knockdown
RNA from FACS-isolated transduced CD34+ cells (n = 8 donors) was extracted using the Qiagen RNeasy Plus Micro Kit according to manufacturer's instructions (Qiagen, #74034), or TRIzol. Stranded libraries were prepared by the University of Miami Sequencing Core using the Illuminia TruSeq Stranded Total RNA Kit (Illumina, #20020596). Libraries were sequenced on the NextSeq-500 with 75 bp paired-end sequencing. Data were aligned and processed as described above. For GSEA, the Wald statistic ranked list was used with the top (ranked by log2 fold change) 500 genes upregulated and downregulated with HSCe aging as gene sets.
RNA-seq of KLF6 Knockdown
RNA from transfected CD34+ cells (n = 3 biological replicates and 1–3 technical replicates for a total of n = 4 sgCTRL and n = 4 sgKLF6) was extracted using the Qiagen Allprep Micro Kit according to the manufacturer's instructions (Qiagen, #80204). Stranded libraries were prepared by the University of Miami Sequencing Core using the Illumina TruSeq Stranded Total RNA kit (Illumina, #20020596). Libraries were sequenced on the HiSeq-3000 with 75 bp paired-end sequencing. Data were aligned and processed as described above. For GSEA, the Wald statistic ranked list was used with the top (ranked by log2 fold change) 500 genes upregulated and downregulated with HSCe aging as gene sets as well as the c2.all.v6.2.symbols gene set using the weighted enrichment score.
Myeloid and Erythroid Liquid Culture
Cells were plated and cultured for 7 days under myeloid-promoting conditions: SCF (100 ng/mL), FLT3 ligands (10 ng/mL), IL3 (20 ng/mL), IL6 (20 ng/mL), GM-CSF (20 ng/mL), and G-CSF (20 ng/mL) and erythroid-promoting conditions: EPO (6 IU/mL) and SCF (100 ng/mL). At day 7 of myeloid and erythroid expansion, cells were stained for CD34 (BD Pharmingen, #304441), myeloid CD11b (BD Pharmingen, #301324), and erythroid CD71 (Pharmingen, #563769) and CD235a (BD Pharmingen, #559943) markers, respectively, and also anti-KLF6 for cells transfected with sgKLF6 (Millipore, #MABN119).
Statistical Analysis
Significance details can be found within the figure legends. For genome-wide sequencing assays, we corrected for multiple testing, and used q < 0.05 as a cutoff to determine significance, for all other statistical tests, we used a P value cutoff of P < 0.05 to define significance. In the figures, *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. For the Western blot analysis, a two-sided parametric unpaired t test was used to calculate P values with Prism software (v7.0).
Data and Software Availability
The ChIP-seq, ERRBS, and bulk RNA-seq data have been deposited in the NCBI Gene Expression Omnibus under ID code GSE104408. The scRNA-seq data has been deposited in the NCBI Gene Expression Omnibus under ID code GSE104379. Code used for analyses is available on Github (https://github.com/exa564/HSCe_Aging).
Disclosure of Potential Conflicts of Interest
R.C. Lindsley is a consultant at Takeda Pharmaceuticals, reports receiving commercial research grants from MedImmune and Jazz Pharmaceuticals, and is a consultant/advisory board member for Jazz Pharmaceuticals. R. Bejar reports receiving a commercial research grant from Takeda and is a consultant/advisory board member for Celgene, Daiichi-Sankyo, Astex, and NeoGenomics. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: E.R. Adelman, N. Salomonis, H.L. Grimes, M.E. Figueroa
Development of methodology: E.R. Adelman, A. Roisman, R.C. Lindsley, N. Salomonis, M.E. Figueroa
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): E.R. Adelman, H.-T. Huang, A. Roisman, A. Olsson, R. Bejar, M.E. Figueroa
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E.R. Adelman, H.-T. Huang, A. Roisman, A. Colaprico, T. Qin, R. Bejar, N. Salomonis, H.L. Grimes, M.E. Figueroa
Writing, review, and/or revision of the manuscript: E.R. Adelman, H.-T. Huang, R. Bejar, N. Salomonis, H.L. Grimes, M.E. Figueroa
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): E.R. Adelman
Study supervision: M.E. Figueroa
Acknowledgments
E.R. Adelman was supported by the T32-AG11431 Biology of Aging grant to the University of Michigan and an award from the University of Michigan Department of Pathology Pilot Training Program in Translational Research. This work was supported by a Quest for Cures grant (0858-15 and R0858-18) as well as a Scholar Award (1357-19) from The Leukemia and Lymphoma Society (to M.E. Figueroa) and R01 HL122661 (to H.L. Grimes). M.E. Figueroa, N. Salomonis, and H.L. Grimes were also supported by R01 CA196658. The authors would like to acknowledge the support of the tissue procurement services of the University of Michigan Comprehensive Cancer Center Tissue Core, the University of Michigan Medical School (UMMS) DNA Sequencing Core, the UMMS Flow Cytometry Core, the University of Miami Sylvester Comprehensive Cancer Center (UM SCCC) Onco-Genomics Shared Resource, and the UM SCCC Flow Cytometry Shared Resource, as well as the Cincinnati Children's Hospital (CCH) Cell Processing Core, which is supported by the CCH Research 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.
References
Supplementary data
Supplementary Figures
Donor information
Differential ChIP-seq peaks
Differentially methylated regions
Gene expression
Gene Set Enrichment Analysis
Overlap with aging signatures
sc-RNAseq
Regions predisposing to AML
Oligonucleotides