Evolution of Multigene Families by Gene Duplication: A Haploid Model

Corresponding author: Hidenori Tachida, Department of Biology, Faculty of Science, Kyushu University, 33, Fukuoka, 812-8581, Japan., htachscb{at}mbox.nc.kyushu-u.ac.jp (E-mail).

Communicating editor: A. G. CLARK

	ABSTRACT

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

Evolution of multigene families by gene duplication and subsequent diversification is analyzed assuming a haploid model without interchromosomal crossing over. Chromosomes with more different genes are assumed to have higher fitness. Advantageous and deleterious mutations and duplication/deletion also affect the evolution, as in previous studies. In addition, negative selection on the total number of genes (copy number selection) is incorporated in the model. First, a Markov chain approximation is used to obtain formulas for the average numbers of different alleles, genes without pseudogene mutations, and pseudogenes assuming that mutation rates and duplication/deletion rates are all very small. Computer simulation shows that the approximation works well if the products of population size with mutation and duplication/deletion rates are all small compared to 1. However, as they become large, the approximation underestimates gene numbers, especially the number of pseudogenes. Based on the approximation, the following was found: (1) Gene redundancy measured by the average number of redundant genes decreases as advantageous selection becomes stronger. (2) The number of different genes can be approximately described by a linear pure-birth process and thus has a coefficient of variation around 1. (3) The birth rate is an increasing function of population size without copy number selection, but not necessarily so otherwise. (4) Copy number selection drastically decreases the number of pseudogenes. Available data of mutation rates and duplication/deletion rates suggest much faster increases of gene numbers than those observed in the evolution of currently existing multigene families. Various explanations for this discrepancy are discussed based on our approximate analysis.

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.

OHTA 1987A , OHTA 1987B , OHTA 1988A , OHTA 1988B theoretically studied evolution by gene duplication in haploid and diploid models using Monte Carlo simulation. She showed among others that positive selection is necessary to acquire genes with new functions and that the variance of the gene number is generally large. Also the ratio (R) of the number of genes with different functions to the number of pseudogenes was proposed to measure effects of positive selection, and an approximate formula to compute it was derived (OHTA 1987B ),

(1)

where v₊ and v_- are rates of mutation to a different gene with a new function and of mutation to a pseudogene, respectively, and u₊ is the fixation probability of an advantageous mutation introduced into a population. Later, CLARK 1994 and WALSH 1995 studied different aspects of two gene systems assuming slightly different models. However, there are questions still to be answered. First, because OHTA's results (OHTA 1987A , OHTA 1987B , OHTA 1988A , OHTA 1988B ) were mostly from simulation except for Equation 1 dependency of gene multiplication on parameters like population size, intensity of selection, mutation rate, and duplication/deletion rate was not so clear. As we do not have good estimates for most of these parameters, deriving analytical expressions describing gene multiplication processes is desirable to quantify speed of gene duplication and gene redundancy. Second, negative selection on the total number of genes per chromosome (copy number selection) was not considered in the previous works. A tendency for DNA loss is observed in Drosophila (PETROV et al. 1996 ), and gene silencing is observed when multiple copies are introduced into plants (HOLLICK et al. 1997 ). Thus, there is a possibility of the action of copy number selection. Third, the rate of gene duplication from a single gene and from multiple genes are considered to be different. However, this difference was ignored to simplify the model in OHTA's simulation.

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

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

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 ),

(2)

where 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 ),

(3)

3. Model III,

(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

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

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 n_l = ^k_i=1 m_i , 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,
   

   (5)

where F(p) = ^p₀ exp[-Ns_ax(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_iNg₊/2, where g₊ is the fixation probability of a deleterious mutant with a selection coefficient -s_d and is expressed as (KIMURA 1962 ),

(6)

4. The number of live genes in the ith allelic class decreases by deletion. The rate for this event in m_iNg_-/2, where g_- is expressed as,

(7)

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

6. The number of pseudogenes decreases by deletion. The rate of decrease is n_pNg_-/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 ^-1, 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 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.01s_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.

View larger version (20K): [in this window] [in a new window] [Download PPT slide]   Figure 1. The probability, p1(t), that the number of different alleles is 1 at t. All combinations of Nsa = 5, 10, 40, and sd = 0, 0.01sa are used for selection coefficients. Other parameters are ua = 0.1 and up = .

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 ₁(t) of n_t being 1, simulations are continued until p₁(t) - ₁(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.

View this table: [in this window] [in a new window]   Table 2. The ultimate probability, P1(), of the population being fixed by chromosomes with one copy

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.

Consider the following process: At the start, m is one. Transition from m = 1 to m = 2 occurs at the rate a = Ng₊/2 by gene duplication. For m > 1, the process moves to m - 1, m + 1, or 1. The transition to 1 is added because we want to know the average Q over all allelic classes, and for this purpose, the process moves to a new allelic class with probability one-half if a different allele is created from this allelic class. The transition to m - 1 occurs by deletion, advantageous mutation, or pseudogene mutation, and its rate is denoted by mb. Because we stay in the same allelic class in one-half of the cases when an advantageous mutation occurs, b = u_p + N(u_af₊/2 + g_-/2). The transition to m + 1 occurs by duplication and its rate is ma. Finally, the transition rate to 1 is mc with c = Nu_af₊/2. Let q_m(t) be the probability of the state being m at time t. Then, q_m(t)'s satisfy the following system of differential equations:

(8)

(9)

With the generating function h(z,t) = _m=1 q_mz^m , the system of differential equations is expressed by

(10)

where M_q = _m=1 mq_m is the mean number. Although the time-dependent solution of this equation is difficult to obtain because of q₁ and M_q in the right-hand side, we can compute the equilibrium solution h(z) by noting _m=1 q_m = 1 and q_m 1 for all n. The resulting solution is

(11)

where

(12)

(13)

The equilibrium distribution of m is easily obtained by expanding the generating function. The mean number M_q is also computed from h(z) and represented as

(14)

It is intuitively clear that increasing b or c decreases and increasing a increases M_q, but quantitatively how it changes is fairly complex, as expressed by (14).

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.

Next, we consider the number (k) of different alleles. Because a new allele is created from a redundant copy, its production rate is proportional to the number of redundant genes and u_aNf₊. Thus, in the second phase, k(t) is expected to be approximated by a pure-birth process with a rate of (M_q - q₁)u_aNf₊. Define k*(t) = E[k(t),k() > 1] as the expectation of k over the event that k eventually enters the second phase. The relationship between k*(t) and E[k(t)] is represented by

(15)

As noted above, because k(t) is a pure-birth process in the second phase, k*(t) takes the form

(16)

where = (M_q - q₁)u_aNf₊. Thus, we can write the expectation of k as

(17)

To check the accuracy of this expression based on the pure-birth process approximation, the values computed from it and E[k(t)] obtained from the Markov chain simulation were compared (Figure 2). To obtain C, the simulation data were fitted by the curve represented by (17). The simulation data are represented by open symbols, and the fitted curve is shown by solid lines in Figure 2. The agreement seems very good except for the very intial phase. Thus, in the second phase, k(t) evolves approximately as a pure-birth process defined as above. Note that the fitted curves are >1 at t = 0. This is probably because we started with a redundant copy (m₁ = 2). In this curve fit, we regard that a certain amount of new alleles are created at t = 0.

View larger version (19K): [in this window] [in a new window] [Download PPT slide]   Figure 2. 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 , Nsa = 5; , Nsa = 10; and , Nsa = 40. The curve fits are obtained assuming Equation 17. See text for more detail. Other parameters are ua = 0.1 and up = .

Finally, we consider the number of pseudogenes, n_p. The average rate of production of new pseudogenes per allelic class is again proportional to the number of redundant genes and represented by ß₁ = (M_q - q₁)u_p. Once a new pseudogene is created, its number can change by duplication/deletion. The rate of decrease is ß₂ = (g_- - g₊)N/2. Thus, the expected number of pseudogenes may take the following form,

(18)

However, recall that the estimated E[k(t)] by (17) is not 1 at t = 0. Because pseudogenes are produced from redundant genes as new alleles but with a different rate ß₁ per allelic class, we need to add the following term due to pseudogenes created at t = 0:

(19)

The term in parentheses is the number of new alleles created at t = 0 in the curve fit mentioned above, and the ratio of the number of pseudogenes to that of new alleles when they are created is ß₁/ (OHTA 1987B ). Combining (18) and (19), we obtain

(20)

Although the computation of this formula requires p₁(), if we want an expectation conditioned on the event of k > 1, it can be computed from the conditional expectation of k(t). Thus, this formula can be used to examine the relationship between conditional expectations of k and n_p. We again checked the accuracy of this equation with the values obtained by the Markov chain simulation, and the results are shown in Table 4. The agreement is generally very good. Without copy number selection (s_d = 0), ß₂ = 0, and if we note k*(t) - (1 - p₁()) = E[k(t)] - 1 [see (15)], Equation 20 is simplified to

(21)

which is equivalent to (1) obtained by OHTA 1987B . With the copy number selection, the following expression is obtained for ß₂t 1,

(22)

Thus, if that increases ß₂ linearly is large, E[n_p]/k* becomes small. This can be seen from Table 4, where even with weak copy number selection (s_d = s_a/100), the number of pseudogenes is drastically reduced.

	SIMULATION

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

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 ₁(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 ₁(t) from the two simulations agree well.

View larger version (26K): [in this window] [in a new window] [Download PPT slide]   Figure 3. 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 (nl), pseudogenes (np), and the proportion of populations with nt = 1 (p1) 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 = up = 0.0001, ua = 0.00001, sa = 0.1, sd = 0.

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. 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

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

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.

First, note that in the SSWMD limit, if we measure time in units of ^-1, the behavior of the process is characterized by the scaled parameters u^*_a = and u^*_p = , in addition to Ns_a and Ns_d. Thus, if we keep these parameters constant, the birth rate is proportional to . Let * = / and consider a simple case of no copy number selection (s_d = 0). In this case, * is expressed as

(23)

where

Here, d is proportional to both u^*_a and Nf₊. Two limiting cases are instructive. If d 1, Equation 23 reduces to

(24)

This shows that * increases linearly as d increases and that the rate is a decreasing function of u^*_p . Thus, if u^*_p is small, * quickly increases as d increases either by an increase of N or u^*_a . However, the increase is saturated, and for d 1,

(25)

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^-5, u_p = 10^-4, = 10^-4, 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.

View larger version (15K): [in this window] [in a new window] [Download PPT slide]   Figure 4. 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, ua = 0.1, up = ; Set 2, ua = 0.02, up = 0.2. Other parameters are sa = 0.1, sd = 0.002.

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 x 10^-5 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 x 10^-6 and 1.7 x 10^-4, respectively (SHAPIRA and FINNERTY 1986 ). In maize, the rate per copy is estimated to be in the range 10^-4 - 10^-3 (SUDPAK et al. 1993; EGGELSTON et al. 1995). Here, we tentatively assume = 10^-4, 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^-6, employing the bandmorph mutation rate 1.03 x 10^-6, estimated by MUKAI and COCKERHAM 1977 in Drosophila. With Ns_a = 40 and no copy number selection, the doubling time is 8 x 10⁴ generations. Because d* 1, the doubling time is shown to be roughly proportional to (u_aNf₊)^-1 from (24), and reducing u_a to one-tenth makes the time about 10 times greater. At any rate, gene 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 x 10^-8 for Adh, 2.8 x 10^-7 for Chs, and 2.7 x 10^-7 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 , OHTA 1987B 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 x 10⁶ (u_a = 10^-7) and 2 x 10⁵ (u_a = 10^-6) 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 Table 4 and Table 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 , 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 , 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.

	ACKNOWLEDGMENTS

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.

Manuscript received October 12, 1997; Accepted for publication May 13, 1998.

	LITERATURE CITED

TOP ABSTRACT MODEL MARKOV CHAIN APPROXIMATION SIMULATION DISCUSSION LITERATURE CITED

BROWN, C. J., C. F. AQUADRO, and W. W. ANDERSON, 1990  DNA sequence evolution of the amylase multigene family in Drosophila pseudoobscura.. Genetics 126:131-138[Abstract].

BUCKLER, E. S., IV, A. IPPOLITO, and T. P. HOLTSFORD, 1997  The evolution ribosomal DNA: divergent paralogues and phylogenetic implications. Genetics 145:821-832[Abstract].

CARROLL, S. B., 1995  Homeotic genes and the evolution of arthropods and chordates. Nature 376:479-485[Medline].

CLARK, A. G., 1994  Invasion and maintenance of a gene duplication. Proc. Natl. Acad. Sci. USA 91:2950-02954[Abstract/Free Full Text].

CLEGG, M. T., M. P. CUMMINGS, and M. L. DURMBIN, 1997  The evolution of plant nuclear genes. Proc. Natl. Acad. Sci. USA 94:7791-7798[Abstract/Free Full Text].

COX, D. R., and H. D. MILLER, 1965 The Theory of Stochastic Processes. Chapman and Hall, London.

CROW, J. F., and M. KIMURA, 1970 An Introduction to Population Genetic Theory. Harper & Row, New York.

DA LAGE, J. L., M. L. LEMEUNIER, M. L. CARIOU, and J. R. DAVID, 1992  Multiple amylase genes in Drosophila ananassae and related species. Genet. Res. 59:85-92[Medline].

EGGLESTON, W. B., M. ALLEMAN, and J. L. KERMICLE, 1995  Molecular organization and germinal instability of R-stippled maize. Genetics 141:347-360[Abstract].

FRYXELL, K. J., 1996  The coevolution of gene family trees. Trends Genet. 12:364-369[Medline].

GILLESPIE, J. H., 1983  Some properties of finite populations experiencing strong selection and weak mutation. Am. Nat. 121:691-708.

HOLLICK, J. B., J. E. DORWEILER, and V. L. CHANDLER, 1997  Paramutation and related allelic interactions. Trends Genet. 13:302-308[Medline].

IWABE, N., K. KUMA, and T. MIYATA, 1996  Evolution of gene families and relationship with organismal evolution: rapid divergence of tissue-specific genes in the early evolution of chordates. Mol. Biol. Evol. 13:483-493[Abstract].

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

KIMURA, M. and J. F. CROW, 1964  The number of alleles that can be maintained in a finite population. Genetics 49:725-738[Free Full Text].

LANGE, B. W., C. H. LANGLEY, and W. STEPHAN, 1990  Molecular evolution of Drosophila metallothionein genes. Genetics 126:921-932[Abstract].

LI, W. H., 1987  Models of nearly neutral mutations with particular implications for nonrandom usage of synonymous codons. J. Mol. Evol. 24:337-345[Medline].

LYCKEGAARD, E. M. S. and A. G. CLARK, 1989  Ribosomal DNA and stellate gene copy number variation on the Y chromosome of Drosophila melanogaster.. Proc. Natl. Acad. Sci. USA 86:1944-1984[Abstract/Free Full Text].

MUKAI, T. and C. C. COCKERHAM, 1977  Spontaneous mutation rates at enzyme loci in Drosophila melanogaster.. Proc. Natl. Acad. Sci. USA 74:2514-2517[Abstract/Free Full Text].

NEITZ, M. and J. NEITZ, 1995  Numbers and ratios of visual pigment genes for normal red-green color vision. Science 267:1013-1016[Abstract/Free Full Text].

OHNO, S., 1970 Evolution by Gene Duplication. Springer-Verlag, New York.

OHTA, T., 1983  Theoretical study on the accumulation of selfish DNA. Genet. Res. 41:1-15[Medline].

OHTA, T., 1987a  A model of evolution for accumulating genetic information. J. Theor. Biol. 124:199-211[Medline].

OHTA, T., 1987b  Simulating evolution by gene duplication. Genetics 115:207-213[Abstract/Free Full Text].

OHTA, T., 1988a  Further simulation studies on evolution by gene duplication. Evolution 42:375-386.

OHTA, T., 1988b  Multigene and supergene families. Oxf. Surv. Evol. Biol. 5:41-65.

OHTA, T., 1991  Multigene families and the evolution of complexity. J. Mol. Evol. 33:34-41[Medline].

PAYANT, V., S. ABUKASHAWA, M. SASSEVILLE, B. F. BENKEL, and D. A. HICKEY et al., 1988  Evolutionary conservation of the chromosomal configuration and regulation of amylase genes among eight species of the Drosophila melanogaster species subgroup. Mol. Biol. Evol. 5:560-567[Abstract].

PETROV, D. A., E. R. LOZOVSKAYA, and D. L. HARTL, 1996  High intrinsic rate of DNA loss in Drosophila. Nature 384:346-349[Medline].

SHAPIRA, S. K. and V. G. FINNERTY, 1986  The use of genetic complementation in the study of eukaryotic macromolecular evolution: rate of spontaneous gene duplication at two loci of Drosophila melanogaster.. J. Mol. Evol. 23:159-167[Medline].

SHIBATA, H. and T. YAMAZAKI, 1995  Molecular evolution of the duplicated Amy locus in Drosophila melanogaster species subgroup: concerted evolution only in coding region and excess of nonsynonymous substitutions in speciation. Genetics 141:223-236[Abstract].

STEPHAN, W., 1986  Recombination and the evolution of satellite DNA. Genet. Res. 47:167-174[Medline].

SUDUPAK, M. A., J. L. BENNETZEN, and S. H. HULLBERT, 1993  Unequal exchange and meiotic instability of disease-resistance genes in the Rp 1 region of maize. Genetics 133:119-125[Abstract].

TAKANO, T. S., S. KUSAKABE, A. KOGA, and T. MUKAI, 1989  Polymorphism for the number of tandemly multiplicated glycerol-3-phosphate dehydrogenase genes in Drosophila melanogaster.. Proc. Natl. Acad. Sci. USA 86:5000-5004[Abstract/Free Full Text].

THEISSEN, G., J. T. KIM, and H. SAEDLER, 1996  Classification and phylogeny of the MADS-box multigene family suggest defined roles of MADS-box gene subfamilies in the morphological evolution of eukaryotes. J. Mol. Evol. 43:484-516[Medline].

WALSH, J. B., 1985  Interaction of selection and biased gene conversion in a multigene family. Proc. Natl. Acad. Sci. USA 82:153-157[Abstract/Free Full Text].

WALSH, J. B., 1995  How often do duplicated genes evolve new functions? Genetics 139:421-428[Abstract].

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar