The underlying mechanism remains to be elucidated. In males, Wolbachia must somehow affect the paternal nucleus before it is shed from maturing sperm, resulting in paternal chromosome loss after fertilization. In females, the bacterium must somehow prevent this loss and thereby rescue the embryo. This conceptual distinction between the male and female sides of CI was formalized by Werren (1997) through the modification/rescue (mod/resc) terminology.

Interestingly, crosses between infected males and infected females can also be incompatible, if the two partners bear different bacteria. The mod and resc functions thus seem to interact specifically: different Wolbachia can harbor different mod resc pairs, that is, different compatibility types. We are interested in the process that allows compatibility types to evolve. In an earlier study (Charlat et al. 2001b), we showed that compatibility types are not constrained by stabilizing selection, if mod and resc are determined by different genes, which, we think, is a reasonable assumption (Poinsot et al. 2003; see also Kose and Karr 1995; Callaini et al. 1997; Tram and Sullivan 2002). We suggested that compatibility types could change through a two-step process: the first involving drift on mod variations and the second involving selection on resc variations. This work highlighted the need for simulation-based models that would incorporate the effects of nondeterministic evolutionary forces. The present article is based on a simulation program developed in an attempt to fulfill this condition, allowing us to follow the evolution of compatibility types under mutation, drift, and selection.

A Wolbachia variant is defined here by five parameters (summarized in Table 1). Two parameters define the compatibility type: the lock and the key. In practice, lock and key are modeled as sequences of S sites, with n possible states for each site (states 1, 2, … , n). In a cross between a male bearing Wolbachia i and a female bearing Wolbachia j, a compatibility score (C_ij) is calculated as the proportion of matching sites in the lock/key alignment. For example, with S = 10, we can have lock_i = 2222211111 and key_j = 2222222222, which gives the compatibility score C_ij = 0.5. Lock and key are comparable to the mod_C and resc_C parameters in our earlier study (Charlat et al. 2001b), but with the notable difference that the current model allows partial compatibility if S > 1. MI (mod intensity) is the efficiency of the mod function; it corresponds to what is often referred to as “CI level,” which can be measured in crosses between infected males and uninfected females. Biologically, MI can be seen as the proportion of infected males' sperm that is actually affected by CI. TE (transmission efficiency) is the proportion of infected eggs among those laid by infected females. Finally, FE (fitness effect) is the fecundity of infected females relative to that of uninfected females.

TABLE 1 Parameter definitions

The evolution of the TE and FE parameters is not in the focus of this study, as their evolutionary trajectories have been well described by the analytical approach (Turelli 1994): selection on Wolbachia always acts to increase infected females fecundity and maternal transmission efficiency. In other words, starting from any initial condition, these parameters will rapidly reach their maximum values under mutation and selection. However, FE and TE are still relevant to the model, as we can investigate the effect of setting upper limits <1 for these two parameters. For example, the model allows us to analyze the evolution of the lock and key parameters in populations where TE cannot exceed 0.9, that is, in populations where uninfected individuals can persist.

Similarly, the evolution of the MI parameter in panmictic populations has been previously worked out (Prout 1994; Turelli 1994): with random mating, MI is not under selection and evolves through drift only. However, keeping MI constant would impede the “realism” of our analysis. Consequently, the model includes the evolution of MI, although we do not focus on this aspect.

Figure 1.— How progeny is determined for any given cross. The male is infected by Wolbachia i, with parameters locki, keyi, MIi, TEi, and FEi (but only MIi and locki are relevant, since the other parameters are not expressed in males). The female is infected by (more ...)

To validate the model, we tested whether it could retrieve earlier results, derived from the analytical approach. We first verified the basic prediction that CI allows Wolbachia to invade uninfected host populations, and more efficiently so if MI is high (Figure 2). We then investigated the combined effects of MI, TE, and FE on invasion dynamics. Caspari and Watson (1959), Fine (1978), and Hoffmann et al. (1990) showed that if TE and/or FE are <1, Wolbachia does not invade unless it first reaches a threshold frequency depending on MI, TE, and FE. Above that point, infection frequency increases toward a stable infection frequency, which is not fixation if MI and TE are <1. As illustrated in Figure 3, the simulation and analytical approaches provide congruent predictions.

Figure 2.— Plot of invasion probability (estimated as the number of times over 1000 runs where Wolbachia finally got fixed) as a function of CI level (MI) and population size, for an initial frequency of 0.1. As expected, we observe that higher mod intensity allow (more ...)

Figure 3.— The curve is a plot of infection frequency at generation i + 1 as a function of frequency at generation i, for a Wolbachia with the following properties: MI = 0.9, TE = 0.8, FE = 0.96. Any point below the x = y line (dashed) indicates that infection frequency (more ...)

Before considering “realistic evolution,” where an initial population can freely change under mutation, selection, and drift, we present here the results of simulations conducted without random mutation: several variants with specific properties are initially introduced and their frequencies followed over 1000 generations. These particular cases will allow the reader to understand which sequence of events can lead to which population state. With these processes in mind, the evolution of populations is analyzed more realistically in the next section.

Figure 4.— The neutrality of MI and lock variations and their possible consequences on infection loss. Population size Ne = 1000. (A) The consequences of MI polymorphism. The three variants (MI = 0.9, 0.6, and 0.3) are initially introduced with respective frequencies (more ...)

Less explicit in earlier analyses is the fact that variations affecting the lock parameter can also lead to infection loss. As detailed elsewhere (Charlat et al. 2001b), mutations affecting the lock sequence are not subject to selection, although they can give rise to self-incompatible, or “suicidal” Wolbachia. Figure 4B illustrates how this can lead to infection loss: as new lock variants (totally or partially self-incompatible) get too frequent by drift, the fitness gain provided by CI to infected females is lowered, which eventually leads to infection loss.

In summary, neutral variations of MI make CI less deleterious to uninfected females, while neutral variations of the lock make CI more deleterious to infected females. The final effect is the same in the two cases: net benefit to infected cytoplasmic lines is reduced.

Figure 5.— The sequence of events leading to the invasion of new compatibility types. Population size Ne = 1000. Two variants (AA, lockA/keyA; BA, lockB/keyA; with lockA = 1111111111, lockB = 2222222222, keyA = 1111111111, and keyB = 2222222222) (more ...)

When such an equilibrium is reached, selection does not favor any Wolbachia variant in particular. However, simulation runs such as that presented in Figure 6 show that not all variants are maintained at stable frequencies. The frequency of BA appears to be stable, but f(AA) and f(AB) vary randomly and symmetrically, leading to the loss of either AB (Figure 6A) or AA (Figure 6B). If AB is lost, the population goes back to the initial neutral polymorphism with AA and BA. In contrast, if AA is lost, the population reaches a stable polymorphism, with BA and AB at equal frequencies.

Figure 6.— How a “suicidal” polymorphism can be maintained by balancing selection. Population size Ne = 1000. Two variants (AA, lockA/keyA; BA, lockB/keyA; with lockA = 1111111111, lockB = 2222222222, keyA = 1111111111, keyB = 2222222222) (more ...)

Notably, the initial conditions (the relative frequencies of BA and AA when AB is introduced) strongly affect the respective likelihoods of the outcomes described in Figure 6, A and B. The higher the initial f(BA)/f(AA) ratio is, the higher is the frequency reached by the mutant AB under selection (in the very first generations following its introduction) and, in turn, the more likely is the outcome illustrated in Figure 6B, which depends directly on the f(AA)/f(AB) ratio at the beginning of the “drift of AA and AB” stage (Figure 6, A and B, left).

A full understanding of the frequency variations during the drift of AA and AB stage is not straightforward (and not necessary for the following sections). The main point here is to show that a balanced suicidal polymorphism (illustrated in Figure 6B) can occur. The sequence of events leading to such balanced polymorphism involves a complex interaction between selection and drift: AA and AB are globally neutral relative to each other (although locally selection takes place) while the frequency of BA is locally and globally stabilized by selection. Figure 7 provides a more detailed explanation to the interested reader.

Figure 7.— How a population bearing AA, BA, and AB reacts to small random variations. The results were obtained with a simple deterministic model, that is, a model where selection takes place only after an initial and controlled “random” variation. (more ...)

• State 1: infection loss, when all individuals are uninfected.
• State 2: population extinction, when no offspring is produced, due to the fixation by drift of a suicidal Wolbachia.
• State 3: stability of compatibility types, when the initial lock/key pair is still predominant. Arbitrarily, we define “predominance” as a frequency >0.9, so that predominance of a type includes cryptic polymorphism, due to recurrent mutation.
• State 4: predominance of a new compatibility type, when a new lock/key type, more compatible with itself than with the initial type, is at a frequency >0.9.
• State 5: balanced suicidal polymorphism, when two types with higher cross-compatibility than self-compatibility (like lock_B/key_A and lock_A/key_B), and only these two, are at frequencies >0.1.
• State 6: this population state is peculiar in that it corresponds to all possible situations that are not described by the five other states. In practice, state 6 will describe mainly (i) populations harboring a neutral polymorphism (like lock_A/key_A and lock_B/key_A) and (ii) populations where a suicidal type (like lock_B/key_A) is at a frequency >0.9. However, we do not rule out that state 6 might include others types of situations, potentially interesting but unidentifiable on the basis of our current understanding of the system.

This typology being defined, the evolution of populations can be depicted as transitions between population states over generations. Starting from an initial condition, different simulations can lead to very different patterns, because of random mutation and drift. To describe general tendencies, we can compile the results of a sufficient number of simulation runs, which will allow us to estimate the probability of the different population states over generations. In the following sections, we use this approach to explore the evolution of compatibility types under various conditions. Preliminary analyses (not shown) have been performed using high mutation rates over a short number of generations, to moderate calculation time. These have allowed us to identify four important factors: (1) the length of the lock and key sequences (parameter S); (2) the mutation rate (parameter Mu); (3) the population size (parameter N_e); and (4) the upper limits for TE, FE, and MI (parameters TE_max, FE_max, and MI_max) that will condition the possible maintenance of uninfected individuals in the populations. The effect of these factors is now examined in detail.

Figure 8.— The evolution of population states over 100,000 generations, for two different S-values (S = 1 and S = 10), two different Mu-values (Mu = 10−6 and Mu = 10−5), and two different Ne-values (Ne = 103 and Ne = 104). The curves plot the frequency (more ...)

When S = 1 (Figure 8A), we observe that the probability of the initial configuration (state 3: predominance of the initial type) decreases very slowly over generations. Indeed, after 100,000 generations, predominance of the initial type is still observed in 90% of the simulations. Furthermore, we observe that among the remaining 10%, most populations have gone extinct due to fixation of a self-incompatible bacterium (state 2: population extinction). Thus, it appears that under these conditions, new compatibility types do not evolve. As illustrated in Figure 8B, things are clearly different when S = 10. Indeed, after 100,000 generations, <40% of the populations are still in the initial configuration, and the remaining 60% are in either state 4 (predominance of a new compatibility type) or state 6 (neutral polymorphism or predominance of partially suicidal Wolbachia). To apprehend the rationale behind these effects of S variations, one must distinguish two aspects of the differences observed between Figure 8A and 8B: (1) the slow vs. rapid decrease of state 3 (the initial configuration) and (2) the replacement of state 3 by state 2 (population extinction) vs. states 4 and 6. The explanation of difference 1 is the following: the mutation rate of the lock and key sequences is defined here as 10⁻⁶ per site; in other words, the overall mutation rate of the lock and key sequences is lower when the sequence is short (Figure 8A) than when it is long (Figure 8B); the initial configuration is lost faster in Figure 8A because the overall mutation rates of the lock and key sequences are higher. The explanation of difference 2 (population extinction vs. other states) is less straightforward. In Figure 8A, where S = 1, compatibility between a given lock and a given key can be only 0 or 1. In other words, if the intensity of mod is maximum (that is, if MI = 1), which is initially the case in these simulations, not a single viable egg is produced in incompatible crosses. Thus, if a suicidal type (with lock = 2 and key = 1) is fixed by drift, the host population goes extinct. In contrast, in Figure 8B, where S = 10, mutations of the lock and key sequences lead only to partial incompatibility. In other words, a suicidal type (partially suicidal) can reach fixation without leading the host population to extinction. In summary, the fact that populations go extinct in Figure 8A suggests that mutated lock sequences often reach fixation before a compatible key occurs by mutation. This readily leads the host population to extinction if S = 1, that is, if mutations of the lock sequence result in complete incompatibility.

We investigated in a similar way the effect of varying the n parameter (the number of possible states at each site in the lock and key sequences). The results are not presented graphically, as only minor quantitative effects were observed. Let us note that simply increasing the value of n tends to slow down the evolution of new compatibility types. This can be illustrated using the following example. Consider a population where a suicidal type (e.g., with lock = 2 and key = 1) is frequent. If n = 2, mutations of the key sequence will give rise to the appropriate key (that is key = 2) with probability 1. In contrast, if n = 10, the appropriate mutation will occur with probability 1/9 only. In other words, increasing the value of n reduces the probability of favorable mutations to occur in the key sequence, resulting in slower evolution of compatibility types. In all simulations presented here, n is arbitrarily set at 10, a value that appeared to be a reasonable compromise between realism and computation time constraints.

Three notable differences can be seen between Figure 8A and 8C: (1) in Figure 8C, <40% of populations are still in state 3 after 100,000 generations, as compared to 90% in Figure 8A; (2) in Figure 8C, >30% of the populations have gone extinct (state 2) after 100,000 generations, as compared to 10% in Figure 8A; and finally, (3) in Figure 8C, >20% of the populations harbor a new lock/key pair after 100,000 generations, as compared to 1% in Figure 8A. Difference 1 (that is, state 3 is less frequent in Figure 8C than in Figure 8A) illustrates that the initial lock/key type is lost faster when mutation rate is higher, as expected. It is notable that the rate of decrease of state 3 in Figure 8C is exactly the same as that observed in Figure 8B. This is consistent with the interpretation we gave when comparing Figure 8A with 8B: state 3 is lost faster when S is larger simply because the overall mutation rate of the lock and key sequences is increased. Difference 2 (that is, state 2 is more frequent in Figure 8C than in Figure 8A) suggests that extinction risk increases with mutation rate. This can be understood by remembering that extinction results from the fixation of a suicidal type (e.g., with lock = 2 and key = 1). Mutations giving rise to such suicidal types being neutral, their rate of fixation depends only on the mutation rate (Kimura 1983). Difference 3 (that is, state 4 is more frequent in Figure 8C than in Figure 8A) suggests that the evolution of compatibility types is facilitated by increased mutation rates, which might seem straightforward. A subtle aspect deserves, however, to be discussed: in Figure 8C, the ratio state 4/state 2 (new type/extinction) is much higher than that in Figure 8A. To understand why, one must remember that spreading of new type (e.g., with lock = 2 and key = 2) can occur only in populations where a suicidal type (e.g., with lock = 2 and key = 1) is sufficiently frequent (see Figure 5). When the mutation rate is increased (as in Figure 8C), mutations giving rise to this new type are more likely to occur before the suicidal type has reached fixation (that is, before populations have gone extinct).

Consider now Figure 8, B and D. Three notable differences can be seen: (1) in Figure 8D, only a tiny remnant of populations are still in state 3 after 50,000 generations, as compared to 60% in Figure 8B; (2) in Figure 8D, the frequency of state 4 reaches a plateau of 60% after 50,000 generations, while no plateau was reached in Figure 8B after 100,000 generations; and finally, (3) in Figure 8D, a small but significant and stable proportion of populations is in state 5 (balanced suicidal polymorphism), while this state was not observed in Figure 8B. Difference 1 (that is, the complete loss of state 3 in Figure 8D) suggests that state 3 (predominance of the initial compatibility type) is unstable in the long term: mutation rate in Figure 8D is high enough for its inescapable loss to be observed after 50,000 generations only. Difference 2 (that is, state 4 reaches a plateau in Figure 8D) suggests that even with high mutation rates, the proportion of populations harboring a new compatibility type can never reach 1. Difference 3 (that is, state 5 is more frequent in Figure 8D) suggests that high mutation rates facilitate the occurrence of balanced suicidal polymorphism. Overall, it is interesting to note that the conditions used in Figure 8D (Mu = 10⁻⁵ and S = 10) are such that 50,000 generations seem to be sufficient for the equilibrium distribution of population states to be observed. At equilibrium, populations are in state 4 (predominance of new type), state 5 (balanced suicidal polymorphism), or state 6 (neutral polymorphism or predominance of partially suicidal Wolbachia). This equilibrium is clearly dynamic: populations themselves are not stable, but the probability of transitions between states 4, 5, and 6 is stable. It is likely that a similar equilibrium would be observed with a mutation rate of 10⁻⁶, or even less, if populations were followed over a sufficiently large number of generations.

The results obtained with N_e = 10² are not presented graphically, as only minor quantitative effects were observed. Overall, reducing the population size from 10³ to 10² has very little effect on the evolution of Wolbachia compatibility types. The conclusions drawn with N_e = 10³ are retrieved: most importantly, the evolution of compatibility type occurs only when S is sufficiently large to prevent population extinction.

More interesting are the simulations performed with N_e = 10⁴ (presented in Figure 8, E–H). Most importantly, we note that when S = 1, increasing population size reduces risks of population extinction. In Figure 8, E and G, state 2 (population extinction) is observed in only a small fraction of all simulations, contrasting with that observed in Figure 8, A and C. To understand this result, one must recall (i) that population extinction can occur when S = 1 due to neutral fixation of a mutated lock (i.e., with lock = 2 and key = 1) and (ii) that the time separating the occurrence of a neutral allele and its fixation by drift is 4N_e on average (Kimura 1983). In other words, in large populations, the two mutation steps necessary for a new compatibility type to evolve often occur before extinction of the host population. This result is important with regard to the parameter space allowing the evolution of compatibility types: with large N_e, the evolution of compatibility types can occur even if S = 1. However, one should not conclude that large populations are protected from extinction in the long term: extinction being an absorbing state, it remains the ultimate fate of all populations if S = 1.

Figure 9.— See Figure 8 legend for details, except for the following: MImax = TEmax = FEmax = 0.9.

To simplify the analysis, note first that the right columns of Figures 8 and 9 (where S = 10) are strikingly similar. In other words, setting upper limits for MI, TE, and FE does not affect the evolution of compatibility types when S is sufficiently large. We therefore focus on comparing the left columns of Figures 8 and 9 (where S = 1).

Let us consider Figure 9A (Mu = 10⁻⁶, S = 1, N_e = 10³, MI_max = TE_max = FE_max = 0.9) and compare it with Figure 8A (same values for Mu, S, and N_e, but MI_max = TE_max = FE_max = 1). In these two figures, the initial population state (state 3: predominance of the initial type) decreases at the same rate. In both cases, after 100,000 generations, predominance of the initial type is still observed in 90% of the simulations. In Figure 8A, state 3 is replaced only by state 2 (population extinction), contrasting with Figure 9A, where state 3 is replaced by state 1 (infection loss). The interpretation is the following: in Figure 8A, suicidal variants (e.g., with lock = 2 and key = 1) can get fixed by drift, which causes population extinction; conversely, in Figure 9A (where MI_max = TE_max = FE_max = 0.9), the unstable equilibrium frequency (below which infection is deterministically lost) is increased by the presence of suicidal Wolbachia variants, so that the infection is always lost before the host population goes extinct. Thus, our model suggests that setting upper limits for MI, TE, and FE changes population extinction to infection loss. A similar conclusion is drawn from the comparison between Figures 8C and 9C (with Mu = 10⁻⁵): in Figure 8C, state 3 is mainly replaced by state 2 (population extinction), while it is replaced by state 1 (infection loss) in Figure 9C.

Shifting to Figure 9, E and G (compared with Figure 8, E and G), we can assess the effect of MI_max, TE_max, and FE_max in large populations (N_e = 10⁴). Most importantly, we note that risks of infection loss are only partially suppressed in large populations. In others words, even large populations are prone to infection loss when S = 1.

We have used here a simulation approach to investigate the evolution of compatibility types under various conditions. Being aware that the flurry of conditions and results is difficult to keep track of, we summarize here our main findings, considering in turn the effect of each parameter (listed in Table 1).

n is the number of possible states at each site of the lock and key sequences; it conditions the potential diversity of compatibility types. We found that the value of n has only minor effects on the evolution of compatibility types. We note, however, that increasing the value of n tends to slow down the evolution of new compatibility types, as it reduces the probability of favorable mutations in the key sequence.

S is the length of the lock and key sequences; it conditions the potential for partial compatibility. If S is small (i.e., S = 1), any mutation in the lock or key sequence leads to complete incompatibility. On the contrary, if S is large (i.e., S = 10), mutation can give rise to different, but partially compatible lock and key. We found that the value of S has profound implications on the evolution of compatibility types: under small S values, population extinction or infection loss is likely to occur before new compatibility types get fixed. On the contrary, if S is large, populations are not prone to extinction or infection loss, allowing the evolution of compatibility types.

Mu is the mutation rate used in the model. We found that higher mutation rates accelerate the evolution of compatibility types, as expected. We also observed that high mutation rates increase risks of population extinction or infection loss when S is small.

N_e is the population size; it conditions the intensity of genetic drift and the mean time from mutation to fixation of neutral alleles. Most importantly, we observed that increasing population size tends to reduce risks of population extinction and, to a lesser extent, of infection loss. However, it must be kept in mind that population extinction and infection loss remain the only ultimate fates of all populations when S is small, since these are stable and absorbing states.

MI_max, TE_max, and FE_max are the upper limits of MI, TE, and FE, respectively; these condition the potential for uninfected individuals to persist in infected populations. We observed that setting these limits below their absolute maximum (that is, <1) prevents populations from going extinct when S is small. However, this does not facilitate the evolution of compatibility types, since risks of population extinction are replaced by risks of infection loss. Overall, we conclude that MI_max, TE_max, and FE_max do not affect the conditions under which the evolution of compatibility types is plausible.

In summary, our analysis points to S as the major parameter: the evolution of compatibility types most readily occurs when the lock and key sequences are sufficiently long, that is, when the lock/key pairs change gradually rather than suddenly. Interestingly, partial compatibility between distinct Wolbachia strains has been observed in Drosophila: the wRi variant from Drosophila simulans is partially compatible with the wMel variant from D. melanogaster (Poinsot et al. 1998). Similarly, the wCer2 variant from Rhagoletis cerasi is partially compatible with wMel and wRi (Charlat et al. 2004; Riegler et al. 2004). Such empirical data support a model with S > 1.

In an alternative model, Dobson (2004) assumed that novel Wolbachia variants arising by mutation can coexist within a host individual together with the original infection. This analysis revealed interesting alternative routes for the evolution of compatibility types. In particular, it suggested that new types might evolve from an initial change in the resc function (mutation in the key sequence). One should note, however, that this remains plausible only if double infection is maintained over a sufficient number of generations. In other words, segregation from double to single infection (due to imperfect vertical transmission) has the potential to impede this particular process. Incorporating these aspects in the present model would allow one to assess the likelihood of this evolutionary trajectory.

In this article, we explored the evolution of compatibility under a wide variety of conditions. This analysis relies on heavy assumptions regarding the molecular mechanism of cytoplasmic incompatibility. Progress in this domain will accelerate now that the Wolbachia genome has been fully sequenced and analyzed (Wu et al. 2004). In turn, this should enrich the models and shed light on the evolutionary processes underlying the diversity of Wolbachia compatibility types.

We are deeply grateful to Thomas Pornin and Matthew Collette for their contribution to debugging. We also thank Frank Jiggins and two anonymous referees for critical reading and constructive comments.