H3K18 lactylation marks tissue-specific active enhancers

Abstract

Background: Histone lactylation has been recently described as a novel histone post-translational modification linking cellular metabolism to epigenetic regulation.

Results: Given the expected relevance of this modification and current limited knowledge of its function, we generate genome-wide datasets of H3K18la distribution in various in vitro and in vivo samples, including mouse embryonic stem cells, macrophages, adipocytes, and mouse and human skeletal muscle. We compare them to profiles of well-established histone modifications and gene expression patterns. Supervised and unsupervised bioinformatics analysis shows that global H3K18la distribution resembles H3K27ac, although we also find notable differences. H3K18la marks active CpG island-containing promoters of highly expressed genes across most tissues assessed, including many housekeeping genes, and positively correlates with H3K27ac and H3K4me3 as well as with gene expression. In addition, H3K18la is enriched at active enhancers that lie in proximity to genes that are functionally important for the respective tissue.

Conclusions: Overall, our data suggests that H3K18la is not only a marker for active promoters, but also a mark of tissue specific active enhancers.

Keywords: Adipocyte; CUT&Tag; ChromHMM; Embryonic stem cell; Enhancer; Epigenetics; H3K18la; Histone post-translational modification; Lactate; Lactylation; Macrophage; Muscle; Promoter.

Fig. 1

H3K18la marks active CGI promoters. A Western blots showing H3K18la and H3 protein expression in all included samples (n = 3). The arrows indicate 15 kDa. B MDS of active hPTMs profiled from various mouse samples, quantified over 3000 bp genome-wide tiles. C Correlation heatmap of genome-wide, quantified H3K18la peak levels (biological replicates (n = 2) for mESC-ser, mESC-2i, ADIPO, GAS, PIM, and MT. For MB and BMDM; n = 1. Pearson’s correlation coefficient R is displayed as color gradient. D Distribution of H3K18la peaks of all mouse samples across genomic features. E H3K18la peak fold enrichment adjusted for genome-wide feature size for all mouse samples. F Venn diagrams depicting the promoter overlaps marked by the active hPTMs in mESC-ser, GAS, and PIM samples. Overlaps are colored according to the absolute number of promoters marked by various combinations of active hPTMs. Percentages indicate the fraction of actively marked promoters belonging to each group. G Scatterplots showing pairwise correlation of promoter H3K18la levels with other hPTM levels (log₂CPM) highlighting the promoters of genes with highest (red, n = 2000) or lowest (cyan, n = 2000) normalized gene expression (RPKM) for mESC-ser, GAS, and PIM. Pearson’s correlation coefficient R and p-values are indicated. H Normalized gene expression (log₂RPKM) per gene category is shown as boxplots. Gene categories are defined by the combination of active hPTMs as F. I Scatter plots showing the correlation between promoter H3K18la levels (log₂CPM, y-axis) and expression of the corresponding gene (log₂RPKM, x-axis) for mESC-ser, GAS, and PIM. Spearman’s correlation coefficient R and p-values are indicated

Fig. 2

H3K18la marks active, tissue-specific enhancers. A Tissue- and cell-type-specific ChromHMM analysis of mESC-ser, GAS, and PIM based on their hPTM profiles. The color scale corresponds to the emission parameter of each hPTM for each state. Fold enrichment of ChromHMM states for total genomic fraction coverage, genomic features, and ENCODE cCREs, scaled from −2 to 2 (see the “Materials and methods” section for details). B Heatmap showing the fold enrichment of hPTM peaks in ENCODE cCREs ; Σ (bp overlap)/[Σ (bp PTM peak)*Σ (bp cCRE)], scaled from −2 to 2. C Bar plots depicting the fraction of published tissue-specific enhancers [–37] that overlap with hPTM peaks. D Venn diagrams depicting the overlap of pELS/dELS marked by the active hPTMs in mESC-ser, GAS, and PIM samples. Overlaps are colored according to the absolute number of ELS marked by various combinations of active hPTMs. Percentages indicate the fraction of actively marked ELS belonging to each group. E ChromHMM analysis of all tissues/cell types based on their H3K18la profiles. The color scale shows the emission parameter of each tissue/cell type for each state. Fold enrichment of ChromHMM states over published tissue-specific enhancer sets [–37], total genomic fraction coverage, genomic features, ENCODE cCREs, house-keeping gene promoters, and house-keeping genes [38], scaled from −2 to 2 (see the “Materials and methods” section for details). F Top 10 GO (category “Biological Process”) terms resulting from the GO enrichment analysis of the genes closest to the top 2000 dELS from ENCODE cCRE with highest H3K18la levels (see the “Materials and methods” section for how dELS were linked to genes)

Fig. 3

Dynamic changes of H3K18la reflect transcriptional adaptations. A Box plots showing H3K18la log₂FC changes from MT versus MB over different genomic features. B Scatterplot showing the correlation between significant (FDR < 0.05) H3K18la log₂FC (>0.5) in promoters and their corresponding gene expression log₂FC (>0.5) based on the overlapping genes from MT versus MB differential analysis. Pearson correlation coefficient R is indicated. C IGV genome browser [50] snapshot of H3K18la profile at the Myh1 promoter region from various mouse samples and the corresponding SEACR-called peak regions. H3K18la level is depicted (rpm). Genomic regions are indicated on the top, as well as RefSeq gene names. D Fold enrichment of significant differentially H3K18la-marked peaks (FDR < 0.05, |log₂FC| > 1.5) in ENCODE cCREs. E Scatterplot showing the correlation between significant (FDR < 0.05) H3K18la log₂FC (>0.5) in dELS and their closest gene expression log₂FC (>0.5) based on the overlapping genes from MT versus MB differential analysis. Pearson correlation coefficient R is indicated. F Top 10 GO terms (category “Biological Process”) based on the GO analysis of the overlapping upregulated genes in MT from E (first quadrant red dots). G Box plots showing H3K18la log₂FC of peaks overlapping with MB- or MT-specific enhancers and of peaks not overlapping with these enhancers. Wilcoxon test p-values are indicated for each pair of groups. H Gene expression changes (log₂FC) after treatment of MB with 10 mM sodium-L-lactate of genes containing a H3K18la peak in MB, in MT, in both or in none of both. Wilcoxon test p-values are indicated for each pair of groups. *p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001

Fig. 4

H3K18la in human muscle marks active promoters and enhancers. A Correlation heatmap of hPTM levels in genome-wide quantified peaks from human muscle. Pearson’s correlation coefficient R is displayed as color gradient. B Distribution of all hPTMs across genomic features. C Peak fold enrichment adjusted for genome-wide feature size. D Venn diagram depicting the promoter overlaps marked by the active hPTMs. Overlaps are colored according to the absolute number of promoters marked by various combinations of active hPTMs. Percentages indicate the fraction of actively marked promoters belonging to each group. E Scatterplots showing pairwise correlation of promoter H3K18la levels with other hPTM levels (log₂CPM) highlighting the promoters of genes with highest (red, n = 2000) or lowest (cyan, n = 2000) normalized gene expression (RPKM). Pearson’s correlation coefficient R and p-values are indicated. F Normalized gene expression (log₂RPKM) per gene category is depicted as boxplots. Gene categories are defined by the combination of active hPTMs as in D. G ChromHMM analysis based on 5 hPTMs profiled in human muscle. The color scale shows the emission parameter of each mark for each state. Fold enrichment of ChromHMM states for total genomic fraction coverage, published skeletal muscle enhancers [57], different genomic features, and ENCODE cCREs, scaled from −2 to 2. H Venn diagrams depicting the overlap of dELS and pELS as marked by the different combination of the active hPTMs. Overlaps are colored according to the absolute number of ELS marked by various combinations of active hPTMs. Percentages indicate the fraction of actively marked ELS belonging to each group. I Heatmap showing the fold enrichment of hPTM peaks in cCREs (ENCODE); Σ (bp overlap)/[Σ (bp PTM peak)*Σ (bp cCRE)], scaled from −2 to 2. J Bar plots depicting the fraction of published human muscle enhancers [57] overlapping with the human hPTM peaks. K Top 10 GO terms (category “Biological Process”) resulting from a GO analysis of the corresponding closest genes to the 2000 dELS with highest H3K18la peaks (see the “Materials and methods” section for details on how enhancer and gene were linked)