Selection Intensity Against Deleterious Mutations in RNA Secondary Structures and Rate of Compensatory Nucleotide Substitutions

Corresponding author: Wolfgang Stephan, Department of Evolutionary Biology, University of Munich, Luisenstrasse 14, 80333 Munich, Germany., stephan{at}zi.biologie.uni-muenchen.de (E-mail)

Communicating editor: G. B. GOLDING

	ABSTRACT

TOP ABSTRACT THEORY COMPUTER SIMULATIONS DISCUSSION LITERATURE CITED

A two-locus model of reversible mutations with compensatory fitness interactions is presented; single mutations are assumed to be deleterious but neutral in appropriate combinations. The expectation of the time of compensatory nucleotide substitutions is calculated analytically for the case of tight linkage between sites. It is shown that selection increases the substitution time dramatically when selection intensity Ns > 1, where N is the diploid population size and s the selection coefficient. Computer simulations demonstrate that recombination increases the substitution time, but the effect of recombination is small when selection is weak. The amount of linkage disequilibrium generated in the process of compensatory substitution is also investigated. It is shown that significant linkage disequilibrium is expected to be rare in natural populations. The model is applied to the mRNA secondary structure of the bicoid 3' untranslated region of Drosophila. It is concluded that average selection intensity Ns against single deleterious mutations is not likely to be much larger than 1.

MODELS of compensatory evolution involve mutations from two or more loci. These mutations are assumed to be deleterious when they occur independently, but, in combination, they (at least partially) compensate each other for their deleterious effects. KIMURA 1985 proposed a two-locus, two-allele model, with alleles A and a at the first locus and B and b at the second locus, as illustrated in Fig 1A. He assumed that the two intermediate haplotypes, Ab and aB, are deleterious, while the wild-type AB and the double mutant ab are neutral. He further assumed that selection intensity against the deleterious intermediates is so strong that their frequencies in a population are very low. Accordingly, he considered only unidirectional mutations from A to a and from B to b and ignored back mutations (Fig 2A).

View larger version (13K): In this window In a new window Download PPT slide   Figure 1. Models of compensatory evolution. (A) KIMURA's (1985) model in which very strong selection is assumed and only unidirectional mutations are considered. (B) The model used in this article. Bidirectional mutations are considered. The fitnesses and frequencies are presented in parentheses.

View larger version (26K): In this window In a new window Download PPT slide   Figure 2. (A) Illustration of the process of compensatory substitution when selection is very strong and the intermediates Ab and aB are maintained in low frequencies. Solid circles present births of ab double mutants. (B) Illustration of the process of compensatory substitution when selection is weak. Ab and aB sometimes fix in the population.

An important example of compensatory evolution is found in RNA secondary structures. In single-stranded RNAs, Watson-Crick (WC) pairing of complementary nucleotide bases is the basic mechanism in the formation of stem-loop structures. It is believed that an individual mutation that breaks up a WC pairing is deleterious and that a second "compensatory" mutation at the complementary site can reestablish pairing and restore fitness. The relatively simple pattern of intramolecular WC base-pairing involved in RNA structures has made them a suitable model for the study of compensatory evolution (STEPHAN and KIRBY 1993 ; GOLDING 1994 ; SCHONIGER and VON HAESELER 1994 ; KIRBY et al. 1995 ; MUSE 1995 ; RZHETSKY 1995 ; TILLIER and COLLINS 1995 ; STEPHAN 1996 ; HIGGS 2000 ; SAVILL et al. 2001 ). Most of these authors analyzed rRNA secondary structures.

Phylogenetic analysis has revealed a number of compensatory nucleotide changes between species. On the other hand, however, mismatches, including not only GU wobble pairs but also other noncanonical pairs, are frequently observed, indicating that selection against deleterious intermediates may not be very strong (ROUSSET et al. 1991 ; PARSCH et al. 2000 ). This suggests that Kimura's compensatory evolution model, which assumes strong selection against deleterious intermediates, may not be generally applicable to the evolution of RNA secondary structures.

In this article, a compensatory evolution model is described in which selection against deleterious single mutations is not necessarily strong but covers a broad range of selection coefficients. Under weak selection, deleterious haplotypes may increase in frequency and even fix in the population as illustrated in a simulation run shown in Fig 2B. Since back mutations play an important role when the frequencies of deleterious haplotypes become large, we use a two-locus, two-allele model in which bidirectional mutations are considered (Fig 1B). This model is different from those used by KIMURA 1985 , IIZUKA and TAKEFU 1996 , and STEPHAN 1996 , who assumed that selection against deleterious intermediates is so strong that the mutation process may be considered unidirectional (see above). Our analysis is also different from that of HIGGS 1998 , who assumed bidirectional mutation but analyzed the model only in the parameter range of strong selection (i.e., Ns >> 1, where N is the diploid population size and s the selection coefficient).

Our goal is to calculate the time to proceed from the wild-type state AB to the fixation of the double mutant ab in this bidirectional mutation model. It should be noted, however, that neither AB nor ab is an absorbing state. That means we do not calculate a "fixation" time sensu stricto but do calculate the time to "flip" back and forth between AB and ab. This is possible as the nucleotide mutation rate is much smaller than 1/N (see below). To describe the transition between AB and ab, it is therefore reasonable to use the term "fixation" (or "substitution"). We present analytical results for the rate of this compensatory substitution event when there is no recombination between the two loci and use computer simulation to obtain this time under the influence of recombination. The theoretical results are applied to DNA sequence data of the bicoid 3' untranslated region (UTR) of Drosophila to estimate the selection intensity against deleterious mutations.

Another purpose of this article is to evaluate the amount of linkage disequilibrium generated by the compensatory evolution model. It is known that epistatic selection may produce significant linkage disequilibrium in natural populations (LEWONTIN 1974 ). SCHAEFFER and MILLER 1993 reported linkage disequilibria in two clusters of DNA polymorphisms in introns of Adh in Drosophila pseudoobscura. These disequilibria are likely due to epistatic selection maintaining pre-mRNA secondary structures (KIRBY et al. 1995 ). Here we examine whether these findings are consistent with our model of compensatory evolution.

	THEORY

TOP ABSTRACT THEORY COMPUTER SIMULATIONS DISCUSSION LITERATURE CITED

Consider a two-locus, two-allele model for a randomly mating population with N diploids. There are alleles A and a at the first locus and B and b at the second locus. The (bidirectional) mutation rate between A and a is given by µ₁ and that between B and b is given by µ₂ (Fig 1B). The mutation rates are assumed to be much smaller than 1/(2N), as is the case for nucleotide substitution rates. The relative fitnesses of the haplotypes AB, Ab, aB, and ab, are given by 1, 1 - s₁, 1 - s₂, and 1, respectively. Effects of fitness are assumed additive within locus and therefore the diploid model is equivalent to a haploid model. The haplotype frequencies are denoted by x₀, x₁, x₂, and x₃, respectively.

In the following, we consider the process of compensatory substitution under the joint action of drift, selection, mutation, and recombination. We assume the process starts at x₀ = 1 at t = 0. First, mutations in AB produce Ab and aB types. Next, mutations in Ab and aB may create ab, and some of them fix in the population. Recombination between Ab and aB may also produce ab. We are interested in the expected time of compensatory substitution, defined as the time from t = 0 (when the system is at state x₀ = 1) to the time point when the double mutant ab is fixed (x₃ = 1). Because the mutation process is bidirectional, the latter state is not an absorption state. The time we are calculating is therefore a first passage time. Because of the assumption µ_i << 1/(2N), double mutants ab either get lost by drift or go to fixation; the proportion of ab haplotypes reverting back to Ab or aB is very small. Recombination can only retard the fixation process (KIMURA 1985 ; STEPHAN 1996 ). The fixation time can be divided into two parts, T₁ and T₂. T₁ is the waiting time for the "successful" double mutant ab (that will eventually get fixed) to appear in the population, and T₂ is the time from the appearance of ab to the fixation event (see Fig 2).

Symmetrical model:
Consider a symmetrical model with s = s₁ = s₂ and µ = µ₁ = µ₂. In the analytical derivations, recombination is neglected. The four haplotypes are divided into two groups: one consists of AB and ab, and the other is the group of the deleterious intermediates Ab and aB. Let X be the frequency of the group of deleterious intermediates (X = x₁ + x₂) and Y the frequency of the other group (Y = x₀ + x₃). Denote the distribution of X by (X). In phase 1 when the system waits for a successful double mutant to appear, it will reach a quasi-equilibrium after a short initial period. At this quasi-equilibrium, its distribution is given approximately by

(1)

where = 4Nµ and C is a constant determined such that ¹₀(X)dX = 1 (WRIGHT 1931 , WRIGHT 1937 ). Equation 1, however, is not valid very shortly after t = 0 because we use the initial condition X = 0 (i.e., Y = 1).

Thus, we assume that when a new ab appears by a mutation in Ab or aB, the distribution of X before the mutation is given by (1). After the mutation, Y (= 1 - X) changes to Y' = Y + . Let p be the fixation probability of the Y group. Since the mutation rate is assumed to be very small, p can be approximated by

(2)

(KIMURA 1962 ), and the probability, p', that the new ab mutant fixes becomes

(3)

Therefore, since a new ab appears with probability 2NµX per generation, the expected number of ab that will fix in the population is given by

(4)

and the waiting time for the appearance of a successful double mutant ab becomes

(5a)

This result suggests that the waiting time for the appearance of a successful double mutant, T₁, is approximately exponentially distributed as

(5b)

In the case of neutrality, since p' is 1/(2N), the time becomes

(6)

because ¹₀X(X)dX = in this symmetrical model. In a one-locus neutral model, a substitution occurs every 1/µ generations on average (KIMURA 1983 ). Thus, T_1neu can be considered as the expected waiting time for two independent neutral substitutions.

When selection is very strong, X is maintained in very low frequency (approximately at the deterministic mutation-selection balance). The expected frequency is then given by

(7)

If 2µ/s << 1, p' is 1/(2N) because the average fitness of the population is nearly one. Therefore, T₁ becomes

(8)

Equation 8 agrees with Equation 8b in STEPHAN 1996 , which was obtained for the expected waiting time in Kimura's model of unidirectional mutation pressure by a different method.

T₂ is the time from the appearance of a successful double mutant haplotype to the fixation event. In the case of neutrality, the expectation of T₂ is 4N (KIMURA and OHTA 1969 ) since we assumed µ_i << 1/(2N). When selection is very strong, T₂ may be close to 4N again. That is expected because x₁ and x₂ are very small and ab has almost no selective advantage in the population. T₂ for moderate selection intensities is expected to be smaller than 4N because ab has a selective advantage when x₁ and x₂ are not very small.

General model:
Consider a general model where µ₁ µ₂ and/or s₁ s₂. Recombination is again neglected in the analytical treatment. In this case, it is necessary to investigate the frequency distributions of Ab and aB separately. Let ₁(x) and ₂(x) be the frequency distributions of Ab and aB, respectively. To obtain these two distributions, we reconsider the symmetrical model, where the distribution of the sum of x₁ and x₂ is given by (1). In a population with N diploids, the probability that i haplotypes are Ab or aB is given by

(9a)

when 0 < i < 2N, and

(9b)

Let P₁(i) and P₂(i) be the probability distributions of the numbers of Ab and aB, respectively. If we assume that µ is so small and Ns is so large that Ab and aB do not coexist frequently, P₁(i) and P₂(i) are given approximately by

(10a)

for i > 0 and by

(10b)

These results indicate that ₁(x) and ₂(x) may be given approximately by

(11)

for x > 1/(4N).

Next we consider ₁(x) and ₂(x) in the general model, where µ₁ µ₂ and/or s₁ s₂. If we assume that Ab and aB do not coexist frequently, the frequencies of Ab and aB follow two independent distributions. From the arguments in Equation 9aEquation 9b HREF="#FD11">Equation 11, it is expected that ₁(x) and ₂(x) are given for x > 1/(4N) by

(12)

where

(13a)

and

(13b)

where ₁ = 4Nµ₁ and ₂ = 4Nµ₂. C₁ and C₂ are constants that are determined such that ¹₀₁(x)dx = 1 and ¹₀₂(x)dx = 1 respectively.

Denote by p^'₁ the fixation probability of a new ab produced by a mutation in Ab given x₁. Since we assume that Ab and aB do not coexist in the population, Y is given by 1 - x₁ before the mutation and Y' = Y + after the mutation. Then, p^'₁ is given by

(14a)

In the same way, the fixation probability of a new ab produced by a mutation in aB given Y' = 1 - x₂ + becomes

(14b)

Therefore, the expected number of ab that will fix in the population per generation is given by

(15)

and T₁ is given by

(16)

When selection is very strong, T₁ becomes

(17)

which agrees with Equation 8b in STEPHAN 1996 . No simple formula for the case of neutrality can be obtained in this way because we assume that Ns_i is so large that Ab and aB do not coexist frequently.

T₂ for the general model is the same as for the symmetrical model. That is, T₂ 4N when selection is very strong and T₂ < 4N when selection intensity is moderate.

	COMPUTER SIMULATIONS

TOP ABSTRACT THEORY COMPUTER SIMULATIONS DISCUSSION LITERATURE CITED

Computer simulations were carried out for the following reasons. The first one is to check the theoretical results for T₁ shown above, because we use the following assumptions in the derivation. We use approximate formulas for the distribution of the haplotype frequencies ignoring the initial condition (X = 0 at t = 0). In the general asymmetric model, we assume that two haplotypes, Ab and aB, do not coexist. The second reason is to examine the effect of recombination on T₁ under a broad range of selection coefficients. In contrast to the strong selection case (KIMURA 1985 ; STEPHAN 1996 ), we were unable to find analytical expressions for T₁ with recombination when selection is weak. T₂ is also investigated by simulations with and without recombination. Another purpose of the simulations is to evaluate the amount of linkage disequilibrium generated in the process of compensatory substitution.

Monte Carlo simulations with mutation, selection, recombination, and random genetic drift were conducted in a constant size population of N diploids as follows. Each replication of the simulations starts from the initial condition (x₀, x₁, x₂, x₃) = (1, 0, 0, 0). In every generation, the frequencies are determined by the pseudosampling method (KIMURA and TAKAHATA 1983 ). The recombination rate between the two loci is assumed to be r per generation. Every 4N generations, the frequencies (x₀, x₁, x₂, x₃) are scored to investigate their frequency distributions. The amount of linkage disequilibrium (D = x₀x₃ - x₁x₂) is also calculated. Each replication ends when (x₀, x₁, x₂, x₃) = (0, 0, 0, 1) is reached for the first time, and T₁ and T₂ are recorded. When 0 t' T₂ (with t' = t - T₁), (x₀, x₁, x₂, x₃) is recorded every N/10 generation to calculate D.

Computer simulations with no recombination were conducted assuming N = 100, s = s₁ = s₂, and µ = µ₁ = µ₂, and the results for T₁ and T₂ are summarized in Table 1. The theoretical results for T₁ in the symmetrical model were compared with the simulation results with no recombination (Fig 3). In Fig 3A, the theoretical expectation calculated by (5a) is shown with the simulation results of = 0.01. Ns was changed from 0 to 5. The theory is in very good agreement with the results of the simulations. Fig 3B shows the results of theory and simulations for a relatively high mutation rate ( = 0.1) although our model assumes very low mutation rates. It is shown that Equation 5a gives a quite good approximation even for = 0.1, unless Ns = 0. The underestimation for Ns = 0 may occur because frequent recurrent mutations reduced the fixation probability of ab. The degree of reduction is expected to be relatively small when T₁ is large. For all parameter sets examined, the standard deviation of T₁ is similar to the mean (Table 1), supporting the exponential distribution of T₁ as suggested by (5b). Similar results were obtained from computer simulations with N = 1000 (data not shown).

View larger version (21K): In this window In a new window Download PPT slide   Figure 3. Relationship between T1 and Ns under the symmetrical model without recombination. The theoretical expectation is obtained from (5a) and represented by a solid curve. The simulation results are based on Table 1 and presented by open circles. (A) Results for = 0.01. (B) Results for = 0.1.

The effect of recombination on T₁ was also investigated by computer simulations with = 0.01 and 0.1 (Table 1). The effect of recombination is very small when selection is weak. The effect, however, can be seen for selection intensity Ns > 1. In Fig 4, a clear positive correlation is observed between 4Nr and T₁ when Ns = 2 and 5. These results are consistent with those of KIMURA 1985 and STEPHAN 1996 . Similar relationships between 4Nr and T₁ were observed in other two-locus models with epistatic interactions (MICHALAKIS and SLATKIN 1996 ; CHRISTIANSEN et al. 1998 ).

View larger version (44K): In this window In a new window Download PPT slide   Figure 4. Effect of recombination on T1 under the symmetrical model with = 0.01. The results of computer simulations are from Table 1.

The results for T₂ are also shown in Table 1. T₂ is much smaller than T₁ in all the parameter sets investigated. First, we consider T₂ with no recombination. In the case of neutrality, T₂ obtained from simulations for = 0.01 is close to 400 as expected from the theory, although T₂ for = 0.1 is a little >4N. When selection is very strong (Ns 5), T₂ is close to 4N again. T₂ for moderate selection intensity is <4N.

T₂ may be negatively correlated with recombination rate (Table 1). The degree of reduction in T₂ is larger when selection is stronger. T₂ for 4Nr = 0 is similar to that for 4Nr = 10 when Ns 1, while T₂ for 4Nr = 10 is much smaller than that for no recombination when Ns 2. The result for the T₂ phase may be understood as follows. Since recombination usually occurs between AB and ab in T₂, it reduces the fixation probability of ab and increases T₁. As the fixation probability is reduced, only ab, which increases its frequency quickly, can successfully fix in the population.

Theoretical results for T₁ in the general model were compared with the results of computer simulations, and the results for ₁ = 0.02 and ₂ = 0.01 are shown in Fig 5 with the theoretical expectations calculated from (16). It is shown that the theoretical expectations are in good agreement with the results of simulations when Ns₁ 0.5 and Ns₂ 0.5. If one of the selection intensities is very small, Equation 16 overestimates T₁ because the assumption used to derive (16) does not hold. In the derivation, we obtained approximate formulas for ₁(x) and ₂(x) under the assumption that Ab and aB do not coexist at the same time. If one of the selection intensities, say Ns₁, is very small, ₁(x) given by (12) does not agree with the frequency distribution of Ab obtained by simulations. Similar results were obtained for other values of ₁ and ₂. Good agreement between the theory and simulations was observed when neither Ns₁ nor Ns₂ is small (data not shown).

View larger version (17K): In this window In a new window Download PPT slide   Figure 5. Relationship between T1 and selection intensities Ns1 and Ns2 under the general model without recombination. The theoretical expectations are obtained from (16) and represented by solid curves. 1 = 0.02 and 2 = 0.01 are assumed. Solid circles, solid squares, open circles, and open squares represent the results of simulations for Ns1 = 0, 0.5, 1, and 2, respectively.

Table 2 shows the amount of linkage disequilibrium (D = x₀x₃ - x₁x₂) in the process of compensatory substitution obtained by simulations under the symmetrical model with no recombination. D was calculated every 4N generations as long as 0 < t < T₁. For this interval of 4N generations, we found almost no correlation between sampling. Mean and variance were calculated for all runs. For = 0.01 the level of linkage disequilibrium is extremely low. Even for = 0.1 D is very small although much larger than that for = 0.01. This indicates that almost no linkage disequilibrium is expected during the waiting time for the appearance of a successful double mutant haplotype. Simulations with recombination showed that recombination reduces the level of linkage disequilibrium (data not shown).

On the other hand, strong positive linkage disequilibrium is generated after a successful double mutant has appeared and is on its way to fixation (0 < t' < T₂). The average of linkage disequilibrium for t' 0 is plotted in Fig 6. Selection increases the level of linkage disequilibrium significantly (Fig 6A and Fig B). Each distribution of D has a peak near t' = 100. As Ns increases, the peak is getting larger and seems to saturate at D 0.2. Note that the theoretical maximum value of D is 0.25, which is reached when x₀ = x₃ = 0.5 and x₁ = x₂ = 0. When selection is weak, the level of linkage disequilibrium is higher for = 0.1 than for = 0.01.

View larger version (75K): In this window In a new window Download PPT slide   Figure 6. Linkage disequilibrium in phase T2. (A) Results of computer simulations without recombination for = 0.01. The effect of selection intensity is investigated. (B) Results of computer simulations without recombination for = 0.1. (C) Results of computer simulations for = 0.01 and Ns = 2. The effect of recombination is investigated.

Fig 6C shows the effect of recombination on the level of linkage disequilibrium when = 0.01 and Ns = 2. As the recombination rate increases, the level of linkage disequilibrium is getting weaker. Similar results were obtained for other mutation rates and selection intensities (data not shown).

	DISCUSSION

TOP ABSTRACT THEORY COMPUTER SIMULATIONS DISCUSSION LITERATURE CITED

Theory:
We analyzed a model of reversible mutations with compensatory fitness interactions; i.e., single mutations are assumed to be deleterious but harmless (neutral) in appropriate combinations. In proceeding under mutation pressure, epistatic selection, and genetic drift from one fitness peak to another, a population must pass through a valley of lower individual fitness. This process of compensatory evolution is investigated by analytical approximation and computer simulation. Our model is more general than KIMURA's (1985), which assumed that mutation pressure is unidirectional and that deleterious intermediates are very strongly selected against. In contrast, our model covers a broad range of selection coefficients, including very small ones, and agrees with analyses of Kimura's model when selection is strong (KIMURA 1985 ; STEPHAN 1996 ).

As in these latter analyses, we focus on the expected time for the compensatory substitution process to go from one fitness peak to another. In addition, we study the structure of variation (linkage disequilibrium) within a population during this transition. The process of compensatory substitution is analyzed by dividing it into two phases defined by the time periods T₁ and T₂. T₁ is the waiting time for a successful double mutant haplotype (ab) that will fix in the population, and T₂ is the time from the appearance of a successful double mutant to the fixation event. The results of theory and computer simulations show that T₁ >> T₂ (Table 1) as suggested by STEPHAN 1996 .

The expectation of T₁ is obtained analytically for the case without recombination. It is shown that selection dramatically increases T₁ (Fig 3). Although the theory is based on the assumption of very low mutation rates, it is shown that it is also applicable to very high nucleotide mutation rates relative to 1/N (e.g., = 0.1) unless Ns is very small. Thus, our analysis is useful for natural populations for which is generally <<0.1 (KIMURA 1983 ; NEI 1987 ; GILLESPIE 1991 ). Our computer simulations demonstrate that there is almost no effect of recombination on T₁ when selection is relatively weak (Ns 1), while T₁ increases significantly with increasing recombination rates when selection is strong (Fig 4).

Parameter estimation:
An important parameter of the compensatory evolution model is the intensity of selection Ns against deleterious intermediates. In the following, we attempt to estimate this parameter for mRNA secondary structures. Our theoretical results predict that T₁ is much larger than T_1neu unless Ns is very small (Fig 3), suggesting that nucleotide substitutions occur very slowly in pairing regions of mRNA secondary structures. It is known that such regions are highly conserved between species in contrast to other (unpaired) regions (MUSE 1995 ; PARSCH et al. 2000 ). Thus, it may be possible to estimate Ns in pairing regions from DNA sequence comparisons.

We compare the rates of substitutions between species in pairing regions with those in regions that are considered selectively neutral. As an example, we analyze the bicoid 3' UTR of Drosophila. It has been shown that bicoid mRNA has a complex secondary structure in the 3' UTR (MACDONALD 1990 ; SEEGER and KAUFMAN 1990 ; FERRANDON et al. 1997 ; MACDONALD and KERR 1998 ). Based on the alignment of DNA sequences from nine Drosophila species, PARSCH et al. 2000 identified eight highly conserved pairing regions, of which seven have been supported by mutational analysis (FERRANDON et al. 1997 ; MACDONALD and KERR 1998 ). We consider these eight stems as pairing regions and the remainder of the 3' UTR as unpaired. It is also assumed that there is no selection in the unpaired regions. The total length of the paired segments is 142 nucleotides, corresponding to 71 bp.

We first compare the nucleotide sequences between D. melanogaster and D. simulans, using the alignment suggested by PARSCH et al. 2000 . In 142 nucleotides of the pairing regions, 3 nucleotide differences are observed and the number of substitutions per site (d_p) is estimated to be 0.0214 by the JUKES and CANTOR 1969 method (Table 3). In the remaining regions of the 3' UTR, 31 nucleotide differences are observed and the number of substitutions per site (d_n) is estimated to be 0.0435. This suggests that the rate of nucleotide substitutions is reduced by roughly a factor of 2 in the pairing regions in comparison with the remainder of the 3' UTR. For the pair of D. melanogaster and D. pseudoobscura, d_p and d_n are estimated to be 0.2099 and 0.5087, respectively. The ratio of d_n to d_p is 2.4, similar to that of the comparison between D. melanogaster and D. simulans.

Since it is assumed that the unpaired regions are selectively neutral, the ratio of d_n to d_p is comparable to T₁/T_1neu. In Fig 7, T₁/T_1neu is plotted as a function of Ns for the symmetrical model without recombination. T₁/T_1neu is calculated by (5a) and (6). The results show that T₁/T_1neu increases rapidly with increasing Ns, as soon as Ns > 1. This is particularly the case for 0.01. T₁/T_1neu depends weakly on when 0.01.

View larger version (18K): In this window In a new window Download PPT slide   Figure 7. Relationship between T1/T1neu and Ns for various values of .

In the bicoid 3' UTR of Drosophila (Table 3), the ratio of d_n to d_p is 2.0–2.4. The estimate of in the unpaired regions for a D. melanogaster population from Zimbabwe is 0.003 (J. F. BAINES, Y. CHEN and W. STEPHAN, unpublished results). Fig 7 suggests therefore that the observed ratio of d_n to d_p can be explained if Ns is 0.6–0.7. If we consider only the number of complete compensatory substitutions (WC/WC in Table 3) for the comparison between D. melanogaster and D. pseudoobscura, the ratio of d_n to d_p becomes 5 and the estimate of Ns 1.

There are, however, some caveats.

1. Ns could be larger than this estimate because d_n is underestimated if selection is acting in the regions that we consider as unpaired. There may be some evidence for weak selection in these regions. First, MACDONALD 1990 suggested the possibility of long-range pairings encompassing almost the entire 3' UTR. However, his suggestion was not supported by a strict phylogenetic analysis (PARSCH et al. 2000 ). Second, average silent divergence in the unpaired segments of the 3' UTR between D. melanogaster and D. simulans is 0.0435, which is about a factor of 2.5 lower than in the rest of the bicoid gene upstream of the 3' UTR (J. F. BAINES, Y. CHEN and W. STEPHAN, unpublished results). On the other hand, even if weak selection is acting in the unpaired regions, the estimate of Ns for the pairing regions does not increase much, as the average time of compensatory substitutions becomes extremely large for Ns > 1 when mutation pressure is relatively weak (i.e., 0.01; Fig 7). Thus, it may be concluded that the selection intensity in the pairing regions of the bicoid 3' UTR is on average not much >1. This estimate is similar to the estimate of average selection intensity for codon usage in Drosophila (AKASHI 1995 ).

2. To estimate the selection intensity, we used the average d_p of eight pairing regions. In other words, Ns is the average selection intensity of these eight pairing regions. PARSCH et al. 2000 found heterogeneity for d_p among pairing regions, caused by variation in both stem length and the physical distance between base-pairing residues. One reason is that long stems are under less selective constraints than short ones (see Fig 3A of PARSCH et al. 2000 ). Another factor is that short-range pairings (hairpins) experience a higher rate of evolution than long-range pairings because of the retarding effects of recombination when selection is sufficiently strong (Fig 4A of PARSCH et al. 2000 ). As a consequence, the estimates of Ns appear to vary substantially among pairing regions.

3. PARSCH et al. 2000 were able to distinguish the effect of stem length from that of physical distance when only pairing regions with covariations were considered. According to their Fig 4A, they found an approximately fivefold drop in the rate of compensatory evolution over a physical distance of nearly 200 bp between base-pairing residues. We have to ask whether such a large decrease of the rate of compensatory evolution is consistent with an estimate of Ns 1. Assuming that a physical distance of 200 bp of the bicoid 3' UTR corresponds to a value of 4Nr in the order of 10 (i.e., using the standard estimates of effective population size and recombination rate for D. melanogaster that are in the order of 10⁶ and 10^-8, respectively), this distance effect can be explained by our model only if Ns is 5 (see Table 1). Thus, it appears that this value is not compatible with our estimate of Ns 1 obtained without taking recombination into account. However, one has to keep in mind that for technical reasons the analysis of PARSCH et al. 2000 is based on pairing regions with covariations only. A much weaker distance effect would presumably result if all pairing regions were considered, including those with no covariations (PARSCH et al. 2000 ). A much larger data set and more sophisticated methods are needed to take the distance effect into account.

Linkage disequilibrium:
Strong linkage disequilibrium is sometimes considered as evidence for epistatic selection (LEWONTIN 1974 ). We investigated the amount of linkage disequilibrium in the process of compensatory substitution. It is demonstrated that the level of linkage disequilibrium is very low during time period T₁. Strong positive linkage disequilibrium, however, is observed in phase T₂ if selection is strong. This suggests that significant linkage disequilibria due to compensatory interactions should be rarely observed in natural populations because T₂ is much smaller than T₁ if selection is strong. On the other hand, if selection is weak, linkage disequilibrium is not very large even in phase T₂.

SCHAEFFER and MILLER 1993 detected two clusters of polymorphisms in the Adh introns of D. pseudoobscura that exhibit significant linkage disequilibrium. In both cases, the disequilibria are due to two highly diverged haplotypes, ha1 and ha2, that have been shown to form different pre-mRNA secondary structures (KIRBY et al. 1995 ). It was also revealed that this structural polymorphism has predated the species split of D. pseudoobscura, D. persimilis, and D. miranda because ha1 and ha2 are similar to D. persimilis and D. miranda haplotypes, respectively (KIRBY et al. 1995 ). This observation is not consistent with our results, which show that a compensatory substitution requires a long waiting time for a successful double mutant to occur and that the fixation event of ab follows relatively quickly. In other words, our model does not predict that the secondary-structure-forming haplotypes of Adh are maintained for such a long time, as observed in these species. This suggests that our model of compensatory evolution is either too simple, as it allows only two sites to undergo base changes, or that some additional form of selection (for instance, balancing selection) may be maintaining the haplotypes ha1 and ha2. While there is no evidence for the latter suggestion, the fact that in both examples the haplotypes were subject to significant rearrangement during evolutionary time (due to insertions and deletions of bases) may indicate that the underlying compensatory process is much more complicated than our two-locus model assumes. Therefore, to model such complex compensatory changes, models need to be developed that include compensatory insertions and deletions in addition to base substitutions.

	ACKNOWLEDGMENTS

We thank an anonymous reviewer for helpful suggestions that improved the manuscript. John F. Baines and Ying Chen kindly provided unpublished results. This research was supported in part by a fellowship from the Japan Society for the Promotion of Science to H.I. and by National Institutes of Health grant GM-58405 to W.S.

Manuscript received February 17, 2001; Accepted for publication June 4, 2001.

	LITERATURE CITED

TOP ABSTRACT THEORY COMPUTER SIMULATIONS DISCUSSION LITERATURE CITED

AKASHI, H., 1995  Inferring weak selection from patterns of polymorphism and divergence at "silent" sites in Drosophila DNA. Genetics 139:1067-1076[Abstract].

CHRISTIANSEN, F. B., S. P. OTTO, A. BERGMAN, and M. W. FELDMAN, 1998  Waiting with and without recombination: the time to production of a double mutant. Theor. Popul. Biol. 53:199-215[Medline].

FERRANDON, D., I. KOCH, E. WESTHOF, and C. NÜSSLEIN-VOLHARD, 1997  RNA-RNA interaction is required for the formation of specific bicoid mRNA 3' UTR-STAUFEN ribonucleoprotein particles. EMBO J. 16:1751-1758[Medline].

GILLESPIE, J. F., 1991 The Causes of Molecular Evolution. Oxford University Press, London.

GOLDING, B., 1994 Using maximum likelihood methods to infer selection from phylogenies, pp. 126–139 in Non-Neutral Evolution: Theories and Molecular Data, edited by B. GOLDING. Chapman & Hall, New York.

HIGGS, P. G., 1998  Compensatory neutral mutations and the evolution of RNA. Genetica 102(103):91-101.

HIGGS, P. G., 2000  RNA secondary structure: physical and computational aspects. Q. Rev. Biophys. 33:199-253[Medline].

IIZUKA, M. and M. TAKEFU, 1996  Average time until fixation of mutations with compensatory fitness interaction. Genes Genet. Syst. 71:167-173.

JUKES, T. H., and D. R. CANTOR, 1969 Evolution of protein molecules, pp. 21–132 in Mammalian Protein Metabolism, edited by H. N. MUNRO. Academic Press, New York.

KIMURA, M., 1962  On the probability of fixation of mutant genes in a population. Genetics 47:713-719[Free Full Text].

KIMURA, M., 1983 The Neutral Theory of Molecular Evolution. Cambridge University Press, Cambridge, UK.

KIMURA, M., 1985  The role of compensatory neutral mutations in molecular evolution. J. Genet. 64:7-19.

KIMURA, M. and T. OHTA, 1969  The average number of generations until fixation of a mutant gene in a finite population. Genetics 61:763-771[Free Full Text].

KIMURA, M. and N. TAKAHATA, 1983  Selective constraint in protein polymorphism: study of the effectively neutral mutation model by using an improved pseudosampling method. Proc. Natl. Acad. Sci. USA 80:1048-1052[Abstract/Free Full Text].

KIRBY, D. A., S. V. MUSE, and W. STEPHAN, 1995  Maintenance of pre-mRNA secondary structure by epistatic selection. Proc. Natl. Acad. Sci. USA 92:9047-9051[Abstract/Free Full Text].

LEWONTIN, R. C., 1974 The Genetic Basis of Evolutionary Change. Columbia University Press, New York.

MACDONALD, P. M., 1990  bicoid mRNA localization signal: phylogenetic conservation of function and RNA secondary structure. Development 110:161-171[Abstract].

MACDONALD, P. M. and K. KERR, 1998  Mutation analysis of an RNA recognition element that mediates localization of bicoid mRNA. Mol. Cell. Biol. 18:3788-3795[Abstract/Free Full Text].

MICHALAKIS, Y. and M. SLATKIN, 1996  Interaction of selection and recombination in the fixation of negative-epistatic genes. Genet. Res. 67:257-269[Medline].

MUSE, S. V., 1995  Evolutionary analysis of DNA sequences subject to constraints on secondary structure. Genetics 139:1429-1439[Abstract].

NEI, M., 1987 Molecular Evolutionary Genetics. Columbia University Press, New York.

PARSCH, J., J. M. BRAVERMAN, and W. STEPHAN, 2000  Comparative sequence analysis and patterns of covariation in RNA secondary structures. Genetics 154:909-921[Abstract/Free Full Text].

ROUSSET, F., M. PÉLANDAKIS, and M. SOLIGNAC, 1991  Evolution of compensatory substitutions through a G · U intermediate state in Drosophila rRNA. Proc. Natl. Acad. Sci. USA 88:10032-10036[Abstract/Free Full Text].

RZHETSKY, A., 1995  Estimating substitution rates in ribosomal RNA genes. Genetics 141:771-783[Abstract].

SAVILL, N. J., D. C. HOYLE, and P. G. HIGGS, 2001  RNA sequence evolution with secondary structure constraints: comparison of substitution rate models using maximum-likelihood methods. Genetics 157:399-411[Abstract/Free Full Text].

SCHAEFFER, S. W. and E. L. MILLER, 1993  Estimates of linkage disequilibrium and the recombination parameter determined from segregating nucleotide sites in the alcohol dehydrogenase region of Drosophila pseudoobscura.. Genetics 135:541-552[Abstract].

SCHÖNIGER, M. and A. VON HAESELER, 1994  A stochastic model for the evolution of autocorrelated DNA sequences. Mol. Phylogenet. Evol. 3:240-247[Medline].

SEEGER, M. A. and T. C. KAUFMAN, 1990  Molecular analysis of the bicoid gene from Drosophila pseudoobscura: identification of conserved domains within coding and noncoding regions of the bicoid mRNA. EMBO J. 9:2977-2987[Medline].

STEPHAN, W., 1996  The rate of compensatory evolution. Genetics 144:419-426[Abstract].

STEPHAN, W. and D. A. KIRBY, 1993  RNA folding in Drosophila shows a distance effect for compensatory fitness interactions. Genetics 135:97-103[Abstract].

TILLIER, E. R. M. and R. A. COLLINS, 1995  Neighbor-joining and maximum likelihood with RNA sequences: addressing interdependence of sites. Mol. Biol. Evol. 12:7-15.

WRIGHT, S., 1931  Evolution in Mendelian populations. Genetics 16:97-159[Free Full Text].

WRIGHT, S., 1937  The distribution of gene frequencies in populations. Proc. Natl. Acad. Sci. USA 23:307-320[Free Full Text].

This article has been cited by other articles:

				J. L. Knies, K. K. Dang, T. J. Vision, N. G. Hoffman, R. Swanstrom, and C. L. Burch Compensatory Evolution in RNA Secondary Structures Increases Substitution Rate Variation among Sites Mol. Biol. Evol., August 1, 2008; 25(8): 1778 - 1787. [Abstract] [Full Text] [PDF]

				K. R. Takahasi and H. Innan The Direction of Linkage Disequilibrium: A New Measure Based on the Ancestral-Derived Status of Segregating Alleles Genetics, July 1, 2008; 179(3): 1705 - 1712. [Abstract] [Full Text] [PDF]

				K. M. Teshima and H. Innan Neofunctionalization of Duplicated Genes Under the Pressure of Gene Conversion Genetics, March 1, 2008; 178(3): 1385 - 1398. [Abstract] [Full Text] [PDF]

				S. A. Shabalina, A. Y. Ogurtsov, and N. A. Spiridonov A periodic pattern of mRNA secondary structure created by the genetic code. Nucleic Acids Res., January 1, 2006; 34(8): 2428 - 2437. [Abstract] [Full Text] [PDF]

				Y. Iwasa, F. Michor, and M. A. Nowak Stochastic Tunnels in Evolutionary Dynamics Genetics, March 1, 2004; 166(3): 1571 - 1579. [Abstract] [Full Text] [PDF]

				J. F. Baines, J. Parsch, and W. Stephan Pleiotropic Effect of Disrupting a Conserved Sequence Involved in a Long-Range Compensatory Interaction in the Drosophila Adh Gene Genetics, January 1, 2004; 166(1): 237 - 242. [Abstract] [Full Text] [PDF]

				T. F. McCutchan, D. Rathore, and J. Li Compensatory Evolution in the Human Malaria Parasite Plasmodium ovale Genetics, January 1, 2004; 166(1): 637 - 640. [Abstract] [Full Text] [PDF]

				Y. Chen and W. Stephan Compensatory evolution of a precursor messenger RNA secondary structure in the Drosophila melanogaster Adh gene PNAS, September 30, 2003; 100(20): 11499 - 11504. [Abstract] [Full Text] [PDF]

				H. Innan A two-locus gene conversion model with selection and its application to the human RHCE and RHD genes PNAS, July 22, 2003; 100(15): 8793 - 8798. [Abstract] [Full Text] [PDF]