Short-chain fatty acid metabolites propionate and butyrate are unique epigenetic regulatory elements linking diet, metabolism and gene expression

Abstract

The short-chain fatty acids (SCFAs) propionate and butyrate have beneficial health effects, are produced in large amounts by microbial metabolism and have been identified as unique acyl lysine histone marks. To better understand the function of these modifications, we used chromatin immunoprecipitation followed by sequencing to map the genome-wide location of four short-chain acyl histone marks, H3K18pr, H3K18bu, H4K12pr and H4K12bu, in treated and untreated colorectal cancer (CRC) and normal cells as well as in mouse intestines in vivo. We correlate these marks with open chromatin regions and gene expression to access the function of the target regions. Our data demonstrate that propionate and butyrate bind and act as promoters of genes involved in growth, differentiation and ion transport. We propose a mechanism involving direct modification of specific genomic regions by SCFAs resulting in increased chromatin accessibility and, in the case of butyrate, opposing effects on the proliferation of normal versus CRC cells.

Fig. 1. Genome-wide H3K18pr distribution.

a, H3K18pr versus H3K18ac differential binding at 10 mM propionate treatment. Sites identified as significantly differentially bound are shown in red; n = 3 technical replicates for each condition. Differential binding was performed using the DiffBind package with DESeq2, using a two-sided test for both increased and decreased binding affinity between conditions followed by multiple hypothesis testing and FDR correction. b, Top ten differentially bound regions associated with H4K12pr, annotated to within 1 kb of the TSS, sorted by FDR-adjusted P value (FDR < 0.05). c, Top GO biological process and molecular function terms of H3K18pr-associated cis-regulatory elements (5 + 1 kb) determined by GREAT against a whole genome background using a binomial test over genomic regions, followed by multiple hypothesis testing using FDR-corrected P values (FDR < 0.05). d, Normalized reads in H4K12pr-associated versus H4K12ac-associated binding sites at 10 mM propionate treatment. Box plots display the minimum, first quartile (Q1, 25th percentile; bottom of the box), median (line within the box, 50th percentile), third quartile (Q3, 75th percentile; top of the box) and maximum. The whiskers extend to the most extreme data points within 1.5 × IQR (interquartile range). e, Differential motif analysis of H3K18pr versus H3K18ac peaks was analysed by HOMER, using a one-sided hypergeometric test for overrepresentation (enrichment) of motifs in the target sequences compared to the background, followed by multiple hypothesis testing and FDR correction. f, Distribution of H3K18pr peaks by gene type was measured by two-sided chi-squared test to assess whether the observed distribution significantly differs from the control distribution, without specifying a direction (enrichment or depletion), followed by multiple hypothesis testing and FDR correction. g, Signal tracks for regions representing BAIAP2. Signal intensity of peaks in 95 kb-spanning BAIAP2 region showing H3K18pr versus H3K18ac binding at 10 mM propionate treatment with input as background. Source data

Fig. 2. Genome-wide H4K12pr distribution.

a, H4K12pr versus H4K12ac differential binding at 10 mM propionate treatment. Sites identified as significantly differentially bound are shown in red; n = 3 technical replicates for each condition. Differential binding was performed using the DiffBind package with DESeq2, using a two-sided test for both increased and decreased binding affinity between conditions, followed by multiple hypothesis testing and FDR correction. b, Top ten differentially bound regions associated with H4K12pr, annotated to within 1 kb of the TSS, sorted by FDR-adjusted P value (FDR < 0.05). c, Top GO biological process and molecular function terms of H4K12pr-associated cis-regulatory elements (5 + 1 kb) determined by GREAT against a whole genome background using a binomial test over genomic regions, followed by multiple hypothesis testing using FDR-corrected P values (FDR < 0.05). d, Normalized reads in H4K12pr-associated versus H4K12ac-associated binding sites at 10 mM propionate treatment. Box plots display the minimum, Q1, median, Q3 and maximum. The whiskers extend to the most extreme data points within 1.5 × IQR. e, Differential motif analysis of H4K12pr versus H4K12ac peaks was analysed by HOMER, using a one-sided hypergeometric test for overrepresentation (enrichment) of motifs in the target sequences compared to the background, followed by multiple hypothesis testing and FDR correction. f, Distribution of H4K12pr peaks by gene type was measured by two-sided chi-squared test to assess whether the observed distribution significantly differs from the control distribution, without specifying a direction (enrichment or depletion), followed by multiple hypothesis testing and FDR correction. g, Signal tracks for regions representing TTC7A. Signal intensity of peaks in 95 kb-spanning TTC7A region showing H4K12pr versus H4K12ac binding at 10 mM propionate treatment with input as background.

Fig. 3. TSS distribution profiles of H3K18ac/pr/bu-associated and H4K12ac/pr/bu-associated ChIP-seq peaks as a function of read coverage.

Upper panels, aggregate read density profile plots of genomic region distributions within ±1 kb of the TSS as a function of log₂(IP/input ratio). Lower panels, read density heatmaps of gene distributions with maximum (z = 4) and minimum (z = −4) values of heatmap intensities. Plots generated by deepTools.

Fig. 4. Propionyl and butyryl differential gene expression by RNA-seq.

a, Volcano plot showing gene upregulation versus downregulation in 10 mM propionate-treated versus control groups; n = 3 technical replicates for each condition. Differential expression analysis performed by DESeq2 b, Hierarchical clustering of GO biological process terms of upregulated versus downregulated pathways in propionate-treated versus control groups. Hierarchical clustering of the pathways was performed using ShinyGO. Pathways were clustered together based on shared genes. Gene enrichment analysis was performed using a two-sided Fisher’s exact test, and FDR correction was applied to adjust for multiple comparisons in the pathway analysis and hierarchical clustering. Size of dots indicates statistically significant FDR-adjusted (FDR < 0.05) P values. c, Heatmaps of the 50 most variable genes in propionate-treated versus control (Cnt) groups. d, Volcano plot showing gene upregulation versus downregulation in 1 mM butyrate-treated versus control groups. e, Hierarchical clustering of GO biological process terms of upregulated versus downregulated pathways in butyrate-treated versus control groups. f, Heatmaps of the 50 most variable genes in butyrate-treated versus control groups.

Fig. 5. Propionyl/butyryl differential RNA-seq.

a, Distribution of log₂(CPM)-transformed expression data for all conditions; n = 3 technical replicates for each condition. Normalization of raw counts performed by ‘cpm’ analysis in edgeR. CPM, counts per million. Box plots display the minimum, Q1, median, Q3 and maximum. The whiskers extend to the most extreme data points within 1.5 × IQR. b, Principal component analysis of log₂(CPM)-transformed expression data for all conditions. c, Volcano plot of propionyl versus butyryl differential expression. d, Heatmaps of the 50 most variable genes for all three conditions. e, Hierarchical clustering of GO biological process terms of differentially expressed pathways in propionyl versus butyryl RNA-seq. Hierarchical clustering of the pathways was performed using ShinyGO. Pathways were clustered together based on shared genes, and gene enrichment analysis was performed using a two-sided Fisher’s exact test. FDR correction was applied to adjust for multiple comparisons in the pathway analysis and hierarchical clustering. Size of dots indicates statistically significant FDR-adjusted (FDR < 0.05) P values.

Fig. 6. Murine cell line butyryl ATAC-seq and Kbu associated targets in mouse intestines.

a, Differential accessibility at 1 mM butyrate supplementation. Sites identified as significantly differentially accessible are shown in red; n = 4 technical replicates for each condition. Differential accessibility was performed using the DiffBind package with DESeq2, using a two-sided test for both increased and decreased binding affinity between conditions followed by multiple hypothesis testing and FDR correction. b, Top ten differentially accessible regions associated with butyrate supplementation sorted by FDR-adjusted P value (FDR < 0.05) and top GO biological process terms associated with positive versus negative fold change determined by GREAT against a whole genome background using a binomial test over genomic regions, followed by multiple hypothesis testing with FDR-corrected P values (FDR < 0.05). c, Normalized reads in binding sites at butyrate supplementation. Box plots display the minimum, Q1, median, Q3 and maximum. The whiskers extend to the most extreme data points within 1.5 × IQR. d, Correlation heatmap showing clustering of replicates from butyrate-supplemented versus untreated group. e, Signal tracks representing differential accessibility between butyrate-supplemented and untreated groups. f, Top GO biological process terms associated with H3K18bu in mouse intestines (5% arabinoxylan (AX)). g, Top GO biological process terms associated with H4K12bu in mouse intestines (5% AX). h, Annotated peak overlap between mouse intestinal H3K18bu (5% AX) targets and butyryl ATAC-seq in murine cells. i, Annotated peak overlap between mouse intestinal H4K12bu (5% AX) targets and butyryl ATAC-seq in murine cells.

Extended Data Fig. 1. Feature distribution of H3K18ac/pr/bu (top panel) and H4K12ac/pr/bu (bottom panel) associated regions (+/- 3 Kb of TSS).

H3K18ac/H4K12ac annotations were taken from ChIP-seq data that was generated without any treatment. H3K18pr/H4K12pr and H3K18bu/H4K12bu ChIP-seq experiments were performed following 10 mM NaPr, and 1 mM NaBu treatments, respectively. The x-axis provides the percentage of sites while colored regions represent distance from TSS. Results obtained using ChIPSeeker R package. Source data

Extended Data Fig. 2. Kac/pr/bu annotated features overlap.

a, d Hierarchical clustering of GO ‘Biological Process’ terms associated with Kac but not Kpr or Kbu. b, e Hierarchical clustering of GO ‘Biological Process’ terms associated with Kpr but not Kac or Kbu. c Hierarchical clustering of GO ‘Biological Process’ terms associated with Kbu but not Kac or Kpr. Hierarchical clustering of the pathways was performed using ShinyGO. Pathways were clustered together based on shared genes and gene enrichment analysis was performed using two-sided Fisher's exact test, and FDR correction was applied to adjust for multiple comparisons in the pathway analysis and hierarchical clustering. Size of dots indicates statistically significant FDR adjusted (FDR < 0.05) P values. Source data

Extended Data Fig. 3. Kpr/bu and ATAC/RNA-seq shared annotated features.

a H3K18pr, propionyl ATAC/RNA-seq shared annotated features. b Hierarchical cluttering of GO ‘Biological Process’ terms associated 782 shared features. c H4K12pr, propionyl ATAC/RNA-seq shared features. d Hierarchical clustering of GO ‘Biological Process’ terms associated with 897 shared features. Hierarchical clustering of the pathways was performed using ShinyGO. Pathways were clustered together based on shared genes and gene enrichment analysis was performed using two-sided Fisher's exact test, and FDR correction was applied to adjust for multiple comparisons in the pathway analysis and hierarchical clustering. Size of dots indicates statistically significant FDR adjusted (FDR < 0.05) P values. e H3K18bu butyryl ATAC/RNA-seq shared features. f Hierarchical cluttering of GO ‘Biological Process’ terms associated with 250 shared features. Source data

Extended Data Fig. 4. Effect of Na⁺ on differential accessibility and expression.

a ATAC-seq MA plot of differential accessibility in 10 mM NaCl-treated vs untreated cells (SW480). Differential analysis was performed by DESeq2 using two-sided Wald tests to identify differentially expressed or accessible genes. The results were then adjusted for multiple comparisons using the FDR correction method. Sites identified as significantly differentially accessible (FDR<0.05, log2 FC>1) are shown in red. n = 3 experimental replicates for each condition. b Volcano plot of differential accessibility in 10 mM NaCl-treated vs untreated cells. c RNA-seq MA plot of differential expression in 10 mM NaCl treated vs untreated cells. Sites identified as significantly differentially expressed (FDR<0.05) are shown blue (Insert: 1 mM vs Cnt). d Volcano plot of differential expression in 10 mM NaCl-treated vs untreated cells. Sites identified as significantly differentially expressed by FDR only (FDR<0.05) are shown blue (Insert: 1 mM vs Cnt). e Correlation heatmap of RNA-seq data showing clustering replicates from 1 and 10 mM NaCl-treated vs untreated groups. Normalized measurement of the covariance between replicates is expressed by Pearson’s correlation coefficient. Source data

Extended Data Fig. 5. H3K18bu vs ac and H4K12bu vs ac differential binding from mouse intestine on HSF + 5 % arabinoxylan diet.

a H3K18bu differential binding in cancer (SW480) vs normal (CCD841) cells and normalized reads in H3K18bu-associated binding sites following 1 mM NaBu treatment. Sites identified as significantly differentially bound are shown in red. n = 3 technical replicates for each condition. Differential binding was performed by DiffBind package with DESeq2 using a two-sided tests for both increased and decreased binding affinity between conditions followed by multiple hypothesis testing and FDR correction (FDR < 0.05). Box plots display: The minimum, first quartile (Q1, 25^th percentile), median, third quartile (Q3, 75^th percentile), and maximum. The bottom of the box is Q1 and the top of the box is Q3. The line within the box represents the median (50^th percentile) value. The whiskers extend to the most extreme data points within 1.5 times the IQR (interquartile range). b Top GO Biological Process terms for H3K18bu-associated cis-regulatory elements in cancer vs normal cells determined by GREAT against a whole genome background using a binomial test over genomic regions, followed by multiple hypothesis testing using FDR corrected P values (FDR < 0.05). c H4K12bu differential binding in cancer vs normal cells and normalized reads in H4K12bu associated binding sites following 1 mM NaBu treatment. d Top GO Biological Process terms for H4K12bu-associated cis-regulatory elements in cancer vs normal cells. e, f CRC differentially bound genes associated with H3K18bu/H4K12bu. g, h Signal tracks showing differential binding in MYC and PDGFA regions in cancer vs normal cells. Source data

Extended Data Fig. 6. H3K18bu vs ac and H4K12bu vs ac differential binding from mouse intestine on HSF + 5 % arabinoxylan diet.

a, b, d, e Sites identified as significantly differentially bound are shown in red. n = 3 biological replicates for each condition. Differential binding was performed by DiffBind package with DESeq2 using a two-sided tests for both increased and decreased binding affinity between conditions followed by multiple hypothesis testing and FDR correction. c, f Chromosomal positions of H3K18bu and H4K12bu-associated regions represented by red dots. Purple lines represent statistically significant enrichment compared to input. The genome was scanned with a sliding window (size 6 Mb) further subdivided into 2 equal-sized steps for sliding. Within each window a hypergeometric test was used to test for enrichment over all protein-coding genes in the genome. FDR adjusted P value cutoff for window was 1E-05. Gene chromosomal mapping performed using ShinyGO. g, h Feature distribution of H3K18bu and H4K12bu bound regions (+/− 3 Kb of TSS). The x-axis provides the percentage of sites while colored regions represent distance from TSS. Results obtained using ChIPSeeker R package. Source data

Extended Data Fig. 7. Overview of SCFAs propionate and butyrate as regulatory elements affecting histone binding, chromatin accessibility and gene expression.

Upper panels: In vitro cellular perturbations and in vivo fiber supplementation experimental designs. Bottom panels: Changes in TSS distribution profiles, differential chromatin accessibility and expression following SCFA supplementation. Top and Bottom volcano plots showing differential accessibility and expression following propionate and butyrate supplementation, respectively. NIAID Visual & Medical Arts. 09/26/2024. Grey laboratory mouse. NIAID BIOART Source. bioart.niaid.nih.gov/bioart/280. Source data