Patterns of genomic divergence and introgression between Japanese stickleback species with overlapping breeding habitats

With only a few absolute geographic barriers in marine environments, the factors maintaining reproductive isolation among marine organisms remain elusive. However, spatial structuring in breeding habitat can contribute to reproductive isolation. This is particularly important for marine organisms that migrate to use fresh‐ or brackish water environments to breed. The Japanese Gasterosteus stickleback species, the Pacific Ocean three‐spined stickleback (G. aculeatus) and the Japan Sea stickleback (G. nipponicus) overwinter in the sea, but migrate to rivers for spawning. Although they co‐occur at several locations across the Japanese islands, they are reproductively isolated. Our previous studies in Bekanbeushi River showed that the Japan Sea stickleback spawns in the estuary, while the Pacific Ocean stickleback mainly spawns further upstream in freshwater. Overall genomic divergence was very high with many interspersed regions of introgression. Here, we investigated genomic divergence and introgression between the sympatric species in the much shorter Tokotan River, where they share spawning sites. The levels of genome‐wide divergence were reduced and introgression was increased, suggesting that habitat isolation substantially contributes to a reduction in gene flow. We also found that genomic regions of introgression were largely shared between the two systems. Furthermore, some regions of introgression were located near loci with a heterozygote advantage for juvenile survival. Taken together, introgression may be partially driven by adaptation in this system. Although, the two species remain clearly genetically differentiated. Regions with low recombination rates showed especially low introgression. Speciation reversal is therefore likely prevented by barriers other than habitat isolation.

organism habitat isolation can act before any other isolating barriers, such as assortative mating and post-zygotic isolation (Coyne & Orr, 2004). Differential use of breeding sites can also drive divergent adaptation and secondarily reduce hybridization (Schluter, 2000).
For example, once two species adapt to different breeding habitats, immigrants from dissimilar environments will have lower survival rates or reproductive success in foreign sites, thereby reducing the chances of heterospecific mating between native residents and immigrants (Nosil, Vines, & Funk, 2005).
Marine environments present a conundrum in this regard; there are seemingly high levels of dispersal and few geographical barriers to gene flow. The factors driving speciation and maintaining reproductive isolation in marine organisms remain elusive (Knowlton, 1993;Palumbi, 1994;Potkamp & Fransen, 2019;Puebla, 2009). Some marine organisms, such as anadromous fishes, use brackish or freshwater environments for breeding and closely related species may diverge in breeding habitat use (Hendry & Stearns, 2004;Kitano et al., 2009;Kume, Kitano, Mori, & Shibuya, 2010). For migratory species with divergent spawning preferences, variation in the spatial distribution and segregation of breeding sites within and among river systems may modulate the extent to which secondary contact occurs. What are the potential consequences of the breakdown of habitat isolation in this scenario? One possible outcome would be speciation reversal, where two species collapse into a single species by extensive hybridization (Elmer, 2019;Seehausen, Takimoto, Roy, & Jokela, 2008;Taylor et al., 2006). A second possibility is that a level of reproductive isolation is maintained in the face of ongoing gene flow. In this case, introgression is likely to be suppressed near loci involved in reproductive isolation (Elmer, 2019;Nosil, 2012;Ravinet et al., 2017), whereas alleles that are beneficial to both populations (Anderson et al., 2009;Pardo-Diaz et al., 2012;Whitney, Randell, & Rieseberg, 2006) or those with heterozygote advantage are likely to introgress elsewhere in the genome (Kim, Huber, & Lohmueller, 2018;Nadachowska-Brzyska, ZieliŃSki, Radwan, & Babik, 2012;Wolfe et al., 2019). Finally, if reproductive isolation is strong enough, breakdown of habitat isolation would neither influence the hybridization rates nor erode overall genomic divergence between species. Such strong reproductive isolation may be achieved by other multiple isolating barriers and/or irreversible isolating barriers, such as complete hybrid sterility and inviability (Nosil, Harmon, & Seehausen, 2009).
The Japanese anadromous stickleback species pair composed of the Japan Sea stickleback (Gasterosteus nipponicus) and Pacific Ocean populations of the three-spined stickleback (G. aculeatus; hereafter referred to as Pacific Ocean stickleback) are an excellent system to address this question. Our previous demographic analyses using whole genome sequence data showed that the Japan Sea and the Pacific Ocean stickleback species diverged about 0.68-1 million years ago (Ravinet et al., 2018), most likely triggered by sea level fluctuation during the Pleistocene (Mix et al., 1995). Furthermore, gene flow between the two species has probably been maintained for most of their evolutionary history (Ravinet et al., 2018). These two anadromous stickleback species share spawning sites in rivers at several locations across the Japanese islands (Ishikawa et al., 2019;Kume et al., 2010;Ravinet, Takeuchi, Kume, Mori, & Kitano, 2014;Yamada, Higuchi, & Goto, 2007). In Bekanbeushi River on Hokkaido Island, where our previous work focused, the most downstream and upstream of stickleback spawning sites are over 5 km apart (Figure 1), contributing a substantial level of habitat isolation (Kitano et al., 2009;Kume, Kitamura, Takahashi, & Goto, 2005;Kume et al., 2010). Isolation of spawning sites contributes to more than 90% of total reproductive isolation between these two species in the Bekanbeushi River (Kitano et al., 2009;Kume et al., 2010). Both species overwinter in the sea but migrate to the river to reproduce.
The Japan Sea stickleback spawns in the brackish estuary, while the Pacific Ocean stickleback spawns further upstream in a more freshwater environment (Kitano et al., 2009;Kume et al., 2010). In the midstream, there is a hybrid zone, where behavioural isolation and hybrid male sterility likely play an important role in reproductive isolation (Kitano, Mori, & Peichel, 2007;Kitano et al., 2009).
Introgression rates have been previously shown to be lower on the neo-X chromosomes than in autosomes (Ravinet et al., 2018).
In contrast to Bekanbeushi River, in the nearby Tokotan River Here, our first aim was to show how the presence and absence of breeding habitat isolation influence the genome-wide patterns of divergence and introgression. To this end, we conducted whole genome resequencing of sticklebacks caught in Tokotan River and compared the data with our previously published whole genome sequences of sticklebacks collected from Bekanbeushi River (Ravinet et al., 2018;Yoshida et al., 2014). We specifically compared the levels of overall interspecies genomic divergence between the two river systems. Our second aim was to test whether the introgression identified may have any functional significance. To this end, we examined whether introgression sites are shared by two systems and overlap with any fitness-related QTLs. Our previous QTL mapping of juvenile survival using an interspecies cross derived from the Bekanbeushi River system found several loci showing heterozygote advantage (Ishikawa et al., 2019). We tested whether these QTL overlapped with genomic signatures of introgression between the two species.
Finally, we asked whether regions with lower recombination rates show less introgression.

| Fish sampling and genome sequencing
Breeding sticklebacks were distributed from the river mouth  (Kume et al., 2005). These fish turned out to be nine females and two males of the Japan Sea stickleback and eight females and 11 males of the Pacific Ocean stickleback, based on morphological characteristics (Higuchi, Sakai, & Goto, 2014;Kitano et al., 2007 (Ravinet et al., 2018), compared to the much lower coverage used for samples from Tokotan (see above). In order to address this, we used samtools 0.1.19 (Li et al., 2009) to downsample all bams to adjust the mean depth to 12× to ensure similar levels of coverage among individuals and to prevent variant calling biases.
Following downsampling, we called variants at all positions across the genome using the bcftools 1.9 multiallelic caller (Danecek & McCarthy, 2017).
Following variant calling, we filtered our data set using vcftools 0.1.17 (Danecek et al., 2011). We removed all sites with >20% missing data, a site Phred quality score <30, a minimum genotype depth of 5×, a maximum genotype depth of 25×, a minimum mean site depth of 5× and a maximum mean site depth of 25×. Further filters, such as minor allele frequency or minimum number of alleles were applied downstream when required for specific analyses and are detailed in the relevant sections.
Prior to analysis, we used the LiftoverVcf module of Picard 2.9 to lift our variant calls to the linkage map-corrected physical positions of the stickleback genome (Roesti, Moser, & Berner, 2013).
We also phased a subset of our data containing only biallelic SNPs using Shapeit2 (O'Connell et al., 2014). Phasing was informed by the recombination map estimated by Roesti et al. (2013) and was conducted using default parameters.
We combined our resequencing data from 20 individuals from when included, we used only females (i.e., X chromosomes) to exclude the confounding effects of intraspecific genetic differentiation between X and Y chromosomes.

| Population structure and admixture
To identify the major axes of differentiation among our sequenced individuals, we used principal components analysis (PCA) implemented in plink 1.9 (Purcell et al., 2007). PCA is a parameter-free approach that has no a priori assumptions of population structure.
Prior to population structure analysis, we only allowed sites with a <5% missing data. Since PCA also requires independence among sites, we also performed linkage pruning, removing all sites with an r 2 > .1. This resulted in 282,612 loci. We did not include P. pungitius for PCA.
In addition to PCA, we used ADMIXTURE, a model-based approach to assessing population structure in a Maximum Likelihood framework (Alexander, Novembre, & Lange, 2009). In short, ADMIXTURE attempts to fit a model of K population clusters to the data, fitting individuals in a way that minimizes deviation from Hardy-Weinberg equilibrium. Like PCA, ADMIXTURE also requires independence among sites, so we reused the linkage-pruned data set for this analysis. We ran ADMIXTURE for K = 1-10 and also performed cross-validation error estimation in order to assess the most suitable value of K (Alexander & Lange, 2011). We did not include P. pungitius in the ADMIXTURE analysis.
We additionally investigated individual variation in admixture and ancestry by estimating hybrid index using the R package introgress (Gompert & Alex Buerkle, 2010). For this, we used Japan Sea and Pacific Ocean individuals from the Bekanbeushi River as 'pure' parental populations. Our rationale for this was that these individuals were previously assigned as pure ancestry based on microsatellite markers, prior to genome resequencing (Kitano et al., 2009).
We first identified autosomal sites with an absolute allele frequency difference of greater than 0.9 between the two populations and then performed linkage pruning (as above) to reduce our data set to 5,044 ancestry informative loci. We then used these sites to calculate hybrid index and interspecific heterozygosity among Tokotan individuals. We additionally used adegenet (Jombart & Ahmed, 2011) to simulate F 1 , F 2 and backcross hybrids in order to assess whether any Tokotan individuals fell close to these classes on the joint distribution of interspecific heterozygosity and hybrid index.

| Genetic differentiation, divergence and introgression
We estimated genetic differentiation (F ST ) and divergence (d XY ) using the popgenWindows.py script (https://github.com/simon hmart in/ genom ics_general). F ST calculated with this script is also known as γ ST (Hahn, 2019). Both statistics were estimated in 10 and 100 kb sliding windows with a 2.5 and 25 kb step respectively. We note here that sliding window estimates were generated for visualization purposes only; we removed overlapping windows to prevent pseudoreplication in downstream statistical analyses. Mean-pairwise genomewide estimates of these statistics were thus derived from the 10 kb nonoverlapping window data using the dplyr package in R (Wickham, François, Henry, & Müller, 2018). We did not include Little Campbell River fish and P. pungitius for analysis of F ST and d XY .
We further examined admixture and introgression using Patterson's D statistic implemented in ADMIXTOOLS (Patterson et al., 2012) and the R package admixr (Petr, Vernot, & Kelso, 2019). For the last test, we set P1 and P2 to JS-B and JS-T respectively in order to investigate site sharing with all possible G. aculeatus populations separately set as P3. By using multiple configurations of the D statistic with different populations with known divergence times (Ravinet et al., 2018), our rationale was that we would be able to broadly infer when introgression likely occurred and also gain some insight into the direction of gene flow. The significance of all D statistics was calculated using block jack-knifing to determine a Z-score distribution, where an absolute value of Z > 3 is equivalent to an approximate p-value of .001.

| Characterization of introgressed regions
To investigate the physical location of introgression in the genome, we used f d , a derivative of ABBA-BABA statistics that is reliable and robust to sliding window estimates (Martin, Davey, & Jiggins, 2014).
f d was calculated in 10 and 100 kb sliding windows with a 2.5 and 25 kb step respectively using ABBABABAwindows.py (https://github. com/simon hmart in/genom ics_general). As for estimates of genetic differentiation and divergence, sliding windows were generated for visualization purposes only; we removed overlapping windows to prevent pseudoreplication in downstream statistical analyses.
In  This enabled us to count the number of unique peaks and to measure the average size of f d peaks observed.
One caveat to both f d and ABBA-BABA statistics is that they assume monophyletic relationships within the groups at the tips of the test phylogeny. An alternative method that attempts to account for this is topology weighting-that is, the weighting of different topologies in a sliding window across the genome. Using our phased, biallelic site data set, we first used phyml (Guindon et al., 2010) to estimate phylogenies in 50 SNP windows across the genome for G. aculeatus Ocean from different river catchments group more closely with one another (Figure 4a). We then used TWISST to weight the contribution of each of these topologies to all the possible subtrees at each window across the genome (Martin & Van Belleghem, 2017).
In order to examine how introgression varied with recombination rate, we used a previously published three-spined stickleback linkage map (Roesti et al., 2013). Following Roesti et al. (2013), we fit a loess model with a span of 0.25 in order to interpolate recombination rate for 10 kb windows across the genome. We then compared mean recombination rate for f d peaks versus the genome background using a one-tailed permutation test, resampling 10,000 times.

| Population structure and admixture
PCA revealed the primary axis of divergence, which explained 33.3% of variance, was between Japan Sea stickleback and all G. aculeatus (Pacific Ocean and Atlantic; Figure 2a and Figure S1). PC2 explained 6.16% of variance and reflected differentiation within G. aculeatus among the marine populations-the Japanese Pacific Ocean, North American Pacific Ocean and the Atlantic lineages (Figure 2a and Figure S1). Notably, Pacific Ocean individuals from Tokotan were closer in PC space to Japan Sea individuals from the same river system, especially in comparison to the same species pair from Bekanbeushi River; suggesting a possibility of a higher level of admixture between the two species in the Tokotan system.
To further investigate this, we used ADMIXTURE to characterize population structure. A value of K = 2 separated all G. aculeatus from G. nipponicus (see Figure S2). However, a value of K = 3 was best supported ( Figure S3) and this further split G. aculeatus marine fish into Atlantic and Pacific Ocean lineages. In both K = 2 and K = 3, ADMIXTURE showed that Pacific Ocean and Japan Sea individuals from Tokotan were more admixed than those sampled from the Bekanbeushi River system (Figure 2b). Most notably, two Tokotan Pacific Ocean females with the highest PC1 scores also had higher proportions of Japan Sea ancestry from the ADMIXTURE analysis (Figure 2a, b;individuals No. 38 and No. 39). Further analysis using ancestry informative markers showed that these two individuals have a higher level of interspecific heterozygosity and a lower hybrid index, indicating that they might be recent backcrosses ( Figure S4).

| Differentiation, divergence and introgression
Pairwise mean F ST was highest between Japan Sea fish from Bekanbeushi and Atlantic Ocean stickleback (0.618, see Figure 2c).
Focusing on pairwise comparisons from our two Japanese river systems, mean F ST was substantially higher between the Pacific Ocean and Japan Sea fish in Bekanbeushi (0.617) than in Tokotan (0.476, one-tailed permutation test: p < 1 × 10 -4 ). Similarly, genomic divergence measured as d XY was lower in Tokotan (0.0182) compared to Bekanbeushi (0.0199: p < 1 × 10 -4 ; Figure 2d). Genomewide F ST in 100 kb sliding windows also clearly demonstrated that interspecies genomic differentiation is lower in Tokotan compared to Bekanbeushi (Figure 3). This pattern was less apparent with d XY , except on chromosome XIX, the Ancestral-sex chromosome (see Figure S5).
To more formally test for the occurrence of introgression, we

| Characterizing the introgression sites
We first investigated introgression regions using topology weighting with TWISST ( Figure 4). Across the genome, the species topology had the highest mean weight (0.926) as expected for species with high levels of genomic divergence. However, the next highest weight was the geographical topology (0.0402) which indicates introgression between the species in one or both river systems (Figure 4b).
We next investigated finer-scale patterns of introgression with the aim to identify where in the genome gene flow might be occurring. We did this by calculating f d , a statistic derived from ABBA-BABA site frequencies that estimates the proportion of a genome window where sites are shared between two focal taxa. By setting P1 to G. aculeatus from the Atlantic Ocean and the outgroup as P. pungitius, higher values of f d suggest an excess of allele shared between the Japanese species pairs. In accordance with population structure fied, 348 occurred in both systems, and the overlap was significantly higher than expected by random sampling (p = 0; hypergeometric test). The mean size of peaks did not differ significantly between the two comparisons (p = .051), but was slightly larger in Bekanbeushi (33.7 kb) compared to Tokotan (29.0 kb). In general, 10 kb windowed estimates of f d from between-species comparisons in both systems were highly correlated (r = .72, p < 2.2 × 10 -16 ; Figure S6), suggesting either a shared history of introgression or multiple independent introgression events in two systems; nonetheless, a substantial number of windows showed high f d in Tokotan but not Bekanbeushi.
Theory and empirical observations suggest introgression more likely to be observed when recombination rate is high as linked neutral or adaptive alleles are able to recombine away from maladaptive alleles. Following this, mean recombination rate for f d peaks versus the genome background was almost two times higher in Bekanbeushi (mean cM/Mb ± SD; peaks = 8.32 ± 19.5; background = 4.57 ± 9.36; permutation test, p < 1 × 10 -4 ; Figure S7a) and 1.5 times higher in Tokotan (peaks = 6.92 ± 14.0; background = 4.58 ± 9.49; p < 1 × 10 -4 ; Figure S7b). We additionally examined the distribution of f d peaks in the genome and found similar numbers of peaks per Mb on autosomes in both Bekanbeushi and Tokotan (Figures 5 and 6). However, in Tokotan, the highest proportion of f d peaks for the entire genome  (Table S3), and two were seen in both river systems. For Tokotan, four QTLs were located close to f d peaks, as opposed to three in Bekanbeushi. Among the five total QTLs, three showed heterozygote advantage. We additionally identified a total of 4,667 genes occurring within 100 kb of f d peaks in both comparisons with 3,411 seen in Bekanbeushi and 3,866 in Tokotan. GO analysis showed that genes in the close vicinity of introgression sites from both sites were enriched for GO terms related to olfactory receptor activity and detection of chemical stimuli (Table S4).

| D ISCUSS I ON
Our results clearly indicate that genomic divergence is lower and introgression rates are higher in the Tokotan system compared to the nearby Bekanbeushi system. In Bekanbeushi, the most extreme breeding habitats of the two species are separated by nearly 5 km, whereas in Tokotan, both species spawn in a 200 m stretch ( Figure 1).
Therefore, our results are consistent with the hypothesis that isolation of breeding habitats contributes to and maintains reproductive isolation between anadromous fishes. Although we observed that these two species are evenly distributed in Tokotan River, further detailed analysis of microhabitat use in Tokotan River is necessary to quantitatively measure the strength of habitat isolation in this river system. Additionally, it is possible that local adaptation to more divergent breeding habitats in Bekanbeushi might further contribute to reduced effective introgression in this river system relative to the Tokotan River system. Analysis of other nearby rivers containing the species pairs and differing in length would also help to test the role of habitat isolation in reproductive isolation of anadromous species.
Although introgression can be adaptive (Arnold, Sapir, & Martin, 2008;Hedrick, 2013;Racimo, Sankararaman, Nielsen, & Huerta-Sánchez, 2015), it is not always straightforward to conclusively show that this is the case. The integration of QTL analysis of tion. In the present study, we found that three QTLs with heterozygote advantage on juvenile survival are located near to introgression sites. This is consistent with the idea that loci with heterozygote advantage can introgress at higher rates than other genomic regions (Kim et al., 2018;Nadachowska-Brzyska et al., 2012;Wolfe et al., 2019), although we caution that the QTL peaks do not overlap with introgressed regions. We further found that the introgression sites we identified were enriched for genes related to olfactory reception and the detection of chemical stimuli. These genes were located on Chromosome I (257,080-372,834 bp) and Chromosome XVI (4,224,301,585 bp and 13,791,682,213 bp).
Interestingly, genes linked to olfactory receptor activity show signatures of overdominance in humans (Alonso, López, Izagirre, & de la Rúa, 2008). A plausible alternative explanation is that olfactory receptors are involved in conspecific mate recognition mechanism and are potentially acting as a one-allele barrier (Felsenstein, 1981;Garner, Goulet, Farnitano, Molina-Henao, & Hopkins, 2018;Smadja & Butlin, 2009). Positive selection for mate recognition where breeding grounds overlap would promote adaptive introgression in this case. Nonetheless, further work is now necessary to distinguish between these two alternative explanations. Furthermore, we found that many introgression sites were shared between both Bekanbeushi and Tokotan although the latter system has a larger number of unique regions. This suggests that introgression between the two species has likely occurred throughout their evolutionary history (Ravinet et al., 2018) but that it is also ongoing. In addition, shared regions of introgression and the presence of introgression close to genes and QTLs of functional significance are consistent with the hypothesis that some introgression in this system may be adaptive.
We additionally found that the strongest signatures of introgression occurred more often in regions with a higher recombination rate. Association between recombination rates and introgression has been reported in both theoretical (Butlin, 2005;Nachman & Payseur, 2012;Noor & Bennett, 2009) and empirical studies (Brandvain et al., 2014;Janoušek et al., 2015;Juric et al., 2016;Martin et al., 2019;Samuk et al., 2017;Schumer et al., 2018). In low recombination regions, selection purging deleterious foreign alleles is more likely to remove nearby physically linked loci, including both neutral and slightly advantageous alleles (Barton & Bengtsson, 1986;Ravinet et al., 2017). In high recombination regions, the probability that linked loci will recombine away from deleterious introgressed alleles is higher and thus signatures of introgression are more likely (Schumer et al., 2018). It should be noted, however, that our present analysis assumed that recombination rates were conserved among stickleback populations and species. The linkage map we used to determine the recombination rate was derived from a cross between sticklebacks from the Atlantic Ocean lineage (Roesti et al., 2013).
Thus, it is necessary to test whether this assumption is valid or not using a recombination rate map estimated using Japanese crosses in the future. Additionally, we also cannot exclude the possibility that regions with signatures of introgression result from introgression that occurred before the split between the Tokotan and Bekanbeushi populations or from ancestral polymorphism (Nelson & Cresko, 2018). The apparent overlap of 'geography' and 'interspecific' topologies in the TWISST analysis ( Figure 4)

rections. A further comparative analysis of introgression patterns is
necessary to distinguish this more clearly.

Although the rate of introgression is likely higher in Tokotan
River, the two species remain clearly genetically differentiated ( Figure 2) and have clearly not merged into a single species, in contrast with several well-known cases of speciation reversal (Seehausen et al., 2008;Taylor et al., 2006). This is likely because there are multiple isolating barriers present in the Japanese stickleback system, including hybrid male sterility, a post-zygotic incompatibility (Kitano et al., 2009). It is predicted that the presence of intrinsic incompatibilities and/or polygenic isolating barriers can lead to irreversible speciation (Barton & Bengtsson, 1986;Nosil et al., 2009;Seehausen, 2006) and our findings support this.
Sex chromosome turnover in the Japanese stickleback has previously been linked to reproductive isolation in this system (Kitano et al., 2009;Yoshida et al., 2014). Here our analysis shows that the neo-X chromosome did not show signs of introgression except at the distal end away from the fusion site ( Figure 5). This could partially be explained by a lack of recombination in males carrying the neo-Y despite the presence of recombination in females, reducing the effective recombination rates in this chromosomal region.
Additionally, the hemizygosity of the neo-X in males can expose recessively deleterious introgressed alleles to negative selection and efficiently purge them (Charlesworth, Coyne, & Barton, 1987;Johnson & Lachance, 2012). Furthermore, a previous QTL analysis showed that the neo-X chromosome contains QTL for the behavioural intensity of dorsal pricking and dorsal spine length, which is more sexually dimorphic in the Japan Sea stickleback and used for male courtship in this species (Kitano et al., 2009). Although these QTL are not located within the nonrecombining region of the neo-X, they are candidates for barrier loci that might act to reduce local introgression in the genome. In contrast to the neo-X, the ancestral-sex chromosome showed a larger number of introgression sites than the genomic background ( Figure 6a) and also a distribution of f d with a higher mean in Tokotan (Figure 6c). We are currently uncertain of what might cause this pattern, but higher introgression on sex chromosomes compared to autosomes has also been observed in other stickleback systems (Dixon, Kitano, & Kirkpatrick, 2019).
In conclusion, we found that the absence of habitat isolation substantially increases introgression and erodes genomic divergence between anadromous species. Additionally, we found that two species pairs with different levels of habitat isolation shared many introgression sites and some regions of introgression were located near QTL with heterozygote advantage on juvenile survival and genes that have previously been linked to heterozygote advantage in other species. Taken together, our findings suggest that at least some of the introgression observed may be adaptive. Despite the presence of introgression, these two species are still clearly differentiated. Other isolating mechanisms such as behavioural isolation and hybrid sterility likely contribute to reproductive isolation. Thus, the combination of multiple isolating barriers maintains reproductive isolation between the anadromous species, and the breakdown of one such barrier, here habitat isolation, is not enough to cause speciation reversal in this species pair. Overall, our data demonstrate that comparative genomic analysis of multiple species pairs with different levels of habitat isolation is informative for understanding how reproductive isolation and genomic divergence can be maintained by habitat isolation and other isolating barriers.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13664.