Coadapted genomes and selection on hybrids: Fisher's geometric model explains a variety of empirical patterns.

Natural selection plays a variety of roles in hybridization, speciation, and admixture. Most research has focused on two extreme cases: crosses between closely related inbred lines, where hybrids are fitter than their parents, or crosses between effectively isolated species, where hybrids suffer severe breakdown. But many natural populations must fall into intermediate regimes, with multiple types of gene interaction, and these are more difficult to study. Here, we develop a simple fitness landscape model, and show that it naturally interpolates between previous modeling approaches, which were designed for the extreme cases, and invoke either mildly deleterious recessives, or discrete hybrid incompatibilities. Our model yields several new predictions, which we test with genomic data from Mytilus mussels, and published data from plants (Zea, Populus, and Senecio) and animals (Mus, Teleogryllus, and Drosophila). The predictions are generally supported, and the model explains a number of surprising empirical patterns. Our approach enables novel and complementary uses of genome-wide datasets, which do not depend on identifying outlier loci, or "speciation genes" with anomalous effects. Given its simplicity and flexibility, and its predictive successes with a wide range of data, the approach should be readily extendable to other outstanding questions in the study of hybridization.

When individuals of different populations mate, the offspring will carry new combinations of alleles. Sometimes the new combinations bring fitness benefits (heterosis). This is often true, for example, when the parental lines are closely related and highly inbred: a fact that can be exploited in artificial breeding programs. Sometimes, the hybrids are much less fit than their parents (hybrid breakdown), suggesting that the populations may be distinct species. These different outcomes depend on the ways in which the alleles interact, and so comparing the outcomes of different types of hybridization can tell us a lot about gene interactions. We developed a general mathematical multigenic model that makes simple predictions for the fitness of hybrids of any type. We show that our model can account for a large number of empirical patterns, including some that were not well explained by alternative theories, developed for specific cases. We tested our predictions with new data from mussels, and published data from plants and animals, and obtained a good fit. Our framework suggests a new and complementary approach to analyzing genomic data from hybrids, which does not rely on identifying small regions of the genome with anomalous effects.

Hybridization and admixture lead to natural selection acting on alleles in different genetic backgrounds. Most classical studies of hybridization can be placed into one of two classes. The first, involves crosses between closely related inbred lines, where there is no coadaptation between the deleterious alleles that differentiate the lines, such that most hybrids are fitter than their parents. Wright's single‐locus theory of inbreeding was developed to interpret data of this kind (Wright 1922, 1977; Crow 1952; Hallauer et al. 2010). The second, involves crosses between effectively isolated species, where viable and fertile hybrids are very rare. Data of this kind are often analyzed by focusing on a small number of “speciation genes,” and interpreted using models of genetic incompatibilities (Dobzhansky 1937; Coyne and Orr 1989; Orr 1995; Gavrilets 2004; Welch 2004; Kalirad and Azevedo 2017).

The differences between these types of hybridization are clear, but it is equally clear that they are extremes of a continuum. Furthermore, the intermediate stages of this continuum are of particular interest, because they include phenomena such as incipient speciation, and occasional introgression between partially isolated populations (Waser 1993; Rosas et al. 2010; Mendez et al. 2012; Fraïsse et al. 2016a; Duranton et al. 2017). However, it can be difficult to model natural selection in this intermediate regime, not least because models require a large number of parameters when they include epistatic effects between many loci. The empirical study of hybrid genotypes in this regime is also difficult. The analysis of lab crosses often focuses on segregation distortions of large effect, and pairwise incompatibilities (Coyne and Orr 2004; Abbott et al. 2013). This QTL‐mapping framework can miss small effect mutations (Noor et al. 2001; Rockman 2012), which are difficult to identify individually, but whose cumulative effect can be substantial (Boyle et al. 2017).

One promising approach is to use Fisher's geometric model, which assigns fitness values to genotypes using a model of optimizing selection on quantitative traits (Fisher 1930; Orr 1998; Welch and Waxman 2003; Martin and Lenormand 2006). The tools of quantitative genetics have often been used to study hybridization (e.g., Melchinger 1987; Lynch 1991; Demuth and Wade 2005; Fitzpatrick 2008), but Fisher's model is fully additive at the level of phenotype, and the “traits” need not correspond in any simple way to standard quantitative traits (Rosas et al. 2010; Martin 2014; Schiffman and Ralph 2017). Instead, the goal is to generate a rugged fitness landscape, which includes a wide variety of mutational effect sizes and epistatic interactions on fitness, with a minimum of free parameters (Orr 1998; Blanquart et al. 2014; Barton 2017; Hwang et al. 2017).

Here, we build on previous studies (Mani and Clarke 1990; Barton 2001; Rosas et al. 2010; Chevin et al. 2014; Fraïsse et al. 2016b; Schiffman and Ralph 2017), and use Fisher's geometric model to study hybridization. We develop a simple Brownian bridge approximation, and show that it can naturally interpolate between previous modeling approaches (Wright 1922, 1977; Orr 1995; Gavrilets 2004), which are appropriate for the two extreme types of hybridization. We then show how the model can account for surprising empirical patterns that have been observed in both regimes (Wright 1977; Moehring 2011; Moran et al. 2017). Finally, we show that the model yields novel predictions for selection on heterozygosity, and test these predictions with a wide range of new and existing datasets (Table 1).

Table 1

Datasets

Hybridization		N	Sex	#Markers	Cross	g_X	Fitness measure	Reference
Zea mays	inbred lines	–	⚥	–	Various	–	Excess yield	See Fig. 3
Mytilus	edulis/galloprovincialis	132	♀/♂	43	F2	–	–	This study
		144	♀/♂	43	BC1	–	–
Drosophila	pseudoobscura/persimilis	1141/1036	♀	2 X; 11 A	BC1	0.37	Motile sperm: present/absent	Noor et al. (2001)
Drosophila	sechellia/simulans	200/200	♀	8 X; 31 A	BC1	0.17	Sperm quantity: 3‐pt. scale	Macdonald and Goldstein (1999)
Drosophila	santomea/yakuba	550/549	♀	10 X; 22 A	BC1	0.17	Motile sperm: 9‐pt. scale	Moehring et al. (2006a, 2006b)
Teleogryllus	oceanicus/commodus	248	♂	–	Various	0.30	Egg and offspring number	Moran et al. (2017)
Populus	alba/tremula	137	♀♂	∼12,000	WH	–	Survival after 4 years	Christe et al. (2016)
Senecio	aethnensis/chrysanthemifolius	64	⚥	966	F2	–	Necrotic/Healthy	Chapman et al. (2016)
Mus musculus	musculus/domesticus	185	♀	14,220	WH	0.039	Testes weight	Turner and Harr (2014)
		305	♀	202 (16 X; 182 A)	F2	0.039	Prop. abnormal sperm	White et al. (2011)

N, The number of individual hybrids, divided by backcross direction where appropriate; #Markers, The number of genetic markers used to estimate p 12 and h, sometimes divided into X‐linked and Autosomal. BC1, First backcross; WH, Wild hybrids. , weight given to X‐linked markers in the calculation of overall genome composition, calculated from the length of annotated coding sequence.

Models and Results

THE MODELS

Notation and basics

We will consider hybrids between two diploid populations, labeled P1 and P2, each of which is genetically uniform, but which differ from each other by d substitutions. Considering all possible combinations of the two homozygotes and the heterozygote, the populations could generate 3 distinct hybrid genotypes, and each might differ in their overall fitness, or in some component of fitness, such as fertility or viability. Below, we will focus on rank order differences between different classes of hybrid (e.g., high vs. low heterozygosity, males vs. females, F1 vs. F2 etc.). As such, following Turelli and Orr (2000), we describe hybrids using a “breakdown score,” S, which is larger for hybrids that have lower values of the fitness component of interest. Breakdown might relate to fitness via, in which case, the parameter α adjusts the rate at which fitness declines with breakdown. This is related to the overall levels of fitness dominance and epistasis (Hinze and Lamkey 2003; Tenaillon et al. 2007; Fraïsse et al. 2016b), and so these can vary independently of other results. We now define the key quantity f, as the expected value of S for a given class of hybrid, scaled by the expected value for an unfit reference class. Here, , is the expected breakdown score for the class of hybrid with the lowest fitness, on the condition that the parental genotypes are optimally fit. In this case, f can vary between zero, for the best possible class of hybrid, and one, for the worst possible class. We will also consider the case where the parental types are themselves suboptimal, with their own level of “breakdown,” denoted f P1 and f P2. In this case, when , then the value of f for the hybrids can be greater than one (see below).

To define classes of hybrid, we also follow Turelli and Orr (2000). We pay particular attention to interpopulation heterozygosity, and define p 12 as the proportion of the divergent sites that carry one allele from each of the parental types. We also define p 1 and p 2 as the proportion of divergent sites that carry only alleles originating from P1 or P2, respectively. Since , it is convenient to introduce the hybrid index, h, which we define as the total proportion of divergent sites originating from P2 (Barton and Gale 1993). This quantity is closely related to measures of ancestry (e.g., Falush et al. 2003; Fitzpatrick 2012; Christe et al. 2016), although it considers only divergent sites. We can now describe each individual genotype via its heterozygosity, p 12, and hybrid index, h. This means that all hybrids can be represented as points in a triangular space, as shown in Figure 1 A (Gompert and Buerkle 2010; Fitzpatrick 2012). Our goal in this article is to represent this space as a kind of fitness landscape, and show how f can vary with p 12 and h.

Figure 1

Diploid hybrid genomes can be represented as points in a triangular space, indicating their hybrid index, h (the total proportion of divergent alleles that originate from parental line P2), and their interpopulation heterozygosity, p 12 (the proportion of divergent sites that carry one allele from each parental line). Panel (A) shows how the standard crosses are placed in this space, with biparental inheritance. The two parental lines P1 and P2 are at the lower corners. The initial F1 cross (P1 × P2) will be heterozygous at all divergent sites, and so be found at the top corner. Individuals from later crosses will vary due to segregation and recombination, but the F2 (F1 × F1) will tend to be found toward the center, while backcrosses (such as or F1 × P2) will be found on the upper edges (see Fig. 3 for more details of these backcrosses). Panel (B) illustrates the selection on heterozygosity predicted by Fisher's geometric model, when the parental lines are well adapted (eq. 9). In backcrosses, heterozygosity will be under diversifying selection, favoring both extreme values. By contrast, in the F2, we predict directional selection toward higher heterozygosity. Panel (C) illustrates some complications introduced by heteromorphic sex chromosomes (see eq. 16 and Appendix 4). With XO sex determination, male F1 carry no heterozygosity on the X, which will tend to reduce their fitness, consistent with Haldane's Rule. For male backcrosses (F1 × P1), the selection acting on (autosomal) heterozygosity, will depend on the alleles carried on the X. When the X carries mostly P2 alleles, fitter individuals will be more heterozygous (darker gray arrow). When the X carries mostly P1 alleles, the fittest individuals will carry no heterozygosity (lighter gray arrow). is the proportion of the divergent sites found on the X, and is set at 37%, as we have estimated for Drosophila pseudoobscura.

Figure 3

Data on hybrid vigor (i.e., increased yield), from crosses of inbred Zea mays. The original data were collated by Wright (1977; see his Table 2.23), and Hallauer et al. (2010; see their Table 9.13), including only data from single crosses, where there was hybrid vigor in the F2, and yield measured in quintals per hectare. Panel (A) plots the excess yield of the F2 (eq. 11). Results are shown for variety crosses (black triangles), as well as crosses of inbred lines in the strict sense (all other points). The dashed line shows the prediction of 0.5 from single‐locus theory (eq. 12). Panel (B) shows the four datasets collated by Wright (1977), which allow us to compare the F2 and various backcrosses. These crosses, chosen to yield different levels of heterozygosity, are the parental type (P1), the second backcross (); the first backcross (), the F2 (F1 × F1), second backcross to the other parent , and the F1 (P1 × P2) (The data of Stringfield (1950), shown as gray triangles, replace BC2* with an F2 between two distinct F1, involving three distinct strains, but the predictions are unchanged). The gray symbols for the four datasets correspond to those used in panel (A). The dotted line in panel (B) shows predictions from Fisher's model, assuming that the between‐strain divergence contains limited coadaptation. The prediction uses equations (13), (14), and (11), with , and . The model predicts both the roughly linear increase in vigor with heterozygosity, and the systematic difference between BC1 and F2.

Fisher's geometric model

Fisher's model is defined by n quantitative traits under selection toward an intermediate optimum (Fisher 1930). We will assume that fitness is always measured in a fixed environment, but we make no assumptions about how this optimum might have moved during the parental divergence (see Appendix 1). If the selection function is multivariate normal, including correlated selection, then we can rotate the axes and scale the trait values, to specify n new traits that are under independent selection of different strengths (Waxman and Welch 2005; Martin and Lenormand 2006; Martin 2014). Examples with traits are shown in Figure 2. We now define the breakdown score of a phenotype as where, for trait i, is its deviation from the optimum and is the strength of selection. By assumption, all mutational changes act additively on each trait, but their effects on fitness can vary with the phenotype of the individual in which they appear, and this yields fitness dominance and epistasis (Appendix 1).

Figure 2

Fisher's geometric model generates a flexible series of fitness landscapes for hybrid genotypes. The left‐hand panels show a schematic representation of the model, with traits, each under optimizing selection of differing strengths. (The orientation of the axes shows that the model allows for correlated selection, although this is ignored in the text, by rotating the axes). The right‐hand panels show how the expected breakdown score (eqs. 2 and 6), varies for hybrids of different types. Panels (A)–(C) show different assumptions about levels of parental maladaptation. Panel (A) shows a scenario where all of the parental divergence is maladaptive, with no tendency for their alleles to be coadapted. In this case, hybrid fitness increases with their heterozygosity, as predicted by Wright's single‐locus theory of inbreeding (eq. 8; Wright 1922, 1977). Panel (B) shows a scenario where the parental lines are optimal (or, at least, very well adapted compared to the worst class of hybrid that can be formed between them). In this case, the hybrid index is under symmetrical diversifying selection, and the form of selection on heterozygosity will vary for different cross types (eqs. 9, 13, and 15). This landscape can also be generated by a general model of genetic incompatibilities (see Appendix 2). Panel (C) shows a situation where only one of the parental lines is maladapted (eq. 10 with ).

The d substitutions that differentiate P1 and P2 can be considered as a chain of vectors, which connect the two parental phenotypes (Fig. S1). While the sizes and directions of these vectors will generally be unknown, in Appendix 1, we show that the chain can be approximated as a constrained random walk, or Brownian bridge (Revuz and Yor 1999, Ch. 1). This approximation relies on the fact that hybrid genomes contain the fixed alleles in effectively random combinations, and it gives accurate results for a wide range of assumptions about the divergence process (Figs. S2–S3), including adaptation of the parental lines to different environments (Appendix 1; Fig. S4).

Under the Brownian bridge approximation, the quantity , that appears in equation (2), is given by where is the expected variance contributed to trait i by a fixed mutation in heterozygous state and we have assumed that any maladaptation has accrued independently in the two parental lines (see Appendix 1 and Table S1). The key quantity f is given by where

(see Appendix 1 for full details). Two features of equation (6) are immediately notable. First, it does not depend on any of the model parameters. For example, the number of traits, n, could affect the accuracy of our approximation (since S will tend to approach normality as n increases); but it does not appear in equation (6), which depends solely on p 12, h and the levels of parental maladaptation. Second, for a given value of h, the fitness of hybrids will tend to increase with their heterozygosity, p 12. This prediction agrees with the widespread finding of heterosis (Crow 1952; Table S1 of Fraïsse et al. 2016b). Indeed, by rearranging equation (6) it is clear that hybrids will tend to be fitter than parental line P1 (such that ) on the condition that .

It is also useful to consider equation (6) in a few special cases. First, let us assume that all of the divergence between the parents is maladaptive, without any tendency for coadaptation between their alleles. In this case, the parental phenotypes can be treated as unconstrained random walks away from the optimum. Each fixed mutation, in homozygous state, will contribute an expected to the variance on each trait, and so, from equations (4) and (5), we have (see also Appendix 1). With this high level of parental maladaptation, the expected breakdown score of the hybrids is The fitness landscape implied by equation (8) is illustrated in Figure 2 A. With these extreme levels of parental maladaptation, hybrid breakdown depends only on the heterozygosity. Indeed, this result is closely related to Wright's (1922) single‐locus theory of inbreeding, which was developed to analyze crosses between closely related inbred lines, where all divergence between the parental lines comprises deleterious alleles. Wright's theory therefore appears as a special case of Fisher's geometric model, when parental alleles show no coadaptation (see below).

Now let us consider the other extreme case, where the parental alleles are completely coadapted, such that P1 and P2 both have optimal fitness, but realized by different combinations of alleles. In this case, we find

The fitness landscape implied by equation (9) is illustrated in Figure 2 B. With well‐adapted parents, the hybrid index is under symmetrical diversifying selection.

As with equation (8), equation (9) can also be derived via an alternative route, using a model of speciation genetics. In particular, we show in Appendix 2 that equation (9) can be obtained from a general model of “Dobzhansky–Muller incompatibilities,” each involving a small number of loci (Orr 1995; Turelli and Orr 2000; Gavrilets 2004; Welch 2004; Fraïsse et al. 2016b). The agreement between the two models depends on a particular set of assumptions about the dominance of incompatibilities, namely (i) partial recessivity on average, and (ii) greater reduction in fitness when they contain homozygous alleles from both parental lines. However, we show in Appendix 2 that these assumptions are biologically realistic, and implied by a number of well‐established empirical patterns (Turelli and Orr 2000). The result is that equation (9) can be interpreted as a model of genetic incompatibilities, but without the large number of free parameters that these models can require.

Equation (9) is expected to apply to many cases of wild hybridization, because it should provide a good approximation even if the parental populations are not truly optimal. The only requirement is that they be much better adapted than the worst possible class of hybrid, such that .

Nevertheless, the general form of equation (6) also allows us to explore intermediate cases. For example, if only P2 is maladapted, then we find In this case, as illustrated in Figure 2 C, selection on hybrid index is skewed toward the fitter parent. Below, we will show how all three of these special cases (Fig. 2 A–C) can be used in data analysis.

TESTING THE PREDICTIONS WITH BIPARENTAL INHERITANCE

Fitness differences between crosses

The simplest predictions from equation (6) assume standard biparental inheritance at all loci. In this case, the standard cross types can be easily located on our fitness landscapes. This is shown in Figure 1 A.

With well‐adapted parental lines (Fig. 2 B) hybrids of high fitness are expected only for the initial F1 cross (P1 × P2) and breakdown is predicted for later crosses, such as the first backcross (BC1: F1 × P1 or F1 × P2) or the F2 (F1 × F1). This pattern of F1 hybrid vigor (high F1 fitness) followed by breakdown in later crosses, has widespread empirical support (see references in Table S1 of Fraïsse et al. 2016b and Rosas et al. 2010).

With highly maladapted parents, by contrast, hybrids of all types can be fitter than their parents (see Fig. 2 A). Plentiful data of this kind come from highly inbred lines of Zea mays (Neal 1935; Wright 1977; Melchinger 1987; Hinze and Lamkey 2003; Hallauer et al. 2010). To analyze these data, a widely used proxy for fitness is the excess yield of a cross (i.e., its increase in yield relative to the parental lines), scaled by the excess yield of the F1. From equations (1)–(2), and using a Taylor expansion, we find where the subscript P refers to both parental lines, which are assumed to have similar yields. For later crosses, these values will vary between individuals within a cross, due to segregation and recombination, but in this section we ignore this variation, and assume that p 12 and h take their expected values for a given cross type (Fig. 1 A). A fuller treatment is outlined in Appendix 3.

Equation (12) confirms that Fisher's model can produce Wright's (1922) single‐locus predictions for inbreeding, but only when all divergence between the parental lines is deleterious (), and when increases in breakdown score act independently on log fitness (). These single‐locus predictions have strong support in Zea mays (Neal 1935; Wright 1977; Melchinger 1987; Hinze and Lamkey 2003; Hallauer et al. 2010). For example, as shown in Figure 3 A, the excess yield of the F2 is often around 50%, which is equal to its expected heterozygosity (Wright 1977; Hallauer et al. 2010). It is also notable that the two outlying points (from Shehata and Dhawan 1975), are variety crosses, and not inbred lines in the strict sense.

Despite this predictive success, Wright (1977) noted a pattern that single‐locus theory could not explain. In Wright's words: “the most consistent deviation from expectation [...] is the low yield of F2 in comparison with the first backcrosses” (Wright 1977, p. 39). Because for both crosses, this cannot be explained without epistasis, or coadaptation between the alleles in the parental lines. In fact, the pattern is predicted by Fisher's model, when there is a small amount of coadaptation, such that . This yields a fitness surface with a small amount of curvature, which is intermediate between those shown in Figure 2 A and 2 B.

Figure 3 B plots the four relevant datasets collated by Wright, and compares the results to predictions from equations (6) and (11), with and . The model predicts the roughly linear increase in yield with mean heterozygosity, as with single locus theory, but also predicts the consistent difference between BC1 and the F2.

Selection on heterozygosity within crosses

In the previous section, we ignored between‐individual variation in heterozygosity within a given cross type. In this section, we show how natural selection is predicted to act on this heterozygosity.

First, let us consider the F2. In this case, we have with relatively little variation between individuals (see Appendix 3 for details). Therefore, if both parents have similar levels of maladaptation, equation (6) is well approximated by The prediction is that p 12 will be under directional selection in the F2, favoring individuals with higher heterozygosity. This is illustrated in Figure 1 B.

Now let us consider a backcross: F1 × P1. In this case, all of the homozygous alleles must come from P1, such that , and so equation (6) becomes So backcrosses will tend to be least fit when they have intermediate levels of heterozygosity. When parents are well adapted, heterozygosity is under symmetrical disruptive selection, favoring heterozygosities that are either higher or lower than (eq. 15). This is illustrated in Figure 1 B .

To test these predictions, we used a new set of genetic data from hybrids of the mussel species: Mytilus edulis and Mytilus galloprovincialis (Bierne et al. 2002, 2006). These species fall at the high end of the continuum of divergence during which introgression persists among incipient species (Roux et al. 2016). We used experimentally bred F2 and BC1, with selection imposed implicitly, by the method of fertilization, and by our genotyping only individuals who survived to reproductive age (Bierne et al. 2002, 2006; see Methods and Fig. S6 for full details).

To estimate heterozygosity in each hybrid individual, we used the 43 markers that were heterozygous in all of the F1 hybrids used to make the subsequent crosses (see Fig. S6). We then asked whether the distribution of p 12 values in recombinant hybrids was symmetrically distributed around its Mendelian expectation of , or whether it was upwardly biased, as would be expected if directional selection were acting on heterozygosity. As shown in the first column of Table 2, Wilcoxon tests found that heterozygosities in surviving hybrids were significantly higher than expected, in both the F2 and backcross. These results may have been biased by the inclusion of individuals with missing data, because they showed higher heterozygosity (see Table S2). We therefore repeated the test with these individuals excluded. As shown in the second column of Table 2, results were little changed, although the bias toward high heterozygosities was now weaker in the backcross.

Table 2

Tests for selection on heterozygosity in F2 and Backcrosses of Mytilus mussels

Markers:	43	43	33	33
Dataset:	All	No missing data	No missing data	Subsampled
Cross	p^_12 (N) P‐value	p^_12 (N) P‐value	p^_12 (N) P‐value	p^_12 (N) P‐value
F2	0.57 (132) 1.5×10^−6 ***	0.56 (88) 6.4×10^−4 ***	0.55 (91) 0.0033**	0.56 (56) 0.0020**
BC	0.57 (144) 1.3×10^−5 ***	0.53 (94) 0.0282*	0.53 (105) 0.0569	0.52 (56) 0.5815

, the estimated median heterozyosity; N, the number of hybrid individuals sampled; P‐value, result of a Wilcoxon test of the null hypothesis median (* ; ** ; *** ). F2, random mating of F1 between M. galloprovincialis and M. edulis; BC, Backcross of the F1 to M. galloprovincialis. No missing data, all individuals with missing data for any of the markers were excluded; Subsampled, for the BC, any individual carrying a marker that was homozygous for the major allele carried by wild M. edulis populations was excluded; for the F2, we downsampled by sequencing order to equalize sample sizes.

Interpreting these results is complicated by the ongoing gene flow between M. edulis and M. galloprovincialis in nature (Fraïsse et al. 2016a). To test for this, we genotyped 129 pure‐species individuals, and repeated our analyses with a subset of 33 markers that were strongly differentiated between the pure species (see Methods, Fig. S6 and Table S3 for details). With these markers, there was evidence of elevated heterozygosity in the F2, but not the backcross (Table 2 third column). We also noticed that many of our backcross hybrids, though backcrossed to M. galloprovincialis, carried homozygous alleles that were typical of M. edulis. We therefore repeated our analysis after excluding these “F2‐like” backcrosses. Results, shown in the fourth column of Table 2, showed that the reduced BC dataset showed no tendency for elevated heterozygosity. However, the bias toward higher heterozygosities remained in the F2, even when we subsampled to equalize the sample sizes.

Despite the problems of interpretation due to introgression and shared variants, the results support the prediction of equations (13)–(15): that directional selection on heterozygosity should act in the F2, but weakly or not at all in the backcross.

PREDICTIONS OF FISHER'S GEOMETRIC MODEL WITH SEX‐SPECIFIC INHERITANCE

Results above assumed exclusively biparental inheritance. But the predictions of Fisher's model are easily extended to include heteromorphic sex chromosomes, or loci with strictly uniparental inheritance, such as organelles (Coyne and Orr 1989; Turelli and Orr 2000; Turelli and Moyle 2007; Fraïsse et al. 2016b). In this section, we introduce the approach, and demonstrate its flexibility, while the full derivations are collected in Appendix 4. The basic method of introducing sex‐specific inheritance is to write p 12 and h as weighted sums of contributions from different types of locus. For example, to focus on the contribution of the X chromosome versus the autosomes, we can write where the subscripts denote the chromosome type (so that is the proportion of divergent sites on the autosomes that are heterozygous), and and are weightings, which sum to one (Turelli and Orr 2000).

We can now predict differences in hybrids of different sexes. For example, previous authors have shown that Fisher's model predicts Haldane's Rule: the tendency of sex‐specific breakdown to appear in the heterogametic sex (Haldane 1922; Turelli and Orr 2000; Barton 2001; Fraïsse et al. 2016b; Schiffman and Ralph 2017; see Fig. 1 C). Using equations (16) and (6), this basic insight can easily be extended to give formal conditions for Haldane's Rule, with different levels of parental maladaptation, and uniparental inheritance (see Appendix 4).

Sex chromosomes also complicate predictions about selection on heterozygosity within crosses. Consider, for example, male backcrosses in an XY species, where the X chromosome is large, the Y chromosome is small or degenerate, and the parental lines are reasonably fit. These backcross males will vary in their autosomal heterozygosity (), but the selection on this heterozygosity will vary with the alleles carried on the X. This is illustrated in Figure 1 C. If the backcross is made to P1, but is large (i.e. divergent X‐linked alleles come mainly from P2), then heterozygosity will tend to be under positive selection (see darker arrow in Fig. 1 C); but if is small, then heterozygosity will tend to be under negative selection (see lighter arrow in Fig. 1 C, Appendix 4 and Fig. S8 for details).

To test this prediction, we used data from Noor et al. (2001) (Fig. S7). These authors generated male backcrosses between Drosophila pseudoobscura and Drosophila persimilis. These species are ideal for our test, because the X carries of the coding sequence, the Y is degenerate, and the hybrids can be fully sterile (see Table 1, Methods, and Noor et al. 2001). Each backcross male was scored for sterility, and multiple genetic markers, including two markers on the X (see Table 1). We characterized the X as low‐ if both markers carried the P1 allele, and as high‐ if both markers carried the P2 allele; we excluded intermediate individuals, for which we have no clear prediction. We then regressed sterility on heterozygosity, , but allowed the intercept and the slope to vary with . Results, shown in Table 3, show that both and are important predictors of sterility in these backcrosses, and that their effects interact. With low , sterility is generally rarer, but positively correlated with . With high , by contrast, sterility is more common, and negatively correlated with . The predictions of Fisher's model are therefore supported, in both cross directions.

Table 3

Regressions of male sterility on autosomal heterospecificity in Drosophila backcrosses

Backcross to	N	Model	intercept (low h_X)	intercept (high h_X)	Best‐fit coefficients for p_12,A	AIC
D. persimilis	577	two intercepts	−1.071	1.797	–	596.79
		two intercepts + one slope	−2.165	0.764	2.147	580.41
		two intercepts + two slopes	−3.033	2.871	3.746 (low‐h_X); −1.973 (high‐h_X)	558.91
D. pseudoobscura	610	two intercepts	−0.197	2.962	–	603.53
		two intercepts + one slope	−2.031	1.219	3.620	558.18
		two intercepts + two slopes	−2.485	4.034	4.505 (low‐h_X); −1.876 (high‐h_X)	545.87

AIC, Akaike Information Criterion; preferred model shown in bold.

The interaction between and may also help to explain a broader set of patterns observed by Moehring (2011) in multiple datasets of Drosophila backcrosses (Macdonald and Goldstein 1999; Noor et al. 2001; Moehring et al. 2006a, 2006b and see Table 1). Moehring found that male fertility correlated strongly and negatively with , while correlations with were weak and inconsistent. This pattern was not consistent with predictions of a simple model of X‐autosome incompatibilities (Moehring 2011), but it is consistent with Fisher's geometric model (see Fig. 1 C, Fig. 2 B, and Appendix 4).

Female backcrosses in Teleogryllus

Fisher's model can help to explain more complex patterns. To illustrate this, we will consider the data of Moran et al. (2017) from the field crickets Teleogryllus oceanicus and T. commodus. These species have XO sex determination, and a large X chromosome (, from Moran et al. 2017). They are also a rare exception to Haldane's Rule, with F1 sterility appearing solely in XX females (Hogan and Fontana 1973). Moran et al. (2017) hypothesized that female sterility might be caused by negative interactions between heterospecific alleles on the X, which appear together in , but not in .

To test this hypothesis, they generated backcross females with identical recombinant autosomes, but different levels of heterozygosity on the X, . The hypothesis of X–X incompatibilities predicts that fertility will decrease with , but this was not observed. This is shown in Figure 4 A.

Figure 4

Data on female fertility, from crosses of the field crickets Teleogryllus oceanicus (P1), and T. commodus (P2), from Moran et al. (2017), compared to theoretical predictions from Fisher's geometric model. Plotting styles denote the level of X‐linked heterozygosity: high (circles and darker shading); intermediate (triangles and lighter shading) or low (squares and no shading). To plot the data (column A), we used the mean number of offspring per pair (upper panel), or mean number of eggs per pair (lower panel), each normalized by the value for P1 (see the supplementary information of Moran et al. 2017 for full details). The theoretical predictions (column B) are listed in Table S4. In the upper panel, we assume , as estimated from the chromosome sizes, and complete silencing of the paternal X (such that in Table S4). In the lower panel, we assume , and incomplete silencing of the paternal X (), to improve fit to the egg data. While predictions apply to the rank order of fitnesses, to aid visualization, we plot (see eq. 1), and set the parameter f P2 via , to reflect the lower fertility of this species under lab conditions.

To try and explain the patterns that were observed, let us note two clear asymmetries in the data. First, there are strong differences between the fitness levels of the two parental species in lab conditions, with T. commodus (labeled P2 in Fig. 4) producing around half the eggs and offspring of T. oceanicus (P1 in Fig. 4). This suggests that equation 10 might apply to these data. Figure 4 shows a second asymmetry in the data: the reciprocal F1 (female–male vs. male–female) have very different fitness (Turelli and Moyle 2007; Fraïsse et al. 2016b). Because the X and autosomes will be identical for both cross directions, this implies some sort of parent‐of‐origin effect on the phenotype. This could be either uniparental inheritance, or selective silencing (Turelli and Moyle 2007; Fraïsse et al. 2016b). One possibility is the speculation of Hoy and collaborators (see e.g., Hoy et al. 1977; Butlin and Ritchie 1989; Dr. Peter Moran pers. comm.), that dosage compensation in Teleogryllus involves silencing of the paternal X.

We can now use the foregoing assumptions to predict the levels of breakdown for each of the crosses produced by Moran et al. (2017). The predictions for each cross are listed in Table S4, and plotted in the upper panel of Figure 4 B. This simple model explains several striking aspects of the Teleogryllus data.

Adjusting the parameter values can further improve the fit. For example, if we increase (as would be case if divergent sites affecting female fecundity were clustered on the X), and assume that paternal X silencing is incomplete (affecting 80% of the divergent sites), then the results, shown in the lower panel of Figure 4 B, agree well with the data on Teleogryllus egg number, as shown in the lower‐left‐hand panel (only the high fitness of P2 × F112 is poorly predicted).

Further adjustments are possible, but these soon become ad hoc, at least without further knowledge of the true nature of parent‐of‐origin effects in Teleogryllus. The important point here is that Fisher's geometric model explains several key features of these hybrid data, while using only a single parameter derived from the data themselves; and even this parameter (f P2), was estimated from the parental control lines.

ESTIMATING THE FITNESS SURFACE

Across a diverse collection of hybrids, equation (6) predicts that the hybrid index will be under disruptive selection, and heterozygosity under directional selection. This prediction can be tested with datasets containing estimates of fitness, h and p 12 for many hybrid individuals. Exactly such an analysis was presented by Christe et al. (2016), for families of wild hybrids from the forest trees, Populus alba and P. tremula (Lindtke et al. 2012, 2014; Christe et al. 2016). These authors scored survival over four years in a common‐garden environment, and fit a generalized linear model to their binary data (binary logistic regression, with “family” as a random effect), and predictors including linear and quadratic terms in p 12 and h. Model selection favored a four‐term model, with terms in p 12, h, and h 2 (see Table S6, and Supplementary information of Christe et al. 2016 for full details). For comparison with our theoretical predictions, we can write their best fit model in the following form: where y is the fitted value for hybrid breakdown. In its general form, equation (6) should hold and Fisher's model predicts that , , and . Additionally if parental maladaptation is not highly asymmetrical then it predicts .

The best‐fit model of Christe et al. (2016) corresponds to , and , which supports the predictions of directional selection toward higher heterozygosity, and near‐symmetrical diversifying selection on the hybrid index.

To obtain confidence intervals on these parameters, we fit equation (17) to the raw data of Christe et al. (2016). We also searched for other data sets, from which we could estimate the hybrid fitness surface. After applying some quality controls (see Methods and Table S2), we identified one other dataset of wild hybrids, from the mouse subspecies Mus musculus musculus/domesticus, where male testes size was the proxy for fertility (Turner and Harr 2014). We also found four datasets of controlled crosses: F2 from the same mouse subspecies (White et al. 2011), and the ragworts Senecio aethnensis and S. chrysanthemifolius (Chapman et al. 2016); and the Drosophila backcrosses discussed above (Macdonald and Goldstein 1999; Moehring et al. 2006a, 2006b). Unlike the data from wild hybrids, these single‐cross datasets leave large regions of the fitness surface unsampled (see Fig. 1); nevertheless, they each contain enough variation in h and p 12 for a meaningful estimation. Details of all six datasets are shown in Table 1, and they are plotted in Figures S9–S11.

Figure 5 shows a summary of the estimated parameters, and full results are reported in Tables S5 and S6, and Figures S9–S11. Taken together, the results show good support for the predictions of equation (6). For all six datasets there was evidence of significant positive selection on heterozygosity ( was preferred in all cases). Furthermore, for all six datasets, we inferred diversifying selection acting on the hybrid index. Estimates of β2, shown in the upper panel of Figure 5, show that this selection was near‐symmetrical in all cases, such that . The poorest fit to the predictions was found for the Drosophila backcrosses, where estimates of β1 were significantly greater than the predicted upper bound of (Fig. 5 lower panel). But these datasets were least suited to our purpose, because estimates of h and p 12 depend strongly on our rough estimate of the relative contributions of the X and autosomes (see Methods), and because they lack F2‐like genotypes, from the center of the fitness surface (Fig. 1 A; Fig. S11). By contrast, results for the Mus musculus F2 (White et al. 2011), are remarkably close to the predictions of equation (9) (Fig. 5; Fig. S10).

Figure 5

Best fit parameter estimates for the GLM of equation (17), with fitness and genomic data from six sets of hybrids (see Table 1 for details). The upper panel shows estimates of the coefficient β2 that determines the form of selection acting on the hybrid index, h. Estimates of are consistent with symmetrical diversifying selection. The lower panel show estimates of the coefficient β1 that determine the relative strength of selection acting on the hybrid index. Estimates of are predicted when the parental types are well adapted (eq. 9), while estimates are predicted when the parental types are maladapted (eq. 6). Confidence intervals are defined as values that reduce the AIC by 2 units. These measures of uncertainty were not obtained for the F2 data, where variation in the hybrid index contributed little to the model fit, as predicted by equation (13). Full details of the model fitting are found in the Methods and Tables S5 and S6.

Two other features of the results deserve comment. First, for the two F2 datasets, it was not possible to provide meaningful confidence intervals for β1 and β2. This is because, for these two datasets, the terms in h and h 2 did not make a significant contribution to model fit, and so the preferred model contained only selection on p 12 (see Table S6). This is consistent with our earlier prediction of equation (13), and stems from the low variation in among F2 hybrids (see Appendix 3 and Figs. S9 and   S10).

Second, for two of the datasets, Populus and Senecio, the estimates of β1 are substantially lower than 4 (Fig. 5; Fig. S9). This is suggestive of parental maladaptation, creating true heterosis in the hybrids (see eq. 6). Consistent with this inference, there is independent evidence of parental load and F1 hybrid vigor in both species pairs (Populus: Caseys et al. 2015; Christe et al. 2017; Senecio: Abbott and Brennan 2014).

Discussion

In this article, we have used Fisher's geometric model to develop predictions for the relative fitness of any class of hybrid. The modeling approach is simple, with few free parameters, and it generates a wide range of testable predictions. We have tested some of these predictions with new and published datasets (Table 1), and the major predictions of the model are well supported.

We emphasize that our approach is designed for coarse‐grained patterns in the data, and typical outcomes of hybridization, without considering the particular set of substitutions that differentiate the parental lines, or the particular combination of alleles in an individual hybrid. The limitations of such an approach are seen in the low r 2 values associated with our model fitting (Table S5). Nevertheless, our approach should enable novel and complementary uses of genomic datasets, which do not depend on identifying individual loci with anomalous effects. In this approach, all alleles, and not just a handful of “speciation genes” might contribute to hybrid fitness.

A second goal of the present work was to show how Fisher's model can interpolate between previous modeling approaches, namely the classical theory of inbreeding (Wright 1922; Crow 1952), and models of genetic incompatibilities, each involving a small number of loci (Dobzhansky 1937; Orr 1995; Gavrilets 2004; Welch 2004). We have also shown that Fisher's model can account for empirical patterns that each approach has struggled to explain (Wright 1977; Moehring 2011; Moran et al. 2017).

Though we have stressed their similarities, we should also stress that Fisher's model and the model of incompatibilities remain different in two ways. First, as shown in Appendix 2, the two models agree only when the dominance relations of incompatibilities take a particular set of values (eq. A36)—albeit a biologically realistic set (Appendix 2; Turelli and Orr 2000). Second, even when predictions are identical for the quantity f (eq. 2), the two approaches still make different predictions for other kinds of data, which were not considered in the present work. The most important difference is the dependency of log fitness on d, the genomic divergence between the species. Under Fisher's geometric model, the log fitness of hybrids declines with (eqs. 1–2 and 5–6). By contrast, with the simplest models of incompatibilities, log fitness declines with where ℓ is the number of loci involved in each incompatibility (Orr 1995; Welch 2004; Appendix 2, eq. A22). This is due to a snowball effect, where the number of incompatibilities grows explosively with . This is a genuine difference between the modeling approaches, although truly discriminatory tests may be difficult. For example, it may not always be possible to distinguish between a snowballing model with a low value of α (equivalent to strong positive epistasis between incompatibilities), or a model where α is higher, but where the number of “incompatibilities” does not snowball, because they appear and disappear as the genetic background changes (Welch 2004; Fraïsse et al. 2016b; Guerrero et al. 2017; Kalirad and Azevedo 2017). Furthermore, Fisher's model can generate an “apparent snowball effect,” when introgressed segments contain multiple divergent sites (Fraïsse et al. 2016b), and this is true of the strongest empirical demonstrations of the effect (Matute et al. 2010; Moyle and Nakazato 2010; Wang et al. 2015).

Finally, given the simplicity and flexibility of the modeling approach explored here, and its predictive successes with a range of data, it should be readily extendable to address other outstanding questions in the study of hybridization. These include the putative role of hybridization in adaptive evolution (e.g., Mendez et al. 2012; Fraïsse et al. 2016b, 2016a; Duranton et al. 2017), the effects of recombination in shaping patterns of divergence (Aeschbacher et al. 2017; Schumer et al. 2018), the analysis of clines (Barton and Gale 1993; Macholán et al. 2011; Baird et al. 2012), and the roles of intrinsic versus extrinsic isolation (Chevin et al. 2014, Fig. S4). Given its ability to interpolate between models of different and extreme kinds, it should also be particularly useful for understanding hybridization in intermediate regimes, where parental genomes are characterized by both maladaptation and allelic coadaptation, or where the architecture of isolation involves many genes of small or moderate effect. Data—including those analyzed here—suggest that such architectures might be quite common (Davis and Wu 1996; Maside and Naveira 1996; Rosas et al. 2010; Morán and Fontdevila 2014; Baird 2017; Buerkle 2017; Boyle et al. 2017).

Methods

Mytilus DATA

Conserved tissues from the mussel species, Mytilus edulis and Mytilus galloprovincialis, and their hybrids, were retained from the work of Bierne et al. (2002, 2006). As reported in those studies, M. edulis from the North of France were crossed with M. galloprovincialis from the French Mediterranean coast to produce F1 hybrids (five males and one female; Bierne et al. 2002). The F1 were then used to produce an F2, and sex‐reciprocal backcrosses to M. galloprovincialis (which we denote here as BC12 and BC21; Bierne et al. 2006). In particular, oocytes from the F1 female were fertilized by the pooled sperm of the five F1 males producing F2 individuals, from which 132 individuals were sampled; oocytes from the F1 female were fertilized by pooled sperm of five M. galloprovincialis males to produce BC12, from which 72 individuals were sampled; and five M. galloprovincialis females were fertilized by pooled sperm from the five F1 males, producing BC21, from which 72 individuals were sampled. In addition to these hybrids, we also genotyped 129 individuals from “reference” populations of the two species, found in regions with relatively little contemporary introgression. In particular, we sampled M. galloprovincialis from Thau in the Mediterranean Sea; and sampled M. edulis from four locations in the North Sea and English Channel (The Netherlands, Saint‐Jouin, Villerville and Réville). Full details of these reference populations are found in Table S3.

In each case, gill tissues were conserved in ethanol at −20°C. DNA was extracted using a NucleaMag 96 Tissue kit (Macherey‐Nagel) and a KingFisher™ Flex (ThermoFisher Scientific). We then genotyped all samples for 98 Mytilus markers that were designed from the data of Fraïsse et al. (2016a). The flanking sequences of the 98 SNPs are provided in Table S7. Genotyping was subcontracted to LGC‐genomics and performed with the KASP™ chemistry (Kompetitive Allele Specific PCR, Semagn et al. 2014). Results are shown in Figure S6. Many of the 98 markers are not diagnostic for M. edulis and M. galloprovincialis, and so we retained only the 43 that were scored as heterozygous in all six of the F1 hybrids. To obtain a reduced set of strongly diagnostic markers, we measured sample allele frequencies in our pure species M. edulis and M. galloprovincialis samples (pooling M. edulis individuals across the four sampling locations; Table S3), and retained only markers for which the absolute difference in allele frequencies between species was > 90%. This yielded the set of 33 markers used for the right‐hand columns in Table 2. The “subsampled” data shown in the fourth column of Table 2, excluded any BC hybrid who carried the major allele typical of M. edulis in homozygous form. This yielded 56 BC hybrids. Because sequencing order was haphazard, to equalize the sample sizes, we retained the first 56 F2 to be sequenced.

COLLATION OF PUBLISHED DATA

We searched the literature for published datasets combining measures of individual hybrid fitness, with genomic data that could be used to estimate p 12 and h. We rejected any datasets where the measure of fecundity or fertility took an extreme low value for one of the pure species, suggesting that it is not a good proxy for fitness (e.g., Orgogozo et al. 2006), data where the fitness proxy correlated strongly with a measure of genetic abnormality such as aneuploidy (Xu and He 2011), or data where the states of many markers could not be unambiguously assigned, for example, due to shared variation. Before estimating the fitness surface, we also excluded any dataset where there was a highly significant rank correlation between the proportion of missing data in an individual, and either their heterozygosity, or fitness. For this reason, we did not proceed with reanalyses of the excellent datasets of Li et al. (2011), or Routtu et al. (2014) (see Table S2 for full details). For our reanalysis of the Teleogryllus data of Moran et al. (2017), we did not consider data from the second backcross, to avoid complications due to selection on heterozygosity during the first backcross, and because of an anomalous hatching rate in the T. oceanicus controls. For our reanalysis of the Mus musculus F2 (White et al. 2011), we used a conservative subset of these data; we excluded any individual where any X‐linked marker was scored as heterozygous (indicative of sequencing errors in heterogametic males; White et al. 2011), and controlled for variation in the uniparentally inherited markers, by retaining only individuals carrying M. m. domesticus mitochondria, and the M. m. musculus Y. However, results were little changed when we used all 304 individuals with fertility data (Table S6). Results were also unaffected when we used alternative proxies for fitness (Table S6). For the wild Mus musculus hybrids data (Turner and Harr 2014), fixed markers and their orientation between species were determined using previously published data of Staubach et al. (2012) and Yang et al. (2009). See details in the Dryad data package.

ESTIMATION OF FROM ANNOTATED GENOMES

To fit the generalized linear models (GLM) to data from species with XY sex determination (Table 1), we needed to account for the different marker densities on the X and autosomes. Accordingly, we estimated the overall values of p 12 and h from equation (16). We estimated the weightings, and , from the total length of coding sequences associated with each chromosome type, ignoring the small contributions from the Y and mitochondria. In each case, we obtained the longest protein product for each unique gene, and then summed their lengths, using a custom R script. The values, shown in Table 1, were calculated as the total length of X‐linked CDS divided by the total CDS length. For Mus musculus, we used the reference genome assembly “GRCm38.p5.” For Drosophila simulans, we used the assembly “GCA_000754195.3 ASM75419v2,” and for Drosophila yakuba “GCA_000005975.1 dyak_caf1.” For Drosophila pseudoobscura, the current annotation was downloaded from FlyBase release 3.04 (Gramates et al. 2017). The .gtf file was then sorted and merged (combining overlapping coding sequences on each chromosome) using BEDTools (Quinlan and Hall 2010). Coding sequence lengths were calculated and summed over each chromosome, using custom awk commands.

GLM METHODS

The linear model results shown in Table 3, Figure 5, Tables S5 and S6, and Figures S9–S11, were all fit in R v. 3.3.2 (R Core Team 2016). For datasets with quantitative fitness measures (White et al. 2011; Turner and Harr 2014; Fig. S10) we used the standard general linear models, with Gaussian errors, and chose data transformations to reduce heteroscedasticity. For binary fitness data (Noor et al. 2001; Christe et al. 2016; Chapman et al. 2016; Table 3; Fig. S9), we used binomial regression with a logit link implemented in the function; and with ordinal fitness data (Macdonald and Goldstein 1999; Moehring et al. 2006b; Fig. S11) we used proportional odds logistic regression (Agresti 2003), implemented in the polr function. In these cases, the P‐values shown in Table S6 were calculated by comparing the t‐value to the upper tail of normal distribution, as in a Wald test. For the non‐Gaussian models, we also report McFadden's pseudo‐r 2, defined as one minus the ratio of log likelihoods for the fitted and null models (McFadden 1974).

DATA ACCESSIBILITY

The newly generated data on Mytilus and collated data are available as a Dryad package (doi:xxx).

Associate Editor: Z. Gompert

Figure S1: A schematic representation of a phenotype evolving over time in two populations, labeled P1 and P2, starting from their most recent common ancestor (MRCA).

Figure S2: The breakdown associated with haploid hybrid genotypes, under simple models of phenotypic divergence.

Figure S3: The breakdown associated with haploid hybrid genotypes, under explicit population genetic simulations of phenotypic divergence.

Figure S4: The breakdown associated with haploid hybrid genotypes, under explicit population genetic simulations of phenotypic divergence, in scenarios involving discrete jumps in the optimal value for one of n = 2 traits.

Figure S5: The effects of an incompatibility on hybrid breakdown score, (eq. 41), when the incompatibility appears in a genotype comprising i loci that are homozygous for alleles from one parental species, j loci that are homozygous for alleles from the other parental species, and k loci that are heterozygous.

Figure S6: Genotype plot for the raw Mytilus data.

Figure S7: Plots of the Drosophila male backcross data reanalyzed here (see Table 1 and Moehring 2011).

Figure S8: Predictions of Fisher's geometric model for heterogametic male hybrids.

Figure S9: Estimation of the fitness surface for interspecific hybrids from plants (Table 1), namely wild hybrids of the genus Populus (row (a); Christe et al. 2016), and an F2 cross of the genus Senecio (row (b); Chapman et al. 2016).

Figure S10: Estimation of the fitness surface for subspecific hybrids from Mus musculus (Table 1, White et al. 2011; Turner and Harr 2014).

Figure S11: Estimation of the fitness surface for backcross male hybrids from Drosophila species pairs (Table 1, Macdonald and Goldstein 1999; Moehring et al. 2006a, b).

Table S1: The contribution of parental maladaptation to hybrid breakdown.

Table S2: Checks for appropriateness for genomic data sets.

Table S3: The reference populations for Mytilus crosses.

Table S4: Expected breakdown scores for homogametic female hybrids with paternal X silencing.

Table S5: Inferring the hybrid fitness surface from genomic data sets.

Table S6: The significance of individual regression coefficients.

Table S7: Information on the 98 markers used for the Mytilus genotyping.

Click here for additional data file.