Introduction

Geographic isolation is the most commonly accepted cause of genetic divergence between populations, and can sometimes lead to new species with distinct phenotypes (Coyne and Orr 2004). In contrast, if reproductive isolation (pre- or post-zygotic isolation) is not completed between populations, hybridization will occur in a secondary zone. The resulting hybrids and gene introgression are subject to natural selective pressures, including endogenous (genetic incompatibility) and exogenous selection (natural selection). Therefore, many genetic studies have been performed on secondary contact zones to investigate their evolutionary outcomes, including genomic and phenotypic features of hybrids (Jiggins and Mallet 2000; Nolte et al. 2005; Hedrick 2013). The outcomes of the resulting interplay between allopatric populations have improved our understanding of the process through which geographic isolation promotes biodiversity (April et al. 2013).

Several studies have shown that genomic features of independent secondary contact zones are variable (Nolte et al. 2009; Aboim et al. 2010), suggesting that the evolutionary outcomes of secondary contact are influenced by various exogenous factors, including local environmental features and past formative history. In this regard, investigations of multiple independent secondary contact zones will provide the answers about the inevitability of evolutionary consequences (Morgan-Richards and Wallis 2003; Mandeville et al. 2017). However, such studies on species inhabiting marine environments are rare. Given the environments spatially varying with oceanographic conditions, the impact of hybridization through secondary contact may play a substantial role in shaping biodiversity in marine environments.

Oscillatory changes in marine environments around the Japanese archipelago were remarkable during the Pleistocene, particularly in the Japan Sea, which is located between the Japanese archipelago and the Asian continent. The Japan Sea was subject to multiple isolation events during glacial periods of the Pleistocene (Tada 1994) and several studies suggest that these isolation events generated contrasting geographic lineages in the Pacific Ocean and the Japan Sea. Accordingly, deep genetic divergences have been identified in several coastal marine species (Kojima et al. 2004; Akihito et al. 2008; Hirase et al. 2012; Hirase and Ikeda 2014a, 2014b) and in anadromous species (Higuchi and Goto 1996; Kokita and Nohara 2011). During Pleistocene glacial periods, environments in the Japan Sea were clearly colder and had lower salinity due to closure of the south Tsushima Strait, with cold fresh water from rivers supplying much of the water. Therefore, the Japan Sea and Pacific Ocean lineages in the Japanese coastal species were subjected to vastly differing selection regimes, offering a useful model system for studies of phenotypic and genomic divergence (Kokita et al. 2013).

There are two potential secondary contact zones between the Japan Sea and Pacific Ocean lineages of coastal marine species: western and eastern coastal areas of Japanese archipelago where these seas are connected. A previous study showed such independent secondary contact zones in the ice goby Leucopsarion petersii (Kokita and Nohara 2011), which is a small paedomorphic bony fish with larval-like forms that include slender translucent and scaleless bodies (Matsui 1986). L. petersii is a temperate-water anadromous species that is distributed around the Japanese archipelago as Japan Sea and Pacific Ocean lineages (JS and PO lineages hereafter), which exhibit substantial genetic divergences that likely date back at least one million years (Kokita and Nohara 2011). Under natural conditions, individuals of the JS lineage have larger body sizes and numbers of vertebra than those of the PO lineage (Fig. 1a; Kokita and Nohara 2011). Given that differences in body sizes have been reproduced in common garden experiments, this phenotypic difference has a genetic basis (Kokita et al. 2017), and a neuropeptide Y (NPY) gene was identified as the candidate gene that contributes to this difference by a limited genome-scan approach (Kokita et al. 2013). Moreover, divergence in female mate choice patterns was detected; the JS lineage females demonstrated preferences for larger and same-lineage males, whereas the PO lineage females did not (Kokita et al. 2017). Nonetheless, reproductive isolation between the two lineages has not been completed, and sharing patterns of divergent haplotypes of mitochondrial (mt) DNA cytochrome b (cyt b) and the nuclear myosin heavy chain 6 (MYH6) gene suggested hybridization in the Seto Island Sea and the Joban-Kashimanada, which are environmentally distinct because of differing ocean currents (Fig. 1b; Kokita and Nohara 2011). Here we obtained genotype and phenotype data from individuals in the western and eastern secondary contact zones and compared these patterns of genomic admixture and phenotypic features as a consequence of secondary contact in different environmental settings.

Fig. 1
figure 1

a Photograph (scale bar: 10 mm) of mature females of the Japan Sea (JS) and Pacific Ocean (PO) lineages of Leucopsarion petersii (reproduced with permission from Kokita et al. 2017). b Sampling locations for this study in the Japanese archipelago. The JS and PO lineages are mainly distributed in each sea region except that JS lineage has spread to the northern part of the Pacific Ocean side (Kokita and Nohara 2011). Closed and open circles indicate sampling locations of pure JS and PO lineages, respectively. Stars show introgressed populations in Seto Inland Sea and Joban-Kashimanada. These symbols were used in the following figures. Tsuruga population of the JS lineage used for the NPY analysis was shown by “z.” c. Pie graphs show proportions of mitochondria genome of JS (black) and PO (white) lineages. ADMIXTURE plots showing individual genotype membership in K clusters based on RAD sequencing. Assuming K = 2, the cluster of PO lineage (black) and JS lineage (white) were estimated. If we assume K = 3, Joban-Kashimanada populations were assigned to the third cluster (gray). Hybrid index was calculated based on divergent SNP loci (hybrid index becomes 0 if all SNP loci have homozygous PO alleles). Horizontal gray lines indicate 0.1 and 0.9 thresholds

Full size image

Materials and methods

Sampling and DNA extraction

Leucopsarion petersii are anadromous and ascend to the lower reaches of rivers to reproduce. We collected 11 L. petersii populations along Japanese coastlines (Fig. 1b and Table 1), and three of these were of the PO lineage, three were of the JS lineage, and five were populations in secondary contact zones where admixtures were detected using mitochondrial and nuclear DNA markers (hereafter we will refer to these five populations as introgressed populations; Kokita and Nohara 2011). Some of the specimens were previously used to sequence mtDNA cyt b and genotype microsatellite markers (Kokita and Nohara 2011; Kokita et al. 2013). DNA was extracted from muscle or fin tissues using phenol/chloroform extraction as described previously (Asahida et al. 1996).

Table 1 Genetic diversity of Leucopsarion petersii populations analysed in this study
Full size table

Mitochondrial DNA (mtDNA) analysis

We analyzed mtDNA of individuals in the five introgressed populations by direct sequencing of an 897-bp region at the 5′-end of cyt b, as described by Kokita and Nohara (2011). Multiple sequence alignments with published cyt b haplotype sequences were performed using the FFT-NS-2 algorithm of MAFFT ver. 7.058 (Katoh and Standley 2013) with manual checking. The sequences were then assigned to haplotype sequences using the DNAcollapser of FaBox (Villesen 2007). Haplotype sequences were submitted to the DNA Data Bank of Japan (DDBJ; LC465830−LC465988; http://www.ddbj.nig.ac.jp). Phylogenetic analyses were conducted using the maximum-likelihood method in MEGA 5.2.2 (Tamura et al. 2011) with 1000 bootstrap replicates in the GTR + I + G substitution model, which was also used in a previous study (Kokita and Nohara 2011).

RAD sequencing, bioinformatics, and genotyping analyses

We identified genome-wide single nucleotide polymorphisms (SNPs) by aligning RAD sequencing reads to L. petersii draft genome contigs, which we roughly constructed. For the draft genome assembly, one JS individual sampled in 2017 from Saga Prefecture, Japan, was chosen. Whole-genome shotgun short-read sequences were generated using a paired-end (100 bp) library with an HiSeq X Ten (Illumina). Library construction and sequencing were performed at Macrogen Inc. We generated 143-Gb raw genome sequence data, representing a 130-fold depth of coverage of the predicted genome size that was estimated to be 1.1 Gb by Kmergenie ver.1.7 (Chikhi and Medvedev 2013). Before assembly, paired-end reads were filtered using platanus_trim with a default parameter included in Platanus ver.1.2.4 (Kajitani et al. 2014). Platanus with a default parameter was used for constructing contigs with filtered paired-end reads, and finally, 5,660,795 contigs (N50 = 372 bp; average length = 257 bp) were obtained.

A double-digest RAD library (Peterson et al. 2012) was prepared according to Peterson’s protocol with slight modifications (Sakaguchi et al. 2015) and sequenced using HiSeq 2500 (Illumina) which generated single-end 51-bp reads. Genomic DNA was digested using EcoRI and BglII (New England Biolabs). RAD sequencing reads were filtered using Trimmomatic ver. 0.32 (Bolger et al. 2014) to remove the Illumina adapter and to eliminate low-quality regions (parameters used: -phred33 ILLUMINACLIP TruSeq3-SE.fa:2:30:10 LEADING:19 TRAILING:19 SLIDINGWINDOW:30:20 AVGQUAL:20 MINLEN:51).

The filtered RAD sequencing reads were aligned to the L. petersii draft genome using BWA-MEM algorithm in BWA ver. 0.7.8 (Li and Durbin 2009) with default settings. Then, SNPs were called using SAMtools ver.1.5 (Li et al. 2009) mpileup (−Q 20) and bcftools, which were filtered using VCFtools ver. 0.1.12 (--minQ 20 --remove-indels --maf 0.05–max-alleles 2 --min-meanDP 15–minDP 10) (Danecek et al. 2011). SNP loci that were not present in 70% of individuals in any one population were then removed using the script pop_missing_filter.sh of dDocent package (Puritz et al. 2014). We searched potential paralogous loci that showed departure from Hardy Weinberg Equilibrium (HWE; p < 0.005) in all populations using GenoDive ver. 2.0b27 (Meirmans and Van Tienderen 2004), but no loci showed departure from HWE.

Population genetic analyses

We calculated genetic diversity of all populations according to effective numbers of alleles and expected and observed heterozygosity using GenoDive. Average site-pi values in all populations were calculated using VCFtools. We then evaluated genetic relationships among populations by constructing a neighbor joining (NJ) tree based on Nei’s standard genetic distances (Ds; Nei 1972). Genetic distances were calculated using adegenet (Jombart 2008) and hierfstat (Goudet 2005) packages for R, and the NJ tree was constructed using the ape package (Paradis et al. 2004) for R. ADMIXTURE was used to detect admixture between the two lineages (Alexander et al. 2009). To eliminate the effects of linkage disequilibrium, this analysis was performed on a dataset of only one SNP per contig. ADMIXTURE was run with K assumptions ranging from 1 to 11 using default parameters and cross-validation scores were used to determine the best K value (Alexander and Lange 2011). In ADMIXTURE analysis, the threshold q values for identifying genomic introgression were set to 0.10 < q < 0.90 because this threshold has reported high efficiency (Vähä et al. 2006). To detect the signature of introgression robustly (Nussberger et al. 2013), hybrid indices (HI) of individuals in introgressed populations were calculated using a set of divergent (diagnostic or semi-diagnostic) SNP loci (FST > 0.9; 37 SNP loci) between PO (Shimizu, Nachikatsuura, and Tosashimizu) and JS (Ajigasawa and Karatsu) lineages using GenoDive. FST values were calculated using VCFtools. HI values between 0.10 and 0.90 were set for identifying individuals with introgression according to previous studies (Trier et al. 2014; Pinheiro et al. 2015).

Morphological analysis and common garden experiments

We examined body sizes and numbers of vertebrae in mature individuals that were ascending to rivers to reproduce or immediately before reproductive migration into rivers in five introgressed populations. Measurement procedures were performed as reported elsewhere (Kokita and Nohara 2011). Given that females usually have larger bodies than males (Matsui 1986), we analyzed their body sizes separately. Vertebral counts did not differ between sexes and the data were pooled (Kokita and Nohara 2011). Data for the JS and PO populations was retrieved from the study by Kokita and Nohara (2011).

Comparisons of body sizes among two lineages and five introgressed populations were performed using analysis of variance (ANOVA) followed by Tukey’s post-hoc honestly significant difference test. Assumptions of homogeneity of variance were confirmed for log-transformed male body sizes and untransformed female body sizes using Levene’s test (P > 0.05). Levene’s test showed significantly differing variances of vertebral counts even after log or square-root transformation. Comparisons of vertebral counts were performed using nonparametric Kruskal–Wallis tests, followed by post-hoc Games–Howell tests. Statistical analyses were performed using R software.

Our morphological analyses showed a tendency of introgressed populations in Joban-Kashimanada (Namie and Hitachi) to have larger body sizes than the two parent lineages (refer to the Results section). Thus, we examined whether this phenotype is heritable using common garden experiments, which were performed as in a previous publication (Kokita et al. 2017). Mature fish that were ascending to rivers for reproduction were collected during May 2006 from the Naraha population (Kido river), which is 10 km away from the Namie population (see Fig. 1b). Two glass tanks (65-L) were used for reproduction in captivity and to obtain newly hatched larvae. In total, 20 males and 20 females were placed in a tank containing temperature-controlled (14°C) freshwater and 10 stones on a sandy bottom. Approximately 1000 newly hatched larvae were obtained from four to six collected egg clutches and were transferred to rearing tanks (100-L) with temperature-controlled flow-through seawater and a natural photoperiod. The rearing tanks were placed in a temperature-controlled water bath (approximately 23°C) and fish were fed ad libitum with live omega-3 highly unsaturated fatty acid (n-3 HUFA)-enriched L-type rotifers (Brachionus plicatilis; approximately 10 individuals mL−1) for the first 50 days after hatching (DAH) and with Artemia nauplii (approximately 10 individuals mL−1) from 40 DAH. At 90 DAH, 15 individuals from each population were fixed in 10% neutralized formalin and body sizes were measured in terms of total lengths. As described by Kokita et al. (2017), common garden experiments were performed using JS and PO populations in the same year; body size data for these populations were analyzed in this previous work. Kokita et al. (2017) also counted vertebral numbers in these individuals and in JS and PO individuals, and trait values of these populations under common garden rearing conditions were compared statistically.

Detection of introgression of NPY gene

We found discordance patterns between nuclear genomic ancestory and body size in most introgressed populations (refer to the Results section). This result suggests that introgression of genes contributes to the body size divergence between the JS and PO lineages. To explore this possibility, we genotyped a minisatellite marker within the intron of the NPY gene (hereafter, we refer to this as the NPY locus). In this NPY locus, a divergent selection occurred between the JS and PO lineages, and unlike in the PO lineage, upregulation of NPY gene expression was shown in the JS lineage in a common garden environment (Kokita et al. 2013). The NPY gene is an important regulator of food intake and is a potent orexigenic agent, and common garden experiments suggested that higher intrinsic growth rates and food intake result in larger body sizes of the JS lineage (Kokita et al. 2017). Therefore, the NPY gene would be a good candidate that is associated with the larger body size of the JS lineage (Kokita et al. 2013). We genotyped the NPY locus for introgressed populations in Joban-Kashimanada (Namie and Hitachi) and Seto Inland Sea (Saiki, Higashihiroshima, and Anan) using the methods described by Kokita et al. (2013). Genotype data for the NPY locus in the four populations was used as reference genotypes (JS, Ajigasawa and Tsuruga; PO, Shimizu and Nachikatsuura; Fig. 1b; Kokita et al. 2013). Genetic diversity at the NPY locus (effective number of alleles and observed heterozygosity) was calculated using GenoDive. Similarities in the allele frequencies between the two lineages and introgressed populations were estimated by constructing a NJ tree based on Goldstein’s δµ2 distance (Goldstein et al. 1995) using Populations (Langella 1999).

Results

Mitochondrial DNA analysis

The ratio of the two mtDNA clades was variable in introgressed populations in the secondary contact zones (Fig. 1c). Genetic diversity within each clade of each population (Supplementary Table 1) was similar to that of the surrounding pure JS or PO populations, as shown in previous phylogeographic studies (Kokita and Nohara 2011). A phylogenetic tree of mtDNA haplotypes from ours and other studies is shown in Supplementary Fig. 1.

RAD sequencing

HiSeq 2500 sequencing yielded 46,978,862,874-base quality-filtered reads. The average number of quality-filtered reads per individual was 3.7 million (range, 1.2–7.6 million). The mean depth per individual was 84.5-fold (range, 22.3–170.2-fold). Finally, 2,262 SNP loci (1,265 first SNP loci in each contig) in 246 individuals were genotyped.

Population genetic analyses

Although genetic diversity was similar among populations, effective allele number and average site-pi values of Hitachi and Namie intogressed populations in Joban-Kashimanada tended to show significantly higher values than that showed by the two parent lineages (Welch’s t test, P < 0.01; Table 1). Overall homozygote excess was observed in all populations (P < 0.005). There were no SNP loci that showed departure from HWE in introgressed populations characteristically. The population tree showed that Hitachi (Joban-Kashimanada), Namie (Joban-Kashimanada), and Anan (Seto Inland Sea) introgressed populations were genetically similar to the PO lineage and that Saiki (Seto Inland Sea) and Higashihiroshima (Seto Inland Sea) populations were genetically similar to the JS lineage (Fig. 2a).

Fig. 2
figure 2

a Neighbor joining tree of 11 populations calculated using Nei’s DS with RAD sequencing data. b Effective numbers of alleles and observed heterozygosity in NPY locus in each population. c Neighbor joining tree of five introgressed populations and PO and JS lineages based on Goldstein’s δµ2 distances, as calculated from NPY data

Full size image

In our clustering analyses with ADMIXTURE, the values of CV error suggested that K = 3 was the best fit (Supplementary Fig. 2). Assuming K = 2, PO and JS lineages were assigned to different clusters, all individuals in Anan introgressed populations showed the sign of introgression (0.1 < q < 0.9) (Fig. 1c). Introgression was detected in the Natori population where introgression of mtDNA and nuclear gene was not detected in a previous study (Kokita and Nohara 2011). However, other introgressed populations showed little introgression except for one individual each in Hitachi and Higashihiroshima introgressed populations (Fig. 1c). Assuming K = 3, Namie and Hitachi introgressed populations in Joban-Kashimanada were clearly assigned to different cluster (Fig. 1c). If we assume K = 4, JS populations were divided in to two genetic clusters (Supplementary Fig. 3).

Although ADMIXTURE analysis did not detect signs of introgression in most individuals of Saiki, Higashihiroshima, Namie, and Hitachi introgressed populations, hybrid indices based on divergent SNP loci showed the signs of introgression (0.1 < HI < 0.9) in 3.1−43.8% individuals in these introgressed populations (3.1% in Namie, 15.6% in Hitachi, 43.8% in Saiki, and 12.5% in Higashihiroshima) (Fig. 1c).

Morphological analyses and common garden experiments

On multiple comparisons using ANOVA, body sizes were observed to significantly differ among the two lineages and five introgressed populations (males, F6,206 = 50.49, P < 0.0001; females, F6,168 = 96.62, P < 0.0001). Interestingly, despite their nuclear genomic ancestry assigned to the PO lineage (Fig. 1c), introgressed populations in Joban-Kashimanada tended to have large bodies that were similar or larger than those of the JS lineage (Namie, mean ± SD = 40.9 ± 1.8 mm for males, 46.2 ± 1.5 mm for females; Hitachinaka, 42.2 ± 1.2 mm for males, 47.0 ± 1.5 mm for females) (Fig. 3). In contrast, body sizes of introgressed populations in Seto Inland Sea (Anan, 39.2 ± 1.4 mm for males, 42.9 ± 1.0 mm for females; Higashihiroshima, 39.6 ± 0.9 mm for males, 42.4 ± 1.2 mm for females; Saiki, 37.3 ± 1.2 mm for males, 40.9 ± 1.2 mm for females) were approximately intermediate between those of JS and PO lineages (Fig. 3). It is notable that body sizes of Saiki and Higashihiroshima populations become smaller (Fig. 3) despite their nuclear genomic ancestry to the JS lineage (Fig. 1c). These trends in body size were almost significant in females but were not in males (Post-hoc Tukey’s honestly significant difference test, P < 0.05; supplementary tables 2 and 3), and body sizes were more variable among females than males. Collectively, body size in introgressed populations was not explained by their nuclear genomic ancestry.

Fig. 3
figure 3

Bubble plot of vertebral numbers and boxplot of body sizes of five introgressed populations, JS, and PO lineages. Nuclear genomic ancestry and NPY genotypic feature of each introgressed population were also shown (black: JS-like, white: PO-like, gray: intermediate)

Full size image

Vertebral numbers also differed significantly among the two lineages and five introgressed populations (Kruskal–Wallis chi-squared = 180.9967, df = 6, P < 0.001; Fig. 3), and vertebral numbers in introgressed populations in Joban-Kashimanada (Namie: 33, 100%; Hitachi: 32, 5.0%; 33, 85.0%; 34, 10.0%) were similar to those in the PO lineage (32, 6.7%; 33, 90.0%; 34, 3.3%), whereas those of Higashihiroshima (34, 65.0%; 35, 35.0%) and Saiki (33, 3.3%; 34, 93.3%; 35, 3.3%) populations from the Seto Inland Sea were similar to those of the JS lineage (33, 1.7%; 34, 81.7%; 35, 16.7%). In contrast, vertebral numbers in individuals of the Anan population were similar to those of both lineages, albeit more commonly to those (33) of the PO lineage (32, 5.0%; 33, 60.0%; 34, 35.0%). These trends were supported by Games–Howell tests (P < 0.05; Supplementary Table 4). Therefore, vertebral numbers in introgressed populations were explained by their nuclear genomic ancestry in contrast to their body size.

Body sizes were significantly greater in the introgressed population in Joban-Kashimanada than in either of the parent lineages after rearing in the common garden (P < 0.05; Supplementary Table 5; Supplementary Fig. 4). Vertebral numbers of the population were similar to those of the PO lineage (P = 0.13; Supplementary Table 6; Supplementary Fig. 4).

Introgression of NPY gene

Body size of the introgressed populations was not explained by their nuclear genomic ancestry. Therefore, we then evaluated the possibility of introgression of NPY gene, which possibly contributes to the body size divergence between the PO and JS lineages. Genetic diversity of the NPY locus varied in the introgressed populations; Saiki and Anan in Seto Inland Sea showed relatively high values close to that of PO lineage, whereas Higashihiroshima in Seto Inland Sea and Hitachi and Namie in Joban-Kashimanada showed relatively low values close to that of JS lineage (Fig. 2b). Surprisingly, in the population tree based on allelic composition in the NPY locus, Saiki population was extremely close to the PO lineage despite its JS-like genomic composition, whereas Hitachi and Namie populations were close to the JS lineage despite their PO-like genomic compositions (Fig. 2c). The Anan population was intermediate despite its PO-like genomic composition. These results suggest that introgression of NPY gene affected body size change in the introgressed populations (Fig. 3).

Discussion

Overall homozygote excess in L. petersii populations identified by RAD sequencing

In this study, overall homozygote excess was observed in all populations of L. petersii, and this tendency was similar even when we preliminarily identified SNP loci using the maximum-likelihood statistical model implemented in the Stacks pipeline (Catchen et al. 2011). Although heterozygotes may have been underestimated because of various factors during RAD library preparation and genotyping methods (Baxter et al. 2011; Davey et al. 2013; Mastretta‐Yanes et al. 2015), homozygote excess in L. petersii populations was also observed in a previous population genetic study using microsatellite DNA (Kokita et al. 2013). Therefore, common factors in L. petersii may have affected the two types of analyses. In the previous study, Kokita et al. (2013) showed high-level heterozygosity in microsatellite loci of L. petersii, suggesting large effective population size in this species. Therefore, allelic dropout due to high-level polymorphism within primer sites for microsatellite DNA or the restriction site for RAD sequencing may make it impossible to observe the associated allele (null allele) in both analyses (Chapuis and Estoup 2006; Gautier et al. 2013). Alternatively, species-specific factors, such as population substructure and inbreeding (Brookfield 1996), may lead to homozygote excess. In any case, the influence on the pattern of genomic admixture in this study is considered to be limited because homozygote excess is detected regardless of the individual being in the JS or PO lineage.

Mitonuclear discordance in both secondary contact zones

Introgressed populations harbored mtDNA of both lineages (JS and PO), although the dominant mtDNA is inconsistent with nuclear genomic clustering: Saiki and Higashihiroshima harbor more of the PO mtDNA in comparison with their JS-like nuclear genomic ancestry, whereas Hitachi and Namie harbor more of the JS mtDNA despite being more genomically similar to the PO lineage. However, Anan is an exception as its mitochondrial and nuclear genomic composition is almost consistent. The conflicting geographic patterns between mitochondrial and nuclear genomes, i.e., mitonuclear discordance, have been commonly observed across a wide range of organisms (Gompert et al. 2008; Currat et al. 2008; Larmuseau et al. 2010; Toews and Brelsford 2012; Morales et al. 2015). Collectively, in comparison with nuclear genomic compositions, our results showed that PO mDNA introgresssed more widely in the Seto Inland Sea and that JS mtDNA introgresssed more widely in the Joban-Kashimanada. These asymmetric patterns observed in both secondary contact zones suggest the inevitability of this pattern in L. petersii and existence of reliable mechanisms that led to the patterns.

Several hypotheses have been proposed to explain asymmetric introgression of mitochondrial genome (Toews and Brelsford 2012). For example, demographic difference between the two lineages, including population or range size, can promote asymmetric introgression of one of the mtDNA lineage (Currat et al. 2008). However, because mtDNA lineages that show asymmetric introgression patterns are different between Seto Inland Sea and Joban-Kashimanada, this possibility was considered to be unlikely. Chan and Levin (2005) showed that incomplete pre-mating barriers, such as strong mating preferences, can promote asymmetric introgression of mtDNA. According to their model, female preferences for hybrid males promoted biased introgression of mtDNA. Kokita et al. (2017) reported that mating preferences of L. petersii are influenced by body sizes and that JS females demonstrated preferences for large males of the same lineage, whereas PO females did not. Our phenotype analyses showed that introgressed populations in Seto Inland Sea have body sizes that are intermediate between those of individuals in the two lineages and introgressed populations in Joban-Kashimanada have larger body sizes than individuals of either parental lineage. Hence, JS lineage females possibly prefer males of populations in Joban-Kashimanada, likely leading to the biased JS mtDNA introgression pattern in these population with PO-like nuclear genomic composition, as observed here. Conversely, JS lineage females may not prefer males of populations in Seto Inland Sea and PO lineage females have no preference and mate freely, likely leading to the biased PO mtDNA introgression pattern in these populations with JS-like nuclear genomic composition, which was also observed here. In addition, given that smaller body size in Seto Inland Sea and larger body size in Joban-Kashimanada are enhanced in females, local selective advantage for females may be associated with the introgression of the mitochondrial genome; PO-lineage females with small body size have selective advantage in Seto Inland Sea, whereas JS-lineage females with large body size have advantage in Joban-Kashimanada. Therefore, mating preference or/and local adaptive pressure may cause asymmetric introgression of mtDNA in both secondary contact zones in theory. At this stage, we cannot discuss the possibilities of other mechanisms, such as adaptive introgression, hybrid zone movement, and genomic incompatibility (Toews and Brelsford 2012). Further detailed analyses are needed to explain discordances between nuclear and mtDNA markers in the ice goby.

The high level of homogeneity in ancestry proportions and mitonuclear discordance in the introgressed populations of ice goby suggest that timing of hybridization was not recent. Previous studies have suggested that secondary contact between the two lineages occurred because of opening of the Japan Sea during any of the recent interglacial periods (Kokita and Nohara 2011; Hirase and Ikeda 2015). Among the two secondary contact zone, given that the Seto Inland Sea was terrestrial during Last Glacial Maximum, it is probable that the secondary contact zone in this sea was formed as a consequence of postglacial colonization and/or range expansion through the Kanmon Strait after this period (Kokita and Nohara 2011). Therefore, asymmetric mitochondrial genome may have occurred in the past several thousand years.

Relationships between genotype and phenotype in Seto Inland Sea introgressed populations

Our analyses of the Anan introgressed population in Seto Inland Sea showed admixture of the genomes of the two parental lineages, and the presence of vertebral number patterns of both PO and JS lineages (33 and 34, respectively). In contrast, the vertebral numbers of Saiki and Higashihiroshima introgressed populations in Seto Inland Sea, which were genetically close to the JS lineage, were predominantly 34, characteristic of the JS lineage. Therefore, expression patterns of vertebral numbers in introgressed populations corresponded with their nuclear genotypes. However, expression patterns of body sizes were not associated with vertebral numbers and nuclear genotypes; individuals of Higashihiroshima and Saiki introgressed populations were significantly smaller than those of the JS lineage.

Expression patterns of the present phenotypes of body size and vertebral numbers have two important implications, and the first of these is that their genetic determinants may differ. Vertebrae numbers and body sizes in fish varies geographically among populations and species and most of these variations can be explained by Jordan’s rule (Jordan 1891), which describes the tendency of fish from higher latitudes and colder waters to have greater numbers of vertebra and larger body sizes than their lower-latitude counterparts. These phenotypes often demonstrate non-adaptive or adaptive plasticity in relation to environmental factors (Billerbeck et al. 1997; Angilletta et al. 2004). However, common garden experiments showed a degree of genetic determinism in the ice goby as shown in some other studies (Billerbeck et al. 1997; Conover et al. 2009), although there were some limitations in these experiments due to nektonic feature and scaleless body of this species (Kokita et al. 2017). Given that body sizes and numbers of vertebrae are correlated (“pleomerism”; Lindsey 1975), it is difficult to assess whether their genetic determinants differ without performing breeding experiments (Billerbeck et al. 1997). The segregation patterns observed here, however, suggest that their genetic bases differ.

The second important implication of our phenotype observations is that they may reflect different gene introgression patterns. In our experiments, expression patterns of body sizes in the introgressed population were unrelated to overall genotypes, whereas vertebral numbers were related to overall genotypes. Hence, genetic variants contributing to body sizes and neutral variants have not moved together, whereas those of vertebral numbers have, reflecting differing adaptive values. Functional relationships between vertebrae numbers and swimming performance were suggested in early studies (Lindsey 1975; Spouge and Larkin 1979), and were supported by observations of selective predation for vertebral numbers in Gasterosteus aculeatus (Swain and Lindsey 1984) and Mylocheilus caurinus (Swain 1988). More recently, differing vertebral phenotypes were correlated with differences in burst swimming speeds (Swain 1992), and recent studies by Yamahira et al. (2006) suggest an adaptive advantage relating to size-dependent winter mortality. Specifically, thermal constraints may favor populations with larger larvae, which tend to have higher numbers of vertebrae and better chances of survival. Whereas the adaptive relevance of vertebral numbers remains unclear (McDowall 2008), body size is a key determinate of physiological and life-history traits of organisms (Peters and Wassenberg 1983; Schmidt-Nielsen 1984; Bonner 2011). Accordingly, body size divergence often promotes pre-mating isolation between populations as a by-product of adaptation to different environments (Schluter 2001; McKinnon et al. 2004). Various selective pressures have also been associated with body size divergences, including abiotic and biotic factors (e.g., food availability, inter-specific competition, predation risk, and river migratory effectiveness) (Hess et al. 2014; Ho et al. 2009; Amarillo‐Suárez et al. 2011). Although we should recognize that the genetic basis of larger body sizes may be obscured by environmental factors that contribute to body size variations (Billerbeck et al. 1997; Angilletta et al. 2004), our results suggest that body size is subject to selective pressures and adaptive introgression in the Seto Inland Sea. We then evaluated the possibility that the previously identified NPY gene contributes the intermediate body size. As expected, we found enhanced introgression of PO alleles in Saiki population that shows body size closest to PO lineage despite its JS-like genomic composition. Therefore, adaptive introgression of NPY gene may contribute to body size variation in Seto Inland Sea.

Seto Inland Sea is a semi-enclosed sea with shallow depth and is characterized by strong seasonal variation in salinity and seawater temperature due to its environmental features (Chang et al. 2009). Kokita et al. (2013) speculated that past and present low seawater temperatures promoted the evolution of higher growth performance, leading to larger body sizes of JS individuals. Assuming that the seawater temperature is also a driving force of body size changes in Seto Inland Sea, intermediate or smaller body size individuals may have advantages in these unstable environments of Seto Inland Sea. Also, a field study has reported that the spawning ground of Higashihiroshima population is locally formed in the downstream area of river closer to the estuary area than other populations (Hasegawa and Shoji 2017). Therefore, river migratory effectiveness may be associated with the body size changes in Seto Inland Sea. Examining the extent to which ecological variables affect selection on body size is difficult, and further detailed field and experimental studies are needed.

Possibility of introgression in the population in Joban-Kashimanada

Despite little admixture being detected in Hitachi and Namie populations in Joban-Kashimanada and their genomic composition being extremely close to that of those in the PO lineage, these populations showed large body size that was similar to that observed in the JS lineage. Given that the Joban-Kashimanada area is characterized by higher biological production due to the intersection of warm Kuroshio and cold Oyashio currents (Yokouchi et al. 2000), environmental contributions to body size may be significant. Although we cannot completely rule out the possibility of maternal effects (Heath et al. 1999) by using F1 individuals, larger body sizes of individuals in this area were also reproduced in common garden experiments, suggesting that this trait is genetically determined. In addition, we found that genotypes of Namie and Hitachi populations in NPY locus were close to those of the JS lineages, which were characterized by reduced genetic diversity (Kokita et al. 2013). Given that the overall genomic diversity of populations in Joban-Kashimanada was not low, low genetic variability of the NPY gene may be locus-specific.

On the possible cause of larger body size in the Joban-Kashimanada populations other than introgression, de novo mutations (Hoekstra et al. 2006) or ancestral polymorphisms (Colosimo et al. 2005) may have provided the allelic variation necessary for the evolution of larger body sizes; however, the processes that have actually resulted in the dominance of these genes in the population are actually either genetic drift or selection or both. Because NPY genotypes of Namie and Hitachi populations were close to those of the JS lineage, common variants near NPY gene derived from ancestral polymorphisms possibly led to their larger body size independently. However, given that mtDNA in Joban-Kashimanada did not suggest any signature of past isolation, such as sub-clades, and that introgression of mtDNA and nuclear MYH6 gene were found (Kokita and Nohara 2011), slight introgression of genes of the JS lineage that contributes to body size may be a reasonable scenario. It is difficult to detect genomic introgression because introgressed individuals share a large part of their genome with one of the parental lineages, similar to the introgressed populations in the ice goby; therefore, more genome-wide markers are required to find the signature of genomic introgression (Nussberger et al. 2013).

Kokita et al. (2013) speculated that past and present low seawater temperatures promoted the evolution of higher growth performance, leading to larger bodies of JS individuals. In fact, the Joban-Kashimanada area would have also been subject to these cooler environments, and a cool Oyashio current often dominated the Joban-Kashimanada area during the Pleistocene glacial periods. Currently, the warm Kuroshio current dominates these areas (Nimura et al. 2006), and the historical dynamics of the Kuroshio and Oyashio currents presumably caused sudden expansion of intertidal goby Chaenogobius annularis populations, and distribution around the Joban-Kashimanada area (Hirase et al. 2012). Moreover, cool Oyashio currents often reach this area even in the present time (Yokouchi et al. 2000). Therefore, past and present cool seawater temperatures might promote the evolution of large body sizes with introgression in the Joban-Kashimanada area, unlike in the case of the Seto Inland Sea, and following isolation from the PO and JS lineages. An alternative explanation for the increased body size in Joban-Kashimanada populations would be competition between introgressed populations and JS individuals as shown in the Neotropical shorebird species by Lipshutz et al. (2018). Further genome-wide studies are needed to gain more insight into the relationship between large body size observed in Joban-Kashimanada and introgression.

Conclusion

We performed integrated analyses of genotypes and phenotypes for two independent secondary contact zones of the ice goby and obtained several important insights. First, we found mitonuclear discordance in introgression patterns in both secondary contact zones, and that body size changes in these zones possibly affect the consequences. Second, we showed that body size changes were not explained by nuclear genomic ancestory in both secondary contact zones and this suggested that introgression of genes regulate body size, such as NPY gene. Third, our findings demonstrated that secondary contact contributes to various genomic and phenotypic consequences in different marine environments. This study provides a valuable empirical example into the processes through which biodiversity is increased following geographic isolation in marine environments.

Data archiving

Mitochondrial DNA sequences were submitted to the DNA Data Bank of Japan (DDBJ; LC465830−LC465988; http://www.ddbj.nig.ac.jp). All RAD sequencing data were deposited in the DDBJ; accession numbers DRR151244–DRR151489.