Abstract

   Nutrition and metabolism are known to influence chromatin biology and
   epigenetics through post-translational modifications, yet how this
   interaction influences genomic architecture and connects to gene
   expression is unknown. Here we consider, as a model, the
   metabolically-driven dynamics of H3K4me3, a histone methylation mark
   that is known to encode information about active transcription, cell
   identity, and tumor suppression. We analyze the genome-wide changes in
   H3K4me3 and gene expression in response to alterations in methionine
   availability in both normal mouse physiology and human cancer cells.
   Surprisingly, we find that the location of H3K4me3 peaks is largely
   preserved under methionine restriction, while the response of H3K4me3
   peak width encodes almost all aspects of H3K4me3 biology including
   changes in expression levels, and the presence of cell identity and
   cancer-associated genes. These findings may reveal general principles
   for how nutrient availability modulates specific aspects of chromatin
   dynamics to mediate biological function.
     __________________________________________________________________

   Methionine availability is known to affect the global levels of histone
   methylation. Here the authors investigate the metabolically driven
   dynamics of H3K4me3 and find that methionine availability influences
   peak width, which is linked to changes in expression of genes
   associated with cell fate and cancer.

Introduction

   Genes interact with environmental factors such as nutrition to shape
   the epigenome that together influences gene activity and organismal
   physiology. Metabolism is also shaped by genes and environment and has
   a substantial contribution to epigenetics^[32]1–[33]4. This nexus is
   essential in numerous biological contexts, including maintaining
   different stages of pluripotency^[34]5–[35]8, mediating an immune
   response^[36]9,[37]10, promoting or suppressing cancer
   progression^[38]11–[39]16, and transducing information about metabolic
   health and longevity from parent to offspring^[40]17–[41]19. The
   molecular foundation of this interaction is in large part determined by
   the modifications on chromatin.

   Chromatin is affected by metabolism through changes in the
   concentrations of metabolites that serve as substrates and cofactors
   for post-translational modifications. These concentrations are dynamic
   and are mediated by changes in metabolic pathway activity or flux that
   arise from transcriptional programs and nutrient availability. For
   example, histone methylation requires S-adenosylmethionine (SAM) as the
   universal methyl donor. SAM is derived from methionine^[42]20 and its
   concentration can fluctuate in physiological conditions around values
   that can limit the activity of histone methyltransferases^[43]21.

   In plasma, methionine is in some reports the most dynamic of the 20
   amino acids and the variation can to a large extent be explained by
   diet^[44]22. Recently work from us and others has shown that dietary
   modulation of methionine concentrations that approach the lower end of
   what can be observed in humans leads to bulk changes in the levels of
   histone methylation^[45]22,[46]23. Other studies have reported similar
   findings in that changes to SAM levels or to the levels of
   alpha-ketoglutarate that modify the activity of demethylase enzymes
   induce global changes in the levels of histone
   modifications^[47]5,[48]12,[49]24–[50]32. When these modifications are
   known to mark key aspects of chromatin status, global changes could
   have broad consequences to epigenomic programs. How these bulk changes
   to the levels of post-translational modifications on chromatin alter
   the genomic architecture of histone marks and relate to gene expression
   is, however, largely unknown.

   One attractive model to investigate this interaction at the genome
   scale is the tri-methylation of histone H3 on lysine 4 (H3K4me3). The
   global (i.e., bulk) levels of this mark are dynamically and reversibly
   responsive to the levels of methionine^[51]22. In addition, there are
   numerous lines of evidence indicating that the structural features of
   H3K4me3, such as the width or breadth of the peak as deposited over a
   genic region, encode information such as gene activity, and gene
   function such as the presence of a developmental program, cell type
   identity, or a tumor suppressor^[52]33–[53]37. Thus, changes in H3K4me3
   may be relevant to developmental transitions and tumor suppression. How
   metabolic dynamics that occur due to differences in nutritional status
   or metabolic pathway activity might affect these programs and gene
   activity related to H3K4me3 is largely unknown.

   We have shown previously that methionine availability modulates bulk
   levels of H3K4me3 by modifying SAM concentrations^[54]22. In this
   present study, we question whether changes in methionine availability
   that are known to affect global levels of H3K4me3 affect specific
   aspects of the genomic architecture and gene expression regulation. We
   consider a mouse model of dietary methionine restriction (MR) and
   focused our analysis on liver. In this organ, this diet results in
   changes to bulk levels of H3K4me3. Similar changes occur in cultured
   human cancer cells (HCT116) subjected to acute MR in culture media,
   that together provide a complementary set of two species, environmental
   conditions, biological statuses (health and cancer), models (in vitro
   and in vivo) and two tissues. We study genome-wide H3K4me3 dynamics
   using a quantitative ChIP-seq analysis that considers peak geometry and
   characterize the connection to gene expression dynamics. We find that
   height and area of the peaks are overall reduced, which account for
   most of the global changes. Strikingly, however, while the most
   conserved feature of H3K4me3 dynamics is the peak width, changes in
   peak width but not other features of peak geometry reflect important
   cellular processes previously linked to H3K4me3, including cell
   identity-related gene expression programs and the dynamics of gene
   expression.

Results

MR reduces H3K4me3 but maintains its genomic distribution

   To begin to study the impact of methionine availability on the genomic
   architecture of H3K4me3, we applied chromatin immunoprecipitation with
   sequencing (ChIP-seq) to map genomic locations enriched with H3K4me3 in
   HCT116 cells cultured under high (100 μM) and low, MR (3 μM) methionine
   availability which is known to lead to several-fold global changes in
   the bulk levels of H3K4me3^[55]22. Aligning the reads to a reference
   genome followed by peak calling identified a set of H3K4me3 peaks,
   i.e., genomic regions significantly enriched (MACS2 enrichment Q-value
   <1e-5) with ChIP-seq reads in comparison to the control (Fig. [56]1a).
   Comparing total peak number (coefficient of variation (CV) = 0.01,
   Fig. [57]1b), genomic location (Jaccard index = 0.88, Fig. [58]1c), and
   the set of genes marked by peaks (Jaccard index = 0.92, Supplementary
   Fig. [59]1a) between high and low methionine conditions showed that the
   distribution of H3K4me3 genomic locations was highly conserved in
   response to MR (Fig. [60]1b, c) while an overall reduction in H3K4me3
   in peak intensity also observed (Supplementary Fig. [61]1b, c). Genes
   with gained or lost H3K4me3 peaks in response to MR tended to be
   observed in smaller peaks (Wilcoxon rank-sum P-value < 1e-4 for five
   out of the six comparisons, Supplementary Fig. [62]1d, e) and were not
   significantly enriched with specific biological functions compared to
   randomly chosen gene sets (median(−log[10](enrichment Q-value)) = 4.25
   for gained peaks compared to 7.76, 8.05, and 8.98 for three random gene
   sets of identical size, Supplementary Fig. [63]1f, and
   median(−log[10](enrichment Q-value)) = 3.00 for lost peaks compared to
   7.12, 6.05, and 5.08 for three random gene sets with identical size,
   Supplementary Fig. [64]1g) and thus likely attributable to technical
   noise. We also analyzed the composition of genomic elements covered by
   the peaks and found, as has been reported with H3K4me3^[65]38, that
   most peaks (11,210 peaks, 81%) covered promoter regions with a smaller
   subset of peaks (2599 peaks, 19%) that appeared on non-promoter regions
   such as intergenic regions and introns (Supplementary Fig. [66]1h). To
   further investigate these changes, we next developed several
   quantitative descriptors for individual peaks (Fig. [67]1d). We used
   height, area, and width to evaluate the size of the peaks and compared
   these metrics in high and low methionine conditions. Genomic regions in
   the high and low methionine conditions were merged (Fig. [68]1a) to
   include less abundant peaks called only in one condition in the
   analysis. All three peak size descriptors showed high correlations
   (Spearman’s rank correlation coefficient >0.95 and random permutation
   test P-value < 1e-323) between high and low methionine conditions
   (Fig. [69]1e), implying that the overall H3K4me3 landscape (i.e.,
   relative size of each individual peak) is robust to methionine
   availability (linear regression coefficient = 0.87 for width, 0.85 for
   height, and 0.70 for area).

Fig. 1.

   [70]Fig. 1
   [71]Open in a new tab

   MR reduces H3K4me3 but maintains its genomic distribution. a MR
   experiment and ChIP-seq data analysis pipeline in human cancer cells
   HCT116. b Number of H3K4me3 peaks called under high (100 μM) and low
   (3 μM) methionine conditions in human cancer cells. c Total length of
   shared and unique peak regions for high and low methionine conditions
   in human cancer cells. d Definition of peak height, area, and width. e
   Density scatter plots comparing H3K4me3 peak heights, areas, and widths
   between high and low methionine conditions in human cancer cells.
   Colors represent for dot density. Dashed lines show identical x and y
   coordinates. Spearman’s rank correlation coefficients between the high
   and low methionine conditions and linear regression coefficients are
   shown in the table at bottom. f MR experiments scheme in vivo. Mouse
   image used with permission from Microsoft. g Density scatter plots
   comparing H3K4me3 peak heights, areas, and widths in mouse liver
   between high (0.84% w/w) and low (0.12% w/w) methionine diets

   To determine whether this conservation in H3K4me3 dynamics under MR
   extends to normal physiology and under a longer-term alteration in
   nutrient availability that has beneficial health effects, we
   investigated the effects of MR on H3K4me3 in vivo. Under these
   conditions, MR has been shown to alter methionine metabolism, improve
   metabolic physiology in liver, extend life-span, and reduce global
   levels of H3K4me3^[72]22,[73]39–[74]41. Thus, this mouse model of liver
   physiology in combination with the human cancer cell model provides a
   breadth of model systems covering both short-term and long-term
   nurtrient alterations, both in vitro and in vivo systems, and both
   pathological and healthy contexts. We focused on liver in profiling the
   epigenomics and transcriptomics because it is the metabolic organ which
   is most responsive to metabolic reprogramming and there are
   liver-related phenotypes associated with MR related to metabolic
   health^[75]40. Livers were obtained from C57BL6 adult mice on a diet
   with either high methionine (0.84% w/w) or low, MR methionine (0.12%
   w/w) for 12 weeks, and ChIP-seq of H3K4me3 was conducted (Fig. [76]1f).
   Quantitation of bulk levels of H3K4me3 under these two conditions
   confirmed as was previously published that this level of dietary MR was
   sufficient to reduce the levels of H3K4me3^[77]22. This physiological
   MR also revealed that the overall distribution of H3K4me3 was conserved
   between the different nutrient conditions as corroborated by conserved
   peak numbers (CV = 0.06, Supplementary Fig. [78]2a), location (Jaccard
   index = 0.82, Supplementary Fig. [79]2b), and marked genes (Jaccard
   index = 0.90, Supplementary Fig. [80]2c). High correlations in the
   values of peak height, area, and width between high methionine and low
   methionine conditions were also observed with peak width as the most
   conserved (Fig. [81]1g, linear regression coefficient = 0.99 for width
   compared to 0.85 for height and 0.80 for area). Consistent with what
   was found in human cancer cells, H3K4me3 peaks gained or lost in
   response to MR in mouse liver also tended to be significantly smaller
   peaks than those conserved in both high and low methionine conditions
   (Wilcoxon rank-sum P-value < 1e-3 for five out of the six comparisons,
   Supplementary Fig. [82]2d, e) and were not enriched with specific
   biological function compared to randomly chosen genes containing
   H3K4me3 (median(−log[10](enrichment Q-value)) = 3.09 for gained peaks
   compared to 4.03, 4.76, and 3.39 for three random gene sets with
   identical size, Supplementary Fig. [83]2f, and
   median(−log[10](enrichment Q-value)) = 14.62 for lost peaks compared to
   16.58, 15.75, and 13.27 for three random gene sets with identical size,
   Supplementary Fig. [84]2g). Reduction in the overall H3K4me3 signal was
   also observed (Supplementary Fig. [85]2h, i) as was a conserved
   distribution of genomic elements (84% promoter peaks and 16%
   non-promoter peaks, Supplementary Fig. [86]2j).

   To assess the robustness of the quantitative peak descriptors to
   variation in ChIP-seq data analysis methodologies, we repeated the
   calculations using several different peak-calling algorithms. Although
   the total number of peaks (CV = 0.12 for human cancer cells,
   Supplementary Fig. [87]3a, and CV = 0.24 for mouse liver, Supplementary
   Fig. [88]3b) and total length of genomic regions covered by peaks
   (Jaccard index = 0.60 for human cancer cells, Supplementary
   Fig. [89]3c, and Jaccard index = 0.43 for mouse liver, Supplementary
   Fig. [90]3d) exhibited some variation among methods, reads mapped to
   peaks (CV = 0.0073 for human cancer cells, Supplementary Fig. [91]3e,
   and CV = 0.0078 for mouse liver, Supplementary Fig. [92]3f), genes
   associated with peaks (Jaccard index = 0.78 for human cancer cells,
   Supplementary Fig. [93]3g, and Jaccard index = 0.66 for mouse liver,
   Supplementary Fig. [94]3h), and dynamics of peak height, area, and
   width (average Spearman’s rank correlation coefficient = 0.92 for human
   cancer cells and 0.86 for mouse liver, random permutation test
   P-values < 1e-323 for all comparisons, Supplementary Fig. [95]4) were
   concordant, implying that quantitation of peak dynamics was robust to
   the ChIP-seq data analysis method used. Taken together, these results
   indicate that changes in metabolism mediated by methionine
   availability, despite altering global levels of histone methylation, do
   not induce a complete redistribution of H3K4me3 marks on the genome.

H3K4me3 width dynamics encode biological information

   Although the genomic positioning of H3K4me3 was conserved in response
   to changes in methionine availability, we observed that peak height,
   area, and width were affected to different extents. To further
   investigate these differences, we computed the Spearman’s rank
   correlation coefficients between fold changes in peak height, area, and
   width for both HCT116 cells and mouse liver and found that fold changes
   in the width and in the other two parameters were less correlated
   (Spearman’s rank correlation coefficient < 0.4) in both models
   (Fig. [96]2a, b) relative to other comparisons. The lower correlation
   between changes in peak width and changes in the other two peak
   features suggested that a change in peak width might encode a different
   dimension of information. To test this hypothesis, we conducted pathway
   enrichment analysis with the Gene Set Enrichment Algorithm
   (GSEA)^[97]42 on the peaks with sensitive (larger reduction under MR)
   and robust (smaller reduction under MR) height, area, and width
   (Fig. [98]2c) using the Molecular Signatures Database (MSigDB)^[99]43.
   For each category of peak dynamics, the number of significantly
   enriched (GSEA FDR Q-value < 1e-5) pathways was used to quantify
   association of this category with biological functions. In cultured
   cancer cells, the change in peak width exhibited the strongest signal
   for enrichment of specific biological processes as measured by the
   number of significantly enriched pathways (298 pathways significantly
   enriched in sensitive width compared to 131 in robust width and less
   than 60 in all other categories, Fig. [100]2d). In mouse liver, the
   consistency of the peak width exhibited the strongest signal (48
   pathways significantly enriched in robust width compared to 14 in
   robust area and 0 in all other categories, Fig. [101]2e). This
   surprising finding appears to indicate that biological information that
   occurs in response to changes in methionine availability is encoded in
   H3K4me3 peak width.

Fig. 2.

   [102]Fig. 2
   [103]Open in a new tab

   H3K4me3 width dynamics encode biological information. a Heat map
   showing Spearman’s rank correlation coefficients among fold changes of
   H3K4me3 peak heights, areas, and widths under MR in human cancer cells.
   The upper right part shows density scatter plots comparing the
   corresponding log2(fold change) values. Colors of dots in the scatter
   plots indicate density of dots. b Same as in a but for mouse liver. c
   Distribution of fold changes (defined by the values in low methionine
   condition divided by the corresponding values in high methionine
   condition) in peak width. The arrow at the bottom denotes the change
   from sensitive (more reduction in peak width) to robust (less reduction
   in peak width). d Number of pathways significantly enriched (GSEA FDR
   Q-value < 1e-5) in peaks with different dynamics under MR in human
   cancer cells. e Same as in d but for mouse liver. f Examples of
   annotated pathways in each category. g Annotation of 298 MSigDB
   pathways enriched in H3K4me3 peaks with sensitive width in human cancer
   cells. h Representative cancer-related pathways enriched in H3K4me3
   peaks with sensitive width in human cancer cells. i Annotation of 48
   MSigDB pathways enriched in H3K4me3 peaks with robust width in mouse
   liver. j Representative liver-specific pathways enriched in H3K4me3
   peaks with robust width in mouse liver

   To further quantify the information contained in peak width, we
   annotated the significantly enriched pathways in each case by
   classifying them into categorical subsets involving fundamental
   cellular function, response to stimulus, cell type (cancer,
   proliferation-related pathways for human cancer cells and metabolism,
   multicellular organ-related pathways for mouse liver), transcriptional
   regulation, differentiation and development, and cancer (Fig. [104]2f).
   We next computed the fraction of each category in the pathways enriched
   in sensitive width in cancer cells and those enriched in robust width
   in liver. Unexpectedly, these two sets were enriched with cell
   type-specific biological functions. Proliferation and cancer-related
   pathways were enriched in peaks with sensitive width in cancer cells
   (154 out of the 298 enriched pathways, Fig. [105]2g, h), while
   metabolism and multicellular organ-related pathways were enriched in
   peaks with robust width in mouse liver (35 out of the 48 enriched
   pathways, Fig. [106]2i, j). Together these findings indicate that the
   peak width dynamics in response to MR encodes information about
   biological function.

H3K4me3 width dynamics encode cell type-specific TF binding

   To further explore this relationship, we probed additional aspects of
   genomic architecture. Using a motif analysis algorithm, HOMER^[107]44,
   we searched for transcription factor (TF) binding motifs enriched in
   the subsets of the different geometrical features of H3K4me3 peaks that
   change or remain consistent during MR. Peak subsets with the 500 most
   sensitive peaks or 500 most robust peaks were considered in this
   analysis (Fig. [108]3a). The complete set of H3K4me3 peaks was used as
   the background model (Fig. [109]3b). In cultured cancer cells, TF
   binding motifs were only enriched (one-sided Fisher’s exact test
   Q-value < 0.05) in peaks with sensitive width (78 motifs enriched in
   sensitive width compared to 1 for robust area, 2 for robust width, and
   0 for all other peak sets, Fig. [110]3c, Supplementary Fig. [111]5a),
   while in mouse liver TF binding motifs were only enriched in peaks with
   robust width (249 motifs enriched in robust width compared to 0 for all
   other peak sets, Fig. [112]3d, Supplementary Fig. [113]5b). These
   findings were also found to be robust to the size of peak subsets
   chosen (Supplementary Fig. [114]5c, d). These two sets of TF binding
   motifs have little overlap (Jaccard index = 0.05 between the two sets
   of 10 motifs with smallest enrichment Q-values), suggesting that the
   role of H3K4me3 dynamics in relation to TF binding has a
   tissue-specific function. The top scoring TF motifs in cancer cells
   that associated with sensitive H3K4me3 peak width dynamics tended to
   involve cancer-associated TFs such as c-MYC and ETS in cancer cells
   (Fig. [115]3e) and liver-specific TFs such as RXR and ESRRB in liver
   (Fig. [116]3f).

Fig. 3.

   [117]Fig. 3
   [118]Open in a new tab

   H3K4me3 width dynamics encode cell type-specific TF binding. a
   Distribution of fold changes (defined by the values in low methionine
   condition divided by the corresponding values in high methionine
   condition) in peak width in all peaks (gray), top 500 peaks with
   sensitive width (orange) and top 500 peaks with robust width (blue) in
   human cancer cells. b Framework for the TF binding motif enrichment
   analysis. c Distributions of TF binding motif enrichment Q-values in
   H3K4me3 peak sets with different dynamics under MR in human cancer
   cells. Box limits are the 25th and 75th percentiles, center lines are
   medians, and the whiskers are the minimal and maximal values. d Same as
   in c but for mouse liver. e Top 10 TF binding motifs enriched in 500
   H3K4me3 peaks with sensitive width in human cancer cells. f Same as in
   e but for robust width in mouse liver. g Number of oncogenes and tumor
   suppressors in top 10 TFs enriched in H3K4me3 peaks with sensitive
   width in human cancer cells and those with robust width in mouse liver

   To investigate the binding of these TFs, we obtained ChIP-seq datasets
   for two TFs putatively associated with H3K4me3 width in each system
   (Supplementary Fig. [119]6a) and computed the overlap between the TF
   binding sites and H3K4me3 peak subsets. We also chose one TF without
   binding motifs associated with peak width for each system as a negative
   control (Supplementary Fig. [120]6b). We found that the cancer-related
   TFs, MYC (Supplementary Fig. [121]6c), and ELF1 (Supplementary
   Fig. [122]6d) indeed had binding sites enriched (one-sided Fisher’s
   exact Q-value < 1e-10) in peaks with sensitive width in human cancer
   cells, while the liver-specific TFs, HNF4a (Supplementary
   Fig. [123]6e), and RXRa (Supplementary Fig. [124]6f–h) were also
   enriched (one-sided Fisher’s exact Q-value < 1e-10) in peaks with
   robust width in mouse liver. Such patterns were not observed in the TFs
   used as negative controls (one-sided Fisher’s exact Q-value = 1,
   Supplementary Fig. [125]6i, j). We also observed the strongest
   correlation between width dynamics and the number cell type-specific
   TFs bound to the genomic region containing the H3K4me3 peak (one-way
   ANOVA P-value = 3.53e-319 for width changes in human cancer cells and
   2.98e-189 for width changes in mouse liver, Supplementary Fig. [126]7).
   This concordance suggests that cell identity is encoded in the H3K4me3
   width dynamics but not in any other structural or locational aspect of
   H3K4me3. A previous study has found that broad H3K4me3 peaks are
   associated with tumor suppressor gene function in normal cells^[127]34.
   We also found that TFs with binding motifs enriched in preserved
   H3K4me3 peak width in liver tended to be tumor suppressors, while those
   enriched in sensitive H3K4me3 peak width in cancer cells tended to be
   oncogenes (chi-squared test P-value = 0.0324, Fig. [128]3g), suggesting
   that H3K4me3 width dynamics under MR may also encode information about
   tumor suppression or progression, depending on the cell or tissue type.
   Altogether, these analyses highlight H3K4me3 peak width dynamics as the
   information carrier under MR in both human cancer cells and mouse
   liver.

H3K4me3 width dynamics predict differential gene expression

   Despite a well-known association of the H3K4me3 mark with active
   promoters, the exact role of H3K4me3 in mediating gene expression is a
   matter of current debate^[129]45,[130]46. We therefore asked whether
   any aspects of genomic alterations in H3K4me3 track with changes in
   gene expression under MR. We first quantified gene expression levels
   using RNA sequencing (RNA-seq) under high and low methionine conditions
   and correlated the gene expression levels with H3K4me3 peak features
   again in both cancer cells and mouse liver (Fig. [131]4a). In cancer
   cells, we observed 7709 expressed genes and 5118 non-expressed genes
   marked by H3K4me3 (Supplementary Fig. [132]8a). There were also 11,412
   expressed genes without H3K4me3 marks, indicating in our models that
   the presence of the H3K4me3 mark is neither necessary nor sufficient
   for gene expression. Nevertheless, expressed genes with no H3K4me3 had
   significantly lower expression levels (Wilcoxon rank-sum
   P-value < 1e-323, Supplementary Fig. [133]8b), supporting the
   well-characterized association between H3K4me3 and active gene
   expression. This observation was further corroborated by significantly
   smaller H3K4me3 peak sizes in non-expressed genes (Wilcoxon rank-sum
   P-values < 1e-26 for height, area, and width, Supplementary
   Fig. [134]8c–e). An evaluation of the correlation coefficients between
   H3K4me3 peak height, area, width, and gene expression levels revealed
   significant positive correlations (Spearman’s rank correlation
   coefficient > 0.2, random permutation test P-value < 1e-50) between all
   three peak size descriptors and gene expression levels in both high and
   low methionine conditions (Fig. [135]4b, Supplementary Fig. [136]8f).
   Consistently, in mouse liver, we identified 9417 expressed and 2449
   non-expressed genes with H3K4me3, as well as 9282 expressed genes
   without H3K4me3 (Supplementary Fig. [137]9a). In accordance with our
   findings in cultured human cells, strong correlations between the
   presence of H3K4me3 and gene expression levels were also observed in
   liver (P-value < 1e-20 for all Wilcoxon rank-sum and random permutation
   tests, Fig. [138]4c, Supplementary Fig. [139]9b–f). Thus, in our
   models, the presence of an H3K4me3 peak, while not a requirement for
   gene expression, predicts overall whether a gene is expressed, and the
   magnitude of the peak does appear to contain information about overall
   gene expression level. This finding is consistent with the well-known
   association between H3K4me3 around transcription start sites (TSS) and
   active transcription^[140]47.

Fig. 4.

   [141]Fig. 4
   [142]Open in a new tab

   H3K4me3 width dynamics predict differential gene expression. a
   Framework of RNA-seq data analysis. Mouse image used with permission
   from Microsoft. b Spearman correlation between peak size descriptors
   and gene expression levels in human cancer cells. c Same as in b but
   for mouse liver. d Expression levels of genes associated with different
   H3K4me3 dynamics under high methionine conditions in human cancer
   cells. The gene set with sensitive width which has been demonstrated to
   associate with cell type-specific biological functions and TF binding
   is highlighted. Box limits are the 25th and 75th percentiles, center
   lines are medians, and the whiskers are the minimal and maximal values.
   e Same as in d but for mouse liver. f Differential gene expression in
   human cancer cells under MR. g Same as in f but for mouse liver. h
   Fraction of differentially expressed genes in genes with or without
   H3K4me3 in human cancer cells. i Same as in h but for mouse liver. j
   Spearman’s rank correlation coefficients between H3K4me3 changes and
   gene expression changes under MR in human cancer cells. k Same as in j
   but for mouse liver

   Having established a baseline for the relationship between H3K4me3 and
   gene expression, we next sought to study whether changes in the
   geometrical features of H3K4me3 are connected to changes in gene
   expression. We compared expression levels of genes associated with
   different H3K4me3 dynamics in cancer cells (Fig. [143]4d) and liver
   (Fig. [144]4e) under high methionine conditions. In addition to our
   previous findings that sensitive H3K4me3 peak width in cancer cells and
   robust H3K4me3 peak width in mouse liver are indicative of biological
   function, we found that genes associated with these peaks also exhibit
   significantly higher expression levels (Wilcoxon rank-sum and
   Kolmogorov–Smirnov P-values < 0.05, Fig. [145]4d, e, Supplementary
   Fig. [146]10). We then conducted differential expression analysis to
   identify the genes with altered expression in response to MR. In human
   cancer cells, we found 1034 genes upregulated (Wald test
   Q-value < 0.05, log[2](fold change) > 0) and 1161 genes downregulated
   (Wald test Q-value < 0.05, log[2](fold change) < 0) under MR
   (Fig. [147]4f) and 287 upregulated (Wald test Q-value < 0.05,
   log[2](fold change) > 0) and 431 downregulated (Wald test
   Q-value < 0.05, log[2](fold change) < 0) genes in mouse liver (Wald
   test Q-value < 0.05) (Fig. [148]4g) in diverse classes of genes. It is
   noteworthy that differential gene expression in vivo is confounded by
   factors including composition of liver by different cell types and
   larger variation between individual mice. Thus, as expected fewer genes
   were found to be differentially expressed in mouse liver.

   We next asked if changes to or consistencies in the geometrical
   features of H3K4me3, especially in peak width, can predict changes in
   gene expression. We first compared the fraction of differentially
   expressed genes that contain or are absent of H3K4me3. We found that in
   both cancer cells and liver, H3K4me3-marked genes were enriched in both
   upregulated and downregulated genes (one-sided Fisher’s exact
   P-value = 1.82e-12 for upregulated genes and 1.89e-98 for downregulated
   genes in human cancer cells, Fig. [149]4h, and 4.22e-4 for upregulated
   genes and 4.43e-16 for downregulated genes in mouse liver,
   Fig. [150]4i), suggesting that the presence of H3K4me3 notes a tendency
   of differential expression during MR. Next, we correlated fold changes
   in H3K4me3 peak height, area, width with changes in gene expression
   levels and found that H3K4me3 peak width dynamics significantly
   correlated (random permutation test P-value < 0.05, Spearman’s rank
   correlation coefficient > 0.1) with alterations in gene expression
   levels, and this correlation was consistent in both models (Spearman’s
   rank correlation coefficient = 0.14 for human cancer cells and 0.17 for
   mouse liver, Fig. [151]4j, k). On the other hand, although H3K4me3 peak
   height and area dynamics were also found to correlate with differential
   gene expression in mouse liver (random permutation test P-value < 0.05,
   Spearman’s rank correlation coefficient > 0.15, Fig. [152]4k), the
   strength of these correlations was smaller in human cancer cells
   (Spearman’s rank correlation coefficient < 0.06, Fig. [153]4j).
   Restricting the analysis to the subset of peaks located at promoters
   (Supplementary Fig. [154]11a) or altering the data analysis pipeline
   (Supplementary Fig. [155]11b) had minimal effects on the resulting
   correlation coefficients (absolute differences in Spearman’s rank
   correlation coefficients < 0.04 for data in Fig. [156]4j, Supplementary
   Fig. [157]11a, b). To further assess the predictability of changes in
   gene expression from changes in peak height, area, and width, we
   performed multiple linear regressions with changes in peak height,
   area, and width as the independent variables and changes in gene
   expression as the dependent variable. Change in peak width is the only
   variable with significant non-zero linear coefficients in predicting
   changes in gene expression in both human cancer cells and mouse liver
   (P-value = 3.75e-15 in human cancer cells and 8.39e-7 in mouse liver,
   Supplementary Fig. [158]11c, d). Interestingly, we also observed a
   stronger correlation between H3K4me3 width changes and gene expression
   changes in genes with H3K4me3 peaks bound by more TFs putatively
   associated with H3K4me3 width in both human cancer cells (Spearman’s
   rank correlation coefficient between width changes and gene expression
   changes = 0.19 in peaks bound by 2 TFs compared to 0.10 for 0 TF and
   0.12 for 1 TF, Supplementary Fig. [159]11e) and mouse liver (Spearman’s
   rank correlation coefficient between width changes and gene expression
   changes = 0.22 in peaks bound by 2 TFs compared to 0.18 for 0 TF and
   0.19 for 1 TF, Supplementary Fig. [160]11f), suggesting that TFs
   associated with H3K4me3 width dynamics regulate expression of their
   target genes thus mediating the connection between H3K4me3 dynamics and
   gene expression dynamics. Moreover, peaks with sensitive width in human
   cancer cells tended to be associated with downregulated genes
   (one-sided Fisher’s exact test P-value = 0.02, Supplementary
   Fig. [161]11g), which were also enriched with more cancer-related
   pathways (305 pathways among which 188 were cancer-related enriched in
   downregulated genes compared to 2 pathways enriched in upregulated
   genes, Supplementary Fig. [162]11h, i), suggesting that the correlation
   between H3K4me3 width dynamics and differential gene expression is
   linked to biological outcomes under MR. Similar pattern was also
   observed in mouse liver, in which upregulated genes were associated
   with peaks with robust width and enriched with liver-specific functions
   (Supplementary Fig. [163]11j–l). Taken together, H3K4me3 width
   dynamics, but not area or height, is the predictor of alterations in
   gene expression in both human cancer cells and mouse liver. Moreover,
   peak width also encodes information of gene expression levels under
   normal methionine conditions in addition to cell type-specific
   biological functions and TF binding preferences.