Differences in gene expression have been observed in a variety of species (1
). However, the extent to which TF binding differences occur both within individuals and closely related species and the global relationship between TF binding and genetic variation are largely unexplored (4
). We used ChIP-Seq to map NFκB and PolII binding sites in ten humans: five are of European ancestry (including a parent-offspring trio), two of eastern Asian ancestry, and three of Nigerian ancestry (Table S1); nine of these have been analyzed by the HapMap (5
) and the 1000 Genomes (http://1000genomes.org
) Project, and one represents an individual for which high resolution SV maps are available (6
). All individuals but one were females; in pair-wise comparisons modest differences in TF binding were observed between the male and nine females; our analyses thus combined results from all ten. For comparison we also analyzed PolII-binding in a female chimpanzee.
We used stringent criteria to identify binding peaks (8
), and clustered them into discrete “binding regions” (BRs) (9
), yielding a total of 15,522 and 19,061 BRs for NFκB and PolII, respectively. Within BRs, most peaks were similar in position and magnitude among individuals (Fig. S1A). However, significant differences in binding were observed (Fig. S1A) and the Spearman correlation coefficients among replicates of different individuals (median values 0.79 and 0.90 for NFκB and PolII, respectively) were less than that of biological replicates of a given individual (median values 0.90 and 0.95, respectively; Fig. S2A, Table S2). 7.5% and 25% of the NFκB and PolII BRs, respectively, differed significantly between two individuals (ANOVA-test (9
), Bonferroni-adjusted P
-value <0.05; SOM and Fig. S3C), and many variable BRs exhibited >2 fold magnitude differences in binding (Fig. S3D). Variable BRs for both NFκB and PolII were often coassociated (P
<1e-4; permutation test; , Fig. S4), a correlation that is particularly strong for BRs <10kb apart (Fig. S4A). Variable NFκB and PolII regions were also often coassociated (P
=2.80E-25, Kolmogorov-Smirnov, Table S3; Fig. S4A), even though the NFκB and PolII data are from TNFα-treated and untreated cells, respectively. These results suggest adjacent binding sites and BRs may influence one another, perhaps through cooperative binding or interactions with other proteins.
Fig. 1 Effect of SNPs on NFκB and PolII binding. (A) Signal tracks of a NFκB motif and a TATA-box demonstrate effects of B-SNPs on TF binding, with correlations in the expected direction (i.e., with “correct trend”). (B) Fold (more ...)
For both NFκB and PolII, BRs within 1 kb of transcription start sites (TSSs) of RefSeq genes showed less variability (6% and 25%, respectively) than intergenic peaks (8% and 28%) (P
<1e-4; permutation-test); TSS BRs also revealed stronger ChIP-Seq signals (1.2 and 2.3 fold, respectively), with many exceptions (Fig. S5). The majority of binding regions (>70%) were occupied in two or more individuals, which argues against cell line artifacts (see Fig. S3B). The signal intensity for 40% and 53% of the BRs absent (i.e., “lost”) in one individual was similar to background for NFκB and PolII (9
), respectively, suggesting complete absence of binding in these cases, rather than “threshold effects”.
BRs differing in TF occupancy among individuals often involve loci of potentially high interest. These include the RPS26, BLK, SP140, and ZNF804A genes for PolII, which have been associated with type 1 diabetes, systemic lupus erythematosus, chronic lymphatic leukemia, and schizophrenia, respectively, and ORMDL3, PTGER4, and LOC253039 for NFκB, associated with asthma, Crohn’s disease, and rheumatoid arthritis (see SOM). Genes with variability in PolII binding showed a slight enrichment with immunity and defense functional gene categories (P
-value=0.045, Benjamini-Hochberg multiple testing correction) among target genes (9
We examined the genetic contribution to binding variation using SNPs from the 1000 Genomes Project. Individual SNPs in NFκB and PolII BRs frequently affected binding (, Fig. S6A), and the number of SNPs in BRs correlated with the frequency of significant binding differences (). SNPs altering the NFκB DNA binding motif had a strong effect, elevating the frequency of significant binding differences by 2.4 fold. ~90% of the binding differences followed the expected trend in which better matches to the consensus yielded higher binding signals (P<1e-3; see , Table S4, Fig. S6B). We call SNPs that putatively affect binding B-SNPs for Binding-SNPs.
We also searched for other associated DNA motifs, such as the Stat1 motif (previously associated with NFκB-binding (10
)), TATA-box, CAAT-box, and GC-box (11
) and also performed de novo
searches for enriched DNA motifs in BRs (9
), which revealed BR enrichments for the NFκB-motif and the GC-box, along with additional motifs (Fig. S7). We assessed the effect of genetic variation on each of the motifs. SNPs in the Stat1 motif markedly elevated the frequency of significant NFκB binding differences (1.3-fold enrichment; P
<1e-3, permutation-test; ), and 71% of the alterations in the Stat1 motif changed NFκB binding in the expected direction; i.e. improved Stat1 motif sequences increased NFκB binding (P
<1e-3; see , Table S4, Fig. S6B). For PolII, SNPs in the CAAT-box had a strong affect on binding (1.6 fold; P
<1e-3), with 63% of cases displaying the correct trend, whereas SNPs in the TATA-box and GC-box had modest effects (1.5 fold and 1.3 fold, with 51% and 52% exhibiting the correct trend). The significant covariance in the Stat1 motif with NFκB binding differences and the NFY CAAT-box with PolII binding suggests a functional interaction of Stat1 with NFκB and NFY and PolII, respectively; the latter has been documented previously (12
). We call this novel approach to examine coassociation of motifs with variable binding regions the Allele Binding Cooperativity test or “ABC test”.
We next analyzed the effect of SVs, >1kb genomic segments displaying copy-number variants (CNVs) or balanced inversions (6
). We probed high-density microarrays to identify CNVs in seven individuals ((9
) Table S5) and combined these with CNVs from another survey (14
). CNVs significantly elevated the frequency of BR differences between individuals by 5.1- and 2.0-fold for NFκB and PolII, respectively (P
<1e-4, permutation-test; , Fig. S8, Table S6). Furthermore, the effect followed the correct trend in 90% and 80% of the respective NFκB and PolII cases (); deletions reduced binding signals, whereas duplications elevated them. A combined set of high-resolution SVs identified by paired-end mapping (6
) also exhibited enrichment in binding differences for deletions intersecting with NFκB and PolII BRs (3.2 fold and 1.7-fold, respectively (P
<1e-4, permutation-test)). Importantly, we found a 2.8-fold significant enrichment for inversions on NFκB BRs (P
<1e-4, permutation-test), and a slight, non-significant enrichment for inversions on PolII BRs (), suggesting that inversions may affect binding. We called SVs associated with binding “Binding-SVs” (B-SVs).
Fig. 2 Effect of SVs on TF binding. (A) Example of a deletion affecting PolII binding. This example also shows a comparison of PolII occupancy in humans and a chimpanzee. A subset of individuals shares the chimpanzee binding phenotype. (B) Effect-sizes for microarray-based (more ...)
The total fraction of significant binding differences coinciding with genetic variations was 35% for NFκB and 26% for PolII (Table S7, Fig. S6C). 34% of the NFκB BRs intersect with SNP-differences between corresponding regions in different individuals (1% intersect with a known TF motif with SNPs falling both in the NFκB or the STAT1 motif; Table S8) and 3% with SVs (note: some SNPs coincide with SVs). Thus, genetic differences affecting the BR can be assigned to many, but not to the majority of, binding differences. Possible reasons for the remaining BR variation include trans-effects, epigenetic variation, as well as B-SNPs and B-SVs that were not ascertained. Some of the binding differences could be related to the different ages of the individuals.
We examined the effect of binding variation on gene expression by generation of deep RNA-Seq data from each cell line (9
) and comparison with binding data (, Fig. S9A). A significant correlation was observed (Spearman correlation coefficients of 0.475 and 0.461 for NFκB and PolII, respectively) (, Fig. S9B, Table S9), suggesting an influence of binding differences on mRNA abundance. Examples of correlated genes include UGT2B17, GSTM1
, and ZNF804A
, encoding glucuronic acid and glutathione transferases and a gene linked to schizophrenia (see SOM). However, a number of BR differences were not associated with differences in gene expression and presumably compensatory (e.g. feedback) mechanisms influence the expression in these cases. We also examined the effect of B-SNPs with differences in both binding and gene expression and found that both NFκB and PolII binding and expression differences correlated with the presence of B-SNPs, including those in the NFκB- and Stat1-motif (for NFκB) and CAAT-, GC-, and TATA-box (for PolII) (Spearman: 0.48–0.82; , Table S9). Copy number differences (i.e., B-SVs) also correlated with gene expression, albeit the correlation was not as strong as that of binding with gene expression (Table S10), indicating a more direct role for genetic variation on TF binding than on gene expression.
Fig. 3 Correlation and effect sizes of TF binding and gene expression. (A) Example showing a correlation of binding and expression. This figure also shows a transgression event, in which the daughter displays a strong increase in binding relative to the parents. (more ...)
The observation that SNPs and SVs are frequently associated with binding differences suggests a crucial role of cis elements in the genetics of TF binding. We thus analyzed the segregation pattern of BR occupancy in the parent-offspring trio, and observed potential Mendelian segregation in >90% of BRs (Fig. S10A), although this was difficult to determine with certainty as not all alleles relevant to TF binding have been ascertained in the parents. Interestingly, 947 and 732 BRs were occupied by NFκB and PolII, respectively, in the child but not in the parents, indicative of transgression, in which a binding event was evident only in the offspring; , Fig. S10B, Tables S11, S12, S13).
We also examined whether some BRs are specific to certain populations. Although the number of individuals analyzed is small, the NFκB data revealed a total of 14 BRs that were specifically occupied or unoccupied in the African or Asian individuals (Table S14). For PolII, the chimpanzee data was used to infer gains and losses relative to the likely ancestral state of binding, and a total of 68 population specific occupancies (gains and losses) were identified in the three population groups (see Table S14). Overall, we found relatively few population-specific events, ~0.1% to ~0.4%, suggesting that most alleles affecting TF binding are shared among different populations.
Since humans and chimpanzees exhibit 5–10% differences in gene expression (15
), we also examined divergence of TF binding among primates by analyzing PolII binding in a single chimpanzee. 15,418 (81%) of human BRs with corresponding syntenic regions in the chimpanzee genome were analyzed. The majority of PolII BRs were occupied both in humans and chimp (Fig. S11A). However, 32% of the BRs exhibited significant differences in binding (corrected P
-value <0.05; e.g. , ), a figure higher than that for human PolII variation (25%). Genes near regions uniquely occupied in the chimp were enriched in (i) nucleoside, nucleotide and nucleic acid metabolism; (ii) steroid metabolism (P
-values: 3.60E-05 and 4.16E-04, respectively). Furthermore, BRs uniquely occupied in humans were significantly enriched in protein modification and mRNA transcription (Fischer Exact test (9
), Benjamini-Hochberg P
-values: 2.22E-89 and 9.08E-139, respectively; Table S15).
Fig. 4 Comparison of PolII binding in humans and a chimpanzee. (A) Signal tracks for a peak found only in the chimpanzee. All ten individuals shown in Fig. S11B. B) Pie charts displaying occupancy by PolII of genomic regions where the chimp and human genomes (more ...)
As in humans, relative differences identified in the chimpanzee were higher in intergenic BRs relative to BRs within 1 kb of a TSS: 33% of the syntenic intergenic PolII BRs differed significantly from the human samples, compared to 31% near TSSs (P<1e-4; permutation test). Consequently, human BRs near TSSs were generally more likely to be scored as “occupied” in chimpanzee (81%) than intergenic BRs (46%; ). Furthermore, human BRs with strong binding signal (i.e., many mapped reads) are more frequently occupied in the chimpanzee than those with weaker signals (Fig. S11C), indicating either divergence of the weaker sites or signals that fell below the threshold at the low signal sites. Finally, we observed a general correlation between polymorphism and divergence in binding: i.e., variable BRs in humans displayed on average more divergence from chimpanzee BRs (in terms of fold-change in normalized read-counts) than non-variable BRs (Spearman 0.68; P=3.9e-07; see Fig. S11D).
Overall our data demonstrate extensive contributions of genetic variations on TF binding, many of which are expected to be functional through their affect on gene expression. Overall, the differences observed here (7.5% and 25% for NFκB and PolII, respectively, for humans; 32% for human/chimpanzee) greatly exceed estimates for sequence variation in coding sequences (estimated as 0.025% for humans (16
) and 0.71% for human/chimpanzee (17
)), suggesting a strong role for binding variation in human diversity. Extending mapping of B-SNPs and B-SVs for additional transcription factors will likely further inform on the genetic underpinnings of phenotypic diversity in humans.