Analysis of Genetically Diverse Macrophages Reveals Local and Domain-wide Mechanisms that Control Transcription Factor Binding and Function

Abstract

Non-coding genetic variation is a major driver of phenotypic diversity and allows the investigation of mechanisms that control gene expression. Here, we systematically investigated the effects of >50 million variations from five strains of mice on mRNA, nascent transcription, transcription start sites, and transcription factor binding in resting and activated macrophages. We observed substantial differences associated with distinct molecular pathways. Evaluating genetic variation provided evidence for roles of ∼100 TFs in shaping lineage-determining factor binding. Unexpectedly, a substantial fraction of strain-specific factor binding could not be explained by local mutations. Integration of genomic features with chromatin interaction data provided evidence for hundreds of connected cis-regulatory domains associated with differences in transcription factor binding and gene expression. This system and the >250 datasets establish a substantial new resource for investigation of how genetic variation affects cellular phenotypes.

Keywords: chromatin structure; cis-regulatory domains; enhancer landscape; gene expression; genetic variation; macrophages; transcription factor binding.

Figure 1. Variation in mRNA expression scales with extent of genetic variation

A. Overview of experimental design and main data sets. B. Comparison of RNA-seq for polyA transcripts in BMDMs derived from the indicated mouse strains under notx. Log2(TPM+1) values are plotted for BALB, NOD, PWK and SPRET vs. C57 (TPM = transcripts per kilobase million). Transcripts exhibiting >2- or >4-fold changes (FDR < 0.01) are light blue and dark blue, respectively. C. WGCNA clustering of differentially expressed genes. Top functional annotations for each cluster are illustrated on the right. (See Figure S1 for modules). D. Ratio-ratio plots of the fold response to KLA in C57 vs. BALB or SPRET BMDMs. Blue dots show genes that are 4-fold reciprocal regulated. Green dots show a 4-fold stronger response to the KLA stimulus in one strain over the other. E. Relationship of differentially expressed genes to number of genetic variation. F. Expression comparison of 46 primary interferon stimulated genes in C57 and SPRET BMDMs under notx and 6h KLA conditions. Right column represents the SPRET/C57 gene expression ratio following KLA treatment.

Figure 2. Effects of genetic variation on nascent transcription are primarily in cis at promoter-distal locations

A. Comparison of GRO-seq gene body tag counts in BMDMs derived from the indicated mouse strains under notx. Log2(tag counts+1) values are plotted for BALB, NOD, PWK and SPRET vs. C57. Colors as in Figure 1A. B. Comparison of GRO-seq, 5′GRO-seq and H3K27ac signal at the Igf1 locus in BMDMs derived from each strain under notx conditions. C. Ratio-ratio plots of GRO-seq tag counts for KLA/notx conditions, comparing C57 vs. BALB or SPRET. Colors as in Figure 1D. D. Relationship of differential RNA-seq expression as a function of mutations between −30 and +20 bp of the TSS defined by 5′GRO-seq signal. E. Ratio-ratio plot of gene body GRO-seq tag counts in BMDMs derived from C57 and SPRET mice versus allele-specific tag counts in BMDMs derived from SPRET x C57 F1 mice. F. GRO-seq expression for Npy and Pla2g7 in BMDMs derived from C57 and SPRET mice and allele-specific tag counts in BMDMs derived from SPRET x C57 F1 mice.

Figure 3. Variation in TF binding greatly exceeds variation in gene expression

A. Scatter plots of log2 tag counts for H3K27ac ChIP-seq regions comparing C57 and BALB or SPRET. Colors as in Figure 1A. B. Scatter plots of log2 tag counts for ATAC-seq peaks passing IDR comparing C57 to BALB or SPRET. Colors as in Figure 1A. C. De novo motif analysis of distal (>3kb from TSS) ATAC-seq peaks associated with H3K27ac signal. Boxes display negative log10 p-values for enrichment of the motif and its rank order in parentheses. D. Pie chart indicating fractions of distal H3K27ac-positive regions of open chromatin occupied by PU.1, C/EBPβ and/or CJUN. E. Heat map of H3K27ac tag density at super enhancers. F. Comparison of log2 ChIP-seq tag counts for PU.1 in BMDMs derived from the indicated mouse strains under notx. Colors as in Figure 1A. G. Ratio-ratio plot of PU.1 ChIP-seq tag counts in BMDMs derived from C57 and SPRET mice versus allele-specific tag counts in BMDMs derived from C57 x SPRET F1 mice. H. SNPs + InDels frequencies in ATAC-seq and PU.1 peaks +/− 150 bp of the peak center for the indicated strain comparisons for stringent VCF filter criteria and more lenient criteria.

Figure 4. Effects of motif mutations enable inferences of a large network of collaborative TFs

A. Heat map of a subset of significant motifs after application of MARGE under notx and KLA treatment conditions (complete listing in Table S4). B. Top 14 of 48 motifs correlated with binding of PU.1 under notx as determined by motif mutation analysis. For highly related motifs (e.g., ETS factor motifs), the motif with the largest effect size is illustrated. Node size is fraction of PU.1 peaks containing the indicated motif and edge thickness is proportional to the effect size of motif mutations. Nodes indicate motifs in which mutations result in reduced PU.1 binding (red) or in which mutations result in increased PU.1 binding (blue). C. Top 15 out of 60 motifs correlated with binding of P65 under KLA treatment as determined by MARGE. Node size and edge thickness are defined in Panel B. D. Integrated network of collaborative TFs. The top 15 of 80 motifs for which motif mutations affected binding of at least one of the three factors are shown. Node sizes are the average fractional overlap of the indicated motif with PU.1, C/EBPβ or CJUN peaks and edges are factor-specific effect sizes. E. Fraction of >4-fold different strain specific binding of TFs explained by mutations in their respective recognition motifs and by all mutations considered by MARGE analysis. F. Overlap of binding of PU.1, RUNX1 and USF2 under notx as determined by ChIP-seq for each factor. G. Fraction of open chromatin marked by H3K27ac and not bound by PU.1, CJUN or C/EBP occupied by RUNX1, USF2 and/or NRF2. H. Relationship of mutations in RUNX motifs on binding of RUNX1 and PU.1 in C57 and SPRET BMDMs.

Figure 5. Clusters of ATAC-seq and ChIP-seq peaks are locally correlated

A–C. Heat maps of Pearson correlation coefficients (PCC) of PU.1, C/EBPβ and CJUN peaks, respectively, across the five strains under notx conditions in a 5 MB window from chromosome 18. Vertical and horizontal lines represent TAD boundaries as defined by C57 Hi-C assays presented in Figure 6. Axes represent sequential locations associated with the indicated feature, with the matrix values corresponding to correlation coefficients defined by the accompanying scale. D. Regional correlation of GRO-seq, ATAC-seq, PU.1, C/EBPβ, CJUN, and H3K27ac signal in the vicinity of the Colec12 gene. E. % of PU.1 CRDs based on minimum peak number and minimum PCC. F. Relationship of strain-specific PU.1 CRDs to enhancer activity measured by 5′GRO-seq and expression of nearest expressed gene measured by RNA-seq (> 16 tag counts). Significance was calculated using a two-sided t-test. G. Heat maps for relative binding and 5′GRO-seq signal at PU.1 CRDs. The ordering of PU.1 signal and corresponding 5′GRO-seq signal is the same for the two plots.

Figure 6. Correlated genomic features are highly connected as determined by proximity ligation assays

A. Hi-C contact frequency maps for chromosome 18 in BMDMs derived from C57 and SPRET mice. The values for the PC1 eigenvector are shown at the bottom (left) with zoomed-in view visualizing TAD boundaries (right). B. RNA-seq, H3K27ac, PC1, and Hi-C contact loops in the vicinity of the Spi1 and Colec12 loci. C Fraction of significant consensus PLAC-seq interactions within, between and outside of ATAC-seq CRDs. D. Example of ATAC-seq notx CRDs highly connected by PLAC-seq consensus interactions.