Accumulation of Deleterious Mutations in Landlocked Threespine Stickleback Populations.

Colonization of new habitats often reduces population sizes and may result in the accumulation of deleterious mutations by genetic drift. Compared with the genomic basis for adaptation to new environments, genome-wide analysis of deleterious mutations in isolated populations remains limited. In the present study, we investigated the accumulation of deleterious mutations in five endangered freshwater populations of threespine stickleback (Gasterosteus aculeatus) in the central part of the mainland of Japan. Using whole-genome resequencing data, we first conducted phylogenomic analysis and confirmed at least two independent freshwater colonization events in the central mainland from ancestral marine ecotypes. Next, analyses of single nucleotide polymorphisms showed a substantial reduction of heterozygosity in freshwater populations compared with marine populations. Reduction in heterozygosity was more apparent at the center of each chromosome than the peripheries and on X chromosomes compared with autosomes. Third, bioinformatic analysis of deleterious mutations showed increased accumulation of putatively deleterious mutations in the landlocked freshwater populations compared with marine populations. For the majority of populations examined, the frequencies of putatively deleterious mutations were higher on X chromosomes than on autosomes. The interpopulation comparison indicated that the majority of putatively deleterious mutations may have accumulated independently. Thus, whole-genome resequencing of endangered populations can help to estimate the accumulation of deleterious mutations and inform us of which populations are the most severely endangered. Furthermore, analysis of variation among chromosomes can give insights into whether any particular chromosomes are likely to accumulate deleterious mutations.

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 (Kimura 1983; Ohta 1992). The accumulation of deleterious mutations results in high mutation load (Agrawal and Whitlock 2012) and can put small populations at the risk of extinction (Lande 1994; Lynch et al. 1995a, 1995b). 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 (Kirkpatrick and Jarne 2000; Hedrick and Garcia-Dorado 2016). 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 (Frankham et al. 2010; Kardos et al. 2016).

Reduction in effective population sizes can occur in populations that colonize isolated habitats (Whittaker and Fernández-Palacios 2007), fragmented populations (Fahrig 2003), populations at the front of a range expansion (Peischl et al. 2013; Peischl and Excoffier 2015), and domesticated strains (Makino et al. 2018). 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 (Frankham et al. 2010). Although deleterious mutation load has been empirically investigated by estimates of inbreeding depression and lethal equivalents (Lynch and Walsh 1998), 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 (Allendorf 2017; Kardos et al. 2016). 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 (Ng and Henikoff 2006; Choi et al. 2012; Choi and Chan 2015). Therefore, whole-genome sequence data can help us to infer the levels of accumulation of putatively deleterious mutations (Kardos et al. 2016; Allendorf 2017).

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 (Vicoso and Charlesworth 2009; Mank et al. 2010), 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 (Smith and Haigh 1974; Charlesworth et al. 1993; Charlesworth 2009; Gossmann et al. 2011).

Genomic analysis has revealed an excess of deleterious mutations in human populations that have recently expanded to new areas (Lohmueller et al. 2008; Tennessen et al. 2012; Fu et al. 2013). Accumulation of deleterious mutations in marginal populations has been also reported in Arabidopsis thaliana (Cao et al. 2011). Accumulation of deleterious mutations has been also revealed in domesticated organisms (Makino et al. 2018), such as rice (Lu et al. 2006), sunflowers (Renaut and Rieseberg 2015), and dogs (Marsden et al. 2016). 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 (Benazzo et al. 2017; Perrier et al. 2017; Ferchaud et al. 2018).

Here, we investigated mutation loads in landlocked freshwater populations of threespine stickleback (Gasterosteus aculeatus) (fig. 1). Pleistocene glacial cycles created many freshwater habitats, where ancestral marine or anadromous threespine sticklebacks colonized and diversified (Wootton 1976, 1984; Bell and Foster 1994). The genetic and genomic basis for adaptation has been extensively investigated in freshwater threespine stickleback populations (Kingsley and Peichel 2007; Kitano et al. 2010; Jones, Grabherr, et al. 2012; Hendry et al. 2013; Savolainen et al. 2013; Peichel and Marques 2017; Ishikawa et al. 2019). For example, selection on standing allelic variation, which is maintained in the ancestral marine populations by migration-selection balance (Schluter and Conte 2009) or balancing selection (Guerrero and Hahn 2017), has often played important roles in adaptation to freshwater habitats (Colosimo et al. 2005; Miller et al. 2007; Kitano et al. 2008, 2010; Jones, Chan, et al. 2012; Jones, Grabherr, et al. 2012; Bassham et al. 2018; Nelson and Cresko 2018; Haenel et al. 2019). There are also cases of freshwater adaptation caused by de novo mutations (Chan et al. 2010; O’Brown et al. 2015; Ishikawa et al. 2019; Xie et al. 2019). In contrast to the genetic basis for adaptation, the mutation loads of freshwater threespine stickleback populations have not yet been investigated.

. 1.

—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.

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 (Wootton 1984), so the extant distribution is limited to several spring-fed ponds and streams (Kitano and Mori 2016), where water temperatures are maintained around 10–20 °C by cold spring water (Kawamata 1980; Hirai et al. 1984; Mori 1985, 1987, 1994). Furthermore, these habitats are geographically isolated from the sea, so no gene flow is likely to occur from marine fish (Kitano and Mori 2016). 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 (Mori 1997). The remaining few populations are endangered and listed as local populations at the risk of extinction by the Ministry of the Environment Japan (Hosoya 2015). Because most of these freshwater populations have unique morphological and physiological characteristics (Goto and Mori 2003; Yamasaki et al. 2019; Ishikawa and Kitano 2020), it is important to conserve them (Mori 1997; Kitano and Mori 2016).

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, 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 (Yoshida et al. 2019). 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 (Yoshida et al. 2014; Ishikawa et al. 2017; Kusakabe et al. 2017; Ravinet et al. 2018). 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 Roesti et al. (2013).

WGRS data of the four marine populations were also reported previously (Yoshida et al. 2014; Ishikawa et al. 2017; Ravinet et al. 2018): 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) (Feulner et al. 2015) and a marine population of Gasterosteus wheatlandi (DRA001086) (Yoshida et al. 2014; Ravinet et al. 2018).

Phylogenetic and Pairwise Sequentially Markovian Coalescent Analyses

Phylogenetic and pairwise sequentially Markovian coalescent (PSMC) analyses (Li and Durbin 2011) were conducted as described previously (Ravinet et al. 2018). 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 (Li et al. 2009; Danecek and McCarthy 2017) using the consensus caller. With a consensus vcf of all positions from across all individuals, we then used a previously published python script (Martin et al. 2013) to calculate maximum likelihood trees in RAxML (Stamatakis 2006) 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. S1A, 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. S1B, 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 2Ne). 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 (Guo et al. 2013).

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 (Yoshida et al. 2014, 2017). We calculated the proportion of genotyped loci that are heterozygous, which is equivalent to multiple-locus heterozygosity in Kardos et al. (2016). 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) (Kitano et al. 2009). 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 (R Core Team 2013). 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 (Wickham 2016), 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 (Yoshida et al. 2017). SIFT is reported to have higher sensitivity, but lower specificity than PROVEAN (Choi et al. 2012), so we used both methods, as conducted previously (Renaut and Rieseberg 2015; Yoshida et al. 2017).

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 (Yoshida et al. 2014; Ravinet et al. 2018). 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 (dN/dS) between genes with and without putatively deleterious mutations. Genes under relaxed selection are predicted to have dN/dS higher than genes under strong purifying selection but not >1 (Ohta 1993). dN/dS analysis was conducted as described previously with a few modifications (Yoshida et al. 2014). 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. dN/dS 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 (Yang 2007) 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, Pcolocalization, was given as , where NG,A and NG,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 Pcolocalization 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), which is consistent with previous phylogenetic analyses using allozyme (Higuchi et al. 1996), microsatellite (Cassidy et al. 2013; Ravinet et al. 2014), and Restriction-site Associated DNA Sequencing data (Ishikawa et al. 2019). 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. S1A, Supplementary Material online).

. 2.

—(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, Supplementary Material online and Ravinet et al. (2018).

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). 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. S1B, Supplementary Material online).

To investigate demographic history of these Japanese freshwater populations, we conducted PSMC analysis (fig. 2). 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 (Ravinet et al. 2018) (supplementary fig. S2, Supplementary Material online).

Reduction in Genetic Diversity in Landlocked Populations

Freshwater individuals have significantly smaller proportion of heterozygous sites than marine populations (fig. 3) (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). 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, 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, Supplementary Material online; see LG3, LG9, and LG10 in Gifu-Tsuya and LG1, LG7, and LG15 in Shiga).

. 3.

—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, 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.

Threespine stickleback has a general trend of recombination rates being lower and higher at the center and both ends of each chromosome, respectively (Roesti et al. 2013; Glazer et al. 2015; Pritchard et al. 2017; Sardell et al. 2018; Shanfelter et al. 2019). 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). 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. S4A, Supplementary Material online). This contrasted with freshwater populations, which often have their highest values at the extreme ends of chromosomes (supplementary fig. S4B, Supplementary Material online). Recent mapping of recombination rates in threespine stickleback has indicated moderately reduced recombination rates at the further ends of the chromosomes (Glazer et al. 2015; Sardell et al. 2018; Shanfelter et al. 2019), which is similar to the patterns of the proportion of heterozygous sites in the marine populations (supplementary fig. S4A, 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, 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, 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). 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), 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).

. 4.

—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, 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 using the corrected number per haploid as a response variable.

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, 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). The frequencies of freshwater-private alleles were higher in deleterious mutations than in neutral mutations (compare the heights of black bars between fig. 5), 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). 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.

. 5.

—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.

To investigate whether negative selection is relaxed in genes with deleterious mutations, we calculated dN/dS on the branch of one freshwater population, Gifu-Tsuya, and that of one marine population, Akkeshi G. aculeatus, and compared dN/dS 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 dN/dS in both the freshwater lineage (supplementary fig. S7A and B, 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 dN/dS 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/d of nearby physically linked genes, we calculated the median d/d 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 dN/dS 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 dN/dS 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. S7C and D, 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, 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, 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, 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, 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 (Watanabe and Takahashi 2009; Hannah 2015). Because the contemporary ambient temperature of central Honshu is too high for sticklebacks (Mori 1997; Kitano and Mori 2016), 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 (Watanabe et al. 2003). Although our phylogenetic analysis is consistent with a split after the divergence between G. aculeatus and G. nipponicus which occurred 0.89 Ma (Ravinet et al. 2018), 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 (Watanabe et al. 2003; Watanabe and Mori 2008), 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.

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 (Withler and McPhail 1985; Mäkinen et al. 2008; Hohenlohe et al. 2010; DeFaveri et al. 2011; Jones, Chan, et al. 2012; Cassidy et al. 2013; DeFaveri and Merilä 2015; Ferchaud and Hansen 2016). Overall reduction in genetic diversity will reduce standing genetic variation, a source for adaptive evolution (Barrett and Schluter 2008), and therefore can increase the risk of extinction of a population when it is faced with environmental change (Frankham et al. 2010).

Furthermore, we observed variation in heterozygosity among genomic regions. First, X chromosomes showed a larger reduction in heterozygosity compared with autosomes (fig. 3). This pattern can be explained by the smaller effective population size of the X chromosome compared with autosomes (Vicoso and Charlesworth 2009; Mank et al. 2010). 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). 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 (Glazer et al. 2015; Sardell et al. 2018; Shanfelter et al. 2019). In contrast, freshwater populations showed a trend of gradual increase in the heterozygosity toward the chromosomal ends (fig. 3 and supplementary fig. S4, Supplementary Material online), although we are unsure what might have caused this difference.

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) and have several unique phenotypic characteristics, such as the predominance of the partially plated morph in the Nasu population (Yamasaki et al. 2019). 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 (Vicoso and Charlesworth 2009; Mank et al. 2010). 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 (Roesti et al. 2013; Glazer et al. 2015; Sardell et al. 2018; Shanfelter et al. 2019); 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 (Haddrill et al. 2007; Bullaughey et al. 2008). 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 (Makova and Hardison 2015; Lynch et al. 2016; Peichel and Marques 2017). 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), 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 dN/dS is higher in genes with putatively deleterious mutations (supplementary fig. S7, 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 and supplementary fig. S8, 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 (Nosil 2012; Ravinet et al. 2017; Hahn 2019). 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 FST (Noor and Bennett 2009; Cruickshank and Hahn 2014; Ravinet et al. 2017). 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 are available at Genome Biology and Evolution online.

Click here for additional data file.