Epigenetic modifications are associated with inter-species gene expression variation in primates

Abstract

Background: Changes in gene regulation have long been thought to play an important role in evolution and speciation, especially in primates. Over the past decade, comparative genomic studies have revealed extensive inter-species differences in gene expression levels, yet we know much less about the extent to which regulatory mechanisms differ between species.

Results: To begin addressing this gap, we perform a comparative epigenetic study in primate lymphoblastoid cell lines, to query the contribution of RNA polymerase II and four histone modifications, H3K4me1, H3K4me3, H3K27ac, and H3K27me3, to inter-species variation in gene expression levels. We find that inter-species differences in mark enrichment near transcription start sites are significantly more often associated with inter-species differences in the corresponding gene expression level than expected by chance alone. Interestingly, we also find that first-order interactions among the five marks, as well as chromatin states, do not markedly contribute to the degree of association between the marks and inter-species variation in gene expression levels, suggesting that the marginal effects of the five marks dominate this contribution.

Conclusions: Our observations suggest that epigenetic modifications are substantially associated with changes in gene expression levels among primates and may represent important molecular mechanisms in primate evolution.

Figure 1

Marks are enriched near transcription start sites. (A) Fold enrichment of the five marks in ±2 kb regions near TSSs in the three primates. Error bars indicate standard deviation calculated across eight individuals in each species. Asterisks indicate significance levels based on Mann-Whitney one-sided tests (*P < 0.05, **P < 0.01, ***P < 0.001). (B) Distribution of normalized peak read counts for five marks around TSSs for each of the three primates. Units are in square root of RPKM (that is, RPKM^0.5) and are averaged across individuals and across genes.

Figure 2

Mark enrichment levels are correlated with gene expression levels in human. (A) Density of enrichment level for five marks around TSSs for genes with low, medium, and high expression levels. Values are averaged across individuals and across genes in each category. (B) Mark enrichment levels plotted against gene expression levels for sliding windows of genes (n = 200) ordered from low to high expression levels. Enrichment levels are obtained in ±2 kb regions near TSSs and scaled to be between 0 and 1. All values are averaged across individuals and across genes in the window. (C) Proportion of variance in gene expression levels explained (R squared) by individual marginal effects (five colored bars), combined mark marginal effects (grey bars), all first-order interaction effects in addition to marginal effects (black bars), and all chromatin state-specific effects in addition to marginal effects (white bars) of the five marks. Results are shown for enrichment levels in TSS regions with increasing length. Error bars indicate standard deviation calculated based on 20 split replicates.

Figure 3

Differentially expressed genes associate with inter-species differences in mark enrichment at transcription start sites. (A) Enrichment level differences for the five marks around TSSs of DE genes (black) and non-DE genes (grey) for each pair of species. Mark differences are considered with respect to the species associated with the lower gene expression level. DE genes are determined based on an FDR cutoff of 5%. (B) TSS regions associated with inter-species differences in any mark are enriched for DE genes. Plotted is the fold enrichment of TSS regions associated with inter-species differences in enriched marks in DE genes across pairs of species, for genes where the mark enrichment levels and the gene expression levels differ in the expected direction (that is, opposite for H3K27me3, same for the other four marks). Both the TSS regions associated with inter-species differences in enriched marks and DE genes are determined based on an FDR cutoff of 5%. Asterisks indicate significance levels from binomial tests (*P < 0.05, **P < 0.01, ***P < 0.001). C, chimpanzee; H, human; R, rhesus macaque.

Figure 4

Differences in mark enrichment level correlate with differences in gene expression level between pairs of primates. (A) Differences in mark enrichment level is plotted against differences in gene expression level for sliding windows of genes (n = 200) ordered based on the differential expression effect size, for all genes. Differences in enrichment level were obtained in ±2 kb regions near TSSs and scaled to be between -1 and 1. All values are averaged across individuals and across genes in the window. (B) Proportion of variance in gene expression level differences explained (R squared) by mark enrichment level differences, for all pairwise comparisons among the three primates. Different linear models are fitted to account for individual marginal effects (five colored bars), combined marginal effects (grey bars), all first-order interaction effects in addition to marginal effects (black bars), and all chromatin state-specific effects in addition to marginal effects (white bars) of the five marks. The DE genes are determined based on an FDR cutoff of 5%. Enrichment level differences are obtained in ±2 kb regions. Error bars indicate standard deviation calculated across 20 split replicates. C, chimpanzee; H, human; R, rhesus macaque.

Figure 5

Importance of marginal and interaction effects from five marks, and their enrichment in different chromatin states, for explaining gene expression level differences between primates. (A) The left panel lists marginal (M) or interaction terms (I2 to I5) among the five marks, where each row represents an interaction term and each column represents the presence (black) or absence (grey) of a particular mark effect for that interaction term. For example, the first row represents the marginal effect of H3K4me1, and the sixth row represents the interaction effect between H3K4me1 and H3K4me3. The right panel lists the corresponding PIP of each term between any pairs of primates for DE genes classified with different FDR cutoffs. (B) The left panel lists marginal (M) or chromatin state-specific terms for 15 chromatin states (S1 to S15) near TSSs, where each column represents the presence (black) or absence (grey) of a particular mark effect for that term. For example, the sixth row represents the state-specific effect of H3K4me1 in chromatin state S1. The right panel lists the corresponding PIP. The PIP measures the importance of each interaction term with higher values indicating higher significance. Mark enrichment level differences and mark enrichment level differences inside chromatin states within ±2 kb regions near TSSs were used for fitting. C, chimpanzee; H, human; R, rhesus macaque. M, marginal effects; I2, interaction term between pairs of marks; I3, interaction term among three marks; I4, interaction term among four marks; I5, interaction term among five marks; S1, active promoter; S2, weak promoter; S3, poised promoter; S4, strong enhancer; S5, strong enhancer; S6, weak enhancer; S7, weak enhancer; S8, insulator; S9, transcription transition; S10, transcription elongation; S11, weak transcription; S12, repressed; S13, heterochroma/lo; S14, repetitive/copy number variation; S15, repetitive/copy number variation.