Massively parallel decoding of mammalian regulatory sequences supports a flexible organizational model

Abstract

Despite continual progress in the cataloging of vertebrate regulatory elements, little is known about their organization and regulatory architecture. Here we describe a massively parallel experiment to systematically test the impact of copy number, spacing, combination and order of transcription factor binding sites on gene expression. A complex library of ∼5,000 synthetic regulatory elements containing patterns from 12 liver-specific transcription factor binding sites was assayed in mice and in HepG2 cells. We find that certain transcription factors act as direct drivers of gene expression in homotypic clusters of binding sites, independent of spacing between sites, whereas others function only synergistically. Heterotypic enhancers are stronger than their homotypic analogs and favor specific transcription factor binding site combinations, mimicking putative native enhancers. Exhaustive testing of binding site permutations suggests that there is flexibility in binding site order. Our findings provide quantitative support for a flexible model of regulatory element activity and suggest a framework for the design of synthetic tissue-specific enhancers.

Figure 1

Synthetic enhancer sequence design and controls. (a) SRESs consist of patterns of 12 consensus binding sequences arranged homotypically (class I) or heterotypically (class II and class III) on 1 of 2 neutral, 168-bp templates. (b) Schematic of massively parallel reporter assay methodology. SRESs were cloned upstream of a minimal promoter in a tagged luciferase library and then assayed in vivo using hydrodynamic tail vein injection. Livers were dissected 24 h after injection, mRNA was generated, and tags were reverse transcribed and sequenced. (c) Bimodal distribution of expression values for 4,966 SRESs. Expression values were calculated using the equation shown. (d) Template-template correlation. Expression values for 2,217 pairs of SRESs (not all SRESs had data for both templates owing to quality control measures) containing the exact same patterns of consensus binding sequences on 2 separate templates are plotted. The red line is a linear regression trace, whereas the dashed line is the diagonal. Template 1, hg19 chr. 9: 83,712,599–83,712,766; template 2, hg19 chr. 2: 211,153,238–211,153,405. (e) Expression values from the three mice in which the SRES library was tested, exhibiting a very high level of correlation.

Figure 2

Homotypic amplification of expression is compatible with a subset of transcription factor binding sites, independent of their spacing. (a) We observed significant correlation between expression and the size of the homotypic cluster for 5 of the 12 transcription factor binding sites (CEBPA, FOXA1, HNF1A, ONECUT1 and XBP1). The PPARA binding site is included as an example of a site that could not be homotypically amplified. Included on the right are box plots for the background expression of all SRESs with a single binding site (B), as well as for the positive (+) and negative (−) controls. Red boxes denote groups of SRESs with significantly higher expression compared to background (Wilcoxon rank-sum test P ≤ 0.05), which is a slightly more stringent test than comparison against negative controls. P values refer to Spearman's correlation coefficients (corrected for multiple testing using FDR). In the box plots, the central rectangle spans the first and third quartiles, the line inside the rectangle is the median, and the lines beyond the box indicate the locations of the minimum and maximum values. (b) In the vast majority of cases, the strength of expression was not dependent on the distance between binding sites, as observed for class I elements. Shown are examples of SRESs, each with two copies of one of the three strongest transcription factor binding sites, including sites for CEBPA, HNF1A and XBP1. P values refer to Spearman's correlation coefficients, and the dashed gray lines are the regression traces.

Figure 3

Heterotypic elements drive stronger expression than homotypic ones. (a) Density of expression by SRES class. The density of expression is plotted for each class of SRES. Red dashed lines denote the mean expression value for each class. (b) Box plot of expression by class and number of sites. Note that there are no class II SRESs with 1 site or class III SRESs with <3 sites. In the box plots, the central rectangle spans the first and third quartiles, the line inside the rectangle is the median, and the lines beyond the box indicate the locations of the minimum and maximum values.

Figure 4

Combinatorial effects in heterotypic clusters of two different transcription factor binding sites. (a) Significant interactions that were identified from a general model of class I and class II SRES–driven expression using the χ² goodness-of-fit test or direct comparisons between class I and class II data sets. Dotted gray lines refer to the 19 combinations (of a possible 121) that were sampled in class II SRESs. Red lines indicate significant synergy, and a black line indicates significant interference (P values ≤ 0.05 in all cases). (b) Direct comparisons between the eight pairs of interacting binding sites and predicted expression based on class I data alone (black lines). Five interactions (NR2F2-XBP1, NR2F2-ONECUT1, NR2F2-FOXA1, RXRA-XBP1 and HNF1A-XBP1) were identified by including combinatorial terms in the model (Wald χ² test P values shown at top right), whereas three (FOXA1-RXRA, FOXA1-PPARA and TFAP2C-RXRA) were identified by directly comparing expression from binding site pairs in combinations with homotypic sequences containing an equal number of the binding sites in isolation (Wilcoxon rank-sum test). A single example is shown for each combination, corresponding to a fixed number of sites (the number of sites in each SRES is given in parentheses, and asterisks to the left of P values indicate those that are significant). In the box plots, the central rectangle spans the first and third quartiles, the line inside the rectangle is the median, and the lines beyond the box indicate the locations of the minimum and maximum values.

Figure 5

Synthetic enhancers mimic mouse liver enhancers. (a) Categorization of 51,850 putative mouse liver enhancers into the 3 SRES classes on the basis of the distribution of the 12 transcription factor binding sites. (b) Frequencies of transcription factor binding site pairs for each of the 8 interactions identified in Figure 4 in putative mouse liver enhancers versus 103,700 matched random genomic controls.

Figure 6

Effects of binding site order in heterotypic enhancers. (a) Depiction of the 10 combinations (out of 39 total combinations) of 3 transcription factor binding sites with a favorable permutation resulting in significantly stronger expression (FDR < 0.05, Wilcoxon rank-sum test). The difference in expression between the best and worst permutation is depicted on the left. (b) Rank-value plot depicting 9 sets of SRESs containing ~49 permutations, each of the same 8 transcription factor binding sites. (c) HepG2 and mouse liver SRES expression strongly agree for the entire set of 441 permutations of 8 transcription factor binding sites. The red line is the regression trace. A specific example of the best and worst permutation for one of these sets (shown in purple in the plots in b,c) appears at the bottom.