A model of evolution with constant selective pressure for regulatory DNA sites

Background Molecular evolution is usually described assuming a neutral or weakly non-neutral substitution model. Recently, new data have become available on evolution of sequence regions under a selective pressure, e.g. transcription factor binding sites. To reconstruct the evolutionary history of such sequences, one needs evolutionary models that take into account a substantial constant selective pressure. Results We present a simple evolutionary model with a single preferred (consensus) nucleotide and the neutral substitution model adopted for all other nucleotides. This evolutionary model has a rate matrix in which all substitutions that do not involve the consensus nucleotide occur with the same rate. The model has two time scales for achieving a stationary distribution; in the general case only one of the two rate parameters can be evaluated from the stationary distribution. In the middle-time zone, a counterintuitive behavior was observed for some parameter values, with a probability of conservation for a non-consensus nucleotide greater than that for the consensus nucleotide. Such an effect can be observed only in the case of weak preference for the consensus nucleotide, when the probability to observe the consensus nucleotide in the stationary distribution is less than 1/2. If the substitution rate is represented as a product of mutation and fixation, only the fixation can be calculated from the stationary distribution. The exhibited conservation of non-consensus nucleotides does not take place if the elements of mutation matrix are identical, and can be related to the reduced mutation rate between the non-consensus nucleotides. This bias can have no effect on the stationary distribution of nucleotide frequencies calculated over the ensemble of multiple alignments, e.g. transcription factor binding sites upstream of different sets of co-regulated orthologous genes. Conclusion The derived model can be used as a null model when analyzing the evolution of orthologous transcription factor binding sites. In particular, our findings show that a nucleotide preferred at some position of a multiple alignment of binding sites for some transcription factor in the same genome is not necessarily the most conserved nucleotide in an alignment of orthologous sites from different species. However, this effect can take place only in the case of a mutation matrix whose elements are not identical.


Background
The controlled expression of genes is the main mechanism responsible for the life cycle and biodiversity [1]. Transcription, which is a crucial process defining the level of gene expression, is regulated by interaction of transcription factor proteins (TFs) with transcription factor binding sites (TFBSs) in a DNA molecule [2]. Thus, adaptable interaction between TFs and TFBSs is one of the main driving forces of biological evolution [3].
Both TFs and TFBSs are subject to mutations and selection affecting their interaction [4]. In this study we focus on mutations in TFBSs. New experimental [5][6][7] and computational [8] methods of TFBS identification produce an increasing amount of data about TFBS sequences, which creates a possibility to study evolutionary events in these regions.
Modelling evolution of regulatory sequences can be useful for understanding both the general mechanisms of gene expression control and the regulatory history of particular genes.
The evolution of regulatory regions has a complex pattern [3,[9][10][11] and is still unclear in many aspects. DNA segments can gain and lose the TFBS function [12], and that can bring new genes under regulation by a particular TF [13][14][15], or divert regulation from other genes [16,17]. One particular type of events is emergence of new or changed sites following displacement of a transcription factor by horizontal transfer [18]. Sometimes this leads to considerable changes in the regulon content [19,20] or even partial or complete rewiring of regulatory cascades [21][22][23][24][25][26], reviewed in [27].
Evolution of functional TFBS sequences is strongly nonneutral [11,28] and under a positive selection [29], which makes it difficult to calculate the rate of TFBS evolution. This rate varies between TFBS positions [30,31]. Moreover, co-evolution of TFBS and TF [16,20] can make the selective pressure vary in different lineages. Indeed, although in some cases the DNA motif bound by orthologous factors may be conserved at surprisingly large evolutionary distances [14,[32][33][34], in other cases not only the motifs themselves may be different [35,36], but even the symmetry of the motif (e.g. palindrome or direct repeat) may change [37,38].
The existing evolutionary models, which were successful in reconstruction of phylogenetic relations, can be applied to evolution of regulatory sequences only with a caution. Such models are historically related to the Jukes-Cantor [39] and Kimura [40] models of molecular evolution. Existing modifications of these models take into account various global characteristics like transition/transversion rate or local GC composition [41][42][43][44]. They are not applicable to the case of strong selection for a specific nucleotide at a particular position.
On the other hand, models developed specifically for the evolution of TFBS are needed to reconstruct the evolutionary origin of a particular TFBS and to evaluate the position-specific mutation rate and selective pressure.
Because of the position-specific variation in the rate of TFBS evolution [30], the rate matrix must also be position-specific. The data produced by mass experiments on TFBS identification or comparative genomic studies produce tens of TFBS for each TF (more exactly, a group of orthologous TFs). That might be sufficient to evaluate the evolutionary rate at each TFBS position.
Here we consider the simplest model of position-specific evolution with one preferred (consensus) nucleotide and three other (minor) nucleotides, the latter considered in a symmetric setting, without any selection or rate preferences [45]. Such a model can be deduced from physical requirements of the TF/TFBS interaction [46] and can explain the observed TFBS fuzziness.
We build a rate matrix, which enhances the model of [45]. We calculate the substitution probability for each finite time and show that the nucleotide conservation in phylogenetic lineages can be non-trivial for some parameter values. Particularly, a non-consensus nucleotide may appear more conserved than the consensus nucleotide, although the latter has a selective preference. This happens when the rate of mutations between non-consensus nucleotides is lower than the rate of mutation into the consensus, or there is selection against any mutation in non-consensus nucleotides.

Model
We start with definitions. We consider an alignment of several sequences; all positions in this alignment are assumed to be independent, and thus may be modelled independently. Consensus nucleotide (or simply consensus) is the most frequent nucleotide in an alignment column (position). Other nucleotides are called non-consensus. The frequency of the consensus nucleotide N c is a fraction of the number of consensus nucleotides in a particular position. Obviously, 1/4 <N c ≤ 1. The consensus is called weak, if 1/4 <N c < 1/2; the consensus is strong, if 1/2 ≤ N c < 1.
Consider the model of nucleotide substitutions given by a Markov process X(t) with four states {g 1 , g 2 , g 3 , g 4 }. Without loss of generality, assume that the state g 1 is the consensus state and the states g 2 , g 3 , g 4 are equiprobable non-consensus states. Suppose that the transition rate matrix A = (q ij ) is given by where q ij is the transition rate from the state g i to g j and α, β, and δ are positive unknown parameters.

Transitional probabilities
Let P ij (t) be the transitional probability from the state g i to the state g j for the time t. From the theory of Markov chains we have for h → 0 For brevity we denote by 'c' the subscript '1' corresponding to the consensus state g 1 and by 'n' and 'm' the subscripts 2, 3, and 4 corresponding to the non-consensus states g 2 , g 3 , g 4 . Thus we have five transitional probabilities P cc , P cn , P nc , P nn , and P nm , where n and m stand for two distinct non-consensus states. The formulas for P ij are derived from simple calculation: For simplicity, introduce a new time-scale, u = tβ, and denote λ = α/β, μ = δ/β. Then and Note that the parameter μ is related to the conservation within the set of non-consensus states {g 2 , g 3 , g 4 }. Simply, λ is a rate of transition from the consensus nucleotide and μ is a rate of transition between non-consensus states up to the time-scale parameter β.
The stationary distribution π = (π c , π n , π n , π n ) of the process X(t) is given by According to our definition of the consensus as the most frequent nucleotide, from π c = 1/(3λ + 1) > 1/4 it follows that 0 <λ < 1. The parameter λ is responsible for transitions from the consensus state to the non-consensus ones. At the same time, the transition from a non-consensus state to the consensus state occurs with the rate 1 (up to the time-scale parameter β). Intuitively, in the model with one selected consensus state the probability of transition to a non-consensus state should be smaller than the probability of transition to the consensus state; thus, λ should be less than 1. Indeed, if λ ≥ 1, then π c ≤ 1/4, π n ≥ 1/4 and we have three equiprobable consensus states.

Conservation and transition of nucleotides
The goal of this section is to compare the transitional probabilities between different states in our model. We will see that the relations we get have clear biological interpretation.

Simple cases
Let a unique consensus exist in our model, i.e. 0 <λ < 1. It is not difficult to obtain the following relations between the transitional probabilities, which hold for any λ ∈ (0, 1), u > 0.
3. P cc (u) > P nm (u) for any μ > 0; 4. P nn (u) > P cn (u) for any positive λ and μ; The interpretation of these results is straightforward. between the transitional probabilities are clearly shown in the figures. The first three inequalities show that the probability of conservation of the consensus state is always higher than the probability of transition between the Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 2/3, μ = 1/6 Figure 3 Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 2/ 3, μ = 1/6. In this case the probability of conservation of the consensus nucleotide is less than the probability of conservation of a non-consensus nucleotide, P cc (u) <P nn (u) for u ∈ (0, u**), since 3λ > 2μ + 1 and λ > μ. The stationary distribution of the consensus state is π c = 1/3. Thus, the consensus is weak. Since μ <λ, we have P cn (u) > P nm (u) for all u > 0. Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 1/6, μ = 1/2 Figure 1 Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 1/ 6, μ = 1/2. This figure describes the case P cc (u) > P nn (u) for all u > 0. The consensus is strong, since π c = 2/3. Since the rate of transition between non-consensus states μ is greater than the rate of transition from the consensus state λ, we have P cn (u) <P nm (u). The relation between the probability of conservation of a non-consensus nucleotide and the probability of transition from the non-consensus state to the consensus one is shown, P nn (u) > P nc (u) for u ∈ (0, u 0 ) and, vice versa, P nn (u) <P nc (u) for u > u 0 . This figure exemplifies all simple cases that hold for any λ ∈ (0, 1), μ > 0, u > 0 (see the corresponding paragraph in the main text). Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 2/3, μ = 3/4 Figure 2 Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 2/ 3, μ = 3/4. Here we have P nn (u) > P nc (u) for u ∈ (0, u 0 ) and, vice versa, P nn (u) <P nc (u) for u > u 0 . In this figure P cc (u) > P nn (u) for all u > 0, as 3λ ≤ 2μ + 1. The consensus is weak, . This figure shows all simple cases that hold for any λ ∈ (0, 1), μ > 0, states of different type. Inequalities 1, 4, and 5 imply that the probability of transition from the consensus state to a non-consensus one is always less than the probability of conservation of a nucleotide or the probability of transition from a non-consensus state to the consensus one. Inequalities 3 and 6 show that the probability of transition between two different non-consensus states is always less than the probability of state conservation. All these results correlate well with our intuitive ideas about an evolutionary model with selective pressure. Note that relations 1-6 hold for any positive μ, λ ∈ (0, 1) and for any time period.

Interesting cases
The most interesting case is the conservation of the consensus or a non-consensus nucleotide. Recall that P cc (h) = Thus, the relation between the probabilities of conservation of the consensus state and conservation of a nonconsensus state during the time u depends on the relation between 3λ and 2μ + 1. We obtained the following result.
The second relation partly describes the case of the weak consensus (Fig. 3). Since 3λ + 1 > 2(μ + 1), the stationary distribution of the consensus state is estimated from above as At the same time, the first relation holds both in the case of strong consensus π c ≥ 1/2 (Fig. 1) and in the case of weak consensus 1/4 <π c < 1/2 (Fig. 2).
The second interesting result is a relation between the probability of conservation of the non-consensus nucleo- ) .
Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 1/6, μ = 2 Figure 4 Transitional probabilities P cn , P cc , P nn , P nc , P nm for λ = 1/ 6, μ = 2. Here we have P cc (u) > P nn (u) for all u > 0, since 3λ ≤ 2μ +1. The consensus is strong and the stationary distribution of the consensus state is π c = 2/3. Since μ > λ, we have P cn (u) <P nm (u). The probability of substitutions between nonconsensus nucleotides P nm (u) is not monotone with a maximum point at u 2 = 2 log 7/11 ≈ 0.354. Since μ > 1, transitions between non-consensus nucleotides occur more frequently than transitions from a non-consensus to the consensus nucleotide, P nc (u) <P nm (u) on the time interval u ∈ (0, u 1 ), whereas P nc (u) > P nm (u) for u > u 1 . This figure also shows the relation between P nn (u) and P nc (u). tide and the probability of transition from a non-consensus state to the consensus state. We have where u 0 ≡ u 0 (λ, μ) is the solution of equation P nn (u) = P nc (u). It can be shown that u 0 is always positive. This relation shows that the probability of conservation of a nonconsensus nucleotide is less than the probability of transition to the consensus nucleotide from a non-consensus one for sufficiently long time u > u 0 . However, on the interval (0, u 0 ) the opposite inequality holds (see Fig. 1, 2).

Less interesting cases
It remains to study the relations between P nm and P cn , P nc . From the practical point of view, these cases are not very interesting, since the events of transition between different non-consensus states can be hardly observed on practice. However, we consider them for the sake of completeness.
The following relations imply that the transitional probability from the consensus to a non-consensus state could be less or greater than the transitional probability between two different non-consensus states depending on the relation between μ and λ (see Fig. 2

and 3):
Indeed, if μ > λ, then δ > α and the transition rate d between the non-consensus states is greater than the transition rate from the consensus state to non-consensus. Clearly, in this case P nm has to be greater than P cn (see Fig.  2 and Fig. 3).
The next case concerns the relation P nm and P nc for λ ∈ (0, 1) (see Fig. 1
An interesting fact is that the probability P nm (u) is not monotone in u for μ > λ (see Fig. 4). If μ > λ, then P nm increases on (0, u 2 ) and decreases for u > u 2 , where Finally, consider the case μ = λ. As it has been shown above, P cc (u) > P nn (u) for u > 0 (see Fig. 5). It is easy to see that in this case P nm ≡ P cn . The function P nm (u) is monotonically increasing for all u. Note that in the degenerate case λ = μ = 1 all states are equiprobable, the stationary distribution of the process π = (1/4, 1/4, 1/4, 1/4) and P nn = P cc , P nc = P cn = P nm (see Fig. 6).

Generalization
The obtained results can be generalized on the following model.  Transitional probabilities for λ = μ = 1/2

Figure 5
Transitional probabilities for λ = μ = 1/2. This case is similar to the one of Fig. 1, since here 3λ ≤ 2μ + 1 and, consequently, the probability of consensus conservation is always greater than the probability of conservation of a nonconsensus nucleotide, P cc (u) > P nn (u) for all u > 0. As λ = μ, we have P nm = P cn . In this case we have two groups of M consensus and N non-consensus states. We use the same notation for subscripts of transitional probabilities . If c and d denote different consensus states, n and m stand for different non-consensus states as before, we have the following transitional probabilities: Thus, this process has the stationary distribution = (π c , ..., π c , π n , ..., π n ) with Clearly, β > γ, since the first M states are the consensus ones. Next, compare the probabilities of conservation of the consensus and non-consensus states P cc and P nn . In a similar way, we obtain that there exists t* > 0 such that P cc (t) <P nn (t) for t ∈ (0, t*). Analyzing the probabilities P cc and P nn we can show that this is possible only for δ (N -1) -α (M -1) + β M -γ N < 0. Then the frequency of the consensus state (stationary distribution of consensus) is estimated from above as This condition coincides with the condition on weak consensus obtained for the case of four states with one consensus (M = 1, N = 3) considered above.

Comparison with the Molecular Evolution Theory
In the framework of the molecular evolution theory, the element a ij of the transition rate matrix is considered to be proportional to the product of the mutation rate p ij and the probability of fixation of a mutation f ij , a ij = kp ij f ij , where k is an arbitrary scaling constant [47]. As in [45] we start from the simplest Jukes-Cantor (1 -p)-scheme, to which we introduce selection. Thus, we ignore the difference between transitions and transversions in p ij .
TFBS regulating different genes in the same genome most likely evolve independently and thus the nucleotide composition π i at the respective positions of different TFBS occurrences approximates the equilibrium frequencies.
With these equilibrium frequencies at hand it is possible to relate p ij and f ij with the equation (see [47]) In our case, for the substitution rate between three nonconsensus positions we obtain π i = π j = π n and p ij = p ji = p nm , which yields f nm ∝ 1 by the l'Hôpital rule as in [47]. Thus, For the substitutions between non-consensus and consensus positions, r nc , both the selection preferences and mutation asymmetry come into consideration. In this case the "asymmetry constant" λ is crucial, which satisfies the inequality r cn /r nc = p n /p c = λ < 1. The following expression is valid: Degenerate case λ = μ = 1
If the mutation rate is symmetric, p nc = p cn , then f cn ∝ λ log(1/λ)/(1 -λ). Conversely, for the non-consensus to consensus substitution f nc ∝ log(1/λ)/(1 -λ) = f cn /λ, the greater flux from a non-consensus state to the consensus maintains a greater consensus frequency. Note that in alignments of sites in a single genome we can observe only the equilibrium constants π n , π c (which actually are rather rough approximations), thus the assumption p nc = p cn may be too strong, and the above general formula for f cn might be more relevant.
The coefficients in the matrix A for the symmetric mutation rates are given by If the background mutation rate is identical for all consensus and non-consensus nucleotides, we obtain p cn = p nc = p nm = p and p may be merged with the constant k. In this most simple case we obtain and, consequently, β > δ > α . This is the simplest generalization of the Jukes-Cantor model for the case with introduced selection.

Conservation of non-consensus nucleotides at the reduced mutation rate
Previously [8] we have observed that non-consensus nucleotides may be highly conserved in the alignments of orthologous TFBS in bacterial genomes. On the other hand, as shown above, the non-consensus nucleotide in the alignment of orthologous sites from different species cannot be more conserved than the consensus nucleotide if we adopt the Kimura model with the identical mutation rate for all pairs.
One way to explain the observation made in [8] is to drop the equivalence of non-consensus nucleotides and assume that different non-consensus nucleotides are under different selection in the sites regulating different rows of orthologous genes. Interestingly, it is possible to observe such conservation pattern in the model with an identical probabilities of mutation and fixation for all non-consensus nucleotides in all orthologous site rows, with a single preferred nucleotide, the consensus, the same in all sites. The only necessary relaxation of the model is to drop the condition p cn = p nc = p nm = p for the condition p nm <p nc = p cn . Doing this it is possible to satisfy the inequalities determining Case 2 of "Interesting Cases": μ <λ, 3λ > 2μ + 1.
The condition p nm <p nc = p cn means that the rate of direct mutations from one non-consensus nucleotide to another one is lower than the mutation rate in pairs involving the consensus nucleotide. A possible example of such specific reduction of the mutation rate comes from correlations of nucleotides occupying different positions of the same binding site. Assume that two positions within the site are not independent and must be occupied by correlated nucleotides. These may be, e.g., adjacent positions in a DNA site or base-paired positions in an RNA structure. Assume also that if any of the two positions is occupied by the consensus nucleotide, the correlating nucleotide may be arbitrary. Conversely, if one position is occupied by a non-consensus nucleotide, the other position should be occupied with some specific nucleotide, e.g. the complementary one in the case of an RNA structure.
In this case, the preferred pathway from a non-consensus nucleotide to another non-consensus nucleotide would become not via a direct mutation, but via an intermediate mutation into the consensus nucleotide. For example, if "C" is both cytosine and consensus, and for some case this position is correlated with another one, so that only A-A, T-T, and G-G pairs involving the non-consensus nucleotides at the first position are allowed, then only mutation A>C is valid, whereas mutations A>T and A>G are forbidden, and may occur via mutation into C and then the compensating mutation in the second position of the site. This simple model to some extent agrees with recent studies demonstrating that protein-DNA interactions are rather complex and probably may not be described by a simple position-independent model such as a positional weight matrix [49].
At the same time, this effect is not in contradiction with the uniform distribution of non-consensus nucleotides obtained from alignments of multiple TFBS regulating dif-  ferent genes in the same genome. It also allows for high conservation of some non-consensus nucleotides in the alignment of orthologous transcription factor binding sites from different species. Indeed, if the context-dependent pattern of conservation is specific to a particular position in a particular set of orthologous TFBSs, exactly this type of behavior may be expected.

Conclusion
The evolutionary model derived here can be generalized to the case of M consensus and N non-consensus nucleotides that are equiprobable, respectively. Thus, for M = 2, N = 2 and α = β, γ = δ. we get the case of neutral evolution with different rates of transitions and transversions [40]. If M = 1, N = 3, then we have the evolution model under constant selective pressure considered in this paper.
One somewhat non-obvious feature of this model is the existence of a combination of the rate of transition between non-consensus states μ and the rate of transition from the consensus state λ, for which the conservation of the consensus nucleotide in a sequence alignment can be lower than the conservation of a neutral (non-consensus) nucleotide. However, this can be observed only for the case of a weak consensus 1/4 <N c < 1/2 and for a relatively short time interval u ∈ (0, u**), and a mutation matrix with different elements, e.g. context-dependent.
Models of this type can be applied not only for the analysis of regulatory sites, but in other situations, e.g. for the analysis of functional sites in proteins [50] or analysis of evolution in the case of nucleotide biases [41][42][43][44]51,52].
In future work, we intend to estimate the parameters of the model from the data based on the method of maximum likelihood trees. Consider a Markov process X(t) describing the evolution at some fixed position. We assume that we observe only the endpoints of k different paths of X(t) depending on the evolution branch. In other words, our data is a set of k nucleotides at a fixed position in k genomes. The parameter λ can be easily estimated by the 'naive' estimator = (1/N c -1)/3, where N c is the frequency of the consensus nucleotide at the fixed position within k genomes. This estimator is obtained from the formula for the stationary distribution of the consensus state π c = 1/(3λ + 1). However, it is a good approximation of the true value of λ. The preliminary analysis of numerically simulated data shows that the performance of the maximum likelihood estimator is also rather good. On the other hand, it is not clear whether it is possible to construct an explicit estimator for the rate of transition between non-consensus states μ from the data, although μ participates in the expression for transitional probabilities.
The inspiration for the constructed model was the example of a set of orthologous transcription factors interacting with the cognate regulatory regions. However, it appears that evolution of sequences under a low selection pressure is a more widespread phenomenon. Indeed, a recent study [53] demonstrated that the force causing conservation of some non-coding genome regions of human, mouse and chimpanzee can be explained by a rather small selective pressure at the genomic level. A similar problem appears in the context of the CG composition of genomic regions [54]. Again, in this case the selective pressure appears to be low, although unlike the previous examples, here two nucleotides become selected for rather than one consensus nucleotide. The appropriate model in this case includes differences in the transition and the transversion rates, which makes the model more complicated, and probably would result in more complex time behavior of the substitution probabilities.
Anyhow, the emerging huge amount of data on orthologous non-coding regions, which has became available recently, brings forward a problem of modelling evolution with a selection pressure at finite times.

Methods
The numerical simulations and figures were produced using Matlab.
Publish with Bio Med Central and every scientist can read your work free of charge