Abstract
Transspecies polymorphism, meaning the presence of alleles in different species that are more similar to each other than they are to alleles in the same species, has been found at loci associated with vegetative incompatibility in filamentous fungi. If individuals differ at one or more of these loci (termed het for heterokaryon), they cannot form stable heterokaryons after vegetative fusion. At the hetc locus in Neurospora crassa and related species there is clear evidence of transspecies polymorphism: three alleles have persisted for ∼30 million years. We analyze a population genetic model of multilocus vegetative incompatibility and find the conditions under which transspecies polymorphism will occur. In the model, several unlinked loci determine the vegetative compatibility group (VCG) of an individual. Individuals of different VCGs fail to form productive heterokaryons, while those of the same VCG form viable heterokaryons. However, viable heterokaryon formation between individuals of the same VCG results in a loss in fitness, presumably via transfer of infectious agents by hyphal fusion or exploitation by aggressive genotypes. The result is a form of balancing selection on all loci affecting an individual's VCG. We analyze this model by making use of a Markov chain/strong selection, weak mutation (SSWM) approximation. We find that transspecies polymorphism of the type that has been found at the hetc locus is expected to occur only when the appearance of new incompatibility alleles is strongly constrained, because the rate of mutation to such alleles is very low, because the number of possible incompatibility alleles at each locus is restricted, or because the number of incompatibility loci is limited.
DISTINGUISHING self from nonself is essential for sexual reproduction, for defense against pathogen invasion, and for the maintenance of individuality. Population genetics theory has been applied successfully to model the evolution of two such recognition systems, the major histocompatibility complex (MHC) of vertebrates and selfincompatibility (SI) loci in plants (Wright 1939; Takahata 1990). The MHC is an array of linked loci that are important in pathogen recognition (Bjorkman and Parham 1990). Sexual reproduction in many plant species is mediated by gametophytic or sporophytic SI loci that elicit recognition and rejection of pollen from individuals with the same or similar genotypes (Clarke and Newbigin 1993).
Although MHC and SI loci encode different gene products, they share common evolutionary features. Numerous alleles at both kinds of loci show evidence of longterm persistence; alleles from one species are often more closely related to alleles in a closely related species than to other alleles in the same species (Figueroaet al. 1988; Ioergeret al. 1990), a phenomenon referred to as “transspecies polymorphism” (Klein 1980; Arden and Klein 1982; Kleinet al. 1998). The observation of transspecies polymorphism, along with the relatively even distribution of allele frequencies found within populations, supports the idea that these loci are under some form of balancing selection, either because of overdominance in fitness or frequencydependent selection favoring rarer alleles or genotypes. There is a welldeveloped population genetics theory of balancing selection at single loci that predicts the occurrence of transspecies polymorphism whenever selection is strong relative to genetic drift (Takahata 1990; Vekemans and Slatkin 1994). Although interactions are known to occur among MHC and SI loci, singlelocus theory provides a reasonable first approximation and has led to a greater understanding of how balancing selection, genetic drift, and mutation together determine allele numbers and allele frequencies at these loci. Models of balancing selection at two closely linked loci in finite populations are more complicated but lead to predictions that are similar to those from onelocus models (Kelly and Wade 2000; Slatkin 2000).
Transspecies polymorphism has also been detected at vegetative incompatibility loci in filamentous fungi. Nonself recognition is thought to be important during vegetative growth. As filamentous fungi grow, fusions between hyphae occur (a process called anastomosis), which yields a network of interconnected hyphae, or mycelium, that makes up the fungal individual (Glasset al. 2000). Filamentous fungi also possess the remarkable attribute of being able to undergo hyphal fusion between different individuals to form vegetative heterokaryons that contain genetically distinct nuclei within a common cytoplasm. Heterokaryon formation has potential benefits of functional diploidy and mitotic genetic exchange (“parasexual cycle”; Pontecorvo 1956) in these haploid organisms. However, the formation of viable heterokaryons is believed to be virtually excluded in nature by the action of genetic differences at heterokaryon incompatibility (het) loci [Mylyk 1976; Pandit and Maheshwari 1996; also termed vegetative incompatibility (vic) loci; Leslie 1993]. Hyphal anastomosis between individuals that have alternative allelic specificities at any het locus usually results in compartmentation and death of the hyphal fusion cell, a phenomenon referred to as either vegetative or heterokaryon incompatibility (Glass and Kuldau 1992; Leslie 1993; Saupe 2000).
Vegetative incompatibility has been described in numerous filamentous ascomycetes and basidiomycetes (Glass and Kuldau 1992; Leslie 1993; Worrall 1997). Genetic studies in several ascomycete species show that there are several effectively unlinked het loci: at least 11 in Neurospora crassa (Perkins and Davis 2000), 9 in Podospora anserina (Saupeet al. 2000), 8 in Aspergillus nidulans (Jinkset al. 1966; Anwaret al. 1993), and 7 in Cryphonectria parasitica (Anagnostakiset al. 1986; Cortesi and Milgroom 1998). Generally, only two or three alternative allelic specificities are found at each het locus. Two types of genetic systems have been described that regulate vegetative incompatibility, allelic and nonallelic (Glasset al. 2000; Saupe 2000). In allelic systems, hyphal anastomosis between individuals that contain alternative specificities at a single het locus triggers vegetative incompatibility. In nonallelic systems, an interaction between specific alleles at two different het loci mediates vegetative incompatibility. In N. crassa, C. parasitica, and A. nidulans allelic systems have been described, whereas in P. anserina both allelic and nonallelic systems have been characterized. In this article, we are concerned with allelic systems such as that found in N. crassa.
In N. crassa,11 het loci have been genetically characterized (Perkins and Davis 2000) and 3 loci have been molecularly characterized. Most loci are biallelic, although hetc and het8 encode three alternative allelic specificities (Howlettet al. 1993; Saupe and Glass 1997). N. crassa populations are polymorphic at het loci (Mylyk 1976; Debetset al. 1994). One of the loci characterized molecularly, hetc, encodes three alternative allelic specificities that show transspecies polymorphisms (Saupe and Glass 1997; Wuet al. 1998). The transspecific lineages are defined by insertion/deletions and are thus unlikely to be the result of convergent evolution of independent mutations. Phylogenetic analysis of other markers showed that species with transspecific hetc lineages were monophyletic (Skupskiet al. 1997) and that genealogies of alleles within a functional hetc class were concordant with the inferred pattern of speciation (Powellet al. 2001), indicating that horizontal transfer is an unlikely explanation for transspecies polymorphisms at hetc. This transspecies polymorphism suggests that some form of balancing selection is operating at hetc, a suggestion further supported by the finding that individuals containing alternative specificities were equally frequent in Neurospora populations (Wuet al. 1998). Genetic analysis of an N. crassa population also showed equal frequencies of molecular polymorphisms associated with the two alternative het6 specificities (MirRashedet al. 2000), although similar patterns are not always found at fungal incompatibility loci (Milgroom and Cortesi 1999).
The role of selection in shaping fungal incompatibility systems is debated. The principal question is whether the incompatibility phenotype is itself favored by selection or whether it is an accidental byproduct of selection for some other phenotype (Bégueretet al. 1994). Explicit population genetic models can help clarify the issue by predicting the patterns of polymorphism expected under different hypotheses about selection. Some features of the het system, such as large numbers of vegetative compatibility groups (VCGs), can be recreated with a model of neutral evolution and a high mutation rate (Nauta and Hoekstra 1996). Transspecies polymorphism, extending >30 million years into the past (Wuet al. 1998), however, cannot be explained by a neutral model without postulating extremely large population sizes, which would imply transspecies polymorphism at loci other than het loci as well. Transspecies polymorphism is one of the most striking outcomes of balancing selection in the known recognition systems MHC and SI. In this article, we assume that balancing selection on het loci results from their effect on vegetative incompatibility and find the conditions under which transspecies polymorphism results from this kind of selection.
The exact mode by which balancing selection might be acting at het loci in nature is not clear. The two observed phenotypes are (i) viable heterokaryon formation between genetically similar individuals (identical at all relevant het loci) and (ii) inviability of heterokaryons resulting from genetic difference at a het locus. If there is balancing selection, the failure to form a viable heterokaryon must be favored over the ability to form a viable heterokaryon. Heterokaryon inviability due to vegetative incompatibility has been associated with a reduced risk of transmission of infectious cytoplasmic elements, such as viruslike dsRNAs and senescence plasmids (Debetset al. 1994; van Diepeningenet al. 1997), which are transmitted via hyphal fusion, and in a reduction in exploitation by aggressive genotypes (Debets and Griffiths 1998). Vegetative incompatibility has also been postulated to play a role in limiting outbreeding in certain species (Esser and Blaich 1994). Some of these hypothesized disadvantages of heterokaryon formation have been incorporated into theoretical models of vegetative incompatibility (Hartlet al. 1975; Nauta and Hoekstra 1994, 1996; Liuet al. 2000). Hartl et al. (1975) considered a deterministic twolocus model, with one locus controlling incompatibility and the other (nuclear) locus providing a means for selection. Nauta and Hoekstra (1994) extended this deterministic model to multiple incompatibility loci and cytoplasmically transmitted selective elements and later to a stochastic model with mutation and loss of new VCGs (Nauta and Hoekstra 1996). They derived an expression for the number of VCG types expected at equilibrium in an asexual population and used simulations to study the dynamics of VCG type in populations that had varying degrees of recombination. Transfer of a nongenetic cytoplasmic element between individuals of the same compatibility type was studied by Liu et al. (2000) in a spatially explicit model based on the infection of chestnut blight fungus (C. parasitica) by hypoviruses.
In this article, we propose a general model of heterokaryon formation in which fusions between individuals identical at all het loci result in a loss in fitness to the individuals involved. This model of selection (“compatibility selection”; Nauta and Hoekstra 1996) is an approximation to more detailed models and is roughly equivalent to replacing the disadvantages experienced by particular fusions by the average of such disadvantage in the population. We assume linkage equilibrium among het loci. Restricted recombination, in the form of asexual reproduction, has previously been shown to have a minor effect on the number of VCGs maintained in a population, unless the percentage of the population reproducing asexually is very large (Nauta and Hoekstra 1996). An additional effect of asexual reproduction is an increase in the number of encounters that a particular haploid genotype has with others during its “lifetime.” Liu et al.'s (2000) simulation study found that virus transmission was increased in asexual, highly spatially structured populations of fungi, in which individuals often encountered others of the same compatibility type. We account for the possibility of multiple vegetative encounters between sexual generations in our model.
THE MODEL
Diallelic loci: We illustrate our method in a simple context first and restrict our analysis initially to loci at which only two alleles are possible. We assume that a total of l unlinked loci can affect vegetative incompatibility and that at a given time i of them are polymorphic. Selection is assumed to be sufficiently strong that the two alleles at each polymorphic locus are equally frequent. In our analysis, i is treated as a random variable that increases when a previously monomorphic locus becomes polymorphic and decreases when one of the polymorphic loci becomes monomorphic. We use the strong selection, weak mutation (SSWM) approximation introduced by Gillespie (1984). The SSWM approximation assumes that mutation is sufficiently infrequent that the time intervals between the appearances of successful mutations are very long and that selection is sufficiently strong that a successful mutation reaches its equilibrium frequency instantaneously on the timescale set by the interval between mutations. The SSWM approximation has proved very useful in the analysis of other models of balancing selection (Sasaki 1989, 1992; Slatkin and Muirhead 1999; Slatkin 2000) because it reduces a highly multiallelic model to a simple Markov chain for the numbers of alleles of different type. In the analysis here, the state variable of the Markov chain is i, the number of dimorphic loci.
Assuming linkage equilibrium, the probability that two randomly chosen haploid individuals have the same genotype is ½^{i} and the probability that they differ in genotype is 1 – ½^{i}. If they have the same genotype, the two individuals undergo vegetative fusion that with probability s causes each individual to die or fail to reproduce (for example, because of the transfer of deleterious agents such as mycoviruses). We assume that each individual encounters E other individuals during its lifetime and that each of the E encounters is with an individual chosen randomly and independently from the population. The probability of an individual surviving a single encounter is 1 – ½^{i} + (1 – s)/2^{i} = 1 – s/2^{i} and the probability of it surviving E independent encounters is w = (1 – s/2^{i})^{E}, which is ∼1 – Es/2^{i} if Es/2^{i} is small, as it will be in applications of this theory. The quantity w represents the average fitness of the population at equilibrium.
To employ the SSWM approximation, we first consider the fate of a new mutation at one of the l – i monomorphic loci. An individual carrying such a mutation will not form a stable heterokaryon with any other individual and hence will have a relative fitness of 1. While such a mutation is in very low frequency, the deterministic rate of increase in frequency is ∼1 – w. Using the SSWM approximation, then, 2(1 – w) = Es/2^{i–1} is the probability in a finite population that the rare mutation will increase in frequency, creating an additional polymorphic locus and increasing i by 1. In a population containing N individuals, the number of mutations at previously monomorphic loci per generation is, on average, Nu(l – i), where u is the per locus mutation rate to alleles affecting compatibility type. Therefore, the net probability per generation that a mutant appears that will result in an increase in i by 1 is
Stochastic loss of one allele at a diallelic locus results in a decrease in i by 1. To find the rate of stochastic loss, we assume that an allele A at one of the i diallelic loci has a frequency x but that the frequencies at the other i – 1 polymorphic loci are ½. We can approximate the expected change in x, a(x), (following the notation of Ewens 1979) under selection alone. The marginal fitness of an individual with A is ∼w_{A} = 1 – xsE/2^{i–1} and the marginal fitness of an individual with the other allele, a, is ∼w_{a} = 1 – (1 – x)sE/2^{i–1}. Therefore the deterministic change in x under one generation of selection is approximately
The integral in Equation 3 cannot be expressed in closed form but it can be approximated very well by
Continuing with the SSWM approximation, we assume that loss or fixation of an allele occurs in a single generation, and hence 1/t̄(½) is the probability per generation that a diallelic locus will become monomorphic. Each of the i diallelic loci can become monomorphic so
Equations 1 and 6 define the Markov chain. This chain has a particularly simple form and its stationary distribution can be found explicitly. The probability that i takes a particular value is α_{i} (0 ≤ i ≤ l), which is given by
Numerical analysis of (7) shows that for given values of S and M, there is a critical value of l that distinguishes cases in which the number of loci is limiting from cases in which a balance is achieved between mutation and stochastic loss. The critical value of l, l*, is the value for which λ(l – 1) = μ(l). Figure 1 shows that l* is approximately a linear function of log_{10}(S) and is nearly independent of M. If l < l*, λ(i – 1) >μ(i) for all i < l, which means that the rate of gain of diallelic loci always exceeds the rate of stochastic loss. In that case, the stationary distribution is piled up at l implying that all available loci are dimorphic. If a locus becomes monomorphic, it will remain so for only a relatively short time. The biological interpretation is that the number of loci at which mutations can create new VCGs is smaller than the equilibrium number that would otherwise be maintained.
If l > l*, then the stationary distribution of i is centered roughly at the value of i that solves the zero flux condition, λ(i) = μ(i), which can be simplified to
More than two alleles per locus: It is straightforward to include loci with more than two alleles. Of the l loci that can affect heterokaryon formation, c_{i} carry i alleles at any one time, so that the population consists of c_{1} monomorphic loci, c_{2} dimorphic, and so on, with a maximum of k alleles at any locus. We construct a Markov chain on the set of integers,
As in the model of diallelic loci, each new mutant has a relative fitness of 1 because an individual carrying it cannot form a stable heterokaryon with any other individual in the population. Hence, each new mutant has a selective advantage of 1 – w over alleles already present at their equilibrium frequencies. Thus, the probability that a new mutation at a locus with i – 1 alleles increases in frequency to become common is ∼2EsF. With c_{i}_{–1} such loci present, the probability that this event occurs is
The rate of stochastic loss of alleles at a locus with i > 2 alleles is slightly different from that for a diallelic locus because only loss and not fixation of an allele must be considered. The mathematical problem is the same as that analyzed previously by Sasaki (1989), Takahata (1990), and Slatkin and Muirhead (1999). We consider an allele A at frequency x and assume that the other (i – 1) alleles are in frequency (1 – x)/(i – 1). We also assume that allele frequencies at other polymorphic loci are at their deterministic equilibria. The marginal fitness of A is w_{A} = 1 – sixEF and the marginal fitness of all of the other alleles is w_{a} = 1 – si((1 – x)/(i – 1))EF. Therefore, the average change in x per generation is
If the coefficient of (1/i – y)^{2} in (11) is large, then the ultimate loss of A is essentially certain; the probability of fixation of A is proportional to e^{–}^{SF}. In that case, the time to loss is most easily found by ignoring the possibility of fixation and using Equations 4.39 and 4.40 from Ewens (1979) for the case in which only loss of A is possible. The result is, approximately,
Equations 9 and 13 provide the transition probabilities for the Markov chain. There is no absorbing state so the chain has a stationary distribution, but the chain is not a continuant, so there is no analytic expression for the stationary distribution, as there is for the model of diallelic loci (Equation 7). The zero flux conditions, λ_{i}(c) = μ_{i}(c), have to be solved numerically.
As in the case of the model of diallelic loci, the results depend on l and the two composite parameters M and S. The maximum number of VCGs possible in a population, assuming equal frequencies of alleles at each locus and linkage equilibrium between loci, is
As in the previous section, there is a critical number of loci, l*, that determines the behavior of the model. The value of l* is the largest value of l for which λ_{i}(0,... 0, 1, l – 1) >μ_{i}(0,...0, l). If l < l*, the rate of appearance of new mutations exceeds the rate of loss even when every locus has the maximum possible number of alleles, which means that the number of loci and hence the total number of VCGs are limited by the number of loci at which such mutations can occur. For given values of M and S, l* decreases as k increases. For example, with M = 0.01 and S = 10^{5}, l* = 10 for the diallelic (k = 2) model, l* = 6 for k = 3, and l* = 2 for k = 10. When l < l*, most loci have the maximum number of alleles. But for k > 3, we find this can occur only when the rate of mutation to new specificities is very high (M > 1).
If M < 1, even extremely strong selection does not necessarily result in highly polymorphic loci if l > l* (Figure 4). For most of Figure 4, there are essentially no loci with more than three alleles. Some multiallelic loci can be maintained with the higher mutation rate shown, but even with M = 0.1, S = 10^{5}, most loci (6.1 out of 10) have three or fewer alleles. With greater numbers of available loci, the conditions for multiallelism become even stricter. Increasing the mutation rate tends to increase overall diversity by increasing the number of polymorphic loci, rather than by increasing the number of alleles per locus.
We conclude that the behavior of the model with multiallelic loci is similar to that of the model with only diallelic loci. Two general features emerge.

Loci will be equally polymorphic to the extent possible. Only when all loci are diallelic are a significant number of triallelic loci maintained, and more than three alleles per locus require very strong selection or very small numbers of available loci.

The behavior of the model depends on whether the number of loci (l) is limiting. If l is small, then all available loci will be polymorphic to the maximum possible extent.
Transspecies polymorphism: For transspecies polymorphism to be detected, alleles must remain in the population for very long times, much longer than alleles at neutral or weakly selected loci. In a haploid species, the average time until loss of alleles initially present in a moderate frequency is ∼N generations. Transspecies polymorphism is possible if t̄/N ≫ 1, where t̄ is the expected time to loss of an allele initially present at its equilibrium frequency (½, ⅓,...). For the diallelic model, an analytical expression for the average age of an allele is given by (3).
Average allele age for the multipleallele model may be calculated by construction of another Markov chain, using the rates of allele gain and loss derived above. In this Markov chain, the state is the current number of alleles at a locus containing an allele of interest. The state changes as alleles are gained and lost at the locus until eventually either the locus becomes monomorphic or the allele of interest itself is lost from the population. The average age of an allele, given that it starts in state i, i > 1, is
The quantity t̄/N depends strongly on the scaled mutation parameter M and on whether the number of loci is limiting. If l > l*, alleles persist for only slightly longer times than neutral alleles but if l < l* persistence times can be dramatically longer. For example, in the diallelic case, with M = 0.1 and S = 10^{3} (l* = 8), t/N = 43,825.9 for l = 6, 7.6 for l = 8, 2.7 for l = 10, and 1.4 for l = 15, which are the values for which the stationary distributions are shown in Figure 2. This pattern is typical and leads us to the conclusion that for transspecies polymorphism to be possible when k = 2, l must be <l*, implying that all or nearly all of the available loci are polymorphic at any time. There is still turnover of alleles but only on a very long timescale. This conclusion is also true if k > 2, but in this case, biologically plausible values of l (5–15 for most species studied) are generally >l* (Table 1). If k > 2 and l is on the order of 10, alleles turn over more slowly than neutral alleles but not slowly enough to account for transspecies polymorphism. In most cases, average allele age is only slightly greater than in the neutral case (Table 1) and is far less than is usually found with balancing selection.
If k > 2, transspecies polymorphism is possible only if M is very small. The equilibrium level of polymorphism is largely determined by S, l, and k, but the turnover of alleles is driven by mutation. Lower values of M result in longer persistence times. We conclude, then, that transspecies polymorphism appears to require severe limitation of mutational opportunity, through a very low mutation rate, or through some constraints on the number of possible alleles at a locus, or both.
DISCUSSION
We have investigated a model that contains the essential features of vegetative incompatibility of the type found in N. crassa and many other species of filamentous fungi. The model's parameters are (i) the number of loci at which mutations creating new incompatibility types can occur (l), (ii) the number of alleles possible at each of these loci (k), (iii) the mutation rate scaled by the effective population size (M = Nu), and (iv) the scaled selection intensity (S = NsE, where E is the number of opportunities an individual has for vegetative fusion).
There are some constraints on the values of these parameters. For the model to represent strong stabilizing selection, S must be ≫1. It is also biologically reasonable that M be ≪S. Genetic studies cited in the Introduction indicate that the biologically interesting range for l is 5–15.
Our results show that, as long as selection is strong (i.e., S ≫ 1), genetic polymorphism will be maintained at loci affecting vegetative incompatibility and numerous incompatibility types will be found under a wide range of values of the other parameters. The model predicts that as many loci as possible will be polymorphic and that, if the number of loci is not limited, most of the polymorphic loci will have only two alleles. If the number of loci is limited, then all available loci will be polymorphic and some may carry more than two functionally different alleles per locus. To the extent possible, loci will tend to have equal numbers of alleles.
The symmetry of our results and in particular the prediction of equal numbers of alleles per locus is a consequence of the symmetry of the underlying assumptions. Mutation rates and numbers of possible alleles could easily differ among loci, leading to differences in the numbers of alleles observed.
Transspecies polymorphism is not a necessary consequence of selection of this type, as it is for strong balancing selection affecting a single locus. The reason is that, although polymorphism can be maintained at loci responsible for vegetative incompatibility, considerable turnover of alleles can occur. As a result, no allele persists long enough for transspecies polymorphism to be commonly found.
Our model predicts that transspecies polymorphism will be found if the rate of turnover of alleles is constrained in some way. One possibility is that the rate of mutation to functionally different alleles at each locus is very small. That would result in long persistence times because only when a new allele appears is there pressure for existing alleles to be lost by genetic drift. Another possibility is that there are so few loci and so few alleles per locus that, at equilibrium, nearly all possible alleles are present. In that case, few if any new functionally different alleles could enter the population and cause turnover of existing alleles.
On the basis of available information, low mutation rate is a more plausible explanation for observations of transspecies polymorphism than is exhaustion of available alleles. At hetc in N. crassa and related species, alleles exhibiting transspecies polymorphism differ by an insertion/deletion (indel) motif, a mutational event that is much rarer than nucleotide substitutions or the insertion or deletion of only one or two nucleotides. That observation is consistent with the idea that the mutation rate to functionally different alleles is very low. Under laboratory conditions, it was possible to generate a limited number of new hetc allelic specificities by altering the indel pattern in the specificity domain of hetc (Wu and Glass 2001). Because synthetic functionally different alleles have been created, it does not appear as if all possible alleles are present, at least at hetc.
If transspecies polymorphisms of the type found by Saupe and Glass (1997) and by Wu et al. (1998) are not the result of their role in vegetative incompatibility, an alternative explanation is that those alleles directly experience balancing selection because of some pleiotropic effect. Additional evidence of transspecies polymorphism from other loci involved in vegetative incompatibility will help resolve this issue.
Acknowledgments
This work was funded in part by grants from the National Institutes of Health to N.L.G. (GM60468) and M.S. (GM40282).
Footnotes

Communicating editor: M. K. Uyenoyama
 Received December 7, 2001.
 Accepted March 13, 2002.
 Copyright © 2002 by the Genetics Society of America