Sample collection
All sticklebacks were collected with seine nets and minnow traps as described previously [20, 30-32, 37, 66, 67] (Fig. 1A). After euthanasia with an overdose of MS-222 (0.5 g/L), the pectoral fins were dissected out and preserved in 99% ethanol until use. Supplementary Tables S1 and S2 provide details of the samples. Morphologically identified species [33] collected at the same locality were denoted as different populations, with the exception of Okinebe (Oki), where G. aculeatus, G. nipponicus, and possibly their hybrids are supposed to be included.
Laboratory experiments and sequencing
Genomic DNA was isolated using a DNeasy Blood & Tissue Kit (QIAGEN, Valencia, CA, USA). Double digest RAD sequencing (ddRAD-seq) was performed as described previously [68]. Briefly, 10 ng of genomic DNA was digested with EcoRI and BglII, followed by adapter ligation and amplification with uniquely barcoded primers. The libraries were run on HiSeq 2000 or 2500 using the 50 bp single-end or a 100 bp paired-end mode at Macrogen (Kyoto, Japan) or the Advanced Genomics Center of the National Institute of Genetics (Shizuoka, Japan). The sequence data are available from DDBJ/EMBL-EBI/NCBI Sequence Read Archive (DRA010673). Some of the ddRAD-seq data has been published previously [15] (see Table S2).
Additionally, we used publicly available whole genome sequence (WGS) data (Table S2). For G. aculeatus collected from PO Akkeshi (P4), G. nipponicus from JS Akkeshi (J4), and Gasterosteus wheatlandi, we used the previously reported whole genome sequences [69]. Sequence data of G. aculeatus from FW Aisaka (FN4) and FW Towada (FN5) were derived from a previous study [32]. For G. aculeatus from northern Europe, the sequences of two randomly selected samples from the marine population reported previously [70], and those of one or two randomly selected samples from each freshwater population reported previously [71] were obtained.
Sequence data processing
The flow of bioinformatic analyses is summarized in Fig. S1. Trimming of ddRAD-seq reads was performed to remove adapter sequences and failed reads using Trimmomatic v0.39 [72] with the following parameters: “ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2 CROP:50 LEADING:3 TRAILING:3 MINLEN:50”. The trimmed reads were mapped to the BROADS S1 stickleback reference genome sequence (soft-masked, Ensembl 99) using NextGenMap v0.5.5 [73]. Variants were called with FreeBayes v1.3.2 [74], skipping sites with the average coverage per sample exceeding 500 and with the options: “--report-monomorphic --use-mapping-quality --use-best-n-alleles 8”. Sites of a sample with a coverage of less than five were discarded with BCFtools v1.9 [75].
We further selected RAD loci with the following criteria using BCFtools and bedtools v2.17.1 [76]. First, the sites genotyped in less than 25% of the samples and located on the mitochondrion were excluded. Next, we searched for the regions consecutively genotyped for at least 40 bp, allowing gaps not longer than 10 bp. The records within the identified RAD regions were extracted and variant representations were normalized with vt v0.5772 [77].
WGS reads were trimmed using Trimmomatic with the following settings: “ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:2:true LEADING:3 TRAILING:3 MINLEN:30”. Overlapped paired-end reads were merged with PEAR 0.9.10 [78] and read pairing was confirmed with fastq-pair v1.0 [79]. Mate-pair reads from Feulner et al. [70] were reversed and complemented using SeqKit 0.10.0 [80]. The reads were mapped to the BROADS S1 stickleback reference sequence using NextGenMap v0.5.5. The maximum insert size for the alignments of the mate-pair reads was set to 6000. Duplicate reads were marked with Picard Tools v2.21.8 [81]. Variants within the selected contiguous RAD loci (see above) were called with FreeBayes using the same settings as that of ddRAD-seq. Sites of a sample with a coverage of less than five were discarded with BCFtool, and normalization of variants was conducted with vt.
The pre-processed variant calls from ddRAD-seq and WGS were merged using BCFtools. Block substitutions were decomposed into their constituent SNPs using vt. Indels, invariant sites, and sites on the sex chromosomes of G. aculeatus and G. nipponicus (Chromosomes IX and XIX) or those in masked regions or on ambiguous nucleotides in the reference sequence were discarded. Samples with excessively missing genotypes (> 80%) were excluded with BCFtools. This process resulted in a dataset of 97145 SNPs genotyped in a total of 310 samples.
Population structure analyses
In order to investigate genetic differentiation and potential introgression among the stickleback populations in the Japanese Archipelago, we first used a model-based likelihood clustering algorithm implemented in ADMIXTURE v1.3.0 [82]. We selected biallelic SNPs that were genotyped in all populations with only one missing population allowed, and that were missing in less than 30% of the overall samples with VCFtools v 0.1.17 [83]. If an allele at a SNP site was found in only one sample, the SNP site was excluded regardless of whether it was identified as “singleton” or “doubleton” with VCFtools. The SNPs were subsampled with VCFtools to maintain a minimum distance of 1 kb to reduce the effect of linkage between SNPs. The input file for ADMIXTURE including 2735 SNPs was created using PLINK v1.90 [84]. ADMIXTURE was run by varying the number of evolutionary clusters K from one to nine. The results were summarized and visualized using CLUMPAK [85] on the web (http://clumpak.tau.ac.il/index.html).
We also conducted principal component analyses (PCA), using the adegenet v2.1.1 package [86, 87] of R [88]. The dataset for the ADMIXTURE analysis was further filtered, keeping SNPs with minor allele frequency ≥ 0.03 and individuals with missing genotypes < 20%. This resulted in a dataset of 813 SNPs.
Phylogenetic analyses
Maximum likelihood (ML) phylogenetic trees were constructed with RAxML-NG v0.9.0 [89] based on two datasets of concatenated SNPs. The first includes 1,919 SNPs of all the samples from Japan excluding putative recent hybrids between G. aculeatus and G. nipponicus, which would violate basic assumptions of phylogenetic reconstruction methods and bias tree topology and branch lengths. Hybrid individuals were identified using the ADMIXTURE analysis described in section 2.4 based on Q values assuming K = 2. When both a Q value for the G. aculeatus cluster and that for the G nipponicus cluster at K = 2 were < 0.875, that individual was classified as a hybrid. The identified hybrids were concordant with those detected by PCA (Supplementary Fig. S2, left panel).
The second tree consists of 3,717 SNPs of non-hybrid individuals from native populations of western and eastern basins of the Pacific and Europe, which were subsampled to include at most six individuals per population, two of which had the least missing genotypes and the rest of which were randomly selected. Hybrids in the East Pacific populations were identified and excluded in a manner similar to that of Japanese samples using ADMIXTURE. We selected biallelic SNPs that were genotyped in all four populations using BCFtools. The SNPs were subsampled with VCFtools to maintain a minimum distance of 1 kb. ADMIXTURE was run by varying the number of evolutionary clusters K from one through four. Identification of hybrid individuals was conducted based on Q values from the ADMIXTURE analysis assuming K = 3. If the Q values of an individual for any cluster did not exceed 0.875, it was classified as a hybrid. As a result, three putative hybrids from Duwamish and Big Soos were removed (Supplementary Fig. S3). The reference sequence obtained from a freshwater stickleback collected at Bear Paw Lake, Alaska [18] was added as a sample in the second dataset. Each dataset included the SNPs genotyped in all but one population and > 70% of the overall samples, keeping the minimum distance between SNPs at 1 kb. We used the general time-reversible model of nucleotide substitution with gamma-distributed rate heterogeneity and ascertainment bias correction [90] using the conditional likelihood method [91]. We conducted bootstrap analyses with 200 replicates and searched for the best scoring trees in each of the two runs. The tree was visualized with FigTree v1.4.4 [92].
Phylogeny and divergence time among stickleback populations was estimated with the multispecies coalescent model using the Bayesian framework of SNAPP v1.5.0 [93] implemented in Beast v2.6.2 [94]. To reduce the computational time, we selected two non-hybrid individuals with the least missing genotypes from 13 representative populations covering the distribution range and distinct lineages of the stickleback. They consisted of G. nipponicus, marine populations of G. aculeatus from the western and eastern basins of the Pacific and Europe, freshwater populations from each of the three regions, including those comprising highly supported clades in the Japanese Archipelago that were revealed by the ML tree analysis. We removed SNPs with missing genotypes, and subsampled SNPs to maintain a minimal distance of 1 kb. This resulted in a dataset of 2,022 biallelic SNPs.
Root divergence was used as the calibration point. We adopted two previously published estimates as the time of divergence between G. aculeatus and G. nipponicus. The first is 680 thousand years (ka) BP following our previous study [37], estimated by a demographic analysis with an Approximate Bayesian Computation approach. The second was the 1.38 Ma BP [43] based on a Bayesian estimation of phylogeny and divergence time with concatenated RAD sequences. Although the potential overestimation of the latter due to incomplete lineage sorting is pointed out [42], we included it to account for uncertainty in the estimation of the divergence time, since it is close to another estimate of 1.22 Ma BP based on an ML-based demographic analysis [37], and within the 95% confidence interval of the former divergence time estimate (0.18–4.1 Ma).
Prior for the divergence time was specified to follow a log-normal distribution with means in real space to the respective divergence times (i.e., 0.68 and 1.38 Ma), and with a standard deviation of 0.18 so that 95% intervals of the two priors do not overlap. We fixed a population parameter theta, which is proportional to the product of effective population size and mutation rate per site, to be equal across lineages with a uniform prior, following Stange et al. [95]. It should be noted that fixed and equal population sizes among all populations could flaw divergence time estimates obtained in the coalescent analysis. Monophyly of G. aculeatus (i.e., all the populations except G. nipponicus) and that of two European populations were set as constraints. We used a script by Matschiner [96] to prepare input files for SNAPP. Three independent runs were performed for each calibration scheme with a chain length of 1.54–2.22 × 106 generations starting from different initial trees. Trees were sampled every 5,000 steps and checked for convergence to the stationary distribution and a sufficient effective sample size (ESS > 200) using Tracer v1.7.1 [97]. The first 10% of the trees were discarded as burn-in and the remaining trees were visualized using DensiTree v2.2.7 [98]. Maximum clade credibility consensus trees of each run after burn-in were summarized with TreeAnnotator v2.6.2 [99] and visualized with FigTree [92].