Evolution of Multigene Families by Gene Duplication: A Haploid Model

Abstract

GENE duplication with subsequent diversification is considered to have played very important roles in the organismal evolution (Ohno  1970, 1988b, 1991). If a gene is duplicated, the selective constraint becomes less for the extra copy, and it can evolve to have a (slightly) different function, while the original function of the gene is kept in the other copy. Thus, gene duplication with subsequent diversification is one of the simplest ways to acquire a new function and is thought to have been employed many times during evolution. For example, Iwabe et al. (1996) suggested that gene duplications contributed to the compartmentalization of a cell and formation of multicellular organisms by creating genes expressed tissue specifically. Also, homeobox genes (Carroll 1995) and MCM1, AGAMOUS, DEFICIENS, and SRF (MADS)-box genes (Thiessen et al. 1996) were apparently created by ancient gene duplications. Gene duplications are not restricted to ancient times. Changes of gene number are observed in closely related species in many taxa, and even examples of polymorphisms with regard to copy number within species are known (e.g., Lyckegaard and Clark 1989; Takano et al. 1989; Lange et al. 1990; Neitz and Neitz 1995). Thus, gene duplications are still ongoing evolutionary processes. How such gene duplications and formation of multigene families occur at the population level is an interesting evolutionary problem.

In the present article, as an extension of Ohta's work, a haploid model of gene duplication is analyzed paying attention to the change of gene number. Copy number selection is explicitly incorporated into the model, and the duplication rate from a single gene was assumed to be 0. First, the model is described. Then, a Markov chain approximation of the model is derived, and its behavior is analyzed assuming low mutation and duplication/deletion rates. Results of extensive computer simulation to check the accuracy of the approximate analysis will be described. Approximate formulas for the rate of increase of the number of different genes and gene redundancy were obtained as a function of population size, mutation rates, duplication/deletion rates, and selection coefficients. Copy number selection was shown to be very effective in reducing the number of pseudogenes. Relationships of various haploid models also will be discussed.

MODEL

Definitions of symbols used are summarized in Table 1. A random-mating haploid population consisting of N chromosomes is assumed. Generations are discrete, and all chromosomes have two identical genes at generation zero. Each chromosome undergoes mutation, gene duplication/deletion, and random sampling with selection, in this order, in one generation.

We call a gene that never has had pseudogene mutation in its descent a live gene. Each live gene mutates to a new different allele or to a pseudogene with rates u_a or u_p, respectively. Here, following the usage of Ohta (1988b), we use the word alleles to denote gene types, although alleles may not be at the same locus. Each new allele is assumed to be different from any other alleles preexisting in the population (Kimura and Crow 1964). Because fitness of a chromosome is a monotone increasing function of the number of different alleles, as described later, u_a is the advantageous mutation rate. Once a gene becomes a pseudogene, it stays a pseudogene.

Next, if a chromosome has more than one gene, duplication

or deletion occurs with a rate γ/2 per gene copy, respectively. If a chromosome has just one gene, no change occurs. Here, we are considering sister-chromatid unequal crossing over as a genetic mechanism for copy number change. If the mechanism for copy number change is duplicative transposition, copy number change can occur even in chromosomes with a single copy.

Finally, N chromosomes of the next generation are sampled from the N chromosomes of the present generation with probabilities proportional to fitness of each chromosome. Let k and n_t be the numbers of different alleles and all genes in a chromosome, respectively. The fitness of the chromosome is zero if k = 0. Otherwise, we consider the following three fitness functions, w(k,n_t).

1. Model I (Ohta 1987b),

   $w(k,nt)={exp(−sdnt)(k≥k¯)exp[−sa(k¯−k)−sdnt](k<k¯)},$

   (2)

   where $k¯$ is the average number of different alleles per chromosome in the population, and s_a and s_d(s_a, s_d > 0) are positive and negative selection coefficients, respectively.
2. Model II (Walsh 1995),

   $w(k,nt)=exp[sak−sdnt].$

   (3)
3. Model III,

   $w(k,nt)={exp[sa(k−k¯)−sdnt](k>k¯)exp(−sdnt)(k≤k¯)}.$

   (4)

   Models I and II are modifications of those by Ohta (1987b) and Walsh (1995), respectively, with addition of deleterious effects determined by the total number of genes (copy number selection). As will be noted in the next section, models I, II, and III correspond to the recessive, semidominant, and dominant advantageous mutation models, respectively, in regard to the increase of the number of different alleles. We mainly show results for Model I. Differences among the models will be discussed where relevant.

Random sampling with selection completes one generation. This cycle is repeated many times. At generation zero, the population starts with all chromosomes having two live genes (k = 1, n_l = 2, n_p = 0), and all genes are identical as mentioned before.

Homologous equal crossing over does not occur in the present haploid model. Thus, this is a single locus infinite-allele model. For large populations, the diffusion approximation can be applied, and, in this approximation, relevant parameters are Nu_a, Nu_p, Nγ, Ns_a, and Ns_d measuring time in units of N (Crow and Kimura 1970).

Because the evolution of the model seems very complex, we first analyze a Markov chain whose behavior approximates that of our haploid model under a certain condition. The validity of the approximation will be examined by computer simulation later.

MARKOV CHAIN APPROXIMATION

We assume that Nu_a, Nu_p, and Nγ are all very small compared to 1. Also, we assume that positive selection is fairly strong, so that the number of different alleles never decreases. Under these conditions, the population is expected to be mostly monomorphic and thus represented by one chromosome type. A mutant chromosome type that is to replace the currently dominant chromosome type in the population is expected to be fixed quickly. For characterization of a chromosome, let k be the number of different alleles in the chromosome as in the previous section. The chromosome is characterized by k numbers, m₁,⋯,m_k, each representing the number of live genes with the same allelic state and the number of pseudogenes, n_p. Thus, in this approximation, the population state is characterized by a vector (m₁,⋯,m_k,n_p). The number of live genes, n_l, is expressed as $n1=Σi=1kmi$⁠, and the total number, n_t, is expressed as n_t = n_l + n_p. Occasionally, transitions among states occur. Because we assume Nu_a, Nu_p, Nγ ≪ 1, one of the following types of changes occurs in a chromosome, and the chromosome is fixed in the population with the probabilities specified below:

1. A live gene in the ith class with m_i > 1 has an advantageous mutation. The rate of this event is u_am_iNf₊, where f₊ is the fixation probability of a mutant chromosome with one more different allele than other chromosomes in the population. Depending on the models of selection, f₊ has a different form and is approximately expressed as,

   $f+={(4sa)∕(πN)(Model I, Ohta1987b)2sa∕(1−exp(−2Nsa))(Model II, Walsh1995)F(1∕N)∕F(1)(Model III)}$

   (5)

   where $F(p)=∫0pexp[−Nsax(2−x)]dx$ (see Crow and Kimura 1970). Note that the fixation probabilities in models I, II, and III correspond to those of recessive, semidominant, and dominant advantageous mutations, respectively. The derivation for Model III can be done in the same way as that in Model I by Ohta (1987b).

2. A live gene in the ith class with m_i > 1 becomes a pseudogene. Because such change is selectively neutral with m_i > 1, the rate for the ith class genes to become a pseudogene is m_iu_p.

3. The number of live genes of the ith class increases by gene duplication. The rate for this event is m_iγNg₊/2, where g₊ is the fixation probability of a deleterious mutant with a selection coefficient −s_d and is expressed as (Kimura 1962),

   $g+=exp(2sd)−1exp(2Nsd)−1.$

   (6)
4. The number of live genes in the ith allelic class decreases by deletion. The rate for this event in m_iγNg₋/2, where g₋ is expressed as,

   $g−=1−exp(−2sd)1−exp(−2Nsd).$

   (7)

5. The number of pseudogenes increases by gene duplication. The rate of increase is n_pγNg₊/2.

6. The number of pseudogenes decreases by deletion. The rate of decrease is n_pγNg₋/2.

The above transition probabilities determine a Markov chain. Because we assumed strong selection (SS) and weak mutation, duplication/deletion (WMD), we call the result of this approximation the SSWMD limit. This type of approximation has been widely used to analyze complex genetic systems (Gillespie 1983; Walsh 1985; Li 1987). In this approximation, further reduction of the number of relevant parameters occurs. If we measure time in units of γ⁻¹, u_a/γ and u_p/γ, in addition to Ns_a and Ns_d, determine the behavior of the process. The process is still a multidimensional Markov chain and is difficult to analyze directly. So we take two approaches. One is computer simulation of the Markov chain and another is further approximations for studying some

Figure 1.

Open in new tabDownload slide

The probability, p₁(t), that the number of different alleles is 1 at t. All combinations of Ns_a = 5, 10, 40, and s_d = 0, 0.01 s_a are used for selection coefficients. Other parameters are u_a = 0.1γ and u_p = γ.

aspects of the process. For each parameter set, 10,000 replications were made in the simulation.

Because k is an increasing function of time with strong positive selection, k is always > 1 once k becomes more than 1. We call the phases with k = 1 and k > 1 the first and second phases, respectively. In the first phase, the total number of genes, n_t, can become 1, and then no change of the gene number occurs thereafter. For this reason, the number of pseudogenes, n_p, affects the evolution of k and n_t. On the other hand, in the second phase, n_l is always larger than 1, and change of the number of genes never stops. The evolution of k and n_l does not depend on n_p. Therefore, these two phases are analyzed separately.

In the first phase, the population starts with k = 1, n_t = 2, and we are interested in when the population enters the second phase (k > 1). Let p₁(t) be the probability that k is 1 at time t. Because evolution in this phase is influenced by n_p, p₁(t) is difficult to compute. We examined p₁(t) by simulation of the Markov chain. Figure 1 shows the time-dependent behavior of p₁(t). Solid and broken lines represent results for cases without (s_d = 0) and with (s_d = 0.01 s_a) copy number selection, respectively, and time is measured in units of 1/γ generations. The mutation rates in this example are 10u_a = u_p = γ. The probability p₁(t) decreases quickly in the first few 1/γ generations, and then the change becomes small. In the later generations, n_t is 1 in most populations with k = 1, and no changes occur in these populations. That is why p₁(t) does not change in the later generations. Positive selection decreases the probability, but very weak copy number selection can increase the probability. The effect of copy number selection is stronger for the cases with stronger positive selection.

Ultimately, the population moves to either the state with n_t = 1 or that with k > 1. The ultimate probability, p₁(∞), was estimated by simulation. Because p₁(t) is a monotone decreasing function of t and bounded from the lower side by the probability $p~1(t)$ of n_t being 1, simulations are continued until p₁(t) − $p~1(t)$ < 0.01 so that we can estimate p₁(∞) with a maximum error of 0.01. The results are shown in Table 2. The ratio u_a/γ and effects of selection are important in determining p₁(∞). For small u_a/γ, most populations become ultimately fixed by chromosomes with a single copy. The ratio u_p/u_a has little effect probably because the duplication/deletion force is strong. As u_a/γ becomes larger, p₁(∞) decreases and the effects of pseudogene mutation become larger. The effect of copy number selection is always significant even though s_d/s_a = 0.01 in the examples.

In the second phase, we are interested in the evolution of k, n_l, and n_p. We first investigate the ratio Q = n_l/k, the average number of genes with the same allelic state and a measure of strict gene redundancy. Because the duplication/deletion process does not stop with n_t > 1, the evolution of the number of genes with the ith allelic state, m_i (we call this the ith allelic class) can be decoupled from n_p or other m_j's in this approximation. Here, we first concentrate on the evolution of one allelic class and denote the number of live genes in the class by m.

Because M_q is the mean number of live genes in one allelic class, its value might be close to the ratio Q = n_l/k, although M_q computed above is the equilibrium value. We can check how well M_q can approximate Q by simulation. Because Q is 1 in populations with n_t = 1, we computed Q conditioned on n_t > 1 from the simulation at t = 10/γ (Table 3). As shown in Table 3, Q estimated from the simulation and M_q computed from (14) agree well. Thus, it seems that the distribution of m is close to equilibrium at t = 10/γ. As Ns_a or u_a becomes larger, the redundancy measured by M_q becomes less because both b and c increase. Increasing u_p increases b and thus decreases M_q. Introduction of the copy number selection also decreases M_q. Only when u_a/γ is very small, is high genetic redundancy observed, but otherwise the genetic redundancy is very small.

Figure 2.

Open in new tabDownload slide

Comparison of the number of different alleles obtained by Markov chain simulation and that expected under a linear pure-birth process. Simulation results are represented by ○, Ns_a = 5; □, Ns_a = 10; and ♦, Ns_a = 40. The curve fits are obtained assuming Equation 17. See text for more detail. Other parameters are u_a = 0.1γ and u_p = γ.

SIMULATION

To examine the validity and limitation of the SSWMD limit obtained by the Markov chain approximation and also to investigate parameter ranges outside that assumed in the approximation, we carried out computer simulation closely following the model. We call this type of simulation the Wright-Fisher simulation. In the simulation, N chromosomes are prepared. In each chromosome, k, m₁, …, m_k, n_p and its fitness are stored. Numbers of mutations and duplication/deletion events are determined by generating Poisson random numbers with means computed by multiplying the number of genes and respective rates. Genes undergoing mutations or duplication/deletion are chosen randomly. After these,

the fitness of each chromosome is determined according to (2), (3), or (4) and is divided by the maximum fitness of the current population. The next generation is constructed in the following way: First, choose a chromosome randomly. Next, draw a uniform random number in [0,1) and, if the fitness of the chosen chromosome is less than the random number, add it to the population of the next generation. This is repeated until the chromosome number becomes N. For each parameter set, 1600 replications are made. We present the simulation results in comparison with the SSWMD limit.

First, we examine time courses of the average k, n_l, n_p and the proportion $p~1(t)$ of populations with n_t = 1 when Nγ, Nu_p, and Nu_a are small and agreement between the Markov chain simulation and Wright-Fisher simulation is expected. Some examples are shown in Figure 3. The values computed by the Markov chain simulation are represented by solid lines and those by the Wright-Fisher simulation are represented by dotted lines. Parameters used are γ = u_p = 0.001, u_a = 0.0001, s_a = 0.1, s_d = 0, and N = 100 (Figure 3A) and N = 400 (Figure 3B). With N = 100, the Markov chain approximation is very close to the Wright-Fisher simulation. For this parameter set, Nγ = Nu_p = 0.01 and Nu_a = 0.001. However, with N = 400, the Markov chain approximation underestimates expected values of k, n_p, n_l by 10 to 20%, especially in the later generations, although they seem to increase exponenetially. The largest discrepancy is found in n_p. As the number of genes per chromosome increases, the population becomes polymorphic because product parameters such as n_lNγ become large. This may be the reason for the discrepancies. Estimates of $p~1(t)$ from the two simulations agree well.

To examine the numbers of genes in a wider parameter range, more simulations were carried out, and some of the results are shown in Table 5 (without copy number selection) and Table 6 (with copy number selection). To clarify the effect of increasing γ, u_a, and u_p, cases with γ = 0.0001 and 0.001 are shown along with the SSWMD limit keeping the u_a/γ and u_p/γ constant. Means and variances of gene numbers at t = 10/γ are shown. With u_a/γ = 0.02 and N = 100, values for γ = 0.001 and 0.0001 are quite close to the SSWMD limit. As population size increases (N = 400) or u_a/γ approaches 0.1, the values with a large γ (0.001) deviate from the SSWMD limit. Especially, the number of pseudogenes becomes large when γ is large, decreasing the ratio of the number of different alleles to those of pseudogenes. With larger N or u_a, populations become polymorphic with regard to gene numbers, thus breaking the assumption of the Markov chain approximation. The increase of pseudogenes is probably due to competition among advantageous alleles and resulting reduction of the probability of fixation of new different alleles. The trend becomes more pronounced when both N and u_a are larger (u_a/γ = 0.1, N = 400). In this parameter set, the numbers of different alleles and live genes are also larger than the SSWMD limit. This is especially true when there is copy number selection. In the Markov chain approximation, fixations of new alleles are assumed to occur after fixations of redundant copies. However, in polymorphic populations, new, different alleles may appear before fixation of redundant copies. This may increase the fixation rate of new genes, especially when there is copy number selection and fixations of redundant copies are difficult to make occur.

The number of different alleles increases as population size increases without copy number selection. However, this is not necessarily true with copy number selection. With u_a/γ = 0.02, the average numbers of different alleles are almost the same for N = 100 and 400. The variance of k is generally very large and in the order of the square of the mean when the mean is large. As u_a increases, keeping u_a/u_p constant, E[k] always increases.

Thus far, we investigated only cases with fairly strong selection (Ns_a = 10, 40), where one of the assumptions of the SSWMD limit is satisfied. We also investigated when selection is weak, and the result is shown in Table 7.

Figure 3.

Open in new tabDownload slide

Comparison of the expected numbers of genes obtained by the Markov chain simulation and Wright-Fisher simulation. Two population sizes are used. A, N = 100; B, N = 400. The average numbers of different alleles (k), live genes (n_l), pseudogenes (n_p), and the proportion of populations with n_t = 1 (p₁) are plotted against generation. Values obtained by the Markov chain simulation (M) are represented by closed symbols with solid lines and those by the Wright-Fisher simulation (W) are represented by open symbols with broken lines. Other parameters are γ = u_p = 0.0001, u_a = 0.00001, s_a = 0.1, s_d = 0.

The result of reducing γ to 0.0001 was indistinguishable from that of γ = 0.001 and is not shown. The average number of different alleles is very close to 1, although pseudogenes are created. Thus, for smaller populations, new alleles are difficult to evolve. Copy number selection has very little effect on E[k] and slight effects on E[n_p]. One notable difference from the strong selection case is that E[k] decreases as u_a increases, keeping u_a/u_p constant. The effect of increasing u_p overrides that of increasing u_a if Ns_a is of the order of 1. Also note that the variance of n_p is very large. Some populations are considered to have a very large number of pseudogenes, whereas n_t = 1 and n_p = 0 in other populations.

DISCUSSION

In the present article, evolution of multigene families by gene duplication was analyzed assuming a haploid population and no interchromosomal recombination. In such a population, redundant copies are created first, and then one of them evolves to become a new gene by adaptive mutation. This system differs from those in which selection directly operates on the number of gene copies (Ohta 1983; Stephan 1986) in that acquiring adaptive genes requires the presence of redundant genes (Ohta 1987a). To analyze this system, we introduced the Markov chain approximation for small mutation and duplication/deletion rates under which it is assumed that each population is fixed by one chromosome type and that evolution occurs by occasional jumps among chromosome types (the SSWMD limit). By this approximation, the extent of gene redundancy was calculated, and the evolution of the average number of different alleles was characterized. Effects of polymorphism expected under higher mutation and duplication/deletion rates were analyzed using computer simulation. The most significant departure from the SSWMD limit with the same ratios of the rates is the increase of the number of pseudogenes. Also, polymorphism elevates the number of different alleles compared to that computed by the SSWMD limit, especially when there is negative selection on the copy number. With these caveats in mind, we consider some of the questions concerning the evolution of multigene families based on the SSWMD limit.

The increase rate with d = ∞ is γ/2 per generation as expected, but the approach to it is affected by the pseudogene mutation rate. Between these two extreme cases, α* was numerically shown to be an increasing function of d in the range 0.01 ≤ u_p/γ ≤ 10 (data not shown). Thus, for example, the rate of increase of gene number is larger in larger populations with s_d = 0. However, this is not always true when there is copy number selection (s_d ≠ 0), as shown in Figure 4. In these cases, as N becomes larger, α* first increases and then decreases. This is because the effect of copy number selection also becomes larger in larger populations. However, as noted in the previous section, the SSWMD limit underestimates the number of different alleles with copy number selection as the change rates (u_a, u_p, and u_γ) become larger. Thus, the decrease of the rate as N increases is not as pronounced, as Figure 4 shows in more realistic situations. Indeed, when u_a = 10⁻⁵, u_p = 10⁻⁴, γ = 10⁻⁴, s_a = 0.1, and s_d = 0.002, corresponding to set 1 in Figure 4, the average numbers of different alleles at the 10/γth generation were 2.048 ± 0.0103 for N = 200 and 1.896 ± 0.081 for N = 800. In summary, the number of different alleles increases as population size becomes large without copy number selection, but this is not necessarily so with copy number selection.

Next, we examine the speed of gene duplication. We measure it by the doubling time of E[k]. It is expressed by log2/[(M_q − q₁)u_aNf₊] and depends on u_a, u_p, γ, Ns_a, and Ns_d. We do not have reliable estimates for most of the parameters, but some rough calculation under

this model is possible. For u_p, we may use the null mutation rate measured in Drosophila by Mukai and Cockerham (1977). Averaging over four loci encoding various enzymes, they obtained 1.03 × 10⁻⁵ for the null mutation rate per loci. Estimates of unequal crossing-over rates vary widely. For example, in Drosophila, the duplication rates from a single gene at the ma and ry loci are estimated to be 2.7 × 10⁻⁶ and 1.7 × 10⁻⁴, respectively (Shapira and Finnerty 1986). In maize, the rate per copy is estimated to be in the range 10⁻⁴ − 10⁻³ (Sudpak et al. 1993; Eggelston et al. 1995). Here, we tentatively assume γ = 10⁻⁴, taking an intermediate value. Currently, we do not have any estimates for the advantageous mutation rate u_a or selection parameters s_a, s_d, but to get a rough idea, let us assume u_a = 10⁻⁶, employing the bandmorph mutation rate 1.03 × 10⁻⁶, estimated by Mukai and Cockerham (1977) in Drosophila. With Ns_a = 40 and no copy number selection, the doubling time is 8 × 10⁴ generations. Because d* ≪ 1, the doubling time is shown to be roughly proportional to (u_aNf₊γ)⁻¹ from (24), and reducing u_a to one-tenth makes the time about 10 times greater. At any rate, gene

Figure 4.

Open in new tabDownload slide

Examples of the rate, α, of increase of the number of different alleles with copy number selection as a function of population size. α/γ is computed from the SSWMD limit as a function of population size. Two examples are shown. Set 1, u_a = 0.1γ, u_p = γ; Set 2, u_a = 0.02γ, u_p = 0.2γ. Other parameters are s_a = 0.1, s_d = 0.002.

duplications can occur very quickly under the present model from the view point of the evolutionary time scale. However, according to Fryxell (1996), the doubling time is of the order of 10 to 100 million years in a particular lineage of a gene family. Clegg et al. (1997) surveyed three multigene families in plants and estimated new gene recruitment rates to be 2.9 × 10⁻⁸ for Adh, 2.8 × 10⁻⁷ for Chs, and 2.7 × 10⁻⁷ for rbcS. Because explosive gene multiplications are not observed in the evolution of multigene families, Fryxell (1996) proposed that simultaneous coevolution of two or more gene families might be necessary for gene multiplication in respective families, which will diminish duplication rates significantly. However, there are other explanations. For example, Ns_a might not be constant. If Ns_a becomes close to 1 or less, the number of different alleles stays near 1 as shown by Ohta (1987a,b) and by the present study. Another explanation is copy number selection. If we introduce weak copy number selection (s_d = 0.01s_a), the doubling times in the above examples become 2 × 10⁶ (u_a = 10⁻⁷) and 2 × 10⁵ (u_a = 10⁻⁶) generations, making the times much longer. Thus, copy number selection might have contributed to retarding the duplication rates. Copy number selection might also explain the general observation of larger multiplicity of genes in plants and mammals compared to those in, for example, Drosophila and lower eukaryotes. The numbers of rDNA genes in plants are very large compared to those in Drosophila and seem to contain many pseudogenes (Buckler et al. 1997). The number of pseudogenes is expected to be very small under strong copy number selection (see Tables 4 and 6).

We usually observe large variation in gene numbers in one multigene family among related species. For example, in Drosophila amylase genes, the numbers of genes are two (Payant et al. 1988; Shibata and Yamazaki 1995), one to three (Brown et al. 1990), and at least four (Da  Lage et al. 1992) in D. melanogaster, D. pseudoobscura, and D. ananassae, respectively. Such variation can be explained by difference in parameters such as s_a and N. However, even with the same parameter set, Clark (1994) showed that multiple equilibrium points concerning gene configuration are possible in two gene systems. Moreover, the variance of gene number is large, an emphasized by Ohta (1987b). As shown here, this is because the increase of gene number behaves like a pure linear-birth process. The coefficient of variation of such a process is asymptotically expressed as $1∕n0$⁠, where n₀ is the initial number (see p. 159 of Cox and Miller 1965). In our case, a further increase of the variance occurs because the gene number remains 1 in a significant proportion of lineages. Thus, we should expect a very large variance of gene numbers in a multigene family.

In the present article, as a fitness function, Ohta's model (Model I; Ohta 1987b) was mostly used. In the Markov chain approximation, the difference in fitness function contributes only through difference in the fixation probability of a chromosome with one more different allele than other chromosomes in the population. As noted, models I, II, and III correspond to recessive, semidominant, and dominant advantageous mutation models in this regard. Thus, the increase in the number of different alleles is more difficult in this order because a mutant chromosome appears singly at first. However, if we consider a decrease of the number of different alleles by chance, models I, II, and III correspond to dominant, semidominant, and recessive deleterious mutant models. Therefore, although it is most difficult to increase the number of different alleles in Model I, it is also most difficult to decrease it. Dependency on parameters also differes among the models. The fixation rate d = u_aNf₊/γ of different alleles is proportional to $Nsa$⁠, in Model I and Ns_a in Model II (Walsh's model; Walsh 1995). Thus, with the same Ns_a > 1, the number of different alleles becomes larger under Model II. For d ≪ 1, α* is proportional to d as shown in (24). Indeed, for Nu_a = 0.04, Nu_p = Nγ = 0.4, Ns_a = 40, and s_d = 0, the Wright-Fisher-type simulation shows that the average number of different alleles at t = 25N was 150 in Model II, whereas it was 6 in Model I. Thus, the selection in Model I is the least efficient in increasing gene number.

We assumed that no duplication/deletion occurs if the total number, n_t, of genes becomes 1 because the duplication rate is much smaller if there is only a single copy. We divided the process into the first and second phases, depending on whether k is 1 or more than 1. In the second phase, the evolution is the same as that in Ohta's model (Ohta 1987b). The analysis of the first phase shows that even if Ns_a is fairly large, nearly half or more of populations starting from chromosomes with two copies ends up in single gene chromosomes unless u_a/γ is not small (Table 2). Thus, we expect to find lineages with only a single gene even if gene number is large in other lineages.

From our approximate analysis, the gene duplication processes in a haploid model without interchromosomal crossing over were characterized, and dependencies of the speed of gene duplication, genetic redundancy, and the number of pseudogenes on various parameters are clarified for strong selection cases. With the accumulation of data on gene families, various discussions on the mechanisms of the evolution of multigene families have been made (e.g., Fryxell 1996; Clegg et al. 1997). Our results will provide a quantitative framework for such discussions by incorporating gene duplication and subsequent divergence into the model. However, higher organisms are mostly diploid, and interchromosomal crossing over does occur. The diploidy and interchromosomal crossing over not included here are known to accelerate the increase of gene number (Ohta 1988a). Incorporating them into the model and analyzing the model in more detail would be necessary for a better understanding of the evolution of multigene families.

Acknowledgement

We thank T. Ohta and M. Iizuka, A. Clark, and an anonymous reviewer for their helpful comments on earlier drafts of the manuscript. This research was partially supported by grants-in-aid to H.T. from the Ministry of Education, Science and Culture of Japan.

Footnotes

LITERATURE CITED