The genomic landscape at a late stage of stickleback speciation: High genomic divergence interspersed by small localized regions of introgression
[['Speciation is a continuous process through which reproductive isolation is established [1–3] .', 'According to the genic view of speciation [4] , when populations are in contact , gene flow is initially restricted at barrier loci ( i . e . loci underlying reproductive isolation ) , leading to the emergence of peaks of genetic differentiation surrounding such barriers; i . e . heterogeneous genomic differentiation [5 , 6] .', 'As speciation progresses , this localised build-up of reproductive isolation spreads to nearby regions due to linkage disequilibrium [4 , 5 , 7] .', 'Once a critical amount of differentiation at multiple barrier loci has accumulated , reduction of the genome-wide effective migration rate will eventually lead to divergence across the entire genome [5 , 7] .', 'This final step of genome-wide congealing may be a rapid and non-linear phase transition under certain conditions , such as when isolating barriers have a polygenic basis or a few strong barrier loci arise [8–10] .', 'Recent empirical genomic studies have revealed regions of high and low differentiation dispersed throughout the genome at early stages of speciation [7 , 11 , 12] .', 'This empirical data has lent strong support to the genic perspective of the speciation process [4] .', 'To-date however , the majority of speciation genomic studies demonstrating heterogeneous genetic differentiation have come from young species or population pairs with low divergence [7 , 11 , 12] .', 'Several thorough genomic studies on old sympatric species pairs exist , including European rabbits [13] , Drosophila species [14] , sunflowers [15] , whitefishes [16] , flycatchers [17 , 18] , wild mice [19] , Mimulus [20] and stick insects [9]; however except in a few cases , such as with Heliconius [21 , 22] , divergence is thought to have occurred during periods of geographical isolation .', 'Distinction between primary and secondary divergence is important for interpreting the patterns of genomic differentiation [12 , 17] .', 'This is because high genome-wide differentiation may have evolved via genetic drift and local adaptation during allopatric isolation , rather than due to divergence with gene flow .', 'Following secondary contact after geographical isolation , heterogeneous genomic differentiation may arise due to introgression .', 'Without a picture of the demographic history , this scenario may be indistinguishable from primary divergence [23] .', 'Despite the fact that the expected pattern of genomic differentiation during speciation is influenced by the timing and duration of geographical isolation [7] , testing different demographic histories has been somewhat neglected by the field [7 , 23] , although this is now changing [17 , 24] .', 'Other factors besides the demographic history of a species pair can also confound patterns of heterogeneous genomic differentiation .', 'For example , variation in recombination rate influences the patterns of genomic differentiation , because local adaptation or background selection in genomic regions where recombination is reduced can elevate differentiation measures and be mistaken for barrier loci [18 , 25 , 26] .', 'Mutation rate variation also influences the patterns of absolute divergence [27] .', 'Regions of low differentiation may be caused by shared ancestral polymorphism rather than gene flow [25 , 28] .', 'Distinction between gene flow and shared ancestral polymorphism is likely easier in more divergent species pairs [27 , 29 , 30] .', 'Furthermore , the use of multiple classical and recently developed methods , such as detection of recent hybrid progeny , ABBA-BABA tests [21 , 31] , model-based inference [32] , and comparisons between allopatric and sympatric pairs [21 , 26] provide a means to distinguish signatures of gene flow from alternative explanations .', 'It is therefore essential to account for factors such as demographic history , recombination rate variation , and shared ancestral polymorphism that can confound the interpretation of genome scan data [7 , 12] .', 'Three-spined stickleback species pairs ( genus Gasterosteus ) span the speciation continuum at varying stages of divergence , making them a model system for speciation research [33 , 34] .', 'To-date genomic research on speciation with gene flow in the stickleback complex has largely focused on weakly divergent species pairs , such as lake-stream ecotypes [35–37] .', 'Such studies have shown that the genomic landscape of differentiation between these recently diverged sympatric or parapatric species pairs is heterogeneous and interspersed with multiple peaks of high differentiation [35 , 37 , 38] .', 'The emerging pattern is consistent with predictions under the genic concept of speciation–i . e . that reproductive isolation is localized in the genome at early stages of divergence [4 , 39] .', 'However , it remains unclear whether such localized differentiation will eventually progress toward genome-wide differentiation in the face of gene flow [40] .', 'Toward the end of the stickleback speciation continuum is a marine species pair in Japan [41 , 42] .', 'The Japan Sea stickleback ( G . nipponicus ) is sympatric with the Pacific Ocean lineage of three-spined stickleback ( G . aculeatus ) ( Fig 1A ) in the waters surrounding the Japanese archipelago ( Fig 1C ) [41 , 43] .', 'Divergence time between the two marine species has been estimated to be 1 . 5–2 million years based on allozyme and microsatellite data [42 , 44] , making it much older than postglacial stickleback species pairs .', 'Divergence between the species may have occurred as a result of the repeated isolation of the Sea of Japan during the Pleistocene , but this divergence scenario remains to be explicitly tested [42 , 44] .', 'A unique feature of the G . nipponicus and G . aculeatus system , relative to postglacial stickleback species pairs , is that a neo-sex chromosome has arisen due to a fusion between a Y chromosome and a previously autosomal chromosome IX ( chrIX ) in the G . nipponicus lineage [41 , 45] .', 'Furthermore , crosses between Japan Sea females and Pacific Ocean males show hybrid male sterility [42] .', 'Previous quantitative trait locus ( QTL ) mapping identified QTL for courtship behaviour on the neo-X and hybrid male sterility on the ancestral-X .', 'However , there are other isolating barriers , such as eco-geographical isolation , temporal isolation , and ecological selection against migrants [42 , 46 , 47] .', 'The combination of these multiple barriers most likely contributes to the strong reproductive isolation in this system [41 , 48] .', 'However , despite such strong divergence , hybrids have been observed where the two species co-occur in Northern Japan [41] and phylogenetic discordance between nuclear and mitochondrial loci suggests some history of introgression during speciation [49 , 50] .', 'Although the Japanese species pair represents one of the furthest points of divergence within the stickleback species complex , speciation remains incomplete .', 'The evolutionary history and genome-wide patterns of genetic differentiation and introgression of this strongly divergent species pair therefore remains an open question .', 'The aim of our study was to address this gap in our knowledge; i . e . to quantify the patterns of genomic differentiation and introgression at a later stage of the stickleback speciation continuum .', 'To this end , we used previously published whole-genome sequences and newly acquired Restriction-site Associated DNA sequencing ( RAD-seq ) data from the Japanese stickleback species pair to determine their evolutionary history and characterise patterns of gene flow between them .', 'Our first aim was to establish how and when divergence took place between G . nipponicus and G . aculeatus .', 'Using thousands of genomic loci and a coalescent modelling approach on the resequence data , we tested a range of divergence scenarios and estimated the timing and duration of isolation , the extent of gene flow and fluctuations in population size .', 'After identifying that the two species have indeed diverged in the face of gene flow , we first used our RAD-seq dataset to investigate patterns of population structure and introgression between the Japanese stickleback species pairs .', 'We then used a comparative genome scan approach with the resequence data , adding G . aculeatus lineage from the Atlantic Ocean [51] as an allopatric control ( Fig 1A , S1 Fig ) .', 'After establishing that gene flow has occurred but that a high level of genomic differentiation has remained , we used two independent measures of gene flow to identify where in the genome introgression has left its mark .', 'We tested whether introgression occurs more frequently in regions of high recombination and whether it occurs in regions with functionally important genes .', 'Our findings suggest a high level of genome-wide divergence can be maintained in the face of gene flow , as introgression is restricted to small , localized genomic regions .'], ['Phylogenetic analysis on 35 , 666 10 kb non-overlapping genome windows on autosomes ( i . e . , excluding chrIX and chrXIX ) using whole genome resequence data on 26 individuals supports a deep split between G . aculeatus ( both Pacific and Atlantic Ocean lineages ) and G . nipponicus ( Japan Sea stickleback ) ( Fig 1A ) .', 'Of all windows , 98 . 8% support the split between species , while only 0 . 51% indicate clustering of fish occurring in Japan ( the Japanese Pacific Ocean G . aculeatus and the Japan Sea G . nipponicus; S1 Table and Fig 1A ) .', 'We calculated genealogical sorting index ( gsi ) [52] on maximum likelihood phylogenies estimated from non-overlapping sliding windows of 10 kb across the autosomes .', 'High gsi indicates monophyly , while low gsi indicates mixed ancestry [52] .', 'Genome-wide averages ( ± SD ) of gsi were high , but not complete , for all three Gasterosteus lineages with that of the Japan Sea stickleback being the highest ( Atlantic gsi = 0 . 45 ± 0 . 10 , Pacific gsi = 0 . 57 ± 0 . 09 , Japan Sea gsi = 0 . 72 ± 0 . 06 ) .', 'This is in stark contrast to the mitogenome phylogeny where sticklebacks from both species occurring in Japan fall into a single clade separate from the clade occurring in the Western Pacific and Atlantic ( Fig 1B , S2 Fig ) .', 'A lack of mitogenome divergence between G . aculeatus and G . nipponicus from the Japanese archipelago suggests mitochondrial introgression might occur where these lineages overlap ( Fig 1C ) .', 'Since the consensus autosomal phylogeny suggests a more recent split between the Japanese Pacific and Atlantic G . aculeatus lineages than the split in the mitochondrial phylogeny , the two mitogenome clades may represent the split between G . aculeatus and G . nipponicus lineages with mitochondrial introgression likely having occurred from the Japan Sea G . nipponicus into the Pacific Ocean G . aculeatus in sympatry .', 'Divergence time estimates between the mitogenome clades are thus informative for dating the divergence time between G . aculeatus and G . nipponicus lineages .', 'Bayesian coalescent analysis using a strict clock model in Bayesian Evolutionary Analysis by Sampling Trees ( BEAST ) suggests a median split date of 1 . 30 million years ( 0 . 15–2 . 41; 95% Highest Posterior Density [HPD] intervals; S2 Table ) for the two major mitogenome clades ( S2 Fig ) , consistent with previous estimates [49] .', 'Divergence between Eastern Pacific and Atlantic haplotypes is more recent at 0 . 39 million years ( 0 . 03–0 . 74; 95% HPD ) but is older than the Most Recent Common Ancestor ( MRCA ) of all haplotypes occurring in Japan ( Fig 1B , S2 Fig ) , suggesting mitochondrial gene flow from G . nipponicus to G . aculeatus may have occurred within the last 0 . 39 million years .', 'To investigate the demographic history of G . aculeatus and G . nipponicus , we first used pairwise sequential Markov coalescent ( PSMC ) on all 26 Atlantic Ocean , Japan Sea and Pacific Ocean resequenced stickleback genomes to examine fluctuations in effective population size .', 'Strikingly , G . nipponicus experienced a severe bottleneck around 0 . 15–0 . 3 million years before present ( BP ) ( Fig 1D ) ; mean Ne fell to 26 , 422 ± 1 , 191 at its lowest point .', 'Subsequently after 0 . 1 million years BP , G . nipponicus underwent a dramatic effective population size expansion ( Fig 1D ) : mean Ne rose to 195 , 974 ± 28 , 832 ( i . e . ~7 . 5 times increase from the bottleneck ) during the late Pleistocene .', 'In contrast , the effective population size of the Japanese Pacific Ocean G . aculeatus has remained relatively stable throughout its history ( mean Ne ± SD = 118 , 150 ± 4 , 330; Fig 1D , see S3 Fig for bootstrap support ) .', 'Although the Atlantic ( Fig 1D ) and Western Pacific lineages of G . aculeatus ( S4 Fig ) also experienced some growth during the late Pleistocene , their effective population sizes remained smaller than that of G . nipponicus .', 'Cryptic population structure in G . nipponicus might explain the disparity in Ne between lineages; however our RAD-sequence dataset confirms substructure is not present in this species ( see below for more details on RAD-seq dataset; S5 Fig and S6 Fig ) .', 'Furthermore , genome-wide averages of Tajima’s D also support a recent demographic expansion for G . nipponicus ( mean ± SD of Tajima’s D = -0 . 82±0 . 45 ) and stable effective population size in the Pacific Ocean ( mean ± SD of Tajima’s D = -0 . 04 ± 0 . 63 ) .', 'To explicitly test whether divergence between G . aculeatus and G . nipponicus occurred in the presence of gene flow , we used an Approximate Bayesian Computation ( ABC ) approach with 1 , 874 2 kb loci randomly sampled from across autosomes .', 'We tested five divergence scenarios–isolation ( I ) , isolation with migration ( IM ) , isolation-with-ancient-migration ( IAM ) , isolation-with-recent-migration ( IRM ) and isolation-with-ancient-and-recent-migration ( IARM ) –i . e . two discrete periods of contact .', 'Since the results of our PSMC analyses indicate Ne has varied throughout divergence ( Fig 1D ) , we performed a hierarchical ABC analysis , first selecting the most appropriate population growth model ( i . e . constant size , population growth and a Japan Sea bottleneck ) within each divergence scenario and then performing final model selection amongst the best supported divergence/growth model scenarios ( see S1 Text for full specification of models , priors , parameters and extensive sensitivity testing ) .', 'Using 20 summary statistics ( see S1 Text for a full list of statistics used ) and a neural-network rejection method with 1% tolerance of simulated datasets , the best-supported divergence scenario was a model of IM with a bottleneck occurring only in the Japan Sea species ( Fig 2A , Table 1 ) .', 'An IARM model was the second best supported model .', 'The use of a standard ABC rejection method gave rise to the qualitatively similar results and we found no evidence of an overrepresentation of introgressed regions in the loci used as the observed data for this analysis ( S1 Text ) .', 'An independent maximum likelihood based demographic analysis using the joint G . aculeatus and G . nipponicus site frequency spectrum ( SFS ) derived from RAD-seq data showed high support for an IARM model ( see S1 Text ) .', 'Parameter estimates from the ABC IM model suggest divergence between G . aculeatus and G . nipponicus occurred 0 . 68 million years ago ( median estimate , 0 . 18–4 . 17 million years , lower & upper 95% HPD; Fig 2B ) .', 'A Japan Sea bottleneck occurred 0 . 3 million years ago ( 0 . 03–2 . 21 million years 95% HPD ) , reducing Ne to about 20% of the contemporary estimate ( Fig 2C , S3 Table ) .', 'Mean migration rates between the two species were low , and migration rate ( expressed as mij−i . e proportion of population i that are migrants from j per generation ) from the Pacific Ocean lineage into the Japan Sea lineage ( m12: median = 1 . 3 x 10−6 , 95% HPD = 8 . 61 x 10−8–5 . 32 x 10−6 ) was slightly greater than in the opposite direction ( m21: median = 1 . 05 x 10−6 , 95% HPD = 4 . 91 x 10−8–6 . 39 x 10−6 , N . B . migration rates are backwards in time; see also Fig 2D & S3 Table ) .', 'In addition to this , the distribution of the migration rate hyperprior suggested that a large number of loci showed some level of gene flow ( S1 Text ) .', 'Contemporary Ne of the Japan Sea lineage is larger than that of the Pacific Ocean , although the Ne estimates differed in magnitude from those estimated by PSMC ( Figs 1D and 2C , S3 Table ) .', 'Given this difference in effective population size , the scaled migration rates , the expected number of migrants per generation , ( 2Nimij ) are higher from Pacific Ocean lineage into the Japan Sea than the alternative ( PO to JS = 0 . 18; JS to PO = 0 . 04 ) in contemporary populations , although still very low .', 'Scaled migration rates were likely more similar during the Japan Sea bottleneck , because lower effective population size of the Japan Sea population ( 1 . 22 x 104 ) at this stage reduces the expected number of migrants from the Pacific Ocean to the Japan Sea ( 0 . 031 ) .', 'Identifying admixture and the presence of backcrossed individuals between species where they co-occur provides strong evidence of on-going introgression [7 , 12] .', 'To address this , we used a RAD-sequencing dataset with a larger sample size of 245 individuals from the Atlantic , Pacific and Japan Sea lineages , including previously published data from Pacific-derived populations in North America [53] .', 'Principal component analysis ( PCA ) of allele frequencies at 3 , 744 high-quality bi-allelic SNPs pruned to remove loci in linkage disequilibrium showed that , consistent with our whole genome data , the main axis explaining 20% of the variance was between G . aculeatus and G . nipponicus ( S5 Fig ) .', 'The secondary axis explaining 9 . 49% of the variance was mainly between the Atlantic and Pacific populations ( S5 Fig ) .', 'Importantly , PCA showed a single individual was intermediate between the Pacific and Japan Sea populations occurring in Akkeshi , the sympatric site in Hokkaido , Japan where our whole genome-sequenced samples were collected ( Fig 1C ) .', 'A separate Bayesian analysis for admixture using STRUCTURE [54 , 55] found greatest support for K = 2 among stickleback populations and also identified the putative F1 hybrid plus individuals with possible recent admixture in Akkeshi ( S6 Fig ) .', 'To further investigate variation in individual ancestry , we identified 5 , 967 ancestry-informative loci i . e . autosomal SNPs with an allele frequency difference of >0 . 8 between the Japan Sea and Pacific Ocean lineages .', 'Using a genomic cline approach , we estimated interspecific heterozygosity ( i . e . , proportion of loci with alleles from both species ) and hybrid index ( i . e . , proportion of alleles from one species ) on simulated hybrid genotypes .', 'This indicates the marker set has high power to detect hybrid ancestry ( S7 Fig ) .', 'Analyses on the observed data suggest the RAD-seq dataset includes one F1 hybrid and several individuals with likely hybrid ancestry in the last few generations ( S7 Fig ) .', 'Taken together , these data indicate that divergence between the Japanese G . aculeatus and G . nipponicus is much older and greater compared to commonly studied postglacial stickleback species pairs .', 'Despite the great extent of divergence between Japanese stickleback species , parameter estimates and observational data suggest that gene flow between them is on-going .', 'Genome-wide differentiation was strikingly high between G . nipponicus and G . aculeatus regardless of their geographical overlap ( Fig 3A & 3B and Fig 4 , and S8 Fig and S9 Fig ) .', 'The genome-wide average of FST between the sympatric species was 0 . 628; this is higher than FST in all other studied stickleback species pairs , which is typically less than 0 . 3 [35–37 , 56] ( see Fig 3C ) .', 'The genome-wide average of absolute divergence ( dXY ) was 0 . 012; which is also high compared to previously calculated dXY values , i . e . less than 0 . 005 , between postglacial parapatric and sympatric stickleback ecotypes [35 , 57 , 58] .', 'Despite consistently high divergence , both FST and dXY values were significantly lower where the two species occur in contact ( Table 2 , Figs 3 & , 4 , S8 Fig and S9 Fig; 10 , 000 replicate permutation tests on 10 kb windows: P < 2 . 2 x 10−16 for both statistics ) , consistent with the presence of gene flow in sympatry .', 'A more fine-scale analysis of genome-wide divergence based on 10 kb non-overlapping windows revealed that the high baseline divergence between G . nipponicus and G . aculeatus is interspersed by regions of low differentiation in both FST and dXY genome scans ( Fig 4 top two panels , S8 Fig and S9 Fig ) , possibly indicating introgression .', 'To identify genomic regions of recent introgression , we calculated two independent measures .', 'The first of these was GMIN , the ratio of the minimum dXY to the average dXY [30] .', 'Under strict isolation , minimum dXY relates to the upper bound of divergence time between two populations , whereas when introgression occurs , minimum dXY reflects the timing of the most recent migration event [30] .', 'The second measure was fd , an estimate of the proportion of introgressed sites in a genome window , calculated using a four population ABBA-BABA test [59] .', 'GMIN is more effective at identifying recent , low level gene flow than either FST or dXY but by definition it is unable to detect genomic regions where complete introgression has occurred [30] , which can however be detected using fd .', 'Importantly , both measures are robust to variation in recombination rate [30 , 59] .', 'Combining these two statistics therefore allows us to identify both low-level ( GMIN ) and strong introgression ( fd ) .', 'Focusing on between species comparisons , mean ( ± SD ) GMIN measured from 10 kb non-overlapping windows was greater in allopatry than sympatry ( Japan Sea vs . Atlantic: 0 . 876 ± 0 . 071; Japan Sea vs . Pacific: 0 . 857±0 . 103; randomization test P < 2 . 2x10-16; Fig 4 ) .', 'Mean fd was also greater when the species overlapped ( JS vs . AT: -0 . 0031 ± 0 . 0540; JS vs . PO: 0 . 0039±0 . 0328; P < 2 . 2x10-16; Fig 4 ) , and both statistics are more strongly negatively correlated in sympatry ( S10 Fig ) supporting gene flow between G . nipponicus and Japanese populations of G . aculeatus .', 'Genomic regions of low GMIN ( i . e . GMIN valleys ) may indicate recent introgression .', 'We identified genome windows with low GMIN values using a Hidden-Markov classification model [60] ( S11 Fig ) .', 'We then clustered 10 kb outlier windows occurring within 30 kb of one another into putative GMIN valleys .', 'GMIN in particular may be susceptible to false positives as a result of shared ancestral polymorphism .', 'However , lower dXY and higher fd in sympatric GMIN valley windows compared to the genomic background suggests shared ancestral polymorphism alone does not explain the patterns observed here ( S12 Fig; randomization test , P < 2 . 2 x 10−16 in both cases ) .', 'These regions of introgression were more common in the genome when the two species overlapped , with 637 valleys in sympatry ( JS-PO comparison ) compared to 337 in allopatry ( JS-AT comparison ) ( randomization test , t = 5 . 35 , P < 2 . 2 x 10−16 ) and a greater number of valleys per chromosome ( Fig 5A ) , although mean valley size did not differ significantly ( 77 . 6 kb and 75 . 4 kb in sympatry and allopatry respectively , P = 0 . 82 ) .', 'Interestingly , 225 valleys were shared between JS-PO and JS-AT comparisons ( Fig 4 ) .', 'These shared valleys may indicate shared ancestral polymorphism but they may also reflect introgression from Pacific Ocean to Japan Sea , where one or a few Japan Sea individuals carry haplotypes derived from Pacific Ocean and therefore are also similar to Atlantic Ocean haplotypes too .', 'However , a larger number of valleys ( 412 valleys ) were unique to the JS-PO comparison , where introgression might occur from Japan Sea to Pacific Ocean .', 'A similar geographical comparison of peaks of fd between species was not possible , due to the fact that fd is much closer to 0 in the comparison between G . nipponicus and the Atlantic G . aculeatus and very few peaks were present ( Fig 4 ) .', 'Nonetheless , Hidden-Markov classification identified 823 fd peaks occurring between G . nipponicus and Pacific G . aculeatus ( S13 Fig ) .', 'If the fd peaks mainly indicate introgression from Pacific Ocean to Japan Sea , dXY between Japan Sea and Atlantic Ocean is expected to be lower in these regions compared to the genome background , as Japan Sea fish carry haplotypes derived from the Pacific Ocean , which in turn are similar to the Atlantic Ocean haplotypes .', 'While JS-AT dXY was lower in fd peaks compared to the genome background ( JS-AT mean dXY ± SD , fd peaks: 0 . 0121±0 . 0035 , genome-background: 0 . 0127±0 . 0026; one-tailed permutation test , P < 2 . 2 x 10−16 ) , this difference was not very clear ( S14 Fig ) .', 'In contrast , if introgression occurred mainly from Japan Sea to Pacific Ocean , dXY in the PO-AT comparison should increase in fd peaks relative to the genome background , as Pacific Ocean fish carry Japan Sea-derived haplotypes , which are divergent from the Atlantic Ocean haplotypes .', 'We clearly observed this pattern ( PO-AT mean dXY ± SD , fd peaks: 0 . 0065±0 . 0035 , genome-background: 0 . 0038±0 . 00182; P < 2 . 2 x 10−16; S14 Fig ) ; suggesting that introgression from Japan Sea to Pacific Ocean may be more predominant than the opposite direction .', 'Importantly , our findings using GMIN , dXY and fd are robust to different missing data thresholds and did not change when phased vs . unphased data is used ( S1 Text ) .', 'To further investigate the direction of gene flow , we used partitioned D statistics ( an extension of the four population test–see S15 Fig ) , which tests the excess of shared derived alleles using five , rather than four populations [61] .', 'To this end , we added an allopatric Japan Sea population ( collected from Lake Shinji , a brackish lake at the Japan Sea coast of southern Honshu ) .', 'A positive D12 statistic is proposed to indicate the predominance of introgression from P3 to P2 ( S15 Fig ) [61] .', 'When P3 was set to Japan Sea ( where P31 is sympatric and P32 is allopatric with the Pacific Ocean ) and P2 to Pacific Ocean ( see S14 Fig ) , D12 was significantly positive in fd peaks ( one-tailed permutation test , P < 2 . 2 x 10−16 ) .', 'In contrast , when we rotated the populations at the tips–i . e . setting P2 to sympatric Japan Sea , P31 to Pacific Ocean , and P32 to Atlantic Ocean ( see S15 Fig ) , D12 was not positive , consistent with the suggestion that introgression is occurring mainly from Japan Sea to Pacific Ocean .', 'However , the resolution of partitioned D statistics has been criticized [62]; positive D12 can also be caused by introgression from the Pacific Ocean ( P2 ) to the common ancestor of the sympatric and allopatric Japan Sea populations ( P31 & P32 ) .', 'To overcome this issue , we calculated DFOIL , which also uses a five-population test but accounts for all possible introgression events [62] .', 'When P1 = sympatric Japan Sea , P2 = allopatric Japan Sea , P3 = Pacific Ocean , and P4 = Atlantic Ocean ( S16 Fig ) , DFOIL clearly indicated the presence of ancestral introgression ( 239 out of 4 , 236 100 kb-windows ) between the Japan Sea ancestor ( P12 ) and the Pacific Ocean ( P3 ) ( see S15 Fig ) .', 'However , we found only a few windows showing unidirectional introgression ( 6 in total ) , meaning we cannot determine the direction of introgression using this analysis ( S16 Fig ) .', 'This low sensitivity may be due to the fact that structuring in the Japan Sea lineage is very low ( S6 Fig ) [63]–i . e . recent divergence time between the sympatric and allopatric Japan Sea populations or high intraspecific gene flow within the Japan Sea species .', 'To investigate whether introgression co-varies with recombination rate , we used a previously published recombination map from an Atlantic G . aculeatus cross [64] to interpolate genome-wide recombination rate variation ( see Methods ) .', 'We detected a negative correlation between recombination rate and GMIN and a positive correlation with fd ( Pearson’s correlation , GMIN: r = -0 . 17 , P < 2 . 2 x 10−16; fd: r = 0 . 08 , P < 2 . 2 x 10−16 , S17 Fig ) .', 'Accordingly , mean recombination rate for putatively introgressed regions was over two times higher than the genome background ( GMIN: valley = 8 . 98 cM/Mb , non-valley = 3 . 99 cM/Mb; fd: peak = 9 . 64 cM/Mb , non-peak = 4 . 16 cM/Mb; randomization test P < 2 . 2 x 10−16 in both cases; Fig 5B ) .', 'Sex chromosomes likely played an important role in speciation between G . aculeatus and G . nipponicus [41 , 45] .', 'A fusion between Y and chrIX means that chrIX segregates as a neo-sex chromosome in G . nipponicus but not G . aculeatus which only carries the ancestral and shared sex chromosome , chrXIX [41 , 45] .', 'The divergent XY ( G . aculeatus ) and X1X2Y ( G . nipponicus ) systems means that recombination is reduced for chrIX and chrXIX in hybrids carrying the neo-Y [45] .', 'Given this recombination rate reduction and previously identified QTL for traits involved in reproductive isolation that map to chrIX and chrXIX [41 , 45] , we tested whether recent introgression ( i . e . measured using GMIN ) was reduced in this part of the genome relative to the autosome .', 'For this , we repeated our analyses using females only ( 5 Japan Sea and 6 Pacific Ocean ) .', 'The number and density of valleys was lowest on the neo-sex chromosome , chrIX ( 16 valleys or 0 . 8 valleys per Mb ) but not on the ancestral sex chromosome ( chrXIX , see S4 Table ) .', 'Finally , we investigated the nature of introgression between the two species .', 'We first asked whether introgression occurs more frequently in genic or non-genic regions .', 'We identified 3 , 261 genes occurring in GMIN valleys and 2 , 958 genes from fd peaks between sympatric G . aculeatus and G . nipponicus; 60% of genes identified were found in both types of introgressed window , whereas 23% occurred only in GMIN valleys and 15% only in fd peaks ( S18 Fig ) .', 'Irrespective of the method used to detect putatively introgressed regions , the number of genes identified was greater than the number expected by chance ( P < 0 . 0001 based on a null distribution generated from 1 , 000 random samples of the genome ) .', 'Mean recombination rate was higher in the genomic windows where genes are present compared to the genomic background ( gene windows = 4 . 92 cM/Mb , genome-background = 4 . 24 cM/Mb; permutation test: P < 2 . 2 x 10−16 ) .', 'This suggests that introgression may be more likely in genic regions of the genome than non-genic regions , which can be partly explained by higher recombination rates in genic regions .', 'To further investigate the functional enrichment of the genes occurring in regions of introgression , we performed gene ontology ( GO ) analysis on 2 , 310 GMIN valley and 2 , 217 fd peak genes with orthologs in the human genome .', 'Enriched GO terms for fd peaks included immune response , metabolic processes and chromatin assembly , while enriched GO terms for GMIN valleys included major histocompatibility complex ( MHC ) protein and metabolic processes ( S5 Table & S6 Table ) .'], ['Determining the demographic and evolutionary history of species pairs is an important first step for understanding how speciation has unfolded in any system [7 , 12] .', 'Our present study has produced several lines of evidence indicating that divergence between the Japanese sticklebacks has occurred in the presence of gene flow .', 'Firstly , our ABC analysis supported a model of isolation with migration .', 'Previously , it has been speculated that the Japan Sea stickleback diverged largely as a result of geographical isolation in the Sea of Japan caused by sea level fluctuation during the early Pleistocene [42 , 44] .', 'Using ABC , we were able to explicitly test several divergence hypotheses in a statistical framework [65]; our findings suggest that gene flow has likely occurred throughout majority of the divergence history .', 'It should be noted that ABC and most established demographic inference methods perf
['Introduction', 'Results', 'Discussion', 'Materials and methods']
['Speciation is a continuous process and analysis of species pairs at different stages of divergence provides insight into how it unfolds .', 'Previous genomic studies on young species pairs have revealed peaks of divergence and heterogeneous genomic differentiation .', 'Yet less known is how localised peaks of differentiation progress to genome-wide divergence during the later stages of speciation in the presence of persistent gene flow .', 'Spanning the speciation continuum , stickleback species pairs are ideal for investigating how genomic divergence builds up during speciation .', 'However , attention has largely focused on young postglacial species pairs , with little knowledge of the genomic signatures of divergence and introgression in older stickleback systems .', 'The Japanese stickleback species pair , composed of the Pacific Ocean three-spined stickleback ( Gasterosteus aculeatus ) and the Japan Sea stickleback ( G . nipponicus ) , which co-occur in the Japanese islands , is at a late stage of speciation .', 'Divergence likely started well before the end of the last glacial period and crosses between Japan Sea females and Pacific Ocean males result in hybrid male sterility .', 'Here we use coalescent analyses and Approximate Bayesian Computation to show that the two species split approximately 0 . 68–1 million years ago but that they have continued to exchange genes at a low rate throughout divergence .', 'Population genomic data revealed that , despite gene flow , a high level of genomic differentiation is maintained across the majority of the genome .', 'However , we identified multiple , small regions of introgression , occurring mainly in areas of low recombination rate .', 'Our results demonstrate that a high level of genome-wide divergence can establish in the face of persistent introgression and that gene flow can be localized to small genomic regions at the later stages of speciation with gene flow .']
['When species evolve , reproductive isolation leads to a build-up of differentiation in the genome where genes involved in the process occur .', 'Spanning the speciation continuum , stickleback species pairs are ideal for investigating how genomic divergence accumulates during speciation .', 'However , much of our understanding of stickleback speciation comes from early stage divergence , with relatively few examples from more divergent species pairs that still exchange genes .', 'To address this , we focused on Pacific Ocean and Japan Sea sticklebacks , which co-occur in the Japanese islands .', 'We established that they are the oldest and most divergent known stickleback species pair , that they evolved in the face of gene flow and that this gene flow is still on going .', 'We found introgression is confined to small , localised genomic regions where recombination rate is high .', 'Our results show high divergence can be maintained between species , despite extensive gene flow .']
['taxonomy', 'fish', 'atlantic ocean', 'japan', 'population genetics', 'geographical locations', 'oceans', 'vertebrates', 'animals', 'osteichthyes', 'phylogenetics', 'data management', 'pacific ocean', 'phylogenetic analysis', 'speciation', 'bodies of water', 'population biology', 'computer and information sciences', 'marine and aquatic sciences', 'evolutionary systematics', 'sticklebacks', 'people and places', 'eukaryota', 'asia', 'heredity', 'earth sciences', 'genetics', 'biology and life sciences', 'gene flow', 'introgression', 'evolutionary biology', 'evolutionary processes', 'organisms']