Dynamics of Transcription Factor Binding Site Evolution

Abstract

Evolution of gene regulation is crucial for our understanding of the phenotypic differences between species, populations and individuals. Sequence-specific binding of transcription factors to the regulatory regions on the DNA is a key regulatory mechanism that determines gene expression and hence heritable phenotypic variation. We use a biophysical model for directional selection on gene expression to estimate the rates of gain and loss of transcription factor binding sites (TFBS) in finite populations under both point and insertion/deletion mutations. Our results show that these rates are typically slow for a single TFBS in an isolated DNA region, unless the selection is extremely strong. These rates decrease drastically with increasing TFBS length or increasingly specific protein-DNA interactions, making the evolution of sites longer than ∼ 10 bp unlikely on typical eukaryotic speciation timescales. Similarly, evolution converges to the stationary distribution of binding sequences very slowly, making the equilibrium assumption questionable. The availability of longer regulatory sequences in which multiple binding sites can evolve simultaneously, the presence of "pre-sites" or partially decayed old sites in the initial sequence, and biophysical cooperativity between transcription factors, can all facilitate gain of TFBS and reconcile theoretical calculations with timescales inferred from comparative genomics.

Fig 1. Biophysics of transcription regulation.

A) TFs bind to regulatory DNA regions (promoters and enhancers) in a sequence-specific manner to regulate transcriptional gene expression (mRNA production) level via different mechanisms, such as recruiting RNA polymerase (RNA-pol). B) A schematic of two types of mutational processes that we model: point mutations (left) and indel mutations (right). C) The mismatch binding model results in redundancy of genotype classes, with a binomial distribution (red) of genotypes in each mismatch class (some examples of degenerate sequences shown) D) The mapping from the TFBS regulatory sequence to gene expression level is determined by the thermodynamic occupancy (binding probability) of the binding site. If each of the k mismatches from the consensus sequence decreases the binding energy by ϵ, the occupancy of the binding site is π _TD(k) = (1 + e ^{β(ϵk−μ))−1}, where μ is the chemical potential (related to free TF concentration). A typical occupancy curve is shown in black (ϵ = 2 k _B T and μ = 4 k _B T); the gray curves show the effect of perturbation to these parameters (ϵ = 1 k _B T, ϵ = 3 k _B T and μ = 6 k _B T); the orange curve illustrates the case of two cooperatively binding TFs (k _c = 0 and E _c = −3 k _B T, see text for details). We pick two thresholds, shown in dashed lines, to define discrete binding classes: strong 𝓢 (π _TD > 2/3) and weak 𝓦 (π _TD < 1/3).

Fig 2. Single TF binding site gain rates at an isolated DNA region.

A) The dependence of the gain rate, 1/⟨t⟩_{𝓢 ← k} shown in units of point mutation rate, from sequences in different initial mismatch classes k (blue: k = 2, red: k = 5), as a function of selection strength. Results with point mutations only (θ = 0) are shown by dashed line; with admixture of indel mutations (θ = 0.15) by a solid line. For strong selection, Ns ≫ n log(2)/2, the rates scale with Ns, which is captured well by the “shortest path” approximation (black dashed lines in the main figure) of Eq (24). The biophysical parameters are: site length n = 7 bp; binding specificity ϵ = 2 k _B T; chemical potential μ = 4 k _B T. Points correspond to Wright-Fisher simulations with Nu = 0.01 where error bars cover ±2 SEM (standard error of mean). Inset shows the behavior of the gain rates as a function of the initial mismatch class k for Ns = 0 and Ns = 100. B, C) Gain rates from redundancy rich classes (k ∼ 3n/4, typical of evolution from random “virgin” sequence) under strong selection, without (B) and with (C) indel mutations supplementing the point mutations. Red crosshairs denote the cases depicted in panel A. Contour lines show constant gain rates in units of Ns u as a function of biophysical parameters n and ϵ. Wiggles in the contour lines are not a numerical artefact but a consequence of discrete mismatch classes.

Fig 3. Convergence to the stationary distribution of TFBS sequences.

A) Evolutionary dynamics of the mismatch classes distribution ψ(k) for an isolated TFBS under point and indel mutations (θ = 0.15), directional selection for stronger binding, and genetic drift is shown for initially well (k = 0, blue) and badly (k = 5, red) adapted populations. At left, no selection (Ns = 0); at right, strong selection (Ns = 100). Different curves show the distribution of genotype classes at different time points (t = 0u ⁻¹, 0.05u ⁻¹, 0.1u ⁻¹ as decreasing opacity); stationary distribution is shown in green. Insets show the time evolution to convergence for initially well (k = 0, blue) and badly (k = 5, red) adapted populations, measured by the Kullback-Leibler divergence D _KL[ψ(t) ∣∣ ψ(t = ∞)]. The biophysical parameters are: n = 7 bp, ϵ = 2 k _B T, μ = 4 k _B T. B) Rate of convergence to the stationary distribution for different ϵ and n values under strong selection (Ns ≫ n log(2)/2; here specifically Ns = 100) and for θ = 0.15. Crosshairs represent the parameters used in a).

Fig 4. TF binding site evolution in a longer sequence of L = 30 base pairs.

The expected number of newly evolved TF binding sites with length n = 7 bp, under strong directional selection (Ns = 100) and both point and indel mutations (θ = 0.15). Time is measured in inverse mutation rates; the number of newly evolved sites is scaled to the selection strength and the sequence length. 1000 replicate simulations were performed with different initial sequences. Average number of sites shown by a solid black line; the gray band shows ±2 SEM (standard error of the mean) envelope. Dashed curves are analytical predictions based on single TFBS gain rates at an isolated DNA region, given by Eqs (27), (28) and (29). Biophysical parameters used: ϵ = 2 k _B T, μ = 4 k _B T. Insets: Expected number of newly evolved sites from a random sequence of length L at t = 0.001u ⁻¹ (left) and t = 0.1u ⁻¹ (right) for different binding length and specificity values, computed using the analytical predictions. Crosshairs denote the values used in the main panel.

Fig 5. Ancient sites and cooperativity can accelerate the emergence of TF binding sites in longer regulatory sequences.

A) The expected number of newly evolved TFBS in the presence (red and brown) or absence (black) of an ancient site, for binding site length n = 10 bp, and specificity, ϵ = 3 k _B T. In this example, the ancient site was a consensus site (k = 0) or two mismatches away from it (k = 2) that evolved under neutrality for t′ = 0.1/u prior to starting this simulation. Dashed lines show the predictions of a simple analytical model, Eq (30). The inset shows how the number of newly evolved TFBS at t = 0.001/u scales with the mismatch of the ancient site k (plot markers: simulation means; error bars: two standard errors of the mean; dashed curve: prediction). B) The expected number of newly evolved TFBS without (black) and with cooperative interactions (for different cooperativity strengths, magenta: E _c = −2 k _B T, yellow: E _c = −3 k _B T, cyan: E _c = −4 k _B T, see Eq (11) in Methods and text) for binding site length n = 7 bp, and specificity, ϵ = 2 k _B T. Both panels use μ = 4 k _B T, strong selection (Ns = 100) and a combination of point and indel mutations (θ = 0.15), acting on a regulatory sequence of length L = 30 bp. Thick solid lines show an average over 1000 simulation replicates, shading denotes ±2 SEM.