CCRA libraries are created by array-based oligonucleotide synthesis (Agilent). Each element of the library is a distinct 230 bp oligonucleotide comprised of five different sequence regions (See Supplemental Note 2 for a diagram of these regions and the specific sequences used in this study). The first region is a 20 bp constant sequence that is homologous to the backbone plasmid to support Gibson cloning. The next (downstream) 11 bp sequence is unique to each sub-library to enable the amplification of subsets of the library elements that are synthesized in each batch. This allows for the analysis of different TFs or the testing of different hypotheses using a single oligonucleotide synthesis. The third region is the 170 bp user-defined variable synthetic promoter sequence. This region is followed by 12 bp ‘promoter’ barcode that identifies the corresponding promoter sequence at Illumina sequencing step. Each promoter barcode is designed to be at least 3 bp different than all other barcodes to control for synthesis, PCR and sequencing errors. The last region of each library element is a constant 17 bp sequence used for PCR amplification. The library pool was synthesized by Agilent as 10 pmol of lyophilized nucleic acid (See Supplemental Note 2 for primer design and additional details for library amplification). To amplify the library, we used 0.15 ng of library DNA template in a final 50 μl PCR reaction. In each 50 μl reaction, we used 0.2 mM dNTP mix, 0.5 μM forward primer, 0.5 μM reverse primer, 1× Herculase II reaction buffer, 1 M Betaine, 0.15 ng DNA template in water, 1 μl of Herculase II polymerase (Agilent). The PCR reaction was cycled as follows: 95° for 1 min, 16 cycles of 95° for 30 s and 58° for 2.5 min and then 72° for 4 min. PCR products were purified by AMPure XP beads from Beckman coulter with 1:1.6 of PCR sample to magnetic particles ratio according to manufacturer's instructions. Typically, we obtained 5 to 10 ng/μl of DNA in a final volume of 15 μl.

Plasmid pRS414 was used as the backbone to create library plasmid pRM1806 (See Supplemental notes 1a for plasmid map, sequence, and Supplemental Table S2 for Addgene accession number). To clone library sequences into the pRM1806 backbone, we linearized the plasmid with high fidelity KpnI and SacI (NEB), and then performed gel extraction using the Qiagen DNA extraction kit. We used 0.03 pmol of the linearized plasmid and 0.12 pmol of purified PCR product in a Gibson assembly reaction (NEB), following the manufacturer's instructions. Nitrocellulose membrane (0.025 μm) was used to filter Gibson assembly product by drop dialysis following the Millipore Sigma protocol. The library was electroporated into 10G SUPREME Electrocompetent cells (Lucigen) using 0.1 cm cuvette and cells were plated on to Kanamycin containing LB plates after 1-hour recovery in SOC. After 16 h of growth, over 50 000 colonies were scraped and the plasmid DNA was extracted using Qiagen Miniprep Kit.

The yeast strain used in this study was yRM1004, which is derived from matA_deltaSir4, and has the following genotype: his3Δ1 leu2Δ0 met15Δ0 ura3Δ0 Δsir4::KanMx Δtrp1::HygMx. Induction of TF directed transposition was performed using a modified calling cards protocol (24). Briefly, plasmid containing a Sir4p (amino acids 951–1200) tagged TF driven by ADH1 promoter with LEU2 auxotrophic marker was transformed into yeast cells (yRM1004) together with the plasmid pRM1804 (see Supplemental notes 1b for plasmid maps, sequences, and Supplemental Table S2 for Addgene accession numbers), which contains the URA3 marker and a galactose inducible Ty5 transposon with an artificial intron inside of His3 gene that is inside of Ty5 gene body for the purpose of selecting transposition positive cells in the next step (28). After transformation, cells were plated onto a Glu-Ura-Leu plate to select for cells carrying both the TF-sir4p fusion plasmid and Ty5 transposon plasmid. Next, a single colony was picked for library plasmid transformation. The library plasmid pRM1806 carries the TRP auxotrophic selection marker, so after the yeast cells were transformed with the library plasmid, they were plated onto a Glu-Ura-Leu-Trp plate to select for all three plasmids. Multiple parallel transformations were performed to obtain a diverse population of library sequences. We typically obtained over 10 000 colonies for each sub library. All colonies were pooled and plated to Gal-Ura-Leu-Trp to induce Ty5 transposition on 10 plates to increase the number of transpositions. Cells were allowed to grow on galactose plates for 4 days at room temperature. After galactose induction, we replica plated cells to Glu-His-Trp to select for yeast with Ty5 transpositions and that carry the library plasmid. After 2–3 days, colonies were scraped, and plasmid extraction was performed using the Yeast Plasmid Mini Kit (Omega).

We performed four independent PCRs to recover transpositions that were inserted into synthetic promoters in either of two possible orientations and upstream or downstream of the barcodes and UMI. We performed an additional PCR to measure the relative abundance of elements in the library for normalization. For these four PCRs, one primer of each pair is specific to either 3′ LTR of Ty5 transposon sequence or 5′ LTR of Ty5 transposon sequence, and the other primer is specific to a constant region either upstream or downstream of the inserted library sequence on the plasmid. For the additional PCR, one primer is specific to an upstream constant region of the inserted library sequence on the plasmid, and the other primer is for the downstream constant region (See Supplemental Table S1 for the primer sequences used). All five PCR products were pooled together for sequencing.

In each PCR reaction, we used 1X RedTaq buffer, 0.2 mM dNTP mix, 1M Betaine, 0.5 μM forward primer, 0.5 μM reverse primer, 4 μl RedTag DNA polymerase (Sigma-Aldrich), 1 μg of the purified plasmid DNA and the corresponding amount of water to reach a final volume of 50 μl. The PCR parameters were set to be 93° for 2 min, 24–28 cycles of 93° for 30 s and 62° for 6 min and 62° for 6 min. The PCR products were then purified with Qiagen PCR purification kit before sequencing.

After transforming the library plasmid into yeast, we divided the cells for either Calling cards or Sort-Seq. For expression measurement, we followed the experimental procedures as well as promoter expression calculation described in (18). We sorted cells into eight bins of 100 000 cells each, and then added yeast culture media to grow the cells for 16 h. Cells from each bin were then pelleted separately and the plasmids were extracted with Yeast Plasmid Mini Kit (Omega) for sequencing.

Next, we performed a separate PCR reaction for each sorted bin. The primer sequences are listed in Supplemental Table S1, and they target the constant regions upstream and downstream of the CCRA library. In each of the eight PCR reactions, the reverse primer was indexed with unique barcode to allow the reactions to be sequenced together. The PCR amplification conditions used were identical to those used for calling cards recovery.

To quantify TF binding to CCRA libraries, we analyze Illumina paired end sequencing reads to count all unique insertions into each library member. An transposition is unique if it can be distinguished by its insertion coordinate relative to the library reference or contains a unique UMI in instances where multiple insertions have landed at the same position across four independent PCRs (See Supplemental notes 3 for examples of CCRA sequencing reads; our python analysis source code and an output sample are available on Gitlab). To identify unique insertions from the sequencing data, we first filter for reads containing the appropriate 12 bp library barcode and 6 bp TF barcode. Filtered reads are then divided into five categories: reads from synthetic promoters where the Ty5 transposon inserted in the forward direction upstream of the promoter barcodes, reads where Ty5 inserted in reverse direction upstream of barcodes, reads where the Ty5 inserted in forward direction downstream of barcodes, reads where the Ty5 inserted in reverse direction downstream of barcodes, and reads from synthetic promoters without insertion. This categorization is achieved by analyzing the first 20 bp of read 1 and read 2. The next 12 bp are used to map the precise location of the transposon insertion into the synthetic sequence. We used the 4 bp UMI to resolve events when multiple calling cards are deposited at the same base pair in a given synthetic sequence. Finally, we use the number of full-length sequences recovered for each library element as a normalization factor to control for the variation in abundance between library members. The total number of independent insertions for each library member is normalized by the relative abundance of each element in the library to compute a normalized binding score (NBS) of TF binding to each synthetic sequence.

For experiments in which changes in binding energies are measured, it is important to measure TF binding strength as accurately as possible. Therefore, we used an expectation maximization algorithm to resolve TF-directed transpositions which occur near TF recognition sites from background transpositions which occur uniformly across the synthetic promoter. Since the distribution of TF directed insertions is approximately Gaussian with the distribution centred at the TF recognition site, we assumed that TF directed insertions can be modelled with this distribution while background insertions follow a uniform distribution. We then used an expectation maximum algorithm to estimate, for each synthetic promoter, the variance of the Gaussian distribution (the mean value is determined by the location of the TF recognition sequence) and the fraction of insertions that were the result of a TF-directed or background transposition. For each library element, we iterate each independent insertion for maximum of 1000 times or until the parameters no longer change. The estimated fraction of TF-directed insertions is used to multiply the raw number of insertions at each promoter to remove insertions due to non-specific transposition (Python script and sample files are provided in Gitlab). This background correction step removes 0–20% of non-specific insertions, which is important for calculating small changes in binding energy; however, incorporating this step does not impact other analysis and should not be used for sequences where the Gaussian assumption is not appropriate (e.g. for sequences with multiple TF sites or for TFs whose recognition sequence is not well-characterized). Therefore, we performed this background correction only for the generation of binding energy landscapes.

To quantitatively compare CCRA with PBM and MITOMI in terms of binding affinity, we calculated the change of binding energy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta G$\end{document}) from consensus site to the alternative site as follows:

Under binding equilibrium, [TF] and [sequence] associate at the same rate that the bound complex [TFS] disassociates:(1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}[{\rm{TF}}] + [{\rm{S}}] \longleftrightarrow [{\rm{TFS}}]\end{equation*}$$\end{document}

The Gibbs free energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta G$\end{document} is related to the binding constant K as follows:(2a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}{{K({\rm S})}} = \frac{{[{\rm TF}][{\rm S}]}}{{[{\rm TFS}]}} = {{\rm e}^{\Delta G/RT}}\end{equation*}$$\end{document}(2b)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}\Delta {{G}} = {{RT\\ {\rm ln}}}({{K}}({\rm{S}}))\end{equation*}$$\end{document}

The binding occupancy on a sequence is defined as the fraction of bound sequence to the total sequence in solution. Replace [TFS] with [TF][S]/K according to 2a, and by approximation that K(S) is much greater than [TF] as the affinity of these sequences are high, we get:(3a)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}{\rm{Occ}}({\rm{S}}) = \frac{{[{\rm TFS}]}}{{[{\rm TFS}] + [{\rm S}]}} = \frac{{[{\rm TF}]}}{{[{\rm TF}] + K({\rm S})}} \approx \frac{{[{\rm TF}]}}{{K\left( {\rm S} \right)}}\end{equation*}$$\end{document}(3b)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}{\rm{K}}({\rm{S}}) = \frac{{[TF]}}{{Occ(S)}}\end{equation*}$$\end{document}

Therefore, the change of binding energy equals:(4)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{eqnarray*}\Delta \Delta G &=& \Delta G({\rm Sconsensus}) - \Delta G({\rm Smutant})\nonumber\\ &=& - {{RT{\rm ln}}}\left( {\frac{{{\rm Occ}({\rm Smutant})}}{{{\rm Occ}({\rm Sconsensus})}}} \right)\end{eqnarray*}$$\end{document}

To determine if Cbf1p binds cooperativity at various synthetic promoters, we compared the observed occupancy to expected occupancy assuming independent binding, and we derived this test by the following:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{eqnarray*}[{\rm{Cbf}}1{\rm{p}}] &+& [{\rm{DNA\,with\,two\,free\,sites}}]\overset {k1} \longleftrightarrow [{\rm{Cbf}}1{\rm{p}} - {\rm{DNA\,with\,one\,free\,site}}] \\ &+& [{\rm{Cbf}}1{\rm{p}}]\overset {k2} \longleftrightarrow [2*{\rm{Cbf}}1{\rm{p}} - {\rm{DNA\,with\,both\,sites\,occupied}}]\end{eqnarray*}$$\end{document}

To simplify: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$[{\rm{P}}] + [{\rm{S}}]\overset {k1} \longleftrightarrow [{\rm{PS}}] + [{\rm{M}}]\overset {k2} \longleftrightarrow \ [{{\rm{P}}_2}{\rm{S}}]$\end{document}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{equation*}{\rm{Occ}}({\rm{P}}) = \frac{{2*K1*P + 2*K1*K2*{P^2}}}{{1 + 2*K1*P + K1*K2*{P^2}}}\end{equation*}$$\end{document}

If Cbf2 binds additively, then K1 = K2 = K;\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{eqnarray*}{\rm{Occ}}({\rm{P}}) &=& \frac{{2*K*P + 2*{k^2}*{P^2}}}{{1 + 2*K*P + {k^2}*{P^2}}} \\ &=& \frac{{2*K*P\left( {1 + K*P} \right)}}{{{{\left( {1 + K*P} \right)}^2}}}\\ &=& 2*\left( {\frac{{K*P}}{{1 + K*P}}} \right)\end{eqnarray*}$$\end{document}

And so, the null expectation for binding occupancy is simple twice the observed binding to a single recognition site.

For yeast TF motifs, we used the recommended PWMs compiled by Spivak and Stormo in the ScerTF database(stormo.wustl.edu/ScerTF). The ScerTF recommended PWM cutoff scores were used to define the presence or absence of TF sites on DNA sequences. The binding motif of MAX, the human bHLH factor, was obtained from factorbook (v1.factorbook.org/mediawiki/index.php/MAX). The NDS sequences used for this study were taken from a study by Raveh-Sadka (29); the NDS1 and NDS2 sequences in this work correspond to the v1 and v37 sequences from that study, respectively.

Cbf1p PBM data was obtained from UniProbe database, and we used dataset UP00397 for calculating free energy changes. We searched for each motif variant in PBM data, all the sequences that contains the same motif variant are grouped together, and the average PBM score was used to reflect the binding affinity for that variant. MITOMI data was obtained from the study by Maerkl (30) and the K_d for each relevant variant reported in the original publication was used for the calculation directly.

The CCRA method is designed to measure both TF binding and gene expression in parallel for hundreds of uniquely barcoded synthetic promoter sequences. To perform CCRA, the TF of interest is C-terminally fused to a short protein tag, so that the TF directs insertion of Ty5 retrotransposons (or ‘calling cards’) (24,25) near its binding sites (Figure 1A and C). For each CCRA assay, TF-directed insertions into the designed promoter library are recovered from yeast cells and the insertion locations and promoter sequence identities are determined via second-generation sequencing (Figure 1C). Each plasmid molecule in a CCRA library has a ‘library barcode’ corresponding to a unique promoter sequence (Figure 1A), as well as a unique molecular identifier (UMI). The library barcode allows each transposon calling card to be assigned to the correct synthetic promoter sequence, and the UMI enables us to determine when multiple transposition have inserted into the same location in distinct copies of the same synthetic promoter sequence. By determining the number of independent transpositions inserted into each synthetic promoter and then normalizing by the promoter's abundance in the library, we generate a normalized binding score (NBS), which is a quantitative measure of TF binding (Figure 1C).

Because the CCRA library is cloned upstream of a yellow fluorescence protein (YFP) reporter gene, it is also possible to measure the transcriptional output of each synthetic promoter in the library using Sort-Seq (17,18) (Figure 1B). To do so, the CCRA library is sorted by flow cytometry into subpopulations according to the ratio of YFP fluorescence to mCherry fluorescence. The mCherry gene is regulated by a constitutive promoter, allowing for normalization of the YFP signal to account for variation due to plasmid copy number, cell size, and other sources of extrinsic expression noise. Next, the sorted subpopulations of yeast cells are sequenced to quantify the abundance of each barcoded sequence in each subpopulation. Relative expression is then calculated by the proportion of each sequence in every binned library as per the standard Sort-Seq protocol (17,18). By combining aspects of both Calling Cards assay and Sort-Seq, CCRA allows us to quantitatively measure the binding of a TF to a library of regulatory sequences, and simultaneously measure the effect of that binding on gene expression.

To determine if CCRA can accurately and reproducibly measure TF binding in parallel, we first analyzed the binding of Cbf1p, a well-studied bHLH protein whose motif is strongly predictive of its in vivo binding pattern (23). To evaluate the sensitivity of the method for the detection of TF binding at weak sites, we created a library of 40 different sequences consisting of 10 synthetic promoters, each with four unique barcodes for replicates. Three of these sequences were taken from different endogenous yeast promoters previously shown to be bound by Cbf1p at a single recognition site (23). We also designed two synthetic promoters with nucleosome disfavouring sequences (29) that flanked a single Cbf1p consensus motif. As negative controls, we included five matched promoters with mutated Cbf1p binding sites. The binding of Cbf1p to a representative promoter, OYE3/DAP1, and its matched control is shown in Figure 2A. Each symbol on the graph represents an independent calling card insertion. Cbf1p-directed transpositions appear to fit a Gaussian distribution centred at Cbf1p motif. Interestingly, the region directly over the motif contains few insertions, likely due to Cbf1p's footprint as binds to its recognition sequence. The wild-type OYE3/DAP1 promoter is bound tightly by Cbf1p (70.1 NBS), but when the Cbf1p binding site is mutated, binding is greatly reduced (7.7 NBS, Figure 2A bottom panel). Cbf1p's binding to all five pairs of promoters is summarized in Figure 2B. In all instances, Cbf1p's binding was significantly stronger at promoters with intact Cbf1p sites than at the mutated promoters, demonstrating that the CCRA method can reliably detect TF binding even at relatively weak sites containing single motifs. It is interesting to note that although Cbf1p binding was significant at all five promoters with intact Cbf1p motifs, the binding was significantly stronger at the two promoters in which the Cbf1p binding sites were flanked by nucleosome disfavouring sequences.

We next investigated the dynamic range of the CCRA assay. Since Cbf1p binding at regulatory elements is known to strongly depend on the number of Cbf1p sites present (23), we designed 183 synthetic promoters containing 0 to 6 sites and measured the binding of Cbf1p to this library. We observed a strong non-linear relationship between the normalized binding score (NBS), and the number of sites present in a given promoter (Figure 2C). Importantly, we were able to measure Cbf1p binding across 3 orders of magnitude. These data demonstrate that CCRA technology can accurately measure TF binding across a large range of binding strengths.

Because the oligonucleotides used to create the synthetic promoters for CCRA are typically 170bp in length, we next sought to determine if TFs still bind in vivo with the same specificity as they do in their native genomic context. Therefore, we designed a 344-element library of genomic promoters derived from endogenous Gcn4p and Gal4p target promoters and used CCRA to measure the binding of these two TFs. We found that Gcn4p directed transpositions almost exclusively to synthetic promoters derived from Gcn4p targets whereas Gal4p directed transpositions to Gal4p targets (Figure 2D), with each TF showing little non-specific binding to the other TF’s set of target sequences. These results indicate that truncated genomic sequences in a plasmid-based system still retain their specificities and are not aberrantly bound by other TFs.

Having established that the CCRA assay measures TF binding with high sensitivity and specificity, we next sought to benchmark the method's reproducibility. To do so, we performed replicate CCRA experiments using a 531-element synthetic promoter library and found that the NBS measured for each library member was highly reproducible (Pearson r = 0.92, P-value = 4.23e–216, Spearman r = 0.63, P-value = 3.21e–59. Figure 2E).

We next sought to establish that the CCRA method could accurately and reproducibly measure expression of the YFP reporter driven by a synthetic promoter library. To determine accuracy, we performed Sort-Seq to measure reporter expression for each member of a library containing sequences derived from Gcn4p and Gal4p promoters. We then cloned 24 of these library members and individually measured their expression levels by flow cytometry. We observed excellent agreement between the two measurements; the Pearson correlation coefficient was 0.95 (Pearson P-value = 6.59e–13, Spearman r = 0.96 and P-value = 7.91e–14), indicating that CCRA methodology accurately measures promoter activities from a library of synthetic sequences (Figure 2F). To further investigate the accuracy of the method using a functional approach, we evaluated reporter expression as a function of the number of TF recognition sites for Gcn4p in an amino acid starvation growth condition and for Gal4p in galactose (31–35). For both factors, reporter expression increased with the number of motifs, as expected from the known mechanism of action for these TFs (Supplementary Figure S2 ). Finally, we also showed that expression measurements are highly reproducible between two biological replicates (Pearson r = 0.97 and P-value = 1.51e–228, Spearman r = 0.95 and P-value = 1.29e–177 Figure 2G).

The CCRA assay requires that the TF of interest be fused to a fragment of the Sir4p protein. This can be achieved by tagging the TF at its endogenous locus or by expressing the fusion from a plasmid, which is more convenient for many experiments. To investigate whether TF fusions expressed from plasmids binds to CCRA libraries in a similar manner as TF fusions expressed at their endogenous loci, we measured the binding for each using the same 531 synthetic promoter library and observed a high concordance (r = 0.84, Supplementary Figure S1). We also confirmed that transcription factors tagged with the Sir4p fragment do not influence Sort-Seq expression measurements as they are highly correlated with measurements made using untagged proteins (r = 0.94 for Gal4p, r = 0.99 for Gcn4p, Supplemental Figure S3A and B). Tagging TFs with Sir4p also does not appear to affect their functions (Supplemental Figure S3C–F). Taken together, these results demonstrate that the CCRA method accurately and reproducibly measures the TF binding and expression consequences to a library of synthetic promoters

Quantitative measurement of TF binding affinities to different DNA sequences is critical for understanding how TFs function in vivo. Because several studies have shown that minute variation in binding site affinity can specify alternative transcriptional or functional programs (36,37), it is important to be able to determine not only a TF’s consensus binding sequence, but also its binding energy landscape (i.e. the TF’s affinity for alternative binding sites). There are several methods that measure binding energy landscapes in vitro, such as MITOMI, PBM, Spec-seq, HT-SELEX, Bind-n-Seq, SPR, CSI and EMSA (30,38–49), and these have proven invaluable for understanding TF-DNA interactions. However, there is currently no method to accurately discriminate the small changes in free energy needed to generate binding energy landscapes in vivo. Such landscapes may differ from those measured in vitro due to the effects of nucleosomes and other chromatin-associated proteins on DNA shape and binding site accessibility. Therefore, we sought to determine whether CCRA could measure binding energy landscapes in vivo. We measured the binding of two basic helix loop helix (bHLH) factors, Cbf1p and MAX, to their consensus motifs and all sequences that differ by 1 bp from the consensus (Figure 3A). The TF binding sites were flanked by two intrinsic nucleosome disfavouring sequences to facilitate comparison to the in vitro binding landscapes previously determined (29). In order to accurately measure small changes in TF affinity, we used an expectation maximization algorithm to distinguish TF-directed transpositions from background insertions by assuming that TF-directed transpositions follow a Gaussian distribution centred at the consensus motif whereas non-specific transpositions follows a uniform distribution across the full synthetic promoter (see Materials and Methods). Cbf1p and MAX occupancies at their consensus binding sites and at all possible one base substitution are shown in Figure 3B and E respectively. As expected, both factors bound most strongly to their consensus sites. The changes in occupancies at non-consensus sites were strongly dependent on the position of the alteration and the identity of the substituted nucleotide. Some positions are crucial, such as the first position of core E-box motif, in the sense that any alternation resulted in completely abolished binding, whereas some positions such as flanking bases next to the core motif are more flexible when changed into other nucleotides. In general, Cbf1p binding appeared to be less tolerant to substitutions in its consensus motif than MAX, in agreement with previous in vitro measurements (30). We calculated the change of binding energy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta G$\end{document}) from consensus site to the alternative site as follows (see Methods for a detailed derivation):\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$$\begin{eqnarray*}\Delta \Delta G &=& \Delta G({\rm Sconsensus}) - \Delta G({\rm Smutant}) \\ &=& - {{RT{\rm ln}}}\left( {\frac{{{\rm Occ}({\rm Smutant})}}{{{\rm Occ}({\rm Sconsensus})}}} \right)\end{eqnarray*}$$\end{document}

To determine whether the measurements performed by CCRA are concordant with the binding energy landscapes of Cbf1p and MAX as measured by well-established in vitro methods, we compared our results to MITOMI and PBM (Figure 3C, D, F and G). Both methods generated energy landscapes that were highly correlated to our CCRA measurements (For Cbf1p the correlation between CCRA and MITOMI: Pearson r of 0.75 and P-value = 1.90 e–5, Spearman r of 0.70 and P-value = 9.40e–5; the correlation between CCRA and PBM: Pearson r of 0.72 and P-value = 5.29e–5, Spearman r of 0.73 and P-value = 3.97e–5. For MAX the correlation between CCRA and MITOMI: Pearson r of 0.72 and P-value = 1.42e–4, Spearman r of 0.74 and P-value = 9.06e–5; the correlation between CCRA and PBM: Pearson r of 0.80 and P-value = 7.24e–6, Spearman r of 0.77 and P-value = 2.36e–5). Since the correlations between the measurements made by the two in vitro methods are similar in magnitude (for Cbf1p, the correlation between MITOMI and PBM: Pearson r = 0.73 and P-value of 3.35e–5, Spearman r = 0.78 and P-value = 4.57e–6. For MAX, the correlation between MITOMI and PBM: Pearson r of 0.79 and P-value = 1.28e–5, Spearman r = 0.79 and P-value = 1.25e–5 Supplementary Figure S4), these results demonstrate CCRA measures binding energy landscapes in vivo with an accuracy comparable to in vitro methods.

The reported binding constant (K) for Cbf1p and MAX is (6.2 ± 1.4) × 10⁷ M⁻¹ at 20 °C (K_d = 1.6 nM) and (7.8 ± 2.6) × 10⁶ M⁻¹ (K_d = 130 nM) respectively (50,51), and therefore the binding energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta G$\end{document} for Cbf1p is about –45 kJ/mol (−18 K_BT) and –39 kJ/mol (–16 K_BT) for MAX. Given the largest \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta G$\end{document} calculated from the consensus to the mutant motif, Cbf1p loses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\frac{1}{4}$\end{document} of its binding energy with one nucleotide difference (e.g. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta G$\end{document} is 9.6 kJ/mol from GTCACGTG to GTCACGTA) and therefore the K_d on the mutated motif GTCACGTA becomes 71 nM, a 40-fold increase relative to the consensus motif. MAX loses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\frac{1}{{12}}$\end{document} of its binding energy with one nucleotide difference in vivo (e.g. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta G$\end{document} is 3.4 kJ/mol from CACGTG to CACTTG), and therefore the K_d on the mutant motif is 500 nM.

Understanding the mechanisms by which TFs select their targets in vivo will likely require more than just a characterization of their cognate DNA binding preferences, since it has been shown that many TFs achieve binding specificity through cooperative interactions with other DNA-binding proteins (10–15). Investigations into the cooperative interactions that occur between TFs are usually performed in vitro, under conditions that may not reflect the actual cellular environment (e.g. the lack of histones). In vivo investigations, which are less common, typically involve genome editing followed by quantitative binding measurement in vivo, which is experimentally challenging and time consuming (23,52,53). Given that CCRA is able to measure small changes in the free energy of TF binding, we sought to extend this approach to analyze TF–TF cooperativity. We focused on a pair of paralogous bHLH proteins, Cbf1p and Tye7p, both of which recognize the E-box motif CACGTG in vitro but bind to two distinct sets of target genes through different types of cooperative interactions.

We first set out to investigate Cbf1p, which has been shown to bind with homotypic cooperativity when two or more sites are present (23). This cooperativity was demonstrated by analyzing Cbf1p binding at mutated versions of the IDH1_NCE103 divergent promoter, which normally contains three Cbf1p binding sites. This study showed that Cbf1p occupancy at the wild-type promoter was much stronger than the sum of the binding occupancies at three mutated promoters, each containing only a single Cbf1p binding site, demonstrating that Cbf1p binding is not additive but instead cooperative at this locus. However, in this study, Cbf1p's cooperativity was investigated at only a single promoter, so it is unclear to what extent this result can be generalized. We therefore sought to use CCRA to determine if this phenomenon occurs at other loci. We selected seven promoters with two or three Cbf1p sites including IDH1_NCE103_pr and designed a CCRA library in which these promoter sequences contained either zero, one or two mutated Cbf1p sites. If Cbf1p binds cooperatively at these loci, we expect that, for each series of synthetic promoters, the sum of the binding scores from sequences with a single Cbf1p site will be significantly less than the binding at the ‘wild type’ promoter sequence with multiple Cbf1p sites. In all seven cases, we found that Cbf1p binding at the wild type promoter was significantly higher than would be expected under an additive binding model (Figure 4A), suggesting that Cbf1p binds cooperatively at all target promoters that contain multiple recognition sites.

We next sought to characterize the relationship between the strength of Cbf1p cooperative binding and the distance between binding sites. Because the DNA double helix is thought to be rigid over length scale less than ∼140 bp due to vertical base-stacking interactions and intra-helix phosphate charge repulsion (54,55), one might expect that Cbf1p dimers would be unable to bind cooperatively at promoters with two recognition sites in close proximity. However, Cbf1p has been shown to sharply bend DNA upon binding (56–58), and, furthermore, DNA is clearly malleable to some proteins, as it is tightly wrapped around nucleosomes and can be twisted and untwisted during replication and transcription (59–61). To investigate the relationship between Cbf1p cooperativity and the distance between recognition sites, we designed synthetic promoters where we varied the distance between two Cbf1p consensus motifs from 9 to 41 bp with 2 bp intervals. We used CCRA to measure Cbf1p binding on these synthetic sequences, and plotted binding occupancy as a function of the distance between two sites. We found that the strength of Cbf1p binding at these synthetic promoters varied periodically with the distance between the binding sites (Figure 4B). We observed strong binding at the shortest distance of 11 bp, and we observed additional peaks at 22, 32 and 41 bp apart. These distances are all shorter than the persistence length of DNA, and at the longest distance investigated, 41 bp, the binding sites are separated by >65 Å, so it seems unlikely that the interaction between Cbf1p dimers could be explained by protein domain flexibility. Therefore, these results suggest that Cbf1p's ability to bend DNA allows the two dimers to interact with one another. We next hypothesized that the observed periodicity could be explained by the fact that Cbf1p makes its base pair contacts in the major groove of DNA so that at some motif distances, contact between Cbf1p dimers would require the rotation of the major groove around the axis of the double helix, incurring an energetic penalty. To test this, we fitted the binding to a cosine function. The calculated period was 10.65 bp, almost exactly the number of base pairs required for DNA to make one complete helical turn about its axis. We evaluated the fit of this model using analysis of variance (ANOVA) and obtained a P-value 1.4e–6, indicating that the data follows the assumed model significantly better than expected by chance. This result suggested to us Cbf1p dimers that are not bound on the same side of the DNA helix must twist the DNA and incur an energetic cost. In contrast, two Cbf1p molecules on the same face of the helix are able to achieve the optimal cooperative binding efficiency. We next sought to compute the free energy cost associated with twisting the DNA double helix. Since we observed a 3.8-fold difference between the highest and the lowest occupancy, we calculated that the free energy lost due to twisting is 3.40 kJ/mol (1.37 K_BT). Compared to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} }{}$\Delta \Delta Gs$\end{document} calculated for the consensus to mutant motif from the previous section, the energic cost of DNA twisting is comparably to a mild nucleotide change in the E-box motif (e.g. from GTCACGTG to GTCTCGTG). Interestingly, over the distance range examined in this experiment, the amplitude of the periodic function did not change appreciably, suggesting that, in contrast to twisting, Cbf1p bends DNA efficiently, with little energetic cost.

We next asked if the phase of Cbf1p binding sites influenced the binding of this transcription factor at native genomic loci. We took published genome-wide Cbf1p Calling Cards data (23) and grouped all intergenic regions with two Cbf1p binding sites within 100 bp according to the relative phase of the two sites. We found that promoters containing two Cbf1p binding sites separated by a multiple of 10.5 bp (i.e. with major grooves on the same side of the DNA helix) were bound significantly more tightly by Cbf1p than promoters with binding sites whose major grooves were on opposite sides of the DNA helix (Figure 4C, P = 0.007). This result demonstrates that the periodicity in cooperative binding that we observed in our CCRA experiments also influences Cbf1p binding in the yeast genome.

Unlike Cbf1p, many of the promoters bound by Tye7p do not encode an E-box, this factor's preferred binding motif (23). It has previously been shown that Tye7p binds cooperatively with the Gcr1p/Gcr2p/Rap1p complex and that by taking into account the DNA binding preferences of these proteins, the in vivo binding of Tye7p can be more accurately predicted (23). However, the biophysical principles that govern the binding of this complex are still unclear. For example, the binding of this complex does not appear to follow either of the two most well-studied models for TF binding, the Enhancesome model or the Billboard model (62,63), because these models both posit a one-to-one correspondence between the binding of a TF and the presence of its recognition site. Instead, Tye7p binding appears to be consistent with the recently described TF collective model, in which a group of TFs bind together, but the motif positioning and composition at target sites is flexible (26,27). However, the TF collective model is ambiguous with regard to the mechanistic details of binding, so important questions about the function of the Tye7p/Rap1p/Gcr1p/Gcr2p collective remain.

We first assessed the predictive power of the collective model by attempting to reprogram yeast promoters that normally bind Cbf1p, a Tye7p paralog, into promoters that bind Tye7p. To do so, we took two promoters, OYE3_DAP1_pr and RPL1_RHO3_pr, that are normally bound by Cbf1p, and removed their E-boxes (i.e. Cbf1p/Tye7p binding sites), and added Gcr1/2p and Rap1p sites with a design based on the TDH3 promoter, which is bound by Tye7p. We then assessed the binding of Tye7p to these reprogrammed promoters using CCRA. Both showed significant decreases in Cbf1p binding (6.1-fold and 2.4-fold respectively) and significant increases in Tye7p (3.3-fold and 2.4-fold respectively) (Figure 5A). We also observed an increase in Gcr1p binding at these reprogrammed promoters. Since neither of these reprogrammed promoters contain a consensus Tye7p binding site, we conclude that Tye7p binding is consistent with the collective model and that this TF can be recruited to promoters via cooperative interactions with Gcr1/2p and Rap1p.

Next, we wanted to better understand the molecular logic by which this collective binds. While Tye7p clearly does not require its motif to be present at a regulatory target, is this true for other members of the collective? When more than one binding site is present for a single TF, do the additional sites contribute to complex stability, or is one site sufficient and the others redundant? How is transcriptional output correlated with binding of each TF member? To answer these questions, we took a Tye7p bound promoter, BMH1_pr, which contains one Tye7p site, three Gcr1/2p sites and one Rap1p site, made every possible combination of mutated sites, and measured Tye7p binding using CCRA. Since Tye7 does not require its recognition sequence for binding, we first wanted to know if its motif made any energetic contribution to stabilize this factor. We divided the mutated sequences into two categories, those with and without a Tye7p motif. Sequences without a recognition site were still significantly bound by Tye7p (Figure 5B, right group), consistent with previous observations, but Tye7p binding at the wild-type BMH1_pr is reduced by 45% when the Tye7p recognition site is mutated (P-value = 0.012). Furthermore, when the 16 pairs of BMH1_pr mutants are compared across groups, we observe a significant reduction in Tye7p when the recognition motif is mutated (P-value = 0.010). These results demonstrate that while the Tye7p motif is not required for Tye7p binding, it makes an energetic contribution when present. Notably, the positional distributions of Tye7p insertions across the BMH1_pr were essentially unaffected by the presence or absence of its cognate motif (Supplementary Figure S6), suggesting that the recruitment of Tye7p may be largely mediated by Gcr1/2p and Rap1p, even though the presence of a Tye7p binding site clearly makes an energetic contribution. Consistent with this hypothesis, we found that Tye7p binding is strongly dependent on Gcr1/2p and Rap1p sites (Figure 5B). In general, we observed a gradual decrease in binding as more collective sites are mutated, and we did not observe large decrease in binding (>2 fold) upon the removal of any one site, suggesting that no single binding site is necessary for Tye7p binding at this promoter, but instead that all sites contribute to the binding affinity of this TF. Based on this observation, we reasoned that Tye7p binding might be predicted by the total free energy from all sites combined on a promoter. Therefore, we performed a regression analysis to understand how well the total sites information explains Tye7p binding (Figure 5C). Given that PWM scores reflect the binding energy of TF to specific DNA sequences, we used the sum of PWM scores for all sites present on the promoters for the analysis and we found that the combined sites information correlates well with Tye7p binding (Pearson r = 0.69 and P-value = 1.07e–5, Spearman r = 0.63 and P-value = 1.04e–4).

We then measured Gcr1p and Gcr2p occupancy on this promoter library. As before, we divided the mutated promoters into two categories based on whether they contained a Gcr1/2p motif. In contrast to what was observed for Tye7p, we found that neither Gcr1p nor Gcr2p was able to bind at any promoters without their shared recognition site (Figures 5D & Supplementary Figure S7A), suggesting that these factors bind independently from the rest of the collective. To confirm this, we regressed Gcr1p and Gcr1p binding against the free energy of binding of Gcr1/2p or the full collective. We found that only Gcr1p/2p sites are required to explain Gcr1p and Gcr2p binding and that incorporating information from the other TF in the collective weakens the predictive power (Figure 5E & Supplementary Figure S7D for Gcr1p and Supplementary Figure S7B & CS for Gcr2p). Thus, the binding of the Gcr1/2p complex appears to be solely dependent on the presence and the number of Gcr1/2p sites. Furthermore, Gcr1/2p binding appears to saturate at two sites. Our Gcr2p binding measurements were more variable and weaker than our Gcr1p measurement, especially at sequences with only one Gcr1/2p motif, which might be due to the fact that Gcr2p is known to bind DNA indirectly through Gcr1p and depends on Gcr1p to function (64,65).

We next sought to investigate the relationship between the binding of the Tye7p collective and its transcriptional output. To do so, we performed Sort-Seq to measure the reporter gene expression from this library. We regressed reporter gene expression against the sum of the free energies of the binding sites (Supplementary Figure S7E). We observed a good correlation, and we found that expression level correlated with the combined TF occupancy (Figure 5F, Pearson r = 0.70 and P-value = 8.34e–6, Spearman r = 0.67 and P-value = 2.95e–9), suggesting that transcriptional output is determined by the whole complex. Similar analysis was done for TDH3 promoter containing two Gcr1/2p sites and two Rap1p sites but no Tye7p site, and again the combined Tye7p, Gcr1p and Gcr2p occupancy correlated well with the expression (Supplementary Figure S7F & GS).

Rap1p binding was not measured in this study due to its inability to be tagged by Sir4p. However, Rap1p has been shown to interact with Gcr1p and Gcr2p as an activating complex (66,67). With expression we measured on both BMH1 and TDH3 promoters, we compared sequence pairs that are with and without Rap1p site (Supplementary Figure S7H). We performed a paired T-test on these sequence in terms of expression, and the P-value is 0.018, indicating Rap1p motif is contributing the genetic regulation.

Taken together, our experiments suggest that Tye7p is recruited to promoters by Gcr1p/Gcr2p/Rap1p complex and that Tye7p binding often occurs in the absence of its recognition site. However, it appears that Tye7p binding is stabilized by the presence of its motif. In contrast, the Gcr1/2p recognition site is necessary and sufficient for the binding of these proteins, suggesting a hierarchy in which these factors can recruit Tye7p but not vice versa (Figure 5G). The transcriptional output at promoters bound by this complex correlates with the combined occupancy of all TFs, suggesting that each TF in the collective aides in the recruitment of the RNA Polymerase II holoenzyme.