Introduction
============

Deleterious mutations arise at each generation in any organism. Although deleterious mutations can be removed from populations by purifying selection, populations with smaller effective population sizes are more likely to accumulate deleterious mutations by genetic drift ([@evaa065-B48]; [@evaa065-B81]). The accumulation of deleterious mutations results in high mutation load ([@evaa065-B1]) and can put small populations at the risk of extinction ([@evaa065-B56]; [@evaa065-B61]a, [@evaa065-B62]b). However, severe reductions in population size can lead to inbreeding and can expose recessive deleterious mutations to negative selection. This may result in purging, rather than accumulation, of deleterious mutations ([@evaa065-B50]; [@evaa065-B35]). Therefore, investigation of deleterious mutation loads in small endangered populations is important for determining which population is the most severely endangered requiring extra caution for conservation efforts ([@evaa065-B24]; [@evaa065-B46]).

Reduction in effective population sizes can occur in populations that colonize isolated habitats ([@evaa065-B106]), fragmented populations ([@evaa065-B20]), populations at the front of a range expansion ([@evaa065-B84]; [@evaa065-B85]), and domesticated strains ([@evaa065-B66]). For example, genetic disorders have been reported in endangered species maintained in captivity for conservation purposes, such as albinism in brown bears and chondrodystrophy in California condors ([@evaa065-B24]). Although deleterious mutation load has been empirically investigated by estimates of inbreeding depression and lethal equivalents ([@evaa065-B63]), recent advances in genomic technologies enable us to conduct whole-genome sequence analysis and bioinformatic analysis to search for putatively deleterious mutations across the genome ([@evaa065-B2]; [@evaa065-B46]). Phylogenetically conserved residues are generally important for protein functions, so any mutations that occur at highly conserved sites are predicted to impair protein functions and be deleterious ([@evaa065-B77]; [@evaa065-B14]; [@evaa065-B13]). Therefore, whole-genome sequence data can help us to infer the levels of accumulation of putatively deleterious mutations ([@evaa065-B46]; [@evaa065-B2]).

In addition, genome-wide data can also inform us which genomic regions and what kinds of genes accumulate deleterious mutations. Because effective population sizes can differ between genomic loci, the levels of accumulation of deleterious mutations are also expected to differ between genomic regions. For example, sex chromosomes generally have smaller effective population sizes than autosomes ([@evaa065-B102]; [@evaa065-B68]), so they may accumulate more deleterious mutations than autosomes. Genomic regions with lower recombination rates also have smaller effective population sizes due to background selection against deleterious mutations and the hitchhiking effects of beneficial alleles ([@evaa065-B98]; [@evaa065-B12]; [@evaa065-B11]; [@evaa065-B27]).

Genomic analysis has revealed an excess of deleterious mutations in human populations that have recently expanded to new areas ([@evaa065-B59]; [@evaa065-B101]; [@evaa065-B25]). Accumulation of deleterious mutations in marginal populations has been also reported in *Arabidopsis thaliana* ([@evaa065-B8]). Accumulation of deleterious mutations has been also revealed in domesticated organisms ([@evaa065-B66]), such as rice ([@evaa065-B60]), sunflowers ([@evaa065-B92]), and dogs ([@evaa065-B69]). These genomic data support the hypothesis that deleterious mutations are likely to accumulate in populations with smaller population sizes. In contrast, genome-wide sequence analysis of mutation loads in natural populations is limited except in a few cases, such as the Apennine brown bear and lake trout ([@evaa065-B6]; [@evaa065-B86]; [@evaa065-B22]).

Here, we investigated mutation loads in landlocked freshwater populations of threespine stickleback (*Gasterosteus aculeatus*) ([fig. 1](#evaa065-F1){ref-type="fig"}). Pleistocene glacial cycles created many freshwater habitats, where ancestral marine or anadromous threespine sticklebacks colonized and diversified ([@evaa065-B109], [@evaa065-B110]; [@evaa065-B5]). The genetic and genomic basis for adaptation has been extensively investigated in freshwater threespine stickleback populations ([@evaa065-B49]; [@evaa065-B54]; [@evaa065-B45]; [@evaa065-B36]; [@evaa065-B95]; [@evaa065-B83]; [@evaa065-B43]). For example, selection on standing allelic variation, which is maintained in the ancestral marine populations by migration-selection balance ([@evaa065-B96]) or balancing selection ([@evaa065-B29]), has often played important roles in adaptation to freshwater habitats ([@evaa065-B15]; [@evaa065-B71]; [@evaa065-B52]; [@evaa065-B44]; [@evaa065-B45]; [@evaa065-B4]; [@evaa065-B76]; [@evaa065-B32]). There are also cases of freshwater adaptation caused by de novo mutations ([@evaa065-B10]; [@evaa065-B80]; [@evaa065-B43]; [@evaa065-B111]). In contrast to the genetic basis for adaptation, the mutation loads of freshwater threespine stickleback populations have not yet been investigated.

![---Map of collection sites in Japan. Freshwater populations are shown by green circles and letters, whereas marine populations are shown by blue circles and letters.](evaa065f1){#evaa065-F1}

Freshwater populations of threespine stickleback in Honshu, the central part of the mainland of Japan, are an ideal system to investigate the accumulation of deleterious mutations in natural systems. Temperatures of these regions during summer are beyond the lethal maximum of most stickleback populations ([@evaa065-B110]), so the extant distribution is limited to several spring-fed ponds and streams ([@evaa065-B51]), where water temperatures are maintained around 10--20 °C by cold spring water ([@evaa065-B47]; [@evaa065-B38]; [@evaa065-B72], [@evaa065-B73]). Furthermore, these habitats are geographically isolated from the sea, so no gene flow is likely to occur from marine fish ([@evaa065-B51]). Because of construction of buildings and river banks since the 1960s, many habitats were directly destroyed and/or became deprived of spring water, resulting in the extinction of many populations ([@evaa065-B75]). The remaining few populations are endangered and listed as local populations at the risk of extinction by the Ministry of the Environment Japan ([@evaa065-B40]). Because most of these freshwater populations have unique morphological and physiological characteristics ([@evaa065-B28]; [@evaa065-B112]; [@evaa065-B41]), it is important to conserve them ([@evaa065-B75]; [@evaa065-B51]).

In this study, we first investigated the phylogeny and demographic history of landlocked freshwater populations in Honshu Island, the central part of the Japanese mainland, using whole-genome resequence (WGRS) data. We next analyzed the level of heterozygosity as an indicator of contemporary effective population size and inferred and characterized putatively deleterious mutations across the genome using two bioinformatics tools, PROVEAN and SIFT. We then tested whether independent colonization events have independently accumulated deleterious mutations and identified regions of the genome with an excess of deleterious mutations.

Materials and Methods
=====================

Whole-Genome Resequencing
-------------------------

Information for all samples used in this study is listed in [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. Briefly, for analysis of deleterious mutations, we investigated nine populations (five freshwater and four marine populations; one female individual per population unless noted). WGRS data of the five freshwater populations (Gifu-keizai1, Shiga, Aizu-tajima, Nasu, Ono; one female for each) were previously reported ([@evaa065-B116]). To confirm that the use of different individuals from the same population gave rise to qualitatively similar results, we additionally sequenced two individuals from the Gifu population, one female collected from a pond in Gifu-keizai University (Gifu-keizai2) and one female collected from a nearby site, a tributary of Tsuya River in Gifu Prefecture (Gifu-Tsuya). In the present study, Gifu-Tsuya was mainly used because the habitat is more natural than a pond in the campus of Gifu-keizai University (recently renamed as Gifu-kyoritsu University). For these two fish, Blood&Tissue DNA extraction kit (Qiagen, Valencia, CA) was used for genomic DNA isolation from female individuals. TruSeq DNA PCR-Free LT Sample Prep kit (Illumina, San Diego, CA) was used for library construction with eight cycles of polymerase chain reaction. HiSeq Rapid SBS Kit v2 (200 Cycle) and HiSeq PE Rapid Cluster Kit v2 were used for sequencing in HiSeq2500. All reads were trimmed and then mapped to the BROADS1, Feb 2006 reference genome sequence using CLC Genomics Workbench 8.5 as described previously ([@evaa065-B115]; [@evaa065-B42]; [@evaa065-B55]; [@evaa065-B91]). Single-nucleotide polymorphism (SNP) calls were conducted as described previously with minimum coverage = 10, maximum coverage = 200, and minimum count = 4. Finally, the coordinates on the genome were corrected following [@evaa065-B93].

WGRS data of the four marine populations were also reported previously ([@evaa065-B115]; [@evaa065-B42]; [@evaa065-B91]): a Japanese Pacific Ocean marine *G. aculeatus* population from Lake Akkeshi (Akkeshi *G. aculeatus*) (DRA001136), a Canadian Pacific Ocean marine population from the Little Campbell River estuary (Canada *G. aculeatus*) (DRA004937), and two marine populations of *Gasterosteus nipponicus* from Lake Akkeshi (Akkeshi *G. nipponicus*) (DRA001136) and from Lake Shinji (Shinji *G. nipponicus*) (DRA005130). For phylogenetic analysis, we additionally included an Atlantic Ocean marine population from the North Sea (PRJEB2954) ([@evaa065-B23]) and a marine population of *Gasterosteus wheatlandi* (DRA001086) ([@evaa065-B115]; [@evaa065-B91]).

Phylogenetic and Pairwise Sequentially Markovian Coalescent Analyses
--------------------------------------------------------------------

Phylogenetic and pairwise sequentially Markovian coalescent (PSMC) analyses ([@evaa065-B57]) were conducted as described previously ([@evaa065-B91]). Briefly, we first exported BAM files using CLC Genomics Workbench 8.5. For phylogenetic analysis, variants were called at all sites with samtools 1.2 and bcftools 1.2 ([@evaa065-B58]; [@evaa065-B17]) using the consensus caller. With a consensus vcf of all positions from across all individuals, we then used a previously published python script ([@evaa065-B70]) to calculate maximum likelihood trees in RAxML ([@evaa065-B99]) using the consensus sequences of 10-, 50-, and 100-kb nonoverlapping windows. To investigate whether the Gifu and Shiga populations belong to either *G. aculeatus* or *G. nipponicus*, we classified trees as either showing clustering of the Gifu and Shiga individuals with *G. nipponicus* or *G. aculeatus* ([supplementary fig. S1](#sup1){ref-type="supplementary-material"}*A*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) using a custom R script. We also calculated the fraction of the genome that shows monophyly among all freshwater populations (see [supplementary fig. S1](#sup1){ref-type="supplementary-material"}*B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) to investigate whether the phylogenetic tree supports a single origin or multiple origins of the freshwater populations.

For PSMC, we again used the samtools/bcftools 1.2 pipeline to produce consensus sequences and converted these into the PSMC format using a 100-bp window, a minimum of 10,000 sites with a Phred quality score \>20, minimum depth of 20×, and a maximum depth of 160×. PSMC was run for 30 iterations with a mean coalescent time of 15 (in units of 2*N*~e~). We estimated effective population size over 45 discrete time intervals using the command "4 + 19 × 2 + 3" and scaled our results with an autosomal mutation rate of 7.1 × 10^−9^ per site per year, assuming a single year generation time ([@evaa065-B30]).

Analysis of Heterozygosity and Deleterious Mutations
----------------------------------------------------

To analyze the five freshwater populations and four marine populations, we selected the whole-genome sequence data of one female individual for each population when multiple individuals were available. All sequence analyses of these nine individuals were performed using the same pipeline. Analysis of proportion of heterozygous sites was conducted as described previously ([@evaa065-B115]). We calculated the proportion of genotyped loci that are heterozygous, which is equivalent to multiple-locus heterozygosity in [@evaa065-B46]. For comparison between X chromosomes and autosomes, we compared between linkage group (LG)19 (the ancestral-X chromosome shared by *G. aculeatus* and *G. nipponicus*) and other linkage groups except LG9 (a neo-X chromosome in *G. nipponicus*) ([@evaa065-B53]). To test whether the ancestral-X chromosome (LG19) has lower proportion of heterozygous sites than autosomes (all linkage groups except LG9 and LG19), we performed statistical analysis using a generalized linear model with a binomial model implemented in the R statistical package ([@evaa065-B88]). A likelihood ratio test was performed between a null model with only habitat (freshwater vs. marine populations) included as an explanatory variable (df = 176) and an alternative test model with habitat and chromosome type (X chromosome vs. autosomes) being included as explanatory variables (df = 177). For sliding window analysis, 100-kb windows with no overlaps were used. For investigation of heterozygosity across the chromosomal positions, the same number of 100-kb windows (152 windows for each chromosome) was randomly selected from each chromosome. The relative chromosomal position of each data point was calculated as the position divided by the chromosomal length. In total 3,192 data points were plotted and fit with Loess curves with 95% confidential intervals using geom_smooth function in a R package, ggplot2 ([@evaa065-B107]), with method = "loess."

For analyses of deleterious mutations, we first searched for derived nonsynonymous mutations in each population using *G. wheatlandi* as an outgroup. We identified nonsynonymous SNPs using the annotation of ENSEMBL Release 83 and selected SNPs at nucleotide positions where the *G. wheatlandi* allele is the same as the reference allele of BROADS1, Feb 2006. Both PROVEAN and SIFT were used to search for putatively deleterious nonsynonymous mutations, as performed previously ([@evaa065-B114]). SIFT is reported to have higher sensitivity, but lower specificity than PROVEAN ([@evaa065-B14]), so we used both methods, as conducted previously ([@evaa065-B92]; [@evaa065-B114]).

In the present study, we defined nonsynonymous mutations as deleterious when they were predicted to be deleterious in both PROVEAN (score \< −2.5) and SIFT (*P *\<* *0.01). Nonsynonymous mutations that were predicted to be neutral in both PROVEAN (score \> −2.5) and SIFT (*P *\>* *0.01) were defined as neutral mutations. For comparison between chromosomes and between positions within chromosomes, we calculated the proportion of putatively deleterious mutations in nonsynonymous mutations. In these comparisons, the frequencies of homozygous and heterozygous SNPs were multiplied by 1 and 0.5, respectively, to reflect the relative allele frequencies per individual.

To analyze whether the nonsynonymous mutations in the freshwater populations are derived from standing genetic variation or not, the whole-genome sequence data previously generated from ten individuals of the marine population of *G. aculeatus* (Akkeshi, Japan) were used ([@evaa065-B115]; [@evaa065-B91]). Only genomic regions with coverage depth between 10 and 200 and without any deletions in all individuals were used. In this analysis, the frequencies of homozygous and heterozygous SNPs were multiplied by 1 and 0.5, respectively, as described above.

To test whether negative selection is relaxed in genes with putatively deleterious mutations, we compared the ratio of nonsynonymous with synonymous substitutions (*d*~N~/*d*~S~) between genes with and without putatively deleterious mutations. Genes under relaxed selection are predicted to have *d*~N~/*d*~S~ higher than genes under strong purifying selection but not \>1 ([@evaa065-B82]). *d*~N~/*d*~S~ analysis was conducted as described previously with a few modifications ([@evaa065-B115]). Briefly, we first generated in silico cDNA sequences using the SNPs called from the whole-genome sequence data. For heterozygous sites, we selected one nucleotide randomly. Nucleotides with insertion--deletion polymorphisms and low mapping qualities and repetitive sequences were masked. Genes with at least ten informative codons were kept for the analysis. *d*~N~/*d*~S~ of all branches of a star phylogeny of three populations, Gifu-Tsuya, Akkeshi *G. aculeatus*, and Akkeshi *G. nipponicus*, were estimated by codeml in PAML v4.5 ([@evaa065-B113]) with runmode = 0 and mode = 1. *Gasterosteus nipponicus* was used just as an outgroup of two populations.

For testing whether independent deleterious mutations occur in the same genes significantly more often than neutral mutations, we calculated the probability that randomly sampled independent mutations from two freshwater populations are colocalized on the same genes and tested whether the probability calculated for deleterious mutations was higher than that for neutral mutations. The probability, *P*~colocalization~, was given as $P_{colocalization} = \frac{\left. {\sum\limits_{G}\left( N_{G,~A} \right.}N_{G,~B} \right)}{{\sum\limits_{G}N_{G,A}}{\sum\limits_{G}N_{G,B}}}$, where *N*~G,A~ and *N*~G,B~ are the number of deleterious or neutral mutations in population A and B, respectively, in a gene G. Only private mutations were used for the analysis. To test the significance, we calculated two-sided *P* values using a permutation test (the number of permutations = 10,000). Briefly, we permutated private deleterious and neutral mutations and calculated *P*~colocalization~ for deleterious mutations and neutral mutations separately. This allowed us to determine the null distribution of odds ratio for both mutation classes. The *P* values were then calculated based on this null distribution.

Results
=======

Multiple Freshwater Colonization Events in the Central Part of Honshu Island
----------------------------------------------------------------------------

Phylogenetic trees using WGRS data showed that all freshwater individuals examined belonged to *G. aculeatus* rather than to its closely related species *G. nipponicus* ([fig. 2*A*](#evaa065-F2){ref-type="fig"}), which is consistent with previous phylogenetic analyses using allozyme ([@evaa065-B37]), microsatellite ([@evaa065-B9]; [@evaa065-B89]), and Restriction-site Associated DNA Sequencing data ([@evaa065-B43]). To confirm the genome-wide closeness between *G. aculeatus* and the Gifu and Shiga fish, whose common ancestor split from the *G. aculeatus* lineage more closely to the *G. aculeatus* and *G. nipponicus* divergence event, we classified the tree topology of phylogenetic trees of genomic fragments with 10-, 100-, and 500-kb windows. For the majority of the genome, the Gifu and Shiga fish clustered with *G. aculeatus* rather than *G. nipponicus* (0.83, 0.98, and 0.99 for 10-, 100-, and 500-kb windows, respectively) ([supplementary fig. S1](#sup1){ref-type="supplementary-material"}*A*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

![---(*A*) A genome-wide phylogenetic tree using autosomal nuclear SNPs excluding sex chromosomes (linkage group \[LG\]9 and LG19). (*B*) PSMC analysis of Japanese mainland freshwater populations. For PSMC of other populations, see [supplementary figure S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online and [@evaa065-B91].](evaa065f2){#evaa065-F2}

Phylogenetic trees also indicated that there are at least two independent freshwater lineages: one lineage composed of Gifu and Shiga populations and another composed of Ono, Aizu_Tajima, and Nasu individuals ([fig. 2*A*](#evaa065-F2){ref-type="fig"}). Tree classifications of the whole-genome-sequence phylogeny confirmed that only small fractions of the genome supported monophyly of these two freshwater lineages (0.07, 0.04, and 0.03 for 10-, 100-, and 500-kb windows, respectively; [supplementary fig. S1](#sup1){ref-type="supplementary-material"}*B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

To investigate demographic history of these Japanese freshwater populations, we conducted PSMC analysis ([fig. 2*B*](#evaa065-F2){ref-type="fig"}). Gifu and Shiga populations increased their population sizes around 200,000--300,000 years ago and then declined, followed by reincrease around 10,000--30,000 years ago. Other three freshwater populations showed clearly different demographic histories. These three decreased their population sizes during the last glacial period (10,000--70,000 years ago). The ancestral marine *G. aculeatus* from Japan maintained relatively constant population sizes as reported previously ([@evaa065-B91]) ([supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Reduction in Genetic Diversity in Landlocked Populations
--------------------------------------------------------

Freshwater individuals have significantly smaller proportion of heterozygous sites than marine populations ([fig. 3*A*](#evaa065-F3){ref-type="fig"}) (Wilcoxon rank sum test, *N *=* *4 and 5 for marine and freshwater populations, respectively, *W* = 0, *P *=* *0.00721). Comparison between the X chromosome and autosomes indicates that the X chromosome has a smaller proportion of heterozygous sites than autosomes in all populations examined ([fig. 3*B*](#evaa065-F3){ref-type="fig"}). A significant reduction on the X chromosome compared with autosomes was confirmed with generalized linear model (see Materials and Methods, coefficient = −0.23; likelihood ratio test statistics = 105.57; *P *\<* *2.2 × 10^−16^). Genome-wide sliding window analysis indicated that freshwater populations show overall reduction in heterozygous sites across the genome ([supplementary fig. S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) except that Gifu and Shiga have a few chromosomal regions with an elevated proportion of heterozygous sites even compared with the marine populations ([supplementary fig. S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online; see LG3, LG9, and LG10 in Gifu-Tsuya and LG1, LG7, and LG15 in Shiga).

![---Analysis of proportion of heterozygous sites. (*A*) Comparison among the populations. Mean values (±S.E.) of autosomes (all linkage groups except LG9 and LG19) are shown for each fish. (*B*) Comparison between an ancestral-X chromosome (LG19) and autosomes (all linkage groups except LG9 and LG19). All comparisons showed significant differences with *P *\<* *2.2 × 10^−16^ (*χ*^2^ test, df = 1). (*C*) Distribution of heterozygous sites across relative chromosomal positions in marine populatons. For each fish, the proportion of heterozygous sites within 100-kb nonoverlapping windows is plotted along the relative chromosomal positions normalized with the chromosomal lengths and is fitted with a Loess smooth line. (See also [supplementary figure S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, where data for the marine and freshwater populations are shown at different scales of *y* axis.)(*D*) Distribution of heterozygous sites across relative chromosomal positions in freshwater populations.](evaa065f3){#evaa065-F3}

Threespine stickleback has a general trend of recombination rates being lower and higher at the center and both ends of each chromosome, respectively ([@evaa065-B93]; [@evaa065-B26]; [@evaa065-B87]; [@evaa065-B94]; [@evaa065-B97]). To investigate the relationship between recombination rates and heterozygosity, we normalized the chromosomal positions of the sliding windows by adjusting the length of each chromosome to 1 and plotted the proportion of heterozygous sites against the normalized chromosomal position ([fig. 3*C* and *D*](#evaa065-F3){ref-type="fig"}). The smooth lines clearly indicate a smaller proportion of heterozygous sites at the center of chromosomes even in the freshwater populations, which exhibit substantially reduced heterozygosity. Interestingly, we observed the peaks of the proportion of heterozygous sites near the peripheries but not at the extreme ends of the chromosomes in the marine populations (see also [supplementary fig. S4](#sup1){ref-type="supplementary-material"}*A*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). This contrasted with freshwater populations, which often have their highest values at the extreme ends of chromosomes ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}*B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Recent mapping of recombination rates in threespine stickleback has indicated moderately reduced recombination rates at the further ends of the chromosomes ([@evaa065-B26]; [@evaa065-B94]; [@evaa065-B97]), which is similar to the patterns of the proportion of heterozygous sites in the marine populations ([supplementary fig. S4](#sup1){ref-type="supplementary-material"}*A*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Accumulation of Deleterious Mutations in Freshwater Populations
---------------------------------------------------------------

To investigate the levels of accumulation of deleterious mutations, we used two deleterious mutation finders, PROVEAN and SIFT. Both PROVEAN and SIFT analyses showed that the proportion of putatively deleterious mutations was significantly larger in the freshwater populations than in the marine populations ([supplementary fig. S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online) (Wilcoxon rank sum test, *N *=* *4 and 5 for marine and freshwater populations, respectively, *W *=* *0, *P *=* *0.00721). Because PROVEAN and SIFT differ in the sensitivity and specificity (see the Materials and Methods), we used conservative criteria for classifying mutations; nonsynonymous mutations were defined as deleterious and neutral mutations when both finders categorized them as deleterious and neutral mutations, respectively. Freshwater populations have significantly larger proportion of deleterious mutations than marine populations in both homozygous and heterozygous mutations ([fig. 4*A* and *B*](#evaa065-F4){ref-type="fig"}, respectively). When a homozygous mutation was counted as 1, a heterozygous mutation was counted as 0.5, and summed to reflect the number per haploid, freshwater populations had significantly larger proportion of deleterious mutations than marine populations ([fig. 4*C*](#evaa065-F4){ref-type="fig"}). The comparison between an X chromosome and autosomes showed the expected trend of more deleterious mutations being accumulated on X chromosomes than on autosomes in the majority of populations analyzed ([fig. 4*D*](#evaa065-F4){ref-type="fig"}), although *G. nipponicus* Shinji and *G. aculeatus* Shiga populations showed no difference between X chromosomes and autosomes. The intrachromosomal pattern of the proportion of deleterious mutations was unclear ([fig. 4*E*](#evaa065-F4){ref-type="fig"}).

![---Analysis of proportion of deleterious mutations in nonsynonymous mutations. Nonsynonymous mutations were predicted to be deleterious when both PROVEAN score is \<−2.5 and *P* value of SIFT is \<0.01 (see the text). (*A--C*) Comparison of the mean values (±S.E.) for homozygous SNPs (*A*), heterozygous SNPs (*B*), and both (the corrected number per haploid; *C*). Only autosomes (excluding LG9 and LG19) were used for the analysis. The results using only either PROVEAN or SIFT are shown in [supplementary figure S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online, which are qualitatively similar to results shown in this figure. (*D*) Comparison between the ancestral-X (LG19) and autosomes (all linkage groups except LG9 and LG19). (*E*) Distribution of deleterious mutations along chromosomal positions. The analysis was performed in the same way as in [figure 3*C*](#evaa065-F3){ref-type="fig"} using the corrected number per haploid as a response variable.](evaa065f4){#evaa065-F4}

Finally, we analyzed two additional individuals collected from Gifu to confirm that our results are not substantially changed by the use of different individuals within a population. Two additional individuals showed qualitatively similar results in all analyses ([supplementary fig. S6](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Characterization of Putatively Deleterious Mutations
----------------------------------------------------

We next investigated whether these putatively deleterious alleles in freshwater populations are segregating in the ancestral source marine populations as standing allelic variation. Analysis of ten Japanese Pacific Ocean marine individuals (20 chromosomes in total) showed that the majority of the putatively deleterious alleles (\>50%) were not found in the marine populations, suggesting that these mutations are either de novo mutations in freshwater populations or rare in the ancestral marine populations (black bars in [fig. 5*A*](#evaa065-F5){ref-type="fig"}). The frequencies of freshwater-private alleles were higher in deleterious mutations than in neutral mutations (compare the heights of black bars between [fig. 5*A* and *B*](#evaa065-F5){ref-type="fig"}), suggesting that deleterious mutations may be efficiently removed by purifying selection in the marine populations. We next investigated the number of freshwater populations sharing each nonsynonymous mutation ([fig. 5*C*](#evaa065-F5){ref-type="fig"}). The number of sharing populations was significantly smaller in deleterious mutations than in neutral mutations: mean ± S.E. of the number of populations sharing deleterious and neutral mutations is 1.35 ± 0.01 and 1.89 ± 0.01, respectively (*χ*^2^ test, df = 1, *P *\<* *2.2 × 10^−16^). This supports the idea that certain portions of deleterious mutations are removed even in the freshwater populations. Because the majority of deleterious mutations (79%) were unique in each freshwater population, the accumulation of deleterious mutations we found here might occur in each freshwater population independently rather than in the common ancestor.

![---Overlap of deleterious mutations among populations. (*A*) Proportion of freshwater alleles that are deleterious (*A*) or neutral (*B*) nonsynonymous mutations in a marine population. Twenty haplotypes of Akkeshi marine population of *G. aculeatus* were used. Black, freshwater alleles not found in the 20 haplotypes; dark gray, alleles that are polymorphic in the marine population; light gray, alleles found in all of the 20 haplotypes. (*C*) Histogram of the number of freshwater populations sharing deleterious (black) or neutral (gray) nonsynonymous mutations.](evaa065f5){#evaa065-F5}

To investigate whether negative selection is relaxed in genes with deleterious mutations, we calculated *d*~N~/*d*~S~ on the branch of one freshwater population, Gifu-Tsuya, and that of one marine population, Akkeshi *G. aculeatus*, and compared *d*~N~/*d*~S~ between genes with and without deleterious mutations that were identified in a freshwater population (Gifu-Tsuya) but not in Akkeshi *G. aculeatus*. Genes with deleterious mutations had significantly higher *d*~N~/*d*~S~ in both the freshwater lineage ([supplementary fig. S7](#sup1){ref-type="supplementary-material"}*A* and *B*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online; Mann--Whitney *U* test, *W *=* *834170, *P *\<* *3.1 × 10^−10^, *N *=* *2811 and 701 for genes without and with deleterious mutations, respectively) and the marine lineage (Akkeshi *G. aculeatus*) (Mann--Whitney *U* test, *W *=* *652280, *P *=* *2.1 × 10^−3^, *N *=* *2357 and 601 for genes without and with deleterious mutations, respectively). Additionally, we found significantly higher *d*~N~/*d*~S~ in the Gifu-Tsuya freshwater lineage compared with the Akkeshi *G. aculeatus* lineage regardless of whether genes are deleterious (Mann--Whitney *U* test, *W *=* *274030, *P *=* *2.2 × 10^−16^, *N *=* *701 and 601 for Gifu-Tsuya and Akkeshi *G. aculeatus*, respectively) or neutral (Mann--Whitney *U* test, *W *=* *4163300, *P *\<* *2.2 × 10^−16^, *N *=* *2811 and 2357 for Gifu-Tsuya and Akkeshi *G. aculeatus*, respectively). Considering the possible violation of statistical independence between *d~N~*/*d~S~* of nearby physically linked genes, we calculated the median *d~N~*/*d~S~* of each category (i.e., genes with and without deleterious mutations; see above) separately for each chromosome under the assumption that genes on different chromosomes are independent. We confirmed that genes with deleterious mutations had significantly higher *d*~N~/*d*~S~ in both the freshwater lineage (Exact Wilcoxon signed-rank test with *N *=* *21 chromosomes, *V *=* *16, *P *=* *1.6 × 10^−4^) and Akkeshi *G. aculeatus* (*N *=* *21, *V *=* *44, *P *=* *0.011) and *d*~N~/*d*~S~ was significantly higher in the Gifu-Tsuya freshwater lineage than in the Akkeshi *G. aculeatus* lineage regardless of whether genes are deleterious (*N *=* *21, *V *=* *231, *P *=* *9.5 × 10^−7^) or neutral (*N *=* *21, *V *=* *230, *P *=* *1.9 × 10^−6^) ([supplementary fig. S7](#sup1){ref-type="supplementary-material"}*C* and *D*, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Despite the fact that the majority of deleterious mutations are unique in each freshwater population (only 21% of deleterious mutations were shared), we found that a relatively large number of genes possess deleterious mutations in multiple freshwater populations: 42% of genes with deleterious mutations were shared by at least two freshwater populations ([supplementary fig. S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). This suggests a possibility that deleterious mutations are likely to occur independently in the same set of genes, even though deleterious mutations themselves are not shared. To test this possibility, we calculated the probability that independent deleterious mutations randomly sampled from two freshwater populations occur on the same genes and compared it with the probability calculated similarly using neutral mutations. Nine out of ten population pairs showed a significantly higher probability for deleterious mutations than that for neutral mutations ([supplementary fig. S9](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). These results suggest that deleterious mutations tend to accumulate on the same set of genes even in different populations.

One hundred forty genes have deleterious mutations in all freshwater populations examined ([supplementary fig. S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Of the 140 genes, only two genes (ENSGACG00000007767 and ENSGACG00000016736) were unique in the freshwater populations (i.e., the remaining 138 genes had deleterious mutations in at least one of the marine populations). ENSGACG00000007767 is a novel gene homologous to integrin alpha subunit. Aizu, Nasu, and Ono, but not other two freshwater populations, shared the same deleterious mutations. ENSGACG00000016736 encodes FAT atypical cadherin 4 (FAT4). All deleterious mutations on this gene had independent origins in each population ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online).

Discussion
==========

Origin of Freshwater Populations in the Central Part of the Japanese Mainland
-----------------------------------------------------------------------------

Our data showed that freshwater colonization in central Honshu, Japan, likely occurred at least twice. Two independent freshwater colonization events may reflect two different waves of southward dispersal of marine threespine stickleback during different glacial periods. Coldwater organisms expanded their distributions southward during glacial periods, and some populations have remained as relic populations during the interglacial periods ([@evaa065-B105]; [@evaa065-B34]). Because the contemporary ambient temperature of central Honshu is too high for sticklebacks ([@evaa065-B75]; [@evaa065-B51]), they are only able to survive in spring-fed habitats where colder temperatures are maintained by groundwater flow. Clearly the Gifu and Shiga "Hariyo" populations split from the main *G. aculeatus* lineage further in the past than the Nasu, Aizu, and Ono populations. Previous estimates of divergence for the Hariyo populations have suggested a split from the marine ancestor around 0.37--0.43 Ma based on the analysis of partial sequences of the mitochondrial cytochrome *b* gene ([@evaa065-B104]). Although our phylogenetic analysis is consistent with a split after the divergence between *G. aculeatus* and *G. nipponicus* which occurred 0.89 Ma ([@evaa065-B91]), more precise divergence time estimation is now necessary using genome-wide sequences of multiple individuals.

Our PSMC analysis showed that the Gifu and Shiga populations showed expansion of population sizes during the last glacial periods (10,000--70,000 years ago). This is consistent with the presence of a large freshwater fluvial environment in the Nobi Plain ([@evaa065-B104]; [@evaa065-B103]), where the Gifu population is distributed today. In contrast, the habitats of Ono, Nasu, and Aizu populations are small basins located at relatively high elevations (\>100--200 m above the present sea levels) and surrounded by high mountains, suggesting that their distributions might have been relatively restricted to small areas even during the glacial periods compared with the Gifu and Shiga populations. Consistent with this, the Ono, Nasu, and Aizu populations substantially declined in effective population size during the last glacial period.

Reduction in Genetic Diversity in Landlocked Populations
--------------------------------------------------------

A reduction in overall genetic diversity in the Japanese freshwater stickleback populations compared with the marine populations is consistent with previous studies on threespine sticklebacks in other geographical regions ([@evaa065-B108]; [@evaa065-B65]; [@evaa065-B39]; [@evaa065-B19]; [@evaa065-B44]; [@evaa065-B9]; [@evaa065-B18]; [@evaa065-B21]). Overall reduction in genetic diversity will reduce standing genetic variation, a source for adaptive evolution ([@evaa065-B3]), and therefore can increase the risk of extinction of a population when it is faced with environmental change ([@evaa065-B24]).

Furthermore, we observed variation in heterozygosity among genomic regions. First, X chromosomes showed a larger reduction in heterozygosity compared with autosomes ([fig. 3*B*](#evaa065-F3){ref-type="fig"}). This pattern can be explained by the smaller effective population size of the X chromosome compared with autosomes ([@evaa065-B102]; [@evaa065-B68]). Second, within each chromosome, chromosome centers tend to have lower genetic diversity compared with the peripheries. This can be explained by the fact that the stickleback exhibits lower recombination rates at chromosomal centers ([fig. 3*C*](#evaa065-F3){ref-type="fig"}). An unexpected result was a difference in the patterns of proportion of heterozygosity near the extreme chromosomal end between the marine and freshwater populations. The moderate reduction of heterozygosity at the further end of chromosome in the marine populations is consistent with the recombination rate map, which shows reduction in recombination rates at the further end of chromosome ([@evaa065-B26]; [@evaa065-B94]; [@evaa065-B97]). In contrast, freshwater populations showed a trend of gradual increase in the heterozygosity toward the chromosomal ends ([fig. 3*C*](#evaa065-F3){ref-type="fig"} and [supplementary fig. S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), although we are unsure what might have caused this difference.

Accumulation of Deleterious Mutations in Freshwater Populations
---------------------------------------------------------------

We also showed that isolated freshwater populations accumulated deleterious mutations, despite the possible purging effects of inbreeding. The Aizu, Nasu, and Ono populations not only have lower genetic diversity but also carry more putatively deleterious mutations than the Gifu and Shiga populations. Aizu, Nasu, and Ono populations belong to a phylogenetic clade different from that of the Gifu and Shiga populations ([fig. 2*A*](#evaa065-F2){ref-type="fig"}) and have several unique phenotypic characteristics, such as the predominance of the partially plated morph in the Nasu population ([@evaa065-B112]). Our present genomic data indicate that these unique freshwater stickleback populations are severely endangered, and any efforts to prevent further reduction in effective population size would be necessary for conserving these populations.

Comparison of deleterious mutations among chromosomes within individuals showed that X chromosomes have higher proportions of deleterious mutations than autosomes in the majority of populations examined. This is consistent with the theoretical prediction that X chromosomes have lower effective population sizes than autosomes ([@evaa065-B102]; [@evaa065-B68]). In contrast, we found no clear intragenomic patterns of mutation loads despite the fact that the recombination rates considerably vary among chromosomal regions with the lowest levels at the center within chromosomes ([@evaa065-B93]; [@evaa065-B26]; [@evaa065-B94]; [@evaa065-B97]); this is also confirmed by our result of the proportion of heterozygous sites. Several previous studies in *Drosophila* and primates have also shown that with the exception of regions with no recombination at all, there is no significant correlation between recombination rate and nonsynonymous mutation rate ([@evaa065-B31]; [@evaa065-B7]). Considering the inconsistency between the pattern of recombination and mutation loads, heterogeneity in deleterious mutation accumulation is likely influenced by confounding factors, such as variation in background mutation rates and the genomic location of adaptive alleles that might drive an increase in frequency of deleterious mutations via hitchhiking. Further investigation of background mutation rates and genomic location of adaptive alleles is necessary and possible using genomic tools and experimental approaches ([@evaa065-B67]; [@evaa065-B64]; [@evaa065-B83]). However, it should also be noted that current recombination rate maps in sticklebacks are based on a small number of individuals and not very precise, which may be a reason why we could not find a significant correlation.

Here, we used only bioinformatic methods for predicting deleterious mutations, so we cannot exclude the possibility that nonsynonymous mutations predicted to be deleterious are actually neutral or even adaptive for freshwater populations. For example, we found that deleterious mutations occurred on the two genes (ENSGACG00000007767 and ENSGACG00000016736) in different freshwater populations but not in the marine population. Mutations of these genes may be adaptive for freshwater residency, or these genes are not important for freshwater residency and therefore under relaxed selection. However, because we showed that deleterious mutations are less likely shared among freshwater populations than nonsynonymous mutations that are predicted to be neutral ([fig. 5](#evaa065-F5){ref-type="fig"}), we suggest that the majority of mutations predicted to be deleterious likely cause reductions in fitness and have been purged to some extent in each population. Furthermore, we have shown that *d*~N~/*d*~S~ is higher in genes with putatively deleterious mutations ([supplementary fig. S7](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), supporting the idea that relaxed negative selection increases the accumulation of deleterious mutations. Interestingly, many individual deleterious mutations are not shared but the genes containing mutations have a higher probability of being shared among freshwater populations ([fig. 5*C*](#evaa065-F5){ref-type="fig"} and [supplementary fig. S8](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Because negative selection is likely to independently purge deleterious mutations on the same genes, genes under relaxed negative selection in freshwater environment may carry independent deleterious mutations that have escaped purging.

In conclusion, we showed that whole-genome sequencing of endangered populations can inform us of the accumulation of deleterious mutations. This information will help to infer which populations are the most severely endangered. Furthermore, information on the regional variations in deleterious mutation loads across the genome can give insight into not only sex chromosome evolution but also genomics of adaptation and speciation. In genomic analysis of adaptation and speciation, regions with low genetic diversity within a population and/or high genetic differentiation between populations are often identified as candidate regions contributing to local adaptation and reproductive isolation ([@evaa065-B79]; [@evaa065-B90]; [@evaa065-B33]). However, negative selection against deleterious mutations can also reduce within population genetic diversity (i.e., via background selection) and also inflate statistics of genetic differentiation such as *F*~ST~ ([@evaa065-B78]; [@evaa065-B16]; [@evaa065-B90]). The chromosomal distribution of deleterious mutations has been relatively underrepresented in genomic analysis of natural populations thus far but could act as a proxy for quantifying the strength of negative background selection across the genome. An important future direction is to investigate how well the distribution of the bioinformatically identified deleterious mutations reflect the strength of negative selection and the level of association with regions that are identified as targets of selection.

Supplementary Material
======================

[Supplementary data](#sup1){ref-type="supplementary-material"} are available at *Genome Biology and Evolution* online.

Supplementary Material
======================

###### 

Click here for additional data file.

We thank all members of Kitano Lab and Dr Akihisa Kitamura for discussion and Asano Ishikawa for providing the map figure of the Japanese archipelago. This research was supported by Ministry of Education, Culture, Sports, Science and Technology Grant-in-Aid for Scientific Research on Innovative Areas to J.K. (JP23113007 and JP23113001) and to A.T. (221S0002), Japan Society for the Promotion of Science Grant-in-Aid for Scientific Research (A) to J.K. (15H02418 and 19H01003), and NIG-JOINT to S.M. (2011-A69, 2012-A63, and 2013-A34).

**Data deposition:** This project has been deposited at Dryad under the accession (dryad.hmgqnk9cs).
