Analysis of ancestry heterozygosity suggests that hybrid incompatibilities in threespine stickleback are environment dependent.

Hybrid incompatibilities occur when interactions between opposite ancestry alleles at different loci reduce the fitness of hybrids. Most work on incompatibilities has focused on those that are "intrinsic," meaning they affect viability and sterility in the laboratory. Theory predicts that ecological selection can also underlie hybrid incompatibilities, but tests of this hypothesis using sequence data are scarce. In this article, we compiled genetic data for F2 hybrid crosses between divergent populations of threespine stickleback fish (Gasterosteus aculeatus L.) that were born and raised in either the field (seminatural experimental ponds) or the laboratory (aquaria). Because selection against incompatibilities results in elevated ancestry heterozygosity, we tested the prediction that ancestry heterozygosity will be higher in pond-raised fish compared to those raised in aquaria. We found that ancestry heterozygosity was elevated by approximately 3% in crosses raised in ponds compared to those raised in aquaria. Additional analyses support a phenotypic basis for incompatibility and suggest that environment-specific single-locus heterozygote advantage is not the cause of selection on ancestry heterozygosity. Our study provides evidence that, in stickleback, a coarse-albeit indirect-signal of environment-dependent hybrid incompatibility is reliably detectable and suggests that extrinsic incompatibilities can evolve before intrinsic incompatibilities.

Introduction

Hybrid incompatibilities—interactions among divergent genetic loci that reduce the fitness of hybrids—are a key component of reproductive isolation between diverging lineages [1]. Incompatibilities have been studied most intensively in the context of sterility and mortality, in part because phenotyping in the laboratory is reliable and because such incompatibilities can have simple genetic architectures involving few loci [2,3]. These sorts of interactions have come to be called “intrinsic” hybrid incompatibilities due to the fact that there are conflicts within the hybrid genome that are expected to impact hybrids in most environmental contexts (although note that the strength of selection against some intrinsic incompatibilities can vary across environments [4,5]). Studies have shown that the number of loci involved in intrinsic incompatibilities tends to increase with genetic divergence between the parent species [6-8] and that incompatibilities can be common throughout the genomes of isolated conspecific populations [9,10]. Collectively, evolutionary biologists have made substantial progress toward identifying generalities about the evolutionary genetics of intrinsic hybrid incompatibilities.

Ecological selection could underpin incompatibilities if particular allele combinations render hybrids unable to function in their ecological environment, such as in avoiding predators or capturing prey. Several recent studies have shown patterns consistent with this effect, wherein hybrids have “mismatched” trait combinations and reduced fitness as a result [11-13]. Such studies have successfully demonstrated that incompatibilities due to trait mismatch exist, but links to the underlying genetics have not been made. Perhaps the most significant barrier to progress in studying the genetics of “ecological” hybrid incompatibilities is the unique difficulty of detecting them. The ability to detect individual incompatibilities depends on the genetic architecture underlying traits—that is, whether quantitative trait loci (QTL) have small or large phenotypic effects [14,15]. For example, Arnegard and colleagues [11] found that combining divergent jaw traits together in F2 threespine stickleback fish (Gasterosteus aculeatus L.) hybrids likely reduced their fitness because these traits interacted in a manner that reduced suction feeding performance. The interacting jaw components map to several regions of the genome that individually explain a small fraction (<10%) of the phenotypic variance (and most variance was unexplained), thus rendering it difficult to study their individual epistatic fitness effects.

Recent theoretical advances, however, suggest ways to test for and measure the net effect of hybrid incompatibilities using experimental crosses. Specifically, selection against hybrid incompatibilities in an F2 hybrid cross causes an increase in ancestry heterozygosity—the number of sites in the genome that carry both parents’ alleles at ancestry informative sites—at the population level [16,17]. This is expected because F2 hybrids have a hybrid index of approximately 0.5—deriving half of their alleles from one parental species and half from the other. Thus, individuals with high heterozygosity relative to their hybrid index have fewer pairs of homozygous loci with opposite ancestry compared to relatively more homozygous individuals with a similar hybrid index. Assuming that most alleles affect the phenotype additively and have noninfinitesimal effect sizes, having many loci with opposite homozygous ancestry can result in hybrids with maladaptive “mismatched” phenotypes, whereas highly heterozygous individuals are expected to have less mismatched phenotypes (Fig 1, S1 Fig). Whether “mismatch” affects fitness, however, ultimately depends on the ecology of the system and the underlying fitness landscape. Such coarse approaches—coarse because they use summary statistics rather than direct mapping—are a promising means to identify the presence of small effect hybrid incompatibilities at the genetic level using field experiments.

Fig 1

Results from simulations illustrating how an ecological mechanism could underlie the heterozygosity–incompatibility relationship in F2 hybrids.

Both panels depict results from a representative simulation run of adaptive divergence and hybridization between 2 populations. We consider an organism with 2 traits that have both diverged as a result of selection. Colored points are individual hybrids, with darker colors indicating higher heterozygosity. Panel (A) depicts the distribution of 500 F2 hybrid phenotypes in two-dimensional trait space. Large black points are the 2 parent phenotypes, which are connected by a black line indicating the “axis of divergence.” Panel (B) depicts the relationship between individual excess ancestry heterozygosity and trait “mismatch” of individual hybrids [13]. Excess ancestry heterozygosity is the observed heterozygosity minus the expected heterozygosity based on ancestry proportion—0 is the expected mean in the absence of selection (approximately observed heterozygosity frequency of 0.5). Mismatch is calculated as the shortest (i.e., perpendicular) distance between a hybrid’s phenotype and the black line connecting parents in (A). Variation parallel to this axis connecting parents in (A) captures variation in the “hybrid index.” The plot shows that trait mismatch is lower in more heterozygous F2 hybrids. Heterozygosity values are discrete because a small number of loci underlie adaptation in the plotted simulation run. Simulations are outlined in the Methods. The “mismatch”–heterozygosity relationship is stronger, although less intuitive, in organisms with greater dimensionality (i.e., more traits; see S1 Fig for a case with 10 traits following [18]). The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

In this study, we compare patterns of selection on ancestry heterozygosity between F2 hybrid stickleback raised indoors in aquaria to those from the same cross types raised in the field in experimental ponds (see Table 1 for overview of data sources). We first consider hybridization between sympatric benthic and limnetic populations of threespine stickleback. These populations, which are reproductively isolated species due to strong assortative mating [19] and experience reduced hybrid fitness due to extrinsic selection pressures [20], have evolved independently in at least 5 watersheds in British Columbia, Canada [21,22]. Although reproductively isolated in the wild, the species pairs have no known intrinsic barriers that reduce fitness in the lab [20]. Second, we consider hybridization between allopatric populations of anadromous and solitary freshwater stickleback. As with the benthic–limnetic species pairs, these populations are recently diverged and can readily hybridize. The benthic × limnetic crosses are between specialist populations that are at extreme opposite ends of the range of trophic phenotypes observed among this region’s native stickleback in fresh water. The marine × freshwater cross involves a limnetic-like marine population and a generalist freshwater population (Cranby Lake) that is intermediate between the limnetic and benthic populations [23,24] and is therefore a less divergent cross with respect to trophic characters—however, the populations differ in other traits involved in marine–freshwater divergence [25].

Table 1

Summary of data sources.

Cross type	Design	Study	Population	Method	Generation	Environment	n fish	n loci ± [1 SD]
Ben × lim	Biparental	unpublished	Priest	Microsatellites	F2	Lab	90	22.9 ± 1.1
Ben × lim	Biparental	Conte and colleagues [26]	Priest	SNP array	F2	Pond	412	89.0 ± 0.0
Ben × lim	Biparental	unpublished	Paxton	Microsatellites	F2	Lab	89	97.1 ± 5.8
Ben × lim	8 × F0	Arnegard and colleagues [11]	Paxton	SNP array	F2	Pond	615	62.5 ± 17.9
Ben × lim	Biparental	Conte and colleagues [26]	Paxton	SNP array	F2	Pond	636	62.0 ± 0.0
Ben × lim	See methods	Bay and colleagues [27]	Paxton	SNP array	F2	Pond	302	183.2 ± 19.7
Ben × lim	4 × biparental	Rennison and colleagues [28]	Paxton	GBS (RAD)	F2 and F3	Pond	649	85.1 ± 34.0
Marine × fresh	Biparental	Rogers and colleagues [29]	LCR* × Cranby	Microsatellites	F2	Lab	374	59.2 ± 4.3
Marine × fresh	Biparental	Schluter and colleagues [30]	LCR* × Cranby	SNP array	F2 and F3	Pond	723	120.3 ± 5.6

*Little Campbell River anadromous.

GBS, genotyping by sequencing; RAD, restriction site–associated DNA; SNP, single nucleotide polymorphism.

If ecological selection on trait mismatch is operating in the field but not in the lab, selection for increased ancestry heterozygosity should be specific to the field (or at least stronger than in the lab). If mismatch is deleterious, the fitness landscape is assumed to be saddle or ridge like (see [11]); hybrids with mismatched phenotypes are displaced along the steep sides orthogonal to the axis of parental divergence and have lower fitness than individuals with relatively “matched” trait values (whether parental or somewhat intermediate). Thus, we predicted that we would observe elevated excess ancestry heterozygosity in samples from the field compared to those from the lab.

Results

In support of our prediction, mean individual excess ancestry heterozygosity—the deviation from Hardy–Weinberg expectations based on the relative frequency of alternative ancestry alleles in the genome—was significantly elevated among pond-raised fish compared to aquarium-raised fish. This was the case in both the benthic × limnetic data ( = 0.021 ± 0.0081 [magnitude of excess ancestry heterozygosity ± SE], z = 2.62, P = 0.009) (Fig 2A—left; also see S2 Fig for plots of individual hybrid index and heterozygosity) and the marine × freshwater data (Fig 2A—right; = 0.038 ± 0.0065 [SE], z = 5.86, P < 0.0001). Patterns were similar for all studies in the dataset—each study that contributed data from pond experiments found significant excess ancestry heterozygosity, and each study that contributed data from aquaria found that excess ancestry heterozygosity did not differ from 0 (Fig 2B). The signal of excess ancestry heterozygosity was variable among chromosomes, although the majority had values exceeding 0 (S3 Fig).

Fig 2

Excess ancestry heterozygosity in recombinant threespine stickleback hybrids in the lab (aquarium) and field (pond).

Panel (A) shows the main test of group differences for (left) benthic × limnetic crosses and (right) marine × freshwater crosses. Drawings above the panels were done by K. Chu and show the first listed population on top. These values are extracted from our statistical model using visreg (colored points and violins; [31]) and emmeans (black estimates of means and CIs; [32]). Panel (B) shows raw data (i.e., not from a statistical model) for each data source (Table 1; LCR × CRN is Little Campbell River × Marine) separately, with colors representing lab (red) versus pond (blue) as in (A). Violin overlays show the full distribution of the data, and small colored points show values for individual fish. Large black points are the group means and 95% CIs. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

Discussion

Ecological selection acting on hybrids is a critical determinant of gene flow between diverging lineages [33]. Yet, detecting how divergent alleles interact to mediate hybrid fitness in ecological contexts has proven difficult due to the effect sizes of interacting loci and the massive experiments required to achieve sufficient power [34]. Here, we tested whether a coarse-grained signal of selection against hybrid incompatibilities—elevated excess ancestry heterozygosity [17]—differed between laboratory and field replicates of genetic crosses between the same populations. We found that excess ancestry heterozygosity was elevated in recombinant stickleback hybrids raised in experimental ponds compared to those from similar crosses raised in aquaria. This result is consistent with the hypothesis that certain ecologically mediated hybrid incompatibilities between recently diverged stickleback populations act more strongly in field settings than in the lab. Our finding implies that individual stickleback with a greater mismatch in parental traits are less likely to survive than those with lesser mismatch. Below, we consider whether other processes could plausibly explain this result, discuss the relevance of our findings for speciation, and highlight opportunities for future research.

The lack of excess ancestry heterozygosity in hybrids raised in aquaria is an expected result given what is known about “intrinsic” hybrid incompatibilities in stickleback. Previous studies of benthic × limnetic hybrids have found no evidence for intrinsic inviability in F2 crosses using measures of embryo development and hatching success [22,20,35] or lifetime fitness [20]. A recent review summarizing the literature on reproductive isolation in threespine stickleback [36] reports that “intrinsic” barriers are typically weak to nonexistent. Both marine (i.e., anadromous) × freshwater and benthic × limnetic crosses had no evidence for intrinsic inviability, whereas the authors found evidence for hybrid ecological inviability in both systems [36]. Thus, it would have been surprising if we had recovered any signal of selection on excess ancestry heterozygosity in lab-raised hybrids.

Alternative causes of excess ancestry heterozygosity

We hypothesize that selection against trait mismatch (i.e., ecological hybrid incompatibility) caused the observed patterns of selection on ancestry heterozygosity among surviving individuals in the ponds, but several other mechanisms could possibly underlie such a pattern. These alternative mechanisms involve processes operating at single loci rather than interactions among loci. Ultimately, the data presented here have limited ability to conclusively distinguish between single-locus processes like heterosis and multilocus processes like incompatibilities, but we discuss the strength of evidence for different possible alternative mechanisms below.

Heterosis refers to a case where, for a given locus, the heterozygote has greater fitness than the parental genotypes. This could obviously lead to selection for elevated ancestry heterozygosity. In the benthic × limnetic crosses, environment-specific heterosis is unlikely based on prior knowledge about hybrid fitness in this system. If heterosis were common, then F1 hybrids should have higher fitness than parents. However, F1 and reciprocal backcross hybrids have lower growth and/or survival than both parent taxa in field experiments [20,37,38]. These patterns are opposite to what would be expected as a result of field-specific heterosis, suggesting that it is unlikely to be acting in the benthic × limnetic crosses. Less is known about heterosis in marine × freshwater crosses, although Hagen [39] found evidence for strong ecology-based selection against hybrids in the hybrid zone between the Little Campbell River anadromous population and the resident freshwater-adapted stream population. In sum, the available evidence from observing selection on F1 hybrids in the field versus lab is inconsistent with there being environment-specific heterozygote advantage. The above patterns also suggest that alleles conferring a fitness advantage are not dominant to lower fitness alleles and thus that environment-specific dominance for fitness does not result in pond-specific selection for increased ancestry heterozygosity.

It is conceivable that heterosis could act in the lab via body condition and general health, but only affects mortality in the field and thus would not affect mean excess ancestry heterozygosity in aquaria. Previous studies have found that the growth rate of F1 benthic × limnetic hybrids in the lab matches the additive expectation of parents [20,40], consistent with there being no intrinsic hybrid vigor. In the present study, we find no relationship between body size and excess ancestry heterozygosity in any of the aquarium-raised crosses (S4 Fig; robust within-family regression analyses are not possible in ponds due to the small number of individuals within families; see [11] for discussion of using body size as a component of fitness). Thus, we conclude that there is no evidence supporting a relationship between ancestry heterozygosity and vigor in the laboratory.

We can also use predictions that are specific to the genetics of heterosis versus hybrid incompatibilities to differentiate between them. Specifically, heterosis depends on interactions within a locus, while incompatibilities depend on interactions between loci. If the benefit of ancestry heterozygosity was due to interactions within loci alone, we would expect to see no relationship between genome-wide admixture proportion (i.e., hybrid index) and excess ancestry heterozygosity—this is because all heterozygosity is beneficial regardless of the genetic background. By contrast, if excess heterozygosity was attributable to interactions among loci, as in the hybrid incompatibility model, we would expect diminishing benefits of excess ancestry heterozygosity as genome-wide ancestry proportions became more parent like (i.e., as hybrid index deviates from 0.5). As expected under the incompatibility model, excess ancestry heterozygosity declines as the hybrid index of pond-raised individuals deviates from 0.5 (Spearman ρ = −0.059; P = 0.0006), while there is no relationship in the lab (ρ = −0.012; P = 0.77) (S5 Fig). While this observation does not eliminate the possibility that some of the observed excess heterozygosity is driven by environment-specific heterosis, it is consistent with our hypothesis that excess ancestry heterozygosity largely results from selection against extrinsic incompatibilities.

Inbreeding depression could also generate our observed pattern if it only affected fitness in ponds [41-43]. If inbreeding depression led to selection on ancestry heterozygosity in our study, we would expect inbred individuals to have higher mean excess ancestry heterozygosity than outbred individuals. To test this specific prediction, we used the data from Arnegard and colleagues [11] where due to the study design, some F2 hybrids resulted from mating between full siblings (n = 74), and some individuals were produced by mating between unrelated individuals (n = 541). We find that these 2 groups do not differ in mean excess ancestry heterozygosity (F1,612 = 0.0734, P = 0.787; S6 Fig). Furthermore, the direction of the nonsignificant difference in mean (estimate of outbred − sibling difference = 0.0035 ± 0.0129 [SE]) is the opposite of what is expected under environment-specific inbreeding depression. We note, however, that this analysis does not account for the possibility of any fixed deleterious alleles in the populations [44,45]. Thus, the data suggest that environment-specific inbreeding depression is not causing heterosis in this system.

Finally, an additional analysis indirectly provides evidence of phenotype-based hybrid incompatibilities in the data from Arnegard and colleagues [11]. Specifically, Arnegard and colleagues [11] classified F2 hybrids into 4 groups (“A,” “B,” “L,” and “O”) based on individual niche use. “B,” “O,” and “L” individuals had benthic-like, intermediate, and limnetic-like diets, respectively. “A” individuals, however, had unusual diets, were smaller, and had a greater extent of mismatched trait combinations compared to the other groups. The authors hypothesize that trait mismatch caused these fish to grow more slowly than more “matched” individuals. Our analysis reveals that “A” group individuals have lower excess ancestry heterozygosity than non-“A” individuals (S7 Fig)—as expected if lower excess ancestry heterozygosity correlates with higher trait mismatch. This reanalysis suggests a link between trait mismatch, ancestry heterozygosity, and fitness in stickleback hybrids raised in a natural environment.

Relation to other studies of incompatibilities and ancestry heterozygosity

Our results contribute to a growing understanding of the biology of environment-dependent hybrid incompatibilities. In natural hybrid populations of swordtail fishes, Schumer and colleagues [46] estimated that dozens of incompatibilities separate parent species [47]. The authors also suggested that many of these are likely subject to natural or sexual selection [46]. Previous studies on hybrid stickleback [11,48] have estimated fitness landscapes that are consistent with the hypothesis that mismatched trait combinations are selected against, and our analysis of genetic data supports this hypothesis. We also note that selection against mismatched combinations of traits has the same genetic basis as selection against single phenotypes that express maladaptive transgressive values after hybridization [18,49,50]. Thus, studies focused on only a single trait under stabilizing selection might still find selection against incompatibilities if hybrids have trait values that are below or above the optimum value [51]. More broadly, our results are consistent with predictions generated from theoretical models of speciation and adaptation [17]. Thus, although stickleback is an excellent system in which to test these predictions, the mechanisms underlying our results are likely general.

Our findings also highlight differences from previous analyses of selection on hybrids in different environments. In yeast, selection for low ancestry heterozygosity is common in hybrids when tested in the lab [52,53]. This difference (which does not occur as a result of aneuploidy) between the yeast studies and our study of stickleback might result from the fact the lab media that yeast were raised in are novel environments, and transgressive traits suited to these environments result from excess ancestry homozygosity (see S1 Fig; also see [54] for how novel environments can result in directional selection on ancestry). Stickleback populations in postglacial lakes are specialists on zooplankton or benthic invertebrates when they coexist with a competing fish species [23,55,56] or are generalist populations that make use of both niches [57]. Freshwater stickleback populations in this region span a range of phenotypes along a limnetic–benthic axis, and no specialists along novel trophic axes are known to occur. Fish in the experimental ponds that we considered herein have diets that are largely representative of what they consume in nature [11], and there is no evidence supporting the hypothesis that ponds contain novel adaptive peaks. Qualitative patterns of selection against hybrids driving excess ancestry heterozygosity might therefore depend on the availability and nature of novel ecological niches.

Outlook, caveats, and conclusions

While we identify genetic signatures consistent with the existence of environment-specific hybrid incompatibilities, we cannot begin to identify their specific mechanisms, including when they arise during ontogeny, without connecting phenotype to genotype. Experiments that directly manipulate individual phenotypes, or manipulate interactions between individuals and their environments, are needed to establish such causality. As a field, we should aim to identify the types of traits that typically underlie ecological hybrid incompatibilities. Integrating field studies of hybrid incompatibility with QTL mapping of ecologically important traits [11,58] represents an exciting new frontier for empirical research into the mechanisms of speciation. Moreover, our data only allow us to scratch the surface of how incompatibilities are spread across the genome—it will be valuable for future studies with higher resolution genomic and phenotypic data to investigate this further.

We also do not know the strength of selection against ecological hybrid incompatibilities. Simple simulations (included in archived R scripts) illustrate that the strength of selection necessary to generate 3% excess ancestry heterozygosity in a population of FFV hybrids (similar to that observed in the present study) can vary by orders of magnitude depending on assumptions about the genetic architecture of selection. Moreover, the degree to which mismatched trait combinations are expressed in hybrids, and thus the ability to detect the coarse signal of ecological incompatibilities in an F2 cross, depends considerably on the underlying genetic architecture of adaptive divergence—in particular, the number, effect sizes, and dominance of QTL. Because the expression of maladaptive trait combinations is predicted to increase with the magnitude of divergence between parent populations [18,24,50], we may predict that the strength of selection against ecological incompatibilities will increase with the magnitude of divergence between parents. However, quantifying the ecological basis of incompatibilities and their genetic structure will remain technically challenging.

The evidence presented here is consistent with the hypothesis that extrinsic hybrid incompatibilities are an important mechanism of postzygotic isolation in this system. Our results imply that selection against ecologically mediated hybrid incompatibilities is active from the earliest stages of divergence. Speciation is largely complete when divergent lineages can stably coexist in sympatry, as is the case for the benthic and limnetic stickleback species pairs. Reproductive isolation between them is primarily thought to have arisen incidentally as a by-product of phenotypic divergence [20,38,59,60], with additional selection favoring the reinforcement of premating barriers as a result of low hybrid fitness [20,61]. Thus, our results are consistent with the idea that selection against trait mismatch, or ecological hybrid incompatibilities, is associated with extrinsic postzygotic isolation in this classic system of rapid and recent adaptive radiation.

Materials and methods

Data sources

We used both previously published and unpublished data in our analyses. Summary information about each data source is listed in Table 1. We base our main inference on a comparison of ancestry heterozygosity in hybrids born and raised in aquaria to hybrids from the same cross types born and raised in experimental ponds, which are seminatural ecosystems. See [11,30] for additional information about ponds. Pond-raised crosses capture both “intrinsic” and “extrinsic” incompatibilities, whereas aquarium-raised crosses are expected to capture “intrinsic” incompatibilities that impact hybrid fitness. In most cases, F2 hybrids were produced via mating between full siblings in both the lab and field.

Studies raising fish from the same cross type (benthic × limnetic or marine × freshwater) in the same environment (aquaria or pond) were combined for analysis (see Fig 2). Two studies [28,30] genotyped both F2 and F3 hybrids, which we analyze together because we found that this choice does not affect our conclusions (S8 Fig). Similarly, data were analyzed together for 4 studies of benthic × limnetic hybrids from Paxton Lake raised in experimental ponds because excess ancestry heterozygosity was statistically indistinguishable among them (F = 0.304; P = 0.82); note that this indicates that potential problems with genotyping by sequencing (GBS), such as improperly inferring heterozygotes [62], likely do not affect our conclusions. This grouping of studies was done only to simplify the presentation of results—patterns are highly repeatable across replicates, and analyses showing results for each pond and/or study separately are shown in S9 Fig. Relevant details of each data source are outlined below, but see the original studies for full details including animal use permits.

Studies genotyped fish using either microsatellites, single nucleotide polymorphism (SNP) arrays, or GBS. All 3 lab studies used microsatellites, whereas the pond studies all used SNP arrays or GBS. We examined potential concerns resulting from different genotyping methods and found no evidence that our results are caused by such differences. First, we only consider loci where parents have no alleles in common and thus can accurately polarize ancestry. Second, we only use loci that were heterozygous for ancestry in F1s, so loci with “null” microsatellites or any difficulties in distinguishing alleles would be filtered out (see section on data filtering below). In the largest microsatellite dataset [29], 100% of loci that were different in parents were accurately called as heterozygous across 8 F1s (288 of 288 loci across all 8 F1 fish). SNP genotypes of the same cross type similarly had 100% heterozygosity in F1s [30]. Finally, in 1,000 simulations resampling our dataset to only a single marker per chromosome, 99.4% of estimates of our statistical main effect were in the same direction as detected in the full dataset. In light of the above, we suggest our the differences detected between lab and pond datasets reflect biology rather than methodology. Allele frequencies and heterozygosity are shown for all genotyped loci (within a given dataset) in S10 Fig.

One additional difference between the pond and aquarium studies is that pond F2 hybrids were a result of natural mating among F1 hybrids, whereas aquarium F2s were produced via artificial crosses. We do not anticipate that this will affect our conclusions, however, because we only consider loci that were fixed differences between F0s, and thus are expected to segregate in a 1:2:1 pattern regardless of the process that united eggs and sperm.

Benthic × limnetic crosses

We obtained data from 4 sources for the pond-raised benthic × limnetic hybrids. The data from aquaria are unpublished. Relevant details of each data source are given below. Final sample sizes from each data source are given in Table 1.

Conte and colleagues [26] generated a single F1 family from each of the Priest and Paxton Lake species pairs. Both were founded by a single wild benthic female and a limnetic male that were collected and crossed in 2009. Thirty-five adult F1 Paxton Lake hybrids and 25 adult F1 Priest Lake hybrids were released into separate ponds where they bred naturally to produce F2 hybrids. F2 adults were collected over 1 year later and were genotyped using a SNP array [63]. A total of 246 SNPs were found in the Paxton cross, and 318 were found in the Priest cross.

Arnegard and colleagues [11] conducted a pond experiment with 8 F0 grandparents from Paxton Lake. Two crosses were between limnetic females and benthic males, and 2 crosses were between benthic females and limnetic males. Five F1 males and 5 F1 females from each family were added to a single pond in March 2008 where they bred naturally. Juvenile F2s were collected in October of that same year and genotyped at 408 SNPs using the SNP array.

Bay and colleagues [27] genotyped F2 hybrid females between Paxton Lake benthics and limnetics. Fish are from several crosses and study designs. One used a cross with 4 unique F0s that were used to produce 2 F1 families—one with a limnetic as dam and the other with a benthic as dam. A second had 8 unique F0s, where 2 F1 crosses were in each direction. These 2 crosses used wild fish collected in 2007. A third set of crosses was done in 2009, 1 in each direction, then the 2 F1 families were released into separate ponds. Since the goal of the authors’ study was to examine the genetics of mate choice, a large number of F2 females were genotyped at a small number of microsatellite markers. A subset of F2 females identified in this parentage analysis were genotyped at 494 SNP markers using the SNP array. A total of 302 females were assigned to families with 10 or more full sibs (which was necessary for linkage mapping).

Finally, Rennison and colleagues [28] conducted a study with 4 unique Paxton Lake benthic × limnetic F1 hybrid families that were each split between 2 ponds. One pond in each pair contained a cutthroat trout predator (ancestry heterozygosity did not differ across pond types and data are pooled across all pairs and pond types). Wild fish were caught in 2011, and F1s were released in 2012. Fish bred naturally and juvenile F2s were sampled in September of that same year. F3 hybrids were collected in September 2013. Approximately 50 fish from each pond and hybrid generation were genotyped at over 70,000 loci using restriction site–associated DNA sequencing (GBS), and after filtering and selection of diagnostic loci, we retained 2,243 SNPs.

The Paxton and Priest Lake laboratory cross data are original to this study. Crosses used a single wild-caught benthic female fish and a single wild-caught limnetic male fish as F0 progenitors. Wild fish were crossed in 2003. Sibling mating of F1 hybrids was used to produce a single F2 hybrid family for analysis, and fish were raised in glass aquaria and fed ad libitum. A total of 92 F2s from the Priest Lake cross were genotyped at 84 microsatellite markers, and 86 F2s from the Paxton Lake cross were genotyped at 216 microsatellite markers following [64]. We constructed a combined linkage map between the 2 families using the map integration function in JoinMap 3.0 [65]. The average centimorgan (cM) distance between the markers in the Priest Lake and Paxton Lake maps was 15.6 ± 1.85 and 4.5 ± 0.51 (mean ± 1 SE), respectively. Use of animals was approved by UBC’s Animal Care Committee (A97-0298).

Marine × freshwater crosses

Schluter and colleagues [30] conducted a pond experiment with anadromous (hereafter “marine”) × freshwater hybrids. This study crossed a marine female from the Little Campbell River, British Columbia, with a freshwater male from Cranby Lake, British Columbia. Over 600 juvenile F2 hybrids were introduced into the ponds directly in August 2006. F2s were produced from 6 F1 families—6 unique females were crossed with 4 males (2 males were crossed twice each). F2s overwintered with an estimated over winter survival rate of approximately 86% (from mark–recapture). In spring 2007, surviving F2s bred and were genotyped at 1,294 bi-allelic SNP markers using a SNP array. A total of 500 of their F3 hybrid offspring were collected in October 2007 and were genotyped with the same methodology.

The data for the laboratory marine × freshwater cross were originally published by Rogers and colleagues [29]. The population used a single Little Campbell River female as the F0 dam and a single Cranby Lake male as the F0 sire. Wild adult fish were captured in 2001 to generate a single F1 hybrid family. Four F2 crosses were made from individuals from this F1 family (8 F1 parents total). F2 hybrid fish were genotyped at 96 microsatellite markers.

Marker filtering and estimating excess ancestry heterozygosity

For each dataset, we restricted our analysis to loci where the F0 progenitors of a given F2 family had no alleles in common (e.g., all “BB” in benthic F0s and all “LL” in limnetic F0s) and where all F1 hybrids were heterozygous for ancestry (e.g., all “BL”). GBS data were filtered to include SNPs with > 20× coverage for a given individual. Final sample sizes of fish and markers are given in Table 1. In all cases, the sex chromosome (chromosome 19) was not analyzed.

Because some studies have more individuals than loci, and others have more loci than individuals, we analyze ancestry heterozygosity both in individuals (averaged across loci) and at loci (averaged across individuals). We retained individuals for which at least 20 loci were genotyped, and retained loci for which at least 20 individuals were genotyped. Differences in genotyping success and/or family structure caused the number of genotyped loci to differ among individuals for a given study (Table 1). For simplicity, we focus on the analysis of individuals in the main text.

Deviation from the expected 50:50 ancestry proportions in F2 hybrids—via variance introduced by the assortment of chromosomes into gametes, recombination, and/or directional selection against one ancestry—reduce the expected heterozygosity below 0.5. To account for this, we base our main inference on estimates of excess ancestry heterozygosity. Excess ancestry heterozygosity was calculated as observed ancestry heterozygosity (pAB) minus expected ancestry heterozygosity (2pApB, where pA and pB are the frequencies of both ancestries at the locus or in the individual’s genome). Our conclusions are unchanged, however, if an uncorrected “observed” heterozygosity is used as the response variable (see S11 Fig) or if the expected heterozygosity is adjusted to account for sample size (i.e., multiplying by following [66]; both analyses included in archived R script).

Data analysis

All data processing and analysis was done in R v4.1.1 [67] with packages included in the tidyverse [68]. We compared excess ancestry heterozygosity between aquarium and pond studies using a linear mixed model [69] where excess ancestry heterozygosity was the response, environment (aquarium or pond) cross (B×L or M×F), their 2-way interaction, and lake (Paxton or Priest; for B×L crosses) were fixed factors, and data source was a random effect. We note that our conclusions are unchanged if we analyze excess ancestry heterozygosity of loci averaged across individuals, rather than ofindividuals averaged across their genotyped loci (S12 Fig). We used the emmeans [32] package to evaluate the statistical significance of between-group differences and visreg [31] to visualize models and extract residuals.

Another possible cause of excess ancestry heterozygosity is genotyping error. Simulations of genotyping error—where all errors are conservatively assumed to have resulted in true homozygotes being called as heterozygotes—indicate that error rates in excess of 5% are necessary to cause the pattern we observe (see archived R script). We believe that our error rate is much smaller than 5%, because we encountered no cases in which a locus was falsely called as homozygous in F1 hybrid individuals (which, barring a rare mutation, are known to be heterozygous in all ancestry informative markers). In addition, we expect that genotyping errors would affect lab and pond data similarly.

Simulations underlying conceptual figures

We used simple simulations in Fig 1 of the main text to illustrate the mechanistic relationship between trait mismatch and heterozygosity. Similar results have been noted elsewhere [17,18], but we give our detailed methods herein. We consider the following life history: A single population adapts to a novel environment and then hybridizes with the ancestral population. The phenotype and genotype (hybrid index and heterozygosity) are recorded for F2 hybrids.

We use the framework of Fisher geometric model [70], wherein the phenotype of an organism is a vector of m traits, z = [z1, z2…,z]. For simplicity and ease of visualization, we only consider m = 2 in the main text (see S1 Fig for m = 10). We assume that mutations are sufficiently rare that they individually sweep through an otherwise monomorphic population. We also assume that mutations all occur at unique loci with free recombination among them (i.e., no linkage) [71,72]. Mutations influence the phenotype additively and are vectors of length m where values are drawn from a random normal distribution with a mean of 0 and a standard deviation of α (α = 0.15). The fitness of a given population is calculated as w = exp(−σ||z−o||2), where ||z−o|| is the Euclidean distance between the populations current phenotype (z) and the optimum (o), and sigma is the strength of selection (σ = 10). The original phenotype of the population is z0 = [0,0], and the optimum phenotype is zopt = [1,1]. The selection coefficient s value when they arise, where s = wmutwwt−1. The probability that a given mutation fixes, π, is calculated as π = 1−exp(−2Nsp)/1−exp(−2Ns), where N is the effective population size (N = 1,000), p is the frequency of the mutation in the population (p = 1N), and s is the selection coefficient.

We allow 1,000 mutations to arise in our simulations, which is sufficient for the adapting population to reach the optimum. After adaptation, we simulate hybridization with the ancestral population (which contains no mutations and has a value of 0 for both traits). We generate 500 F2 hybrids, which inherit homozygous ancestral, heterozygous, or homozygous derived ancestry at each locus with probabilities 0.25:0.5:0.25. The genotype of these hybrids determines their phenotype for both traits; these phenotypes are plotted in Fig 1A. We consider 2 orthogonal properties of the phenotype: hybrid index and heterozygosity. Hybrid index is the fraction of alleles an individual inherited from the derived parent, and heterozygosity is the fraction of loci that are heterozygous. We finally calculated the Euclidean phenotypic distance from each hybrid’s phenotype to the line connecting parent phenotypes. This distance is the individual hybrid’s “mismatch” as calculated elsewhere [13].

All data and analysis code have been deposited in the Dryad repository: https://doi.org/10.5061/dryad.h18931zn3 [73].

Simulation model illustrating the incompatibility–heterozygosity relationship.

The model is as in Fig 1 in the main text except there are 10 traits instead of 2 and the optimum of the adapting population is “0” for traits 2 to 9. Plots and model are inspired by Fig 1 in [18]. Both panels depict results from a representative simulation run of adaptive divergence and hybridization between 2 populations. Colored points are individual hybrids, with darker colors indicating higher heterozygosity. The left panel depicts the distribution of 500 F2 hybrid phenotypes where the x-axis depicts the value of the selected trait, and the y-axis depicts the Euclidean distance from the optimum for traits 2–9 (i.e., y = ). Large black points are the 2 parent phenotypes. The right panel depicts the relationship between excess ancestry heterozygosity and the maladaptive distance from the optimum for individual hybrids [13]. Points are slightly jittered horizontally. The plot shows that this maladaptive trait expression is lower in F2s with greater excess ancestry heterozygosity. Heterozygosity values are fairly discrete because a small number of loci underlie adaptation in the plotted simulation run. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

de Finetti ternary diagrams for genotyped individuals.

Each point represents an individual hybrid and shows each individual’s hybrid index (frequency of benthic or marine alleles in its genome) and its mean heterozygosity. Hybrid index and heterozygosity are used because many loci are being considered simultaneously. These graphs are not used for analysis, but rather are shown to allow readers to visualize the structure of the raw data that underlies our analysis. Specifically, the shapes of the distributions of heterozygosity and hybrid index values are similar between environments and crosses—the means are just subtly different. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Mean excess ancestry heterozygosity across chromosomes.

Each point is the average excess ancestry heterozygosity for all loci on a given chromosome (raw values; not residuals from a statistical model). Linkage group XIX contains the sex determining region and is plotted or analyzed. Error bars are 1 SE. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

No relationship between individual mean heterozygosity and growth (standard length) in the aquarium-raised biparental benthic–limnetic F2 hybrids.

Results are residuals from visreg [31]. Each point is an individual F2 hybrid. Standard length is standardized within family (1 family each for Paxton and Priest lakes for B×L and 4 families for M×F. The interaction between lake-of-origin mean heterozygosity was nonsignificant so we plot the main effect across both lakes of origin (Paxton and Priest). Mean heterozygosity was not significantly associated with standard length for either cross ( = −0.05 ± 0.67 [SE], F1,174 = 0.0075, P = 0.93; M = 0.65 ± 0.59 [SE], F1,372 = 1.24, P = 0.26). Analyses considering body depth (either individually or in a combined metric of “overall size”) give the same qualitative result (see archived R script). The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

The benefit of excess heterozygosity declines with deviations from a hybrid index of 0.5 in ponds but not in the lab.

Each point is an individual recombinant hybrid and data are pooled across cross types. Excess ancestry heterozygosity declines as the hybrid index of pond-raised individuals deviates from 0.5 (Spearman ρ = −0.059; P = 0.0006), whereas there is no relationship in the lab (ρ = −0.012; P = 0.77). Bootstrap tests indicate that these 2 correlations are statistically indistinguishable, so we consider this analysis to be interesting and consistent with our hypothesis, but not conclusive. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Comparing selection on ancestry heterozygosity between individuals of different relatedness.

There is no difference in mean excess ancestry heterozygosity between F2 hybrids whose parents were unrelated and those whose parents were full siblings. Data from [11]. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn.

(PDF)

Click here for additional data file.

Fish assigned a priori as “phenotypically mismatched” (“A” group) have lower excess ancestry heterozygosity than nonmismatched fish.

Assignments are from [11] and methods are described therein. Each point is an individual F2 hybrid. This result implies that phenotypically mismatched individuals have lower excess ancestry heterozygosity than nonmismatched individuals. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Excess ancestry heterozygosity does not differ between F2 and F3 hybrids.

The plots show individual excess ancestry heterozygosity from the 2 studies that genotyped both the F2 and F3 generations [28,30]. The means (black dots, ± 95% CI) do not differ between generations in either study (Rennison group difference = 0.014 ±0.0095 [SE], F1,667 = 2.35, P = 0.13; Schluter group difference = 0.0016 ±0.0078 [SE], F1,721 = 0.042, P = 0.84). The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Estimates of mean (± 95% CI) excess ancestry heterozygosity for individuals and loci across “replicates”.

We consider a replicate to be a unique biparental F0 cross for aquarium studies and a unique pond for pond studies. Mean excess ancestry heterozygosity is shown for each such replicate for both individuals (upper) and loci (lower). In each panel, the horizontal line indicates no excess ancestry heterozygosity. Red points are “lab” replicates, and blue points are “pond” replicates. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

de Finetti ternary diagrams for genotyped loci.

Each point represents genotyped locus within a given study (i.e., line in Table 1 in the main text) and shows the frequency of either benthic or marine alleles on the x-axis and its heterozygosity on the y-axis. These graphs are not used for analysis, but rather are shown to allow readers to visualize the structure of the raw data that underlies our analysis. Specifically, the shapes of the distributions of heterozygosity and hybrid index values are similar between environments and crosses—the means are just subtly different. The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Main analysis with observed ancestry heterozygosity rather than excess ancestry heterozygosity as the response.

For full details, see caption of Fig 2 in the main text. Qualitative conclusions of statistical models are identical to those of the main analysis (see archived R script). The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

Test of main hypothesis with loci instead of individuals.

Figure is as in Fig 2 in the main text (and S11 Fig) and qualitative conclusions of statistical models are identical to those of the main analysis (see archived R script). The data and code required to recreate this figure may be found at https://doi.org/10.5061/dryad.h18931zn3.

(PDF)

Click here for additional data file.

15 Jul 2021

Dear Dr Thompson,

Thank you for submitting your manuscript entitled "Genetic evidence for environment-dependent hybrid incompatibilities in threespine stickleback" for consideration as a Short Report by PLOS Biology.

Your manuscript has now been evaluated by the PLOS Biology editorial staff and I am writing to let you know that we would like to send your submission out for external peer review.

Please accept my apologies for the delay incurred while we sought external advice; because of Academic Editor absences, we were unable to secure advice on your paper, and so we will be looking for some enthusiasm from the reviewers to judge whether you paper is ultimately better suited to PLOS Biology or PLOS Genetics.

However, before we can send your manuscript to reviewers, we need you to complete your submission by providing the metadata that is required for full assessment. To this end, please login to Editorial Manager where you will find the paper in the 'Submissions Needing Revisions' folder on your homepage. Please click 'Revise Submission' from the Action Links and complete all additional questions in the submission questionnaire.

Please re-submit your manuscript within two working days, i.e. by .

Login to Editorial Manager here: https://www.editorialmanager.com/pbiology

During resubmission, you will be invited to opt-in to posting your pre-review manuscript as a bioRxiv preprint. Visit http://journals.plos.org/plosbiology/s/preprints for full details. If you consent to posting your current manuscript as a preprint, please upload a single Preprint PDF when you re-submit.

Once your full submission is complete, your paper will undergo a series of checks in preparation for peer review. Once your manuscript has passed all checks it will be sent out for review.

Given the disruptions resulting from the ongoing COVID-19 pandemic, please expect delays in the editorial process. We apologise in advance for any inconvenience caused and will do our best to minimize impact as far as possible.

Feel free to email us at plosbiology@plos.org if you have any queries relating to your submission.

Kind regards,

Roli Roberts

Roland Roberts

Senior Editor

PLOS Biology

rroberts@plos.org

24 Aug 2021

Dear Dr Thompson,

Thank you very much for submitting your manuscript "Genetic evidence for environment-dependent hybrid incompatibilities in threespine stickleback" for consideration as a Short Reports at PLOS Biology. Your manuscript has been evaluated by the PLOS Biology editors, an Academic Editor with relevant expertise, and by three independent reviewers.

In light of the reviews (below), we are pleased to offer you the opportunity to address the comments from the reviewers in a revised version that we anticipate should not take you very long. We will then assess your revised manuscript and your response to the reviewers' comments and we may consult the reviewers again.

IMPORTANT: You'll see that the reviewers are all very positive about your study, but each has a number of requests for improvement. The Academic Editor supplied the following additional advice that you might find useful when prioritising your revisions: "I think reviewer #2 raises some interesting points about interpretation that need to be addressed. In particular his/her suggestion that inbreeding might lead to similar patterns seems valid to me - I think the authors need to address this in the discussion. I would be less concerned about some of the points raised by reviewer #3, such as using a different set of genetic markers, which would represent a significant additional burden of work. I think this goes beyond the scope of this paper, but the implications of possible biases in RAD-seq should be discussed."

We expect to receive your revised manuscript within 1 month.

Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension. At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may end consideration of the manuscript at PLOS Biology.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point by point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Related" file type.

*Resubmission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this resubmission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Roli Roberts

Roland Roberts

Senior Editor

PLOS Biology

rroberts@plos.org

*****************************************************

REVIEWERS' COMMENTS:

Reviewer #1:

I really appreciated this paper addressing the interesting question of coarse effects of incompatibilities across the genome.

I found it concise and well written.

I like the approach and the use of the existing data on sticklebacks to explore the interplay between incompatibilities and the environment.

I only have a few minor questions/comments.

L73: "relatively homozygous" -> relatively more homozygous?

L98 "interactions between individual hybrids": do you mean competition and/or sexual selection? It might be worth adding some detail.

Additionally, biotic interactions can also happen with parents don't you think?

I do not grasp here what you want to convey by separating individual hybrids from other abiotic and biotic factors.

L339 about yeast: I am not a yeast expert at all, but I remember there can be frequent aneuploidies in some crosses. Is it possible it could explain this pattern of low heterozygosity if parts of the genome end up in one copy?

L348: How does ponds and lake environments differ in all those stickleback experiments?

While the difference compared to marine environments for example might not be too high,

I expect the environment of a small pond might still be quite different from the one of a big lake.

In such cases there might be some room for positive selection of some transgressive phenotypes.

Reviewer #2:

In their manuscript titled "Genetic evidence for environment-dependent hybrid incompatibilities in threespine stickleback" Thompson et al. explore the genomic evidence for sets of trait mismatches in hybrids that reduce fitness (what the authors refer to as 'ecological incompatibilities'). This paper presents a creative follow-up to a previously published meta-analysis of trait distributions (and their potential fitness implications) in hybrids (Thompson et al. 2021), and integrates (mostly) previously published genetic/genomic datasets to assess if hybrids reared in semi-wild conditions experience selection against mixed ancestry and subsequent increased heterozygosity more so than hybrids reared in the lab.

Overall I thoroughly enjoyed the manuscript and feel this work provides much needed perspective. I also applaud the authors for a creative use of previously published data. The work makes a novel and, in my opinion, very important contribution to the speciation literature. It is also well written and extremely thorough in analysis. I applaud the authors (again!) for the transparency of the data and methodology. Although my review is overall very positive, I do have some concerns that I think could be addressed with written caveats or minor additional analyses.

1. The authors make the argument that environment-dependent heterosis is unlikely to explain the small increase in heterozygosity in pond populations versus lab reared populations of stickleback hybrids (lines 278-327), although they do acknowledge that "Ultimately, the data presented here have limited ability to conclusively distinguish between single-locus processes like heterosis and multi-locus processes like incompatibilities" (lines 281-284). While I mainly find their arguments against environment-dependent heterosis compelling, there is one scenario that I think is unaccounted for in the current manuscript, which is environment-specific inbreeding depression (given full sibling mating from wild collected grandparents). For example, imagine a scenario in which one or both wild-caught parents harbors recessive deleterious alleles in heterozygous form (perhaps at many loci). Given that these alleles are predicted to be relatively common overall, but rare for a given locus, F1s are also likely to be heterozygous for any given recessive deleterious allele. However, if full sibling F1s are mated, then F2s have a high probability of inheriting both recessive deleterious alleles, and expressing some manifestation of inbreeding depression. In this case, the lines of evidence against environment-specific heterosis are not as compelling: (1) we do not expect that F1s should show heterosis/ increased fitness (recessive deleterious alleles are equally masked in both F1s and in wild parents). (2) Condition might not correlate with overall heterozygosity (just heterozygosity at specific regions of the genome) and/or the inbreeding depression in question might not affect overall condition. (3) I find the argument that excess heterozygosity should decline with further deviations from a hybrid index of 0.5 (and the data presented) fairly weak and unconvincing overall, particularly as these models are indistinguishable when accounting for the vastly different sample size via bootstrapping. Furthermore, if there is any epistatic variation for environment-specific inbreeding depression, one might expect the same results as under a model of 'ecological incompatibilities'. (4) The data from re-analyzed from Arnegard et al. is quite compelling, although I think this particular argument might benefit from subsampling the 'not A' ground to be the same size as the 'A' group, and asking how often these two classes differ from each other (and from 0).

2. This is not an issue per se, but it seems the authors are really wanting to detect whether hybrids have a dearth of pairs of alleles with alternate homozygous ancestry, with an increase in heterozygosity being a natural side effect of this. I wonder if a complimentary analysis would be to calculate some modified hybrid index using only loci that are homozygous within an individual and ask if surviving individuals are skewed to one parent or the other at homozygous regions in the pond versus lab. To avoid any issues caused by directional selection, the authors would need to find skew in both directions (i.e. for both parents) within their dataset(s). I think this type of analysis might also quell some worries about the role of inbreeding depression (and other explanations for environment-specific heterosis) in their dataset, if it's possible to do.

3. I think the method presented by the authors is, overall, of great value to the field, and as previously stated, this work has creatively re-used previously collected data to say something very new, and in my mind, of utmost importance. However, while the authors elucidate throughout the manuscript that ability to detect differences in heterozygosity between field vs lab reared hybrids at these genome-wide coarse scales will depend largely on the polygenic nature of the traits of interest and the scale of additivity/dominance for the alleles controlling those traits. I think the authors should perhaps highlight this caveat under the 'Outlook, caveats, and conclusions' section.

I have many, relatively minor and surely picky points that I think would improve the clarity of the manuscript:

- Lines 43-45: For the severity of 'intrinsic' incompatibilities varying across environments, I also suggest referencing Bundus et al. 2015. Gametic selection, developmental trajectories, and extrinsic heterogeneity in Haldane's Rule. Evolution.

- Line 46: "the number of intrinsic incompatibilities" feels ambiguous to me, because it is unclear whether the authors are referring to the number of interacting loci (as the references would suggest) or the number of unique incompatibilities (which may or may not influence the same hybrid phenotype).

- Line 48: I would also reference Chae et al. 2014. Species-wide genetic incompatibility analysis identifies immune genes as hot spots of deleterious epistasis. Cell.

- Lines 51-55: I might include the example of extreme transgressive trait combinations and developmental mismatches outlined in Coughlan et al. 2021. The genetic architecture and evolution of life history divergence among perennials in the Mimulus guttatus species complex. Proc B.

- Lines 59-60: I find this sentence a bit odd, as the discussion of how polygenic traits are depends, in part, on what type of selection is acting on them. While there is good reason (and evidence) to suggest that Rockman's 2012 piece describes quantitative traits that under stabilizing selection within populations (or even within species), there are also good theoretical reasons to think that the genetic architecture of divergence should be much less complex (see Remington 2015. Alleles versus mutations: Understanding the evolution of genetic architecture requires a molecular perspective on allelic origins. Evolution). Of course, the dearth of genetic studies quantifying the genetic architecture of species divergence and its relation to fitness make these two world views somewhat difficult to assess, but to date there are many studies which have shown that species divergence can be explained by a fairly small handful of loci. This has implications for how effective this differential heterozygosity test will be (see above).

- Lines 210: is this really variation in recombination, or is it variation in independent assortment of chromosomes and/or random chance in which product of meiosis is destined to be the egg (in the case of female meiosis).

- Lines 214-216: Could the authors provide a supplemental figure with the raw, observed heterozygosity data plotted? Similarly, could the authors add the data from the three lab studies and Schluter et al. 2021 to Figure S3? I think this would better allow readers to look at the variance among lab studies and marine studies as well. (The authors could even add each as a separate facet/panel so as to avoid confusion by comparing between different crosses and/or environments).

- Lines 218-230 (Data Analysis section)- do the authors think it would be beneficial to do some sort of bootstrapping of the pond data so that sample sizes are matched between pond and aquaria (which overall have much lower sample sizes) to see how often pond heterozygosity is statistically elevated from 0?

- Throughout the manuscript the authors suggest that mismatched homozygous ancestry is resulting in trait mismatches (which in turn have fitness consequences). Although the some of the authors have previously shown the importance of trait mismatches in hybrid fitness (for example, Thompson et al. 2021), I would be cautious to so definitely say that the reduction in fitness is due explicitly to trait mismatches vs transgressive values of single traits (due to epistatic interactions, for example).

- I just want to reiterate that the link to phenotype, as described in lines 317-327 is absolutely lovely- huge congrats to the authors for thinking of this.

Lines 339-342: the authors might also find the following reference useful: Smukowski Heil et al. 2019. Temperature preference can bias parental genome retention during hybrid evoluition. PLoS genetics.

Reviewer #3:

The authors present a new analysis of numerous previous pond and lab experiments measuring the survival of stickleback hybrid crosses. They make the fascinating observation that heterozygosity is about 3% higher in pond-raised stickleback than lab-raised stickleback in both benthic x limnetic crosses and marine x freshwater crosses. They conclude that this pattern is consistent with extrinsic genetic incompatibilities, i.e. field selection against incompatible interacting homozygous loci from different parental lines which does not occur in the lab environment. They reject the hypothesis that heterosis in field environments better explains the observed pattern, i.e. hybrids with more heterozygous loci are more fit in ponds.

I think this is a fascinating observation that will inspire numerous field biologists to conduct similar studies or re-examine their old data to look for this pattern. It is also open to numerous interpretations - the strongest alternative being heterosis/balancing selection. The authors are careful to discuss this alternative explanation that heterozygous loci are directly selected for in field environments and offer several arguments in favor of genetic incompatibilities. Indeed, their conclusion is careful in stating that the observed heterozygosity "is consistent with" extrinsic genetic incompatibilities. I agree that it is consistent, but I think the authors could be more balanced throughout their presentation - particularly in the title, abstract, and introduction. In my view, they do not present explicit evidence for extrinsic genetic incompatibilities - the pattern could also be explained by simple heterosis in field environments.

The analyses presented are relatively short, although they do explore some simulations. I was also left wondering - how many genetic incompatibilities should we expect in these crosses and are they enough to be predominant in the genome-wide signal of 3% increased heterozygosity? The incredible replication across different crosses, genotyping designs, and studies is impressive and some of these supplemental results should be included in the main manuscript if possible.

I would expect many more incompatibilities in the more divergent marine x freshwater crosses than within-lake benthic x limnetic crosses, particularly given the snowball effect, however excess heterozygosity doesn't appear to increase in more divergent crosses? Why not?

I also think the authors could compare different types of heterozygous sites to better test their hypotheses. Is excess heterozygosity still observed in gene desert regions or only genic/near-gene regions? Do gene regulatory regions and/or coding regions exhibit more excess heterozygosity? The authors would have to postulate additional reasons for different site patterns between heterosis and genetic incompatibilities, but this could still point in useful directions for sorting alternative explanations.

Finally, all the pond hybrids were the result of natural assortative mating events in the field whereas all the lab hybrids were artificial crosses? As I note below, are the authors worried about the effects of increased inter-F1 family assortative mating resulting in increased heterozygosity in ponds? More detail and controls for lab artificial crossing schemes would be useful.

Similarly, models exploring the effect of different genetic markers on heterozygosity would also be useful. RADseq has a well known bias against detection of heterozygous sites due to restriction-site dropout (Arnold et al. 2013).

**************

Line 26: Could the authors qualify this a bit more? I suppose they don't consider tests of postzygotic extrinsic isolation to be tests of hybrid incompatibilities?

Line 32: Space is limited, but it would be good to know why heterosis is an unlikely explanation? The authors could be much more balanced throughout in presenting heterosis as a plausible additional explanation.

Also, I think the authors should clarify in line 30 that this excess heterozygosity could also be a signature of heterosis, not just selection against incompatibilities as they claim, correct?

Line 35: incompatibilities or environment-dependent heterosis?

Line 76: How does this scale with the genetic architecture of the phenotype? What happens for more polygenic traits if opposite homozygous ancestry is not the predominant pattern?

Line 88: The authors are also assuming some disruptive selection on trait-matched hybrids to get the saddle shape? As stated, this suggests a ridge, not a saddle.

Line 97: I'm not convinced this is 'compelling' support for extrinsic incompatibilities because there are alternative explanations for the observed pattern of increased heterozygosity in field samples.

Line 126: This looks thorough, but can the authors also run a glm with the effects of experiment type, genotyping method, and their interaction on heterozygosity levels? It would be good to know if the magnitude of increased heterozygosity observed in the field environments differs across different studies?

Line 136: So, all field experiment hybrids were the result of natural assortative mating whereas all lab 'control' hybrids were random crosses with no assortative mating? Since there are multiple F1 families in each experiment, was there increased inter-family assortative mating?

I agree that only considering fixed differences between parents in the F2 generation probably makes this a moot point, but could assortative mating also potentially explain an increase in the survival of more heterozygous fish in fish enclosures? E.g. increased assortative pairings between F1 families in pond experiments could result in more outbred F2 offspring compared to random lab pairings?

Line 222: An additional model should include the effect of genotyping method (GBS, etc.). Could also include models with effects of sex, presence/absence cutthroat trout, assortative mating/or artificial crosses in the lab to see which model best explains heterozygosity patterns?

Line 295: Is this true - there was zero mortality in the lab before juveniles were released into field enclosures? Heterosis could also act on increased viability of eggs - were there zero egg deaths in lab reared clutches?

Line 305: remove extra 'which'

Line 312: This is a very interesting argument in support of the author's extrinsic genetic incompatibilities argument. Should we really expect the fitness effects of genetic incompatibilities to decline linearly from 0.5? Also, why isn't a quadratic fit used as fitness should decline in either direction away from 0.5 ancestry?

Line 367: Citations would be useful here, too?

Fig. 1 could probably be expanded with at least on additional panel. E.g. what are the effects of selection on trait mismatch and heterozygosity? How does heterozygosity relate to fitness effects of genetic incompatibilities vs. heterosis (see my previous comment)?

Looking at Fig. 2, I wonder if this pattern is mostly/entirely driven by the tails of the distribution in pond experiments? Would ranked analyses produce similar results to parametric results?

Some fish are approaching 40% heterozygous ancestry? What is the predicted vs. observed frequency of these individuals in an F1 intercross (the design of these studies)?

15 Oct 2021

Submitted filename: GEEI_R1_R2R_FINAL.pdf

Click here for additional data file.

1 Nov 2021

Dear Dr Thompson,

Thank you for submitting your revised Short Reports entitled "Genetic data suggest that hybrid incompatibilities in threespine stickleback are environment-dependent" for publication in PLOS Biology. The Academic Editor and I have assessed your revisions and responses to the reviewers.

Based on this assessment, we will probably accept this manuscript for publication, provided you satisfactorily address the following editorial and policy-related requests.

IMPORTANT:

a) We suggest something more specific for the title: perhaps "Genetic signatures suggest that hybrid incompatibilities in threespine stickleback are environment-dependent" or "Ancestry heterozygosity suggests that hybrid incompatibilities in threespine stickleback are environment-dependent" or "Analysis of ancestry heterozygosity suggests that hybrid incompatibilities in threespine stickleback are environment-dependent"

b) Your financial statement currently says “The authors received no specific funding for this work.” Can you confirm that this is indeed correct?

c) Please address my Data Policy requests below; specifically, please supply numerical values and/or the code required to recreate Figs 1AB, 2AB, S1, S2, S3, S4, S5, S6, S7, S8 S9, S10, S11, S12, and cite the location of the data clearly in each relevant Fig legend.

As you address these items, please take this last chance to review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the cover letter that accompanies your revised manuscript.

We expect to receive your revised manuscript within two weeks.

To submit your revision, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' to find your submission record. Your revised submission must include the following:

-  a cover letter that should detail your responses to any editorial requests, if applicable, and whether changes have been made to the reference list

-  a Response to Reviewers file that provides a detailed response to the reviewers' comments (if applicable)

-  a track-changes file indicating any changes that you have made to the manuscript.

NOTE: If Supporting Information files are included with your article, note that these are not copyedited and will be published as they are submitted. Please ensure that these files are legible and of high quality (at least 300 dpi) in an easily accessible file format. For this reason, please be aware that any references listed in an SI file will not be indexed. For more information, see our Supporting Information guidelines:

https://journals.plos.org/plosbiology/s/supporting-information

*Published Peer Review History*

Please note that you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*Early Version*

Please note that an uncorrected proof of your manuscript will be published online ahead of the final version, unless you opted out when submitting your manuscript. If, for any reason, you do not want an earlier version of your manuscript published online, uncheck the box. Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us as soon as possible if you or your institution is planning to press release the article.

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Please do not hesitate to contact me should you have any questions.

Sincerely,

Roli Roberts

Roland G Roberts, PhD,

Senior Editor,

rroberts@plos.org,

PLOS Biology

------------------------------------------------------------------------

DATA POLICY:

You may be aware of the PLOS Data Policy, which requires that all data be made available without restriction: http://journals.plos.org/plosbiology/s/data-availability. For more information, please also see this editorial: http://dx.doi.org/10.1371/journal.pbio.1001797

We note that the raw data are provided in your supplementary data files. However, we also need the individual numerical values that underlie the figures to be made available in one of the following forms:

1) Supplementary files (e.g., excel). Please ensure that all data files are uploaded as 'Supporting Information' and are invariably referred to (in the manuscript, figure legends, and the Description field when uploading your files) using the following format verbatim: S1 Data, S2 Data, etc. Multiple panels of a single or even several figures can be included as multiple sheets in one excel file that is saved using exactly the following convention: S1_Data.xlsx (using an underscore).

2) Deposition in a publicly available repository. Please also provide the accession code or a reviewer link so that we may view your data before publication.

Regardless of the method selected, please ensure that you provide the individual numerical values that underlie the summary data displayed in the following figure panels as they are essential for readers to assess your analysis and to reproduce it: Figs 1AB, 2AB, S1, S2, S3, S4, S5, S6, S7, S8 S9, S10, S11, S12. NOTE: the numerical data provided should include all replicates AND the way in which the plotted mean and errors were derived (it should not present only the mean/average values).

IMPORTANT: Please also ensure that figure legends in your manuscript include information on where the underlying data can be found, and ensure your supplemental data file/s has a legend.

Please ensure that your Data Statement in the submission system accurately describes where your data can be found.

------------------------------------------------------------------------

DATA NOT SHOWN?

- Please note that per journal policy, we do not allow the mention of "data not shown", "personal communication", "manuscript in preparation" or other references to data that is not publicly available or contained within this manuscript. Please either remove mention of these data or provide figures presenting the results and the data underlying the figure(s).

------------------------------------------------------------------------

4 Nov 2021

Dear Ken,

On behalf of my colleagues and the Academic Editor, Chris Jiggins, I'm pleased to say that we can in principle accept your Short Reports "Analysis of ancestry heterozygosity suggests that hybrid incompatibilities in threespine stickleback are environment-dependent" for publication in PLOS Biology, provided you address any remaining formatting and reporting issues. These will be detailed in an email that will follow this letter and that you will usually receive within 2-3 business days, during which time no action is required from you. Please note that we will not be able to formally accept your manuscript and schedule it for publication until you have any requested changes.

IMPORTANT:

a) Many thanks for supplying the Dryad link. This enabled me to check it and sign off on the MS editorially. However, you will need to include the Dryad URL/DOI in each relevant main and supplementary figure legend (e.g. "The data and code required to recreate this Figure may be found at https://datadryad.org...."). I've left a note with my Production colleagues to contact you about this.

b) Those fish pics in Fig 2 are fantastic!

Please take a minute to log into Editorial Manager at http://www.editorialmanager.com/pbiology/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production process.

PRESS: We frequently collaborate with press offices. If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximise its impact. If the press office is planning to promote your findings, we would be grateful if they could coordinate with biologypress@plos.org. If you have not yet opted out of the early version process, we ask that you notify us immediately of any press plans so that we may do so on your behalf.

We also ask that you take this opportunity to read our Embargo Policy regarding the discussion, promotion and media coverage of work that is yet to be published by PLOS. As your manuscript is not yet published, it is bound by the conditions of our Embargo Policy. Please be aware that this policy is in place both to ensure that any press coverage of your article is fully substantiated and to provide a direct link between such coverage and the published work. For full details of our Embargo Policy, please visit http://www.plos.org/about/media-inquiries/embargo-policy/.

Thank you again for choosing PLOS Biology for publication and supporting Open Access publishing. We look forward to publishing your study.

Best wishes,

Roli

Roland G Roberts, PhD

Senior Editor

PLOS Biology

rroberts@plos.org