De Novo Identification and Biophysical Characterization of Transcription Factor Binding Sites with Microfluidic Affinity Analysis

Recent evidence suggests that knowledge of both strongly-and weakly-bound sequences and their interaction affinities is required for an accurate understanding of transcriptional regulation. Weak-affinity sites are evolutionarily conserved, make significant contributions to overall transcription^1,2, and may allow closely related TFs to mediate different transcriptional responses³. In addition, quantitative models require both strongly-and weakly-bound sequences and their binding affinities to recapitulate transcriptional responses^4-7.

Unfortunately, quantitative data detailing TF binding are often lacking, even for model organisms. In vivo immunoprecipitation-based methods (e.g. ChIP-chip⁸ and ChIP-SEQ⁹ provide genome-wide information about promoter occupancy. However, these techniques require knowledge of physiological states under which TFs are bound to promoters, cannot distinguish whether a TF contacts DNA directly or is tethered via another DNA-binding protein, and do not measure affinities.

In vitro methods complement in vivo data by measuring binding affinities, distinguishing whether TFs directly bind DNA, and allowing manipulation of post-translational modifications and buffer conditions. Furthermore, in vitro methods can be used without knowledge of conditions under which TFs are active. However, current in vitro methods cannot simultaneously discover both high-and low-affinity target sequences and measure their affinities. Electromobility shift assays (EMSAs)¹⁰ DNAse footprinting¹¹ and surface plasmon resonance¹² require prior knowledge of potential binding sites, precluding motif discovery. Conversely, selection techniques (i.e. SELEX) and one-hybrid systems¹³ discover motifs from a large sequence space, but recover only the most strongly bound sequences, without affinity information. Protein binding microarrays (PBMs)^3,14-18 can discover both strongly-and weakly-bound sequences but cannot measure reactions at equilibrium, preventing affinity measurements. PBMs also suffer from reduced sensitivity: a recent study using PBMs to probe TF binding in S. cerevisiae failed to recover consensus motifs for 49 of 101 TFs with previous evidence of direct DNA binding¹⁵. Embedding immobilized DNA in hydrogels¹⁹ extends the PBM technique to allow affinity and kinetic measurements, but limits available DNA sequences to ∼ 100.

An alternative approach is Mechanically-Induced Trapping of Molecular Interactions (MITOMI), a technique that uses a microfluidic device to measure binding interactions at equilibrium, allowing construction of detailed maps of binding energy landscapes. The first-generation MITOMI device measured 640 parallel interactions and required TF-specific DNA libraries²⁰.

Here, we report a second-generation MITOMI device (MITOMI 2.0) capable of measuring 4,160 parallel interactions. Devices were fabricated in polydimethylsiloxane (PDMS) using multilayer soft lithography; each device had 4,160 unit cells and approximately 12,555 valves to control fluid flow (Fig. 1a). Each unit cell contained a DNA chamber and a protein chamber, controlled by micromechanical valves: a ‘neck’ valve, ‘sandwich’ valves, and a ‘button’ valve (Fig. 1a, Supplementary Fig. 1). Unit cells were programmed with particular DNA sequences by aligning and bonding the device with a non-covalently spotted DNA microarray containing a library of 1457 double-stranded Cy5-labeled oligonucleotides. To accommodate all 65,536 DNA 8-mers, each 70-bp oligonucleotide contained 45 overlapping, related 8-mer de Bruijn sequences²¹ (Fig. 1b). Each oligonucleotide sequence appeared in at least 2 unit cells.

To evaluate the performance of this technique, we measured DNA binding for 28 S. cerevisiae TFs from 10 different families (Supplementary Table 1). Of these, 26 TFs had prior evidence of direct, sequence-specific DNA binding and 2 TFs had no previously annotated literature motifs, despite multiple previous attempts^14,15,22.

All TF protein was produced by in vitro transcription/translation. PCR-generated linear expression templates were added directly to rabbit reticulocyte lysate off-chip in the presence of a small fraction of BODIPY-labeled lysine charged tRNA to produce BODIPY-labeled, His-tagged TFs (Fig. 1c, Supplementary Fig. 2). In each experiment, ∼ 50 μL of extract (∼ 100 ng of protein) was loaded into the device.

Following alignment to DNA microarrays, slide surfaces within the protein chamber were derivatized with anti-pentaHis antibodies beneath the button valve and passivated elsewhere (Fig. 1d). Introduction of His-tagged TFs into both chambers solubilized spotted DNA, allowing TFs and DNA to interact. TF-DNA complexes were captured on the surface beneath the button valve during a ∼ 1 hour incubation; rapid closure of the button valve trapped interactions at equilibrium concentrations prior to a final wash to remove unbound material before imaging²⁰.

BODIPY intensities under the button valve reflect the number of surface-bound protein molecules; Cy5 intensities under the button valve reflect the number of DNA molecules bound by surface-immobilized protein (Fig. 1d,e,f). Therefore, the ratio of Cy5 to BODIPY fluorescence is linearly proportional to the number of protein molecules with bound DNA, or protein fractional occupancy. Cy5 intensities within the DNA chamber reflect the amount of soluble DNA available for binding.

All 28 TFs showed oligonucleotide-specific variations in bound Cy5 intensities, demonstrating marked preferences for individual oligonucleotides (Fig. 2a, Supplementary Fig. 3). By contrast, the distribution of intensities for rabbit reticulocyte extract alone was well-fit by a Gaussian (reduced χ² = 1.0, p = 0.47), establishing that binding is due to expressed TFs and not components of the in vitro transcription/translation system (Fig. 2a).

Variations in fluid flow between channels can lead to differences in the number of protein molecules beneath each button valve. To account for these differences and generate a quantity proportional to fractional occupancy, Cy5 intensities were normalized by BODIPY intensities to yield a dimensionless intensity ratio (Cy5 intensity/BODIPY intensity) (Fig. 1e). Intensity ratios also showed strong preferences for individual oligonucleotide sequences, with no clear preference detected for rabbit reticulocyte lysate alone (Fig. 2b; Supplementary Fig. 4, Supplementary Table 2). Intensity ratios were well correlated both between measurements of the same 70-mer oligonucleotide at different locations within a given device (Fig. 2c; Supplementary Table 3) and between experiments (Supplementary Fig. 5).

Binding affinity can be described by a single-site binding model relating intensity ratio (r) to DNA concentration ([D]); K_d, the DNA concentration at which measured intensities reach half their maximum value (r_max) provides a quantitative measure of binding affinity (Eqn. 1):

At low DNA concentrations, measured intensity ratios are approximately inversely proportional to K_d. Calibrated measurements of DNA chamber intensities in our experiments establish that soluble DNA concentrations are indeed low (150 ± 25 nM, mean ± s.e.m.) (Supplementary Fig. 6), suggesting it might be possible to accurately estimate interaction affinities from intensity ratios measured at a single, low DNA concentration.

To test this hypothesis, we first measured concentration-dependent binding for 4 TFs from 2 different families (Cbf1p, Cin5p, Pho4p, and Yap1p), each interacting with 10 oligonucleotides from the 8-mer DNA library. We then globally fit Eqn. 1 over all oligonucleotides at all concentrations to get accurate K_d measurements (Fig. 3a,b; Supplementary Fig. 7, Supplementary Fig. 8, Supplementary Fig. 9).

Next, we calculated K_d values for the exact same oligonucleotides from single-concentration measurements. The low DNA concentration used for these measurements prevented direct determination of r_max, a parameter that depends on quantities that vary between experiments (e.g. amount and intensity of BODIPY and Cy5 dyes incorporated during protein and DNA library production, respectively), and must be empirically determined. K_d values from concentration-dependent binding can be used to “calibrate” the appropriate r_max value (Supplementary Information). Single-concentration K_d values calculated using calibrated r_max values were in excellent agreement with those derived from concentration-dependent binding (r² = 0.90, p = 2.1 × 10^-19) (Fig. 3b). Furthermore, once calibrated, r_max values can be used to calculate K_d values for all oligonucleotides with signals above background, providing absolute affinities for all 1457 oligonucleotides with only a few additional measurements (Fig. 3c,d; Supplementary Fig. 10). The range of K_d values calculated here for Pho4p and Cbf1p agree with those measured in previous studies (∼ 10 nM to 10 μM)²⁰, validating our approach. Relative differences in binding affinities between oligonucleotides (the Gibbs free energy upon binding, ΔΔG) can also be calculated using these calibrated r_max values (Supplementary Fig. 11).

Even in the absence of additional information to calibrate r_max values, however, measured intensity ratios provide accurate information about binding affinity. To demonstrate this, we assumed an r_max value of 1 for all TFs and again compared measured and calculated K_d values. K_d measurements were well correlated (r² = 0.67, p = 1.8 × 10^-10), although precise values were somewhat offset (Supplementary Fig. 12a). ΔΔG describes relative affinity differences between oligonucleotides and is therefore less sensitive to these offsets, with stronger correlations (r² = 0.76, p = 8.0 × 10^-13) (Supplementary Fig. 12b).

Measured intensity ratios reflect interaction affinities between a given TF and a 70-bp oligonucleotide. Identifying TF target sites requires determination of the precise subsequences responsible for TF binding within each oligonucleotide. Traditionally, analysis of TF binding requires designation of sequences into bound and unbound populations, followed by a search for sequences overrepresented in the bound population, which ignores relative strengths of binding interactions, and can be sensitive to the precise threshold used to delineate populations. Here, we used a pipeline that incorporates all intensity information for all oligonucleotides to generate a position-specific affinity matrix (PSAM)²³ describing the change in binding affinity upon mutation of a specific position within a consensus sequence (Supplementary Fig. 13). Notably, PSAMs describe actual binding affinities for any combination of nucleotides and can used to calculate predicted affinities to arbitrary sequences.

First, we analyzed all measured intensity ratios using fREDUCE, an enumerative algorithm that searches for sequences whose occurrence within oligonucleotides correlates strongly with their measured signal²⁴. For all 28 proteins, fREDUCE returned sequences whose appearance within an oligonucleotide correlated strongly with measured intensity ratios (Fig. 5, Supplementary Table 6, Supplementary Fig. 14).

Next, the highest-correlated 7-and 8-bp fREDUCE sequences were converted to PSAMs using MatrixREDUCE²³, an algorithm that fits all measured intensity ratios with a statistical mechanical model assessing the effects of individual base pair substitutions on binding affinity. Investigations of MatrixREDUCE performance have recommended the use of initial seed sequences derived from enumerative analysis to ensure optimization of global minima²⁴ therefore, the fREDUCE sequences were used as seeds. MatrixREDUCE assumes that the free energy contributions of each position in the binding site are independent; although this is known to be false in some instances, we employ linear motifs here to compare our results with the largest possible set of previous literature.

To choose the single PSAM that best explains measured binding, we compared occupancies predicted by each PSAM for all oligonucleotides in the DNA library with measured intensity ratios (Supplementary Fig. 15). Predicted and measured values were well-correlated for almost all TFs (Supplementary Table 7). For all 26 TFs with described motifs, the final recovered motif was in agreement with those previously reported in the literature (Fig. 4) ^14,15,22. We also derive PSAMs for two TFs, Msn1p and Nrg2p, that were previously resistant to characterization, establishing significantly enhanced sensitivity over both ChIP-based and PBM techniques.

Two well-characterized basic helix-loop-helix (bHLH) proteins (Pho4p and Cbf1p) provide a test of the ability to detect both high-and low-affinity target sequences. Pho4p binds both high-affinity (5′-CACGTG-3′) and low-affinity (5′-CACGTT-3′) sites²⁵ Cbf1p binds to a degenerate ‘5-RTCACRTG-3’ motif^20,26. For both proteins, we recover the expected motif variants (Fig. 4, Supplementary Fig. 15).

Detailed analysis of differences between measured and calculated binding profiles can provide additional information about binding preferences. For example, oligonucleotides with high measured intensity ratios but low predicted occupancies, could indicate binding to additional motifs. In addition, this comparison allows investigation of whether free energy contributions at each position within the sequence are truly independent.

For most TFs, optimized PSAMs successfully described gross binding properties (e.g. Pho4p, Cin5p, Msn2p, and Sko1p; Supplementary Fig. 16), albeit with outliers at weak binding energies that may represent cooperative interactions between base pair substitutions. For a few transcription factors (Rpn4p, Cup9p, Cad1p, Matα2p, and Pdr3p), correlations between measured and predicted binding were much weaker (r² < 0.25). To determine if low correlations resulted from binding to additional target sequences, we used BioPROSPECTOR²⁷, MDScan²⁷, MEME²⁸, and WEEDER²⁹ to scan for overrepresented sequences within oligonucleotides with high measured intensity ratios (Z-score > 25 or 75) but low predicted occupancies (Z-score < 3).

For Rpn4p, although both PBM studies and our first analysis identified binding to a 5′-GCCACC-3′ motif, ChIP and expression data suggest a T-rich 5′ extension of this motif upstream of Rpn4p target genes. Strikingly, analysis of the 13 oligonucleotides with discordant measured and predicted binding returned this precise extension, establishing that unexpected binding data can yield biologically relevant results (Supplementary Fig. 17).

The Cup9p optimized PSAM also agreed with previous PBM¹⁵ results (Fig. 4); however, 14 sequences showed stronger-than-predicted binding (Supplementary Fig. 18). Analysis of these sequences yielded motifs similar to the optimized PSAM, but with an ‘ACGT’ core (Supplementary Fig. 18, grey box). To assess the affinity of Cup9p for this candidate alternate motif, we measured concentration-dependent binding of Cup9p to the primary motif, candidate secondary motif, and several related motifs (Supplementary Fig. 19a). A random 2 bp substitution abolished binding, but mutating these bases or the entire second half of the motif to the candidate secondary motif reduced affinity only ∼ 20-fold (Supplementary Fig. 19b), confirming weak-affinity binding. Interestingly, this motif is found only 29 times in the genome outside of coding regions, primarily at the boundary of subtelomeric repeats and upstream of genes regulated by iron depletion, metal toxicity, or oxidative stress (Supplementary Table 8). While the physiological role of these putative binding sites is unknown, these results demonstrate the ability of MITOMI 2.0 to detect weak but potentially biologically relevant TF binding sites.

For the remaining 3 TFs (Cad1p, Matα2p, and Pdr3p), low correlations between predicted and measured binding likely result from experimental variability and not binding to additional motifs. Correlations between technical replicates across the device were relatively low (Supplementary Table 3) due to either binding to a limited number of oligonucleotides (Cad1p, Supplementary Fig. 3) or large variations in protein coverage (for Matα2p and Pdr3p). Consistent with this, these TFs do not bind any oligonucleotides with stronger-than-expected affinity.

The data presented here demonstrate increased sensitivity over current state-of-the-art techniques, detecting sequence-specific binding for several proteins that have failed to yield results in multiple experiments (Cad1p, Msn1p, Nrg2p, Sko1p, Yap7p, and Pdr3p). Moreover, these data represent the most comprehensive investigation of biophysical binding affinities to date, including ΔΔG values for 28 TFs and K_d values for 4 TFs from 2 different families (Cbf1p, Cin5p, Pho4p, and Yap1p) binding to 1457 individual sequences. These data can be used to test basic assumptions underlying current models of TF-DNA specificity and more accurately model cooperativity between nucleotide binding sites (‘non-additivity’).

The DNA library used here is not organism-specific, making this technique useful for a wide range of organisms, including higher eukaryotes and pathogens. In addition, the programmable nature of MITOMI 2.0 allows subsequent detailed examination of unexpected binding phenomena or systematic mutational analysis of candidate motifs through direct observations of concentration-dependent binding. Although these experiments probed TF binding to double-stranded DNA, MITOMI 2.0 can be used with only minimal changes to investigate single-stranded DNA binding and RNA binding. When paired with advances in rapid whole genome sequencing, we anticipate that MITOMI 2.0 characterization of all recognizable TFs in a given proteome will allow transcriptional networks and regulons to be quickly identified and ultimately modeled.