Interplay of cis and trans mechanisms driving transcription factor binding and gene expression evolution

Abstract

Noncoding regulatory variants play a central role in the genetics of human diseases and in evolution. Here we measure allele-specific transcription factor binding occupancy of three liver-specific transcription factors between crosses of two inbred mouse strains to elucidate the regulatory mechanisms underlying transcription factor binding variations in mammals. Our results highlight the pre-eminence of cis-acting variants on transcription factor occupancy divergence. Transcription factor binding differences linked to cis-acting variants generally exhibit additive inheritance, while those linked to trans-acting variants are most often dominantly inherited. Cis-acting variants lead to local coordination of transcription factor occupancies that decay with distance; distal coordination is also observed and may be modulated by long-range chromatin contacts. Our results reveal the regulatory mechanisms that interplay to drive transcription factor occupancy, chromatin state, and gene expression in complex mammalian cell states.

Fig. 1

F1 mice were used to interrogate the regulation of TFBS variation. a In vivo binding of liver-specific TFs FOXA1, HNF4A, and CEBPA were profiled in the livers of male mice from inbred strains C57BL/6J (BL6), CAST/EiJ (CAST), and their F1 crosses: C57BL/6J × CAST/EiJ (BL6xCAST) and CAST/EiJ × C57BL/6J (CASTxBL6). Six biological replicates were generated for each TF and genetic background combination. b The number of TFBS that could be classified with associated number of SNVs. c Venn diagram illustrates the numbers of classifiable SNVs that overlap between TFs. Each variant is at least 250 bp from any other SNV. Numbers shown are the final numbers of regulatory loci used for downstream analyses. d Heatmap confirming overall accuracy of regulatory class assignments. BL6 (black) vs. CAST (brown) binding intensity ratios for different regulatory categories for CEBPA. A subset of variants from each class was randomly sampled to match the overall distribution. Sparkline in key shows the number of observations at each color category, where density is increasing from left to right

Fig. 2

Differences in TF binding intensities strongly affected by variation acting in cis and are additively inherited. a Mean F0 vs. F1 TF binding intensity ratios (BL6 vs. CAST) for CEBPA are plotted in the left panel. The right panel shows mean F0 vs. F1 gene expression ratios for liver-expressed protein-coding genes. The correlation coefficient reflects the extent of cis-directed regulatory mechanisms. b Proportion of CEBPA binding locations at promoters and enhancers. The width of the bar is proportional to the overall number of TFBSs in the “All” category. Binomial tests were used to test for enrichment at promoters and enhancers for each regulatory class based on the overall numbers of TFBSs (“All”). ***P < 0.0001, **P < 0.001,* P < 0.05. c Most allele-specific TFBSs are affected by variation acting in cis. Lineage-specific TFBSs were defined as TFBSs where binding occurs either in BL6 or CAST in F0 individuals and in an allele-specific manner in F1 individuals based on a cut-off (F0_B6/(B6+CAST) > 0.95, F1_B6/(B6+CAST) > 0.95, F0_B6/(B6+CAST) < 0.05, F1_B6/(B6+CAST) < 0.05). These TFBSs can be sorted into the three categories described. d Mean CEBPA log2 F0 total read counts were plotted against mean log2 F1 read count (BL6 + CAST allele) multiplied by 2. For the scatterplot, we used averages across biological replicates. TFBSs affected by variation acting in cis are thus expected to fall along the diagonal and these have been colored blue (see c). Categories shown were determined by maximal likelihood estimation. e The majority of cis-directed TFBSs are inherited additively. TFBSs affected by variation acting in trans may show additive or dominant inheritance patterns in TF binding intensities. Different modes of inheritance were defined by comparing overall peak binding intensities between F0 and F1 individuals. Total F1 counts were individually scaled to 1 (yellow). Red indicates TFBSs where F1 > F0; blue indicates TFBSs where F1 < F0. CEBPA data are shown

Fig. 3

Effect of genomic distance on cis-acting inter-peak correspondence. a Strategy for measuring the span of cis regulatory effects. Successive 1 kb bins were taken from each TFBS affected by variation acting in cis starting 400 bp from the location of the SNV and extending in both directions. For each bin, Spearman’s ρ was calculated using the BL6:CAST allelic ratio between queried TFBSs against TFBSs assigned as anchorages for the analysis. b Spearman’s ρ values for each bin were plotted for each TF. The linear regression line (solid red) calculated from these values is shown. Red dashed lines mark the 90% confidence intervals of the true slope of the line. Gray dots represent the null background distribution of correlation values constructed by the random subsampling of TFBSs to anchor TFBSs (see Methods section). The numbers of TFBSs in each randomly sampled bin were matched to those in the observed bins. The gray line is the linear regression line for the correlation values derived from sampled points. c TFBSs are enriched at regions of chromatin contact. Enrichment values were calculated compared with expected rate of chromatin contact given the general enrichment for contact in each regulatory dataset (i.e., cons, trans, cis, cis–trans). “Any” denotes the null background set of randomly chosen locations in the genome

Fig. 4

Genetic and epigenetic influences that change TF binding have parallel consequences for gene expression and chromatin. a, b Coordination between the regulatory categories of variation in TF binding occupancy variation and chromatin (a) and gene expression (b). Locations of the considered TFBSs are noted in the cartoons on the left. Separate logistic regressions were performed for each chromatin regulatory class (see Methods section). Odds ratios were mean-centered for comparison across chromatin regulatory classes. Absolute values of Z-scores >2 (α < 0.05) are denoted by a black border. c Direct association between chromatin and gene expression. Genes were linked to H3K4me3 modifications if the mark was located within 5 kb upstream of the TSS. Binomial tests were performed based on the expected background probability of observing the same regulatory mechanism underlying both expression and histone enrichment change. d High diversity in regulatory mechanisms of TF binding variation is associated with gene expression influenced by cis–trans-acting variation. Calculations are on a gene-by-gene basis for TFBSs 20 kb upstream and 10 kb downstream of TSSs. These scores were compared between genes grouped by transcriptional regulatory class. Significant P-values for Mann–Whitney U-tests are shown. The surface area of the violin plot is proportional to the number of genes in each class