Genomic landscape of reproductive isolation in Lucania killifish: The role of sex loci and salinity.

Adaptation to different environments can directly and indirectly generate reproductive isolation between species. Bluefin killifish (Lucania goodei) and rainwater killifish (L. parva) are sister species that have diverged across a salinity gradient and are reproductively isolated by habitat, behavioural, extrinsic and intrinsic post-zygotic isolation. We asked if salinity adaptation contributes indirectly to other forms of reproductive isolation via linked selection and hypothesized that low recombination regions, such as sex chromosomes or chromosomal rearrangements, might facilitate this process. We conducted QTL mapping in backcrosses between L. parva and L. goodei to explore the genetic architecture of salinity tolerance, behavioural isolation and intrinsic isolation. We mapped traits relative to a chromosome that has undergone a centric fusion in L. parva (relative to L. goodei). We found that the sex locus appears to be male determining (XX-XY), was located on the fused chromosome and was implicated in intrinsic isolation. QTL associated with salinity tolerance were spread across the genome and did not overly co-localize with regions associated with behavioural or intrinsic isolation. This preliminary analysis of the genetic architecture of reproductive isolation between Lucania species does not support the hypothesis that divergent natural selection for salinity tolerance led to behavioural and intrinsic isolation as a by-product. Combined with previous studies in this system, our work suggests that adaptation as a function of salinity contributes to habitat isolation and that reinforcement may have contributed to the evolution of behavioural isolation instead, possibly facilitated by linkage between behavioural isolation and intrinsic isolation loci on the fused chromosome.

INTRODUCTION

Local adaptation often plays a critical role in speciation when the force of natural selection is channelled into the generation of reproductive isolation (i.e. ecological speciation, Nosil, 2012; Schluter, 2009). This occurs in two ways: (1) natural selection may directly generate reproductive isolating barriers such as habitat isolation, immigrant inviability, or extrinsic post‐zygotic isolation; or (2) the force of natural selection may contribute to reproductive isolation indirectly (i.e. as a by‐product) (Cruz, Carballo, Conde‐Padin, & Rolan‐Alvarez, 2004; Kay, 2006; Lowry et al., 2008; Moser et al., 2016; Nosil, 2012; Nosil et al., 2009; Nosil et al., 2005; Ramsey et al., 2003; Rundle & Nosil, 2005; Schemske, 2010; Schluter, 2001). The strength of divergent natural selection can be indirectly funnelled into the creation of other forms of reproductive isolation via physical linkage between adaptive loci and loci that contribute to other barriers (but are not under natural selection themselves) or through pleiotropy where a given locus is under natural selection via its effect on one trait, with indirect effects on traits influencing reproductive isolation. Environmental factors requiring adaptation in different traits at multiple life stages may play a larger role in generating both direct and indirect reproductive isolation as they affect many loci in the genome and provide stronger barriers to migration, increasing divergence between the species (Feder & Nosil, 2010; Nosil et al., 2009). Studies of genomic architecture are particularly useful in determining contributions (direct or linked) of natural selection to reproductive isolation by testing for signatures of natural selection on barrier loci (Presgraves et al., 2003) or genomic signatures of selection on hybrids along ecological clines (Grahame et al., 2006; Schaefer et al., 2016; Sung et al., 2018). Genome‐wide association studies often find co‐localization of loci contributing to ecological traits and loci contributing to reproductive isolation (Bay et al., 2017; Renaut et al., 2011). However, linkage disequilibrium between loci coding for ecological traits and those contributing to reproductive isolation is consistent with ecological speciation but can also be caused by spatial subdivision among species (Beirne, Welch, Loire, Bonhomme, & David, 2011). Genetic mapping studies can help differentiate between these two possibilities. Much of this prior work has focused on terrestrial or aquatic freshwater systems. The indirect effects of natural selection on reproductive isolation in marine and coastal species remain less well studied.

An important environmental factor driving local adaptation in coastal aquatic environments is salinity. The transition from fresh water to salt water is marked by a rapid shift in communities in nearly all groups of aquatic organisms, indicating the strength of this boundary (Godfrey & Wooten, 1979, 1981; Gunter, 1945, 1950a, 1961). Multiple studies have suggested salinity has played a significant role in speciation in teleost fish (Betancur et al., 2015; Hrbek & Meyer, 2003; Huyse, Van Houdt, & Volckaer, 2004; Lee & Bell, 1999; Whitehead, 2010). In particular, salinity often contributes to habitat isolation among closely related species (Torres‐Dowdall et al., 2013). Environmental salinity requires complex physiological adaptation because, in high salinity environments, organisms are subject to ion influxes and loss of water from tissues. Conversely, in low salinity environments, fluxes of water into tissues and loss of ions to the environment occurs (Evans, 2008; Evans et al., 2005). This complex adaptation causes physiological divergence within species in many tissues and life stages (Seehausen & Wagner, 2014; Taylor, 1999; Whitehead et al., 2011). Previous work suggests the genomic basis of this important trait may be dispersed across the genome on multiple chromosomes or linkage groups (Berg et al., 2015; Brennan et al., 2018). Whereas it is clear that salinity tolerance directly contributes to habitat isolation, extrinsic isolation, and immigrant inviability, the contribution of salinity to other reproductive isolating barriers via linked selection remains unknown.

Linked selection can only generate reproductive isolation if traits that contribute to isolating barriers are coupled to traits under divergent natural selection (Butlin & Smadja, 2018). In general, the coupling of different isolating barriers is necessary for the completion of speciation (i.e. the full cessation of gene flow). However, coupling is particularly difficult when speciation occurs with gene flow as recombination directly opposes divergence by homogenizing allelic combinations, thus breaking down linkage disequilibrium between reproductive barriers (Butlin, 2005; Felsenstein, 1981). Chromosomal rearrangements have been hypothesized to facilitate speciation with gene flow because recombination is reduced between homologous chromosomes with different arrangements. Both theoretical and empirical work has shown that the genes underlying reproductive isolation are less likely to become decoupled and homogenized via gene flow when they co‐occur in areas of low recombination (Charlesworth & Barton, 2018; Faria & Navarro, 2010; Hoffmann & Rieseberg, 2008; Kirkpatrick & Barton, 2006; Lowry & Willis, 2010; Noor et al., 2001; Wellenreuther & Bernatchez, 2018; Wellenreuther et al., 2019). Thus, genomic features such as rearrangements may allow the environment to play a larger role in speciation by facilitating tight physical linkage between loci under divergent natural selection and loci conferring other reproductive isolating barriers. Rearrangements may also contribute to speciation by coupling multiple reproductive isolation barriers unrelated to the environment (Butlin & Smadja, 2018). If a chromosomal rearrangement has facilitated divergence, then the expectation would be that ecologically important traits and reproductive isolating barriers should map to the rearranged region.

In this study, we examine the contribution of divergent natural selection for salinity tolerance to the evolution of behavioural and intrinsic post‐zygotic isolation between Lucania goodei and L. parva and ask whether a chromosomal rearrangement has facilitated the speciation process. Lucania goodei and L. parva are closely related species (Duggins et al., 1983; Whitehead, 2010) that differ in salinity tolerance. Lucania goodei is found primarily in freshwater sites (restricted mainly to Florida and southern Georgia), whereas L. parva can be found in fresh, brackish, and marine habitats as far west as central Mexico and as far north as Massachusetts (Lee, 1980). Within Florida, the two species ranges overlap by ~ 12%–19%. The species co‐occur in either brackish pools along the coast or freshwater portions of rivers that are tidally influenced (Fuller & Noa, 2008).

Differential adaptation to salinity between the two species is present at multiple life stages (Dunson & Travis, 1991; Fuller, 2008; Fuller et al., 2007). In the early life‐history stage, L. parva has high survival across all salinities (0.2–30 ppt). In contrast, L. goodei has high survival in fresh and slightly brackish salinities (0–10 ppt) but has low survival at higher salinities (20–30 ppt) (Fuller et al., 2007). These differences in salinity tolerance are reflected in genetic differences in coding sequence and expression differences between the two species (Berdan & Fuller, 2012b; Kozak et al., 2014).

Hybrids between L. parva and L. goodei can be found in nature at low levels (Hubbs, Walker, & John, 1943) in the portions of the range where the two species co‐occur. Several reproductive isolating barriers other than habitat isolation exist. There is reduced viability of hybrid offspring at high salinities and reduced overall fitness of F2 offspring and backcrosses to L. goodei (Fuller, 2008; Fuller et al., 2007). Previous work on Lucania indicates that a large genetic incompatibility is segregating between the two species that results in some hybrid males having drastically reduced fitness (Fuller, 2008). Offspring from hybrid F1 males (from crosses between L. parva females and L. goodei males) have reduced viability and nearly half die during the first few days of development. Besides these post‐zygotic isolating barriers, assortative mating due to male and female preferences causes behavioural isolation between the two species (Berdan & Fuller, 2012a; Fuller et al., 2007; Kozak et al., 2015; St. John & Fuller, 2019). A pattern consistent with reinforcement is present where this behavioural isolation is heightened in areas where the two species co‐occur (Gregorio et al., 2012; Kozak et al., 2015; St. John & Fuller, 2019).

In addition to phenotypic divergence, previous work has established that a chromosomal fusion has occurred in the lineage leading to L. parva. Karyotypes and genetic mapping indicate that two acentric chromosomes have been fused into a single metacentric chromosome in L. parva (a centric fusion), leading to 1N = 23 in L. parva compared to 1N = 24 in L. goodei (Berdan et al., 2014; Uyeno & Miller, 1971). The sex‐determining region is currently unmapped in Lucania, but karyotypes of both species do not differ between males and females, suggesting that the sex chromosomes are homomorphic (Berdan et al., 2014; Uyeno & Miller, 1971). Homomorphic sex chromosomes are common among many fish and reptile species (Bachtrog, 2013; Bachtrog et al., 2014; Mank, Promislow, & Avis, 2005). Chromosomes containing the sex‐determining region often have important effects in speciation because they can play a disproportionate role in intrinsic isolation and also tend to diverge faster than autosomes (Bachtrog et al., 2011; Coyne & Orr, 2004; Meisel & Connallon, 2013; Presgraves, 2018; Saether et al., 2007). Therefore, we sought to determine the nature of sex determination in this system and to map loci physically linked to the sex‐determining region.

We used crosses between L. goodei and L. parva to map loci linked to the sex‐determining region, early life stage survival at different salinities (hereafter referred to as ‘salinity tolerance’), behavioural isolation (male attractiveness for each species), and intrinsic post‐zygotic isolation (reduced male fertilization success and offspring viability). We wanted to determine the distribution of these traits in the genome and whether behavioural and post‐zygotic isolation loci are physically linked to salinity loci, the chromosomal fusion, or the sex‐determining region. To do this, we created a series of backcrossed pedigrees, phenotyped the backcrossed offspring for sex, salinity tolerance, male attractiveness, male fertilization success and offspring viability, and genotyped the offspring to conduct QTL mapping. We predicted that if salinity adaptation contributes to reproductive isolation indirectly via physical linkage, then loci that contribute to survival at high salinities should map to similar genomic regions as loci for other reproductive isolating traits. We predicted that if genomic features have facilitated coupling of adaptive loci and other reproductive isolating loci, then these traits should co‐localize to the chromosomal fusion and/or the sex chromosome. If the fusion by itself contributes to the coupling of multiple reproductive isolating barriers (independent of salinity), then traits for other reproductive barriers, but not salinity tolerance, should map there.

METHODS

Field collections, animal husbandry, and crosses

Our goal was to QTL map sex determination, salinity tolerance and multiple components of reproductive isolation. To do this, we created three separate crosses among Lucania populations and genotyped the offspring (Crosses 1–3; see Table 1). Cross 1 was the main mapping backcross to L. goodei in which we genotyped all traits (Table 1). We chose to examine reproductive isolation and salinity tolerance from backcrosses into L. goodei because previous work had shown that significant genetic incompatibilities and reductions in survival in saltwater were present in these backcrosses. In contrast, backcrosses to L. parva show limited evidence for reduced survival and genetic incompatibilities (Fuller, 2008; Fuller et al., 2007). Hence, backcrosses into L. goodei were expected to produce offspring in which salinity tolerance and reproductive isolating barriers were segregating. For mapping of sex‐linked loci, we used two additional mapping crosses (Crosses 2 and 3; see Table S1).

Table 1

Crosses used. For sample sizes for each assay see Figure S2, S3, and S4. For location see Figure 1

Cross	Population/Species 1	Population/Species 2	Backcross direction	Phenotypes
1	Delks L. goodei	Delks L. parva	L. goodei	Male Fertilization Success, Male Mating Success, Offspring Viability, Sex, Salinity Tolerance
2	Blue Springs L. goodei	Indian R. L. parva	Both	Sex
3	Pecos R. L. parva	Indian R. L. parva	Both	Sex

Figure 1

Location of the sex determination locus. LOD scores across theL. parvalinkage map for cross 2 interspecies backcrosses to L. parva and L. goodei (a, b) and cross 3L. parva interpopulation cross (c, d). Red dashed line indicates LOD score = 3 (approximately equivalent top < .05). Inset shows LOD score across the fused chromosome LG1 (b, d). Black dashed line separates the putative ancestral 1A and 1B portions of the chromosome (1B starts at 30cM)

Individuals were collected using dip nets and seines between 2009 and 2011 (map in Figure S1) and housed as described in supplemental methods. We crossed fish by pairing a male and a female in a 38‐litre (10 gallon) aquarium. For each cross, F1 families were generated in both directions from parental populations and pairs were checked for eggs every 2–3 days (see supplemental methods). Fry were raised to adulthood and then backcrossed to one or both parental populations (Table 1). All F1 and backcrosses were performed in fresh water. In all experiments, our freshwater source was dechlorinated city water treated with Start Right (Jungle Laboratories, Cibolo, Texas).

For the main mapping cross 1, parental adult L. goodei and L. parva were collected from a population where they co‐occur: The Boat Ramp at Delk's Bluff near Ocala (Marion County, Florida). Families of all four possible types of backcrosses to L. goodei were made (Figure S2). A randomly selected portion of the offspring was raised to one month of age in either fresh or salt water, euthanized and later genotyped. Genotypes were compared between fresh and salt water in order to map the location of salinity tolerance (see below). The remaining freshwater offspring were raised to adulthood to assay male fertilization success and offspring viability (i.e. survival to hatching for offspring of the backcrossed males), behavioural isolation and sex determination (Figure S2).

For mapping crosses 2 and 3, F1s were created in both directions, raised to adulthood, and used to create all possible backcrosses in both directions. These offspring were genotyped to map the location of the sex‐determining region. For cross 2, the parental populations were L. goodei from Blue Springs along the Santa Fe River, Florida, and L. parva from Indian River along Merritt Island, Florida. For cross 3, the parental populations were L. parva from Indian River, Florida and L. parva from the Pecos River, Texas (Figure S1).

Predictions for loci linked to the sex‐determining locus

Our goal was to map loci linked to the sex‐determining region. We assumed that loci that are closely linked to the sex‐determining region would not segregate in Mendelian ratios. Instead, the SNPs should exhibit classic signatures of a given sex determination system (Staelen et al., 2008; Star et al., 2016). As sexually dimorphic loci were not present on our linkage map, we scanned for linked loci that showed the expected species/population of origin patterns (see Figure S3–S4). In the case of an XX‐XY system (male determining locus), F1r (L. parva ♀ × L. goodei ♂) males should pass on an L. goodei allele to male offspring and an L. parva allele to female offspring at loci linked to the sex‐determining region. When backcrossed to L. goodei, we expected the resulting female offspring to be heterozygous and the male offspring to be homozygous for L. goodei alleles for loci linked to the sex‐determining region (Figure S3). Similarly, F1 (L. goodei ♀ × L. parva ♂) males should pass on an L. parva allele to male offspring and an L. goodei allele to female offspring. When backcrossed to L. goodei, we expected the resulting male offspring to be heterozygous and the female offspring to be homozygous for L. goodei alleles. Predictions for a ZZ‐ZW system are in the supplemental material (methods; Figure S3).

We assigned the offspring from diagnostic crosses with a sex genotype consistent with the predictions from a given sex determination system and mapped this as a trait. Adult males and females are visually identifiable by sexually dimorphic coloured anal fins (males have these, females do not; Fuller et al., 2007). Hence, we can assign them genotypes predicted from a hypothesized sex determination system, the identity of their parents and their sex (Figure S3). Our results were consistent with an XX‐XY sex determination system (see results below), which is in keeping with the sex determination system of other members of Fundulidae (Chen & Ruddle, 1970). However, our sample size was low for offspring from diagnostic crosses for the XX‐XY system (i.e. backcrosses involving F1r males).

To map the sex‐determining region more finely, we used backcrossed offspring from two additional crosses (Cross 2: between L. goodei and L. parva; Cross 3: between L. parva populations; Table 1; Figures S3–S4; Supplemental Methods). To map the sex‐determining region, we used the pure species maps from Berdan et al. (2014) (L. goodei map = 915 SNP markers, L. parva map = 766 SNP markers).

Salinity tolerance assays

For the salinity tolerance assay, we used the main mapping backcross to L. goodei (Cross 1). We obtained clutches of backcrossed eggs from multiple families and divided the eggs between fresh water (0.2 ppt salinity, 32 families, 1,152 eggs total) and salt water (15 ppt salinity, 23 families, 521 eggs total). Eggs were placed in groups in 177‐mL tubs treated with methylene blue (anti‐fungal agent) in fresh water (~0.2 ppt) or salt water (15 ppt; see supplemental methods). We chose 15 ppt because it creates a significant osmoregulatory challenge without inducing excessive mortality. Fifteen ppt is greater than the isosmotic point (~10 ppt). Previous studies show that survival to fry to the eating stage is ~ 13% for L. goodei offspring raised in salinities 20–30 ppt compared to ~ 87% at 0.2 ppt (Fuller, 2008). We recorded the number of eggs that hatched and the number of fry that survived to one month of age. At one month of age, a subset of offspring was euthanized and preserved for genotyping (freshwater: 61 offspring from 18 families; saltwater: 84 offspring from 17 families). Additional offspring from the freshwater treatment were raised to adulthood to obtain adults for measurements of behavioural isolation and male fertilization success.

Male behavioural isolation

We measured behavioural isolation between the backcrossed males (cross 1) and both L. goodei and L. parva females. We used a no‐choice mating assay, which has been used successfully in previous studies of behavioural isolation in Lucania (Berdan & Fuller, 2012b; Fuller et al., 2007; Kozak et al., 2012; St. John & Fuller, 2019). Adult backcross males were placed in a 38‐litre (10 gallon) aquarium with a stimulus female: either a female L. goodei or a female L. parva. We checked pairs every 2nd day for eggs over 20 days, recording the day on which eggs were first found and the total number of eggs produced. We obtained stimulus females from two populations – Delk's Bluff (6 trials with L. goodei females and 6 trials with L. parva females) and the Wakulla River (23 trials with L. goodei females and 27 L. parva females). Thirty‐three males were tested in total; 29 males were tested with both L. goodei and L. parva females and another 4 were tested only with L. parva females. We randomized the order in which males were paired with L. goodei and L. parva females.

From these data, we recorded whether or not males mated as a binary variable (latency < 20 days = 1; latency ≥ 20 days = 0) and daily egg production rate, which served as indices of male attractiveness to and willingness to mate with stimulus females. We did not check fertilization status or survival of the eggs for this measure. St. John and Fuller (2019) found egg laying to be more robust than time measures of preference in Lucania. These fish are external fertilizers, so there is little opportunity for male seminal fluids to alter female mating rates. Furthermore, females lay ~ 1–2 eggs per spawning event (Breder & Rosen, 1966), so the number of eggs is approximately equal to the number of spawns. The disadvantage of no‐choice mating assays is that it is difficult to determine precisely what attributes of males and females are being measured. In the Lucania system, both males and females have mating preferences (Kozak et al., 2015; St. John & Fuller, 2019). Hence, the assays may be detecting loci that affect both hybrid male attractiveness and hybrid male preference of females.

Reduced male reproductive success

We assayed both the fertilization success (proportion of eggs fertilized) and the offspring viability (proportion of viable eggs that hatched) from the backcross males mated to females in the mating assays (N = 23 males). We checked all collected eggs under a light microscope to assess fertilization. We considered eggs that were already dead upon collection to be unfertilized. We saved the fertilized eggs and measured their survival until hatching as our metric of ‘offspring viability’. After mating and fertilization trials, males were subsequently euthanized with MS‐222 and stored in ethanol at −20°C.

Genotyping

DNA isolation protocols are described in the supplemental methods. For cross 1, DNA was extracted from 173 offspring from the salinity tolerance assay (61 freshwater, 84 saltwater), 33 males from the male behavioural isolation and intrinsic isolation assays, and 27 females. For cross 2 and 3, DNA from 50 adults and 36 adults were isolated, respectively (Figure S4c). Samples were diluted to a concentration of 75 ng/μl prior to genotyping.

All genotyping was done using a single custom designed Illumina Infinium Beadchip as described in Berdan et al. (2014). The Bead Chip contained probes for genotyping of 4,545 SNPs: 1,679 putatively species‐specific loci (used in crosses 1 and 2), 1,369 segregating among L. parva populations (used in cross 3) and 1,497 segregating among L. goodei populations (not used in this study).

All DNA samples were spread across twelve 96‐well plates and genotyped on the Illumina Infinium Bead Chip, scanned using the iScan System (Illumina) at the Keck Center for Comparative and Functional Genomics at the University of Illinois in 2011–2013. There was no indication that backcross type (Figure S2) influenced the general segregation of markers using a quasi‐binomial glm model with the proportion of homozygous genotype calls as a dependent variable and a likelihood ratio test for difference between models with and without cross type as an explanatory variable ( deviance = 0.100, p = 0.57).

Backcross linkage map

L. goodei and L. parva linkage maps have been previously constructed and genome size in L. goodei was previously estimated to be 1.3 Gb (Berdan et al., 2014). We created a hybrid map because (a) map length differs between L. goodei and L. parva (605 cM and 392 cM respectively) and (b) we did not know whether recombination in hybrids differs from the parental species. We constructed a backcross linkage map from the backcrossed offspring from cross 1 (N = 189 individuals). We constructed the map in Joinmap 4.0 (Van Ooijen, 2006) following methods described in Berdan et al. (2014). We considered all hybrids to come from a single family, and we only constructed the map using species‐specific SNPs markers that were in Hardy‐Weinberg equilibrium (222 markers total) and that showed the correct inheritance pattern (heterozygous in all F1 parents and homozygous in all L. goodei parents).

QTL mapping—sex determination

All QTL mapping, other loci association tests, and visualization were performed in R v.3.5 (R Core Team, 2018), using the rQTL (Broman & Sen, 2009) and ggplot (Wickham, 2016) packages. We tested whether our sex‐specific genotypes were consistent with either an XX‐XY or a ZZ‐ZW sex determination system. To do this, we used the predicted genotypes for sex‐linked loci to map the sex‐determining region on our linkage maps (Figures S3–S4). We used scanone with a binary model to calculate LOD scores and determined significance and the 95% Bayesian credible interval using 5,000 permutations. For cross 1, we used 60 genotyped individuals and the new backcross map. For cross 2, we used 50 genotyped backcross offspring and species‐specific SNPs which had been previously mapped on the L. parva and L. goodei linkage maps (Berdan et al., 2014) (353 markers L. goodei map, 355 L. parva map). For Cross 3, we used 36 genotyped backcrossed individuals (Figure S4), filtered SNP data, identified alleles that were fixed between parental population (Kozak et al., 2014) and used the 821 markers out of these that had a position on the L. parva map.

Salinity tolerance genotype testing

We sought to determine the location of loci associated with salinity tolerance. To do this, we compared the frequency of the different genotypes across the genome among offspring that survived in fresh and salt water. We considered only offspring that were raised to 1 month of age and then euthanized so as to directly compare the same window of mortality between fresh and saltwater (inclusion of adults raised in fresh water does not alter the results). Along a given linkage group, some markers were in complete linkage and were removed (28 of 172 markers removed). Previous work indicates that juveniles of both L. goodei and L. parva survive well in fresh water. We therefore used the frequency of the SNP genotypes among the 61 freshwater offspring as the expected frequency and asked whether the frequencies in saltwater differed from expected proportions using a chi‐square test (df = 1). We corrected for multiple testing by using the Benjamini and Hochberg (1995) method as implemented in R with ‘p.adjust’.

Mapping gametic disequilibrium—interactions among loci

Genetic incompatibilities among loci are predicted to generate distortions in genotype frequencies in surviving backcrossed individuals. Incompatible allele combinations will be reduced or absent among the survivors. The goal here was to determine whether backcrossed offspring differed in their probability of survival due to interactions among genotypes located on different linkage groups. To do this, we tested for nonrandom patterns of genotypes, using a chi‐square analysis. We only included backcrossed individuals that had been raised in fresh water (N = 121 individuals, including offspring and adults) to avoid the distorting effects of differential survival in salt water. We considered both offspring that were raised until one month of age (and then killed) as well as offspring that were raised to adulthood. Analyses that excluded adult backcrossed offspring (n = 61) had little effect on the qualitative results, but reduced power. Along a given linkage group, some of the markers were in complete linkage, so we used one representative marker from each linkage block. We also only considered patterns among loci located on different linkage groups. We performed a total of 10,675 tests. For each test, we calculated chi‐square statistic (df = 1), the associated p‐value and the frequencies of the four combinations of genotypes (homozygous at both locus 1 and 2, heterozygous at both locus 1 and 2, homozygous at locus 1/heterozygous at locus 2 and vice versa). Again, we used the false discovery rate correction in ‘p.adjust’.

QTL mapping—phenotypes of backcrossed males

We sought to determine the location of QTL related to multiple aspects of reproductive isolation: male fertilization success, offspring viability as a function of male identity and behavioural isolation of backcross males (egg production and whether or not mating occurred) with both L. goodei and L. parva females (sample size: Figure S2; phenotypic distribution: Figure S5–S7). Backcrossed males were paired separately with L. goodei and L. parva stimulus females. We performed QTL analyses separately for all traits. In total, we had six traits that we mapped in backcrossed males: whether or not a male mated with L. goodei, whether or not a male mated with L. parva, egg production with L. goodei, egg production with L. parva, male fertilization success and offspring viability. For each of these traits, the QTL mapping was done in rQTL using the backcross linkage map and scanone with standard interval mapping (‘em’) (Broman & Sen, 2009). We calculated the significance of LOD scores using 5,000 permutations and the 95% Bayesian credible interval for any significant QTL identified. Inclusion of backcross types (BC i‐iv) as a covariate or estimating LOD thresholds separately for the sex chromosome (Broman et al., 2006) did not qualitatively change any of our findings (Table S2).

RESULTS

Backcross map

The backcross linkage map contained a total of 172 placed markers on 23 linkage groups, which were numbered based on syntenic markers shared with the L. goodei and L. parva maps (see Berdan et al., 2014). In both the L. parva map and the backcross map, linkage group 1 represented a fusion of two linkage groups (1A and 1B) from L. goodei. The total length of our map was 448 cM, which is intermediate between the lengths of our previously published maps for L. goodei (605 cM) and L. parva (392 cM) (Berdan et al., 2014).

Loci linked to the sex‐determining region

There was support for an XX‐XY sex determination system that mapped to the fused chromosome. There was little support for a ZZ‐ZW system in cross 1, as the predicted ZZ‐ZW genotypes failed to map to the hybrid map (LOD < 1.32, p > .53; N = 44 informative individuals). In contrast, the XX‐XY predicted genotypes in cross 1 mapped to chromosome 1 at 0 cM near marker 05836 (LOD = 3.35, p = 0.02, 95% Bayesian Credible Interval 0–12 cM; N = 16 individuals; Figure S8a,b). Likewise, both crosses 2 and 3 provided strong support for an XX‐XY sex determination system that mapped to one the end of the fused chromosome (Figure 1; Figure S8). In cross 2 (second L. parva — L. goodei interspecies cross), the predicted XX‐XY genotypes mapped to chromosome 1A at 2 cM between markers 13121 and 14413 on the L. goodei map (LOD = 5.21, p < .001, N = 50; 95% Bayesian Credible Interval 0.5–3 cM; Figure S8c,d). Using these same data and the L. parva linkage map, the sex‐determining region mapped to chromosome 1 at 10.5 cM near marker 13,005 (LOD = 6.82, p < .001, N = 50; 95% Bayesian Credible Interval 9–11 cM; Figure 1a,b; a location of 12cM on chromosome 1 in the L. parva corresponds to a location 1cM on the backcross map; Table S3). Using cross 3 (L. parva population backcrosses) and the L. parva map, the predicted XX‐XY genotypes mapped to chromosome 1 at marker 11321 at 20.8 cM (LOD = 5.54, p < .001, N = 36; 95% Bayesian Credible Interval 13–44 cM; Figure 1c,d). Thus, the data supported an XX‐XY sex determination system that consistently mapped to the chromosome 1A portion of the fused chromosome in all 3 crosses.

Salinity tolerance

Survival in salt water was approximately half of that in fresh water (Figure S5a, Table S4, freshwater = 40.3 ± 0.037% SE, n = 29; saltwater = 24.4 ± 0.046%, N = 19; for families where 5 or more eggs were tested). This survival rate is similar to that found previously in Fuller et al. (2007) and Fuller (2008) for backcrosses to L. goodei (freshwater survival: 65%, saltwater survival: 32%).

In our experiment, salinity tolerance loci are diagnosable by an overabundance of heterozygotes (i.e. individuals with an L. parva allele) in salt water compared to fresh water. We compared the proportion of homozygous (L. goodei) and heterozygous genotypes at each marker between fresh and saltwater rearing conditions (Table 2). We examined QTL at the linkage group level (i.e. a maximum of one QTL per linkage group) and only considered linkage groups with more than one significant locus as being involved in adaptation to salinity as a way to avoid potential false positives that could be caused by genotyping errors. Linkage groups with multiple markers where heterozygotes were under‐represented in fresh water and over‐represented in salt water were as follows: 3, 6, 7, 12 and 17 (Table 2, Table S5, Figure 2). The effects were particularly strong for linkage group 7, where the heterozygotes were 1.9 times as abundant in salt water (~0.65) as they were in fresh water (~0.34). Loci at linkage group 16 showed the opposite pattern to our expectations, where heterozygous individuals were common among freshwater and rare among saltwater offspring. An additional 5 linkage groups (11, 18, 19, 21, 22) had a single significant marker implicated in salinity tolerance (Table 2, Table S5). Of these, linkage group 18 showed the expected pattern with an overrepresentation of heterozygotes in salt water. Linkage group 11,19, 21 and 22 showed the opposite pattern with an overrepresentation of heterozygotes in fresh water relative to salt water. Only linkage group 21 was implicated in intrinsic isolation, whereas linkage group 11 was implicated in behavioural isolation (see below).

Table 2

Salinity associated loci: loci showing significant different proportions of heterozygotes in freshwater versus saltwater

Marker	Linkage Group	Position (cM)	Proportion of heterozygotes in freshwater	Proportion of heterozygotes in saltwater	Chi‐square p‐value	FDR p‐value
09418	7	16.96	0.34	0.65	2.10E−09	3.01E−07
00141	7	17.49	0.34	0.64	8.42E−09	5.88E−07
14667	7	0.00	0.18	0.41	1.23E−08	0.0000
14398	17	24.83	0.44	0.68	0.0000	0.0005
18723	17	25.90	0.45	0.67	0.0000	0.0011
11877	21	0.00	0.54	0.33	0.0001	0.0027
00137	17	25.89	0.46	0.67	0.0001	0.0027
13073	6	0.00	0.48	0.68	0.0002	0.0034
13872	18	9.84	0.33	0.51	0.0003	0.0052
11937	3	0.00	0.34	0.52	0.0006	0.0079
15386	12	6.62	0.38	0.57	0.0006	0.0079
10789	3	0.00	0.35	0.52	0.0008	0.0092
05062	12	6.66	0.39	0.57	0.0008	0.0092
11514	3	0.08	0.34	0.51	0.0009	0.0094
14634	3	3.34	0.34	0.51	0.0012	0.0106
11023	12	6.11	0.41	0.58	0.0012	0.0106
10999	16	0.00	0.62	0.45	0.0013	0.0106
12180	19	7.45	0.57	0.40	0.0017	0.0132
11538	16	3.18	0.61	0.44	0.0018	0.0132
05635	3	0.00	0.36	0.52	0.0018	0.0132
09531	3	0.00	0.37	0.53	0.0020	0.0136
16266	11	34.49	0.59	0.43	0.0026	0.0162
11737	12	6.87	0.41	0.57	0.0026	0.0162
03290	3	0.00	0.36	0.52	0.0028	0.0168
15226	12	6.88	0.41	0.57	0.0038	0.0215
12650	6	16.40	0.40	0.55	0.0041	0.0227
07450	16	2.65	0.59	0.44	0.0053	0.0280
11889	22	10.68	0.57	0.43	0.0071	0.0364
06697	6	16.41	0.40	0.54	0.0100	0.0491

Results from chi‐square test with df = 1.

Abbreviations: cM, CentiMorgan; FDR, False discovery rate.

Figure 2

Salinity tolerance loci. Difference in proportion of heterozygous individuals in salt versus freshwater plotted for loci across all 23 linkage groups. Linkage group numbers listed above, position of loci in centiMorgans (cM) on the hybrid map shown. Different linkage groups separated by white partitions. Red lines indicate FDR cut‐offs. LG 3, 6, 7,12, 17 showed multiple outliers in the predicted positive direction. See Table 2for loci names

Gametic disequilibrium—interactions among loci

The chromosomal fusion was implicated in genetic incompatibilities. Survival of backcrossed offspring in fresh water was low (~40%), and the surviving backcrossed offspring were a nonrandom subset that had favourable combinations of alleles at different loci. Twenty‐six of 10,675 tests for interactions among genotypes at loci on different linkage groups remained significant after correcting for multiple tests. All of these interactions among loci involved an overrepresentation of offspring that had concordant genotypes (i.e. either homozygous for the L. goodei specific marker at both loci or heterozygous at both loci). Individuals with mixed genotypes (i.e. homozygous at one locus, but heterozygous at another) were either absent or under‐represented. Although there were 26 significant interactions, these involved loci on only five pairs of linkage groups (Table 3). There were multiple significant interactions involving loci on linkage group 1 and both linkage groups 13 and 16. One interaction between linkage group 1 and linkage group 13 involved a marker very close to the sex‐determining region (marker 13005). There were also significant interactions between linkage groups 13 and 16, linkage groups 21 and 22, and linkage groups 23 and 2. Due to our somewhat low sample size (N = 121), a two‐fold difference in the abundance of concordant genotypes versus mixed genotypes was required for a marker to be inferred as highly significant after FDR correction (Table S6).

Table 3

Genotypes between markers on different linkage groups with significant frequency distortion

Linkage groups		p‐value	# genotypes				Marker ID
LG Marker 1	LG Marker 2	FDR	AABB	AABb	AaBB	AaBb	Marker 1	Marker 2
1 ‐ sex	13	1.53E−23	57	0	0	64	13005	02161
2	23	1.53E−23	55	0	0	66	03425	11531
21	22	2.10E−23	57	0	0	62	02541	06333
2	23	2.10E−23	54	0	0	65	17258	11531
2	23	3.67E−23	55	1	0	65	15948	11531
2	23	3.67E−23	54	0	1	66	14340	11531
21	22	8.33E−23	57	1	0	61	23270	06333
21	22	4.27E−22	57	0	2	60	03323	06333
21	22	4.27E−22	57	0	2	60	03556	06333
21	22	6.34E−22	57	0	2	59	03555	06333
22	21	2.46E−21	56	2	1	60	06333	06712
1*	13	1.26E−18	44	0	0	51	10924	02161
1 ‐ sex	13	5.93E−11	46	8	8	52	11211	02161
1 ‐ sex	13	8.93E−07	40	7	14	45	11521	02161
13	16	0.0066	38	16	19	48	02161	12642
1 ‐ sex	16	0.0066	38	19	16	48	13005	12642
13	16	0.012	37	16	20	48	02161	04992
13	16	0.012	37	16	20	48	02161	17027
1 ‐ sex	16	0.012	37	20	16	48	13005	04992
1 ‐ sex	16	0.012	37	16	20	48	13005	17027
16	13	0.025	36	21	16	48	01506	02161
1 ‐ sex	16	0.025	36	21	16	48	13005	01506
13	16	0.036	36	17	19	45	02161	13269
1 ‐ sex	16	0.036	36	17	19	45	13005	13269
13	16	0.046	35	16	22	48	02161	11538
1 ‐ sex	16	0.046	35	22	16	48	13005	11538

Genotypes refers to the number of individuals that are homozygous for the L. goodei marker at both loci (AABB), are heterozygous at both loci (AaBb), or are homozygous at one locus but heterozygous at another (AABb and AaBB). 1‐sex indicates a marker on linkage group1 located within the sex determining region (marker ID in bold); 1* indicates a marker on linkage group 1 that is adjacent to the sex determining region (marker ID in italics). p‐value from a chi‐square test with df = 1.

Abbreviation: FDR, False discovery rate.

Behavioural isolation

As expected, behavioural isolation was higher when backcrossed males were paired with L. parva females than with L. goodei females. We found that 75% of males mated at least once with L. goodei females during the 20‐day period, but only 48% of males mated with L. parva females. Similarly, backcrossed males had higher daily egg production with L. goodei (1.78 eggs/day) than they did with L. parva (0.61 eggs/day) (Figure S6).

The number of eggs produced in pairings with L. goodei females mapped to chromosome 1 at 32 cM (LOD = 2.89, p = .01, N = 29; Figure 3a). Males that were homozygous for the L. goodei allele had higher mating success than heterozygotes that carried an allele from L. parva (homozygotes: 3.41 eggs/day, heterozygotes: 0.77 eggs/day; Figure S9). There were no QTL identified for the probability of a male mating with L. goodei females, which is not surprising as most males successfully mated with L. goodei.

Figure 3

LOD scores from QTL mapping of behavioural isolation. (a) number of eggs produced when male mated toL. goodei,(b) number of eggs produced when male mated toL. parva, C) male probability of mating withL. parva.Red dashed line indicates LOD score = 2.5 (approximately equivalent top < .05)

Mating success with L. parva females mapped to two linkage groups, but the genotypic pattern did not match our expectation that heterozygotes should have increased success with L. parva (backcrosses could only be heterozygous for L. parva alleles, not homozygous). Whether or not a male mated with L. parva females mapped to linkage group 1, marker 13870 at 57 cM (LOD = 2.87, p = .021, N = 33; Figure 3b) and heterozygotes were less likely to mate with female L. parva than were homozygotes (90% of homozygotes mated with L. parva compared to 27% of heterozygotes; Figure S9). Likewise, the number of eggs laid when males were mated to L. parva females mapped to chromosome 11 at 16.5 cM (LOD = 6.2, p = .01, N = 33) (Figure 3c), but the pattern was one where heterozygotes had lower egg production with L. parva females than individuals that were homozygous for the L. goodei allele (Figure S9).

Fertilization success and offspring viability as a function of male genotype

Fertilization success and offspring viability did not differ between males when mated L. goodei or L. parva females (fertilization success: L. goodei females = 0.52 ± 0.12 SE, L. parva females 0.57 ± 0.12 SE; offspring viability: L. goodei females = 0.70 ± 0.09 SE, L. parva females = 0.77 ± 0.11). Fertilization success (proportion of fertilized eggs) of backcrossed males was generally bimodal with males having high or low success (Figure S7). Fertilization success of backcrossed males mapped to a single QTL located on linkage group 7 at 25 cM (LOD = 4.15, p = .03, N = 23; Figure 4a). Male fertilization success was 2.6 times higher in heterozygotes (0.64) than L. goodei homozygotes (0.24) (Figure S9). Offspring viability was defined as the proportion of fertilized eggs surviving to hatching. Offspring viability as a function of paternal genotype mapped to linkage group 1 at 9 cM (LOD = 3.47, p = .04, N = 20; Figure 4b) in the same region where sex determination maps and offspring viability was 2.4 times higher in homozygotes (0.78) than in heterozygotes (0.32) (Figure S5, Figure S7).

Figure 4

LOD scores from QTL mapping of male fertilization and offspring viability. (a) Male fertilization success, (b) male offspring viability. Red dashed line indicates the LOD score = 2.5, approximately equivalent top < .05

DISCUSSION

In this study, we investigated the role of salinity adaptation in the evolution of behavioural and intrinsic isolation in the killifish Lucania goodei and L. parva by genetically mapping salinity tolerance, sex determination, behavioural isolation and intrinsic isolation. We found that salinity tolerance has a polygenic basis, mapping to multiple linkage groups, but with little co‐localization with other reproductive isolating barrier loci. This pattern provides little evidence that selection via physical linkage or pleiotropy has contributed much to the evolution of behavioural or intrinsic isolation among Lucania species. This study is consistent with adaptation to salinity contributing to reproductive isolation through habitat isolation, extrinsic isolation and/or immigrant inviability. We also found evidence that a fusion between an autosome and the chromosome with the sex‐determining region has occurred in L. parva. Several of the genetic incompatibilities we mapped were located on the fused sex chromosome, with L. parva alleles causing gametic disequilibrium (segregation distortion) and offspring inviability (reduced hatching success of offspring). These results provide preliminary evidence that the chromosomal fusion may harbour loci related to reproductive isolation. Below, we discuss these results in more detail.

Adaptation to salinity and speciation in lucania

Salinity tolerance mapped to numerous locations in the Lucania genome. Admittedly, some of these loci may be false positives. Still, salinity tolerance in Lucania is almost certainly polygenic as decades of research have revealed that salinity tolerance in teleosts is a complex trait that involves multiple tissues (e.g. gills, kidneys) and physiological pathways (Evans, 2008; Evans et al., 2005; Larsen et al., 2011; Laverty & Skadhauge, 2012). Indeed, our own previous work examining both gene expression and coding changes in Lucania indicates that salinity tolerance is a multi‐locus trait (Berdan & Fuller, 2012b; Kozak et al., 2014). Other studies of the genomic basis of salinity tolerance in teleosts have revealed similarly distributed genetic architectures. Lucania is a member of the Fundulidae family (Whitehead, 2010). In Fundulus heteroclitus, a genome‐wide association study found salinity tolerance loci on at least 5 different chromosomes (Brennan et al., 2018). In salmonids, between 3 (in Oncorhynchus mykiss) and 10 (Salmo salar and Salvelinus alpinus) linkage groups are involved in salinity tolerance (Norman et al., 2012). Loci associated with salinity map to 11 different linkage groups in Atlantic cod (Gadus morhua) (Berg et al., 2015). It is not yet clear if this type of genetic architecture will facilitate or hinder divergence during speciation with gene flow. If many targets of divergent selection are spread across the genome, this increases the potential for the coupling of adaptive loci with other isolating barriers; however, the effect of divergent selection at each locus will be diluted. Therefore, whether polygenic adaptation leads to the evolution of other reproductive isolating barriers via linked selection should depend on the number of loci, their physical positions and the strength of selection.

Our current study provides little support for the hypothesis that the force of divergent natural selection has played a major role in the evolution of behavioural and intrinsic isolation. Of the five linkage groups where heterozygotes with L. parva alleles survived better at high salinity, only linkage group 7 affected both salinity tolerance and male fertilization success. However, the pattern in fertilization success was not consistent with this linkage group increasing reproductive isolation. Rather than L. parva alleles conferring low fertility, backcrossed individuals that were homozygous for L. goodei alleles on linkage group 7 were more likely to be infertile. Our finding that areas of the genome associated with salinity tolerance differ from those contributing to other forms of reproductive isolation is in keeping with prior evidence that sympatry between L. goodei and L. parva led to increased behavioural isolation among L. parva populations, regardless of salinity habitat (Kozak et al., 2015). Novel adaptations can become coupled with genetic incompatibilities due to spatial subdivision of populations as indicated by theoretical models (Bierne et al., 2011). Given that karyotypes differ among species and intrinsic isolation maps to the fused chromosome, our current data cannot rule out an incompatibility‐first model, with coupling of salinity tolerance and incompatibilities in Lucania occurring later.

Two caveats deserve mention before the hypothesis of reproductive isolation as an indirect outcome of salinity adaptation is rejected. First, this study had a relatively low sample size for a QTL study. Hence, only loci with very large effects could be detected and small effects can be important in reproductive isolation (Ellison et al., 2011). Furthermore, sparse marker coverage may have led to missed QTL, broader intervals, and the potential of considering multiple small effect QTL as a single large effect QTL. Second, this study mapped salinity tolerance only at the earliest life‐history stage of egg/larval survival. We chose this life stage because previous work has shown that both species can tolerate fresh water well in the egg/larval fish stage and that L. goodei have reduced survival in salt water in this stage relative to L. parva. However, juvenile and adult L. parva have reduced survival relative to L. goodei when they are forced to overwinter in fresh water (Fuller et al., 2007). We did not map salinity tolerance at later life stages, although these may conceivably be important to the evolution of other reproductive isolating barriers. Given these caveats, our current understanding of the Lucania system is that differential adaptation to salinity appears to contribute to habitat isolation (Fuller & Noa, 2008), but does not result in correlated responses in other traits that promote behavioural or intrinsic isolation.

Sex determination, chromosomal rearrangements and speciation

The male sex‐determining region (XX‐XY) mapped to the fused chromosome in both hybrid (L. goodei × L. parva) and pure L. parva crosses. This suggests that the centric fusion in L. parva occurred between the Lucania chromosome with the sex‐determining region and an autosome. Combined with the pattern of reduced fertility in F1r males (with L. parva mothers), but not F1r females (Fuller, 2008), this finding is consistent with Haldane's Rule of the heterogametic sex being more likely to exhibit sterility or inviability (Turelli & Moyle, 2007). Although all karyotypic evidence points to the Lucania sex chromosomes being homomorphic (Berdan et al., 2014; Uyeno & Miller, 1971), it is possible that F1r males are potentially hemizygous for some loci in the male determining region itself (Peichel et al., 2020), particularly if there has been some change in the genetic content of the chromosome accompanying the fusion. In fishes, sex chromosomes are often involved in fusions because of sexually antagonistic selection or male‐mutation bias (Bachtrog, 2013; Bachtrog et al., 2011; Kitano & Peichel, 2012; Pennell et al., 2015). However, unlike many other known fusions in fish (Kitano & Peichel, 2012; Kitano et al., 2009), our fusion does not appear to represent a neo‐Y system with unfused × chromosomes, because both males and females possess fused chromosomes (Berdan et al., 2014). Future work can investigate the differences that have arisen in between the autosome in L. goodei and the homologous portion in L. parva that is now sex‐linked, similar to previous work with Drosophila species (Counterman, Ortíz‐Barrientos, & Noor, 2004).

The chromosomal fusion may contribute to reproductive isolation between L. goodei and L. parva as several genetic incompatibilities map there (Figure 5; Table 4). Chromosomal fusions often differentiate populations and have been shown both theoretically and empirically to facilitate adaptation (Dobigny et al., 2017; Franchini et al., 2010; Guerrero & Kirkpatrick, 2014; Wellband et al., 2019). QTLs for behavioural isolation mapped to the fused chromosome, including number of eggs laid with L. goodei and the probability of mating with L. parva females. We found that the fused chromosome contained a locus associated with offspring viability. Additionally, loci on the fused chromosome appeared to interact with loci on linkage groups 13 and 16 to cause an incompatibility (segregation distortion). Overall, of 6 reproductive isolating barriers mapped, 4 mapped to the fused chromosome. However, this may be a result of a mapping bias caused by regions of low recombination when QTL loci are spread across the genome (Noor et al., 2001). Further work is needed to clarify the role of the fusion to reproductive isolation in Lucania. We also identified an incompatibility segregating between two other pairs of linkage groups (2 with 23 and 21 with 22). The interaction between linkage group 21 and 22 is interesting because it involves markers that mapped to a putative translocation. These markers map to linkage group 21 in one species and linkage group 22 in the other (Berdan et al., 2014). Our data suggest there are at least several segregating incompatibilities among multiple linkage groups that contribute to the reduced fitness previously observed in backcrosses to L. goodei (Fuller, 2008; Fuller et al., 2007).

Figure 5

Sex‐determining and isolating loci mapping to Linkage group 1. Bayesian credible intervals for sex determination and isolating traits (solid rectangles) mapped relative to position (in cM) along linkage group 1 (the fused chromosome) from the hybrid linkage map. Blue indicates hybrid incompatibility loci; green indicates behavioural isolation; sex determination region shown in red. The ancestral autosomal portion (chromosome 1B) is located at ~ 40–57 cM on this hybrid map. This figure was made with LinkageMapView (Ouellette et al., 2017).

Table 4

Summary of locations of isolating barriers. Grey shading indicates that a QTL for the specified trait localized on the specified linkage group.

Linkage Group	Intrinsic Isolation		Behavioural Isolation				Adaptation	Sex determination
	Gametic Disequilibrium (Incompatibilities)	Male Fertilization Success	Male Offspring Viability	Eggs Laid with L. parva ♀	Probability of Mating with L. parva ♀	Eggs Laid with L. goodei ♀	Differential Survival in Fresh versus Saltwater	Sex‐Determining Locus
1					*
2
3
4
5
6
7		*
8
9
10
11				*
12
13
14
15
16							*
17
18
19
20
21
22
23

A * indicates that the expected pattern was reversed.

Finally, we note that intrinsic isolation is asymmetric in Lucania, in line with Darwin's corollary to Haldane's rule (Turelli & Moyle, 2007). Additionally, in the F1 generation, F1r (L. parva ♀ × L. goodei ♂) males suffer a 50% reduction in early egg survival relative to parental males and F1 (L. goodei ♀ × L. parva ♂) males (Fuller, 2008). We did not map this incompatibility in this study, although the incompatibility was expressed in our F1 generation, with some F1r males having no fertility and being excluded as sires (Table S4). It is possible that the incompatibility was selected against in backcross iv (F1r ♂ × L. goodei ♀), as males from this cross tended to have higher male fertilization success and offspring viability than male offspring from the other three backcrosses. Unfortunately, we do not have the power to determine if these reductions in offspring viability represent the same incompatibilities across generations here. The major takeaway is that we have mapped loci important to reproductive isolation preventing gene flow from L. parva into L. goodei, but the reproductive barriers preventing gene flow from L. goodei into L. parva remain unknown.

Implications for reinforcement

Potential linkage of behavioural and intrinsic isolation on the fused chromosome is intriguing given prior work in Lucania, suggesting that reduced hybrid fitness may have caused increased species‐specific preferences in areas where the two species co‐occur (Gregorio et al., 2012; Kozak et al., 2015). This process is often referred to as reinforcement, where costly mating with heterospecifics generates selection for altered mating traits and preference (Servedio & Noor, 2003). Reinforcement can only occur when post‐zygotic isolation is present. Lucania goodei and L. parva are fairly genetically diverged with an estimated genome‐wide F of 0.38 and intrinsic post‐zygotic incompatibilities are present (Fuller, 2008; Kozak et al., 2014), suggesting that these two species have been isolated for a significant amount of time. We hypothesize that the initial stages of divergence likely occurred via local adaptation and habitat isolation, after which many of the reproductive isolating barriers that we mapped in this study may have evolved.

Our findings combined with our previous work (Gregorio et al., 2012; Kozak et al., 2015) add to the growing support for the importance of chromosomal rearrangements involving the sex‐determining region as facilitators of reinforcement. In order for reinforcement to occur, incompatibilities and behavioural isolation must be coupled with one another (Butlin & Smadja, 2018; Ortíz‐Barrientos, Grealy, & Nosil, 2009). Physical linkage, with subsequent reduced recombination between loci, is one of the strongest ways to generate coupling. Fusions can generate linkage disequilibrium in two ways: by bringing previously unlinked loci together and by reducing recombination (Dumas & Britton‐Davidian, 2002; Franchini et al., 2010). An example of this occurs in Japan Sea sticklebacks (Gasterosteus nipponicus), where the Y‐chromosome has fused to an autosome containing a behavioural isolation locus (Kitano et al., 2009). We found preliminary evidence that major effect QTL for behavioural isolation and hybrid incompatibilities may be located on the fused sex chromosome in Lucania. However, additional studies with larger sample sizes and an enhanced ability to detect QTL of smaller effect are needed to determine the true effect of our detected QTL, as Beavis effects may have upwardly biased these (Xu, 2003), and to rule out false positives. Such studies are needed to definitively determine whether incompatibility loci are linked to behavioural isolation loci in Lucania.

Both of our QTL for mating success with L. parva were in the opposite direction than anticipated. Individuals homozygous for L. goodei alleles had higher success than heterozygous individuals. There are several possible explanations for this. First, these QTL could be false positives caused by a small sample size and the difficulties of properly measuring behavioural responses. Second, an L. parva allele may potentially only increase mating success with L. parva females when mating in saltwater environments. However, Berdan and Fuller (2012b) investigated the effect of salinity on behavioural isolation between L. goodei and L. parva and found little effect of salinity. Third, behavioural sterility could be occurring. Behavioural sterility occurs when admixed individuals fail to mate due to a defect (neurological, pheromonal, etc.) that renders them incapable of correctly courting the parental species (Coyne & Orr, 2004; Naisbit et al., 2001; Stratton & Uetz, 1986). Notably, our crossing design allowed backcross individuals to be homozygous for L. goodei alleles but only heterozygous for L. parva alleles. One possibility is that being heterozygous for L. parva alleles in an L. goodei genetic background is detrimental for effectively courting any female; a defect that would only manifest in certain backcrosses. Preferences of males and females show evidence for character displacement in L. parva and L. goodei (Gregorio et al., 2012; Kozak et al., 2015); however, the traits upon which these preferences are based have not yet been identified. Additional investigation of the male trait differences associated with L. parva genotypes and preferences of L. parva females against admixed males will be needed to clarify this result.

In summary, previous work in Lucania indicates that differences in salinity tolerance have profound effects on habitat distribution (Fuller & Noa, 2008), ecological interactions (Dunson & Travis, 1991), survival across different life stages (Fuller, 2008; Fuller et al., 2007), physiology (Berdan & Fuller, 2012b) and genomic differentiation (Kozak et al., 2014) in Lucania. Salinity tolerance in L. parva from the egg stage to small juvenile is polygenic and distributed across the genome. However, there is currently no evidence that adaptation to salinity contributes to the evolution of reproductive isolating barriers that prevent L. goodei and L. parva from collapsing into a hybrid swarm in areas where they co‐occur. Instead, a chromosomal rearrangement has created a fused chromosome in L. parva that contains both the sex‐determining region and loci that contribute to some forms of reproductive isolation. Theory predicts reinforcement is more likely to occur when incompatibility loci and mate preference loci are located on sex chromosomes than when they are on autosomes, because of reduced recombination and the greater efficiency of natural selection on sex chromosomes (as recessive alleles are exposed in the heterogametic sex) (Hall & Kirkpatrick, 2006; Lemmon & Kirkpatrick, 2006; Servedio & Saetre, 2003). Our results suggest that ecological adaptation as a function of salinity may be important in the early stages of differentiation in allopatry, but that chromosomal rearrangements involving the sex chromosome may contribute to reinforcement when species come back into secondary contact.

Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13725.

Fig S1‐S9

Click here for additional data file.

Table S1‐S6

Click here for additional data file.

Supplemental methods

Click here for additional data file.