Toward genome-wide identification of Bateson-Dobzhansky-Muller incompatibilities in yeast: a simulation study.

The Bateson-Dobzhansky-Muller (BDM) model of reproductive isolation by genetic incompatibility is a widely accepted model of speciation. Because of the exceptionally rich biological information about the budding yeast Saccharomyces cerevisiae, the identification of BDM incompatibilities in yeast would greatly deepen our understanding of the molecular genetic basis of reproductive isolation and speciation. However, despite repeated efforts, BDM incompatibilities between nuclear genes have never been identified between S. cerevisiae and its sister species S. paradoxus. Such negative results have led to the belief that simple nuclear BDM incompatibilities do not exist between the two yeast species. Here, we explore an alternative explanation that such incompatibilities exist but were undetectable due to limited statistical power. We discover that previously employed statistical methods were not ideal and that a redesigned method improves the statistical power. We determine, under various sample sizes, the probabilities of identifying BDM incompatibilities that cause F1 spore inviability with incomplete penetrance, and confirm that the previously used samples were too small to detect such incompatibilities. Our findings call for an expanded experimental search for yeast BDM incompatibilities, which has become possible with the decreasing cost of genome sequencing. The improved methodology developed here is, in principle, applicable to other organisms and can help detect epistasis in general.

Introduction

Speciation, the “mystery of mysteries” in Darwin’s words (Darwin 1859), is one of the most important processes in evolution, responsible for the generation of the tremendous biodiversity on Earth. Important as it is, speciation is not well understood at the genetic level. For example, it is unknown how many genetic changes underlie the formation of a new species in nature, and the relative roles of natural selection and genetic drift in causing these changes are debated (Schluter 2009; Nei and Nozawa 2011). A key step in speciation is the establishment of reproductive isolation, which can occur pre- or postzygotically (Coyne and Orr 2004). Genetic incompatibility is thought to be the major cause of postzygotic isolation. Specifically, the Bateson–Dobzhansky–Muller (BDM) model asserts that a genetic change at locus A in one population and a genetic change at locus B in another population may be incompatible when residing in the same genome upon the hybridization between individuals of the two populations, which could result in postzygotic incompatibility and lead to inviability, infertility, or inferiority (Orr 1996). Although this model is generally accepted, only a small number of genes in a few species pairs have been identified to be genetically incompatible (Wu and Ting 2004; Maheshwari and Barbash 2011; Nosil and Schluter 2011). One classical example involves the melanoma formation in the hybrids of Xiphophorus species. Normally, the Tu locus controls the formation of spots composed of black pigment cells. In interspecific hybrids between the platyfish X. maculatus and swordtail X. helleri, these spots sometimes spontaneously develop into malignant melanomas (Wittbrodt et al. 1989). A two-locus BDM model can explain this phenomenon: overexpression of Tu, which has been identified to be Xmrk on the X chromosome, causes melanomas to form (Adam et al. 1993), whereas an autosomal repressor gene mapped near cdkn2a/b negatively regulates Tu (Schartl et al. 2013). The hybrids that have Tu but not the repressor will develop melanomas (Meierjohann et al. 2004). There is, however, much disagreement on the existence of such major BDM incompatibilities and their role in speciation in general (Liti et al. 2006; Maheshwari and Barbash 2011). Identifying such genes and studying their functions and evolution could help settle this debate and uncover the molecular genetic basis of reproductive isolation and speciation. Because BDM incompatibilities are expected to accumulate with the divergence of two species, identifying such incompatibilities from closely related species is most relevant to understanding speciation (Nosil and Schluter 2011).

For four reasons, the budding yeast Saccharomyces cerevisiae (Sc) and its sister species S. paradoxus (Sp) appear to be ideal for identifying BDM incompatibilities and studying their roles in speciation. First, S. cerevisiae is one of the best studied eukaryotes, with abundant information on its genetics, genomics, physiology, cell biology, and molecular biology. There are also numerous genetic tools and methods for studying yeast. Its short generation time allows rapid genetic analysis and its small genome (∼12 million bases) makes genotyping and fine genetic mapping easier than in other species. Second, separated approximately 10 Ma (Kawahara and Imanishi 2007) and with approximately 85% genome sequence identity (Kellis et al. 2003), Sc and Sp are relatively closely related. The two species can readily mate with each other (Murphy et al. 2006); yet, their postzygotic isolation is strong, with Sc–Sp hybrids producing only approximately 1% viable spores (Hunter et al. 1996). Third, the genomes of the two species are essentially collinear with no gross chromosomal rearrangements and no reciprocal translocation; only four inversions and three segmental duplications exist (Kellis et al. 2003). This fact eliminates chromosomal rearrangement as a major contributor to their postzygotic isolation. Fourth, the genotypes and phenotypes of yeast haploids can be directly analyzed, avoiding the need to generate homozygotes from the spores produced by F1 hybrids. Note that F1 hybrids are not suitable for identifying genetic incompatibilities unless they are dominant, but a previous study has excluded the existence of dominant genetic incompatibilities underlying the infertility of the hybrid between Sc and Sp (Greig et al. 2002). One complication of the yeast system is that a large fraction of spores produced by Sc–Sp hybrids are killed by aneuploidy (Hunter et al. 1996). At least one recombination is usually required for correct segregation of homologous chromosomes during meiosis. In the Sc–Sp hybrid, the sequence differences between homologous chromosomes cause the mismatch repair system to suppress recombination, resulting in a high frequency of aneuploidy (Chambers et al. 1996). Deleting the mismatch repair gene MSH2 increases the recombination rate in the hybrid from 5.4 to 35.6 crossovers per meiosis (Kao et al. 2010). Consequently, F1 spore viability rises to approximately 10% (Kao et al. 2010).

Research in the last decade has focused on understanding the genetic basis of Sc–Sp F1 hybrid infertility, which is equivalent to F1 spore inviability. Despite the multiple advantages of the study system and repeated efforts (Greig et al. 2002; Greig 2007; Kao et al. 2010; Xu and He 2011), no nuclear–nuclear genetic incompatibilities have been identified for Sc–Sp F1 infertility, although a mitochondrial-nuclear incompatibility has been reported for F2 hybrid infertility (Chou et al. 2010). Two general strategies have been used to identify nuclear–nuclear genetic incompatibilities between Sc and Sp. The first approach is to replace chromosomes in Sc with their Sp homologs one at a time. If interchromosomal incompatibilities exist, one would observe a reduction in strain fertility, viability, or growth rate upon a chromosomal replacement. The fact that such replacements were made for at least 9 of the 16 chromosomes demonstrates the lack of BDM incompatibility for F1 spore viability in the 9 chromosomes (Greig 2007). This result, however, does not exclude the possibility of incompatibilities for F1 spore growth rate or higher order incompatibilities for viability. Note that even when an interchromosomal incompatibility is detected using this approach, further work is needed to localize the incompatible genes.

The second approach is to identify genetic incompatibilities in F1 spores by linkage analysis. Briefly, if the Sc allele at locus A (A) is incompatible with the Sp allele at locus B (B), spores of the genotype A may have reduced viability and thus may be underrepresented among viable F1 spores. This decrease in frequency also applies to pairs of markers closely linked to A and B, respectively. Thus, it is possible to use existing genetic markers such as single nucleotide differences (SNDs) between the two species to map BDM incompatibilities. This approach is virtually identical to mapping genetic interaction or epistasis. Because of the large number of marker pairs to be tested, the statistical power is expected to be low.

Two groups have used the above second approach to look for incompatibilities between Sc and Sp that kill F1 spores with 100% penetrance, but with no success (Kao et al. 2010; Xu and He 2011). The negative result has led to the suggestion that two-locus BDM incompatibilities do not exist in yeast and are unimportant to yeast speciation (Kao et al. 2010). However, for two reasons, genetic incompatibility need not have 100% penetrance. First, an incompatibility may only increase the probability of spore inviability rather than killing the spore deterministically, because spore viability is likely to be a complex trait controlled by multiple genes. Second, a high-order incompatibility behaves like a two-locus incompatibility with incomplete penetrance. For instance, a three-locus incompatibility with 100% penetrance behaves exactly as a two-locus incompatibility with 50% penetrance. Given the possibility of incomplete penetrance, one wonders what conclusion about the genetic incompatibility between Sc and Sp can be drawn from the existing data of the linkage analysis. To answer this question, it becomes necessary to understand the properties of this linkage analysis. Here, we use computer simulation to inspect the statistical properties of the linkage analysis, under the scenario that two-locus genetic incompatibilities cause F1 spore inviability with incomplete penetrance, which, as aforementioned, includes the possibility of multiple-locus incompatibility. We show that the previously designed statistical method is not ideal and propose a modified method that improves the statistical power. We find previously used sample sizes too small to detect genetic incompatibilities and offer guidelines for future experimental searches of the BDM incompatibilities between Sc and Sp. These results may apply to the study of BDM incompatibilities in other species and more generally to epistasis mapping.

Materials and Methods

General Strategy of Simulating the Identification of BDM Incompatibilities

Based on theoretical predictions and experimental results (Welch 2004; Wu and Ting 2004; Lee et al. 2008), we assume that genetic incompatibilities are asymmetric. That is, if A and B are incompatible, A and B can still be compatible (fig. 1A). We define I as the probability that an F1 spore dies due to an incompatible allelic pair. We consider the use of msh2 mutants of both Sc and Sp in this study (Kao et al. 2010) such that spore deaths have three potential causes: random death, aneuploidy, and genetic incompatibility. Random death refers to spore death caused by deleterious mutations, meiotic errors, or environmental factors, and is assumed to have the same rate in the parental species and their hybrid.

F

General strategy of simulating the identification of BDM incompatibilities between Saccharomyces cerevisiae (Sc) and S. paradoxus (Sp). (A) The Sc allele at locus A and the Sp allele at locus B are incompatible, leading to reduced viability when in the same spore. (B) Procedure for detecting BDM incompatibility between Sc and Sp. (C) A 2 × 2 table for spore counts of each marker pair. Several statistics for genetic incompatibility are computed using these counts.

The following steps outline the procedure of simulating spore production (fig. 1B). First, to simulate the hybridization between the two yeast species, we set the in silico genome to contain 16 chromosomes with lengths following those of Sc. SND density was set to be one per seven nucleotides based on the 85% sequence identity between the two species. We assume N pairs of incompatibilities and randomly assign them to the existing SNDs. The effects of these N pairs of incompatibilities on F1 spore inviability were either set to be equal or set to follow a certain distribution. The number of crossovers generated during the meiosis of F1 hybrids followed a Poisson distribution with a mean of 35.6 per meiosis (Kao et al. 2010) and the crossovers were randomly assigned to the genome. Meiotic gene conversion and variable recombination rates across the genome are not considered. After meiosis, four spores are generated. We then calculate spore viability as described in the next section and stochastically determine viable spores based on their viabilities.

In actual experiment, the viable spores may be genotyped by restriction enzyme digestion (Xu and He 2011), microarray-based SND typing (Kao et al. 2010), or genome sequencing. Here, we use 1,207 SNDs (1 per 10 kb) as markers in linkage analysis. Using more markers does not improve the precision or power of identifying BDM incompatibilities because of limited recombination in msh2 Sc–Sp hybrids: 10,000 nucleotides correspond to 1.5 cM. Using 1 marker per 10 kb means that the expected mapping resolution is at best 2.5 kb.

Our preliminary analysis revealed that any BDM incompatibility between two intrachromosomal loci is difficult to detect due to strong linkage. Hence, we examine the frequencies of spores for every pair of interchromosomal SND markers. That is, for markers A and B that are located on different chromosomes, we obtain the numbers of spores with the genotypes of A (a), A (b), A, (c), and A (d), respectively. These numbers form a 2 × 2 table (fig. 1C), from which three statistics are calculated: chi-squared value, G test statistic, and odds ratio (OR) (discussed later). Because of viability differences among the four genotypes, the incompatible genotype should have a reduced frequency, compared with its expected value.

In theory, when the sample size is sufficiently large, we should be able to recover the pre-assigned incompatible allelic pairs. After acquiring a statistic of genetic incompatibility for each pair of markers, we determine statistical significance using a familywise 5% type I error rate (discussed later). We then attempt to estimate the chromosomal segments encompassing the incompatibility genes (discussed later).

Calculating Spore Viability

In our simulation, random death, aneuploidy, and BDM incompatibility are three causes of F1 spore inviability. We set the random death rate to be R = 1–0.804 = 0.196, based on the fact that S. cerevisiae and S. paradoxus msh2 mutants have spore viabilities of 84.0% and 80.4%, respectively (Hunter et al. 1996). It has been estimated that aneuploidy occurs at a frequency of 0.29 per viable msh2 Sc–Sp hybrid spore (Kao et al. 2010), but it is unknown what the corresponding fraction is in dead spores. The impact of aneuploidy on spore viability is complicated. Although loss of a chromosome is lethal, gain of an extra chromosome could be beneficial if it masks the deleterious effect of genetic incompatibility. We set the probability of spore inviability due to aneuploidy to be either U = 0% or 50% to obtain a minimal and a more realistic estimate of the required sample size for identifying BDM incompatibilities, respectively. Inviability caused by aneuploidy is applied to pairs of sister spores because nondisjunction typically occurs in meiosis I of the hybrid (Hunter et al. 1996). We assume no epistasis among incompatible gene pairs. Let T be the fraction of viable spores produced by F1 hybrids, N be the number of BDM incompatibility pairs between Sc and Sp, and I be the probability of spore death caused by the kth pair of incompatibility or penetrance. We have In the simple case of I = I for all N incompatible pairs, we have

Statistics Characterizing Genetic Incompatibility

Genetic incompatibility between A and B leads to a reduction in the frequency of A, compared with its expected value. This signal can be detected in multiple ways. Because of strong linkage within a chromosome, we only evaluate pairs of markers that reside on different chromosomes. In a previous study (Kao et al. 2010), a chi-squared test was used to test whether the frequency of a recombinant equals the product of corresponding allele frequencies. For example, if the A and B frequencies among viable F1 spores are 0.3 and 0.5, respectively, the expected frequency of viable A spores is 0.3 × 0.5 = 0.15. Chi-squared is then calculated by summing over all genotypes the squared difference between the expected and observed numbers of a genotype divided by the expected number. This test is nondirectional in the sense that it does not distinguish whether the recombinants are overrepresented or underrepresented. Besides the chi-squared test, the G test of independence may be used to test the goodness of fit of the observed genotype frequencies to their expected values. G test is designed for cases where the margins of a 2 × 2 table are not fixed by investigators whereas the total number in the four cells of the table is fixed (Sokal and Rohlf 1995). We conduct the G test with Williams’s correction (Sokal and Rohlf 1995). In addition, we calculate an OR by dividing the product of the numbers of the two parental genotypes by that of the two recombinant genotypes: OR = (a × d)/(b × c) (fig. 1C).

Because multiple pairs of markers are tested in an experiment, we evaluate the significance of the earlier statistics by controlling the familywise type I error rate. We first randomly shuffle each of the 16 chromosomes among spores and then find the highest statistic among all pairs of markers. We conduct this shuffling 100 times and rank the resulting 100 highest statistics. The 5th largest number among these 100 numbers is chosen as the critical value corresponding to a familywise type I error rate of 5%.

After applying the cutoff, we group statistically significant pairs of markers as follows. Let us use the OR as an example, but the same procedure applies to the other statistics used. First, we find the maximal OR, and take a step of seven markers on each side of each focal marker to obtain the initial square of close linkage. The number seven is chosen by considering the tradeoff between grouping markers showing signals of different incompatibilities and dividing markers showing the signal of the same incompatibility. Second, we keep expanding the square with a step size of one marker until it is no longer significant or it reaches an end of the chromosome. Third, if two squares overlap with each other, we ignore the square with the lower maximal OR. Fourth, we repeat these steps until all significant pairs of markers are included in the squares. Fifth, the marker pair of the maximal OR of each square is recorded. If two adjacent marker pairs in the same square tie for the maximal OR, we record the locations of their midpoints.

A preassigned BDM incompatible pair is considered to be correctly identified when both causal SNDs are within seven markers from the maximum in an aforementioned square. Sensitivity is calculated as the fraction of true incompatible pairs identified. False discovery rate is calculated as the total number of false discoveries divided by the total number of discoveries. When no discovery is made in all simulations, false discovery rate is defined as 0. Genomic distance is calculated as the average distance between the two identified markers and their respective causal SNDs. Standard errors of sensitivity, false discovery rate, and genomic distance estimates are estimated using 1,000 bootstrap samples.

Results

OR Outperforms Other Statistics in Identifying Genetic Incompatibility

Following Kao et al. (2010), we use msh2 mutants of Sc and Sp in our simulation of identifying BDM incompatibilities, unless otherwise noted. Based on theoretical predictions and experimental results (Welch 2004; Wu and Ting 2004; Lee et al. 2008), we assume that genetic incompatibility is asymmetrical. That is, if A and B are incompatible, A and B can still be compatible (fig. 1A). It is difficult to detect BDM incompatibility between two loci that reside in the same chromosome because of limited recombination in the hybrid. Hence, we only examine pairs of markers located on different chromosomes. That is, for markers A and B that are located on different chromosomes, we obtain the numbers of spores with the genotypes of A (a), A (b), A (c), and A (d), respectively, which form a 2 × 2 table (fig. 1C). Because of viability differences among the four genotypes, the incompatible genotype should have a reduced frequency, compared with its expected value (fig. 1B). In theory, when the sample size is sufficiently large, we should be able to detect such incompatible allelic pairs.

We calculate three statistics using the 2 × 2 table: chi-squared, G test statistic, and OR = (a × d)/(b × c) (see Materials and Methods), and evaluate their relative performances in identifying preassigned incompatibilities by simulation. The chi-squared statistic was previously used in this context (Kao et al. 2010), but this statistic does not differentiate between overrepresentation and underrepresentation of a genotype relative to its expectation and thus may be less specific. Because chi-squared test is an approximation of the G test, they have similar properties, although G test may be more precise. By contrast, a lower-than-expected OR indicates overrepresentation of A and/or A, whereas a higher-than-expected OR indicates depletion of these genotypes, which is predicted under genetic incompatibility. After acquiring a statistic of genetic incompatibility for each interchromosomal marker pair, we determine statistical significance using a familywise 5% type I error rate to control multiple testing. We then identify the chromosomal segments that are likely to encompass the incompatibility genes (see Materials and Methods).

Because the incompatible marker pairs are preassigned in the simulation, we can evaluate how well the three statistics perform in terms of the 1) sensitivity, 2) false discovery rate, and 3) mean genomic distance between the identified markers and the preassigned incompatible SNDs. For each parameter set, we conduct 400 simulation replications and pool the data in our analysis. Sensitivity is the fraction of all preassigned incompatible pairs that are recovered by the analysis. False discovery rate is the number of false discoveries divided by the total number of discoveries. The standard errors of these estimates are estimated by bootstrapping the pooled data 1,000 times. There are 12.07 million nucleotides × 15% = 1.8105 million SNDs between Sp and Sc. We randomly assigned N pairs of SNDs to form N incompatibility pairs. In mapping these incompatibilities, however, we use only 1,207 markers, or 1 marker per 10,000 nucleotides, because the use of more markers does not increase mapping resolution due to limited recombination (see Materials and Methods).

We start the simulation with the following parameters. We assume no contribution of aneuploidy to spore inviability, and set N = 10 pairs of incompatibilities that have equal effects on inviability. Given the known viability of msh2 hybrid spores, the 10 pairs each contribute I = 0.75 to spore inviability. That is, a spore with one pair of incompatibility is 25% as viable as a spore without any incompatibility. The 10 pairs of incompatibilities (i.e., 20 causal SNDs) are randomly distributed in the 16 yeast chromosomes. The number of viable spores genotyped is M = 200. When OR is used, the sensitivity is 40%, significantly greater than that of chi-squared (28%) or G test statistic (30%) (fig. 2A). The false discovery rate under OR is 24%, not significantly different from that under the other two statistics (22% and 23%, respectively) (fig. 2B). The mean genomic distance between the identified marker and the preassigned incompatibility loci is 18.3 kb under OR, significantly smaller than that under the other two statistics (19.3 and 19.1 kb, respectively) (fig. 2C).

F

Performances of OR, chi-squared, and G test statistic for detecting BDM incompatibilities. Data shown are from 400 simulations of 10 incompatible pairs with equal I and no contribution of aneuploidy to spore inviability. The sample size is 200 viable spores. OR, χ2, and G represent odds ratio, chi-squared, and G test statistic, respectively. χ2* and G* respectively consider χ2 and G only when OR > 1. Standard error, shown by error bars, is estimated by 1,000 bootstrap replications. (A) Sensitivity of the five tests. P values are from paired t test. (B) False discovery rates of the five tests. (C) Average genomic distance between preassigned incompatibilities and the identified significant markers.

If the differences among the three methods are simply due to the fact that chi-squared and G test statistic cannot distinguish whether parental or nonparental types are in excess, we could use the directional information from OR and consider only those chi-squared or G test statistic values when OR > 1. Although such modified chi-squared and G test statistic outperform their original versions in sensitivity, they are still worse than OR (fig. 2A). In terms of the false discovery rate, the modified versions appear worse than the original versions (fig. 2B). In terms of the genomic distance, the modified versions are similar to the original versions (fig. 2C). We subsequently confirmed the advantage of OR over chi-squared and G test statistic in multiple conditions, by varying N, M, and the influence of aneuploidy (U) (table 1). When genetic incompatibility is symmetrical, however, the advantage of OR over chi-squared and G test statistic disappears (supplementary table S1, Supplementary Material online).

Table 1

OR Outperforms Other Statistics in Detecting Asymmetrical Genetic Incompatibilities

Parameters				Sensitivity (%)					False Discovery Rate (%)					Genomic Distance (kb)
Ua	Nb	Ic	Md	ORe	χ2f	Gg	χ2*h	G*i	OR	χ2	G	χ2*	G*	OR	χ2	G	χ2*	G*
0	8	0.92	100	45.4	22.8**	30.6**	27.8**	35.7**	19.4	26.3**	23.5*	30.4**	27.1**	18.7	20.2*	20.1**	20.4**	20.3**
0	8	0.92	200	85.8	77.3**	81.2**	79.4**	82.2**	15.4	15.2	14.7	17.8**	17.2**	14.3	15.9**	15.3**	15.8**	15.3**
0	8	0.92	400	88.9	88.8	88.9	88.8	88.8	17.7	12.9**	13.3**	15.7**	16.0*	11.4	12.4**	12.0**	12.4**	12.0**
0	8	0.92	800	88.4	88.7*	88.6	88.5	88.5	23.1	18.5**	18.5**	20.4**	20.6**	10.3	10.5**	10.4*	10.5**	10.4*
0	10	0.75	100	8.0	4.1**	4.9**	5.7**	6.7**	37.1	45.8**	43.4**	47.6**	45.7**	25.0	24.1	24.2	24.0	24.1
0	10	0.75	200	40.3	28.2**	30.2**	33.4**	35.6**	24.4	22.0	23.2	26.8*	25.9	18.3	19.3**	19.1**	19.3**	19.1**
0	10	0.75	400	81.8	75.8**	76.9**	77.5**	78.3**	17.2	17.0	16.5	19.0**	18.6*	15.3	16.4**	16.2**	16.4**	16.2**
0	10	0.75	800	88.1	87.9	88.0	87.9	88.1	15.2	12.9**	12.9**	15.3	15.2	12.0	12.7**	12.6**	12.7**	12.6**
0	15	0.52	100	0.7	0.2**	0.3**	0.4*	0.4*	62.3	82.4	78.9	78.5	78.5	30.2	43.4	38.7	35.9	33.9
0	15	0.52	200	3.8	2.4**	2.6**	3.2**	3.2**	48.4	49.3	47.7	54.5	55.1	25.2	26.2	26.3	26.5	26.1
0	15	0.52	400	17.1	11.8**	11.9**	14.5**	14.8**	35.7	35.4	35.4	37.9**	37.5**	21.8	21.8	21.7	22.1	22.1
0	15	0.52	800	56.8	49.3**	50.0**	53.0**	53.4**	23.1	21.2**	21.1**	24.0**	23.8*	18.3	18.8**	18.7**	18.7**	18.7**
0.5	5	0.97	100	77.5	36.0**	48.9**	42.2**	53.9**	18.0	25.3**	24.0**	27.9**	26.3**	18.1	19.4**	20.2**	19.7**	20.6**
0.5	5	0.97	200	90.8	87.9**	90.0**	88.2**	90.0**	16.6	14.0*	13.4**	16.6	17.2	13.2	14.7**	14.3**	14.7**	14.3**
0.5	5	0.97	400	89.6	89.8	89.8	89.7	89.8	17.4	13.1**	13.5**	15.9*	16.3	10.7	11.5**	11.1**	11.5**	11.1**
0.5	5	0.97	800	90.2	90.3	90.3	90.3	90.2	20.8	17.0**	17.0**	19.7	19.8	9.8	9.9	9.8	9.9	9.8
0.5	7	0.72	100	7.2	3.1**	3.6**	4.2**	4.6**	40.8	48.2*	46.0*	53.2**	53.1**	23.9	27.8*	27.2*	25.8	25.5
0.5	7	0.72	200	33.5	18.8**	20.8**	23.2**	25.0**	26.7	31.8**	30.6**	33.9**	32.8**	19.6	20.3	20.3	20.3*	19.9
0.5	7	0.72	400	77.7	68.6**	70.1**	71.8**	72.6**	18.7	17.9	17.9	20.8**	20.7**	16.6	17.4**	17.1**	17.4**	17.2**
0.5	7	0.72	800	89.2	88.6**	88.8*	88.5**	88.7**	14.4	13.2	13.2	15.5*	15.3*	12.8	13.6**	13.4**	13.6**	13.4**
0.5	10	0.52	100	0.9	0.4**	0.5**	0.7	0.7	66.7	75.0	76.2*	79.7	79.1	27.2	32.9	31.8	30.0	29.3
0.5	10	0.52	200	3.4	1.8**	1.8**	2.5**	2.7**	53.8	53.9	55.7	62.0*	60.5	23.3	25.3*	25.4*	24.3	24.3
0.5	10	0.52	400	17.2	11.2**	11.6**	14.2**	14.6**	35.6	37.7	37.1	39.9**	39.4*	20.6	21.1	21.0	20.9	21.0
0.5	10	0.52	800	57.4	49.4**	50.0**	53.2**	53.6**	24.2	23.5	23.4*	26.5**	26.4**	18.1	18.5*	18.4	18.6**	18.5**

Note.—The results are from 400 simulations for each parameter set.

aProbability of aneuploidy-induced inviability.

bNumber of pre-assigned BDM incompatibility pairs.

cProbability of spore death caused by each pair of incompatibility.

dTotal number of genotyped spores.

eOdds ratio.

fχ2 statistic.

gG test statistic.

hχ2 statistic only when OR > 1.

iG test statistic only when OR > 1.

*P < 0.05 when comparing the performance of a statistic with that of OR by a paired t test.

**P < 0.005 when comparing the performance of a statistic with that of OR by a paired t test.

Previous Studies Were Underpowered

To understand why previous experimental searches of nuclear BDM incompatibilities between Sc and Sp were unsuccessful, we perform a simulation following the scheme of a previous experiment study, which genotyped 58 spores from F1 with MSH2 and 48 spores from F1 lacking MSH2 (Kao et al. 2010). Before we started the simulation, we confirmed that no pair of markers in that study (Kao et al. 2010) showed significant OR using our methodology. The simulation parameters used for msh2 spores are the same as described earlier. For mismatch repair proficient spores, random death rate is set to be R = 0.05 (Greig et al. 2002). Given the observed viability of 1% among these spores, the contribution of aneuploidy to spore inviability (U) is calculated using equation (2) to be 91.54% and 95.77%, for the corresponding numbers of 0% and 50% in msh2 spores, respectively. To be consistent with the previous study (Kao et al. 2010), we used the density of 1 marker per 2 kb. Using 1 marker per 10 kb yielded similar results.

Assuming different pairs of incompatibilities in the simulation, we calculate the corresponding probabilities of nondiscovery, which is the probability that no marker pair has an OR that deviates significantly from the expectation at the familywise 5% level. We first assume equal effects on spore viability from all pairs of incompatibilities. When aneuploidy does not reduce msh2 spore viability, at least 8 pairs of incompatibilities are required to explain the observed spore inviability. We found the probability of nondiscovery to exceed 5% in all cases except when N = 8 (fig. 3A). If aneuploidy reduces msh2 spore viability by 50% and correspondingly reduces the viability of MSH2 spores, there should be at least 5 pairs of incompatibilities. Under this assumption, we found the probability of nondiscovery to exceed 0.05 in all cases except when N = 5 (fig. 3B). Thus, it is possible for the previous experiment to have missed all incompatibilities. Our analysis tends to overestimate the power of the previous study, because segments in spores with aneuploidy were ignored in the experimental study (Kao et al. 2010) such that the actual sample size is smaller than the number of sampled spores. Furthermore, we have not considered genotyping errors, which would further decrease the statistical power. It might seem counter-intuitive that the more pairs of genetic incompatibility there are, the more difficult it is to identify any of them. The underlying reason is that the total contribution of all incompatibility pairs on inviability is fixed in this simulation and that all pairs are assumed to contribute equally. Thus, having a larger number of incompatible pairs means a smaller contribution from each pair.

F

Sample size in Kao et al. (2010) is too small to detect BDM incompatibilities with incomplete penetrance. Data shown are from 200 simulations for each parameter set used. (A) Probability of nondiscovery in a study by Kao et al. (2010) when aneuploidy is assumed to cause no msh2 spore inviability (U = 0). White bars show the results for incompatibilities with equal effects (i.e., equal-penetrance), whereas the gray bar shows the result for 150 incompatibility pairs with unequal effects as described in (C). (B) Probability of nondiscovery in the study by Kao et al. when aneuploidy is assumed to cause U = 50% inviability to msh2 spores. White bars show the results for incompatibilities with equal effects, whereas the gray bar shows the result for 100 incompatibility pairs with unequal effects as described in (D). (C) Distribution of the effect sizes (i.e., penetrances) of 150 BDM incompatibility pairs (under U = 0) used for the simulation of the gray bar of (A). (D) Distribution of the effect sizes of 100 BDM incompatibility pairs (under U = 50%) used for the simulation of the gray bar of (B). Error bars in (A) and (B) are standard errors estimated from 1,000 bootstrap samples.

Because multiple pairs of genetic incompatibility are unlikely to have equal effect sizes on spore viability, it would be more realistic to consider unequal effect sizes. The difficulty, however, is that there is no prior knowledge on the effect size distribution. Because BDM incompatibilities may be similar to loss-of-function mutations (Maheshwari and Barbash 2011), we assume that the effect size distribution follows the distribution of the deleterious fitness effects of single-nonessential-gene deletions in yeast (Qian et al. 2012). We randomly sample I from this distribution until the total incompatibility explains the observed spore inviability. The mode of the number of incompatible pairs required to explain the observed spore inviability is 150 (fig. 3C) and 100 (fig. 3D) when the contribution of aneuploidy to msh2 spore inviability is 0% and 50%, respectively. The corresponding distributions of I under the two scenarios used in this simulation study are presented in figure 3C and D, respectively, and the probability of nondiscovery is 79% (fig. 3A) and 77% (fig. 3B), respectively.

Because the study by Kao et al. (2010) was the largest experiment for identifying BDM incompatibilities between Sc and Sp, our results suggest that none of the previous studies on the subject were sufficiently powerful to detect BDM incompatibilities between the two yeasts.

Sample Sizes Required for Identifying BDM Incompatibilities

How many viable spores should be genotyped to identify BDM incompatibilities with a reasonable success rate? Here, we again assume the exclusive use of msh2 strains in the experiment. Under the assumption of no effect from aneuploidy on viability, we examine the sceneries of N = 8, 10, and 15 incompatible pairs with equal effects, respectively. We use the sample size of M = 100, 200, 400, and 800 spores, respectively. In the case of N = 8, the probability of nondiscovery is negligible even when M = 100 (fig. 4A). In the case of N = 10 and 15, the probability of nondiscovery declines quickly as M increases from 100 to 200 and 400 (fig. 4A). As expected, the total number of discoveries increases with the sample size M (fig. 4B), so does the sensitivity (fig. 4C). By contrast, the false discovery rate (fig. 4D) and the mean genomic distance between the causal SNDs and the identified markers (fig. 4E) generally decline with M. We also examined the situation when the probability of msh2 spore inviability due to aneuploidy is 50% and obtained overall similar results (fig. 4F–J). Figure 5 shows randomly picked examples of our simulation results under various M when N is fixed at 10 and U at 0. Because one incompatibility pair happens to reside on the same chromosome, the maximal number of pairs detectable is 9. It is clear how increasing the sample size increases the power of detection. Similar patterns can be seen when U = 0.5 (supplementary fig. S1, Supplementary Material online).

F

Genotyping more F1 spores improves the efficiency of identifying BDM incompatibilities with equal effects. (A) Probability of nondiscovery, (B) number of total discoveries, (C) sensitivity, (D) false discovery rate, and (E) mean genomic distance between the preassigned and identified incompatibilities, when aneuploidy is assumed to have no impact on spore inviability. (F) Probability of nondiscovery, (G) number of total discoveries, (H) sensitivity, (I) false discovery rate, and (J) mean genomic distance between the preassigned and identified incompatibilities, when aneuploidy is assumed to cause a 50% probability of spore inviability. Data shown are from 200 simulations per parameter set. Error bars show standard errors estimated from 1,000 bootstrap samples.

F

An example showing the benefit of using large samples in identifying genetic incompatibilities. (A) Genomic positions of 10 pairs of randomly placed equal-effect genetic incompatibilities in the simulation. Genomic positions are defined by marker numbers on both axes. Note that one pair of incompatibility near marker 1,200 on both axes are located in the same chromosome and therefore are undetectable in our study because only interchromosomal marker pairs are examined. Color shows the expected OR. Spore viability is assumed to be immune to aneuploidy. (B, D, F, H) ORs for all interchromosomal marker pairs when the sample size (number of viable msh2 spores genotyped) is (B) 100, (D) 200, (F) 400, and (H) 800, respectively. Color shows the observed OR (OR < 1 is not shown). (C, E, G, I) Interchromosomal marker pairs whose OR values are significant, when the sample size is (C) 100, (E) 200, (G) 400, and (I) 800, respectively. The identified incompatibilities are circled, with the correct identifications in green and incorrect identifications in red. Note that an incompatible pair is considered to be correctly identified only when both loci of a preassigned pair are within 7 markers (i.e., 70 kb) from an identified OR peak. X and Y labels in (B–I) are the same as in (A).

To obtain a more realistic estimate of the required sample size for detecting incompatibilities, we use the aforementioned unequal effect sizes depicted in figure 3C and D, respectively. Because, under this model, most incompatibilities have small effects, which are hard to detect, we focus on incompatibilities with I > 0.2 and its subset that has I > 0.4, respectively, when evaluating sensitivity, false discovery rate, and genomic distance. The probability of nondiscovery, however, is evaluated as originally defined. As aforementioned, when there is no contribution of aneuploidy to msh2 spore inviability, 150 incompatibility pairs are required to explain the observed spore inviability. Among them, 10 pairs have I > 0.2, four of which have I > 0.4 (fig. 3C). When there is a 50% contribution of aneuploidy to msh2 spore inviability, 100 incompatibility pairs are required to explain the observed spore inviability. Among them, six pairs have I > 0.2, two of which have I > 0.4 (fig. 3D). Our simulation (fig. 6) shows that a much larger sample is required for successful detection of BDM incompatibilities under unequal effect sizes than under equal effect sizes. For example, when M = 1,600, the probability of nondiscovery becomes negligible (fig. 6A and E). With such a large sample, the sensitivity is approximately 40% for I > 0.2 and approximately 80% for I > 0.4 (fig. 6B and F) and the false discovery rate is approximately 30% for I > 0.2 and approximately 50% for I > 0.4 (fig. 6C and G). The mean genomic distance is between 15 and 20 kb for both I > 0.2 and 0.4, respectively (fig. 6D and H).

F

Genotyping more F1 spores improves the efficiency of identifying BDM incompatibilities with unequal effect sizes. (A) Probability of nondiscovery, (B) sensitivity, (C) false discovery rate, and (D) mean genomic distance between the preassigned and identified incompatibilities, when aneuploidy is assumed to have no impact on spore inviability. The effect sizes of the 150 incompatibility pairs are shown in figure 3C. We only show results for incompatibilities with I > 0.2 and I > 0.4, respectively. Probability of nondiscovery refers to the probability of no significant marker pair regardless of effect size. (E) Probability of nondiscovery, (F) sensitivity, (G) false discovery rate, and (H) mean genomic distance between the preassigned and identified incompatibilities, when aneuploidy is assumed to cause a 50% probability of spore inviability. The effect sizes of the 100 incompatibility pairs are shown in figure 3D. Data shown are from 200 simulations per parameter set. Error bars show standard errors estimated from 1,000 bootstrap samples.

Discussion

In this study, we demonstrate that OR outperforms chi-squared and G test statistic in detecting asymmetrical BDM incompatibility through linkage analysis. Our simulation suggests that the existence of two-locus BDM incompatibility between Sc and Sp cannot be excluded and its nondiscovery in previous yeast experiments could be due to the limited sample size and low statistical power. Our study provides important guidelines for designing experiments for identifying yeast BDM incompatibilities and for interpreting potential experimental outcomes. More generally, it highlights the importance of understanding the statistical properties of an experimental method (e.g., sensitivity and false discovery rate) to use it efficiently and interpret the result correctly.

We made several assumptions in our simulation that are worth discussion. First, for simplicity, we assumed that recombination rates are equal throughout the genome and ignored recombination hot/cold spots and interferences between crossovers (Mancera et al. 2008). This assumption should not affect the overall results because of the relatively low marker density used (1 per 10 kb). But recombination rate variation would make the genomic distances between the causal SNDs and the identified markers more variable across the genome. Second, due to the lack of prior knowledge on the distribution of I, we assumed either equal I values for different incompatibility pairs or unequal I values that follow a specific distribution mimicking the fitness effects of gene deletions. We believe that the result from the unequal I are closer to the truth than that from the equal I. Third, we assumed that BDM incompatibility is asymmetrical, which is in accordance with the theory and most of the incompatible pairs identified so far (Wu and Beckenbach 1983; Meierjohann et al. 2004; Welch 2004). Nevertheless, our test still works even when it is symmetrical (supplementary table S1, Supplementary Material online). Fourth, it is unclear how much aneuploidy affects viability in msh2 spores, and we used 0% and 50%, respectively, in our study to have a sense of the range of possible outcomes. Fifth, we assumed no error in genotyping the spores. Although genotyping errors would reduce the statistical power, we expect the genotyping error rate to be low, especially when high-coverage next-generation DNA sequencing is used. Moreover, due to low recombination rates, nearby SNDs can be used for correction of sequencing errors at specific positions. Sixth, we did not explicitly study high-order incompatibility, but because high-order incompatibility is equivalent to two-locus incompatibility with incomplete penetrance, our results apply to high-order incompatibility. For example, I = 0.5 in a two-locus incompatibility (fig. 3) is equivalent to a three-locus incompatibility with 100% penetrance.

In our simulation, we used 1 marker per 10 kb to look for BDM incompatibility. Although next-generation sequencing-based genotyping will offer much more markers, the extra markers do not enhance the mapping resolution, because the low recombination rate in msh2 F1 makes all markers within a 10 kb segment almost completely linked. Because of this property, pairs of incompatible genes that are located in the same chromosome are difficult to detect and therefore are not examined in our simulation. Intrachromosomal incompatible gene pairs are expected to constitute only 7.54% of all incompatible pairs if incompatibility genes are uniformly distributed in the genome.

We found that, by the current method, much larger samples than previously used are required for identifying yeast BDM incompatibilities with incomplete penetrance. Given the rapid increase in DNA sequencing capacity and decline in sequencing cost, genotyping approximately 1,000 spores is no longer out of reach. In fact, a recent study sequenced the genomes of 1,000 F2 individuals from a genetic cross between two yeast strains in order to map quantitative traits (Bloom et al. 2013). Our simulation shows that by genotyping 800 to 1,600 F1 spores, there is a reasonable chance of identifying genetic incompatibilities with relatively high penetrance (>20%).

Given the power of today’s DNA sequencing capacity, an alternative strategy of identifying BDM incompatibility may be used. This strategy involves two steps. First, because an incompatibility allele (e.g., A in fig. 1A) has a fitness of 1–0.25I, relative to its alternative (e.g., A), it is relatively easy to identify it by sequencing a pool of viable F1 spores en masse. Second, after identifying low-fitness alleles, one can then look for their incompatible partners by sequencing individual spores. Because of the reduced number of marker pairs to be tested, the sample size required in the second step will be much smaller. A critical requirement in this design is to minimize the competition among spores in mitotic growth before sequencing them en masse, because allelic differences in growth rate between Sc and Sp that are unrelated to the incompatibility for spore viability may be common.

Although Sc and Sp are used here to parameterize our simulation study, our methodology and results are useful for mapping recessive genetic incompatibilities in other species when the haploid stage can be assayed, including species with haplontic or haploid–diploid life cycles and diplontic species that can undergo homozygous diploidization. Because BDM incompatibility is a type of epistasis, our methods and results also apply in genomic detection of epistasis.

Supplementary Material

Supplementary figure S1 and table S1 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).